大数の計算バグ

[michidora]

式(赤字＝正答,青字＝記入式,緑字=途中式,黒字=表記数字)
C18=(10^14)-1
C18=99,999,999,999,999
や
C18=(10^15)-3
C18=999,999,999,999,997
のように、これはいいのですが…
C18=(10^15)-2
C18=1,000,000,000,000,000
正答=999,999,999,999,998
や
C18=(10^15)-1
C18=1,000,000,000,000,000
正答=999,999,999,999,999
という答えが出てしまうバグを見つけました。
他にも、
C18=(10^16)-22
C18=9,999,999,999,999,980
正答=9,999,999,999,999,978

C18=(10^16)-21
C18=10,000,000,000,000,000
正答=9,999,999,999,999,979

C18=1000000000000000-1
C18=1E+16-1　←最終的にこの表記になっている
C18=10,000,000,000,000,000
正答=9,999,999,999,999,999

概数や四捨五入にする機能だったら納得いきますがそれをするには明らかに中途半端な数字なのでバグだと思います。
この際はどうすればいいのですか?

ちなみに、10の14乗までなら普通に計算できますが、15乗からの計算がおかしくなります。

[akionnpu]

課題トラッカーに投稿すれば良いと思います。
http://user.services.openoffice.org/ja/ ... m.php?f=42
↑
詳細はこちらです。

原因は内部的にlog10で桁数を求めているようですが
桁数１５桁以上の場合うまく桁数をとれないようです。

なお、１６桁ではそもそも有効桁数１５桁のため正しい値は出ません。
それぞれ下記の結果部分の青字があるべき値です。

コードは参考程度に課題トラッカーに提出するときにでも
つけてみてください。（無責任に言ってみます
時間に余裕があれば書いてみます

【参考】
opt_bugが想定される内部ロジック、opt_updateが修正候補

コード: 全て選択

#include<stdlib.h>
#include<stdio.h>
#include<math.h>
double opt_update(double v)
{
        if(v!=0)
        {
                int keta=0;
                double vv = v;
                while(log10(vv)>14)
                {
                        keta +=14;
                        vv /= pow(10, 14);
                }
                keta += floor(log10(vv));
                double s = pow(10, 14-keta);
                v = floor((v + .5 / s)*s)/s;
        }
        return v;
}
double opt_bug(double v)
{
        if(v!=0)
        {
                double s = pow(10, 14-floor(log10(v)));
                v = floor((v + .5 / s)*s)/s;
        }
        return v;
}
int main(int argc, char *argv[])
{
        double k ;
        int i;
        printf("case 1\n");
        k = 1000000000000000.;
        for(i=0;i<30;i++)
        {
                volatile double v = k - i;
                printf("%d : %.14e : %.14e : %.14e : %.lf\n", i,  v, opt_bug(v), opt_update(v), v);
        }

        printf("case 2\n");
        k = 10000000000000000.;
        for(i=0;i<30;i++)
        {
                volatile double v = k - i;
                printf("%d : %.14e : %.14e : %.14e : %.lf\n", i,  v, opt_bug(v), opt_update(v), v);
        }
}

結果：３列目がバグ値、４列目が訂正値

$ gcc -O2 -lm -o doubletest test.c && ./doubletest
case 1
0 : 1.00000000000000e+15 : 1.00000000000000e+15 : 1.00000000000000e+15 : 1000000000000000
1 : 9.99999999999999e+14 : 1.00000000000000e+15 : 9.99999999999999e+14 : 999999999999999
2 : 9.99999999999998e+14 : 1.00000000000000e+15 : 9.99999999999998e+14 : 999999999999998
3 : 9.99999999999997e+14 : 9.99999999999997e+14 : 9.99999999999997e+14 : 999999999999997
4 : 9.99999999999996e+14 : 9.99999999999996e+14 : 9.99999999999996e+14 : 999999999999996
5 : 9.99999999999995e+14 : 9.99999999999995e+14 : 9.99999999999995e+14 : 999999999999995
6 : 9.99999999999994e+14 : 9.99999999999994e+14 : 9.99999999999994e+14 : 999999999999994
7 : 9.99999999999993e+14 : 9.99999999999993e+14 : 9.99999999999993e+14 : 999999999999993
8 : 9.99999999999992e+14 : 9.99999999999992e+14 : 9.99999999999992e+14 : 999999999999992
9 : 9.99999999999991e+14 : 9.99999999999991e+14 : 9.99999999999991e+14 : 999999999999991
10 : 9.99999999999990e+14 : 9.99999999999990e+14 : 9.99999999999990e+14 : 999999999999990
11 : 9.99999999999989e+14 : 9.99999999999989e+14 : 9.99999999999989e+14 : 999999999999989
12 : 9.99999999999988e+14 : 9.99999999999988e+14 : 9.99999999999988e+14 : 999999999999988
13 : 9.99999999999987e+14 : 9.99999999999987e+14 : 9.99999999999987e+14 : 999999999999987
14 : 9.99999999999986e+14 : 9.99999999999986e+14 : 9.99999999999986e+14 : 999999999999986
15 : 9.99999999999985e+14 : 9.99999999999985e+14 : 9.99999999999985e+14 : 999999999999985
16 : 9.99999999999984e+14 : 9.99999999999984e+14 : 9.99999999999984e+14 : 999999999999984
17 : 9.99999999999983e+14 : 9.99999999999983e+14 : 9.99999999999983e+14 : 999999999999983
18 : 9.99999999999982e+14 : 9.99999999999982e+14 : 9.99999999999982e+14 : 999999999999982
19 : 9.99999999999981e+14 : 9.99999999999981e+14 : 9.99999999999981e+14 : 999999999999981
20 : 9.99999999999980e+14 : 9.99999999999980e+14 : 9.99999999999980e+14 : 999999999999980
21 : 9.99999999999979e+14 : 9.99999999999979e+14 : 9.99999999999979e+14 : 999999999999979
22 : 9.99999999999978e+14 : 9.99999999999978e+14 : 9.99999999999978e+14 : 999999999999978
23 : 9.99999999999977e+14 : 9.99999999999977e+14 : 9.99999999999977e+14 : 999999999999977
24 : 9.99999999999976e+14 : 9.99999999999976e+14 : 9.99999999999976e+14 : 999999999999976
25 : 9.99999999999975e+14 : 9.99999999999975e+14 : 9.99999999999975e+14 : 999999999999975
26 : 9.99999999999974e+14 : 9.99999999999974e+14 : 9.99999999999974e+14 : 999999999999974
27 : 9.99999999999973e+14 : 9.99999999999973e+14 : 9.99999999999973e+14 : 999999999999973
28 : 9.99999999999972e+14 : 9.99999999999972e+14 : 9.99999999999972e+14 : 999999999999972
29 : 9.99999999999971e+14 : 9.99999999999971e+14 : 9.99999999999971e+14 : 999999999999971
case 2
0 : 1.00000000000000e+16 : 1.00000000000000e+16 : 1.00000000000000e+16 : 10000000000000000
1 : 1.00000000000000e+16 : 1.00000000000000e+16 : 1.00000000000000e+16 : 10000000000000000
2 : 1.00000000000000e+16 : 1.00000000000000e+16 : 1.00000000000000e+16 : 9999999999999998
3 : 1.00000000000000e+16 : 1.00000000000000e+16 : 1.00000000000000e+16 : 9999999999999996
4 : 1.00000000000000e+16 : 1.00000000000000e+16 : 1.00000000000000e+16 : 9999999999999996
5 : 1.00000000000000e+16 : 1.00000000000000e+16 : 1.00000000000000e+16 : 9999999999999996
6 : 9.99999999999999e+15 : 1.00000000000000e+16 : 1.00000000000000e+16 : 9999999999999994
7 : 9.99999999999999e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999992
8 : 9.99999999999999e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999992
9 : 9.99999999999999e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999992
10 : 9.99999999999999e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999990
11 : 9.99999999999999e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999988
12 : 9.99999999999999e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999988
13 : 9.99999999999999e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999988
14 : 9.99999999999999e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999986
15 : 9.99999999999998e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999984
16 : 9.99999999999998e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999984
17 : 9.99999999999998e+15 : 1.00000000000000e+16 : 9.99999999999999e+15 : 9999999999999984
18 : 9.99999999999998e+15 : 1.00000000000000e+16 : 9.99999999999998e+15 : 9999999999999982
19 : 9.99999999999998e+15 : 1.00000000000000e+16 : 9.99999999999998e+15 : 9999999999999980
20 : 9.99999999999998e+15 : 1.00000000000000e+16 : 9.99999999999998e+15 : 9999999999999980
21 : 9.99999999999998e+15 : 1.00000000000000e+16 : 9.99999999999998e+15 : 9999999999999980
22 : 9.99999999999998e+15 : 9.99999999999998e+15 : 9.99999999999998e+15 : 9999999999999978
23 : 9.99999999999998e+15 : 9.99999999999998e+15 : 9.99999999999998e+15 : 9999999999999976
24 : 9.99999999999998e+15 : 9.99999999999998e+15 : 9.99999999999998e+15 : 9999999999999976
25 : 9.99999999999998e+15 : 9.99999999999998e+15 : 9.99999999999998e+15 : 9999999999999976
26 : 9.99999999999997e+15 : 9.99999999999998e+15 : 9.99999999999998e+15 : 9999999999999974
27 : 9.99999999999997e+15 : 9.99999999999997e+15 : 9.99999999999997e+15 : 9999999999999972
28 : 9.99999999999997e+15 : 9.99999999999997e+15 : 9.99999999999997e+15 : 9999999999999972
29 : 9.99999999999997e+15 : 9.99999999999997e+15 : 9.99999999999997e+15 : 9999999999999972

[tmacchant]

コードをみるとdouble型浮動小数点の有効桁数が大体16桁から考えてこうなるのは，まあ当然ですね。

いま，MS-Officeが手持ちにないのでKingsoft　OfficeやGnumericでためしました。
=(10^16)-22
9999999999999980
倍精度浮動小数点の演算としてはまちがった結果ではなく正しい計算です。

OpenOfficeと同じ結果です。
これに対応するにはすべてのできなくはないですが多分かなり難しいです。

実数演算ではある意味仕方ないことで，このような計算をするには，多倍長計算ができる特別なソフトウェアを使用する必要があります。
MS-Office　ではどうなるのかはわかりませんが。

[tmacchant]

>OpenOfficeと同じ結果です。
>これに対応するにはすべてのできなくはないですが多分かなり難しいです。

すみません。先ほどの返信に推敲ミスがありました。

>これに対応するにはすべてのできなくはないですが多分かなり難しいです。
は
これに対応するに理論的にはできなくはないですが多分かなり難しいです。

が正しいです。

すいませんでした。<m(__)m>

[akionnpu]

tmacchantさんこんにちは＾＾

文言訂正は[編集]ボタンからできるみたいですよ♪

私のところのgnumericは結果の5列目と同じっぽいです＾＾；

基本、gnumericは未制御のため

1E16-21=9999999999999980
1E16-22=9999999999999978

となっているようです。。

デバッグ重視なのでしょうか

弊害として次の値を入力し表示できません。
999999999999995000000000000000

999999999999994953335044046848
となります。（内部的な値保持は正しい模様）

エンドユーザ向けにはCalcの端数処理が正しい気がします。
有効桁数１５桁の処理をしようとしているので。。

あとは、桁数さえ正確にとれれば
問題ないかなと。。

[tmacchant]

tmacchantさんこんにちは＾＾
文言訂正は[編集]ボタンからできるみたいですよ♪

ありがとうございます。今度から使わせていただきます。

余談をすこし

任意桁数の計算ができる簡単なツールとしてUnixでは，当たり前についてくるツールである
ｂｃ があります。

ｃｙgwin上で動かしてみました。
*******************
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'.

10^30-22
999999999999999999999999999978
****************

Gnuwin32 版もあります。これならcygwin不要です。
http://gnuwin32.sourceforge.net/packages.html

ただし， bc はコンソールプログラムなので，GUIしか使ったことのない人には使いづらいかも知れません。
しかし，逆にUNIX的なパイプとかリダイレクトを駆使して，OOのBasicとくみあわせればうまく使うことができるかもしれません。

OOoBasicはまったくいじったことないので保証はできませんが<m(__)m>