はじめまして

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 で実行できるよう、一部、スクリプトを付加しています）

有効桁数の指定は、この問題での最低限にしています。
計算の内容によっては、もっと大きな桁数で計算する必要があります。

ru88.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

ru88.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

ru88.py

#!/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

ru88.rb

#!/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

ru88.prl

#!/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

ru88.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))

ru88.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));
}

ru88.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

ru88.f95

      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（４倍）

ru88-4.py

#!/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（４倍）

ruby-4.rb

#!/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（４倍）

ru88-4.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

ru88-bc.php

#!/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

ru88-decimal.py

#!/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

ru88-decimal.rb

#!/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

ru88-boost.cpp

#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

ru88-big.go

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

ru88-FM.f95

      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

ru88-bigrat.prl

#!/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)

ru88-rational.rb

#!/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"

【計算結果が正しくない！？】案外知らない、計算誤差の話