異なる浮動小数点数を異なるように表示する
浮動小数点演算の解説をしようとして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)"); }
参考文献
- フォーマット指定子一覧
- Manpage of PRINTF
- 浮動小数点演算ではまった話 - bkブログ
- http://0xcc.net/blog/archives/000164.html
- お、奥が深い!!(gkbr
- http://0xcc.net/blog/archives/000164.html
- IEEE 754 - Wikipedia
- IEEE 754-1985 - Wikipedia, the free encyclopedia
- Manpage of STRCMP
- 浮動小数点演算について
- 非正規化数 - Wikipedia
-
-
-
- -
-
-
追記
printf("%.20e", ...)とかで実用上は十分。