double版 powについて

先日powのint版を作成して、通常のpowを用いた場合と比較した。
今回はint版と同じアルゴリズムを用いたdouble版powと通常のpowを比較してみる。
ただし乗数は整数に限るとする。

double pow( double x,  int exp );

かけ算のルーチンは以下の通り

double pow( double x, int exp )
{
  double a = 1.0;

  if( x == 0.0 ) return ( exp == 0 )? 1.0 : 0.0;
  if( exp < 0 ) return pow( 1.0/x, -exp );

  for( ; exp >= 1; --exp )
  {
    a *= x;
  }
  return a;
}

かけ算2のルーチンは以下の通り

double pow( double x, int exp )
{
  double a = 1.0;

  if( x == 0.0 ) return ( exp == 0 )? 1.0 : 0.0;
  if( exp < 0 ) return pow( 1.0/x, -exp );

  for( ; exp >= 1; exp >>= 1 )
  {
    if( exp & 1 )
      a *= x;
    x *= x;

  }
  return a;
}

コンパイル環境はC++ Builder XE 10.2.3で、 32bit版はbcc32とbcc32c(CLang版)、64bit版はbcc64(CLang版)を使った。
実行環境は Windows10 Pro、Intel Core i5 7600で行っている。なお、実行環境は前回のint版でも同じだ。

速度の傾向はint版と同じだ。10乗位まではかけ算とかけ算2がほぼ同等で、それ以上だとかけ算2が最も早い。bcc32版だとなぜか3乗あたりからすでにかけ算2の方が早くなる。

32ビット(bcc)版

32bit(CLang版)

64bit(CLang版)

さて、int版と違う点は浮動小数点の場合かけ算の繰り返しによる桁落ちがどの程度あるかが問題になる。
-10乗から50乗の範囲で、基数を1.0～2.0の範囲でランダムに生成し、オリジナルのpowの値との差を比較した。
各乗数毎に100万個のサンプルで計算した結果、50乗くらいになるとある程度誤差が発生するようだ。
この1例で自作関数の速度と精度について完全な証明をしているわけでもないし、なにか見落としがあるかもしれないが、実用上は今回挙げたような簡易版でも大丈夫だと考える。
乗数がintに限られるpowは割とニーズがあると思うのだが、標準化されないのだろうか?

2020.11.9追記

C03で"double pow( double, int )"も追加されたのだが、C++ Builderでは"double pow( double, double)"のエイリアスで定義されていて実益がないことが判明。残念!

以下計算結果一覧表

表. 32bit(bcc32) double

乗数	pow	かけ算	かけ算2
0	640	1172	1172
1	2125	1406	1422
2	2188	1672	1422
3	2204	1921	1687
4	2219	2219	1484
5	2250	2422	1656
6	2219	2703	1703
7	2282	2984	1937
8	2265	3188	1547
9	2375	3391	1625
10	2297	3640	1641
15	2532	4937	2140
20	2406	6172	1829
25	2516	7437	1875
30	2640	8797	2125
35	2656	10000	2079
40	2531	11265	2063
45	2765	12516	2156
50	2656	13813	2094
75	2907	20109	2359
100	2781	26391	2390

表. 32bit(bcc32c:CLang) double

乗数	pow	かけ算	かけ算2
0	610	125	93
1	2063	125	156
2	2125	156	282
3	2140	203	188
4	2156	219	375
5	2172	328	344
6	2188	359	359
7	2250	375	281
8	2157	406	453
9	2187	453	407
10	2188	468	453
15	2453	703	344
20	2313	953	500
25	2437	1125	438
30	2562	1453	407
35	2578	2000	531
40	2438	2156	593
45	2719	2375	547
50	2579	2859	609
75	2828	5266	593
100	2735	7297	687

表. 64bit(bcc64:CLang) double

乗数	pow	かけ算	かけ算2
0	1600	125	94
1	7800	156	156
2	7800	125	266
3	7800	203	219
4	7800	219	375
5	245300	297	359
6	235900	344	282
7	242200	343	282
8	228100	406	438
9	242200	406	406
10	250000	453	422
15	246900	562	360
20	253100	781	500
25	242200	875	484
30	250000	1141	375
35	248400	1422	594
40	248400	1641	594
45	248400	1906	610
50	240600	2328	578
75	248400	4547	657
100	243700	6453	688

表. 誤差確認 64bit(bcc64:CLang) double

乗数	かけ算	かけ算2
-10	0.00E+00	0.00E+00
-9	0.00E+00	0.00E+00
-8	0.00E+00	0.00E+00
-7	0.00E+00	0.00E+00
-6	0.00E+00	0.00E+00
-5	0.00E+00	0.00E+00
-4	0.00E+00	0.00E+00
-3	0.00E+00	0.00E+00
-2	0.00E+00	0.00E+00
-1	0.00E+00	0.00E+00
1	0.00E+00	0.00E+00
2	0.00E+00	0.00E+00
3	0.00E+00	0.00E+00
4	0.00E+00	0.00E+00
5	0.00E+00	0.00E+00
6	0.00E+00	0.00E+00
7	0.00E+00	0.00E+00
8	0.00E+00	0.00E+00
9	0.00E+00	0.00E+00
10	0.00E+00	0.00E+00
15	0.00E+00	0.00E+00
20	0.00E+00	0.00E+00
25	0.00E+00	0.00E+00
30	0.00E+00	0.00E+00
35	0.00E+00	0.00E+00
40	0.00E+00	0.00E+00
45	0.00E+00	0.00E+00
50	4.00E-06	7.02E-03