printfのa指定子が便利だ

マルチポスト元→http://exth.net/~tgbt/wordpress/2009/08/28/2497/


CPUとGPUの演算精度について調査していた（まだ途中）際に気がついたことを書いておく．
もしかしたらものすごく今更で誰でも知っている内容なのかも知れないけど，そんなことは知ったことか！あと，もしかしてもしかしたら，間違ったことを書いているかも知れないので，そのときはごめんなさい．


浮動小数型の変数の中身を(C言語の)printfで表示するための指定子としては，fやeが有名だと私は思っている．（fとeから適切な方を選んでくれるgもあるけど，ここでは割愛．）
例えば，doubleの0.1をprintfのfとeで表示すると，

  • f 0.100000
  • e 1.000000e-01

といった感じになる．
これだとあまり差を感じないかも知れないが，
0.00001234を表示しようとすると

  • f 0.000012
  • e 1.234000e-05

となるし，
0.000000001234を表示しようとすれば

  • f 0.000000
  • e 1.234000e-09

となる．


そんなわけで，とても小さな値を表示するならfよりeを使った方が良い．
（gを使えば適切な方を使ってくれるけどね．）


以上が前置き１．


演算精度について調べはじめて気がついたんだけど，最初に例に挙げた0.1という数は，実はdouble型で表現することができない．なぜなら，0.1を2進数で表現すると，
0.000110011001100110011001100110011001...
といように，循環小数になってしまうためだ．
その割にeやfの出力は綺麗に0.1を出してくれるわけだが，これはprintfの丸めの仕様らしい．（偶数方向への丸めが働いた結果か？）

ここでプログラム内で0.1同士のかけ算をしてみることにした．人間の脳的には0.01が答えになるはず．
結果は

  • f 0.010000
  • e 1.000000e-02

となり，正しい．だがよく考えてみるとこれはおかしくて，「0.1に近い循環小数」同士をかけ算したら「0.1同士の積」になってしまっている．人間的に意図した結果は得られているんだけど，計算機的には妙だ．


ここまでが前置き２．


さて，ここでa指定子．
C99の拡張のようだが，普通のGNU/Linuxでgccを使う限りは有効だろう．

double引き数を(abcdefの文字を使って) [-]0xh.hhhhp±d; 形式の16進表記に変換する。

らしい．
何はともあれ，まずは0.1を出力させてみた．

  • f 0.100000
  • e 1.000000e-01
  • a 0x1.999999999999ap-4

謎の結果が得られた．
これは途中で出した
0.000110011001100110011001100110011001...
に対応している．
始めに1が現れた位置で区切り，その後の値を4つずつ区切ると，
0.0001 1001 1001 1001 1001 1001 1001 1001 1001 ...
となる．
ここで4つずつのまとまりを16進表記にすると，
0.0001 9 9 9 9 9 9 9 9 ...
となる．あとは小数点を移動した数を末尾につけてやればOKだ．（4つ右に動かしたので-4．）


それならば，というわけで，a指定子で「0.1」，「0.01」，「0.1*0.1」それぞを出力してみた．

  • 0.1 = 0x1.999999999999ap-4
  • 0.01 = 0x1.47ae147ae147bp-7
  • 0.1*0.1 = 0x1.47ae147ae147cp-7

確かに0.01と0.1*0.1の値が別になった．
64bit doubleは仮数部が52bitあるのだが，0.1を見ると9が12個+aが1つの合計13個の文字がならんでおり，13*4=52bitがフルに見えていることになる．
ここまで見えればCPU内部での値の保持の仕方が見えて非常に面白い．特に，0.01よりも0.1*0.1の方が大きいあたりが興味深い．


というわけで，ちょっとa指定子で遊んでみた結果をまとめてみた．誰かの参考になればいいね．


おまけ：
参考のために（a形式の出力の詳細を確認するために），a形式っぽい出力を行うサンプルプログラムを書いてみた．需要はゼロだと思うけど，自分のメモのために貼っておく．

void printHex(double d)
{
  const int INTMAX = 64;
  const int DECMAX = 64;
  int i;
  char out[INTMAX*DECMAX];
  char *iout = out;
  char *dout = &out[INTMAX];
  int integer = (int)floor(d);
  long double decimal = d - (long double)integer;

  // integer
  for(i=0;i<INTMAX;i++)iout[i]=0;
  i = 1;
  while(i<=INTMAX){
    iout[INTMAX-i] = integer&1;
    integer >>= 1;
    i++;  }
  for(i=0;i<INTMAX;i++)putchar('0'+iout[i]);

  putchar('.');

  // decimal
  for(i=0;i<DECMAX;i++)dout[i]=0;
  i = 0;
  long double dcheck = 1.0;
  while(i<DECMAX){
    dcheck *= 0.5;
    if(decimal >= dcheck){dout[i] = 1; decimal-=dcheck;}
    else{dout[i] = 0;}
    i++;
  }
  for(i=0;i<DECMAX;i++)putchar('0'+dout[i]);

  putchar('\n');

  // a type format
  int found=0;
  int count=0;

  printf("printf a correspond 1: ");
  for(i=0;i<INTMAX+DECMAX;i++){
    if(found==0){
      if(out[i]==1){
        found++;
        putchar('0'+out[i]);
        putchar(' ');
      }
    }else{
        putchar('0'+out[i]);
        if(++count==4){putchar(' ');count=0;}
    }
  }
  putchar('\n');

  char hex[4];
  int pos = 0;
  found = 0;
  count = 0;
  printf("printf a correspond 2: ");
  for(i=0;i<INTMAX+DECMAX;i++){
    if(found==0){
      pos++;
      if(out[i]==1){
        found++;
        putchar('0'+out[i]);
        putchar('.');
      }
    }else{
      hex[count] = out[i];
      if(++count==4){
        printf("%x", hex[0]*2*2*2+hex[1]*2*2+hex[2]*2+hex[3]);
        count=0;
      }
    }
  }
  printf("%+d", INTMAX-pos);
  putchar('\n');

}