FMT(高速剰余変換)の応用
2002/12/20 後 保範 ( Ushiro Yasunori )
2001/06/09 作成のものをweb化
2003/05/20 更新(自然数巡回FMT及び分割乗算等を追加) 後 保範
1. 概要と目的
FMT原理 |
で示した方法の多数桁演算への応用を述べる。
多数桁係数を持つ行列計算(行列乗算、連立一次方程式、固有値計算)や暗号解読への応用
はここでは述べない。
FMTを利用したソースプログラムは多数桁計算プログラム | 、
π計算プログラム | を参照。
2. 整数FMTと複数個のFMTの重ね合せ
FFTは原始n乗根
として複素数ω=e-2πi/nを使用したものといえる。一方、FMTは
FFTと同様に複素数でも定義できるが、整数でも可能である。ここでは、説明を簡単にするため
FMTの要素数nは偶数とする。
2.1 概要
FMTは1の原始n乗根となるωの上で定義される。そのため、mod(ωn,P)=1で、
mod(ωk),P)≠1,
k=1,2,…,n-1となる整数n,ω,Pでも可能となる。
ω=2,P=ωn/2+1とするとωはmod(P)上で原始n乗根となることが知られている。
一般に、ωが2
以上の整数ならP=ωn/2+1上でωは原始n乗根となり整数FMTが構成できる。
また、適当な大きさの素数P上で原始n乗根となるωが多数存在する。下記に現在の計算機で
扱い易い素数Pとn及びωの値を示す。この場合、複数個のPの上での多数桁乗算を
百五減算
で合成し、より大きい値の上での多数桁乗算が構成できる。
2.2 整数ωとmod(ωn/2+1)を利用
nを偶数とし、ωが2以上の整数とすると、P=ωn/2+1上でωは1の原始n乗根となる。
これを利用すると、要素数nのFMTが構成できる。
nが大きくなるとPは指数的に増大し、一般に計算機で扱いづらいと言われているが、2段階FMT
による多数桁の乗算等で利用できる。
2.3 倍精度計算可能な原始根
Pが224の近くの原始n乗根となる整数n,ω,Pの例を表1に示す。
表1. Pが224の近くのFFTは原始n乗根となる整数n,ω,Pの例
番号 |
n = α = 217 |
n = β = 3・215 |
n = γ = 32・214 |
NO. |
ω |
P |
ω |
P |
ω |
P |
1 |
770 |
149 α + 1 |
26 |
173 β + 1 |
175 |
114 γ + 1 |
2 |
201 |
153 α + 1 |
866 |
187 β + 1 |
434 |
132 γ + 1 |
3 |
70 |
155 α + 1 |
562 |
194 β + 1 |
12 |
136 γ + 1 |
4 |
20 |
164 α + 1 |
637 |
204 β + 1 |
604 |
147 γ + 1 |
5 |
1498 |
165 α + 1 |
100 |
220 β + 1 |
459 |
150 γ + 1 |
2.4 倍精度2個で計算可能な原始根
Pが248の近くの原始n乗根となる整数n,ω,Pの例を表2に示す。
表2. Pが248の近くのFFTは原始n乗根
となる整数n,ω,Pの例
番号 |
n = α = 235 |
n = β = 3・232 |
n = γ = 32・231 |
NO. |
ω |
P |
ω |
P |
ω |
P |
1 |
360 |
8319 α + 1 |
9375 |
22100 β + 1 |
8426 |
14636 γ + 1 |
2 |
933 |
8334 α + 1 |
5549 |
22226 β + 1 |
226 |
14771 γ + 1 |
3 |
1714 |
8549 α + 1 |
9520 |
22254 β + 1 |
4987 |
14896 γ + 1 |
4 |
474 |
8755 α + 1 |
3991 |
22314 β + 1 |
9607 |
14913 γ + 1 |
5 |
744 |
8759 α + 1 |
6955 |
22316 β + 1 |
7806 |
14994 γ + 1 |
2.5 複数個のFMTの重ね合せ
表1及び表2の複数個のPによる、FMTを用いた畳み込み演算結果は、複数個のPは互いに素な
なため、百五減算(中国剰余定理)で畳み込み演算結果の各要素の桁数を拡大することができる。
3. 正巡回FMTと負巡回FMTによる畳込み演算
3.1 概要
n個の要素の畳込み演算は、両方の入力にn個のゼロ要素を追加して、FMT変換又はFFT変換
を行い、互いの2n要素を対応する要素単位に乗算し逆変換することで2n個の結果が得られる。
もし、入力要素を2n桁に拡張しないで、そのままn桁のまま順変換、要素単位乗算及び逆変換
を行うと畳込み演算結果のn個の要素はどのようになるであろうか。
一般に行なわれているFFTで計算すると、2n要素の結果がn要素に後半のn要素が前半のn要素に
加算されて求まる。計算方法を変えると後半のn要素が前半の要素から減算されて求まる。
FFTの代りに、FMTを利用した場合は、前半と後半の要素がFFTと逆転する。
加算されて求まるとき正巡回FMT(FFT)、減算のとき負巡回FMT(FFT)と呼ぶ。
3.2 畳込み演算の方法
下記のように表現されるf,gの畳み込み演算結果hを正巡回と負巡回FMTで計算する方法を述べる。
f = an-1En-1 + …… + a1E + a0 ,
g = bn-1En-1 + …… + b1E + b0
--- (1)
h = f・g = c2n-1E2n-1 + c2n-2E2n-2 +
…… + c1E + c0 --- (2)
正巡回h+と負巡回h-のf・g畳込み演算結果を下記で表す。
h+ = mod (h , En - 1) = dn-1En-1 +
…… + d1E + d0
h- = mod (h , En + 1) = en-1En-1 +
…… + e1E + e0
--- (3)
このとき、dj, ejとcj, cj+nは(3)式の定義から
次の関係があることが分かる。
dj = cj + cj+n
ej = cj - cj+n , j = 0, 1, …… ,n-1
--- (4)
正巡回順FMTをFMT()で、正巡回逆FMTをFMT*()で表し、負巡回順FMTをFMTI()で、
負巡回逆FMTをFMTI*()で表すとh+, h-は下記のように
計算できる。
ここで、×印は要素別乗算を示す。
h+ = FMT*( FMT(f)×FMT(g) )
h- = FMTI*( FMTI(f)×FMTI(g) )
--- (5)
h+, h-の係数から求める2n個の畳込み演算結果hの係数cj
, cj+nは下記で求められる。
cj = ( dj + ej ) / 2
cj+n = ( dj - ej ) / 2 , j = 0, 1, …… ,n-1
--- (6)
3.3 負巡回FMTの計算方法
正巡回FMTはn要素のFMTそのままの計算であり、負巡回FMTの計算は下記のように正巡回FMT
に変更を加えて求める2種類の方法がある。
これは、En-1及びEn+1が下記のように原始n乗根ωを使用して分解できる
ためである。
( En - 1 ) = (E - ωn-1)(E - ωn-2) ……
(E - ω1)(E - ω0)
( En + 1 ) = (E - ωn-1/2)(E - ωn-3/2) ……
(E - ω3/2)(E - ω1/2)
(1) 原始n乗根ωのテーブルを変更する方法
下記の様に正巡回のωのテーブルを負巡回のωのテーブルに変更する。
逆変換を行うには、テーブル上のωをω-1に変更する。
ωのテーブルの与え方がFMTとFFTでは逆順になることに注意を要す。
正巡回のωのテーブル : ( ωn-1, ωn-2,
…… , ω1, ω0 )
負巡回のωのテーブル : ( ωn-1/2, ωn-3/2,
…… , ω3/2, ω1/2 )
(2) 入力係数及び出力係数に変換係数を掛ける
(1)式の入力係数ak,bkに(3)式の出力係数akに下記を
施して正巡回の計算を適用すると負巡回
FMTの結果が得られる。
入力係数変換: ak = akωk/2 ,
bk = bkωk/2
出力係数変換: ek = dkω-k/2 ,
k = 0,1,...,n-1
4. 2段階FMTを使用した畳込み演算
FMTやFFTを利用した時、メモリ的に問題となる多数桁の乗算を、2段階FMTを使用することにより、
FMTの計算効率を損なわずに、使用メモリ量を減少させる効果がある。
4.1 概要
1要素が多数桁を持つN要素の乗算を、P=ωN/2+1を法とする整数FMTを用いて実施する。
ここでは、Pを法とするFMTを上位FMT、要素ごとの項別乗算で必要たなるFMTを下位FMTと呼ぶ。
上位FMTは整数演算を使用し、ωkの乗算は桁上げ処理を除き整数の加減算だけで構成できる。
8バイト整数が利用できる場合はω=260(2進計算)又はω=1016
or 1018(10進計算)とし、4バイト整数
の場合にはω=228(2進計算)又はω=108 or 109に選んで
計算するのが、メモリ的にも速度的にも
効率がよい。
4.2 利用機能
2段階FMTでは、1要素が多数桁で構成される上位FMTと、その1要素の乗算に使用される下位FMT
がある。その各々に必要なFMTの機能を示す。
(1) 上位FMT
・ ωが2以上の整数で、P=ωN/2を法とする要素数Nの整数FMT
・ 正巡回FMTと負巡回FMTによる畳み込み演算
・ 整数8バイト又は整数4バイトの加減算処理(桁上げ処理を除く)で構成
・ 要素ごとの項別乗算は下位FMTを利用
・ ノード間並列化が必要な場合は上位FMTだけ並列化
・ スーパーコンでは1要素が百万桁に達する演算が可能
(2) 下位FMT
負巡回FMTを利用し、多数桁乗算結果が自動的にmod(P)を満たすようにする
・ 上位FMTによる多数桁乗算の要素単位に行う項別乗算に利用する
・ 実数FMTによる多数桁の乗算が現在のところ一番高速である
・ 計算は複素数16バイト又は実数8バイトを利用する
・ 整数乗算に実数値を使用するため、丸め誤差が発生し結果は整数化が必要
・ ノード間並列化ではこの部分の並列化は不要
4.3 計算方法
FMTを利用し多数桁演算の基底は2進数である必要はなく、10進数でも他の基底数でも可能で
あるが、説明を簡単にするため2進数を仮定して説明する。基底を変更するときは、ωの値を
基底数のべき乗にする必要がある。
FMTの順変換及び逆変換の計算方法はFMT原理 | を参照。
4.3.1 記号の定義
ここで使用する記号を下記のように定義する。
N : 上位FMTの要素数 n : 下位FMTの要素数
FP, FP-1 : 上位正巡回FMTの順変換及び逆変換
FN, FN-1 : 上位負巡回FMTの順変換及び逆変換
fN , fN-1 : 下位負巡回FMTの順変換及び逆変換
・ : 通常の乗算(畳込み演算による全要素乗算) × : 項別乗算
4.3.2 上位FMTの計算
上位FMTの要素数をNとしたとき、上位FMTの1要素に2進m・N桁を記憶できるようにする。
このとき、上位FMTはω=22mでP=ωN/2+1として構成する。
共にN個の要素を持つAとBの多数桁乗算を下記の様に表す。乗算結果は2N個の要素になる。
A ・ B = ( c2N-1 , c2N-2 ,……, cN ,……, c1
, c0 )
AとBの多数桁乗算を2N個の要素のFMTではなく、N個の要素の正巡回FMTによる乗算と、
負巡回FMTによる乗算を利用して下記の様に計算する。
正巡回乗算 : FP-1 ( FP ( A ) × FP
( B ) ) = ( dN-1 , dN-2 ,……, d1, d0 )
負巡回乗算 : FN-1 ( FN ( A ) × FN
( B ) ) = ( eN-1 , eN-2 ,……, e1, e0 )
AとBの多数桁乗算の2N個の要素ck+N, ckは下記で計算できる。
ck+N = ( dk - ek ) / 2
ck = ( dk + ek ) / 2
k = 0, 1,……, N-1
正巡回乗算と負巡回乗算で現われる要素単位の項別乗算(×)は下位FMTを利用して行う。
4.3.3 下位FMTの計算
上位FMTの要素単位に実施される項別乗算に下位FMTは利用される。そのため下位FMT
を利用する多数桁の乗算はmod(P)の上で実行される必要がある。
P=ωN/2+1=2mN+1のため、下位FMTの要素数をn 、下位要素の1要素に詰める
2進桁数をhとすると、
下記関係が必要となる。
h・n = m・N
この条件を満たせば、下位FMTに実数FMTを使用するか、整数FMTを使用するかは自由である。
要素数nでこの条件の下に、多数桁の乗算結果がmod(P)を満たすためには、負巡回FMTを下位
FMTに利用すればよい。即ち、aとbの2進h桁を1要素とするh・n桁の乗算は下記の様にn要素の
負巡回FMTを使用することにより、mod(P)のもとでh・n桁となり求まる。
mod( a・b , P ) = fN-1 ( fN ( a ) ×
fN ( b ) )
4.4 具体的計算サイズ
入力A,Bは0.5兆桁で出力Cが1兆桁(10進換算)となる3種類の多数桁乗算の使用メモリ量を
比較する。2段階FMT(実際は2段階FFT)は2002年の1兆2411億桁の
円周率世界記録 | 更新に
使用した多数桁乗算方式である。それまでの世界記録には実FMT(FFT)+Karatsuba法が使用
されていた。
(1) 2段階FMT(整数8バイトに2進60桁詰めで計算可能)
FFT(A)及び計算結果 : 0.5TB (FFT用ワーク、テーブルは不要)
FFT(B) : 0.5TB
正巡回 結果の保存 : 0.25TB
Work(下位FFT及び通信): 0.05TB
合計 : 1.3TB
(2) 実FMT(実数8バイトに2進9桁詰めが限度)
FFT(A)及び計算結果 : 3.0TB (1要素2進9桁、入力の12倍)
FFT(B) : 3.0TB
FFTワーク : 3.0TB (MATRIX/MPP)
合計 : 9.0TB
(3) 実FMT+Karatsuba法(Karatsuba法の適用回数でメモリ量と演算量が変わる)
1回カラツバ法を使用すると、FFT領域は半減し、Karatsuba法のワークが必要。
FFT ワーク 合計 演算量比
1回適用: 4.5TB + 1.0TB = 5.5TB 1.5 倍
2回適用: 2.3TB + 1.5TB = 3.8TB 2.3 倍
3回適用: 1.2TB + 1.7TB = 2.9TB 3.4 倍
4回適用: 0.6TB + 1.8TB = 2.4TB 5.1 倍
5. 複素FMTによる実数の畳込み演算
FMT原理 |
の「5. 畳込み演算とFMTとの関係式」で示した(10)式の応用でq=1/4のケースである。
FMT原理の(5),(6)式と(10),(11)式にq=1/4を代入した式を下記に示す。
f = an-1En-1 + …… + a1E + a0 ,
g = bn-1En-1 + …… + b1E + b0
h = f ・ g = c2n-1E2n-1 + c2n-2E2n-2
+ …… + c1E + a0
h(1/4) = mod ( h, En - ωn/4 ) =
c(1/4)n-1 E n-1 + …… +
c(1/4)1 E + c(1/4)0
c(1/4)j = cj + i ・ cj+n ,
j = 0, 1, …… , n-1
ここで、a j, b j, c j及びEを実数とし、ω = e 2πi/n
とすると、ωn/4 = i (虚数単位)となり、c(1/4)jは
cjを実数部に持ち、cj+nを虚数部にもつ複素数となる。
即ち、n要素(E進n桁)の実数乗算は、本複素
FMTを適用したn要素の複素数畳込み演算結果の虚数部に上位桁が、実数部に下位桁が求まる。
更に、下記の様にn/2-1次に畳込む、すなわちn/2要素のFMTにq=1/4を適用すると負巡回
多数桁乗算が直接得られる。
h(1/4,n/2) = mod ( h, En/2 - ωn/4 ) =
c(1/4,n/2)n/2-1 E n/2-1 + …… +
c(1/4,n/2)1 E + c(1/4,n/2)0
c(1/4,n/2)j = ( cj - cj+n )
+ i ・ ( cj+n/2 - cj+3n/2 ) ,
j = 0, 1, …… , n/2-1
これは、n要素(E進n桁)の実数乗算は、本複素FMTを適用したn/2要素の複素数畳込み演算結果
の虚数部に負巡回の上位桁が、実数部に負巡回の下位桁が求まる。
6. αで巡回する多数桁乗算
FMT原理 | の
(11)式を使用する。再度、この式を下記に示す。
c(q)j = cj + cj+n ω qn
c(q+n/2)j = cj - cj+n ω qn
q : 有理数 , j = 0, 1, …… , n-1
この式は、要素数nに対して法Pと自然数αに対してα=mod(ωqn,P)なる有理数q
と自然数P,ω,ωq
が存在すれば、+α及び-αで巡回する多数桁乗算となる。
表3にn=220で+-αで巡回するFMT乗算の係数を示す。
表3. n=220で+-αで巡回するFMTの係数の例
α |
P |
q |
ω |
ωq |
2 |
280,049,478・n+1 |
1/93,349,826 |
24,328,030,273,313 |
196,084,934,801,470 |
3 |
322,504,032・n+1 |
1/80,626,008 |
27,808,6974,322,725 |
45,280,667,942,168 |
4 |
283,546,518・n+1 |
1/94,515,506 |
237,787,923,135,329 |
113,109,652,731,064 |
5 |
205,984,108・n+1 |
1/205,984,108 |
43,678,380,610,349 |
192,726,477,692,165 |
6 |
282,686,616・n+1 |
1/282,686,616 |
20,159,304,809,371 |
288,630,981,189,601 |
7 |
222,921,552・n+1 |
1/111,460,776 |
161,161,411,978,886 |
163,047,718,620,128 |
8 |
28,0049,478・n+1 |
1/93,349,826 |
243,728,030,273,313 |
70,308,039,835,214 |
9 |
322,504,032・n+1 |
1/40,313,004 |
93,762,856,436,902 |
259,915,309,872,599 |
10 |
314,771,758・n+1 |
1/314,771,758 |
149,271,167,055,472 |
9,402,286,612,793 |
7. 分割畳込み演算
多数桁乗算のFMTを利用した分割乗算に付いて記述する。FMT原理の(11)式を応用する。
n要素(1要素はE進1桁)同士で結果が2n要素の多数桁乗算を、m=2l個に分割した乗算は
s=2n/m要素FMTを利用し下記の原理で行なう。qは有理数である。
h(q,s) = mod(h, Es - ωqn)
= mod(f ・ g, Es -ωqn)
= c(q,s)s-1 + …… +
c(q,s)1 + c(q,s)0
このとき、c(q,s)jと多数桁乗算の係数cjの関係式は次式の
様になる。
c(q,s)j = cj + cj+sωqn +
…… + cj+(m-1)s ω(m-1)qn
j=0,1, …… , s-1
FMTによる4分割乗算(m-4)の例を以下に示す。
m=4(s=n/2)でqを0,1/4,1/2及び3/4として上式に代入し整理すると得られる。
cj = [ ( c(0,n/2)j + c(1/2,n/2)j )
+ ( c(1/4,n/2)j + c(3/4,n/2)j ) ] / 4
cj+n/2 = [ ( c(0,n/2)j - c(1/2,n/2)j )
+ ( c(1/4,n/2)j - c(3/4,n/2)j ) ・
ω-n/4 ] / 4
cj+n = [ ( c(0,n/2)j + c(1/2,n/2)j )
- ( c(1/4,n/2)j + c(3/4,n/2)j ) ] / 4
cj+3n/2 = [ ( c(0,n/2)j - c(1/2,n/2)j )
- ( c(1/4,n/2)j - c(3/4,n/2)j ) ・
ω-n/4 ] / 4