DEIM Forum 2016 D8-5

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

上土井陽子†

若林 真一†† 広島市立大学大学院情報科学研究科 〒 731–3194 広島市安佐南区大塚東三丁目 4 番 1 号
E-mail: †shotah@lcs.info.hiroshima-cu.ac.jp, {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 暗号，楕円曲線暗号のような公開鍵暗号システムの出現によって，乗算の計算量の低減に関する研究は大きな重要性をもつことになった．
上記のような暗号技術で用いる桁数の大きな整数に対する高′A = a0 × B 0 + a1 × B 1 + · · · · · · + aN ′ × B N −1

(1)′
次に ai = 0 (N ′ <= i < 2N ，i は整数) とし，式 (1) の係数

系列を式 (2) のように表現する．

速乗算アルゴリズムとして，Fürer の乗算アルゴリズム [3] の( a0 , a1 , · · · , aN ′ −1 , 0, 0, · · · , 0)

中でも高速フーリエ変換 (fast Fourier Transform, FFT) と畳さらに，係数系列の各要素を時間の関数 f (t) の関数値とすみ込み定理に注目し，基本となっているアルゴリズムを実装し(2)

ると，式 (3) のように表現できる．

てシミュレーションによる実験的性能評価を行った．その結果，多倍長整数乗算を通常の筆算で行った場合に比べて，大幅に計算時間の高速化ができていることがわかった．しかし，円周率

π や三角関数 sin, cos などの無理数の計算を必要とする複素数体上のフーリエ変換では乗算結果に誤差が生じる場合があることがわかった．乗算結果を暗号技術で用いる場合には高速に計( x(0), x(1), · · · , x(N ′ − 1), · · · , x(2N ′ − 1) )式 (2)，式 (3) の間には，次のような関係が成り立っている．

atx(t) =
0

算できることよりも，高精度な計算ができる方が重要度が高いと考えられる．
本稿では，誤差が生じることがない有限体上での変換である
数論変換と同じ計算フレームワークである高速数論変換 (fastNumber 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] を利()*)+! ",)(+ !,#%+ $!(%+ "*%,*+ %!!)+ ),,"+ "*,,+,-./! ")!

一般に，離散フーリエ変換 (DFT) を用いて時間領域の信号

x(n) を周波数領域の信号 X(k) に変換するためには膨大な計算

(!!

#%!

"!!

図2

筆算による乗算)"!

,,!

用する．FFT は，M = 2i (i は自然数) 点 DFT の計算量を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 Fouriertransform, IFFT) という．

,-"!

1 ∑
−nk
X(n) WMM
M −1x(n) =;<"!

;<.!

!"#$%&'!
("#$%)"#&%)***)"#'(&%+!

!"#$%&'!
(.#$%).#&%)***).#'(&%+!

(8)

))*!

k=0

))*!

(+#$%)+#&%)***)+#'(&%+!

./!
"!

1234
-"!

#$'!#$!%!

#$(!#$!!!

(!+#$%,#$%)+#&%,#&%)***)+#'(&%,#'(&%!+!

#$!!#$&!

/))*!

#$%!#$$!
01!

#$)!#$!'!

56789:%!"#$%&'!
(#"-.%#$%)#"-.%!#&%)***)#"-.%!#'(&%+!

!!

-!!

#$*!#$!(!

,-"!

=>?@A!

#$+!#$!*!
-"!

(,#$%),#&%)***),#'(&%+!

#$./%0!

BCDEF"!-!.!

#$,!#$!)!

図 6 FFT と畳み込み定理を用いた乗算手法

図 4 回転因子 W8
X(0)!

x(0)!

W 80 !
x(2)!

W 80!

X(1)!

W 84!

W 82!
W 80 !

x(6)!

W84!

W 82!

W86!

W83 !

W 80 !

W 84!

W8 2 !

W85 !

W8 4 !

W8 6 !

x(5)!
x(7)!

W8 6 !

W87!

X(5)!

W 84!
W80 !

図5

ミュレーション実験を行った．実験で用いた整数は，10 進数でみなし生成した整数である．本稿では，N = 2i (10 <
=i<
= 13)
となる各桁数 N の分割桁数 sp = 1, 2, 4, 8 に対して，異なる 5
つの整数の組 Data1〜Data5 に対して乗算結果を計算する実験

X(6)!
X(7)!

W 84!

て，CPU:3.33GHz 6-Core Intel Xeon の計算機上に実現しシ

1 から 9 の数をランダムに N 個並べた数値を N 桁の 10 進数と

X(4)!

W 80!

理を用いた乗算手法 (図 6) の 2 つの乗算手法を C 言語を用いて実装し，gcc 4.2.1 コンパイラ，最適化オプション-O3 を用い

X(3)!

W 84!

x(1)!
x(3)!

W 81!
X(2)!

x(4)!

3. 4 実験的性能評価 (1)
畳み込み演算を用いた乗算手法 (図 3)，FFT と畳み込み定

W 80!

を行った．

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)](10)

式 (10) から，畳み込みをフーリエ変換と逆フーリエ変換を使

算時間が約半分になっていることから，sp を増加させることによる高速化の可能性がある．
表 1 各手法の乗算平均計算時間 (表中の灰色の部分は結果に誤差有)
桁数

計算時間 [ms](10 進)

畳み込み演算 (筆算)FFT と畳み込み定理

N

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. 高精度な多倍長整数乗算手法の非再帰的実装に乗算を計算できる可能性がある．

4. 2 数論変換のためのパラメータ P ，α の導出

FFT は，実験結果より桁数 N や分割桁数 sp によって計算

式 (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

乗算を行うためのパラメータ一覧 (位数 p − 1 の原始根を求める一般的な計算方法による算出)

桁数(10 進)

乗算を行うためのパラメータ一覧
分割桁数 sp = 1

分割桁数 sp = 2

N

素数 P

素数 P

1024

176129153 10039297

3201

2048

331777386 20078593

4096

737281

96 40148993

8192 1376257

73 80330753

原始根 α

分割桁数 sp = 4

原始根 α

素数 P

分割桁数 sp = 8

原始根 α

素数 P

原始根 α

2559999948800007169△

47486 102379527169 244863846

5119999897600035841△

15170 204759054337

55931003

10239999795200005121△

46283 409518112769

24724426

64bit で表現できない△

51189761537 182252995※ △ は 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) の位桁からの桁上がりである．

モンゴメリリダクションは，

a 1 = ( a & 0 x f f f f f f f f 0 0 0 0 0 0 0 0 ) > >32
MR(T ) = T R−1(modP )

(23)

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 Mit 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 ) → α248Mi≡1248−i Mi→ (α48) ≡1

t 1 = ( ( ( ( ( T [1] < <1) > >1)+ ( ( T [5] < <1) > >1)+C a r r y ) > >63)

i
→ αM
≡1
i

+ ( T [1] > >63)+(T [5] > >63)) > >1

(mod P48 )(mod P48 )(mod P48 )上式より，素数 P48 上の位数 Mi の原始根は αi = α2t2 = temp1+Carry

48−i

と

計算できる．例えば，位数 M9 = 29 = 512 のとき，位数 M9

最後に，t か t − P を解とするが，t1 ∗ 264 + t2 は 2P 以下であることより，以下のコードに示す条件分岐により計算する．

の原始根 α9 は，
48−92
α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 の範

文では，N = 2i (10 <= i <= 13) となる各桁数 N の分割桁数

囲に高々M 個しかないため，4.2 節で述べたように一般的な位

sp = 1, 2, 4, 8 に対して乗算結果を計算する実験を行った．本論

数 P − 1 の原始根の計算方法では計算時間が莫大となり，α を文で示す実験では，積剰余の実装方法とパラメータの導出方法

求めることが難しい．本論文では sp = 8 の場合では，文献 [7]を分割桁数によって適する方法を選択した．使用した積剰余のに示唆されている特定の P の範囲での位数 M の原始根の計算

実装，高速数論変換のパラメータの素数 P と原始根 α は次のを一般化し，α を計算しやすくするために本来計算したい位数通りである．モンゴメリ乗算を用いる場合は，モンゴメリ乗算

M よりも大きくして，素数 P を求めておき，それを使って位で用いるパラメータ R2，P ′ を事前に計算したが実際の値につ′′数 M の原始根 α を計算する方法を実装し，実用的な時間で′位数 M の原始根 α を計算する．
本 論 文 で は ，位 数 M48 = 248 と し ，P48 = sM48 + 1

いては省略する．
・分割桁数 sp = 1，2 のとき
積剰余：ライブラリ関数の ∗，%をそのまま利用（s：自 然 数 ）の 形 で 表 現 す る こ と が で き る 素 数 P48 =パラメータ：表 2 で示した各桁数 N に対応する素数 P と原始6195545712378249217（s = 22011）を用いることにする．こ根 α の値表 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 乗算の計算時間の比較
桁数定したとき，高速数論変換を用いた乗算手法では計算できたす

計算時間 [ms]FFT を用いた乗算手法

べての桁数 N ，分割桁数 sp において正しい乗算結果を得るこ

FNTT を用いた乗算手法a ∗ b(modP )モンゴメリ乗算実装

N

sp = 1

sp = 2

sp = 4

sp = 1

sp = 2

sp = 4

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

数に対応可能にすれば，高速かつ正確な乗算ができると推測で(10 進)とができた．

6. まとめと今後の課題

きると考えられる．また，標準ライブラリで精度が保証されて表 5 から，高速数論変換と畳み込み定理を用いた手法の場合でも高速フーリエ変換を用いた手法の場合と同様に分割桁数をいない計算の範囲を超えて非再帰的実装する際の 2 つの課題を克服した．

大きくとると高速に計算できた．しかし，分割桁数 sp = 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. Fürer, “Faster Integer Multiplication,” Proceedings ofthe 39th Annual ACM Symposium on Theory of Computing, pp. 57-66, 2007.
[4] A. Schö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 年．