数値演算法

(4) 高速フーリエ変換

乗算処理を高速化するための手法として、前回は Karatsuba 法と Toom-Cook 法を紹介しました。どちらの方法も、数値を多項式として表現してその各係数を求めることによって、乗算結果を得ることができるということを利用しています。今回は、多項式を利用した手法としてまずは一般的な解法を考え、最終的には高速フーリエ変換を利用した乗算処理の高速化について説明したいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) 連立方程式を使った解法

Toom-Cook 法で用いた手段をもう一度おさらいします。
N 桁の二数 u, v をそれぞれ N 次の多項式 u(x) = Σi{0→N}( uixi ), v(x) = Σi{0→N}( vixi ) で表し、その積を w(x) = Σi{0→2N}( wixi ) としたとき、u(x) と v(x) に任意の値 xi を代入してその積 w(xi) を求める操作を 2N 回繰り返すことにより、係数 wi を未知数とする 2N 個の連立方程式を作ることができます。よって、連立方程式を解くことによって wi を求め、それによって得られた多項式 w(x) に基数を代入して二数の積を求めることができます。Toom-Cook 法は、xi が 0 から 2N までの自然数である場合に利用できるアルゴリズムですが、xi として任意の数を利用しても、連立方程式を解くことができれば wi を求めることは可能です。例えば、1251 と 2211 の積を計算するとき、これらを多項式と見なして次のように計算することが可能です。

u(x)=12x + 51
v(x)=22x + 11
w(x)=u(x)v(x)
=Σi{0→2}( wixi ) として
u(-1) = 39u(0) = 51u(1) = 63
v(-1) = -11v(0) = 11v(1) = 33
w(-1) = -429w(0) = 561w(1) = 2079

よって、連立方程式は

w2 - w1 + w0=-429
w0=561
w2 + w1 + w0=2079

となって、w0 = 561、w1 = 1254、w2 = 264 が得られ、w(x) = 264x2 + 1254x + 561 より w(100) = 2640000 + 125400 + 561 = 2765961 と求めることができます。

連立方程式は、行列を利用して解くことができます。そのことを説明する前に、行列に関して知っておかなければならない内容を整理しておきます。


N 行 N 列の正方行列 ( 行列数が等しい行列 ) A に対して、第 r 行と第 c 列を除いた N - 1 行 N - 1 列の行列の行列式を、A の第 r 行と第 c 列成分 arc に関する「余因子 ( Cofactor )」といいます。行列 A の行列式を detA、arc に関する余因子を cofArc としたとき、以下の公式が成り立ちます。

detA = Σi{1→n}( ( -1 )r+i・ari・cofAri )

detA = Σi{1→n}( ( -1 )i+c・aic・cofAic )

上側の式を第 r 行に関する「余因子展開」、下側の式を第 c 列に関する余因子展開といいます。余因子展開を利用することで、任意の大きさの正方行列に対する行列式を計算することができます。例えば、4 行 4 列の行列に対する行列式は以下のように計算できます。

第 4 列に対して余因子展開

|4,3,2,1|
|3,2,1,0|
|2,1,0,0|
det|1,0,0,0|
 
|3,2,1||4,3,2||4,3,2||4,3,2|
|2,1,0||2,1,0||3,2,1||3,2,1|
=-1 x det|1,0,0| + 0 x det|1,0,0| - 0 x det|1,0,0| + 0 x det|2,1,0|
 
|3,2,1|
|2,1,0|
=-det|1,0,0|

3列に対して余因子展開

|2,1||3,2||3,2|
=-1 x det|1,0| + 0 x det|1,0| - 0 x det|2,1|
 
|2,1|
=-1 x det|1,0|
= -( 2 x 0 - 1 x 1 ) = 1

なお、最後にある 2 行 2 列の行列式の計算部分は、以下の公式を利用しています。

   | a, b |
det| c, d | = ad - bc

行列 A に対し、a~rc = (-1)r+c・cofAcr を第 r 行第 c 列の要素に持つ、A と同じ大きさの行列を、行列 A の「余因子行列 ( Adjugate Matrix )」と言います。ここで、cofAcr は行列 A の第 c 行と第 r 列の成分 acr ( arc ではない ) に関する余因子であることに注意してください。

この時、以下の関係式が成り立ちます。
Σi{1→n}( ari・a~ic ) = detA・δrc

上式にある δrc は「クロネッカのデルタ (Kronecker Delta)」と呼ばれ、r = c のとき 1、それ以外は 0 になります。

上記関係式から、行列 A の余因子行列を A~ としたとき、以下の関係式が成り立ちます。
A A~ = detAE

ここで E は「単位行列 ( Identity Matrix )」であり、各要素を δrc で表すことができます ( 対角要素が 1 で、他の要素が 0 になります )。

行列 A に対し、AX = E を満たす行列 X を、行列 A の「逆行列 (Inverse Matrix)」といい、A-1 と表します。上式から、detA ≠ 0 ならば、逆行列は次のように表せることになります ( 補足 1 )。
A-1 = A~ / detA

例として前に示した連立方程式を行列を使って表すと、

|-429|=|1,-1,1||w2|
|561||0,0,1||w1|
|2079||1,1,1||w0|

になります ( 右辺を計算すると、連立方程式の形に戻るのが分かると思います )。
右辺の行列に対する逆行列が分かっているとき、両辺にその逆行列を左側から掛ければ、右辺の行列は単位行列になります。

|1,-1,1| -1|-1,2,-1|
|0,0,1| = -(1/2)|1,0,-1|
|1,1,1||0,-2,0|

より、

|-1,2,-1||-429||1,0,0||w2||w2|
-(1/2)|1,0,-1||561| = |0,1,0||w1| = |w1|
|0,-2,0||2079||0,0,1||w0||w0|

となるため、左辺を計算することで各係数を求めることができます。
以上の内容から、余因子展開によって行列式を求め、それを利用して余因子行列を求めれば、任意の連立方程式を解いて乗算処理を行うのが可能であることが分かります。

上記説明の中では、行列式の定義そのものやその値が持つ意味については触れていません。ここでは、n 個の未知数を持つ連立一次方程式を行列を使って表したとき、この連立方程式が唯一の解を持つための必要十分条件は、行列式の値が 0 でないことであることを知っておけば十分です。なお、行列式の値が 0 である場合、その連立方程式は解が存在しないか(不能)、解が無数に存在する(不定)状態になります。

まずは、行列演算用のサンプル・プログラムを示します。

/*
  Determinant : 正方行列 mat の行列式を求める
*/
template< class T >
T Determinant( const SquareMatrix< T >& mat )
{
  if ( mat.size() == 1 ) return( mat[0][0] );
  if ( mat.size() == 2 )
    return( mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0] );

  T ans = 0;
  for ( typename SquareMatrix< T >::size_type i = 0 ; i < mat.size() ; ++i ) {
    T num = mat[i][0];
    num *= ( ( ( i & 1 ) == 0 ) ? 1 : -1 );
    if ( num == 0 ) continue;
    SquareMatrix< T > minor = MinorMatrix( mat, i, 0 );
    if ( minor.size() > 0 ) {
      T detVal = Determinant( minor );
      num *= detVal;
    }
    ans += num;
  }

  return( ans );
}

/*
  MinorMatrix : 行列 mat から 行 r と列 c を除いた小行列を返す
*/
template< class T >
SquareMatrix< T > MinorMatrix( const SquareMatrix< T >& mat, typename SquareMatrix< T >::size_type r, typename SquareMatrix< T >::size_type c )
{
  if ( mat.size() == 0 ) return( mat );
  SquareMatrix< T > minorMat( mat.size() - 1 );

  typename SquareMatrix< T >::size_type r1 = 0;
  for ( typename SquareMatrix< T >::size_type r0 = 0 ; r0 < mat.size() ; ++r0 ) {
    if ( r0 == r ) continue;
    typename SquareMatrix< T >::size_type c1 = 0;
    for ( typename SquareMatrix< T >::size_type c0 = 0 ; c0 < mat.size() ; ++c0 ) {
      if ( c0 == c ) continue;
      minorMat[r1][c1] = mat[r0][c0];
      ++c1;
    }
    ++r1;
  }

  return( minorMat );
}

/*
  Cofactor : 正方行列 mat の行 r 列 c の要素に関する余因子を返す
*/
template< class T >
T Cofactor( const SquareMatrix< T >& mat, typename SquareMatrix< T >::size_type r, typename SquareMatrix< T >::size_type c )
{
  if ( mat.size() == 1 ) return( 1 );

  return( Determinant( MinorMatrix( mat, r, c ) ) );
}

/*
  AdjugateMatrix : 正方行列 mat の余因子行列を返す(余因子行列の各要素を行列式で割ると逆行列になる)
*/
template< class T >
SquareMatrix< T > AdjugateMatrix( const SquareMatrix< T >& mat )
{
  SquareMatrix< T > adj( mat.size() );

  for ( typename SquareMatrix< T >::size_type r = 0 ; r < mat.size() ; ++r ) {
    for ( typename SquareMatrix< T >::size_type c = 0 ; c < mat.size() ; ++c ) {
      adj[c][r] = Cofactor( mat, r, c );
      adj[c][r] *= ( ( ( ( r + c ) & 1 ) == 0 ) ? 1 : -1 );
    }
  }

  return( adj );
}

SquareMatrix は正方行列を表すクラスですが、この中では実装されていません。行列の各要素は二次元配列と同様の方法でアクセス方法が可能であり、メンバ関数としては、行列数 ( サイズ ) を表す size があります。また、サイズを表す型として size_type が定義されています。
行列式は Determinant で計算しています。また、行列式の計算で必要となる、指定した行と列を除いた小行列は MinorMatrix にて作成します。余因子の計算は Cofactor が行い、余因子行列を求める AdjugateMatrixCofactor を順番に呼び出しているだけの関数です。なお、逆行列を求めるメンバ関数が存在しないのは、逆行列を求める場合に行列式による除算処理を含むため、各要素が整数とならない可能性があるためです。その他の関数では除算処理がないため、行列の要素が整数値である場合は必ず処理結果も整数値となります。

行列を使った計算では値が負数になる場合があります。前回までは、正数用の多倍長整数までしか用意されていなかったので、ここで符号を含めた多倍長整数クラスを用意しておきます。

namespace BigNum
{
  class Signed
  {
    Unsigned number_;
    bool isPos_; // 符号(正,0=true,負=false)

    // t による初期化
    template< class T >
      void init( T t );

   public:

    // コンストラクタ
    Signed( int i = 0 );
    explicit Signed( const Unsigned& num )
      : number_( num ), isPos_( true ) {}

    bool div( const Signed& v, Unsigned* r );
    bool isNaN() const { return( number_.isNaN() ); }

    bool operator!() const { return( ! number_ ); }

    /* 比較演算子 */

    bool operator==( const Signed& num ) const;
    bool operator<( const Signed& num ) const;
    bool operator>( const Signed& num ) const;
    bool operator!=( const Signed& num ) const;
    bool operator<=( const Signed& num ) const;
    bool operator>=( const Signed& num ) const;

    /* 計算ルーチン */

    Signed& operator+=( const Signed& num );
    Signed& operator-=( const Signed& num );
    Signed& operator*=( const Signed& num );
    Signed& operator/=( const Signed& num );
    Signed& operator%=( const Signed& num );

    Signed& operator>>=( UnsignedBuff num ) { number_ >>= num; if ( number_.isZero() ) isPos_ = true; return( *this ); }
    Signed& operator<<=( UnsignedBuff num ) { number_ <<= num; return( *this ); }

    Unsigned toUnsigned() const { return( number_ ); }
  };

  /* 演算子の多重定義 */

  Signed operator+( const Signed& n1, const Signed& n2 );
  Signed operator-( const Signed& n1, const Signed& n2 );
  Signed operator*( const Signed& n1, const Signed& n2 );
  Signed operator/( const Signed& n1, const Signed& n2 );
  Signed operator%( const Signed& n1, const Signed& n2 );
}

演算処理には、BigNum::Unsigned 用のルーチンをそのまま流用することができます。乗算・除算処理は符号が同一ならば結果を正、そうでなければ結果を負とするだけで対応が可能です。加算・減算処理では、符号をチェックして適切な処理を行うようにします。例えば加算処理の場合、符号が異なっていれば減算処理に切り替える必要があります ( ここでは具体的な実装は示していません )。ビットシフトは符号ビットの変更はなく、数の部分だけをシフトします。但し、右ビットシフトでは演算結果がゼロになる可能性があるので、ゼロであれば符号を正にする処理を加えています。

以上のルーチンが揃えば、逆行列を使った乗算処理を行うことができます。

/*
  CreateMatrix : a[r][c] = pow( r, c - 1 ) となる行列を mat に作成する

  | 1 0 0 0 ... |
  | 1 1 1 1 ... |
  | 1 2 4 8 ... |
  | : : : :     |
  |             |
*/
template< class T >
void CreateMatrix( SquareMatrix< T >& mat )
{
  for ( typename SquareMatrix< T >::size_type r = 0 ; r < mat.size() ; r++ ) {
    mat[r][0] = T( 1UL );
    T num( r );
    for ( typename SquareMatrix< T >::size_type c = 1 ; c < mat.size() ; c++ ) {
      mat[r][c] = num;
      num *= r;
    }
  }
}

/*
  MulMatVec : 行列 mat とベクトル vec の乗算処理

  要素の大きさが一致しない場合、ベクトルの大きい要素は無視され、小さい場合は末尾を 0 とする
*/
template< class T, class U >
vector< T > MulMatVec( const SquareMatrix< T >& mat, const vector< U >& vec )
{
  vector< T > ans; // 演算結果

  for ( typename SquareMatrix< T >::size_type r = 0 ; r < mat.size() ; r++ ) {
    ans.push_back( T() );
    for ( typename SquareMatrix< T >::size_type c = 0 ; c < mat.size() ; c++ ) {
      if ( c >= vec.size() ) break;
      ans.back() += mat[r][c] * U( vec[c] );
    }
  }

  return( ans );
}

namespace BigNum
{
  /*
    Unsigned::matrixMul : 逆行列を使った乗算処理

    const UBigNum& v : 乗数
  */
  Unsigned& Unsigned::matrixMul( const Unsigned& v )
  {
    if ( isNaN( v ) ) return( *this );

    // 結果のベクトルのサイズを計算
    size_t sz = ( ( size() > v.size() ) ? size() : v.size() - 1 ) * 2 + 1;

    // 連立方程式の行列表現
    SquareMatrix< Signed > mat( sz );
    CreateMatrix( mat );

    // u(i), v(i)を求める
    vector< Signed > uAns = MulMatVec( mat, number_ );
    vector< Signed > vAns = MulMatVec( mat, v.number_ );

    // w(i) = u(i) * v(i) を求める(結果はuAnsに格納)
    for ( size_t i = 0 ; i < uAns.size() ; i++ )
      uAns[i] *= vAns[i];

    SquareMatrix< Signed > adjMat = AdjugateMatrix( mat ); // 余因子行列

    // 余因子行列との積を求めて積の各係数w(i)を求める
    vector< Signed > wVec = MulMatVec( adjMat, uAns );

    // 係数を数値に変換
    clear();
    for ( unsigned int i = wVec.size() ; i > 0 ; --i )
      partialAdd( wVec[i-1].toUnsigned(), i - 1 );

    operator/=( Determinant( mat ).toUnsigned() ); // 行列式による除算は最後に

    return( *this );
  }
}

CreateMatrixでは、0 から n - 1 までの値を多項式に代入してできる連立方程式の、行列による表現を作成するために使います。この部分は、他の値を利用しても ( 例えば -n / 2 から n / 2 まで ) 成り立ちますが、異なる行に対して同じ値を利用した場合、行列式が 0 となり解が求められなくなります。
partialAddは、多倍長整数として表現されている数列に対して、途中の桁から加算処理をするためのメンバ関数です。この具体的な実装は「数値演算法 (2) 多倍長整数の演算」のサンプル・プログラムの中に示してあります。

このように、逆行列を利用して連立方程式を機械的に解き、求めた係数から積を計算することができることができるわけですが、問題なのは、逆行列を解く操作自体が非常に重い処理であるということです。余因子行列を求めるときに各要素ごとで行列式を計算する必要があり、その計算量は莫大なものになります。そもそも逆行列や行列式を計算するときは、余因子行列を計算するのではなく、連立方程式を解いてその結果から導くのが通常の方法です。連立方程式を機械的に解く手法としては、「ガウスの消去法」や「ガウス・ジョルダンの消去法」、「LU分解」があります。しかし、これらも各要素の乗算回数は少なくはありません ( 「ガウスの消去法」や「ガウス・ジョルダンの消去法」を使った場合、乗算回数は桁数の三乗にほぼ比例します )。


2) 離散フーリエ変換

少なくとも逆行列が既知のものであれば、余因子行列を求めたり連立方程式を解くことなく、機械的に係数を求めることができます。そのような便利なものに「離散フーリエ変換 ( Discrete Fourier Transform ; DFT )」があります。

離散データ { fl } から周波数成分 { Fk } への離散フーリエ変換式は以下のようになります。

Fk = (1/N)Σl{0→N-1}( fle-i2πkl/N )

また、逆変換は次のようになります。

fl = Σk{0→N-1}( Fkei2πkl/N )

{ fl } の添字 l ( = 0, 1, ..., N ) をデータ番号、{ Fk } の添字 k ( = 0, 1, ..., N ) を周波数番号と言います。

ここで、

ωN = ei2π/N

とすると、上式は次のように表すことができます。

fl = Σk{0→N-1}( FkωNkl )
Fk = (1/N)Σl{0→N-1}( flωN-kl )

さらに、行列とベクトルを使って表すと、次のようになります。

|f0||1,1,1,...1||F0|
|f1||1,ωN,ωN2,...ωNN-1||F1|
|f2| = |1,ωN2,ωN4,...ωN2(N-1)||F2|
|:||::::||:|
|fN-1||1,ωNN-1,ωN2(N-1),...ωN(N-1)(N-1)||FN-1|
|F0||1,1,1,...1||f0|
|F1||1,ωN-1,ωN-2,...ωN-(N-1)||f1|
|F2| = (1/N)|1,ωN-2,ωN-4,...ωN-2(N-1)||f2|
|:||::::||:|
|FN-1||1,ωN-(N-1),ωN-2(N-1),...ωN-(N-1)(N-1)||fN-1|

この関係式を見ると、上側の式にある行列の逆行列が下側の行列であることは明らかです。また、

fl = Σk{0→N-1}( FkωNkl ) = Σk{0→N-1}( FkNl)k )

より、

f(x) = Σk{0→N-1}( Fkxk ) = F0 + F1x + F2x2 + ... + FN-1xN-1

としたとき fl = f(ωNl) であり、これは Fk を未知数とする連立方程式に対して、ωNl を代入して未知数を求める解法そのものになります。これなら、逆行列を計算することなく機械的に係数を求めることができます。

離散フーリエ変換・逆変換を使って、実際に数値を計算してみます。最初の章で例に示した二数、1251 と 2211 をもう一度使います。

u(x)=12x + 51
v(x)=22x + 11
w(x)=u(x)v(x)
=Σi{0→2}( wixi ) として
u(1)=63u(ω)=12ω + 51u(ω2)=12ω2 + 51
v(1)=33v(ω)=22ω + 11v(ω2)=22ω2 + 11
w(1)=2079w(ω)=264ω2 + 1254ω + 561w(ω2)=264ω4 + 1254ω2 + 561
=1254ω2 + 264ω + 561

上式で、ω は ω3 = ei2π/3 を表しています。

フーリエ変換式を使って、

|w0|= (1/3)|1,1,1||2079|
|w1||1,ω-1,ω-2||264ω2 + 1254ω + 561|
|w2||1,ω-2,ω-4||1254ω2 + 264ω + 561|

より、

w0=(1/3)( 2079 + ( 264ω2 + 1254ω + 561 ) + ( 1254ω2 + 264ω + 561 ) )
=506ω2 + 506ω + 1067
w1=(1/3)( 2079 + ( 264ω2 + 1254ω + 561 )ω-1 + ( 1254ω2 + 264ω + 561 )ω-2 )
=275ω2 + 275ω + 1529
w2=(1/3)( 2079 + ( 264ω2 + 1254ω + 561 )ω-2 + ( 1254ω2 + 264ω + 561 )ω-4 )
=605ω2 + 605ω + 869

困ったことに、計算結果は複素数を含んでいます。しかし、ω3 に対しては ω2 + ω + 1 = 0 が成り立つので ( 理由は次の章にて説明します )、さらに次のように変形することができます。

w0=506( ω2 + ω + 1 ) + 561 = 561
w1=275( ω2 + ω + 1 ) + 1254 = 1254
w2=605( ω2 + ω + 1 ) + 264 = 264

よって、前と同じ結果の w(x) = 264x2 + 1254x + 561 が得られました。

以下は、離散フーリエ変換・逆変換処理と、それを使った乗算処理用のサンプル・プログラムです。

/*
  IDft : 離散フーリエ逆変換

  N : データ数
  vec : 逆変換を行うベクトル
  ans : 結果を格納する変数へのポインタ
*/
void IDft( BigNum::UnsignedBuff n, const vector< Digit >& vec, vector< complex< long double > >* ans )
{
  for ( BigNum::UnsignedBuff i = 0 ; i < n ; i++ ) {
    ans->push_back( complex< long double >( 0, 0 ) );
    unsigned int k = 0;
    for ( vector< Digit >::const_iterator s = vec.begin() ; s != vec.end() ; ++s ) {
      long double re = cosl( k * 2 * M_PIl / n );
      long double im = sinl( k * 2 * M_PIl / n );
      (*ans)[i] += static_cast< long double >( *s ) * complex< long double >( re, im );
      k = ( k + i ) % n;
    }
  }
}

/*
  Dft : 離散フーリエ変換

  n : データ数
  vec : 変換を行うベクトル
  ans : 結果を格納する変数へのポインタ
*/
void Dft( BigNum::UnsignedBuff n, const vector< complex< long double > >& vec, vector< complex< long double > >* ans )
{
  for ( BigNum::UnsignedBuff i = 0 ; i < n ; i++ ) {
    ans->push_back( complex< long double >( 0, 0 ) );
    BigNum::UnsignedBuff k = 0;
    for ( vector< complex< long double > >::const_iterator s = vec.begin() ; s != vec.end() ; ++s ) {
      long double re = cosl( k * 2 * M_PIl / n );
      long double im = sinl( k * 2 * M_PIl / n );
      (*ans)[i] += *s * complex< long double >( re, im );
      k = ( k < i ) ? k + n - i : k - i;
    }
  }
}

namespace BigNum
{
  /*
    Unsigned::dftMul : 離散フーリエ変換を使った乗算処理

    v : 乗数
  */
  Unsigned& Unsigned::dftMul( const Unsigned& v )
  {
    UnsignedBuff i;
    if ( isNaN( v ) ) return( *this );

    // 分割数Nを決定
    UnsignedBuff n = ( ( ( size() > v.size() ) ? size() : v.size() ) - 1 ) * 2 + 1;

    vector< Digit > uVec = number_;   // u(i)
    vector< Digit > vVec = v.number_; // v(i)

    // 離散フーリエ逆変換で u(i), v(i)を求める
    vector< complex< long double > > uAns;
    IDft( n, uVec, &uAns );
    vector< complex< long double > > vAns;
    IDft( n, vVec, &vAns );

    // w(i) = u(i) * v(i) を求める(結果はuAnsに格納)
    for ( i = 0 ; i < uAns.size() ; i++ )
      uAns[i] *= vAns[i];

    // w(i)を離散フーリエ変換する
    vector< complex< long double > > wVec;
    Dft( n, uAns, &wVec );

    // 係数を数値に変換
    clear();
    for ( i = n ; i > 0 ; --i ) {
      if ( wVec[i - 1].real() < 1 ) continue;
      Unsigned n( wVec[i - 1].real() );
      partialAdd( n, i - 1 );
    }
    operator/=( n );

    return( *this );
  }
}

処理の流れは、最後の乗算処理ルーチン dftMul から順に見れば理解できると思います。乗数と被乗数をそれぞれ離散フーリエ逆変換により 1 から ωNN-1 までの解に変換し、両者を掛け合わせて結果に対する解を求めた後、それを使って離散フーリエ変換を行い係数を得ます。処理の中で逆行列を計算する部分が省略できた分、前のルーチンよりは早くなっています。
離散フーリエ変換・逆変換処理は DftIDftで 行っています。ここでは coslsinl を使って正弦・余弦を求め、その値を使って複素数を得ています。このあたりの内容については次の章で後述します。

離散フーリエ変換・逆変換の結果は複素数ですが、最終的に求められた係数は正の整数値となります ( 二つの多項式を乗算して得られる多項式の係数は、元の多項式が正の整数であれば必ず正の整数になります )。よって係数から数値を求めるとき、虚数部は相殺されると見なして実数部だけ取り出し、それを整数に変換しています。浮動小数点数から整数値への変換は、コンストラクタを別に用意することで対応しています。

以下に、浮動小数点数から BigNum::Unsigned 型に変換するコンストラクタを示します。

namespace BigNum
{
  // 浮動小数点数のビットフィールド宣言
  template< class D >
  struct FP_BitField;

  // double 用ビットフィールド
  template<>
  struct FP_BitField< double >
  { union ieee754_double d; };

  // long double 用ビットフィールド
  template<>
  struct FP_BitField< long double >
  { union ieee854_long_double d; };

  /*
    コンストラクタ(浮動小数点数 d を代入)

    D は浮動小数点数の型、eLen は符号ビット長、fLen は仮数ビット長をそれぞれ表す

    S ... 符号ビット
    E ... 指数ビット
    F ... 仮数ビット
    I ... Integer part (整数部)

    double 型
    SEEEEEEEEEEEFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF FFFFFFFFFFFFFFFF
     <- 11bit -><------------------ 52bit(6.5byte) ------------------->

    long double 型
    SEEEEEEEEEEEEEEE IFFFFFFFFFFFFFFF ... FFFFFFFFFFFFFFFF
     <--- 15bit --->  <------------- 63bit -------------->
  */
  template< class D, unsigned int eLen, unsigned int fLen >
  void Unsigned::initDbl( D d )
  {
    assert( sizeof( Digit ) * 2 <= sizeof( UnsignedBuff ) );

    const Digit eCen = ( 1 << ( eLen - 1 ) ) - 1; // 指数に対するバイアス値
    UnsignedBuff fByte = fLen / DIGIT_BIT; // 仮数部最上位の有効ビット数
    UnsignedBuff msFLen = fLen % DIGIT_BIT; // 仮数部最上位の有効ビット数
    if ( msFLen != 0 ) ++fByte;

    nan_ = false;

    // 特殊な値ならNaN扱い
    if ( isinf( d ) || isnan( d ) ) {
      push_upper( 0 );
      nan_ = true;
      return;
    }

    // 負数の場合はNaN扱い
    if ( d < 0 ) {
      push_upper( 0 );
      nan_ = true;
      return;
    }

    FP_BitField< D >* field = reinterpret_cast< FP_BitField< D >* >( &d );

    // 指数部の取得
    SignedBuff e = ( field->d ).ieee.exponent - eCen;

    if ( e < 0 ) {
      if ( e == -1 ) number_[0] = 1; // 0.5以上なので四捨五入して1とする
      return;
    }

    // 仮数部の取得
    unsigned int mantissa = ( field->d ).ieee.mantissa1;
    for ( size_t i = 0 ; i < sizeof( mantissa ) / sizeof( Digit ) ; ++i )
      push_upper( ( mantissa >> i * DIGIT_BIT ) & DIGIT_MAX );
    mantissa = ( field->d ).ieee.mantissa0;
    for ( size_t i = 0 ; i < sizeof( mantissa ) / sizeof( Digit ) ; ++i )
      push_upper( ( mantissa >> i * DIGIT_BIT ) & DIGIT_MAX );
    if ( msFLen != 0 ) {
      number_.back() |= 1 << msFLen; // 小数 + 1 の処理
    } else {
      push_upper( 1 );
    }

    // 桁合わせ
    if ( static_cast< UnsignedBuff >( e ) >= fLen ) {
      operator<<=( e - fLen );
    } else {
      operator>>=( fLen - e - 1 );
      if ( ( number_[0] & 1 ) != 0 ) operator+=( 2 ); // 小数部が0.5以上なら四捨五入
      operator>>=( 1 );
    }

    negZero();
  }

  // double 型から変換
  Unsigned::Unsigned( double d )
  { initDbl< double, 11, 52 >( d ); }

  // long double 型から変換
  Unsigned::Unsigned( long double d )
  { initDbl< long double, 15, 63 >( d ); }
}

浮動小数点数のビット構成は「IEEE 浮動小数点数演算基準 ( IEEE Standard for Floating-Point Arithmetic ; IEEE 754 )」という標準規格で定義されています。GNUC ライブラリでは ieee754.h というヘッダ・ファイルがあり、この中でビット・フィールドとして定義されています。指数部や仮数部のビットを取得するため、サンプル・プログラムではこのビット・フィールドを利用しています。
気を付けるべき点として、少数部の四捨五入をしている箇所があります。10 進数の 0.5 は、2 進数で表すと 0.1 となるため ( 1 を 2 で割ると 0.5 になります。これは 2 進数の場合 1 を 10 で割る操作と同等なので、0.1 と表されるわけです )、小数点第一位が 1 であれば、四捨五入により整数部に 1 を加えます。また、浮動小数点数は、指数部が 0 でない場合 1.nn... x 2N の形になるよう指数部を調整しており ( これを正規化といいます )、しかも整数部分の 1 は小数部には含まれず省略されています ( 例外もあります。これについては後述します )。よって、小数部を取得した後で、省略されたビットを必ず立てる必要があります。
最後の注意事項として、Intel CPU の場合、バイトが下位部分から順に並んでいます ( リトル・エンディアン Little Endian )。よって、配列から要素を取得する方向に気を付ける必要があります ( 補足ですが、配列の要素の中身まで気を付ける必要はありません。要素を取得する時は自動で並べ替えられています。さもないと、多バイトの要素をポインタや配列から取得するときにわざわざ並べ替えなければならなくなってしまいます )。 ieee754.h では CPU がリトル・エンディアンかビッグ・エンディアン ( Big Endian ) かも考慮してあり、サンプル・プログラムのように取得をすればどの場合でも正しく修得ができるようになっています。バイトの配置順については、バイト単位の他にワード単位で異なる場合もあるので、移植性を考慮するのであれば ieee754.h の利用はほぼ必須となります。

double 型と long double 型のビット構成は上記プログラムの中のコメントに記載してありますが、詳細な構成を以下に示しておきます。

■ 倍精度 ( double ) の記憶形式

図 2-1. double 型のビット構成
double形式
 
表 2-1. double 型のビット・パターン例
ビット・パターン(16進)数値(10進)
000000000 000000000.0
13FF00000 000000001.0
-1FFF00000 00000000-1.0
最大正規数7FEFFFFF FFFFFFFF1.7976931348623157e+308
正の最小正規数00100000 000000002.2250738585072014e-308
最大非正規数000FFFFF FFFFFFFF2.2250738585072009e-308
正の最小非正規数00000000 000000014.9406564584124654e-324
+∞7FF00000 00000000+∞
−∞FFF00000 00000000−∞
非数7FF80000 00000000NaN

■ 拡張倍精度 ( long double ) の記憶形式 ( Intel )

図 2-2. long double 型のビット構成
long double形式
 
表 2-2. long double 型のビット・パターン例
ビットパターン(16進)数値(10進)
00000 00000000 000000000.0
13FFF 80000000 000000001.0
-1FFFF 80000000 00000000-1.0
最大正規数7FFE FFFFFFFF FFFFFFFF1.18973149535723176505e+4932
正の最小正規数0001 80000000 000000003.36210314311209350626e-4932
最大非正規数0000 7FFFFFFF FFFFFFFF3.36210314311209350608e-4932
正の最小非正規数0000 00000000 000000013.64519953188247460253e-4951
+∞7FFF 80000000 00000000+∞
−∞FFFF 80000000 00000000−∞
非数7FFF C0000000 00000000NaN

浮動小数点数は、符号ビット ( s )・指数部 ( e )・小数部 ( f ) ( 小数部は仮数部ともいいます ) から構成されます。指数部が 0 でない場合、小数部に対して正規化が行われ、数値は ( -1 )s 1.f x 2e-b の形を取ります。ここで、b は double 型で 1023、long double 型で 16383 になります。このように正規化が行われた時の数を正規数といいますが、指数部が 0 の場合は正規化は行われず、数値は ( -1 )s 0.f x 21-b となります。この時の数は非正規数といいます。なお、上記サンプル・プログラムで非正規数のことを考慮していないのは、正の非正規数が必ず 0 と判断されるためです。
無限大や非数を表す場合、指数部の全ビットが 1 になります。このときの小数部が 0 ならば無限大、それ以外なら非数として扱われます。サンプル・プログラムでは最初に isnanisinf でチェックをしているため、後の処理においてはこの部分も考慮していません。
long double 型の構成にある j ビットは明示的先行仮数ビットと呼ばれるものです。double型では、正規化されるときに先頭に付加される 1 は暗黙的なのに対し、long double 型ではそれが j ビットで明示されます。なお、これは無限大や非数の場合も 1 になります。
最後に、上記 long double 型の構成は IntelCPU に対するもので「拡張倍精度浮動小数点数 ( Double Extended Precision Floating Point Number )」と呼ばれ、SPARCPowerPC では 128 ビットで構成される「四倍精度浮動小数点数 ( Quadruple Precision Floating Point Number )」が使われています。


3) 高速フーリエ変換 ( Fast Fourier Transform )

離散フーリエ変換によって、多項式に対する ωNk の解から、逆行列を計算することなく多項式の係数を得ることが可能となったものの、処理速度は全く改善されてはいません。しかし、ωN の性質を利用することにより、処理量を劇的に減らすことができます。まずは、ωN の正体が何であるかを説明し、その性質を利用したアルゴリズムの「高速フーリエ変換 ( Fast Fourier Transform )」を紹介します。

以下の式は、「オイラーの公式 ( Euler's formula )」といいます。

e = cosθ + i sinθ

フーリエ変換は、任意の関数を正弦波の重ね合わせで表すという発想から誕生しました。前に示した、離散フーリエ変換と逆変換の式に三角関数が含まれていないのは、オイラーの公式によって指数関数で表現されているためです。三角関数を使って表した離散フーリエ変換・逆変換の式は次のようになります ( なお、離散フーリエ変換・逆変換用のサンプル・プログラムの中で正弦・余弦を求めているのは、各 fl, Fk 値に対してこの計算を行うためだったことになります )。

Fk = (1/N)Σl{0→N-1}( fl[ cos ( -2πkl / N ) + isin ( -2πkl / N ) ] )

fl = Σk{0→N-1}( Fk[ cos ( 2πkl / N ) + i sin ( 2πkl / N ) ] )

θ = 2π をオイラーの公式に代入すると、以下の関係式が得られます。

ei2π = cos 2π + i sin 2π = 1

ここで、前節にて定義した ωN = ei2π/N を使うと

ωNN=( ei2π/N )N
=ei2π/N ・ N
=ei2π = 1

なので、ωN は 1 の N 乗根になります。また、ωN は次の関係を満たします。

ωNN = 1
ωNk ≠ 1 ( k = 1, 2, ..., N - 1 )

つまり、ωN は N 乗されて初めて 1 になります。このような条件を満たす数を「1 の原始 N 乗根 ( Primitive Root of Unity )」と言います。このことから、離散フーリエ変換を利用した方法では、1 の原始 N 乗根を使って連立方程式を解いていたことになります。

オイラーの公式を使って ωNk を三角関数で表すと、以下のようになります。

ωNk = ei2πk/N = cos ( 2πk / N ) + i sin ( 2πk / N )

つまり複素平面上において、ωNk は点 ( cos ( 2πk / N ), sin ( 2πk / N ) )を取ることになります。原点からの距離は cos2( 2πk / N ) + sin2( 2πk / N ) = 1 なので、原点を中心とする単位円を描いたとき、ωNk は必ず単位円周上に位置します。その位置は、実軸(横軸)から 2πk / N の角度を取ることから、ωN のべき乗は、原点を中心とする複素平面上の単位円を、( 1, 0 ) から始めて N 等分した位置に並ぶことになります。
下図は、複素平面の単位円周上に並んだ ωNk を示したものです。それぞれ、2, 3, 4, 8 乗根を表しています。

図 3-1. 複素平面の単位円周上に並んだ ωNk ( N = 2, 3, 4, 8 )
複素平面上における1の原始N乗根のべき乗

この周期性から、以下の関係式が成り立ちます。

ωNmN+k = ωNk ( 但し、m, k は任意の整数 )

さらに、N が偶数の場合は、以下の関係式も成り立ちます。

ωNN/2 = ( ei2π/N )N/2 = e = -1

ωNN/2+k = ωNN/2・ωNk = -ωNk ( 但し、k は任意の整数 )

ωN2 = ( ei2π/N )2 = ei2π/(N/2) = ωN/2

ωN2k = ωN/2k

上側の式は、単位円上の下側にあたる部分が上半分の符号を換えたものと等しくなることを示しています。また、下側の式は、ωN2 をべき乗すると、単位円上の N 等分点を一つおきに進み、N / 2 回で一周することを意味しています。

Xn - 1 の解には必ず 1 が含まれるので、次のように因数分解できます。

Xn - 1 = ( x - 1 )( xn-1 + xn-2 + ... + x + 1 )

N > 1 ならば ωN ≠ 1 なので、ωN が Xn - 1 = 0 の解であるためには次の関係式が必ず成り立つ必要があります。

ωNn-1 + ωNn-2 + ... + ωN + 1 = 0

前節で、離散フーリエ変換を利用して乗算処理を実際に行った例の中で、ω3に対して ω2 + ω + 1 = 0 が成り立つことを利用したところがありましたが、これは上記関係式を利用していたことになります。


ωN が持っている性質を利用すると、離散フーリエ変換処理を高速化することができます。

まず、離散フーリエ逆変換の式の一つを、Σを使わずに書き表してみます。

fl = F0 + F1ωNl + F2ωN2l + ... + FN-2ωN(N-2)l + FN-1ωN(N-1)l

ここで、N は偶数であるとき、ω の指数から l を除いた数が奇数のものと偶数のものとでちょうど半分ずつに項を分けることができます。

fl=( F0 + F2ωN2l ... + FN-2ωN(N-2)l ) + ( F1ωNl + F3ωN3l + ... + FN-1ωN(N-1)l )
=( F0 + F2ωN2l ... + FN-2ωN(N-2)l ) + ωNl( F1 + F3ωN2l + ... + FN-1ωN(N-2)l )

後半部分については、さらに ωNl を括り出していることに注意してください。こうすることによって、両括弧内にある各項の ωN の指数は全て偶数となります。ここで、ωN が持つ次の性質

ωN2k = ωN/2k

を利用すると、式を次のように変形することができます。

fl = ( F0 + F2ωN/2l ... + FN-2ωN/2(N-2)l/2 ) + ωNl( F1 + F3ωN/2l + ... + FN-1ωN/2(N-2)l/2 )

括弧の中にある二つの式は、どちらもデータ数が N / 2 個の離散フーリエ変換なので、それぞれに対して個別に処理を行っても結果が得られることを意味しています。さらに N / 2 が偶数であれば、同じ操作を行うことによって計 4 つの離散フーリエ変換を個別に行う形に変形することができます。N が 2 のべき乗であれば、最終的にはデータ数が一つの離散フーリエ変換となり、それは何も処理をしないことを意味します。

上の処理を、各 fl 毎に並べてみます。

f0=( F0 + F2 + ... + FN-2 ) + ( F1 + F3 + ... + FN-1 )
f1=( F0 + F2ωN/2 + ... + FN-2ωN/2(N-2)/2 ) + ωN( F1 + F3ωN/2 + ... + FN-1ωN/2(N-2)/2 )
f2=( F0 + F2ωN/22 + ... + FN-2ωN/2(N-2)2/2 ) + ωN2( F1 + F3ωN/22 + ... + FN-1ωN/2(N-2)2/2 )
:
fl=( F0 + F2ωN/2l + ... + FN-2ωN/2(N-2)l/2 ) + ωNl( F1 + F3ωN/2l + ... + FN-1ωN/2(N-2)l/2 )
:
fN/2=( F0 + F2ωN/2N/2 + ... + FN-2ωN/2(N-2)(N/2)/2 ) + ωNN/2( F1 + F3ωN/2N/2 + ... + FN-1ωN/2(N-2)(N/2)/2 )
fN/2+1=( F0 + F2ωN/2N/2+1 + ... + FN-2ωN/2(N-2)(N/2+1)/2 ) + ωNN/2+1( F1 + F3ωN/2N/2+1 + ... + FN-1ωN/2(N-2)(N/2+1)/2 )
fN/2+2=( F0 + F2ωN/2N/2+2 + ... + FN-2ωN/2(N-2)(N/2+2)/2 ) + ωNN/2+2( F1 + F3ωN/2N/2+2 + ... + FN-1ωN/2(N-2)(N/2+2)/2 )
:
fN/2+l=( F0 + F2ωN/22/N+l + ... + FN-2ωN/2(N-2)(N/2+l)/2 ) + ωNN/2+l( F1 + F3ωN/2N/2+l + ... + FN-1ωN/2(N-2)(N/2+l)/2 )
:

ここで、次の性質を利用します。

ωNmN+k = ωNk ( 但し、m, k は任意の整数 )

ωNN/2+k = -ωNk ( 但し、k は任意の整数、N は偶数 )

上記関係式から、次のように変形することができます。

f0=( F0 + F2 + ... + FN-2 ) + ( F1 + F3 + ... + FN-1 )
f1=( F0 + F2ωN/2 + ... + FN-2ωN/2(N-2)/2 ) + ωN( F1 + F3ωN/2 + ... + FN-1ωN/2(N-2)/2 )
f2=( F0 + F2ωN/22 + ... + FN-2ωN/2(N-2)2/2 ) + ωN2( F1 + F3ωN/22 + ... + FN-1ωN/2(N-2)2/2 )
:
fl=( F0 + F2ωN/2l + ... + FN-2ωN/2(N-2)l/2 ) + ωNl( F1 + F3ωN/2l + ... + FN-1ωN/2(N-2)l/2 )
:
fN/2=( F0 + F2 + ... + FN-2 ) - ( F1 + F3 + ... + FN-1 )
fN/2+1=( F0 + F2ωN/2 + ... + FN-2ωN/2(N-2)/2 ) - ωN( F1 + F3ωN/2 + ... + FN-1ωN/2(N-2)/2 )
fN/2+2=( F0 + F2ωN/22 + ... + FN-2ωN/2(N-2)2/2 ) - ωN2( F1 + F3ωN/2N/2+2 + ... + FN-1ωN/2(N-2)2/2 )
:
fN/2+l=( F0 + F2ωN/2l + ... + FN-2ωN/2(N-2)l/2 ) - ωNl( F1 + F3ωN/2l + ... + FN-1ωN/2(N-2)l/2 )
:

上半分と下半分を見比べると、fl と fN/2+l において、括弧内の式は同じで、それを加算するか減算するかが異なっているだけです。よって、上半分の括弧内の式を計算してしまえば、下半分の処理は不要になります。この計算方法は、図を用いて表現すると蝶の形に似た図形が現れることから「バタフライ ( Butterfly )」と呼ばれています。
以上のことから、N が 2 のべき乗であれば、最初に行った分解処理とバタフライを繰り返すことによって計算を行うことができます。この手法を「高速フーリエ変換 ( Fast Fourier Transform )」といいます。

下の図は、N = 8 の場合の分解処理とバタフライを示しています。Dn で表された部分が分解処理、Bn で表された部分がバタフライで、バタフライにおいて演算の対象を示す矢印が蝶のような形をしていることがわかると思います。

図 3-2. バタフライの概念図
高速フーリエ変換

分解処理を繰り返し実施した結果、N = 8 の場合は 0, 4, 2, 6, 1, 5, 3, 7 の順に並べ替えられています。最初の分解処理において、要素は偶数番目と奇数番目で分割されており、次の分解処理ではそれぞれのグループにおいて第 2 ビットの値が 0 のものと 1 のもので分割されます。これを繰り返すと先ほど示した順に並べ替えられることになり、それはちょうど各要素の番号に対してビットの並びを反転し、小さい順に並べる処理と同等になります。この置換方法をビット反転といいます。以下にその内容を示します。

表 3-1. ビット反転の例
番号番号(2進表現)並びを反転した結果(2進表現)並びを反転した結果
00000000
41000011
20100102
61100113
10011004
51011015
30111106
71111117

高速フーリエ変換は再帰的な処理であり、プログラムもそのように記述することができます。しかし、ビット反転とバタフライを利用すれば、プログラムを再帰処理にすることは不要になります。なお、高速フーリエ変換の手法は上記以外にも存在し、今回紹介した方法は特に「Cooley-Tukey 法」と呼ばれています。Cooley-Tukey 法は、1965 年に「ジェイムズ・クーリー ( J.W.Cooley )」と「ジョン・テューキー ( J.W.Tukey )」によって発表されました。しかし、この計算法は 1805 年頃、ドイツの数学者「ガウス ( Karl Gauss )」によってすでに発見されていたそうです。

以下に、高速フーリエ変換による乗算処理のサンプル・プログラムを示します。

/*
  InvertBit : 0 から順にビット反転した配列を作成する

  bitCnt : 作成するビット列のビット数
  index : 結果を格納する配列へのポインタ
*/
void InvertBit( size_t bitCnt, vector< size_t >* index )
{
  size_t size = 1;
  for ( size_t i = 0 ; i < bitCnt ; i++ )
    size *= 2;

  index->clear();
  for ( size_t i = 0 ; i < size ; i++ ) {
    size_t j = i;
    index->push_back( 0 );
    size_t bitMask = size >> 1;
    for ( size_t k = 0 ; k < bitCnt ; k++ ) {
      if ( ( j & 1 ) != 0 ) (*index)[i] |= bitMask;
      j >>= 1;
      bitMask >>= 1;
    }
  }
}

/*
  IFft : 高速フーリエ逆変換

  n : データ数
  blockCnt : 処理を行うブロックの数
  vec : 逆変換対象データ
  index : 間引き処理後のデータの並び
*/
void IFft( BigNum::UnsignedBuff n, BigNum::UnsignedBuff blockCnt, vector< complex< long double > >& vec, const vector< size_t >& index )
{
  BigNum::UnsignedBuff h = n / 2;
  for ( BigNum::UnsignedBuff c = 0 ; c < blockCnt ; c++ ) {
    for ( BigNum::UnsignedBuff s = c * n ; s < c * n + h ; s++ ) {
      complex< long double > r = vec[index[s]];
      long double re = cosl( ( s - c * n ) * 2 * M_PIl / n );
      long double im = sinl( ( s - c * n ) * 2 * M_PIl / n );
      vec[index[s+h]] *= complex< long double >( re, im );
      vec[index[s]] += vec[index[s+h]];
      r -= vec[index[s+h]];
      vec[index[s+h]] = r;
    }
  }
}

/*
  Fft : 高速フーリエ変換

  n : データ数
  blockCnt : 処理を行うブロックの数
  vec : 変換対象データ
  index : 間引き処理後のデータの並び
*/
void Fft( BigNum::UnsignedBuff n, BigNum::UnsignedBuff blockCnt, vector< complex< long double > >& vec, const vector< size_t >& index )
{
  BigNum::UnsignedBuff h = n / 2;
  for ( BigNum::UnsignedBuff c = 0 ; c < blockCnt ; c++ ) {
    for ( BigNum::UnsignedBuff s = c * n ; s < c * n + h ; s++ ) {
      complex< long double > r = vec[index[s]];
      long double re = cosl( ( n - ( s - c * n ) ) * 2 * M_PIl / n );
      long double im = sinl( ( n - ( s - c * n ) ) * 2 * M_PIl / n );
      vec[index[s+h]] *= complex< long double >( re, im );
      vec[index[s]] += vec[index[s+h]];
      r -= vec[index[s+h]];
      vec[index[s+h]] = r;
    }
  }
}

namespace BigNum
{
  /*
    Unsigned::fftMul : 高速フーリエ変換を使った乗算処理

    v : 乗数
  */
  Unsigned& Unsigned::fftMul( const Unsigned& v )
  {
    if ( isNaN( v ) ) return( *this );

    // 分割数 N を決定 ( 2 のべき乗にする )
    UnsignedBuff i = ( ( ( size() > v.size() ) ? size() : v.size() ) - 1 ) * 2 + 1;
    UnsignedBuff n = 2;
    size_t bitCnt = 1;
    while ( i > n ) {
      n *= 2;
      ++bitCnt;
    }

    vector< complex< long double > > uAns; // u(i)
    uAns.resize( n );
    for ( i = 0 ; i < size() ; i++ ) {
      uAns[i] = complex< long double >( number_[i], 0 );
    }
    vector< complex< long double > > vAns; // v(i)
    vAns.resize( n );
    for ( i = 0 ; i < v.size() ; i++ ) {
      vAns[i] = complex< long double >( v.number_[i], 0 );
    }

    // 間引き処理 ( 実際にはビット反転した番号を作成 )
    vector< size_t > index;
    InvertBit( bitCnt, &index );

    // 高速フーリエ逆変換で u(i), v(i) を求める
    for ( i = 2 ; i <= n ; i *= 2 ) {
      UnsignedBuff blockCnt = n / i;
      IFft( i, blockCnt, uAns, index );
      IFft( i, blockCnt, vAns, index );
    }

    // w(i) = u(i) * v(i) を求める
    vector< complex< long double > > wAns;
    for ( i = 0 ; i < n ; i++ ) {
      wAns.push_back( uAns[index[i]] );
      wAns[i] *= vAns[index[i]];
    }

    // w(i) を高速フーリエ変換する
    for ( i = 2 ; i <= n ; i *= 2 ) {
      UnsignedBuff blockCnt = n / i;
      Fft( i, blockCnt, wAns, index );
    }

    // 係数を数値に変換
    clear();
    for ( i = n ; i > 0 ; --i ) {
      if ( wAns[index[i-1]].real() < 1 ) continue;
      Unsigned m( wAns[index[i-1]].real() );
      partialAdd( m, i - 1 );
    }
    operator>>=( bitCnt );

    return( *this );
  }
};

ビット反転は InvertBit にて行っています。実際の要素を並べ替えるのではなく、ビット反転した結果を別の配列に用意して、それを使ってデータ要素に間接的にアクセスする形式を取っています。
フーリエ変換・逆変換処理 ( Fft / IFft ) では、ブロックの前半と後半に分けて、それぞれを上側から順に加算・減算しています ( 計算前に、後半部分に対して ωNk を掛けています )。ここがバタフライ処理にあたる箇所になります。

離散フーリエ変換と高速フーリエ変換による乗算処理の性能評価結果を以下に示します。各グラフの X軸は処理する数の桁数を表し、ほぼ同じ桁の数どうしを演算処理させています。また Y軸は処理時間(sec.)を表し、処理時間は演算を 25回繰り返した合計時間になります。計測は Windows7 上の VirtualBox でゲスト OS を Xubuntu 16.04 LTS として行っています。CPU は Intel Core i7-2600 3.40GHz で、ゲスト OS 上ではコア数を 1 としています。

表 3-2.・図 3-3. DFT 法の処理時間計測結果
被乗数桁数乗数桁数結果桁数演算時間(sec)
55100.000095
2472474930.161937
4884879740.393537
72972814570.864114
97096919391.51814
1210121024203.10833
1454145129044.04207
1694169333865.28991
1934193338666.69977
2175217543508.94490
24172417483410.7603
26572656531313.3994
29002898579815.5912
31413141628118.2423
33803380676021.3117
36223621724224.4596
38643863772627.5337
41044103820631.5307
43444343868645.2617
45864585917139.9662
DFT法

DFT の場合、筆算処理と比較してもかなり処理時間が長くなります ( X 軸の桁数に注意してください )。演算は二重のループとなっていることから、処理時間は桁数の二乗のオーダーで増加しています。この結果を見る限り、DFT を使って乗算を行うのは無理があるということになります。

表 3-3. FFT 法の処理時間計測結果
被乗数桁数乗数桁数結果桁数演算時間(sec)
55100.000047
2375237247460.175960
4742474294840.338378
71117110142210.391403
94799478189570.361585
1184911846236940.770190
1421714214284300.776966
1658316581331641.31860
1895618955379100.792095
2133121326426572.42663
2369323693473861.65691
2606126056521171.67226
2842928426568541.67029
3080230798616001.67791
3316433162663261.70338
3553535534710691.69340
3790737901758082.26130
4027440267805413.79886
4264442642852863.85685
4501045006900153.79873
 
図 3-4. FFT 法の処理時間計測結果 ( 筆算処理との比較 )
筆算処理による乗算FFT法
筆算処理FFT法

FFT 法は今までの高速演算法と比較しても最も最速となりました。FFT法は Karatsuba法と同じように、処理時間が突然大きくなる箇所がいくつか見られます。これは、データ量が 2 のべき乗となる値に変換されるために発生します。FFT 法は、データ量 N が 2 のべき乗である場合 Nlog2N の計算量となるため、それを反映して処理時間はほとんどデータ量に比例して増加しています。

高速な乗算処理が可能で、実装も非常にシンプルにできるということで、いいことばかりのように見える FFT ですが、サンプル・プログラムには致命的な欠陥があり、あるサイズを越えると演算時の精度により正しくない結果が得られるという問題があります。具体的には、演算回数が 65536 ( = 216 )を越えると変換結果の最大ビット数が 32 ビットとなり、これらを乗算処理したときに ( = 64 ビット ) long double 型が持つ有効桁数の 263 ビットを超えてしまうため下位の桁が削られてしまう場合があり、この時は乗算結果も正しくなくなります。正しく計算することができる乗数・被乗数の限界の桁数は約 70000 になります。


乗算処理の高速化をテーマに二章に分けて紹介してきましたが、桁数に限界がない実装方法となると、Toom-Cook 法が最有力候補となります。桁数 ( 値の大きさ ) をチェックして、ある範囲内なら FFT を利用する手段もあり、演算対象の数にある程度限りがあるのならそれが最も高速な手段になり得ると思います。


補足1) 行列式の性質

行列式の定義については「確率・統計 (5) 正規分布」の「補足 4) 行列の積の行列式」でも紹介していますが、ここで再度定義内容を確認すると

|A| = Σ( i1 i2 ... iN )( σ( i1 i2 ... iN ) a1,i1・a2,i2・...・an,iN )

と表されます。但し、( i1 i2 ... iN ) は「順列」を表し、1 から N までの数値を任意の順番で並べたものになります。Σ( i1 i2 ... iN ) は全ての順列に対する和を表し、ar,c は r 行 c 列の要素を表します。σ( i1 i2 ... iN ) は、順列 ( i1 i2 ... iN ) が「偶置換」ならば 1、「奇置換」ならば -1 になるとします(これを「順列符号」といいます)。全ての順列は 1 から N までならんだ数列に対して任意の回数分要素の入れ替え ( これを「互換」と言います ) を行うことによって得られ、互換を偶数回行って得られる順列を偶置換、奇数回の場合は奇置換といいます。

行列式の性質として次のようなものがあります。

  1. 行列 A の任意の一行 ( または一列 ) を k 倍して得られる行列 B に対し、|B| = k|A| が成り立つ
  2. 行列 A の任意の一行 ( または一列 ) を二つの行ベクトル ( または列ベクトル ) に分解し、その行(列) をそれぞれのベクトルに置き換えた二つの行列を A1, A2 としたとき、|A| = |A1| + |A2| が成り立つ
  3. 行列の任意の二行 ( または二列 ) を入れ替えると行列式は符号を換える ( 絶対値は等しい )

1 については行列式の定義から和の項に必ず k 倍された要素がただ一つだけ含まれることになるので明らかです。また、2 に関してもほぼ同様の考え方で成り立つことが示せます。3 については、隣り合った行・列を入れ替えると和の各項において互換が一回余分に発生することから奇置換は偶置換に、偶置換は奇置換に換わります。この、隣り合った行・列の入れ替えを繰り返すことで任意の場所の行・列を入れ替えることが可能ですが、 k 行(列) だけ離れた行・列を入れ替えるためには 2k - 1 回の入れ替えが必要です。従って、必ず符号を換えることになります。

これらの性質を使って、余因子展開を証明することができます。

|A| = Σi{1→n}( ( -1 )r+i・ari・cof Ari )

について見ると、A の r 列目の列ベクトル ar は次のように展開することができます。

ar=ar1|1| + ar2|0| + ... + arn|0|
|0||1||0|
:::
|0||0||1|

よって |A| は、第 r 列が i 行目のみ 1 でその他は 0 になるような行列の行列式に ari を掛けたものに分解することができます。第 r 列を第一列と入れ替え、その i 行目にある 1 の要素が第一行になるようにさらに入れ替えを行うと、係数は ( -1 )r+i だけ掛けあわせた値になります。このとき、一列目の列ベクトルは ( 1, 0, ... 0 )T になるので、一列目からは一行目以外の要素を選択すると積がゼロになり、結果として r 列目 i 行目を除いた行列の行列式、すなわち余因子 cof Ari と等しくなります。これにより、余因子展開が証明されました。

行列 B は、r' 行目以外は行列 A と等しく、r' 行目は A の r 行目と等しい要素を持つものとします。但し、r' ≠ r とします。このとき、B の余因子 cof Br'iA の余因子 cof Ar'i ( i = 1, 2, ... n ) は明らかに等しくなります。よって、行列 B の r 行 c 列目の要素を brc とすると、余因子展開の式から

Σi{1→n}( ( -1 )r'+i・ari・cof Ar'i ) = Σi{1→n}( ( -1 )r'+i・br'i・cof Br'i ) = |B|

が成り立ちます。B が要素が等しい二つの行ベクトルを持っています。これらの行を入れ替えた行列を B' としたとき、明らかに B = B' なので |B| = |B'| となります。ところが、行を入れ替えた行列に対しては符号が入れ替わるので |B| = -|B'| です。これが成り立つためには |B| = 0 でなければなりません。従って、

Σi{1→n}( ( -1 )r'+i・ari・cof Ar'i ) = 0

が成り立つことになります。r' = r のときは明らかに余因子展開と同じ式になるので、

Σi{1→n}( ( -1 )r'+i・ari・cof Ar'i ) = δr'r|A|

となります。但し、δr'r は「クロネッカのデルタ ( Kronecker Delta )」です。この結果から、余因子行列 A~ を使って

AA~ = |A|E

と表されることになり、

A-1 = A~ / |A|

が示されたことになります。


<参考文献>
  1. 「これなら分かる応用数学教室」 金谷健一著 (共立出版)
  2. Sophere 多倍長演算について
  3. PIT Lab 高速フーリエ変換(FFT)
  4. Wikipedia

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system