高速な疑似六倍精度・疑似八倍精度計算法

高速な疑似六倍精度・疑似八倍精度計算法
山中脩也
(早稲田大学)
大石進一
(早稲田大学)
多倍長精度計算フォーラム@工学院大学
2014 年 3 月 7 日
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
はじめに
近年の高精度計算・多倍長計算の発展は凄まじい
ハードの変更が必要なく，現存の資産で動作
.
.
様々な関数が用意され，非常に便利
任意精度計算が容易
いろいろなソフトウェア：
apfloat, ARPREC, bbnum library, CLN, exflib, FMLIB, GMP,
HugeCalc, mpmath.... and so on.
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
はじめに
（一般的な）高精度計算・多倍長計算の現状
浮動小数点数形式を構成
仮数部 × 基数
指数部
事前に仮数部の桁数を決めて実行
倍々精度の仮数部：2 進 106bit
仮数部や指数部を長く保持するための方法
浮動小数点数を組み合わせて実装（ARPREC, DD/QD など）
整数を組み合わせて実装（exflib，GMP/MPFR など）
.
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
はじめに
高精度計算・多倍長計算の現状
浮動小数点数を組み合わせた実装（ARPREC, DD/QD など）
指数部と仮数部を分けて考える必要がない
比較的低い精度での多倍長で高速
オーバーフロー・アンダーフローの制約を受ける
桁数を増やすと，実行時間の増加が激しい
整数を組み合わせた実装（exflib，GMP/MPFR など）
仮数部をそのままビットで区切ればよく比較的分かりやすい
指数部・仮数部をそれぞれ多倍長化する必要がある
オーバーヘッドが大きい
桁数を増やしてもそこまで実行時間が変化しない
.
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
実行時間比較
乱数同士各 1 億回の計算（単位は秒，括弧内は sec per bit）
Mac OS 10.8, Intel Core i5, Memory 4GB, gcc 4.7
掛け算
DD/QD (浮動)
exflib (整数)
高速な疑似六倍精度・疑似八倍精度計算法
106 bit
0.58 (5.4e-3)
4.05 (3.8e-2)
159 bit
1.63 (1.0e-2)
5.05 (3.1e-2)
212 bit
4.59 (2.0e-2)
6.30 (2.9e-2)
. 山中脩也 (早稲田大学)
はじめに
本講演での焦点
高々 200 bit （≈ 10 進 64 桁）程度の精度を想定
浮動小数点数を利用した高精度計算法を使用
高速性に注目
精度保証付き数値計算もしたい！
八倍精度の高速計算フォーマットを提案
高速な疑似六倍精度・疑似八倍精度計算法
.
.
山中脩也 (早稲田大学)
200bit の応用例
CR-LIBM
A library of correctly rounded elementary functions in
double-precision
.
“99% correct rounding or 0.501 ulp accuracy.”
“The parameters are set up so that all the operators offer a
. require
relative error better than 2−208 . The worst cases never
more than 158 bits of accuracy.”
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
倍精度のみで高精度化
x1 , x2 ∈ F とする．このとき DD（Doubled-Double）の数 x は
x ' x1 + x2
(1)
x1 = fl (x1 + x2 )
(2)
と書け，
|x2 | 5 2−53 |x1 |
(3)
が成立する．
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
Key となるテクニック
無誤差変換
.
.
Knuth（1969）
.
[a, b] = TwoSum (x, y)
Dekker（1971）
[a, b] = TwoProduct (x, y)
Theorem
F : 浮動小数点数の集合，◦ ∈ +, −, ×,
x ◦ y = a + b,
.
∀x, ∀y ∈ F ⇒ ∃a, ∃b ∈ F
|b| 5 2−53 |a|
.
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
Key となるテクニック
[a, b] = TwoSum (x, y)
a = fl (x + y)
c = fl (a − x)
b = fl ((x − (a − c)) + (y − c))
function
end
function
.
.
function ( [xh , xt ] = Split
) (x)
c = fl (227 + 1) · x
xh = fl (c − (c − x))
xt = fl (x − xh )
.
無誤差変換
end
[a, b] = TwoProduct(x, y)
a = fl (x · y)
[x1 , x2 ] = Split(x)
[y1 , y2 ] = Split(y)
b = fl (x2 · y2 − (((a − x1 · y1 ) − x2 · y1 ) − x1 · y2 ))
end
高速な疑似六倍精度・疑似八倍精度計算法
. 山中脩也 (早稲田大学)
Key となるテクニック
無誤差変換
.
.
Knuth（1969）
.
[x, y] = TwoSum (a, b)
Dekker（1971）
[x, y] = TwoProduct (a, b)
.
擬似的な多倍長化
⇒ 倍精度だけを使って擬似的な多倍長を達成できる
誤差見積もりの簡便化
⇒ 誤差を明示的に知れ，精度保証付き数値計算に役立つ
.
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
QD とは
QD は 4 つの浮動小数点数を並べたもの
a = a0 + a1 + a2 + a3
QD には次の関係が成立する
a0 = nearest53 (a)
a1 = nearest53 (a − a0 )
a2 = nearest53 (a − a0 − a1 )
a3 = nearest53 (a − a0 − a1 − a2 )
QD の表現方法は唯一である． [Hida ら]
ただし，広く利用されている sloppy アルゴリズムは，
.
必ずしもこれを満たさない．
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
QD はもっと高速化できないか？（掛け算の例）
x × y ≈ x0 y0
+ x0 y1 + x1 y0
+ x0 y2 + x1 y1 + x2 y0
+ x0 y3 + x1 y2 + x2 y1 + x3 y0
= z0 + b0
+ b1 + c0
+ b2 + c1
+ c2 + d0
+ c3 + d1
+ c4 + d2
+ x0 y3
+ x1 y2
+ x2 y1
+ x3 y0
高速な疑似六倍精度・疑似八倍精度計算法
O(1) term
O(ε) terms
O(ε2 ) terms
O(ε3 ) terms
O(1) term
O(ε) terms
O(ε) terms
O(ε2 ) terms
O(ε2 ) terms
O(ε2 ) terms
O(ε3 ) terms
O(ε3 ) terms
O(ε3 ) terms
O(ε3 ) terms
山中脩也 (早稲田大学)
QD はもっと高速化できないか？（掛け算の例）
[a0 , a1 , a2 , a3 ] = multiplication(x0 , x1 , x2 , x3 , y0 , y1 , y2 , y3 )
[z0 , b0 ] = TwoProduct (x0 , y0 )
[b1 , c0 ] = TwoProduct (x0 , y1 )
[b2 , c1 ] = TwoProduct (x1 , y0 )
[z1 , z2 ] = ThreeSum (b0 , b1 , b2 )
[c2 , d0 ] = TwoProduct (x0 , y2 )
[c3 , d1 ] = TwoProduct (x1 , y1 )
[c4 , d2 ] = TwoProduct (x2 , y0 )
[z2 , z3 ] = SixSum (z2 , c0 , c1 , c2 , c3 , c4 )
z3 = fl (z3 + d0 + d1 + d2 + x0 · y3 + x1 · y2 + x2 · y1 + x3 · y0 )
[a0 , a1 , a2 , a3 ] = Renormalize (z0 , z1 , z2 , z3 )
function
end
.
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
提案フォーマット
4 つの浮動小数点数を並べたもの
a = a0 + a1 + a2 + a3
次の関係が成立する
( faithfuln (a) は 1 番目と 2 番目に近い n bit 浮動小数点数)
a0 ∈ faithful50 (a)
a1 ∈ faithful100 (a) − a0
a2 ∈ faithful150 (a) − a0 − a1
|a3 | < 2ε350 · ufp (a0 )
このとき次が成立する
|ak | < 2εk50 ufp (a0 )
表現方法はいくつかある（情報ロスはない）
高速な疑似六倍精度・疑似八倍精度計算法
.
山中脩也 (早稲田大学)
提案フォーマット
縦の足し算を，浮動小数点加算で計算が行なえるようにする
縦の足し算を，誤差の発生がないように定める
従来の QD では発生していた過剰な高精度を抑制し，高速化
表現方法はいくつかあるが，情報ロスはない
⇒ 53 bit から 50bit へ
.
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
提案フォーマットでの高速化
[a0 , a1 , a2 , a3 ] = multiplication(x0 , x1 , x2 , x3 , y0 , y1 , y2 , y3 )
[z0 , b0 ] = TwoProduct50 (x0 , y0 )
[b1 , c0 ] = TwoProduct50 (x0 , y1 )
[b2 , c1 ] = TwoProduct50 (x1 , y0 )
z1 = fl (b0 + b1 + b2 )
[c2 , d0 ] = TwoProduct50 (x0 , y2 )
[c3 , d1 ] = TwoProduct50 (x1 , y1 )
[c4 , d2 ] = TwoProduct50 (x2 , y0 )
z2 = fl (c0 + c1 + c2 + c3 + c4 )
z3 = fl (d0 + d1 + d2 + x0 · y3 + x1 · y2 + x2 · y1 + x3 · y0 )
[a0 , a1 , a2 , a3 ] = Renormalize (z0 , z1 , z2 , z3 )
end
.
function
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
提案フォーマットの新しい点
無誤差変換に加え任意の場所でビット列を切る関数を多用
function [a, b] = cut(x)
a = const[n] − (const[n] − x);
b = x − a;
end
.
function( [a, b] )= Split(x)
p = 227 + 1 ∗ x;
a = p − (p − x);
b = x − a;
end
高速な疑似六倍精度・疑似八倍精度計算法
.
山中脩也 (早稲田大学)
実行時間比較
掛け算
提案手法
DD/QD
exflib
203 bit
2.00 (9.8e-3)
4.54 (2.2e-2)
212 bit
3.81 (1.7e-2)
4.75 (2.2e-2)
.
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)
まとめ
高速八倍精度計算フォーマットの提案
従来の QD より高速な八倍精度フォーマットを提案した
.
ただし精度は QD より若干落ちる
.
これにより桁数を上げた際の実行時間の低下を抑制
高速化を達成
高速な疑似六倍精度・疑似八倍精度計算法
山中脩也 (早稲田大学)

Download Report