FMT(高速剰余変換）の原理

　　　　　　　　2002/12/20 日立製作所 & 早稲田大学　後 保範 ( Ushiro Yasunori )
　　　　　　　　　2001/6/9　作成のものをweb化
　　　　　　　　2003/05/20 更新（5章)　後 保範

1. 概要と目的

　FMT(Fast Modulo Transform, 高速剰余変換)はFFT(高速フーリエ変換)の原理に剰余(mod)
の拡張を追加して作成したものである。FFTとFMTの計算式は同一であるが、与える係数の
順が異なる。FMTは順変換の結果が剰余として意味つけられて多数桁計算への応用ではFFT
より見通しが良い。
　円周率世界記録 の計算に使用した多数桁乗算はFFTで構成していると記載したが、これは
まだFMTの理論を確立する前にプログラムを作成したためである。このため、円周率世界記録
の計算に使用した多数桁乗算は、与える係数順がFFTのものを使用したが、整数上で定義した
ものを使用しており、本FMTの原理の応用と言える。
　ここでは、計算原理を説明する。原理の概観は FMTを、 応用の説明は FMTの応用を参照。
FMTを利用したソースプログラムは多数桁計算プログラム、 π計算プログラムを参照。

2. FMTとFFT及び百五減算の関係

　n個の素数P1,P2,...,Pnを使用して、各剰余を求め、 各剰余から元の値を復元するのが百五減算
(3,5,7の剰余から105の剰余が求まることに由来、中国剰余定理)である。
ここでは、剰余modを整数だけでなく、一般の記号にも拡張し、記号Eと数ω及び整数kに対して
mod(f,E-ωk)をEの多項式f中のEをωkに置き換える処理と定義する。 一方ωをFFTの基本原理である、
１の原始n乗根とし、素数Pkの代わりにE-ωk, k=1,2,...,nを利用して、 畳み込み演算に利用可能な
アルゴリズムを構築する。このアルゴリズムをFMTと呼ぶ。
　FMTの順変換が百五減算の剰余を求めることに、逆変換が剰余から元の値に復元する処理に対応
する。この相似関係を利用して、百五減算で利用可能なアルゴリズム はすべて、FMTに適応できる。
FFTではωは複素数でω=e-2πi/nと定義する。これに対して、 FMTではωを複素数以外にも拡張して、
mod(ωn,P)=1となる整数n,ω,Pの組でも良い。見かけの数式上はFFTとFMTは同じ表現となる。
但し、FFTとFMTでは変換に利用するωの配列がω0,ω1,...,ωn-1 の順からωn-1,...,ω1,ω0と逆転する。
これは、入力する要素位置を逆にしたのと同じである。多数桁の乗算ではFMTの方がFFTより
応用上便利である。

3. 基本式

　FMTは拡張mod関数の上で定義し、En-1がn個に分解できる原理を活用する。
　ここで、nは説明を簡単にするため偶数とする。

3.1 En-1の因数分解

　ωが1のFFTは原始n乗根 の性質を利用すると、En-1は下記のように因数分解できる。
　　　En-1 = (E - ω0)(E - ω1) …… (E - ωn-1)
　理由：ωが１の原始n乗根で下記の関係が成立するため。
　　　ω0 + ωj + ω2j + …… + ω(n-1) j = 0 , 　　j=1, …… ,n-1
　 　　　 　　　 　　　 　　　 　　　 　　= n ,　　j=0 　　　　　　　　　　 　　　- - - (1)
　　　ω0 ･ ω1 ･ ω2 ･ …… ･ ω(n-1) = ωn(n-1)/2 = ω-1/2 = -1　　　　　- - - (2)

3.2 順変換

　　記号Eに関するn-1次の下記の多項式関数
　　　f = an-1En-1 + …… + a1E + a0
　に対して、
　　　fk = mod( f, E - ωk )
　となるfk, k=0,1,...,n-1 を求めるのがFMTの順変換である。
　　ここで、mod( f, E - ωk )は記号Eのベキと係数akで表わされる 多項式fの記号を
　数値ωkで置きかえることで定義する。
　順変換結果のfkは下記のようになる。
　　　fk = an-1ω(n-1) k + …… + a1ωk + a0　　　　　　　　　　　　　　　　　- - - (3)
　　　　　　k = 0, 1, …… , n-1　　

3.3 逆変換

　　k=0,1,...,n-1なる下記のn個の関数値
　　　fk = mod( f, E - ωk )
　から元の記号Eに関するn-1次の多項式関数
　 　　　f = an-1En-1 + …… + a1E + a0
　のn個の係数ajを求めるのがFMTの逆変換である。
　　n個の係数ajは下記で計算される。
　　　aj = ( fn-1ω-(n-1) j + …… + f1ω-j + f0 ) / n　　　　　　　　　　　　　- - - (4)
　　　　　　j = 0, 1, …… , n-1　
　　(4)式はωが1のFFTは原始n乗根の性質である(1)式から容易に導出できる。

3.4 FMTとFFTの関係

　　(3)式及び(4)式からFMTとFFTはωを複素数で定義すると全く同一と見える。
　ただし、多数桁の乗算に利用する観点から見ると、係数akの順序が異なる。
　FFTでは下位桁から順に与えるのに対して、FMTでは上位桁から順に与える。これは、
　実際の計算でFFTのωに対するテーブルが逆順になることを意味する。また、FMTとFFT
　では、乗算結果の桁が１要素分異なる。また、負巡回の場合の符号が異なる。
　このため、多数桁の乗算を考えるときはFFTよりFMTの方が桁の位置や負巡回等の符号を
　含めて便利である。
　　FMTではωを複素数(共役関係を利用した実数を含む)だけでなく、整数上で定義する
　ものへ拡張される。また、記号上の剰余定理(mod関数)により、FMTの途中経過が全て、
　どの位置の剰余に対応するか分かるため、百五減算が適用される高精度解法にそのまま
　適用できる。

4. 畳込み演算への適用方法

　　記号Eに関するn-1次の2個の関数f,gから畳込み演算で2n-1次の関数hの各係数ckを
　求める方法を述べる。
　　　f = an-1En-1 + …… + a1E + a0 , 　　g = bn-1En-1 + …… + b1E + b0　　- - - (5)
　　　h = f ･ g = c2n-1E2n-1 + c2n-2E2n-2 + …… + c1E + a0　　　　　　　　　　- - - (6)
　関数hの各係数ckは下記3種類の変換を順に行って求める。
　(1) FMT順変換
　　　(5)式のf,gにFMT変換に原始2n乗根ωを用いて(3)式に従いFMT順変換を行い、下記の
　　各2n個の係数fk,gkを求める。
　　　　fk = mod( f, E - ωk) = a2n-1ω(2n-1) k + a2n-2ω(2n-2) k + …… + a1ωk + a0
　　　　gk = mod( g, E - ωk) = b2n-1ω(2n-1) k + b2n-2ω(2n-2) k + …… + b1ωk + b0　
　　　　　　k = 0, 1, …… , 2n-1　　　　　　　　 　　　　　　　　　　　　　　　　　　- - - (7)
　(2) 項別乗算
　　　(7)式のfk,gkを対応する項別に乗算して、下記のhk を求める。
　　　　hk = fk ･ gk 　　　　　　k = 0, 1, …… , 2n-1　　　　　　　　　　　　　　　- - - (8)
　　ここで、hk = mod(h , E - ωk )なる関係がある。
　(3) FMT逆変換
　 　　hk = mod(h , E - ωk )なる関係を利用して、FMT逆変換で(8)式で 求めたhk から
　　Eに関する2n-1次の関数hの2n個の係数cjを(4)式に従い下記で求める。
　　　cj = ( f2n-1ω-(2n-1) j + f2n-2ω-(2n-2) j + …… + f1ω-j + f0 ) / ( 2n )
　　　　　　j = 0, 1, …… , 2n-1　　　　　　　　 　　　　　　　　　　　　　　　　　　- - - (9)

5. 畳込み演算とFMTとの関係式

　　(5)式に示す記号Eに関するn-1次の2個の関数f,gの畳込み演算結果は(6)式に示す2n-1次
　の関数hとなる。これに対して、適当な有理数qと mod(h,En-ωqn)を導入し結果も入力と同じn-1次
　の関数にすることを考える。このときの関係式と計算方法を示す。

5.1 関係式

　(6)式hに対して、適当なqに対してmod(h,En-ωqn)を行ったものh(q) と表すと次のようになる。
　　　h(q)　　= mod ( h, En - ωqn ) = c(q)n-1 E n-1 + …… + c(q)1 E + c(q)0
　　　h(q+1/2) = mod ( h, En + ωqn ) = c(q+1/2)n-1 E n-1 + …… + c(q+1/2)1 E + c(q+1/2)0
　　　　　q : 有理数　　　　　　　　　　　 　　　　　　　　　　　　　　　　- - - (10)
　 (5)式の畳込み演算hの係数cjと(10)式の係数c(q)jは下記の 関係を持つ。
　　　c(q)j　　　= cj + cj+n ω qn
　　　c(q+1/2)j = cj - cj+n ω qn
　　　　　q : 有理数 ,　　　j = 0, 1, …… , n-1　　　　　　　　　　　　- - - (11)
　逆にcj, cj+nはc(q)jを使用して下記の様に 計算できる。
　　　cj　　= ( c(q)j + c(q+1/2)j ) / 2
　　　cj+n = ( c(q)j - c(q+1/2)j ) ω -qn / 2
　　　　　q : 有理数 ,　　　j = 0, 1, …… , n-1　　　　　　　　　　　　- - - (12)

5.2 計算方法

　　(10)式のFMTの計算は下記のように通常のFMTに変更を加えて求める2種類の方法がある。
　これは、En-ωqn及びEn+ωqnが下記のように 原始n乗根ωを使用して分解できるためである。
　　　( En - ωqn ) = (E - ωn-1+q)(E - ωn-2+q) …… (E - ω1+q)(E - ωq)
　　　( En + ωqn ) = (E - ωn-1/2+q)(E - ωn-3/2+q) …… (E - ω3/2+q)(E - ω1/2+q)

　(1) 原始n乗根ωのテーブルを変更する方法
　　下記の様に通常のFMTのωのテーブルをmod付きωのテーブルに変更する。
　　逆変換を行うには、テーブル上のωをω-1に変更する。
　　ωのテーブルの与え方がFMTとFFTでは逆順になることに注意を要す。
　　　　通常(mod(h,En-1))のテーブル ： ( ωn-1, ωn-2, …… , ω1, ω0 )
　　　　mod(h,En-ωqn)のωのテーブル ： ( ωn-1+q, ωn-2+q, …… , ω1+q, ωq )
　　　　mod(h,En+ωqn)のωのテーブル ： ( ωn-1/2+q, ωn-3/2+q, …… , ω3/2+q, ω1/2+q )
　　
　(2) 入力係数及び出力係数に変換係数を掛ける
　　入力係数aj出力係数とcjに下記を施して通常のFMT計算を適用するとmod変換されたFMT
　　の結果が得られる。
　　　　入力係数変換： a(q)j = aj ωqj , 　 　　 a(q+1/2)j = aj ωqj+1/2
　　　　出力係数変換： c(q)j = cj ω-qj/n , 　　　 c(q+1/2)j = cj ω-(qj+1/2)
　　　　　　　q ： 有理数 ,　　　j = 0, 1, …… , n-1