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