誤差を出さずに整数計算をしたいです

質問

いつもお世話になっております。

0～18446744073709551615までの整数値を計算して出力するプログラムを作成しております。
しかし、計算に誤差が出てしまい意図した出力結果が得られません。
試してみた方法と出力は以下のとおりです。
----------------------------------------------------------------------
[方法１]
/* 18446744073709551615に0.0～1.0までの乗算 */
#define RANGE_MAX 18446744073709551615ULL
#define MAGNIFICATION 1.0 /* 0.0から1.0 */
printf("%0.lf", RANGE_MAX * MAGNIFICATION);
[出力１]
0～18446744073709552000
----------------------------------------------------------------------
[方法２]
/*
562949953421311に0.0～1.0までの乗算した結果と
562949953421312に0.0～1.0までの乗算した結果を
32767回加算した結果を加算する
*/
#define RANGE_MAX 562949953421312ULL
#define MAGNIFICATION 1.0 /* 0.0から1.0 */
#define CYCLE 32767
unsigned long long int y = (RANGE_MAX - 1) * MAGNIFICATION;
int loop;
for(loop = 0; loop < CYCLE; loop++){
y += RANGE_MAX * MAGNIFICATION;
}
printf("%llu \n", y);
[出力２]
0～9223372036854775808
----------------------------------------------------------------------

方法１は丸め誤差が原因だと考え方法2を試してみました。
しかし、方法2では情報落ちになってしまうのか、途中から計算されない
模様です。
誤差を出さずに計算によって0～18446744073709551615までの整数を
得るにはどのようにすれば良いのかご教授願います。

回答No.7

sakusaker7 · ベストアンサー

あーすいません。誤差出てますね。
まあ *1.0 なら取り除いてもいいじゃないかと思うんですが。

#include <stdio.h>
#include <stdlib.h>

#define RANGE_MAX 562949953421312ULL
#define MAGNIFICATION 1.0 /* 0.0から1.0 */
#define CYCLE 32767

int
main()
{
unsigned long long int y = (RANGE_MAX - 1);
unsigned long long int y2 = y;
int loop;

for(loop = 0; loop < CYCLE; loop++){
y += RANGE_MAX;
y2 +=
#ifdef USE_CAST
(unsigned long long)
#endif
(RANGE_MAX * MAGNIFICATION);
if (loop % 2048 ==0 || loop > 32760)
printf("%05d %21llu %21llu \n", loop, y, y2);
if (y != y2) {
printf("y!=y2\n%05d %21llu %21llu\n", loop, y, y2);
break;
}
}
printf("%llu \n", y);
printf("\n\n%llu \n", ~0ULL);
printf("%llu \n", ~0ULL - RANGE_MAX);

return 0;
}

gcc でコンパイルしたものをアセンブリ言語レベルでみると

L2:
cmpl $32766, -28(%ebp)
jg L3
leal -16(%ebp), %eax
addl $0, (%eax)
adcl $131072, 4(%eax)
fildq -24(%ebp)
fstpt -56(%ebp)
cmpl $0, -20(%ebp)
jns L5
fldt LC0
fldt -56(%ebp)
faddp %st, %st(1)
fstpt -56(%ebp)
L5:
fldt -56(%ebp)
fstpl -40(%ebp)
fldl -40(%ebp)
fldl LC1
faddp %st, %st(1)
fstpl (%esp)
call ___fixunsdfdi
movl %eax, -24(%ebp)
movl %edx, -20(%ebp)
movl -28(%ebp), %eax
andl $2047, %eax
testl %eax, %eax
je L7
cmpl $32760, -28(%ebp)
jg L7
jmp L6
L7:

加算で浮動小数点命令を使ってます。そして
00000 1125899906842623 1125899906842623
y!=y2
00015 9570149208162303 9570149208162304
9570149208162303
loop が15のところで食い違いがでてます。
つまり、RANGE_MAX を16個足しこんだところで double の精度を超えちゃっているので誤差が出たと。

で、上のソースの、MAGNIFICATIONを使わない加算と、
加算の前にMAGNIFICATIONを乗算した結果を unsigned long long で
キャストするようにしたものは期待通りの値になってるんじゃないですか?
つまりは浮動小数点演算が出てくるようにしちゃだめってことです。
あるいは long double がdoubleよりも精度のよい型であるような処理系なら
キャストなどしなくてもうまく行くかもしれません。

>#define MAGNIFICATION 1.0 /* 0.0から1.0 */
→
#define MAGNIFICATION 1.0L /* 0.0から1.0 */

お礼

お礼が遅くなりまして申し訳ございません。

>9570149208162303
>loop が15のところで食い違いがでてます。
>つまり、RANGE_MAX を16個足しこんだところで double の精度を超えちゃっているので誤差が出たと。
"y += RANGE_MAX * MAGNIFICATION;"と書いてしまうと
"y = y + (unsigned long long int)(RANGE_MAX * MAGNIFICATION);"
と解釈してもらえず
"y = (double)y + (RANGE_MAX * MAGNIFICATION);"
という感じになっていたのですね。

以下のようにキャストを使用するとうまく計算できるようになりました。
y += RANGE_MAX * MAGNIFICATION;
→
y += (unsigned long long int)(RANGE_MAX * MAGNIFICATION);

ありがとうございました。

回答No.6

sakusaker7

結局はどういう結果がほしいのかという#5での補足要求になるわけですけど、
方法2について MAGNIFICATION が1であれば演算誤差は出てないですよ。
#0も誤差がないといえばないけど意味がないので除外

#include <stdio.h>
#include <stdlib.h>

#define RANGE_MAX 562949953421312ULL
#define MAGNIFICATION 1.0 /* 0.0から1.0 */
#define CYCLE 32767

int
main()
{
unsigned long long int y = (RANGE_MAX - 1) * MAGNIFICATION;
int loop;

for(loop = 0; loop < CYCLE; loop++){
y += RANGE_MAX * MAGNIFICATION;
if (loop % 2048 ==0 || loop > 32760)
printf("%05d %21llu \n", loop, y);
}
printf("%llu \n", y);
printf("\n\n%llu \n", ~0ULL);
printf("%llu \n", ~0ULL - RANGE_MAX);

return 0;
}

00000 1125899906842623
02048 1154047404513689600
04096 2306968909120536576
06144 3459890413727383552
08192 4612811918334230528
10240 5765733422941077504
12288 6918654927547924480
14336 8071576432154771456
16384 9224497936761618432
18432 10377419441368465408
20480 11530340945975312384
22528 12683262450582159360
24576 13836183955189006336
26624 14989105459795853312
28672 16142026964402700288
30720 17294948469009547264
32761 18443929323942445056
32762 18444492273895866368
32763 18445055223849287680
32764 18445618173802708992
32765 18446181123756130304
32766 0
0

18446744073709551615
18446181123756130303

0.xという数値をMAGNIFICATIONに使った場合は小数部分が掛け算で発生するので
そこで誤差が生じるというのはすでに述べたとおりです。

> この値を32768（2^15乗）回だけ加算すれば、
> 18446744073709551616になると考えたものが方法2です。

とのことですけど、MAGNIFICATIONを変化させることによってなにを得ようとしたのでしょうか?
あと念のため、18446744073709551616 だと 2^64 なのでお使いの処理系の
unsigned long long を使っても最大値を飛び越して0になっちゃいますよ。
64bit 長の unsigned な数値は 0 ～ 2^64-1 です。～2^64じゃないです。

補足

回答して頂いた皆様に誤解を与える質問をしてしまった為、
質問内容を訂正いたします。

[質問内容]
"(2^64-1) * 1.0"の計算結果が意図した結果（2^64-1）にならない。
"1.0"で乗算した時に誤差が出ないように計算したい。

"(2^64-1) * 0.x"の計算で誤差が出ないのは無理なので質問にミスがあったことをお詫びいたします。

>00000 1125899906842623
この時の変数yは(562949953421312 - 1) + (562949953421312 * 1)
と同じ値であり、意図通りです。
>02048 1154047404513689600
この時の変数yは(562949953421312 - 1) + (562949953421312 * 2049)
と違う値（1だけ大きい）であり、意図通りではないです。

"y += RANGE_MAX * MAGNIFICATION;"という計算では誤差が出る時と
出ない時があるという事なのでしょうか？

回答No.5

Oh-Orange

★補足要求です。
・そもそも何がしたいのでしょうか？
　unsigned long long型に 0.0～1.0 の小数を掛けた時点で double 型へ変換されます。
　その後に unsigned long long型に代入しても double の誤差をそのまま代入します。
・MAGNIFICATION という定数は 0.1 刻みで 1.0 までを掛けるのなら 10 で割れば。
　0.0001 刻みで 1.0000 なら 10000 で割って処理すれば良いが何がしたいのかが
　分からないので多倍長演算を組むのが適切なのかこちらでは不明です。
　どうのような目的で 0.0～1.0 を掛け、unsigned long long型の整数を何に
　利用したいのでしょうか。詳細を補足にどうぞ。

補足

>・そもそも何がしたいのでしょうか？
６４bit対応プログラムのデバッグに使用します。
0～18446744073709551615の値を正しく処理できるか判断する為に、
0～18446744073709551615の値を返すデバッグ用のプログラムを
作成しております。

極端に良いますと、0～18446744073709551615の値が返れば良いので小数の計算をする必要はありません。
しかし、意図した結果がどうしても得られなかったので質問させて頂いた次第です。

回答No.4

yphkz4063

コンパイラの性能を上回る処理をしたければ、加減乗除ですら自前で実装するのが正しい方法です。

回答No.3

jacta

一番安直な方法は、Visual C++ 2008を止めて、GCCかBorland C++ Compilerにした上でlong double型を使えば解決します。

Visual C++ 2008にこだわるのであれば、C/C++の範囲では、自分で演算処理を実装せざるを得ません。

回答No.2

sakusaker7

> /* 18446744073709551615に0.0～1.0までの乗算 */
> #define RANGE_MAX 18446744073709551615ULL

2^64 - 1 ですか。
であれば、1以下の浮動小数点数を掛けた時点でdouble にしているでしょうから、
その有効桁数(ゲタ含めて53bit)に収まらなくなるのでその分
答えが合わなくなっているというオチでしょう。

方法2も結局doubleを経由してダメですね。

補足

double型の有効桁数が10進数で約15桁と書いてあった為、2^49の"562949953421312"を使用しました。
この値を32768（2^15乗）回だけ加算すれば、
18446744073709551616になると考えたものが方法2です。

有効桁数以内でも無理なようなのでお手上げです。

参考URL
http://www.cc.kyoto-su.ac.jp/~yamada/pB/float.html

回答No.1

asuncion

直接の回答ではありません。
方法１、２とも、コードの断片でなく、
main関数の定義などを含めた全体を見せていただけますか？

それから、

> #define MAGNIFICATION 1.0 /* 0.0から1.0 */

0.0から1.0まで変化している形跡がありませんね。
まあ、コード全体を見せてくださればわかる話かもしれませんけれど。

補足

以下が方法1と方法2のソースコードです。
「#include "stdafx.h"」はVisualC++2008が自動で付けており、
内容の変更はしておりません。
----------------------------------------------------------------------
[方法１]
----------------------------------------------------------------------
#include "stdafx.h"
#include <stdio.h>
#define RANGE_MAX 18446744073709551615ULL
#define MAGNIFICATION 1.0 /* 0.0から1.0 */

int main(void)
{
printf("%0.lf", RANGE_MAX / MAGNIFICATION);
return 0;
}
----------------------------------------------------------------------
[方法２]
----------------------------------------------------------------------
#include "stdafx.h"
#include <stdio.h>
#define RANGE_MAX 562949953421312ULL
#define MAGNIFICATION 1.0 /* 0.0から1.0 */
#define CYCLE 32767

int main(void)
{
unsigned long long int y = (RANGE_MAX - 1) * MAGNIFICATION;
int loop;
for(loop = 0; loop < CYCLE; loop++){
y += RANGE_MAX * MAGNIFICATION;
}
printf("%llu \n", y);
return 0;
}
----------------------------------------------------------------------

なお、現在はMAGNIFICATIONについては手動で変更しております。