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

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

3. 基本式

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

3.1 Eⁿ-1の因数分解

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

3.2 順変換

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

3.3 逆変換

　　k=0,1,...,n-1なる下記のn個の関数値
　　　f_k = mod( f, E - ω^k )
　から元の記号Eに関するn-1次の多項式関数
　　　　f = a_n-1E^n-1 + …… + a₁E + a₀
　のn個の係数a_jを求めるのがFMTの逆変換である。
　　n個の係数a_jは下記で計算される。
　　　a_j = ( f_n-1ω^{-(n-1) j} + …… + f₁ω^-j + f₀ ) / n　　　　　　　　　　　　　- - - (4)
　　　　　　j = 0, 1, …… , n-1　
　　(4)式はωが1のFFTは原始n乗根の性質である(1)式から容易に導出できる。

3.4 FMTとFFTの関係

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

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

　　記号Eに関するn-1次の2個の関数f,gから畳込み演算で2n-1次の関数hの各係数c_kを
　求める方法を述べる。
　　　f = a_n-1E^n-1 + …… + a₁E + a₀ , 　　g = b_n-1E^n-1 + …… + b₁E + b₀　　- - - (5)
　　　h = f ･ g = c_2n-1E^2n-1 + c_2n-2E^2n-2 + …… + c₁E + a₀　　　　　　　　　　- - - (6)
　関数hの各係数c_kは下記3種類の変換を順に行って求める。
　(1) FMT順変換
　　　(5)式のf,gにFMT変換に原始2n乗根ωを用いて(3)式に従いFMT順変換を行い、下記の
　　各2n個の係数f_k,g_kを求める。
　　　　f_k = mod( f, E - ω^k) = a_2n-1ω^{(2n-1) k} + a_2n-2ω^{(2n-2) k} + …… + a₁ω^k + a₀
　　　　g_k = mod( g, E - ω^k) = b_2n-1ω^{(2n-1) k} + b_2n-2ω^{(2n-2) k} + …… + b₁ω^k + b₀　
　　　　　　k = 0, 1, …… , 2n-1　　　　　　　　　　　　　　　　　　　　　　　　　　- - - (7)
　(2) 項別乗算
　　　(7)式のf_k,g_kを対応する項別に乗算して、下記のh_k を求める。
　　　　h_k = f_k ･ g_k 　　　　　　k = 0, 1, …… , 2n-1　　　　　　　　　　　　　　　- - - (8)
　　ここで、h_k = mod(h , E - ω^k )なる関係がある。
　(3) FMT逆変換
　　　h_k = mod(h , E - ω^k )なる関係を利用して、FMT逆変換で(8)式で求めたh_k から
　　Eに関する2n-1次の関数hの2n個の係数c_jを(4)式に従い下記で求める。
　　　c_j = ( f_2n-1ω^{-(2n-1) j} + f_2n-2ω^{-(2n-2) j} + …… + f₁ω^-j + f₀ ) / ( 2n )
　　　　　　j = 0, 1, …… , 2n-1　　　　　　　　　　　　　　　　　　　　　　　　　　- - - (9)

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

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

5.1 関係式

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

5.2 計算方法

　　(10)式のFMTの計算は下記のように通常のFMTに変更を加えて求める2種類の方法がある。
　これは、Eⁿ-ω^qn及びEⁿ+ω^qnが下記のように原始n乗根ωを使用して分解できるためである。
　　　( Eⁿ - ω^qn ) = (E - ω^n-1+q)(E - ω^n-2+q) …… (E - ω^1+q)(E - ω^q)
　　　( Eⁿ + ω^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,Eⁿ-1))のテーブル： ( ω^n-1, ω^n-2, …… , ω¹, ω⁰ )
　　　　mod(h,Eⁿ-ω^qn)のωのテーブル： ( ω^n-1+q, ω^n-2+q, …… , ω^1+q, ω^q )
　　　　mod(h,Eⁿ+ω^qn)のωのテーブル： ( ω^n-1/2+q, ω^n-3/2+q, …… , ω^3/2+q, ω^1/2+q )
　　
　(2) 入力係数及び出力係数に変換係数を掛ける
　　入力係数a_j出力係数とc_jに下記を施して通常のFMT計算を適用するとmod変換されたFMT
　　の結果が得られる。
　　　　入力係数変換： a^(q)_j = a_j ω^qj , 　　　 a^(q+1/2)_j = a_j ω^qj+1/2
　　　　出力係数変換： c^(q)_j = c_j ω^-qj/n , 　　　 c^(q+1/2)_j = c_j ω^-(qj+1/2)
　　　　　　　q ：有理数 ,　　　j = 0, 1, …… , n-1　　　