Kahan summation algorithmを理解する

はじめに

Kahan summation algorithmの考え方を説明します。
浮動小数点数列の総和を精度よく求めることができます。

Kahan summation algorithm

具体例で説明します。
浮動小数点数の和は桁数を合わせて足し算を行います。このときに下位桁が失われるのですが
この下位桁をうまく取り出して次の加算に反映することで精度を向上させます。

sample

サンプルプログラムを示します。
0~1.0をloop数で割って一定幅に分割します。分割した幅dxを合計します。
合計値と真の値を比較して誤差を計算します。

オプションkにeを与えるとKahan summation algorithmを有効にします。
e以外を与えると無効にして通常の合計を行います。

test.c

# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define FLOAT_TYPE  double

void usage()
{
  printf("usage:./a.out [loop] [k] <q>\n");
  printf(" e:enable kahan !e:disable kahan\n");
}

void kahan_sum(FLOAT_TYPE* d, int size, FLOAT_TYPE* sum, FLOAT_TYPE* c)
{
  int i=0;
  for (i=0; i<size; i++) {
    FLOAT_TYPE tmp_sum;
    if (fabs(d[i]) <= fabs(*sum)) {
        tmp_sum = *sum + (d[i] + *c);
        *c = (d[i] + *c) - (tmp_sum - *sum);
    } else {
        tmp_sum = (*sum + *c) + d[i];
        *c = (*sum + *c) - (tmp_sum - d[i]);
    }
    *sum = tmp_sum;
  }
}

int main(int argc, char* argv[])
{
  int i, loop, q, ek = 0; /* ek:enalbe kahan */
  FLOAT_TYPE max = 1.0;
  FLOAT_TYPE dx, x = 0;
  FLOAT_TYPE e = 0; /* error */
  FLOAT_TYPE  k  = 0; /* kahan c */

  if (!(argc == 3 || argc == 4)) {
    usage();
    return -1;
  }

  loop = atoi(argv[1]);
  dx = max / loop;
  if (argc == 4 && argv[3][0] == 'q') {
    q = 1;
    if (argv[2][0] == 'e') {
        ek = 1; /* enable kahan */
    } else {
        ek = 0; /* disable */
    }
  } else if (argc == 3) {
    q = 0;
    if (argv[2][0] == 'e') {
        ek = 1; /* enable kahan */
    } else {
        ek = 0; /* disable */
    }
  } else {
    usage();
    return -1;
  }

  for (i=0; i<loop; i++) {
    FLOAT_TYPE c, d; /* c:correct, d:diff */
    if (ek == 1) {
        kahan_sum(&dx, 1, &x, &k);
    } else {
        x += dx;
    }
    c = dx*(i+1);
    d = c-x;
    e = fmax(e, d);
    if (!q) {
      printf("%+.15e\t%+.15e\t%+.15e\t%+.15e\n", x, c, d, k);
    }
  }
  printf("dx\t%+.15e\te\t%+.15e\tloop\t%d\n", dx, e, loop);

  return 0;
}

コンパイルします。

console

gcc -O2 -Wall -Werror test.c -lm

実行します。
通常の合計とKahan summation algorithmそれぞれの誤差を比較します。
Kahan summation algorithmでは誤差が小さくなっていることが分かります。

console

$ cat exec.sh 
# !/bin/bash

echo "disable kahan"
./a.out 10 d q
./a.out 100 d q
./a.out 1000 d q
./a.out 10000 d q
./a.out 100000 d q
./a.out 1000000 d q
./a.out 10000000 d q

echo "enable kahan"
./a.out 10 e q
./a.out 100 e q
./a.out 1000 e q
./a.out 10000 e q
./a.out 100000 e q
./a.out 1000000 e q
./a.out 10000000 e q

$ ./exec.sh 
disable kahan
dx	+1.000000000000000e-01	e	+1.110223024625157e-16	loop	10
dx	+1.000000000000000e-02	e	+2.775557561562891e-17	loop	100
dx	+1.000000000000000e-03	e	+0.000000000000000e+00	loop	1000
dx	+1.000000000000000e-04	e	+9.381384558082573e-14	loop	10000
dx	+1.000000000000000e-05	e	+1.916244940503020e-12	loop	100000
dx	+1.000000000000000e-06	e	+6.459610624176548e-12	loop	1000000
dx	+1.000000000000000e-07	e	+2.498300455400226e-10	loop	10000000
enable kahan
dx	+1.000000000000000e-01	e	+0.000000000000000e+00	loop	10
dx	+1.000000000000000e-02	e	+0.000000000000000e+00	loop	100
dx	+1.000000000000000e-03	e	+0.000000000000000e+00	loop	1000
dx	+1.000000000000000e-04	e	+0.000000000000000e+00	loop	10000
dx	+1.000000000000000e-05	e	+0.000000000000000e+00	loop	100000
dx	+1.000000000000000e-06	e	+0.000000000000000e+00	loop	1000000
dx	+1.000000000000000e-07	e	+0.000000000000000e+00	loop	10000000

references

• https://en.wikipedia.org/wiki/Kahan_summation_algorithm