FMT(高速剰余変換）の応用

　　　　　　　　2002/12/20 日立製作所 & 早稲田大学　後保範 ( Ushiro Yasunori )
　　　　　　　　　2001/6/9　作成のものをweb化
　　　　　　　　2003/05/20 更新（自然数巡回FMT及び分割乗算等を追加) 後保範

1. 概要と目的

　FMT原理で示した方法の多数桁演算への応用を述べる。
　多数桁係数を持つ行列計算（行列乗算、連立一次方程式、固有値計算）や暗号解読への応用
はここでは述べない。
　FMTを利用したソースプログラムは多数桁計算プログラム、 π計算プログラムを参照。　　

2. 整数FMTと複数個のFMTの重ね合せ

　FFTは原始n乗根 として複素数ω=e^-2πi/nを使用したものといえる。一方、FMTは
FFTと同様に複素数でも定義できるが、整数でも可能である。ここでは、説明を簡単にするため
FMTの要素数nは偶数とする。

2.1 概要

　　FMTは１の原始n乗根となるωの上で定義される。そのため、mod(ωⁿ,P)=1で、 mod(ω^k),P)≠1,
　k=1,2,…,n-1となる整数n,ω,Pでも可能となる。
　　ω=2,P=ω^n/2+1とするとωはmod(P)上で原始n乗根となることが知られている。一般に、ωが2
　以上の整数ならP=ω^n/2+1上でωは原始n乗根となり整数FMTが構成できる。
　　また、適当な大きさの素数P上で原始n乗根となるωが多数存在する。下記に現在の計算機で
　扱い易い素数Pとn及びωの値を示す。この場合、複数個のPの上での多数桁乗算を 百五減算
　で合成し、より大きい値の上での多数桁乗算が構成できる。　　

2.2 整数ωとmod(ω^n/2+1)を利用

　　nを偶数とし、ωが2以上の整数とすると、P=ω^n/2+1上でωは1の原始n乗根となる。
　これを利用すると、要素数nのFMTが構成できる。
　nが大きくなるとPは指数的に増大し、一般に計算機で扱いづらいと言われているが、2段階FMT
　による多数桁の乗算等で利用できる。

2.3 倍精度計算可能な原始根

　　Pが2²⁴の近くの原始n乗根となる整数n,ω,Pの例を表1に示す。
　

表1. Pが2²⁴の近くのFFTは**原始n乗根**となる整数n,ω,Pの例
番号	n = α = 2¹⁷		n = β = 3・2¹⁵		n = γ = 3²・2¹⁴
NO.	ω	P	ω	P	ω	P
1	770	149 α + 1	26	173 β + 1	175	114 γ + 1
2	201	153 α + 1	866	187 β + 1	434	132 γ + 1
3	70	155 α + 1	562	194 β + 1	12	136 γ + 1
4	20	164 α + 1	637	204 β + 1	604	147 γ + 1
5	1498	165 α + 1	100	220 β + 1	459	150 γ + 1

2.4 倍精度2個で計算可能な原始根

　　Pが2⁴⁸の近くの原始n乗根となる整数n,ω,Pの例を表2に示す。
　

表2. Pが2⁴⁸の近くのFFTは**原始n乗根** となる整数n,ω,Pの例
番号	n = α = 2³⁵		n = β = 3・2³²		n = γ = 3²・2³¹
NO.	ω	P	ω	P	ω	P
1	360	8319 α + 1	9375	22100 β + 1	8426	14636 γ + 1
2	933	8334 α + 1	5549	22226 β + 1	226	14771 γ + 1
3	1714	8549 α + 1	9520	22254 β + 1	4987	14896 γ + 1
4	474	8755 α + 1	3991	22314 β + 1	9607	14913 γ + 1
5	744	8759 α + 1	6955	22316 β + 1	7806	14994 γ + 1

2.5 複数個のFMTの重ね合せ

　　表1及び表2の複数個のPによる、FMTを用いた畳み込み演算結果は、複数個のPは互いに素な
　なため、百五減算(中国剰余定理)で畳み込み演算結果の各要素の桁数を拡大することができる。

3. 正巡回FMTと負巡回FMTによる畳込み演算

3.1 概要

　　n個の要素の畳込み演算は、両方の入力にn個のゼロ要素を追加して、FMT変換又はFFT変換
　を行い、互いの2n要素を対応する要素単位に乗算し逆変換することで2n個の結果が得られる。
　もし、入力要素を2n桁に拡張しないで、そのままn桁のまま順変換、要素単位乗算及び逆変換
　を行うと畳込み演算結果のn個の要素はどのようになるであろうか。
　　一般に行なわれているFFTで計算すると、2n要素の結果がn要素に後半のn要素が前半のn要素に
　加算されて求まる。計算方法を変えると後半のn要素が前半の要素から減算されて求まる。
　FFTの代りに、FMTを利用した場合は、前半と後半の要素がFFTと逆転する。
　　加算されて求まるとき正巡回FMT(FFT)、減算のとき負巡回FMT(FFT)と呼ぶ。

3.2 畳込み演算の方法

　　下記のように表現されるf,gの畳み込み演算結果hを正巡回と負巡回FMTで計算する方法を述べる。
　　　　f = a_n-1E^n-1 + …… + a₁E + a₀ , 　　g = b_n-1E^n-1 + …… + b₁E + b₀ 　　　　　　--- (1)
　　　　h = f・g = c_2n-1E^2n-1 + c_2n-2E^2n-2 + …… + c₁E + c₀　　　　　　　　　　　　　　--- (2)
　　正巡回h⁺と負巡回h^-のf・g畳込み演算結果を下記で表す。
　　　　h⁺ = mod (h , Eⁿ - 1) = d_n-1E^n-1 + …… + d₁E + d₀
　　　　h^- = mod (h , Eⁿ + 1) = e_n-1E^n-1 + …… + e₁E + e₀　　　　　　　　　　　　　　　--- (3)
　　このとき、d_j, e_jとc_j, c_j+nは(3)式の定義から次の関係があることが分かる。
　　　　d_j = c_j + c_j+n
　　　　e_j = c_j - c_j+n , j = 0, 1, …… ,n-1　　　　　　　　　　　　　　　　　　　　　　　　　　--- (4)
　　正巡回順FMTをFMT()で、正巡回逆FMTをFMT^*()で表し、負巡回順FMTをFMTI()で、
　　負巡回逆FMTをFMTI^*()で表すとh⁺, h^-は下記のように計算できる。
　　ここで、×印は要素別乗算を示す。
　　　　h⁺ = FMT^*( FMT(f)×FMT(g) )
　　　　h^- = FMTI^*( FMTI(f)×FMTI(g) )　　　　　　　　　　　　　　　　　　　　　　　　　--- (5)
　　h⁺, h^-の係数から求める2n個の畳込み演算結果hの係数c_j , c_j+nは下記で求められる。
　　　　c_j = ( d_j + e_j ) / 2
　　　　c_j+n = ( d_j - e_j ) / 2 , j = 0, 1, …… ,n-1 　　　　　　　　　　　　　　　　　　　　　--- (6)

3.3 負巡回FMTの計算方法

　　正巡回FMTはn要素のFMTそのままの計算であり、負巡回FMTの計算は下記のように正巡回FMT
　に変更を加えて求める2種類の方法がある。
　これは、Eⁿ-1及びEⁿ+1が下記のように原始n乗根ωを使用して分解できるためである。
　　　( Eⁿ - 1 ) = (E - ω^n-1)(E - ω^n-2) …… (E - ω¹)(E - ω⁰)
　　　( Eⁿ + 1 ) = (E - ω^n-1/2)(E - ω^n-3/2) …… (E - ω^3/2)(E - ω^1/2)

　(1) 原始n乗根ωのテーブルを変更する方法
　　下記の様に正巡回のωのテーブルを負巡回のωのテーブルに変更する。
　　逆変換を行うには、テーブル上のωをω^-1に変更する。
　　ωのテーブルの与え方がFMTとFFTでは逆順になることに注意を要す。
　　　　正巡回のωのテーブル： ( ω^n-1, ω^n-2, …… , ω¹, ω⁰ )
　　　　負巡回のωのテーブル： ( ω^n-1/2, ω^n-3/2, …… , ω^3/2, ω^1/2 )
　(2) 入力係数及び出力係数に変換係数を掛ける
　　(1)式の入力係数a_k,b_kに(3)式の出力係数a_kに下記を施して正巡回の計算を適用すると負巡回
　　FMTの結果が得られる。
　　　　入力係数変換： a_k = a_kω^k/2 , 　 b_k = b_kω^k/2
　　　　出力係数変換： e_k = d_kω^-k/2 , 　　　 k = 0,1,...,n-1　

4. 2段階FMTを使用した畳込み演算

　FMTやFFTを利用した時、メモリ的に問題となる多数桁の乗算を、2段階FMTを使用することにより、
FMTの計算効率を損なわずに、使用メモリ量を減少させる効果がある。

4.1 概要

　　１要素が多数桁を持つN要素の乗算を、P=ω^N/2+1を法とする整数FMTを用いて実施する。
　ここでは、Pを法とするFMTを上位FMT、要素ごとの項別乗算で必要たなるFMTを下位FMTと呼ぶ。
　　上位FMTは整数演算を使用し、ω^kの乗算は桁上げ処理を除き整数の加減算だけで構成できる。
　8バイト整数が利用できる場合はω=2⁶⁰(2進計算)又はω=10¹⁶ or 10¹⁸(10進計算)とし、4バイト整数
　の場合にはω=2²⁸(2進計算)又はω=10⁸ or 10⁹に選んで計算するのが、メモリ的にも速度的にも
　効率がよい。

4.2 利用機能

　　2段階FMTでは、1要素が多数桁で構成される上位FMTと、その1要素の乗算に使用される下位FMT
　がある。その各々に必要なFMTの機能を示す。
　(1) 上位FMT
　　・ ωが2以上の整数で、P=ω^N/2を法とする要素数Nの整数FMT
　　・正巡回FMTと負巡回FMTによる畳み込み演算
　　・整数8バイト又は整数4バイトの加減算処理(桁上げ処理を除く)で構成
　　・要素ごとの項別乗算は下位FMTを利用
　　・ノード間並列化が必要な場合は上位FMTだけ並列化
　　・スーパーコンでは１要素が百万桁に達する演算が可能
　(2) 下位FMT
　　負巡回FMTを利用し、多数桁乗算結果が自動的にmod(P)を満たすようにする
　　・上位FMTによる多数桁乗算の要素単位に行う項別乗算に利用する
　　・実数FMTによる多数桁の乗算が現在のところ一番高速である
　　・計算は複素数16バイト又は実数8バイトを利用する
　　・整数乗算に実数値を使用するため、丸め誤差が発生し結果は整数化が必要
　　・ノード間並列化ではこの部分の並列化は不要

4.3 計算方法

　　FMTを利用し多数桁演算の基底は2進数である必要はなく、10進数でも他の基底数でも可能で
　あるが、説明を簡単にするため2進数を仮定して説明する。基底を変更するときは、ωの値を
　基底数のべき乗にする必要がある。
　　FMTの順変換及び逆変換の計算方法はFMT原理を参照。

4.3.1 記号の定義

　　ここで使用する記号を下記のように定義する。
　　　N ：　上位FMTの要素数　　　　　　　　　　　　　　　　　　n ：下位FMTの要素数
　　　F_P, F_P^-1 ：上位正巡回FMTの順変換及び逆変換
　　　F_N, F_N^-1 ：上位負巡回FMTの順変換及び逆変換
　　　f_N , f_N^-1 ：下位負巡回FMTの順変換及び逆変換
　　　・：通常の乗算(畳込み演算による全要素乗算）　　　　　× ：項別乗算

4.3.2 上位FMTの計算

上位FMTの要素数をNとしたとき、上位FMTの１要素に2進m・N桁を記憶できるようにする。
　このとき、上位FMTはω=2^2mでP=ω^N/2+1として構成する。
　共にN個の要素を持つAとBの多数桁乗算を下記の様に表す。乗算結果は2N個の要素になる。
　　　A ・ B = ( c_2N-1 , c_2N-2 ,……, c_N ,……, c₁ , c₀ )　
　AとBの多数桁乗算を2N個の要素のFMTではなく、N個の要素の正巡回FMTによる乗算と、
　負巡回FMTによる乗算を利用して下記の様に計算する。
　　　正巡回乗算： F_P^-1 ( F_P ( A ) ×　F_P ( B ) ) = ( d_N-1 , d_N-2 ,……, d₁, d₀ )
　　　負巡回乗算： F_N^-1 ( F_N ( A ) ×　F_N ( B ) ) = ( e_N-1 , e_N-2 ,……, e₁, e₀ )
　AとBの多数桁乗算の2N個の要素c_k+N, c_kは下記で計算できる。
　　　c_k+N = ( d_k - e_k ) / 2
　　　c_k　　= ( d_k + e_k ) / 2
　　　　　　k = 0, 1,……, N-1
　正巡回乗算と負巡回乗算で現われる要素単位の項別乗算(×)は下位FMTを利用して行う。

4.3.3 下位FMTの計算

　　上位FMTの要素単位に実施される項別乗算に下位FMTは利用される。そのため下位FMT
　を利用する多数桁の乗算はmod(P)の上で実行される必要がある。
　P=ω^N/2+1=2^mN+1のため、下位FMTの要素数をn 、下位要素の１要素に詰める 2進桁数をhとすると、
　下記関係が必要となる。
　　　h・n = m・N
　　この条件を満たせば、下位FMTに実数FMTを使用するか、整数FMTを使用するかは自由である。
　要素数nでこの条件の下に、多数桁の乗算結果がmod(P)を満たすためには、負巡回FMTを下位
　FMTに利用すればよい。即ち、aとbの2進h桁を１要素とするh・n桁の乗算は下記の様にn要素の
　負巡回FMTを使用することにより、mod(P)のもとでh・n桁となり求まる。
　　　mod( a・b , P ) = f_N^-1 ( f_N ( a ) ×　 f_N ( b ) )

4.4 具体的計算サイズ

　　入力A,Bは0.5兆桁で出力Cが1兆桁(10進換算)となる3種類の多数桁乗算の使用メモリ量を
　比較する。2段階FMT(実際は2段階FFT)は2002年の1兆2411億桁の円周率世界記録更新に
　使用した多数桁乗算方式である。それまでの世界記録には実FMT(FFT)+Karatsuba法が使用
　されていた。
　(1) 2段階FMT(整数8バイトに2進60桁詰めで計算可能)
　　　FFT(A)及び計算結果　：　0.5TB (FFT用ワーク、テーブルは不要)
　　　FFT(B) 　　　　　　　　　　：　0.5TB
　　　正巡回結果の保存　　：　0.25TB
　　　Work(下位FFT及び通信)：　0.05TB
　　　合計　　　　　　　　　　　　：　1.3TB
　(2) 実FMT(実数8バイトに2進9桁詰めが限度)
　　　FFT(A)及び計算結果　：　3.0TB (1要素２進9桁、入力の12倍)
　　　FFT(Ｂ) 　　　　　　　　　：　3.0TB
　　　FFTワーク　　　　　　　：　3.0TB (MATRIX/MPP)
　　　合計　　　　　　　　　　：　9.0TB
　(3) 実FMT+Karatsuba法(Karatsuba法の適用回数でメモリ量と演算量が変わる)
　　　1回カラツバ法を使用すると、FFT領域は半減し、Karatsuba法のワークが必要。
　　　　　　　　　　FFT　ワーク　合計　　　演算量比
　　　1回適用： 4.5TB + 1.0TB = 5.5TB　　　1.5 倍
　　　2回適用： 2.3TB + 1.5TB = 3.8TB　　　2.3 倍
　　　3回適用： 1.2TB + 1.7TB = 2.9TB　　　3.4 倍
　　　4回適用： 0.6TB + 1.8TB = 2.4TB　　　5.1 倍

5. 複素FMTによる実数の畳込み演算

　　FMT原理の「5. 畳込み演算とFMTとの関係式」で示した(10)式の応用でq=1/4のケースである。
　FMT原理の(5),(6)式と(10),(11)式にq=1/4を代入した式を下記に示す。
　　　f = a_n-1E^n-1 + …… + a₁E + a₀ , 　g = b_n-1E^n-1 + …… + b₁E + b₀
　　　h = f ･ g = c_2n-1E^2n-1 + c_2n-2E^2n-2 + …… + c₁E + a₀　
　　　h^(1/4)　= mod ( h, Eⁿ - ω^n/4 ) = c^(1/4)_n-1 E^n-1 + …… + c^(1/4)₁ E + c^(1/4)₀
　　　c^(1/4)_j　= c_j + i ･ c_j+n , 　　　j = 0, 1, …… , n-1　
　ここで、a_j, b_j, c_j及びEを実数とし、ω = e^2πi/nとすると、ω^n/4 = i (虚数単位)となり、c^(1/4)_jは
　c_jを実数部に持ち、c_j+nを虚数部にもつ複素数となる。即ち、n要素(E進n桁)の実数乗算は、本複素
　FMTを適用したn要素の複素数畳込み演算結果の虚数部に上位桁が、実数部に下位桁が求まる。
　　　更に、下記の様にn/2-1次に畳込む、すなわちn/2要素のFMTにq=1/4を適用すると負巡回
　多数桁乗算が直接得られる。
　　　h^(1/4,n/2)　= mod ( h, E^n/2 - ω^n/4 ) = c^(1/4,n/2)_n/2-1 E^n/2-1 + …… + c^(1/4,n/2)₁ E + c^(1/4,n/2)₀
　　　c^(1/4,n/2)_j　= ( c_j - c_{j+n )} + i ･ ( c_j+n/2 - c_j+3n/2 )　, 　　　j = 0, 1, …… , n/2-1　
　これは、n要素(E進n桁)の実数乗算は、本複素FMTを適用したn/2要素の複素数畳込み演算結果
　の虚数部に負巡回の上位桁が、実数部に負巡回の下位桁が求まる。　

6. αで巡回する多数桁乗算

　FMT原理の (11)式を使用する。再度、この式を下記に示す。
　　　c^(q)_j　　　= c_j + c_j+n ω^qn
　　　c^(q+n/2)_j = c_j - c_j+n ω^qn
　　　　　　　　　q ：有理数 ,　　　j = 0, 1, …… , n-1
　この式は、要素数nに対して法Pと自然数αに対してα=mod(ω^qn,P)なる有理数q と自然数P,ω,ω^q
　が存在すれば、+α及び-αで巡回する多数桁乗算となる。
　表3にn=2²⁰で+-αで巡回するFMT乗算の係数を示す。

表3. n=2²⁰で+-αで巡回するFMTの係数の例
α	P	q	ω	ω^q
2	280,049,478･n+1	1/93,349,826	24,328,030,273,313	196,084,934,801,470
3	322,504,032･n+1	1/80,626,008	27,808,6974,322,725	45,280,667,942,168
4	283,546,518･n+1	1/94,515,506	237,787,923,135,329	113,109,652,731,064
5	205,984,108･n+1	1/205,984,108	43,678,380,610,349	192,726,477,692,165
6	282,686,616･n+1	1/282,686,616	20,159,304,809,371	288,630,981,189,601
7	222,921,552･n+1	1/111,460,776	161,161,411,978,886	163,047,718,620,128
8	28,0049,478･n+1	1/93,349,826	243,728,030,273,313	70,308,039,835,214
9	322,504,032･n+1	1/40,313,004	93,762,856,436,902	259,915,309,872,599
10	314,771,758･n+1	1/314,771,758	149,271,167,055,472	9,402,286,612,793

7. 分割畳込み演算

　　多数桁乗算のFMTを利用した分割乗算に付いて記述する。FMT原理の(11)式を応用する。
　n要素(1要素はE進１桁)同士で結果が2n要素の多数桁乗算を、m=2^l個に分割した乗算は
　s=2n/m要素FMTを利用し下記の原理で行なう。qは有理数である。
　　　h^(q,s) = mod(h, E^s - ω^qn) = mod(f ･ g, E^s -ω^qn)
　　　　　　= c^(q,s)_s-1 + …… + c^(q,s)₁ + c^(q,s)₀
　このとき、c^(q,s)_jと多数桁乗算の係数c_jの関係式は次式の様になる。
　　　c^(q,s)_j = c_j + c_j+sω^qn + …… + c_j+(m-1)s ω^(m-1)qn
　　　　　　　　　　　　j=0,1, …… , s-1
　FMTによる4分割乗算(m-4)の例を以下に示す。
　m=4(s=n/2)でqを0,1/4,1/2及び3/4として上式に代入し整理すると得られる。
　　　c_j　　　= [ ( c^(0,n/2)_j + c^(1/2,n/2)_j ) + ( c^(1/4,n/2)_j + c^(3/4,n/2)_j ) ] / 4
　　　c_j+n/2 = [ ( c^(0,n/2)_j - c^(1/2,n/2)_j ) + ( c^(1/4,n/2)_j - c^(3/4,n/2)_j ) ･ ω^-n/4 ] / 4
　　　c_j+n　　= [ ( c^(0,n/2)_j + c^(1/2,n/2)_j ) - ( c^(1/4,n/2)_j + c^(3/4,n/2)_j ) ] / 4
　　　c_j+3n/2 = [ ( c^(0,n/2)_j - c^(1/2,n/2)_j ) - ( c^(1/4,n/2)_j - c^(3/4,n/2)_j ) ･ ω^-n/4 ] / 4