たくさんの数の平均を求める方法について

どうもこんにちは
研究でシミュレート用のプログラムを書いています

大量の数を入力し、その平均値を求めるコードを書いているのですが、
誤差ができるだけ小さくなる方法はないでしょうか
入力する数はdouble型の実数値あるいはint型の整数値で、
個数は1億程度です。
最初は1つずつ足していたのですが、整数型の場合はオーバーフローしてしまい、実数型の場合も徐々に加算する値が相対的に小さくなり、誤差が大きくなっていきました。
100万個ずつに区切って平均を求め、それを後で合計する方法も考えましたが、あまりきれいな方法になりません

なにかいい方法はないでしょうか

No.1 回答者： episteme 回答日時： 2009/10/03 16:21

処理系によっては倍長(たいてい64bit)整数をサポートしています。

# long long とか __int64 とか

この回答への補足

int64 は使用できないのです。

補足日時：2009/10/04 21:22

この回答へのお礼

締め切りました。

お礼日時：2010/05/12 21:58

No.2 回答者： usokoku 回答日時： 2009/10/03 16:45

C言語での変数の書式を知らないので、変なことを書いている場合があります。

ひとつの場合は、指数部が存在しない変数の場合には、
求めた和を保存する配列を作っておいて、
値を加える前の数値を保存しておいて、オーバーフローの割り込みが発生したときに、それまでの値を配列に保存します。そして加える値を最初の値として計算して行きます。
最後の計算として、配列の値を2のn乗で割って、オーバーフローしない桁に直して足します。

もうひとつの場合、指数部が存在する場合には、
求めた和を保存する配列を複数個作っておいて、
指数部が同じ値同士をたして行きます。指数部で桁落ちが発生するような状態になったときに、配列に保存します。
桁落ちの発生する範囲が指数部の値によって異なりますので、指数部の値によって配列を変える必要があります。ただ、この指数部の範囲が大きく変わるような状態での和を求めることは、「平均値」と呼べるものではなくなってしまう可能性が大きいです。

別法としては、指数部の持たない整数等に型変換をして、指数部が存在しない方法で計算して行く方法があります。

MS-DOSの時代ですと、多倍精度演算が各種公開されていました。当時の公開ライブラリを探すとよいものがあるかもしれません。また、多倍精度演算専用言語もあるようです。
http://search.yahoo.co.jp/search?p=%E5%A4%9A%E5% …
http://www.vector.co.jp/vpack/filearea/dos/prog/c/
どのような言語系をお使いかわかりませんが、ちょっと手直しをする程度で使えるはずです。

この回答への補足

計算量が極端に多くなる処理は使えないのです。
変数もintやdoubleなど基本のものしか使えません。
すいません。

補足日時：2009/10/04 21:23

この回答へのお礼

締め切りました。

お礼日時：2010/05/12 21:58

No.3 ベストアンサー 回答者： Tacosan 回答日時： 2009/10/03 17:57

情報落ちを避ける方法は存在します.

誰が考案した方法なんだったっけ....

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

この回答へのお礼

回答ありがとうございます。

この処理も試してみようと思います。
桁数に極端な差がある場合は厳しいかもしれませんが。。。

お礼日時：2009/10/04 21:25

No.4 回答者： usokoku 回答日時： 2009/10/04 17:31

２番です。

３番の方の「情報落ち」の意味で「桁落ち」という言葉を使っています。

加算だけならば、整数の固定長100桁１０進演算専用ライブラリをご自身で作っても、半日ぐらいでできるでしょう。１０進ですとそのままテキストかできますので、計算時間が多少かかってもよい場合(せいぜい１０回程度しか使わないで捨てるソフト)なら、考慮の対象に入れてもよいでしょう。

この回答への補足

計算時間が多少かかるのは避けたいのです。

補足日時：2009/10/04 21:25

この回答へのお礼

締め切りました。

お礼日時：2010/05/12 21:59

No.5 回答者： ICE_FALCON 回答日時： 2009/10/04 21:42

例えば入力がint型 indtという変数で

int型が32bitなら

double a[32],b;
for(i=0;i<32:i++) a[i]=0.;

if(indt&0x80000000) a[0] += indt;
else if(indt&0x40000000) a[1] += indt;
else if(indt&0x20000000) a[2] += indt;
else if(indt&0x10000000) a[3] += indt;
else if(indt&0x08000000) a[4] += indt;
else if(indt&0x04000000) a[5] += indt;
else if(indt&0x02000000) a[6] += indt;
else if(indt&0x01000000) a[7] += indt;
...
else if(indt&0x00000002) a[30] += indt;
else a[31] += indt;

for(i=0;i<32:i++) b+=a[i];

でどうにかなるんでは？
つまり入力の精度によって、合計を求める変数を変えれば桁おちしにくいかも。

この回答への補足

a[0] に1億個とか値を加算すればa[0]の値が相対的にindtより大きくなっていき、精度が落ちていくと思います。
それに高々1つの値を加算するのに条件分岐を大量に使うのは避けたいです。

補足日時：2009/10/04 23:52

この回答へのお礼

締め切りました。

お礼日時：2010/05/12 22:00

No.6 回答者： akayoroshi 回答日時： 2009/10/05 16:33

　ほぼ同じ値の実数値が多数個あるなら、最初、単純に合計を求めて、それが大きな誤差を含んでいても、それを使って平均値を求め、次に、元の各データとこの平均値との差の合計を求める。

それに先ほどの平均値をデータの個数倍したものを加えれば、もう少し精度が上がります。
　別のやり方としては、最初データを隣り合った2つずつの組にして、それぞれの組の合計を、もとのデータ数の半分の大きさの配列の各要素に入れる、というのを繰り返せば、配列の長さが1回ごとに半分になり、最後は1になって、元の全データの合計が得られる、というのはどうですか。

この回答へのお礼

回答ありがとうございます。

一度平均を求めてからあらためて平均値を求めるやり方は確かに精度が上がりました。
これも試してみようと思います。

2つずつ加算していく方法は考えましたが、処理がややこしくなりそうなので試してませんでした。
こちらも試してみようと思います。

お礼日時：2009/10/07 20:46

No.7 回答者： usokoku 回答日時： 2009/10/05 17:08

＞計算量が極端に多くなる処理は使えないのです。

だから、障害発生時のみの例外処理のある計算方法を答えたでしょう。
配列のアドレス(ポインター)の計算が増えますが、主記憶を内部レジスターに使っていたCPUが昔ありましたので、それぼとむちゃくちゃな計算方法ではありません。

＞変数もintやdoubleなど基本のものしか使えません。
１ライン
アセンブラ
は使えませんか。補助レジスター(自由に使えるのは2つか3つですが)を使えば、今のCPUは32bitなので64-98bit演算が可能。
インテル系ならばSP, BP, IPを保護して、SI, DIを入出力ポインターに割り振り、EAXからEDXの４レジスターを使って128bit演算が可能なはず。

この回答へのお礼

締め切りました。

お礼日時：2010/05/12 21:59

No.8 回答者： ICE_FALCON 回答日時： 2009/10/05 22:11

No.5 の回答者です。

１億個足すと精度が下がるとのことですが、
No.5のやりかたは、普通に足すのにくらべて31bit多く精度があります。
でも今きずいたけど、unsignedにしか対応してないな・・・。
signedに対応すると・・・。こんな感じ。
int i,j;
double a[32],b;
int indt,abs_dt;
int msk =0x4000000;

for( i=0;i<32;i++) a[i]=0.;

for (i=0;i<100000000;i++){
indt = rand();　//入力
abs_dt = (indt<0) ? ~indt : indt;

for(j=0;j<31;j++){
if(abs_dt & (msk>>j)){
a[j] += indt;
break;
}
}
}
for(i=0;i<32;i++) b+=a[i];
out = b/100000000.;

条件分岐が多いのなら、１ビット毎に分けているのを２ビットごととかにすれば、少なくなります。

この回答へのお礼

回答ありがとうございます。

もう少し検討してみようと思います。

お礼日時：2009/10/07 20:48

No.9 回答者： tatsu99 回答日時： 2009/10/06 17:22

>int64 は使用できないのです。

お使いのコンパイラとＯＳは何でしょうか。
たぶん使えると思いますが・・・・
それとも、int型かdouble型を使用しなさいという制約があるのでしょうか？

この回答への補足

> それとも、int型かdouble型を使用しなさいという制約があるのでしょうか？

そうです。

補足日時：2009/10/07 20:48

この回答へのお礼

締め切りました。

お礼日時：2010/05/12 21:59

No.10 回答者： tatsu99 回答日時： 2009/10/09 16:20

>> それとも、int型かdouble型を使用しなさいという制約があるのでしょうか？

>そうです。
それでは、以下のような方法はいかがでしょうか。
合計を求める領域をint型とdouble型で用意します。（最初に０クリアしておきます）
int型に加算を繰り返していきますが、もし加算した結果がオーバーフローする場合は、加算する前のint型の値を、double型に加算します。
そして、int型を０クリアのち、int型に加算します。
上記の処理を繰り返していきます。
全て加算した後に、最後にint型の合計をdouble型に加算します。
そうするとdouble型に全ての合計が格納されています。
尚、オーバーフローしたかどうかは、加算前と加算後の値の大小を比較すれば判ります。加算前＞加算後の場合、オーバーフローが起こっていると判断します。

この回答へのお礼

締め切りました。

お礼日時：2010/05/12 21:59