はじめまして
63歳ですが、プログラム書いてます。
30年以上昔に知った話を、最近、また耳にしましたので、ちょっとご紹介させていただきます。
実数計算が正しいか、９種類の言語で試しました。
今回試したのは、linux 標準の bc（計算機）、PHP、python3、ruby、perl、node.js、C++、go、fortran95 の９つの言語です。
C++ のコンパイルは gcc ru88.c -lm -o ru88.out、gcc ru88-boost.cpp -lstdc++ -o ru88-boost.out、
fortran95 のコンパイルは gfortran -O3 -o ru88.exe ru88.f95、
gfortran -O3 fmsave.o fm.o fmzm90.o fm_interval.o fm_doubleint.o fm_quadint.o fm_quadreal.o fm_rational.o -o ru88-FM.exe ru88-FM.f95、
go は go run ru88.go、go run ru88-big.go でおこないました。
事前に、必要なライブラリ等はインストールしています。
■実数の計算には、誤差がつきものです。
有名な、Rump の例題 があります。
${\displaystyle f(a,b)=333.75b^{6}+a^{2}(11a^{2}b^{2}-b^{6}-121b^{4}-2)+5.5b^{8}+{\frac {a}{2b}}}$
の、$a=77617.0,b=33096.0$ のときの値を計算するのですが、そのままプログラムを書くと正しくない値になります。
各言語における計算結果です。みごとに、バラバラですね・・・
本当の値は、最初に紹介した論文にも書いてありますが、 -0.8273960599468213～4（最後の桁が３と４の間）になります。
| 言語 | 結果 | 正誤 |
|---|---|---|
| bc | -.8273960599468214 | ○ |
| PHP | -1.1805916207174E+21 | × |
| python | 1.1805916207174113e+21 | × |
| ruby | -1.1805916207174113e+21 | × |
| perl | -1.18059162071741e+21 | × |
| JavaScript(node.js) | -1.1805916207174113e+21 | × |
| C | 576460752303423489.187500 | × |
| go | -1.1805916207174113e+21 | × |
| fortran | 6.33825300E+29 | × |
■なぜ、正しく計算できないのでしょう？
まず、最初の $333.75b^{6}$ を手計算してみましょう。
$b^{6}$ は、$1,314,174,534,371,215,466,459,037,696$ なので、$438,605,750,846,393,161,930,703,831,040$ です。
次の項のカッコの中を計算していきます。
$11a^{2}b^{2}$ は、$72,586,759,116,001,040,064$、
$b^{6}$ は、$1,314,174,534,371,215,466,459,037,696$、
$121b^{4}$ は、$145,173,518,207,904,485,376$ なので、
カッコの中は $-1,314,174,606,957,974,558,362,483,010$。
それに$a^{2}$ を掛けて $-7,917,111,779,274,712,207,494,296,632,228,773,890$ となります。
$5.5b^{8}$ は、$7,917,111,340,668,961,361,101,134,701,524,942,848$ で、
ここまでの（$333.75b^{6}+a^{2}(11a^{2}b^{2}-b^{6}-121b^{4}-2)+5.5b^{8}$）を計算すると、$-2$ になります。
途中の桁が多くて目が回りませんか？実は、コンピュータも目を回すんです。
ほとんどの CPU は、IEEE754 という規格で数値を扱っているので、16桁以上の数は正しく扱えません。
最初の $333.75b^{6}$ が、$4.3860575084639\times10^{29}$ という数になってしまいます。
次の項は $-7.9171117792747\times10^{36}$ で、ここまでを計算すると $-7.917111340669\times10^{36}$ になります。
$5.5b^{8}$ は、$7.917111340669\times10^{36}$ なのですが、これって上の計算結果と符号が違うだけで同じに見えますよね？
そして、ここまでの計算結果は $-1.1805916207174\times10^{21}$ と、全然違ったものになります（本当は $-2$）。
大きな似た数の引き算（この場合は、符号が逆な足し算）の場合、桁落ちという現象がおき、答えが正しくなくなります。
■では、どうすれば正しい結果が計算できるでしょう？
◎なるべく整数の演算にする。
python や ruby、JavaScript(node.js) は、整数であれば大きな数でも正しく計算できます。
そこで、上式を４倍して、
${\displaystyle g(a,b)=1335b^{6}+4a^{2}(11a^{2}b^{2}-b^{6}-121b^{4}-2)+22b^{8}+{\frac {2a}{b}}}$
とすると、最後の$2a\div{}b$ だけが、実数の計算になります。
| 言語 | 結果 ${\displaystyle g(a,b) \div 4}$ |
|---|---|
| python | -0.8273960599468213 |
| ruby | -0.8273960599468213 |
| JavaScript(node.js) | -0.8273960599468213 |
◎多倍長演算のライブラリを利用する
PHP には BC Math、python には decimal、ruby には dcimal、C++ には boost、go には math/big、fortran95 には FMLIB など、色々な多倍長計算のライブラリがあります。
また、他にも GMP(GNU Multiple Precision) を使ったライブラリが色々あります。
問題は、ライブラリにあわせたかたちに、プログラムを修正する必要があることです。
オブジェクト指向な言語なら演算子のオーバライドで済むのですが、定数などは多倍長であることを明示する必要がでてくる場合があります。
| 言語 | 結果 | 修正方法 |
|---|---|---|
| PHP + BC Math | -0.8273960599468214 | ２項演算を関数に書き換え |
| python + decimal | -0.827396059946821368141165095479816292 | 引数、定数を多倍長指定 |
| ruby + decimal | -0.827396059946821368e0 | 引数を多倍長指定 |
| C++ + boost | -0.8273960599468213681411650954798162920003 | 階乗関数を定義 |
| go + math/big | -0.827396059946821368141165095479816292 | ２項演算を関数に書き換え |
| fortran95 + FMLIB | -.827396059946821368141165095479816290 | 型定義が必要 |
◎有理数（分数）で計算する。
ruby には Rational、perl には bigratという、分数で計算するライブラリがあります。結果の分数は実数に変換できます。
| 言語 | 結果 | 修正方法 |
|---|---|---|
| perl + bigrat | -54767/66192 = -0.827396059946821 | bigrat 指定のみ |
| ruby (rational) | -54767/66192 = -0.8273960599468214 | 定数を Rational 指定 |
■ということで、一番簡単なのは ruby や python、perl を使うことでした。
ruby + decimal や python + decimal は、引数や定数を多倍長指定するだけでの手軽さです。
また、ruby(rational) や perl + bigrat だと、分数のまま表示することができるので、誤差０な結果が得られます。
●試したプログラム（linux で実行できるよう、一部、スクリプトを付加しています）
有効桁数の指定は、この問題での最低限にしています。
計算の内容によっては、もっと大きな桁数で計算する必要があります。
bc
#!/usr/bin/env -S bc -q
define rump(a, b)
{
return 333.75 * (b ^ 6) + (a ^ 2) * (11 * (a ^ 2) * (b ^ 2) - (b ^ 6) - 121 * (b ^ 4) - 2) + 5.5 * (b ^ 8) + a / (2 * b)
}
scale=16
rump(77617.0, 33096.0)
quit
PHP
#!/usr/bin/env php
<?php
function rump($a, $b)
{
return 333.75 * ($b ** 6) + ($a ** 2) * (11 * ($a ** 2) * ($b ** 2) - ($b ** 6) - 121 * ($b ** 4) - 2) + 5.5 * ($b ** 8) + $a / ( 2 * $b);
}
echo rump(77617.0, 33096.0), PHP_EOL;
python
#!/usr/bin/env python3
def rump(a, b):
return 333.75 * (b ** 6) + (a ** 2) * (11 * (a ** 2) * (b ** 2) - (b ** 6) - 121 * (b ** 4) - 2) + 5.5 * (b ** 8) + a / ( 2 * b)
print(rump(77617.0, 33096.0))
ruby
#!/usr/bin/env ruby
def rump(a, b)
333.75 * (b ** 6) + (a ** 2) * (11 * (a ** 2) * (b ** 2) - (b ** 6) - 121 * (b ** 4) - 2) + 5.5 * (b ** 8) + a / (2 * b)
end
puts rump(77617.0, 33096.0)
perl
#!/usr/bin/env perl
sub rump
{
($a, $b) = @_;
return 333.75 * ($b ** 6) + ($a ** 2) * (11 * ($a ** 2) * ($b ** 2) - ($b ** 6) - 121 * ($b ** 4) - 2) + 5.5 * ($b ** 8) + $a / (2 * $b);
}
print &rump(77617.0, 33096.0), "\n";
node.js
#!/usr/bin/env nodejs
function rump(a, b)
{
return 333.75 * (b ** 6) + (a ** 2) * (11 * (a ** 2) * (b ** 2) - (b ** 6) - 121 * (b ** 4) - 2) + 5.5 * (b ** 8) + a / (2 * b)
}
console.log(rump(77617.0, 33096.0))
C
#include <stdio.h>
#include <math.h>
long double rump(long double a, long double b)
{
return 333.75 * powl(b, 6) + powl(a, 2) * (11 * powl(a, 2) * powl(b, 2) - powl(b, 6) - 121 * powl(b, 4) - 2) + 5.5 * powl(b, 8) + a / (2 * b);
}
int main(int argc, char **argv)
{
printf("%Lf\n", rump((long double)77617.0, (long double)33096.0));
}
go
package main
import("fmt"; "math")
func main() {
fmt.Println(rump(77617.0, 33096.0))
}
func rump(a float64, b float64) float64 {
return 333.75 * math.Pow(b, 6) + math.Pow(a, 2) * (11 * math.Pow(a, 2) * math.Pow(b, 2) - math.Pow(b, 6) - 121 * math.Pow(b, 4) - 2) + 5.5 * math.Pow(b, 8) + a / (2 * b)
}
fortran
PROGRAM RU88
WRITE(*, *) RUMP(77617.0, 33096.0)
END PROGRAM RU88
FUNCTION RUMP(A, B)
RUMP=333.75*(B**6)+(A**2)*(11*(A**2)*(B**2)-(B**6)-121*(B**4)-2)+5.5*(B**8)+A/(2*B)
END FUNCTION RUMP
python（４倍）
#!/usr/bin/env python3
def rump4(a, b):
return 1335 * (b ** 6) + 4 * (a ** 2) * (11 * (a ** 2) * (b ** 2) - (b ** 6) - 121 * (b ** 4) - 2) + 22 * (b ** 8) + (2 * a) / b
print(rump4(77617, 33096)/4)
ruby（４倍）
#!/usr/bin/env ruby
def rump4(a, b)
1335 * (b ** 6) + 4 * (a ** 2) * (11 * (a ** 2) * (b ** 2) - (b ** 6) - 121 * (b ** 4) - 2) + 22 * (b ** 8) + (2.0 * a) / b
end
puts rump4(77617, 33096)/4
node.js（４倍）
#!/usr/bin/env nodejs
function rump4(a, b)
{
return Number(1335n * (b ** 6n) + 4n * (a ** 2n) * (11n * (a ** 2n) * (b ** 2n) - (b ** 6n) - 121n * (b ** 4n) - 2n) + 22n * (b ** 8n)) + Number(2n * a) / Number(b)
}
console.log(rump4(77617n, 33096n)/4)
PHP + BC Math
#!/usr/bin/env php
<?php
function rump($a, $b)
{
return bcadd(bcadd(bcadd(bcmul(333.75, bcpow($b, 6)), bcmul(bcpow($a, 2), bcsub(bcsub(bcsub(bcmul(bcmul(11, bcpow($a, 2)), bcpow($b, 2)), bcpow($b, 6)), bcmul(121, bcpow($b, 4))), 2))), bcmul(5.5, bcpow($b, 8))), bcdiv($a, bcmul(2, $b)));
}
bcscale(16);
echo rump(77617.0, 33096.0), PHP_EOL;
python + decimal
#!/usr/bin/env python3
import decimal
def rump(a, b):
return decimal.Decimal(333.75) * (b ** 6) + (a ** 2) * (11 * (a ** 2) * (b ** 2) - (b ** 6) - 121 * (b ** 4) - 2) + decimal.Decimal(5.5) * (b ** 8) + a / (2 * b)
decimal.getcontext().prec = 37
print(rump(decimal.Decimal(77617.0), decimal.Decimal(33096.0)))
ruby + decimal
#!/usr/bin/env ruby
require 'bigdecimal'
def rump(a, b)
333.75 * (b ** 6) + (a ** 2) * (11 * (a ** 2) * (b ** 2) - (b ** 6) - 121 * (b ** 4) - 2) + 5.5 * (b ** 8) + a / (2 * b)
end
puts rump(BigDecimal(77617.0, 5), BigDecimal(33096.0, 5))
C++ + boost
#include <iostream>
#include <boost/multiprecision/cpp_dec_float.hpp>
namespace mp = boost::multiprecision;
typedef mp::number<mp::cpp_dec_float<17>> mpf;
mpf power(mpf a, int b)
{
mpf value = (mpf)1;
for (int i = 0; i < b; i ++)
value *= a;
return value;
}
mpf rump(mpf a, mpf b)
{
return 333.75 * power(b, 6) + power(a, 2) * (11 * power(a, 2) * power(b, 2) - power(b, 6) - 121 * power(b, 4) - 2) + 5.5 * power(b, 8) + a / (2 * b);
}
int main(int argc, char **argv)
{
std::cout << rump((mpf)77617.0, (mpf)33096.0).str() << std::endl;
}
go + math/big
package main
import("fmt"; "math/big")
func main() {
fmt.Println(rump(77617.0, 33096.0))
}
func Pow(d *big.Float, e uint64) *big.Float {
result := new(big.Float).SetFloat64(1)
for i:=uint64(0); i<e; i++ {
result = new(big.Float).Mul(result, d)
}
return result
}
func rump(x float64, y float64) *big.Float {
a := new(big.Float).SetPrec(122).SetFloat64(x)
b := new(big.Float).SetPrec(122).SetFloat64(y)
return new(big.Float).Add(new(big.Float).Add(new(big.Float).Add(new(big.Float).Mul(new(big.Float).SetFloat64(333.75), Pow(b, 6)), new(big.Float).Mul(Pow(a, 2), new(big.Float).Sub(new(big.Float).Sub(new(big.Float).Sub(new(big.Float).Mul(new(big.Float).Mul(new(big.Float).SetFloat64(11), Pow(a, 2)), Pow(b, 2)), Pow(b, 6)), new(big.Float).Mul(new(big.Float).SetFloat64(121), Pow(b, 4))), new(big.Float).SetFloat64(2)))), new(big.Float).Mul(new(big.Float).SetFloat64(5.5), Pow(b, 8))), new(big.Float).Quo(a, new(big.Float).Mul(new(big.Float).SetFloat64(2), b)));
}
fortran95 + FM
PROGRAM RU88
USE FMZM
TYPE (FM) :: V, A, B
TYPE (FM), EXTERNAL :: RUMP
CHARACTER(80) :: TEXT
CALL FM_SET(26)
A=77617.0
B=33096.0
V=RUMP(A, B)
CALL FM_FORM('F40.36', V, TEXT)
WRITE(*, *) TEXT
END PROGRAM RU88
FUNCTION RUMP(A, B) RESULT(V)
USE FMZM
TYPE (FM) :: A, B, V
V=333.75*(B**6)+(A**2)*(11*(A**2)*(B**2)-(B**6)-121*(B**4)-2)+5.5*(B**8)+A/(2*B)
RETURN
END FUNCTION RUMP
perl + bigrat
#!/usr/bin/env perl
use bigrat;
sub rump
{
($a, $b) = @_;
return 333.75 * ($b ** 6) + ($a ** 2) * (11 * ($a ** 2) * ($b ** 2) - ($b ** 6) - 121 * ($b ** 4) - 2) + 5.5 * ($b ** 8) + $a / (2 * $b);
}
$value = &rump(77617.0, 33096.0);
no bigrat;
print $value, ' = ', eval($value), "\n";
ruby (rational)
#!/usr/bin/env ruby
def rump(a, b)
Rational(333.75) * (b ** 6) + (a ** 2) * (11 * (a ** 2) * (b ** 2) - (b ** 6) - 121 * (b ** 4) - 2) + Rational(5.5) * (b ** 8) + a / (2 * b)
end
value = rump(Rational(77617.0), Rational(33096.0))
print value, " = ", value.to_f, "\n"