数論変換による多倍長整数乗算アルゴリズムの非再帰的実装の解析

DEIM Forum 2016 D8-5
数論変換による多倍長整数乗算アルゴリズムの非再帰的実装の解析
橋本 翔太†
上土井陽子†
若林 真一†
† 広島市立大学大学院情報科学研究科 〒 731–3194 広島市安佐南区大塚東三丁目 4 番 1 号
E-mail: †[email protected], {yoko,wakaba}@hiroshima-cu.ac.jp
あらまし
桁数の大きな整数のための高速乗算アルゴリズムに，高速フーリエ変換と畳み込み定理を用いたアルゴリ
ズムがある．このアルゴリズムでは，乗算の対象となるそれぞれの整数を定数桁ずつの成分に分割し，各定数桁成分
を時間成分とみなして高速フーリエ変換を適用した後，畳み込み定理により 2 つの整数の乗算結果を得る．ただし，
高速フーリエ変換は複素数体上での変換のため，数値計算に誤差が生じる場合がある．そこで，高精度な計算を行う
ため，高速フーリエ変換のかわりに有限体上での変換である高速数論変換と畳み込み定理を用いた乗算アルゴリズム
が知られている．本稿では，高速数論変換と畳み込み定理を用いた乗算アルゴリズムを非再帰的に実装し，数千桁の
整数の乗算に対して分割成分の桁数と計算効率の関係を解析する．
キーワード
多倍長整数乗算，数論変換，モンゴメリ乗算
1. は じ め に
2. 準
備
整数の乗算は，あらゆるアルゴリズムの中でよく用いられる
2. 1 多倍長整数の表現方法
基本的な算術命令の一つである．そのため乗算は，アルゴリズ
本稿で用いる多倍長整数を時間成分として表現する方法につ
ムに計算の基本的なステップとして遍在する．一般的なアルゴ
いて述べる．桁数が N の非常に大きな 10 進数の整数 A を定数
リズムの中で乗算の対象となる整数は，ハードウェアにより高
桁 sp で分割し，桁数 N ′ = N/sp，基数 B = 10sp として式 (1)
速に処理できる程度のビット長が小さな整数であるため，乗算
のように多項式表現する．ここで N は 2 のべき乗で表現可能
の計算量はアルゴリズムの全体効率にほとんど影響しない．し
な十分に大きな整数，sp は 2 のべき乗で表現できる数とする．
かし，ビット長の大きな整数同士の乗算を必要とする RSA 暗
号，El-Gamel 暗号，楕円曲線暗号のような公開鍵暗号システ
ムの出現によって，乗算の計算量の低減に関する研究は大きな
重要性をもつことになった．
上記のような暗号技術で用いる桁数の大きな整数に対する高
速乗算アルゴリズムとして，Fürer の乗算アルゴリズム [3] の
中でも高速フーリエ変換 (fast Fourier Transform, FFT) と畳
み込み定理に注目し，基本となっているアルゴリズムを実装し
A = a0 × B 0 + a1 × B 1 + · · · · · · + aN ′ × B N
′
−1
(1)
′
次に ai = 0 (N ′ <
= i < 2N ，i は整数) とし，式 (1) の係数
系列を式 (2) のように表現する．
( a0 , a1 , · · · , aN ′ −1 , 0, 0, · · · , 0)
(2)
さらに，係数系列の各要素を時間の関数 f (t) の関数値とす
ると，式 (3) のように表現できる．
てシミュレーションによる実験的性能評価を行った．その結果，
多倍長整数乗算を通常の筆算で行った場合に比べて，大幅に計
算時間の高速化ができていることがわかった．しかし，円周率
π や三角関数 sin, cos などの無理数の計算を必要とする複素数
体上のフーリエ変換では乗算結果に誤差が生じる場合があるこ
とがわかった．乗算結果を暗号技術で用いる場合には高速に計
( x(0), x(1), · · · , x(N ′ − 1), · · · , x(2N ′ − 1) )
式 (2)，式 (3) の間には，次のような関係が成り立っている．

at
x(t) =
0
算できることよりも，高精度な計算ができる方が重要度が高い
と考えられる．
本稿では，誤差が生じることがない有限体上での変換である
数論変換と同じ計算フレームワークである高速数論変換 (fast
Number Theoretic Transform, FNTT) と畳み込み定理に注目
し，筆算，高速フーリエ変換と畳み込み定理を用いた乗算アル
ゴリズムと計算時間，計算精度について実験的評価を行った．
また，高速数論変換と畳み込み定理を用いた乗算アルゴリズム
の非再帰的実装における課題とそれに対する解決法について考
察する．
(3)
′
(0 <
=t<N )
′
(N ′ <
= t < 2N )
(4)
本稿では，N 桁の整数を分割桁数 sp で分割し，0 を詰めた
あとの桁数を M (= 2N ′ ) とする．
例えば，N = 8 の 10 進数の整数 31589182 を分割桁数 sp = 2
のとき，桁数 N ′ = 8/2 = 4，基数 B = 102 = 100 として，
31589182 = 82 × 1000 + 91 × 1001 + 58 × 1002 + 31 × 1003
と表せるので，係数系列は
(82, 91, 58, 31, 0, 0, 0, 0 )
である．また，この係数系列を時間成分の系列とみなすと，図
45#!
45%!
()*+',-!
.#"'&/#"(&/000/#")*(&1!
()*+',-!
.%"'&/%"(&/000/%")*(&1!
,-!.!
&!"#$%&"'(/
)*+*!,-."+/
"/01234!
+#!
*$!
'*!
+!
%#!
"#$%&"+&!,!!
"! "! "! "!
"! #! $! %! &! '! (! )!
.!,!'!
23)!
!!
-!#".&$%"+*.&!
!"#$%&'()*+',-!
."#$%&"'&/"#$%&!"(&/000/"#$%&!")*(&1!
図 1 整数を波形で表現
6789:!
1 のように係数を波としても表現できる．
;<=>?#!$!%!
一般に桁数 N において，分割桁数 sp を大きくとると，整数
A を表す時間成分の系列数は減少し，それぞれの時間成分の値
図 3 畳み込み演算による乗算手法 (筆算)
は大きくなる．このことが，計算時間や計算精度にどのような
影響を与えるかについては後述することにする．
ることとほぼ同じであると述べた．本節では，畳み込み演算に
2. 2 畳み込み演算
畳み込み演算とは，式 (4) で示したような時間に関する関数
x(t) を平行移動しながら関数 y(t) を重ね足し合わせる二項演
算 (x ∗ y)(t) である．離散値で定義された関数 x(t), y(t) に対す
る畳み込みは，式 (5) で定義される．
比べて畳み込み値を高速に得られる可能性がある循環畳み込み
性質 (Cyclic Convolution Property, CCP) をもった変換であ
る FFT と畳み込み定理に注目する．FFT と畳み込み定理を用
いて畳み込み値を計算し，計算時間，計算精度について筆算で
計算した場合と比較を行う．
(x ∗ y)(t) =
t
∑
x(u) y(t − u)
(5)
3. 1 高速フーリエ変換とその逆変換
乗算を行う整数を図 1 のように波で表現し，その波に対し高
u=0
式 (5) は巡回畳み込みになるので，本研究で扱う整数の係数
系列では N ′ 桁の整数に対して桁数を 2 倍にし，上位の桁の部
分に 0 をつめて巡回しないようにし，畳み込み演算を利用する．
速フーリエ変換を適用し，実空間からフーリエ空間へ移動する．
またフーリエ変換は可逆変換であるため逆変換も存在し，フー
リエ空間から実空間に移動する場合は，高速逆フーリエ変換を
適用する．
2. 3 筆算と畳み込み演算との関係
通常の筆算による乗算を，8 桁の 10 進数の整数 31589182 と
整数 54177913 を分割桁数 sp = 2 とし，4 桁の 100 進数の整数
信号 x(n) に対するつぎの計算を，M 点離散フーリエ変換
(discrete Fourier transform, DFT) という．
としたときを例に図 2 に示す．図 2 から，乗算結果は畳み込み
X(k) =
M
−1
∑
x(n) e−j2πnk/M
(6)
n=0
!"!
#$!
%"!
'!"#(!
")!
)%!
$&!
"!! ./"01/*0!
(*"! )#(! ""$!! +"*,,!
&((%! (#$&! )"$%! ,()$!
#&)! %$,! "#()! "!%(!
()*)+!
./*01/"0!
/.210/"03./*01/"04./"01/*0!
")!
を必要とする．
本稿では，DFT に比べて少ない計算量で計算可能である
",)(! !"!&! (%"(! ((&$!
Cooley と Tukey によって開発された時間間引き FFT [1] を利
",)(+ !,#%+ $!(%+ "*%,*+ %!!)+ ),,"+ "*,,+
用する．FFT は，M = 2i (i は自然数) 点 DFT の計算量を
-#$%&'!
,-./!
一般に，離散フーリエ変換 (DFT) を用いて時間領域の信号
x(n) を周波数領域の信号 X(k) に変換するためには膨大な計算
""!
(!!
#%!
#(!
"!!
図2
筆算による乗算
)"!
,,!
O(M 2 ) から O(M log M ) に削減したものである．
式 (6) に対して，
値を桁上げ処理にすることで得られるということがわかる．桁
WM = e−j2π/M = cos
を利用すると，
数 M の畳み込み演算の計算量は O(M 2 ) であり，乗算結果を
X(k) =
得るための計算時間のほとんどは畳み込み値の計算 (畳み込み
演算) であり，畳み込み値を高速に計算することが可能になれ
ば，乗算結果を高速に得ることができると考えられる．
畳み込み演算を用いた乗算手法の流れを図 3 に示す．
3. 高速フーリエ変換に基づく多倍長整数乗算
2π
2π
− j sin
M
M
M
−1
∑
nk
x(n) WM
(7)
n=0
となる．ここで，式 (7) は回転因子とよばれ，図 4 のような
単位円を M 分割したものであり，周期性と対称性をもつ．図
4 は M = 8 のときの回転因子 W8 である．
以上の性質とバタフライ演算の繰り返し使用により計算量の
大幅な低減が可能となる．M = 8 のときの FFT の計算フロー
前節で乗算結果を高速に得ることは，畳み込み値を高速に得
を図 5 に示す．
また，次の計算を高速逆フーリエ変換 (inverse fast Fourier
transform, IFFT) という．
x(n) =
,-"!
M −1
1 ∑
−nk
X(n) WM
M
;<"!
;<.!
!"#$%&'!
("#$%)"#&%)***)"#'(&%+!
!"#$%&'!
(.#$%).#&%)***).#'(&%+!
(8)
))*!
k=0
))*!
(+#$%)+#&%)***)+#'(&%+!
./!
"!
1234
-"!
#$'!#$!%!
#$(!#$!!!
(!+#$%,#$%)+#&%,#&%)***)+#'(&%,#'(&%!+!
#$!!#$&!
/))*!
#$%!#$$!
01!
#$)!#$!'!
56789:%!"#$%&'!
(#"-.%#$%)#"-.%!#&%)***)#"-.%!#'(&%+!
!!
-!!
#$*!#$!(!
,-"!
=>?@A!
#$+!#$!*!
-"!
(,#$%),#&%)***),#'(&%+!
#$./%0!
BCDEF"!-!.!
#$,!#$!)!
図 6 FFT と畳み込み定理を用いた乗算手法
図 4 回転因子 W8
X(0)!
x(0)!
W8
x(2)!
0!
W8
0!
W8
X(1)!
W 84!
W 82!
W 81!
W8
x(6)!
0!
4!
W 82!
W86!
W83 !
W8
X(3)!
W 84!
W8
x(3)!
W8
0!
W8
4!
X(5)!
W 84!
x(5)!
x(7)!
0!
W8 2 !
W8
0!
W8
4!
W85 !
W8
W8 6 !
図5
ミュレーション実験を行った．実験で用いた整数は，10 進数で
みなし生成した整数である．本稿では，N = 2i (10 <
=i<
= 13)
となる各桁数 N の分割桁数 sp = 1, 2, 4, 8 に対して，異なる 5
つの整数の組 Data1∼Data5 に対して乗算結果を計算する実験
X(6)!
6!
X(7)!
W 84!
て，CPU:3.33GHz 6-Core Intel Xeon の計算機上に実現しシ
1 から 9 の数をランダムに N 個並べた数値を N 桁の 10 進数と
X(4)!
x(1)!
理を用いた乗算手法 (図 6) の 2 つの乗算手法を C 言語を用い
て実装し，gcc 4.2.1 コンパイラ，最適化オプション-O3 を用い
X(2)!
x(4)!
3. 4 実験的性能評価 (1)
畳み込み演算を用いた乗算手法 (図 3)，FFT と畳み込み定
0!
W87!
を行った．
3. 4. 1 計算時間について
上記の 2 つのアルゴリズムに，整数 x，y の組み合わせ Data1
∼Data5 に対して乗算時間の平均値を表 1 にまとめた．表中の
バタフライ演算
3. 2 畳み込み定理
N は計算に用いた 10 進整数の桁数，sp は整数の分割桁数を
関数 x(t), y(t) の畳み込み演算とフーリエ変換には，式 (9)
表す．
の関係が成り立つ．ただし，F [z(t)] は関数 z(t) をフーリエ変
表 1 から，FFT と畳み込み定理を用いた手法は桁数増加に
対して計算時間の増加は緩やかで，多倍長整数に対する計算に
換したものであるとする．
適していることがわかる．また，sp を増加させることにより計
F [(x ∗ y)(t)] = F[x(t)] · F[y(t)]
(9)
式 (9) から，関数 x(t), y(t) に対する畳み込みをフーリエ変換し
たものは，関数 x(t), y(t) をそれぞれフーリエ変換したもの同
士の積と等しいということが分かる．さらに，次の式も導ける．
[
]
(x ∗ y)(t) = F −1 F [x(t)] · F [y(t)]
算時間が約半分になっていることから，sp を増加させることに
よる高速化の可能性がある．
表 1 各手法の乗算平均計算時間 (表中の灰色の部分は結果に誤差有)
桁数
N
(10)
式 (10) から，畳み込みをフーリエ変換と逆フーリエ変換を使
計算時間 [ms]
(10 進)
畳み込み演算 (筆算)
FFT と畳み込み定理
sp = 1
sp = 2
sp = 4
sp = 8
sp = 1
sp = 2
sp = 4
sp = 8
1024
2.180
0.679
0.297
0.190
0.604
0.338
0.224
0.182
2048
7.783
2.186
0.739
0.364
1.125
0.619
0.384
0.287
4096
30.064
7.963
2.343
0.880
2.301
1.226
0.713
0.502
8192
117.201
30.322
8.225
2.610
4.651
2.404
1.409
0.946
うことで，高速に計算できる可能性がある．
3. 3 FFT と畳み込み定理を用いた乗算手法
FFT と畳み込み定理を用いた乗算手法の流れを図 6 に示す．
また，FFT は円周率や三角関数などの浮動小数点数を使用し
ているため，計算された畳み込み値を小数第一位を四捨五入し
て桁上げ処理を行う．
3. 4. 2 計算精度について
畳み込み演算による乗算結果が正しいとしたとき，実験を
行った全ての桁数 N に対して sp = 1, 2, 4 について FFT と
畳み込み定理を用いた手法で正しい結果が得られた．しかし，
sp = 8 のとき最悪の場合ですべての畳み込み値に誤差が生じて
おり，正しい乗算結果を得ることができなかった．理由として，
sp = 8 の場合，64 ビット浮動小数点数 double 型変数は 10 進
数で表すことのできる精度の桁数が 15(実験で用いたコンパイ
き，FFT と同様にバタフライ演算等により，高速に変換がで
ラのヘッダファイル float.h の DBL_DIG の値) であり，精度の
きるので以降，高速数論変換 (fast Number Theorem Trans-
保障された範囲を超えた計算が行われたからと考えられる．sp
form, FNTT) と呼ぶことにする．
の増加により高速化は期待できるが，計算精度が低下してしま
うことがわかる．
FNTT と畳み込み定理を用いて，四則演算が定義され閉じて
いる有限集合の中で整数を扱い，誤差が生じることもなく高速
4. 高精度な多倍長整数乗算手法の非再帰的実装
FFT は，実験結果より桁数 N や分割桁数 sp によって計算
に乗算を計算できる可能性がある．
4. 2 数論変換のためのパラメータ P ，α の導出
式 (14) の FNTT を行うためには，素数 P と位数 M の原
精度が低下する場合があった．そこで，複素数体上で変換して
始根 α の 2 つのパラメータを事前に計算しておく必要がある．
いた数を有限個の元からなる体，つまり有限体上で行うフーリ
まず，初等整数論の性質，用語についてまとめ，具体的なパラ
エ変換とみなされる高速数論変換を用いた乗算手法について考
メータの算出方法について説明する．
察する．また，標準ライブラリで用意されている演算の精度保
4. 2. 1 整数論の性質・用語 [5] [8]
証範囲を超えて分割桁数 sp を大きくすることによる高速化を
p は素数とする．一般に，
目標として，高速数論変換を用いた乗算手法の非再帰的実装を
ap ≡ a
試みる．
数論変換では合同式を多用する．整数 a と b に対して，a − b
が自然数 m で割り切れるとき，a と b は m を法として合同で
(mod m) ⇐⇒ m|(a − b)
割って，
ap−1 ≡ 1
(11)
(16)
となる．もし，a が p の倍数でなければ，式 (16) の両辺を a で
あるといい，記号で a ≡ b (mod m) と書き，つまり
a≡b
(mod p)
(mod p)
(17)
が得られる．これをフェルマーの小定理という．式 (17) にお
いて，p − 1 乗しなくても 1 と合同になることがある．例えば，
とする．
a = 3，p = 13 のとき
4. 1 数 論 変 換
a3 = 27 ≡ 1
式 (7) の FFT と同じ構造
X(k) =
M
−1
∑
x(n) αkn
(12)
n=0
で表現され，循環畳み込み性質を有する変換を考える．ここ
で，M が αM = 1 を満たす最小の整数となれば，式 (12) の
α=e
(18)
となる．ae = 1 (mod p) となる最小の自然数 e を， mod p で
の a の位数という．
a の位数が e のとき，ak ≡ 1 ならば k は e の倍数となる．な
ぜなら，k = q · e + r ，0 <
= r < e とすると，
1 ≡ ak = (ae )q · ar ≡ ar
変換は循環畳み込み性質をもつことが知られている．特に
−j2π/M
(mod 13)
(mod n)
(19)
= WM とした変換は前述した FFT である．
これに対して，整数の世界で，ある整数 P を法とする演算の
となり，0 < r なら e の最小性に反するからである．つまり位数
もとでは循環畳み込み性質を有する変換がいくつかあることが
e は，p − 1 の約数である．a の位数が p − 1 の約数である e と
知られている [6]．式 (12) の変換式で，
なるとき，a を mod n での位数 e の 1 の原始根という．また，
αM ≡ 1
(mod P )
(13)
を満たす最小の整数が M であるような α を 1 の原始 M 乗
j≡k
(mod e)
→
aj ≡ ak
(mod n)
(20)
である．
根，または，位数 M の原始根といい，この α を用いると循環
自然数 n に対して 1, 2, 3, · · · , n の中で n と互いに素な数の個
畳み込み性質をもつ．このようにして定義される変換を数論変
数を ϕ(n) で表すことにする．ϕ(n) はオイラー関数といわれて
換（Number Theoretic Transform, NTT）といい，式 (14) に
いる．ϕ(1) = 1 であり，素数 p に対しては，1, 2, 3, · · · , p の中
示す．
で p と互いに素な数は 1, 2, 3, · · · , p−1 であるから ϕ(p) = p−1
X(k) =
M
−1
∑
x(n) αkn
となる．その他の n についてのオイラー関数については本稿で
(mod P )
(14)
n=0
また，FFT と同様，数論変換も可逆変換であるため，次の
式 (15) のように逆数論変換 (Inverse Number Theorem Trans-
は扱わないので説明は省略する．
p を素数とするとき， mod p の位数 p − 1 原始根が存在し，そ
の個数は ϕ(p − 1) 個存在する (原始根存在定理)．さらに，一
般に p の大きさに関わらず，modp の位数 e の原始根の個数
form, INTT) を定義できる．
は ϕ(e) で表される．また，オイラー関数 ϕ(n) の値は高々n で
x(n) =
1
M
M
−1
∑
X(k) α−kn (modP )
(15)
k=0
上記の数論変換において，桁数 M を 2 のべきとしたと
あり，かつ平均して Θ(n) の大きさをもつことも重要な性質で
ある．
表2
乗算を行うためのパラメータ一覧
桁数
(10 進)
(位数 p − 1 の原始根を求める一般的な計算方法による算出)
乗算を行うためのパラメータ一覧
分割桁数 sp = 1
原始根 α
分割桁数 sp = 2
素数 P
分割桁数 sp = 4
原始根 α
素数 P
1024
176129
153 10039297
2048
331777
386 20078593
4096
737281
96 40148993
15170 204759054337
55931003
8192 1376257
73 80330753
46283 409518112769
24724426 64bit で表現できない
△
3201
素数 P
分割桁数 sp = 8
N
原始根 α
素数 P
原始根 α
2559999948800007169
△
5119999897600035841
△
10239999795200005121
△
51189761537 182252995
47486 102379527169 244863846
※ △ は 24 時間では計算できなかった場合
αM/2 ≡ P − 1
4. 2. 2 計 算 方 法
(mod P )
(22)
具体的な数論変換で用いるパラメータ，素数 P と位数 M の
を計算すればよい．
原始根 α の計算方法を説明する．
FNTT と畳み込み定理を用いて乗算する整数の桁数を N ，分
割桁数 sp とし，位数 M（ = 2N/sp）として数論変換を行いた
いとする．
まず，分割桁数 sp のとき，基数 B = 10sp となるので，桁数
表 2 に桁数 N と分割桁数 sp に対応した一般的なライブラリ
で実用的に算出できた素数 P と位数 M の原始根 α をまとめ
た．表中の △ がついているところは位数 M に対して素数 P
が非常に大きく，一般的な位数 p − 1 の原始根を求める方法と
同様な方法では 24 時間では計算できず，実用的な時間では計
M の畳み込み値の最大値 convM ax は，
算できない場合である．
convM ax = (B − 1) · (B − 1) · M
(21)
である．位数 M の原始根 α を求め，数論変換を用いてすべて
の計算を mod P 上で行うためには次の 2 つの条件を同時に満
たす素数 P を選択する必要がある．
条件 1 素数 P は畳み込み値の最大値 convM ax より大きい素
数である
4. 3 非再帰的実装における課題
本研究では数論変換による乗算手法を特殊なライブラリを用
いず，64bit の整数 long 型に対応する標準ライブラリ関数によ
り，図 6 で高速フーリエ変換により変換していた部分を高速数
論変換に置き換えた高精度な多倍長整数乗算手法を非再帰的実
装する． 非再帰的な実装のため，素数 P を大きくとらなけれ
条件 2 P − 1 が M の倍数である，つまり P − 1 = sM （s は
自然数）と表現できる
012"
UV2"
UV4"
92'(!:;$
#2*+,$2*-,$%%%$2*.'/0-,&"
92'(!:;$
#4*+,$4*-,$%%%$4*.'/0-,&"
条件 1 を満たさなければいけない理由は，convM ax よりも
小さい素数を P として採用すると真の値 x なのか P − x なの
!"#$%!&'(("
#)*+,$)*-,$%%%$)*.'/0-,&"
かが判別できなくなるためである．
条件 2 は，位数 M の原始根 α を求めるために必要な条件で
!"#$%,
-!./"
<=>?@A*BCDE$
FGHIJBKL"
M の 1 の原始根 α が存在することがわかる．しかし，実際に乗
!"#$%!5&'(("
345678!92'(!:;$
#*234,*+,$*234,$*-,$%%%$*234,$*.'/0-,&"
算を行うときにこのような素数 P と位数 M の差が 1 となる P
と M の組み合わせで計算することはない．本稿で用いる数論
#1*+,$1*-,$%%%$1*.'/0-,&"
'()*!+"
#$)*+,1*+,$)*-,1*-,$%%%$)*.'/0-,1*.'/0-,$&"
ある．フェルマーの小定理 (式 (17)) より，素数 P としたとき，
P − 1 乗して 1 になる．つまり M = P − 1 としたとき，位数
!"#$%!&'(("
012"
変換，特に乗算で用いる原始根は一般的な暗号計算で用いられ
M,NOG"
PQRST2$3$4$
る素数 P の位数 p − 1 の原始根（そのような原始根は ϕ(p − 1)
個ある）ではなく，位数 M ，つまり mod 上で M 乗して初め
て 1 となる原始根を計算しなければならない．
図 7 FNTT と畳み込み定理を用いた乗算手法
ばならず，数論変換を用いた乗算計算では sp = 1, 2 のみしか
そこで，位数 M の原始根を計算できる素数として，P =
標準ライブラリ関数では対応できなかった．また，64bit で表
sM + 1（s：自然数）で表現できる素数を用いる．前述したよう
現できる素数 P 上で変換や計算を行っても，位数 M の原始根
に，求めたい α の位数が M のとき，素数 P に対して αP −1 ≡ 1
α の計算や剰余計算において以下に示す課題があった．
(mod P ) は成り立ち，さらに，P − 1 が M の倍数である．よっ
4. 3. 1 mod の計算
て，素数 P = sM + 1 としたとき，位数 M の原始根 α を計算
数論変換では 64bit 整数 x，y の積剰余 x ∗ y(modP ) を計算
することができる．この計算方法により，計算したい乗算の桁
を多用しており，sp が 4 以上の場合，素数 P が 64bit で表現
数や畳み込み値に対して適切な素数 P と位数 M の 1 の原始根
できても積 x ∗ y が 64bit を超える場合，積をとって余りを取
α を求めることができる．
る一般的な方法では計算時間が大きい．そこで，除算を用いず
また，素数 P 上で位数 M の原始根 α を計算する際，αM ≡ 1
に乗算・加減算・シフト演算のみで，すべての演算を 64bit で
(mod P ) を計算する必要はなく，α は原始根であることから，
行える積剰余を求めるモンゴメリ乗算 [2] の実装を組み込んだ．
αM/2 ̸≡ 1 (mod P ) なので，αM/2 ≡ −1 (mod P ) つまり，
モンゴメリ乗算とは，P > 0 を法とした合同算術に関して，
演算したい値を，ある定数 R を掛けた表現であるモンゴメリ表
下位 32bit で表現し，a ∗ b を次のように上位 64bit，下位 64bit
現に変換し，この表現によってすべての計算を行った後，最後
に分けて表現することにした．つまり，
に元の領域での表現に逆変換することである．
T1 ∗ 264 + T2 = a ∗ b
モンゴメリ表現での乗算では R が余分に残るので，R−1 を
(25)
掛けて P による剰余を求める処理をモンゴメリリダクションと
となるような T1，T2 を求めるようにする．T1 と T2 の具体
いう．モンゴメリリダクションは，モンゴメリ乗算の基本とな
的な計算方法は，次のコードである．コード中の Carry は下
る演算であり，P を法とする R に関する T (0 <
= T < P R) の
モンゴメリリダクションは，
位桁からの桁上がりである．
MR(T ) = T R
−1
(modP )
(23)
a 1 = ( a & 0 x f f f f f f f f 0 0 0 0 0 0 0 0 ) > >32
a 2 = a & 0 x f f f f f f f f > >32
b 1 = ( b & 0 x f f f f f f f f 0 0 0 0 0 0 0 0 ) > >32
と定義される．ここで，R は R > P および P と互いに素とな
b 2 = b & 0 x f f f f f f f f > >32
任意の整数である．本稿では 64bit の整数の積剰余を対象とし
T 1 = a 1 * b 1 + ( a 1 * b2 > >32)+(b 1 * a2 > >32)+C a r r y
ているため，R を R = 264 とし，2 のべきで表現できる定数に
T 2 = a 2 * b 2 + ( a 1 * b2 < <32)+(b 1 * a2 < <32)
した．また，R−1 は RR−1 ≡ 1( mod P ) を満たすモジュラー逆
その後，モンゴメリリダクションの手続きで，
数である．
モンゴメリリダクションは次の図 8 に示す手続きで計算でき
′
′
る．図 8 中の P は P P ≡ −1(modR) を満たす定数で事前に
t = ((T1 ∗ 264 + T2) + ((T1 ∗ 264 + T2) ∗ P ′ mod R)P(26)
)/R
求めておく必要がある．図 8 より，モンゴメリリダクションの
を計算する．まず，(T1 ∗ 264 + T2) ∗ P ′ を計算することになる
過程で modR と /R の計算がモンゴメリリダクションでは行
が，今 modR の R が 264 以下の値であり，かつ，2 の冪であ
われているが，R を 2 のべきを選ぶことによりビットシフトや
るとすると，(T1 ∗ 264 + T2) ∗ P ′ のうち 2 進数で 64 桁までの
ビットマスクにより効率的に計算できることがわかる．
値しか modR で残らないので，64 桁以下の計算だけすればよ
い．よって，
t = (T+(TP’mod R)P)/R;!
(T1 ∗ 264 + T2) ∗ P ′ mod R = (T2 ∗ P ′ ) mod R
if(t >= P){!
return t-P;!
(27)
となる．つまり，次に計算するのは，T2 ∗ P ′ である．これも
}else{!
a ∗ b と同様に上位 32bit と下位 32bit に分けて計算し，
return t;!
}!
T 2 ∗ P ′ = T3 ∗ 264 + T4
(28)
図 8 モンゴメリリダクション
モンゴメリリダクションを用いてある整数 x をモンゴメリ
表現 Y へ変換するためには，R を掛けて P による剰余を求
めればよいので，あらかじめ R2 = R2 (modP ) を用意して，
X = MR(xR2 ) を求めればよい．モンゴメリ表現から元の整数
となるような T4 を求める．上位 64bit の T3 はその後 modR
を取るので，正確な計算をしなくても modR で求める範囲で
はないので計算をしなくてよい．具体的な計算方法は以下の
コードである．
表現に戻す逆変換は，モンゴメリリダクションそのものであり
T 2 1 = ( T 2 & 0 x f f f f f f f f 0 0 0 0 0 0 0 0 ) > >32
x = MR(X) である．
T 2 2 = T 2 & 0 x f f f f f f f f > >32
以上より，整数 a(0 <
= a < P ) と整数 b(0 <
= b < P ) の乗算剰
余 c = ab(modP ) は，モンゴメリリダクションを用いた変換，
P ’1 = ( P ’& 0 x f f f f f f f f 0 0 0 0 0 0 0 0 ) > >32
逆変換を用いて，
T 4 = T 2 2 * P ’2 + ( T 2 1 * P ’2 < <32)+(P ’1 * T22 < <32)
c = MR(MR(ab)R2 )
(24)
により計算できる．本論文では，すべての計算が 64bit ででき
るように工夫してモンゴメリ乗算を実装した．その実装の手続
きを以下で説明する．
モンゴメリ乗算の実装方法
64bit の整数 a，b の乗算剰余 a ∗ b (mod P ) を計算したいとす
る．まず，式 (24) の中の MR(ab) を計算する．図 8 において，
T を ab に置き換えたものである．a，b はともに 64bit で表さ
れる整数なので a ∗ b は最大で 128bit になってしまう．それで
は 64bit で計算ができないため，a，b をそれぞれ上位 32bit，
P ’2 = P ’& 0 x f f f f f f f f > >32
よって，モンゴメリリダクションの手続きは，
t = ((T1 ∗ 264 + T2) + (T4 mod R)P )/R
(29)
となる．次の (T4 mod R)P の計算も a ∗ b と同様に正確に計
算する必要がある．この計算結果を仮に，(T4 mod R)P =
T5 ∗ 264 + T6 とすると，次に計算するのは，
t = ((T1 ∗ 264 + T2) + (T5 ∗ 264 + T6))/R
(30)
である．ここで，R = 264 とすると，式 (30) は，
t = ((T1 + T5 ∗ 264 + T2) + (T5 ∗ 264 + T6))/R
= (T1 + T5) + (T2 と T6 の加算結果の桁上がり) (31)
となるので，T2 の 63bit 目と T6 の 63bit 目の論理積を T1 と
のとき位数 M48 の原始根は，α48 = 150730 であった．この素
T5 の加算結果に加算することで t を求めることができ，以下
数 P は，N = 2048，sp = 8 のときの畳み込みの最大値に近い
のように実装した．
素数よりも大きい素数である（表 2 参照）．
上記の素数 P48 と位数 M48 の原始根 α48 を用いて他の桁数，
temp1 = T[1]+T[5]
言い換えると他の位数 Mi（i：48 より小さい自然数）における
t e m p 2 = ( ( ( T [2] < <1) > >1)
原始根 αi を，素数 P48 をそのままの値にして計算することがで
+ ( ( T [6] < <1) > >1)) > >63
きる．M は 2 のべきで表現できる数なので，M48 = 248−i Mi
t e m p 3 = T [2] > >63
と表せる．このとき，
t e m p 4 = T [6] > >63
48
αM
≡1
48
C a r r y = ( t e m p 2 + t e m p 3 + t e m p 4 ) > >1
48−i
(mod P48 ) → α248
→
t 1 = ( ( ( ( ( T [1] < <1) > >1)
+ ( ( T [5] < <1) > >1)+C a r r y ) > >63)
Mi
≡1
248−i Mi
(α48
)
i
→ αM
≡1
i
+ ( T [1] > >63)+(T [5] > >63)) > >1
≡1
(mod P48 )
(mod P48 )
(mod P48 )
上式より，素数 P48 上の位数 Mi の原始根は αi = α2
t2 = temp1+Carry
48−i
と
計算できる．例えば，位数 M9 = 29 = 512 のとき，位数 M9
最後に，t か t − P を解とするが，t1 ∗ 264 + t2 は 2P 以下で
あることより，以下のコードに示す条件分岐により計算する．
の原始根 α9 は，
48−9
2
α9 = α48
39
(modP ) = α248 (modP )
549755813888
= α48
(modP ) = 4011842735377954027
i f ( t1 > 0 ) {
t = 0 xff ffffffffffffff - P + 1 + t 2
と計算できる．
} e l s e i f ( t2 > = P ) {
このように大きめの位数 M （原始根の数が多くなる位数）で
t = t2 - P ;
の素数 P を決めておき，位数 M の原始根 α を求めておくこと
} else{
で，1/2，1/3，· · · 倍の M の場合の原始根は素数 P 上のもとで
t = t2;
α2 ，
α3 ，· · · を計算して求めることができる．
}
本論文で用いる整数乗算のためのパラメータの例を上記のよう
ここで，積剰余 a ∗ b(modP ) の計算方法として，標準ライブ
ラリ関数を用いた方法，2 進分解法実装，モンゴメリ乗算実装
を C 言語により実装し，32bit より小さい整数 a，b を決定し，
100000 回繰り返して積剰余を計算した際の計算時間を表 3 に
まとめた．表 3 よりモンゴメリ乗算実装は標準ライブラリ関
数による実装と比べて，64bit 整数同士の計算ができるものの，
に計算し表 4 にまとめた．表 4 の位数 M の各原始根は，実用的
な計算時間（数ミリ秒）で計算できた．表中の N = 4096，8192
のときの P ，α は畳み込み値の最大値が 6195545712378249217
より小さいと分かっているとして，計算時間の推測をするため
に算出したものである．
5. 高速数論変換を用いた乗算手法の実験的評価
計算時間は約 7 倍となっていることがわかる．
5. 1 実 験 環 境
表 3 積剰余計算の比較
実装方法
a*b( mod P) 2 進分解法実装 モンゴメリ乗算実装
×
64bit 整数同士の計算
計算時間 [ms]
15.180
⃝
571.768
高速数論変換と畳み込み定理を用いた乗算手法 (図 7) を C 言語
を用いて実装し，gcc 4.2.1 コンパイラを用いて，CPU:3.33GHz
110.695
6-Core Intel Xeon の計算機上に実現しシミュレーション実験
4. 3. 2 位数 M の原始根 α の計算
を行った．実験で用いた整数は，3.4 節で行った高速フーリエ
乗算対象となる整数の桁数，分割桁数から素数 P を計算する
変換と畳み込み定理を用いた乗算手法のシミュレーション実
ことは簡単である．しかし，桁数や分割桁数を大きくとると素
験で用いた整数と同じ整数の組 Data1∼Data5 である．本論
数 P が大きくなり，小さい位数 M の原始根 α は素数 P の範
囲に高々M 個しかないため，4.2 節で述べたように一般的な位
文では，N = 2i (10 <
= i <
= 13) となる各桁数 N の分割桁数
sp = 1, 2, 4, 8 に対して乗算結果を計算する実験を行った．本論
数 P − 1 の原始根の計算方法では計算時間が莫大となり，α を
文で示す実験では，積剰余の実装方法とパラメータの導出方法
求めることが難しい．本論文では sp = 8 の場合では，文献 [7]
を分割桁数によって適する方法を選択した．使用した積剰余の
に示唆されている特定の P の範囲での位数 M の原始根の計算
実装，高速数論変換のパラメータの素数 P と原始根 α は次の
を一般化し，α を計算しやすくするために本来計算したい位数
通りである．モンゴメリ乗算を用いる場合は，モンゴメリ乗算
M よりも大きくして，素数 P を求めておき，それを使って位
で用いるパラメータ R2，P ′ を事前に計算したが実際の値につ
′
′
数 M の原始根 α を計算する方法を実装し，実用的な時間で
′
位数 M の原始根 α を計算する．
本 論 文 で は ，位 数 M48 = 248 と し ，P48 = sM48 + 1
（s：自 然 数 ）の 形 で 表 現 す る こ と が で き る 素 数 P48 =
6195545712378249217（s = 22011）を用いることにする．こ
いては省略する．
・分割桁数 sp = 1，2 のとき
積剰余：ライブラリ関数の ∗，%をそのまま利用
パラメータ：表 2 で示した各桁数 N に対応する素数 P と原始
根 α の値
表 4 実験で用いた乗算を行うためのパラメータ一覧 (大きめの素数 P をとったとき)
桁数 N
位数 M
畳み込み値の最大値に近い素数 P
1024
256
2559999948800007169
2048
512
5119999897600035841
4096
1024
10239999795200005121
8192
2048
64bit で表現できない
大きめの素数 P
原始根 α
1495320057967618216
6195545712378249217
4011842735377954027
37690190028693632
2128855843988593852
※ N = 4096，8192 のときの P ，α は畳み込み値の最大値が 6195545712378249217 より小さいと分
かっているとして算出
・分割桁数 sp = 4 のとき
精度な計算ができない．また，畳み込み演算においても分割桁
積剰余：64bit モンゴメリ乗算による実装
数 sp = 8 の次に大きい分割桁数 sp = 16 では標準ライブラリ
素数 P ，原始根 α：表 2 に示す一般的な計算方法により算出
の対応範囲を越えるため計算できない．高速数論変換では，実
した P ，α
験で行った最大の分割桁数 sp = 8 のときの，畳み込み値の最
・分割桁数 sp = 8 のとき
大値より非常に大きい素数 P で計算したのにも関わらず，畳み
積剰余：64bit モンゴメリ乗算による実装
込み演算よりも高速に計算可能であった．sp = 1, 2 の結果の傾
素数 P ，原始根 α：表 4 に示す位数 M を大きめにとって算出
向から，畳み込み値に近い素数 P ，原始根 α が計算できたとき
した P ，α
でも，桁数 N = 4096，8192 において畳み込み演算よりも高速
5. 2 実験結果と考察
に計算できる可能性が高いと考えられる．よって，高速数論変
5. 2. 1 計 算 時 間
換ではこれまでに述べた mod 計算，素数 P や原始根 α の計算
高速数論変換と畳み込み定理を用いた乗算手法による乗算計
の課題を克服することができれば，実験結果の傾向からさらな
算時間の平均値と畳み込み演算，FFT を用いた手法の平均計
5. 3 計 算 精 度
算時間 (表 1) を合わせた表を表 5 に示す．
畳み込み演算を用いた手法で求まった乗算結果が正しいと仮
表 5 乗算の計算時間の比較
桁数
(10 進)
N
定したとき，高速数論変換を用いた乗算手法では計算できたす
計算時間 [ms]
FFT を用いた乗算手法
べての桁数 N ，分割桁数 sp において正しい乗算結果を得るこ
FNTT を用いた乗算手法
a ∗ b(modP )
sp = 2
る分割桁数を大きくとることによる高速化の可能性がある．
モンゴメリ乗算実装
sp = 4
とができた．
sp = 1
sp = 2
sp = 4
sp = 1
sp = 8
1024
0.604
0.338
0.224
0.736
0.387 0.453
0.322
2048
1.125
0.619
0.384
1.385
0.731
0.789
0.501
4096
2.301
1.226
0.713
2.788
1.461
1.652
0.944
高速数論変換と畳み込み定理を用いた手法で，大きな分割桁
8192
4.651
2.404
1.409
5.749
2.946
3.365
1.874
数に対応可能にすれば，高速かつ正確な乗算ができると推測で
6. まとめと今後の課題
きると考えられる．また，標準ライブラリで精度が保証されて
表 5 から，高速数論変換と畳み込み定理を用いた手法の場合
でも高速フーリエ変換を用いた手法の場合と同様に分割桁数を
大きくとると高速に計算できた．しかし，分割桁数 sp = 2 か
いない計算の範囲を超えて非再帰的実装する際の 2 つの課題を
克服した．
今後の課題として，96bit で素数を表現，モンゴメリ乗算の拡
ら sp = 4 への計算時間の関係が sp = 1 と sp = 2 の実験結果，
張を行うことでさらに大きな分割桁数でも対応可能にし，非再
および高速フーリエ変換の実験結果でも見られるような分割桁
帰的実装が高速化につながる限界を調査することが挙げられる．
数を 2 倍にすると計算時間が約 1/2 倍とならなくなった．この
ような結果となった理由として，64bit の整数同士の積剰余を
計算できるようにモンゴメリ乗算実装を組み込んだためである
ことが挙げられる．モンゴメリ乗算実装により 64bit 整数同士
の計算が可能になったが計算時間が約 7 倍もかかっており，分
割桁数を増加させることによる高速化はできなかった．また，
高速数論変換は高速フーリエ変換と同じ計算フレームワークで
計算しているため計算時間も桁数増加に対して計算時間の増加
も緩やかである．
次に，高速数論変換と畳み込み定理を用いた手法と，高速
フーリエ変換を用いた手法と畳み込み演算を用いた手法を比較
する．実験で扱った最大桁数 N = 8192 のとき最も高速に乗算
結果を得られるのは，分割桁数 sp = 4 にし高速フーリエ変換
と畳み込み定理を用いた手法である．高速フーリエ変換を用い
た手法の非再帰実装では，sp = 4 より大きい分割桁数では高
文
献
[1] J. W. Cooley and J. W. Tukey , “An algorithm for the machine calculation of complex Fourier series,” Math. of Comput., 19, pp. 297-301, 1965.
[2] P. L. Montgomery, “Modular multiplication without trial division,” Mathmatics of Computation, vol.44, No.170, pp.
519-521, 1985.
[3] M. Fürer, “Faster Integer Multiplication,” Proceedings of
the 39th Annual ACM Symposium on Theory of Computing, pp. 57-66, 2007.
[4] A. Schönhage and V. Strassen,“Schnelle Multiplikation
großer Zahlen,” Computing 7, pp. 281-292, 1971.
[5] 高木貞治，
“ 初等整数論講義 第 2 版 ”，共立出版，1994 年．
[6] 辻井重雄，鎌田一雄，
“ ディジタル信号処理 ”，昭晃堂，1993 年．
[7] 後保範，
“ 高速剰余変換による多数桁乗算 ”，情報処理学会論文
誌，Vol.44，No.12，2003 年．
[8] 和田秀男，
“ コンピュータと素因子分解 ”，遊星社，2001 年．