高速乗算法の設計と実装(1)
梅谷 武
作成:2001-08-25 更新:2005-04-20
ショーンハーゲ-ストラッセン法を俯瞰することによって、多倍長整数の高速乗算法の構造についてより抽象的な視点で検討し、その応用として32bit算術演算器向け高速乗算法を設計し、Pentiumへの実装を試みる。
IMS:20010825001; NDC:412.1; keywords:離散Fourier変換,数論変換;
ショーンハーゲ-ストラッセン法の概略を少し改良した形で述べ、その全体構造を俯瞰する。
数のP進表現と積の評価
K=2k(k>0), N=2Kとし、正の整数
a,bが基数
P=2Kにより次のように
P進表現されるものとする。
この場合
0≦a,b≦2N-1であるが、この積について
0≦ab≦2N-1となるように
Nを定める。このとき
abは
2N-1を法として求めればよい。
多項式としての積とその係数の評価
PK≡1(mod 2N-1)より、
となる。この係数
cr=s+t≡r(mod K)asbtの大きさを評価すると、
となる。したがって
crは
K(22K+1)を法として求めればよい。
孫子の剰余定理による係数の計算
(K, 22K+1)=1であるから、孫子の剰余定理から
crを、
c'r及び
c''rが既知であるとして、連立合同式
を解くことによって求めることができる。
離散フーリエ変換による係数の計算
c'rは直接計算し、
c''rは、
(K, 22K+1, 24)型の離散Fourier変換を利用して次のように計算する。
| | |
| | |
| | K-1 | K-1 Σ i=0 | F(a)iF(b)iζ-ir | (mod 22K+1) |
|
桁上げ処理
最後にcr, r=0,…,K-1に桁上げ処理を施すことで、c=abのP進表現が得られる。
ショーンハーゲ-ストラッセン法の構造からわかることは、高速乗算法が次のような
乗算変換を組み合わせてできあがっているということである。(→[
N5],[
D2],[
D3])
定義1.1 (乗算変換)
可換環
Rから可換環
R'への準同型写像を
- ある元pが生成するイデアルにより剰余環R'=R/(p)と自然な準同型写像R → R'=R/(p)
- 何らかの同型定理に基づく同型写像R → R'
という2つの方法のどちらかで作成し、可換図式
によって
Rの乗算を
R'の乗算に帰着させることを
乗算変換と呼ぶ。
ショーンハーゲ-ストラッセン法を構成する乗算変換を列挙してみる。
乗算変換2
自然な準同型写像
による乗算変換
| | |
| | |
(Z/K(22K+1)Z)[x] | × | (Z/K(22K+1)Z)[x] |
| | |
乗算変換3
孫子の剰余定理による同型写像
による乗算変換
| | |
| | |
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)[G] | × | (Z/(22K+1)Z)[G] |
| | |
ショーンハーゲ-ストラッセン法が4つの乗算変換に分解することができるということをみたが、すべての高速乗算法はこのような乗算変換へ分解できることがわかっている。(→[
N5])
高速乗算法が乗算変換へ分解できるという性質を意識することによって、全体の見通しがよくなると同時に、部品となる乗算変換をいくつか用意して、それらをうまく組み合わせることによって高速乗算法を組み立てていくという設計法が可能になる。
ショーンハーゲ-ストラッセン法を変形して、実際に高速乗算法を設計することについて考えてみる。まず多倍長整数を計算機の内部でどのように表現するかであるが、これは現在の計算機の多くが
8=23bitの倍数の長さの2進算術演算器をもっているということから、
P=2K, K=2k(k>0)として
P進表現するとしてよいであろう。
ションハーゲ-ストラッセン法においては長さ
Nと
Kには
N=2Kという関係があったが、これは理論上の計算量を最小にするための工夫であって、このままでは実用には適さない。実際にはこの関係をもっと自由度の高いものにしなければならない。
とするとき積
abを
2N-1を法として計算すると
となる。この係数
cr=s+t≡r(mod N)asbtの大きさを評価すると、
となる。この右側の評価式
cr ≦ N(2K-1)2を変形して、右辺が孫子の剰余定理でうまく分解できるようにすることが設計の第1段階である。
長さNのFFT算法を用意することが設計の第2段階である。このNの選択にあたっては離散フーリエ変換の条件を満たさなければならないという制約がある。この制約は第1段階とも関係し、また1の原始N乗根の値が計算しやすい形かどうかということも検討すべき項目である。
次節以降、
- 適当な素数pを探して、係数の評価式としてcr < Npが成り立つようにする。
- N=2n, n>0であるとし、FFT算法としてCooley-Tukey型算法を使用する。
という方針で高速乗算法を設計し、そのPentiumへ実装法について検討する。
我々はすでに可換環上で離散フーリエ変換を定義し、その諸性質をまとめている。(→[
N3]) 工学分野では整数の剰余環上で定義された離散フーリエ変換のことを
数論変換(NTT, Number Theoretic Transform)と呼ぶことがある。(→[
S3],[
D1])
ここでは法mと長さnが与えられたとき、離散フーリエ変換が定義できるための条件について考える。長さの逆元が存在することは(n,m)=1と同値である。法mが素数のとき、1の原始n乗根の存在については次のことがいえる。
証明 既約剰余類群(Z/pZ)*が位数p-1の巡回群になることからわかる。
■
法mが合成数のときも既約剰余類群(Z/mZ)*の構造を考えれば、1の原始n乗根の存在する条件がわかる。
離散フーリエ変換の
1の原始
n乗根
ζが存在する場合に条件3
について考える。
( | n-1 Σ i=0 | ζik | )(1 | - | ζk) | = | 1 | - | ζkn | = | 0 |
が可換環であるという条件だけで成立するが、ここで
1≠ζkであるから、すべての
kについて
(1-ζk)が単元、すなわち
のとき条件3が満たされることがわかる。法
mが素数のとき、この条件はつねに満たされる。
以下、素数
pを法とする場合についてのみ考える。この場合剰余環
Z/pZは有限体
Fpとなる。この単元群
(Fp)*は位数
p-1の巡回群であり、その生成元
rは
pを法とする
原始根と呼ばれる。
n|p-1なることと
1の原始
n乗根
ζが存在することは同値であり、
と表現される。原始根は一意的には定まらないのでこの表現も一意的ではない。上で示したように
Fpにおいては
1の原始
n乗根
ζが存在すれば
(n,p,ζ)は離散フーリエ変換の条件を満たす。これを定理の形で整理しておく。
係数の評価式
cr < Npにおける素数
pを探す。
cr | ≦ | N(2K-1)2 | = | N(22K-2K+1+1) |
が成り立っているので、
なる素数
pであって、
Fpの演算、すなわち法
pの剰余演算が高速に計算できるものを探さなければならない。我々はションハーゲ-ストラッセン法の実装実験(→[
N4])において、数論変換を使う高速乗算法は、剰余演算が高速実行できなければ実用にはならないということをすでに経験している。
とおくと
が成り立つ。
T(n),n=0,1,2,3,4,5,6を素数判定してみる。
| | |
| | |
| | |
| | 28(28-1)+1 | = | 65281 | = | 97*673:合成数 |
|
| | 216(216-1)+1 | = | 4294901761 | = | 193*22253377:合成数 |
|
| | 232(232-1)+1 | = | 18446744069414584321:素数 |
|
| | |
T(5)、すなわち
K=25=32の場合に素数になっているが、これは32bitアーキテクチャのCPUには大変都合がよい。各係数は
と表現できるが、このとき、
より除算を使わずに剰余演算ができる。(→[
S5])
p = T(5) = 232(232-1)+1とする。これは素数であり、
という係数の評価式が成り立つ。さらに
Fp上に長さ
Nの離散フーリエ変換が存在するための必要十分条件は、
であるが、
はこの条件を満たし、さらにこのとき
pの最小原始根
7を使って
Fp上の
1の原始
N乗根を
と表現することができる。
計算数論
[
N2] 和田 秀男, 高速乗算法と素数判定法, 上智大学数学教室, 1983
算法
[
S1] 野下 浩平, 高岡 忠雄, 町田 元, 基本的算法, 岩波書店, 1983
ディジタル信号処理