高速乗算法の設計と実装(1)

梅谷 武
作成:2001-08-25 更新:2005-04-20
ショーンハーゲ-ストラッセン法を俯瞰することによって、多倍長整数の高速乗算法の構造についてより抽象的な視点で検討し、その応用として32bit算術演算器向け高速乗算法を設計し、Pentiumへの実装を試みる。
IMS:20010825001; NDC:412.1; keywords:離散Fourier変換,数論変換;
目  次
1. 高速乗算法の構造
2. 高速乗算法の設計
参考文献
1. 高速乗算法の構造
1.1 ショーンハーゲ-ストラッセン法の構造
 ショーンハーゲ-ストラッセン法の概略を少し改良した形で述べ、その全体構造を俯瞰する。
数のP進表現と積の評価
 K=2k(k>0), N=2Kとし、正の整数a,bが基数P=2Kにより次のようにP進表現されるものとする。
a
aK-1PK-1++a1P+a0,0≦ai<P
b
bK-1PK-1++b1P+b0,0≦bi<P
この場合0≦a,b≦2N-1であるが、この積について0≦ab≦2N-1となるようにNを定める。このときab2N-1を法として求めればよい。
多項式としての積とその係数の評価
 PK≡1(mod 2N-1)より、
cabK-1
Σ
r=0
(
Σ
s+t≡r(mod K)
asbt)Pr
(mod 2N-1)
となる。この係数cr=s+t≡r(mod K)asbtの大きさを評価すると、
0crK(2K-1)2K(22K+1)
となる。したがってcrK(22K+1)を法として求めればよい。
孫子の剰余定理による係数の計算
(K, 22K+1)=1であるから、孫子の剰余定理からcrを、c'r及びc''rが既知であるとして、連立合同式
cr
c'r(mod K)
cr
c''r(mod 22K+1)
を解くことによって求めることができる。
離散フーリエ変換による係数の計算
 c'rは直接計算し、c''rは、(K, 22K+1, 24)型の離散Fourier変換を利用して次のように計算する。
F(a)i
K-1
Σ
s=0
asζsi(mod 22K+1)
F(b)i
K-1
Σ
t=0
btζti(mod 22K+1)
c''r
K-1K-1
Σ
i=0
F(a)iF(b)iζ-ir(mod 22K+1)
桁上げ処理
 最後にcr, r=0,…,K-1に桁上げ処理を施すことで、c=abP進表現が得られる。
1.2 高速乗算法の構造
 ショーンハーゲ-ストラッセン法の構造からわかることは、高速乗算法が次のような乗算変換を組み合わせてできあがっているということである。(→[N5],[D2],[D3])
定義1.1 (乗算変換)
 可換環Rから可換環R'への準同型写像を
  1. ある元pが生成するイデアルにより剰余環R'=R/(p)と自然な準同型写像R → R'=R/(p)
  2. 何らかの同型定理に基づく同型写像R → R'
という2つの方法のどちらかで作成し、可換図式
R×R
R
R'×R'
R'
によってRの乗算をR'の乗算に帰着させることを乗算変換と呼ぶ。
 ショーンハーゲ-ストラッセン法を構成する乗算変換を列挙してみる。
乗算変換1
 自然な準同型写像
ZZ/(2N-1)Z
による乗算変換
Z×Z
Z
Z/(2N-1)Z×Z/(2N-1)Z
Z/(2N-1)Z
乗算変換2
 自然な準同型写像
Z[x](Z/K(2N-1)Z)[x]
による乗算変換
Z[x]×Z[x]
Z[x]
(Z/K(22K+1)Z)[x]×(Z/K(22K+1)Z)[x]
(Z/K(22K+1)Z)[x]
乗算変換3
 孫子の剰余定理による同型写像
Z/K(22K+1)Z
Z/KZ×Z/(22K+1)Z
による乗算変換
Z/K(22K+1)Z×Z/K(22K+1)Z
Z/K(22K+1)Z
Z/KZ×Z/(22K+1)Z×Z/KZ×Z/(22K+1)Z
Z/KZ×Z/(22K+1)Z
乗算変換4
 G={ζi|i=0,1,…K-1}としたとき離散フーリエ変換の構造定理(→[N3])による同型写像
(Z/(22K+1)Z)[x]/(xK-1)(Z/(22K+1)Z)[G]
による乗算変換
(Z/(22K+1)Z)[x]/(xK-1)×(Z/(22K+1)Z)[x]/(xK-1)
(Z/(22K+1)Z)[x]/(xK-1)
(Z/(22K+1)Z)[G]×(Z/(22K+1)Z)[G]
(Z/(22K+1)Z)[G]
1.3 高速乗算法の設計法
 ショーンハーゲ-ストラッセン法が4つの乗算変換に分解することができるということをみたが、すべての高速乗算法はこのような乗算変換へ分解できることがわかっている。(→[N5])
 高速乗算法が乗算変換へ分解できるという性質を意識することによって、全体の見通しがよくなると同時に、部品となる乗算変換をいくつか用意して、それらをうまく組み合わせることによって高速乗算法を組み立てていくという設計法が可能になる。
 ショーンハーゲ-ストラッセン法を変形して、実際に高速乗算法を設計することについて考えてみる。まず多倍長整数を計算機の内部でどのように表現するかであるが、これは現在の計算機の多くが8=23bitの倍数の長さの2進算術演算器をもっているということから、P=2K, K=2k(k>0)としてP進表現するとしてよいであろう。
a
aN-1PN-1++a1P+a0,0≦ai<P
ションハーゲ-ストラッセン法においては長さNKにはN=2Kという関係があったが、これは理論上の計算量を最小にするための工夫であって、このままでは実用には適さない。実際にはこの関係をもっと自由度の高いものにしなければならない。
a
aN-1PN-1++a1P+a0,0≦ai<P
b
bN-1PN-1++b1P+b0,0≦bi<P
とするとき積ab2N-1を法として計算すると
cabN-1
Σ
r=0
(
Σ
s+t≡r(mod N)
asbt)Pr
(mod 2N-1)
となる。この係数cr=s+t≡r(mod N)asbtの大きさを評価すると、
0crN(2K-1)2
となる。この右側の評価式cr ≦ N(2K-1)2を変形して、右辺が孫子の剰余定理でうまく分解できるようにすることが設計の第1段階である。
 長さNのFFT算法を用意することが設計の第2段階である。このNの選択にあたっては離散フーリエ変換の条件を満たさなければならないという制約がある。この制約は第1段階とも関係し、また1の原始N乗根の値が計算しやすい形かどうかということも検討すべき項目である。
 次節以降、
  1. 適当な素数pを探して、係数の評価式としてcr < Npが成り立つようにする。
  2. N=2n, n>0であるとし、FFT算法としてCooley-Tukey型算法を使用する。
という方針で高速乗算法を設計し、そのPentiumへ実装法について検討する。
2. 高速乗算法の設計
2.1 数論変換
 我々はすでに可換環上で離散フーリエ変換を定義し、その諸性質をまとめている。(→[N3]) 工学分野では整数の剰余環上で定義された離散フーリエ変換のことを数論変換(NTT, Number Theoretic Transform)と呼ぶことがある。(→[S3],[D1])
 ここでは法mと長さnが与えられたとき、離散フーリエ変換が定義できるための条件について考える。長さの逆元が存在することは(n,m)=1と同値である。法mが素数のとき、1の原始n乗根の存在については次のことがいえる。
命題2.1
 pが素数のとき、
1の原始n乗根∈Z/pZn|p-1
証明 既約剰余類群(Z/pZ)*が位数p-1の巡回群になることからわかる。 ■
 法mが合成数のときも既約剰余類群(Z/mZ)*の構造を考えれば、1の原始n乗根の存在する条件がわかる。
 離散フーリエ変換の1の原始n乗根ζが存在する場合に条件3
n>k>0n-1
Σ
i=0
ζik=0
について考える。
(n-1
Σ
i=0
ζik)(1-ζk)=1-ζkn=0
が可換環であるという条件だけで成立するが、ここで1≠ζkであるから、すべてのkについて(1-ζk)が単元、すなわち
((1-ζk),m)1,n>k>0
のとき条件3が満たされることがわかる。法mが素数のとき、この条件はつねに満たされる。
 以下、素数pを法とする場合についてのみ考える。この場合剰余環Z/pZは有限体Fpとなる。この単元群(Fp)*は位数p-1の巡回群であり、その生成元rpを法とする原始根と呼ばれる。n|p-1なることと1の原始n乗根ζが存在することは同値であり、
ζ=r(p-1)/n
と表現される。原始根は一意的には定まらないのでこの表現も一意的ではない。上で示したようにFpにおいては1の原始n乗根ζが存在すれば(n,p,ζ)は離散フーリエ変換の条件を満たす。これを定理の形で整理しておく。
定理2.2 (有限体上の離散フーリエ変換)
 素数pと自然数nについて、
n|p-1∃ζ∈Fp:(n,Fp,ζ)は離散フーリエ変換の型
2.2 数列T(n)
 係数の評価式cr < Npにおける素数pを探す。
crN(2K-1)2=N(22K-2K+1+1)
が成り立っているので、
22K-2K+1+1p
なる素数pであって、Fpの演算、すなわち法pの剰余演算が高速に計算できるものを探さなければならない。我々はションハーゲ-ストラッセン法の実装実験(→[N4])において、数論変換を使う高速乗算法は、剰余演算が高速実行できなければ実用にはならないということをすでに経験している。
T(n)=22t-2t+1,t=2n,n0
とおくと
22K-2K+1+1T(k),K=2k
が成り立つ。T(n),n=0,1,2,3,4,5,6を素数判定してみる。
T(0)
=
2(2-1)+1=3:素数
T(1)
=
22(22-1)+1=13:素数
T(2)
=
24(24-1)+1=241:素数
T(3)
=
28(28-1)+1=65281=97*673:合成数
T(4)
=
216(216-1)+1=4294901761=193*22253377:合成数
T(5)
=
232(232-1)+1=18446744069414584321:素数
T(6)
=
264(264-1)+1:合成数
T(5)、すなわちK=25=32の場合に素数になっているが、これは32bitアーキテクチャのCPUには大変都合がよい。各係数は
A264+B,0≦A,B≦264-1
と表現できるが、このとき、
A264+B
=
(264-232+1)A+(232-1)A+B
(232-1)A+B(mod T(5))
より除算を使わずに剰余演算ができる。(→[S5])
2.3 離散フーリエ変換の型
 p = T(5) = 232(232-1)+1とする。これは素数であり、
crNp
という係数の評価式が成り立つ。さらにFp上に長さNの離散フーリエ変換が存在するための必要十分条件は、
N|232(232-1)
であるが、
N=2n,0<n≦32
はこの条件を満たし、さらにこのときpの最小原始根7を使ってFp上の1の原始N乗根を
7t,t=232-n(232-1),0<n≦32
と表現することができる。
参考文献
計算数論
[N1] 和田 秀男, コンピュータと素因子分解, 遊星社, 1999
[N2] 和田 秀男, 高速乗算法と素数判定法, 上智大学数学教室, 1983
[N3] 梅谷 武, 離散Fourier変換
[N4] 梅谷 武, ストラッセン-ショーンハーゲ法
[N5] Daniel J. Bernstein, Multidigit multiplication for mathematicians
算法
[S1] 野下 浩平, 高岡 忠雄, 町田 元, 基本的算法, 岩波書店, 1983
[S2] D. E. Knuth(中川 圭介訳), 準数値算法―算術演算, サイエンス社, 1986
[S3] H. J. Nussbaumer(佐川雅彦他訳), 高速フーリエ変換のアルゴリズム, 科学技術出版社, 1989
[S4] David H. Bailey, The computation of π to 29,360,000 decimal digits using Borweins' quartically convergent algorithm
[S5] Mikko Tommila, Number theoretic transforms
ディジタル信号処理
[D1] 電子情報通信学会, ディジタル信号処理の基礎, コロナ社, 1988
[D2] G. A. Jullien, Number theoretic techniques in digital signal processing
[D3] G. A. Jullien, Residue arithmetic with application in digital signal processing