ときどきの雑記帖 RE* (新南口)
dtoa
Python で浮動小数点数を文字列化するのに使われている _Py_dg_dtoa に 次のような謎めいた数値が出てくる。
ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 +
i*0.301029995663981;
k = (int)ds;
じつはこの行の少し前に以下のようなコメントがあって、要は十進変換したときの桁数を求めたいのだけど log10 を使うのは重いので変形して近似値を求めることで代用しようということらしい。
/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
* log10(x) = log(x) / log(10)
* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
*
* This suggests computing an approximation k to log10(d) by
*
* k = (i - Bias)*0.301029995663981
* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
*
* We want k to be too large rather than too small.
* The error in the first-order Taylor series approximation
* is in our favor, so we just round up the constant enough
* to compensate for any error in the multiplication of
* (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
* and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
* adding 1e-13 to the constant term more than suffices.
* Hence we adjust the constant term to 0.1760912590558.
* (We could get a more accurate k by invoking log10,
* but this is probably not worthwhile.)
*/
このコメントから直接は読み取れないけど(だよね?)、コード作者の書いた論文 (Correctly rounded Binary-Decimal and Decimal-Binaey Conversion)を読むと コメント中の d2 に当たる演算対象の値は 1.0 ≦ d2 < 2.0 なのが前提らしい。 実際、この計算を行う前に
dval(&d2) = dval(&u);
word0(&d2) &= Frac_mask1;
word0(&d2) |= Exp_11;
のようにして指数部を調整しているのだった。
ということでdtoa.cから適当に切り出して次のようなプログラムを作って確認してみた。
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <inttypes.h>
typedef uint32_t ULong;
typedef union { double d; ULong L[2]; } U;
#define word0(x) (x)->L[1]
#define word1(x) (x)->L[0]
#define dval(x) (x)->d
#define Exp_shift1 20
#define Exp_mask 0x7ff00000
#define Exp_shift1 20
#define Frac_mask1 0xfffff
#define Exp_11 0x3ff00000
#define Bias 1023
int
main()
{
U d2;
U u;
double ds;
for (int i=0; i<64; i+=1) {
double v = 1.0 + i / 64.0;
u.d = v;
int t = (int)((word0(&u)>>Exp_shift1) & (Exp_mask>>Exp_shift1)) - Bias;
dval(&d2) = dval(&u);
word0(&d2) &= Frac_mask1;
word0(&d2) |= Exp_11;
ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + t*0.301029995663981;
printf("val=1.0+%d/64 log10=%17.15f approx=%17.15f k=%d\n", i, log10(u.d), ds);
}
}
実行結果
val=1.0+0/64 log10=0.000000000000000 approx=0.031326431754716 k=0
val=1.0+1/64 log10=0.006733382658968 approx=0.035850332607875 k=0
val=1.0+2/64 log10=0.013363961557982 approx=0.040374233461034 k=0
val=1.0+3/64 log10=0.019894828716939 approx=0.044898134314193 k=0
val=1.0+4/64 log10=0.026328938722349 approx=0.049422035167351 k=0
val=1.0+5/64 log10=0.032669116753368 approx=0.053945936020510 k=0
val=1.0+6/64 log10=0.038918066030370 approx=0.058469836873669 k=0
val=1.0+7/64 log10=0.045078374735188 approx=0.062993737726828 k=0
val=1.0+8/64 log10=0.051152522447381 approx=0.067517638579987 k=0
val=1.0+9/64 log10=0.057142886136569 approx=0.072041539433146 k=0
val=1.0+10/64 log10=0.063051745747089 approx=0.076565440286305 k=0
val=1.0+11/64 log10=0.068881289407813 approx=0.081089341139464 k=0
val=1.0+12/64 log10=0.074633618296904 approx=0.085613241992623 k=0
val=1.0+13/64 log10=0.080310751188595 approx=0.090137142845781 k=0
val=1.0+14/64 log10=0.085914628706593 approx=0.094661043698940 k=0
val=1.0+15/64 log10=0.091447117306554 approx=0.099184944552099 k=0
val=1.0+16/64 log10=0.096910013008056 approx=0.103708845405258 k=0
val=1.0+17/64 log10=0.102305044894763 approx=0.108232746258417 k=0
val=1.0+18/64 log10=0.107633878399830 approx=0.112756647111576 k=0
val=1.0+19/64 log10=0.112898118392187 approx=0.117280547964735 k=0
val=1.0+20/64 log10=0.118099312077994 approx=0.121804448817894 k=0
val=1.0+21/64 log10=0.123238951730406 approx=0.126328349671052 k=0
val=1.0+22/64 log10=0.128318477259681 approx=0.130852250524211 k=0
val=1.0+23/64 log10=0.133339278634731 approx=0.135376151377370 k=0
val=1.0+24/64 log10=0.138302698166281 approx=0.139900052230529 k=0
val=1.0+25/64 log10=0.143210032661026 approx=0.144423953083688 k=0
val=1.0+26/64 log10=0.148062535455438 approx=0.148947853936847 k=0
val=1.0+27/64 log10=0.152861418337206 approx=0.153471754790006 k=0
val=1.0+28/64 log10=0.157607853361668 approx=0.157995655643165 k=0
val=1.0+29/64 log10=0.162302974570048 approx=0.162519556496323 k=0
val=1.0+30/64 log10=0.166947879615811 approx=0.167043457349482 k=0
val=1.0+31/64 log10=0.171543631304961 approx=0.171567358202641 k=0
val=1.0+32/64 log10=0.176091259055681 approx=0.176091259055800 k=0
val=1.0+33/64 log10=0.180591760282358 approx=0.180615159908959 k=0
val=1.0+34/64 log10=0.185046101708608 approx=0.185139060762118 k=0
val=1.0+35/64 log10=0.189455220613663 approx=0.189662961615277 k=0
val=1.0+36/64 log10=0.193820026016113 approx=0.194186862468436 k=0
val=1.0+37/64 log10=0.198141399798755 approx=0.198710763321594 k=0
val=1.0+38/64 log10=0.202420197778030 approx=0.203234664174753 k=0
val=1.0+39/64 log10=0.206657250721285 approx=0.207758565027912 k=0
val=1.0+40/64 log10=0.210853365314893 approx=0.212282465881071 k=0
val=1.0+41/64 log10=0.215009325086051 approx=0.216806366734230 k=0
val=1.0+42/64 log10=0.219125891280883 approx=0.221330267587389 k=0
val=1.0+43/64 log10=0.223203803701322 approx=0.225854168440548 k=0
val=1.0+44/64 log10=0.227243781503063 approx=0.230378069293707 k=0
val=1.0+45/64 log10=0.231246523956736 approx=0.234901970146865 k=0
val=1.0+46/64 log10=0.235212711174338 approx=0.239425871000024 k=0
val=1.0+47/64 log10=0.239143004802770 approx=0.243949771853183 k=0
val=1.0+48/64 log10=0.243038048686294 approx=0.248473672706342 k=0
val=1.0+49/64 log10=0.246898469499533 approx=0.252997573559501 k=0
val=1.0+50/64 log10=0.250724877352585 approx=0.257521474412660 k=0
val=1.0+51/64 log10=0.254517866369725 approx=0.262045375265819 k=0
val=1.0+52/64 log10=0.258278015243031 approx=0.266569276118978 k=0
val=1.0+53/64 log10=0.262005887762274 approx=0.271093176972136 k=0
val=1.0+54/64 log10=0.265702033322238 approx=0.275617077825295 k=0
val=1.0+55/64 log10=0.269366987408644 approx=0.280140978678454 k=0
val=1.0+56/64 log10=0.273001272063738 approx=0.284664879531613 k=0
val=1.0+57/64 log10=0.276605396332563 approx=0.289188780384772 k=0
val=1.0+58/64 log10=0.280179856690861 approx=0.293712681237931 k=0
val=1.0+59/64 log10=0.283725137455511 approx=0.298236582091090 k=0
val=1.0+60/64 log10=0.287241711178348 approx=0.302760482944249 k=0
val=1.0+61/64 log10=0.290730039024169 approx=0.307284383797407 k=0
val=1.0+62/64 log10=0.294190571133676 approx=0.311808284650566 k=0
val=1.0+63/64 log10=0.297623746972070 approx=0.316332185503725 k=0
ふむ。
まだ論文の記述と実際のコードでどう対応しているのかよくわからない部分があるので もうちょっと粘る。
続く?