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

# はじめまして
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 の例題](https://www.tuhh.de/ti3/paper/rump/Ru88a.pdf) があります。
　　${\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](https://ja.wikipedia.org/wiki/IEEE_754) という規格で数値を扱っているので、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](https://www.php.net/manual/en/book.bc.php)、python には [decimal](https://docs.python.org/3/library/decimal.html)、ruby には [dcimal](https://docs.ruby-lang.org/ja/latest/class/BigDecimal.html)、C++ には [boost](https://www.boost.org/)、go には [math/big](https://pkg.go.dev/math/big)、fortran95 には [FMLIB](http://dmsmith.lmu.build/) など、色々な多倍長計算のライブラリがあります。
また、他にも [GMP(GNU Multiple Precision)](https://gmplib.org/) を使ったライブラリが色々あります。

問題は、ライブラリにあわせたかたちに、プログラムを修正する必要があることです。
オブジェクト指向な言語なら演算子のオーバライドで済むのですが、定数などは多倍長であることを明示する必要がでてくる場合があります。

| 言語 | 結果 | 修正方法 |
|----|----|----|
| 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](https://docs.ruby-lang.org/ja/latest/class/Rational.html)、perl には [bigrat](https://perldoc.perl.org/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 で実行できるよう、一部、スクリプトを付加しています）
有効桁数の指定は、この問題での最低限にしています。
計算の内容によっては、もっと大きな桁数で計算する必要があります。

<details><summary>bc</summary>

```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
```
</details>

<details><summary>PHP</summary>

```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;
```
</details>

<details><summary>python</summary>

```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))
```
</details>

<details><summary>ruby</summary>

```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)
```
</details>

<details><summary>perl</summary>

```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";
```
</details>

<details><summary>node.js</summary>

```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))
```
</details>

<details><summary>C</summary>

```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));
}
```
</details>

<details><summary>go</summary>

```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)
}
```
</details>

<details><summary>fortran</summary>

```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
```
</details>

<details><summary>python（４倍）</summary>

```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)
```
</details>

<details><summary>ruby（４倍）</summary>

```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
```
</details>

<details><summary>node.js（４倍）</summary>

```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)
```
</details>

<details><summary>PHP + BC Math</summary>

```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;
```
</details>

<details><summary>python + decimal</summary>

```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)))
```
</details>

<details><summary>ruby + decimal</summary>

```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))
```
</details>

<details><summary>C++ + boost</summary>

```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;
}
```
</details>

<details><summary>go + math/big</summary>

```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)));
}
```
</details>

<details><summary>fortran95 + FM</summary>

```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
```
</details>

<details><summary>perl + bigrat</summary>

```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";
```
</details>

<details><summary>ruby (rational)</summary>

```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"
```
</details>