84
Help us understand the problem. What are the problem?

posted at

updated at

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

はじめまして

63歳ですが、プログラム書いてます。
30年以上昔に知った話を、最近、また耳にしましたので、ちょっとご紹介させていただきます。

実数計算が正しいか、9種類の言語で試しました。

今回試したのは、linux 標準の bc(計算機)、PHP、python3、ruby、perl、node.js、C++、go、fortran95 の9つの言語です。
C++ のコンパイルは gcc ru88.c -lm -o ru88.outgcc 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.gogo 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(最後の桁が3と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) は、整数であれば大きな数でも正しく計算できます。
そこで、上式を4倍して、
  ${\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 2項演算を関数に書き換え
python + decimal -0.827396059946821368141165095479816292 引数、定数を多倍長指定
ruby + decimal -0.827396059946821368e0 引数を多倍長指定
C++ + boost -0.8273960599468213681411650954798162920003 階乗関数を定義
go + math/big -0.827396059946821368141165095479816292 2項演算を関数に書き換え
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 だと、分数のまま表示することができるので、誤差0な結果が得られます。


●試したプログラム(linux で実行できるよう、一部、スクリプトを付加しています)

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

bc
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))
C
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));
}
go
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(4倍)
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(4倍)
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(4倍)
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"

Register as a new user and use Qiita more conveniently

  1. You can follow users and tags
  2. you can stock useful information
  3. You can make editorial suggestions for articles
What you can do with signing up
84
Help us understand the problem. What are the problem?