浮動小数型で大量のデータの合計を求める際に気をつけないといけないこと

##浮動小数型の合計値が正しく出力されない

浮動小数型で大量のデータ(数千万個程度)の合計を計算する際に、正しい値が出力されず詰まったのでその理由を考えてみました。
まず問題のプログラムがこちら。

sum.c

#include <stdio.h>
#include <stdlib.h>
#define N 100000000 //　= 10^8
int main() {
  float sum = 0.0;
  float *a = (float *)malloc(sizeof(float) * N);
  unsigned int i;
  // initialization
  for (i = 0; i < N; i++)
    a[i] = 1.0;

  // summation
  for (i = 0; i < N; i++)
    sum += a[i];

  printf("total: %e\n", sum);
  free(a);
  return 0;
}

このプログラムでやりたいことは、配列aのすべての要素の合計を計算すること。
今回aの中身はすべて1.0のためは$1.0\times10^8$が出力されれば正解です。
（この時点でこのプログラムの問題点に気づいた方は、この記事を読む必要はありません(笑)。）
いざ実行してみると…

$ gcc sum.c
$ ./a.out
total: 1.677722e+07

えっ。。。
想定していた結果とは異なり、$1.677722\times10^7$が出力されてしまいました。
単精度浮動小数型(float型)は約 $-1.75\times10^{-38}$ から $3.40\times10^{38}$までの値を表現できるため、オーバーフローしたわけではありません。
また$N$が小さい場合は、このプログラムでも正しい値を出力します。
では今回はなぜこのような結果を出力してしまったのでしょうか。

##原因は情報落ち
浮動小数型は、int型などの整数型とは異なり、符号部、指数部、仮数部から構成されます。

float型の場合は仮数部が23ビットであるため、有効桁数は2進数で24桁、10進数では7桁程度になります。
今回のプログラムでは、
sum←0+1
sum←1+1
sum←2+1
...
と加算され、終盤では
sum←1.677722e+07+1
という演算がなされます。
加減算の有効数字は大きい方の値が基準となり、1.677722e+07+1の計算において1は有効桁の範囲外になります。
そのため、この演算を行ってもsumの中身の値は変化していません。
また有効桁の範囲外とならなくても、大きい値と小さい値の加減算では誤差が発生し得ます。
このときに発生する誤差を情報落ちと呼びます。
このように非常に大きい値と小さい値の加減算をする際は、誤差が発生することを認識する必要があります。

##解決策は計算順序の変更
ではどうすればこの問題が解決されるでしょう。
簡単な解決策として、変数sumを単精度から倍精度に変更することが挙げられます。
倍精度浮動小数型では有効桁数が10進数で約15桁あり、$N=10^8$でも余裕を持って計算できるでしょう。
しかし$N=10^{16}$ほどのデータの合計を計算しようとしたらどうでしょう。(そもそもメモリに収まりきらないですね)
また今回のようにデータがすべて同じ場合なら良いのですが、データにムラがある場合は情報落ちが発生する可能性も考えられます。

そこで行うべきなのがそれぞれの部分和を計算し総和を求めていく方法です。
図で例を見てみましょう。
ここでは簡単のため、N=8として示しています。

このようにN=8の場合、$\mathrm{log}_28=3$ステップで計算を行い、最終的にa[0]に合計値が格納されます。
このようなアルゴリズムで合計値を計算することにより、桁数が大きく異なる数値同士の加減算が発生しにくくなります。(データに偏りがあると発生し得るが)
また同じレベルにおいてデータの依存性がないため、並列処理が可能です。
ただし、この図を見てもわかるように、これを行うと最終的に配列aがもとの形ではなくなってしまいます。
配列を壊したくない人はコピーをしてから行うと良いでしょう。

合計計算を図の形で求めたときのプログラムを見てみましょう。

sum_reduction.c

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

#define N 100000000 // 10^8
int main() {
  float *a = (float *)malloc(sizeof(float) * N);
  int i, j;
  const int L = (int)log2(N - 1) + 1; //ステップ数
  // initialization
  for (i = 0; i < N; i++)
    a[i] = 1.0;

  // summation

  int M = 1 << L; // M=2^L

  //最初のステップだけ特異処理
  for (i = 0; i < (M >> 1); i++)
    if (i + (M >> 1) < N)
      a[i] += a[i + (M >> 1)];

  for (j = 1; j <= L; j++)
    for (i = 0; i < (M >> (j + 1)); i++)
      a[i] += a[i + (M >> (j + 1))];

  printf("total: %e\n", a[0]);
  free(a);
  return 0;
}

少し複雑になりましたが、図のアルゴリズムをそのまま実装しているだけです。
また、図ではNが2の累乗のため特異処理をする必要はありませんが、現実ではデータ数は2の累乗になっていることはほぼないでしょう。
データ数以上の配列にはアクセスしないように、最初のステップだけif文を入れて特異的に処理をしています。(if文はコストが大きいので、できるだけ省略したい)
また変数iについての繰り返しの並列化もこのアルゴリズムでなら可能です。(今回はめんどくさいので並列化はしていません。)

計算結果を見てみましょう。

$ gcc sum_reduction.c
$ ./a.out
total: 1.000000e+08

すべての値の合計がちゃんと出力されていることがわかりました。
このように、数千万個程度の大量のデータの合計を計算するときには、注意しないといけないよって話でした。

めでたしめでたし。