異なる浮動小数点数を異なるように表示する

c

浮動小数点演算の解説をしようとしてCのprintfの%gでは1.0と1.0 / 6.0を6回足し合わせた数(1.0ではない)が両方「1」と表示されることに気づいて愕然としました。異なる値が同じように表示されるのは教育上よろしくないですね。
というわけでdouble専用のsprintfを作ってみました。目指すところは「とにかくビット表現が異なるdoubleは同じ見かけにしない」こと。どういう結果がかえってくるのかはテストケースを見るのが一番。

  // Zero
  TEST(0.0, "+0.0");
  TEST(-0.0, "-0.0");
  // Infinity
  TEST(1.0 / 0.0, "+Inf");
  TEST(-1.0 / 0.0, "-Inf");
  // Not a Number
  TEST(0.0 / 0.0, "+NaN(mantissa=2251799813685248)");
  // Normalized Number
  TEST(1.0, "1(+1.0 * 2^0)");
  TEST(4.25, "4.25(+1.0625 * 2^2)");
  TEST(1.1, "1.1(+1.100000000000000088817841970012523233890533447265625 * 2^0)");

  double very_small = 1.0;
  for(i=0; i<16; i++){
    very_small /= (1ULL << 63);
  }
  TEST(very_small, "3.64556e-304(+1.0 * 2^-1008)");
  very_small /= (1ULL << 14);
  TEST(very_small, "2.22507e-308(+1.0 * 2^-1022)");
  very_small /= 2;
  TEST(very_small, "1.11254e-308(+0.5 * 2^-1022)"); // Denormalized number

  double sum = 0.0;
  for(i=0; i<6; i++){
    sum += 1.0 / 6.0;
  }
  TEST(sum, "1(+1.9999999999999997779553950749686919152736663818359375 * 2^-1)");

1.1が正確に表現できない数であること(1/5が2進法では循環小数になるから)とか、1.0 / 6.0を6回足しても1.0には足りないこととかがわかります。あとNaN(Not a Number)には符号とmantissa(仮数部)が異なる複数のNaNが存在しうることも表現しています。なんで2251799813685248になっているのかとか、他の環境や別のシチュエーションではどういう値になるのか、に関してはよくわかりません。

#include <stdio.h>
#include <math.h>
#include <stdint.h>
#include <assert.h>

void sprint_double(char* result, double f)
{
  uint64_t i = *(uint64_t*)&f;
  int flag = i >> 63;
  int exp = i >> 52 & 0x7FF;
  uint64_t MANTISSA_MASK = 0xFFFFFFFFFFFFFULL;
  uint64_t mantissa = i & MANTISSA_MASK;
  int j;
#ifdef DEBUG
  // double -> bits
  char digits[64];
  for(j=0; j<64; j++){
    digits[63 - j] = '0' + (i >> j & 1);
  }
  for(j=0; j<64; j++){
    printf("%c", digits[j]);
    if(j == 0 || j == 11){
      printf(",");
    }
  }
  printf("\n");
  printf("flag: %d, exp: %d, mantissa: %llu\n", flag, exp, mantissa);
#endif

  char sign = (flag == 0) ? '+' : '-';

  if(exp == 0x7ff){
    if(mantissa){
      sprintf(result, "%cNaN(mantissa=%lld)", sign, mantissa);
    }else{
      sprintf(result, "%cInf", sign);
    }
    return;
  }

  char int_part = '1';
  if(exp == 0){
    if(mantissa == 0){
      sprintf(result, "%c0.0", sign);
      return;
    }else{
      // denormalized number
      int_part = '0';
      exp++;
    }
  }


  // mantissa -> char*
  char smant[52 + 1];
  smant[52] = '\0';
  uint64_t m = mantissa;
  for(j=0; j<52; j++){
    m *= 10;
    smant[j] = '0' + (m >> 52);
    m &= MANTISSA_MASK;
    if(m == 0){
      smant[j + 1] = '\0';
      break;
    }
  }

  exp = exp - 1023;
  printf("%g(%c%c.%s * 2^%d)\n", f, sign, int_part, smant, exp);
  sprintf(result, "%g(%c%c.%s * 2^%d)", f, sign, int_part, smant, exp);
}

int main(){
  int i;
  char buf[1024];
  #define TEST(f, s) sprint_double(buf, f); assert(strcmp(buf, s) == 0);
  // Zero
  TEST(0.0, "+0.0");
  TEST(-0.0, "-0.0");
  // Infinity
  TEST(1.0 / 0.0, "+Inf");
  TEST(-1.0 / 0.0, "-Inf");
  // Not a Number
  TEST(0.0 / 0.0, "+NaN(mantissa=2251799813685248)");
  // Normalized Number
  TEST(1.0, "1(+1.0 * 2^0)");
  TEST(2.0, "2(+1.0 * 2^1)");
  TEST(4.0, "4(+1.0 * 2^2)");
  TEST(0.5, "0.5(+1.0 * 2^-1)");
  TEST(0.25, "0.25(+1.0 * 2^-2)");
  TEST(1.25, "1.25(+1.25 * 2^0)");
  TEST(4.25, "4.25(+1.0625 * 2^2)");
  TEST(100.125, "100.125(+1.564453125 * 2^6)");
  TEST(1.1, "1.1(+1.100000000000000088817841970012523233890533447265625 * 2^0)");

  double very_small = 1.0;
  for(i=0; i<16; i++){
    very_small /= (1ULL << 63);
  }
  TEST(very_small, "3.64556e-304(+1.0 * 2^-1008)");
  very_small /= (1ULL << 14);
  TEST(very_small, "2.22507e-308(+1.0 * 2^-1022)");
  very_small /= 2;
  TEST(very_small, "1.11254e-308(+0.5 * 2^-1022)"); // Denormalized number

  double sum = 0.0;
  for(i=0; i<6; i++){
    sum += 1.0 / 6.0;
  }
  TEST(sum, "1(+1.9999999999999997779553950749686919152736663818359375 * 2^-1)");
}


参考文献

        • -

追記
printf("%.20e", ...)とかで実用上は十分。