FFTと多数桁乗算 No.2 多数桁乗算の各種方法 - 東京工芸大学

FFTと多数桁乗算 No.2
多数桁乗算の各種方法
2007年5月
後保範 (東京工芸大学)
1
目次
1. 筆算方式
2. Karatsuba法
3. 中国剰余定理（百五減算）
4. FFT(高速フーリエ変換)による乗算
4.1 整数FFTによる乗算
4.2 実数FFTによる乗算
4.3 複素FFTによる乗算
4.4 FFTとFMT(高速剰余変換)の関係
2
1. 筆算方式
35
73
26
87
×)
56
43
90
81
------------------------------------------------------------------------2835 5913 2106 7047
3150 6570 2340 7830
1505 3139 1118 3741
n2回乗算
+) 1960 4088 1456 4872
------------------------------------------------------------------------1960 5593 7745 15395 11994 9936 7047 ：桁上
桁上
20 16
72
00
15
94
06 47 ：結果
結果
3
2. karatsuba法
半分の桁の乗算を4回から3回に削減
通常FFT方式と合わせて使用される。
A = a1 E + a 2 ,
B = b1 E + b 2
C = A × B = c1 E
c 1 = a 1 × b1 ,
2
+ c2 E + c3
c3 = a 2 × b2
c 2 = a 1 × b 2 + a 2 × b1
d = ( a 1 + a 2 ) × ( b1 + b 2 ) = c 1 + c 2 + c 3
c 2 を c 2 = d − ( c 1 + c 3 ) で計算
4
3. 中国剰余定理(百五減算)
C = A×Bを下記で求める。
(1) m個の素数P1,P2,...,Pmを決める
(2) Ak=A mod(Pk), Bk=B mod(Pk)
(3) Ck=Ak×Bk mod(Pk)
(4) C=A×BをCkから中国剰余定理で求める
注) P1×P2×...×Pm > A×B
k=1,2,...,m
高速化のため筆算方式と組み合わせる
5
3.1 中国剰余定理の計算法
P1,P2,...,PmでのC1,C2,...,Cmから結果Cを復元
(1) C=C1とする。
(2) 下記をi=2からmまで繰り返す
w=(Ci - C)×ri mod(pi)
C = w×qi + C
注) 前もって下記を計算(i=2,3,...,m)
qi=P1P2 ...Pi-1 , ri=1/qi=qiPi-2 mod(Pi)
6
3.2 中国剰余定理の具体例(1/5)
• 素数Pi(i=1,2,3,4)とqi,ri(i=2,3,4)の値
P1= 1000003, P2 = 1000033
P3= 1000037, P4 = 1000039
q2 = 1000003, q3 = 1000036,000099
q4 = 1000073,001431,003663
r2 = 233341,
r3 = 551491
r4 = 682897
7
3.2 中国剰余定理の具体例(2/5)
• 入力(A,B)と乗算結果(C)
C=A×Bの計算
A = 123456,789012
B = 987654,321098
C = 121932,631136,585886,175176
8
3.2 中国剰余定理の具体例(3/5)
• Ci = Ai×Bi mod(Pi), i=1,2,3,4の計算
C1=418644×358145=805578 mod(1000003)
C2=715096×729605=400320 mod(1000033)
C3=221288×779269=498340 mod(1000037)
C4=974423×804113=644714 mod(1000039)
9
3.2 中国剰余定理の具体例(4/5)
• i=2の計算
w = (C2 - C1)×233341 = 813535 mod(P2)
C = w×q2 + C1= 813538,246183
• i=3の計算
w = (C3 - C)×551491= 730556 mod(P3)
C = w×q3 + C = 730583,113632,571227
10
3.2 中国剰余定理の具体例(5/5)
• i=4の計算(最終)
w = (C4 - C1)×682897 = 121923 mod(P2)
C = w×q4 + C
= 121932,631136,585886,175176
= A×B
11
4. FFT(高速フーリエ変換)による乗算
A=(a1,...,an), B=(b1,...,bn) : n個の要素
C=A×B, C=(c1,...,c2n) : 2n個の要素
(1) A=(a1,...,an,0,...0) : n個のゼロ追加
B=(b1,...,bn,0,...,0) : n個のゼロ追加
(2) A=FFT(A), B=FFT(B)
(3) C=A*B ; FFTの直交性で項別乗算
(4) C=FFT-1(C) ; C=A×Bとなる
12
4.1 整数FFTによる乗算
• 整数FFTによる乗算の特徴
(1) 計算方法は乗算原理そのまま
(2) 法Pの元での1の原始n乗根ωが必要
ω1/2 = -1 mod(P) , ωn = 1 mod(P)
(3) 下記条件要(hは入力要素の最大値)
n×h2 < P
(4) P1,P2,...,Pmを用いて重ね合わせが可能
13
4.1.1 整数FFTの係数例1
• 倍精度計算可能な1の原始n乗根
n
17
=1
mod(P)
,
n=2
の例
ω
NO
P(224程度)
ω
770
149・n + 1
1
2
201
153・n + 1
3
70
155・n + 1
4
20
164・n + 1
5
1498
165・n + 1
14
4.1.2 整数FFTの係数例2
• 倍精度計算可能な1の原始n乗根
n
35
=1
mod(P)
,
n=2
の例
ω
NO
P(248程度)
ω
360
8319・n + 1
1
2
933
8334・n + 1
3
1714
8549・n + 1
4
474
8755・n + 1
5
744
8759・n + 1
15
4.2 実数FFT方式
B=FFT(A), A=FFT -1(B) のとき
Aが実数 ==> Bは共役複素数を利用
A : 2n個の実数
B : 共役性を利用してn個の複素数
FFT, FFT-1 ; n要素の複素FFTと同等
項別乗算(*) ; 共役性でn個の複素乗算
==> 要素数が半分の複素FFTと同等
16
4.3 複素数FFT方式
•拡張FMT(入出力に変換を加える)から導出．
ωｊ/4(ｊ=0,...,n-1)の変換で下記の式が得られる．
Bk = Ak + i･Ak+n , k=0,1,...,n-1
Bk : 複素拡張FFT(FMT)による乗算結果
Ak : 実FFTによる乗算結果, i: 虚数単位
==> 入力値の虚部をゼロにした複素拡張FFT
の乗算の実部が元の多数桁乗算の上位
に、虚部が下位に対応する。
17
4.4 FFTとFMT(高速剰余変換)の関係
•FFT: 高速フーリエ変換 (フーリエが発見)
Fast Fourier Transformation ( ω = e 2π i / n )
信号処理や偏微分方程式解析用に開発
乗算結果の巡回が通常乗算と異なり上位から
•FMT: 高速剰余変換 (後が提案)
Fast Modulo Transformation (ωは原始n乗根)
多数桁乗算用に開発(変換中の意味明確)
FFTと多数桁乗算の演算量は同じ
剰余(mod)の定義に従い、巡回は下位から
18

Download Report