高速な倍精度指数関数expの実装

高速な倍精度指数関数expの実装
Cybozu Labs
2011/10/1 光成滋生
x86/x64最適化勉強会2

内容
fmath.hppの紹介
最終実装コード
ベンチマーク
expの性質
パラメータの決めかた
2011/8/6 /202

fmath.hpp
floatに対する高速なlogとexp
VCやgccのライブラリに比べて数倍速い
http://github.com/herumi/fmath/
http://homepage1.nifty.com/herumi/soft/fmath.html
exp(float)のアルゴリズムは情報科学苦手の会で
話したことがある
今回のはそれのexp(double)版で，かつ改良版
2011/8/6 /203

最終実装コード
今回の目標はこのコードの意味を理解すること
2011/8/6 /204
union di {
uint64_t i;
double d;
};
double expd(double x) {
di di;
di.d = x * 2954.639443740597 + 6755399441055744ULL;
uint64_t iax = tbl[di.i & 2047];
double t = (di.d - 6755399441055744ULL)
* 0.0003384507717577858 - x;
uint64_t u = ((di.i + 2095104) >> 11) << 52;
double y = (3.0000000027955394 - t) * (t * t)
* 0.16666666685227835 - t + 1;
di.i = u | iax;
return y * di.d;
}

ベンチマーク
3種類用意した
単純和
配列版
SSE版（配列版fmath_expdをSSEで実装したもの）
fmath::expd_v
2011/8/6 /205
double sum = 0;
for (double x = 0; x < 1; x += 1e-8) {
sum += fmath::expd(x);
}
void fmath_expd(double *px, int n){
for (int i = 0;i < n; i++) {
px[i] = fmath::expd(px[i]);
}
}

単純和のベンチマーク
gcc 4.4.5 gcc 4.5.1 VC10 icc 12.0.5
std::exp 77.584 80.432 35.060 31.352
fmath::expd 12.680 13.067 12.842 11.582
2011/8/6 /206
 Core i7-2600 3.4GHz + 64-bit (Linux/Windows7)
 結果はどれも同一のsum=171828182.101727を出力
一つの関数あたりの処理クロック（clk）
gccは相変わらず遅い

配列版のベンチマーク(1/2)
誤差と速度の評価(gcc 4.4.5)
平均0, 標準偏差1の正規分布から10^6個生成した乱数
Taylor, Remezは岡崎直観さんの実装
cf. http://www.chokkan.org/blog/archives/352
expdの精度はRemez9次より良くて11次より悪い
とはいえdoubleの精度は約16桁なのでかなりよい
2011/8/6 /207
関数 clk 最大誤差二乗平均誤差
std::exp 67.17 0 0
fmath::expd 11.70 6.575812e-16 1.117645e-16
Taylor 11th 36.97 8.792464e-15 1.372958e-15
Remez 9th 16.29 1.909579e-14 9.909547e-15
Remez 11th 19.35 2.571683e-16 6.146323e-17
fmath::expd_v 8.27 6.575812e-16 1.117645e-16
配列版
（速度不利）
SSE版
（速度有利）

配列版のベンチマーク(2/2)
コンパイラによる差@Core i7-2600 3.4GHz(64bit)
結論
計算主体なら，とりあえずicc使っとけ
2011/8/6 /208
関数 gcc 4.4.5 VC10 icc 12.0.5
std::exp 67.17 35.45 11.83
fmath::expd 11.70 14.89 13.43
Taylor 11th 36.97 42.40 34.73
Remez 9th 16.29 15.34 15.63
Remez 11th 19.35 20.91 18.93
fmath::expd_v 8.27 7.20 8.24
iccすげー！
内部でベクトル版expを利用

exp(𝑥)の性質
和のexpはそれぞれのexpの積
exp 𝑥 + 𝑦 = exp 𝑥 exp(𝑦)
テイラー展開
exp 𝑥 = 1 + 𝑥 +
𝑥2
2
+
𝑥3
6
+
𝑥4
24
+ …
戦略その1
入力値𝑥を，テーブル引きするための大雑把な値𝑠とテイ
ラー展開するための小さい値𝑡にわける
𝑥 = 𝑠 + 𝑡
exp 𝑥 = exp 𝑠 exp(𝑡)
2011/8/6 /209
テーブル引きテイラー展開

テイラー展開の項数
どこまで計算するか?
1 + 𝑥, 1 +
𝑥2
2
, 1 +
𝑥2
2
+
𝑥3
6
, 1 +
𝑥2
2
+
𝑥3
6
+
𝑥4
24
, …
2次ならaddが1回, mulが2回
3次はaddが3回, mulが3回
4次はaddが4回, mulが4回, 5次は...
精度からの要請
doubleの精度は53bit(1e-16のオーダー)
𝑥が12bitなら𝑥4/24は12*4+4=52なので3次でよさそう
4次までするなら𝑥は9bitでよい(𝑥5/120→9*5+7=52)
今回は3次を選択
2011/8/6 /2010

𝑥を整数と小数部分にわける
𝑥 = round(𝑥) + (𝑥 – round(𝑥))
round(𝑥) = 𝑥に一番近い整数
𝑠 = round(𝑥)とおけば 𝑡 = 𝑠 − 𝑥で 𝑡 ≤ 1/2
 𝑡 ≤ 1/212
にするには
𝑥の代わりにx' = x * 2048を考える
𝑥′ = 𝑠′ + 𝑡′とすると|𝑡′| ≤ 1/2で
𝑥 = round(𝑥 ∗ 2048) / 2048 + 𝑡′/2048
2011/8/6 /2011
input : x
x' = x * 2048
s' = round(x')
s = s' / 2048
t = x – s // このとき x = s + tで|t|<=1/4096

テーブル引き
𝑥 = exp 𝑠 exp 𝑡 の exp 𝑠 の部分
𝑠は整数𝑛 / 2048の形で𝑛はどれぐらい?
exp(𝑥) = 0となる𝑥はlog(DBL_MIN)=-708.396
exp(𝑥) = Infとなる𝑥はlog(DBL_MAX) = 709.782
とするとテーブルエントリサイズは
(709 + 708) ∗ 2048 = 2.9 ∗ 106
でかすぎ
2011/8/6 /2012

敗因を考える
exp(𝑛/2048), 𝑛は整数を求めるのは難しかった
2の巾ならシフトと組み合わせれば簡単なのに
無理やり2の巾に持っていこう
𝑠と𝑡への分割方法を見直すと2巾である必然は無い
𝑥や𝑡が整数ではないし
exp 𝑛/′2000に近い値′ = 2^(𝑛/2^整数 )
の形にもっていけないか
exp(𝑛/𝛼) = 2^(𝑛/2 𝛽
) → 2 𝛽
= 𝛼 log 2
𝛼 = 2048のとき𝛼 log 2 = 1419.5
ということは𝛼 = 2048/(log 2)にすれば𝛽 = 12
2011/8/6 /2013

テーブル引き再度
𝛼 = 2048/log2とすると
exp 𝑛 ⁄ 𝛼 = 2^(𝑛/2048)
こうすれば
𝑛 = 2048(𝑛/2048) + (𝑛 % 2048) = 2048𝑛1 + 𝑛2
として
exp(𝑛/𝛼) = 2^(𝑛1)2^(𝑛2/2048))
テーブルエントリは2.9*10^6から2048に減った
これなら使える
2011/8/6 /2014

ここまでのおさらい
𝑥をα倍して整数部と小数部に分けて整数部は表
引き，小数部は3次までのテイラー展開
残りはround(𝑥)の処理方法
2011/8/6 /2015
input : x
x' = x * α; // α = 2048 / log 2 = 2954.639...
s' = round(x')
s = s' * (1/α);
t = x – s;
// exp(x) = exp(s) exp(t)

浮動小数点数を整数に丸める
FPU命令を使う
fistp
コントロールレジスタに依存（場合によっては変更の必要あり）
SSE命令を使う
cvtsd2si(double→int), cvtpd2dq（doublex2→intx2）
コントロールレジスタに依存
fisttp, cvttsd2si, cvttpd2dq
切り捨て専用なので直近に丸めるのとは違う
SSE4.1命令を使う
roundpd, roundps
CPU判別が必要で今回は保留
2011/8/6 /2016

浮動小数点数のトリックを使う（1/2）
doubleのフォーマット
𝑥 = −1 𝑠
× 2 𝑒−1023
× (1 + ⁄𝑓 252
)
すごく大きい値を足すと仮数部が丸められる
たとえば12.25は
12.25 = 23
∗ 1.53125 = 23
∗ (1 + 17/25
)
252
を足すと
12.25 + 252
= 252
∗ 1 + 2.72. . 10−15
= 252
∗ (1 + ⁄12 252
)
2011/8/6 /2017
符号ビットs（1bit）指数部e（11bit）仮数部f（52bit）
仮数部が丸められて整数になってる

浮動小数点数のトリックを使う（2/2）
今回の𝑥は負の数にもなるので補正項が必要
252
ではなく252
+ 251
を足す
doulbeのビットパターンを取り出すにはunionを使う
strict aliasingルールを破らないように（厳密には?）
2011/8/6 /2018
union di {
double d;
uint64_t i;
};
uint64_t roundC(double x) {
di di;
di.d = x + (3LL << 51);
return di.i & mask(52); // 52bitマスク
}

最終コード
マジックナンバーの意味
2011/8/6 /2019
double expd(double x) {
di di;
di.d = x * (2048 / log(2)) + (3ULL << 51);
uint64_t iax = tbl[di.i & 2047];
double t = (di.d – (3ULL << 51)) * (log(2) / 2048) - x;
uint64_t u = ((di.i + 2095104) >> 11) << 52;
double y = (3.0 - t) * (t * t) * (1.0/6) - t + 1;
di.i = u | iax;
return y * di.d;
}
α倍 round処理
テイラー展開
もとに
戻して1/α倍
2巾の簡単な方
2巾の表引き

おまけ（少しでも誤差を減らす）
テイラー展開の係数
1 + 𝑥 + ⁄𝑥2 2 + ⁄𝑥3 6 = 1 + 𝑥 + ⁄1 6 𝑥2 × 3 + 𝑥 で
1/6の代わりに0.16666666685227835
3.0の代わりに3.0000000027955394を利用
誤差がほんの少し減る（RMSが3%減）
係数の決め方 : 区間[0, t]でのexp(x)との二乗平均誤差
𝐼 = න
0
𝑡
exp 𝑥 − 𝑎 + 𝑏𝑥 + 𝑐𝑥2
+ 𝑑𝑥3 2
𝑑𝑥
が最小になるように

𝜕𝐼
𝜕𝑎
=
𝜕𝐼
𝜕𝑏
=
𝜕𝐼
𝜕𝑐
=
𝜕𝐼
𝜕𝑑
= 0を解き, 𝑡 = log(2)/4096を代入
2011/8/6 /2020

高速な倍精度指数関数expの実装

高速な倍精度指数関数expの実装

More Related Content

What's hot

Similar to 高速な倍精度指数関数expの実装

More from MITSUNARI Shigeo

高速な倍精度指数関数expの実装