数値演算法

(12) 素因数分解 -2-

前章では、比較的単純な素因数分解の手法と、もっと複雑な二次ふるい法を紹介しました。今回も引き続き、素因数分解の手法について紹介したいと思います。

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

1) 複数多項式二次ふるい法 ( Multiple Polynomial Quadratic Sieve ; MPQS )

二次ふるい法 ( Quadratic Sieve ; QS )」は、素因数分解対象の数 N に対して f(r) = r2 - N の中から小さな素因数で分解できるものをふるいにかけて、f(r) の積がちょうど平方数になるような組み合わせをガウスの消去法で見つけるという処理を行っていました。f(r) の積が平方数 F2 となったとき、ri2 ( i = 1, 2, ... ) の積はやはり平方数なので、F2 ≡ (Πiri)2 ( mod N ) を満たすことになり、F - Πiri と N の最大公約数が求める素因数となる可能性があります。

r の値は、f(r) がなるべく小さくなるよう中央値を √N に近い整数としました。求める r の幅を M としたとき、r は √N - M / 2 から √N + M / 2 までの値をとることになります。N が M に対して十分に大きい数であれば √N - M / 2 > 0 なので、f(r) は √N のときゼロになり、f( √N - M / 2 ) が最小値、f( √N + M / 2 ) が最大値となります。

f( √N ± M / 2 )=( √N ± M / 2 )2 - N
=N ± M√N + ( M / 2 )2 - N
=M( √N + M / 4 ), -M( √N - M / 4 )

より、N が M に対して十分に大きい数であれば f(r) の大きさはおよそ ±M√N の幅を持ちます。

図 1-1. f(r) = r2 - N の取りうる幅
f(r)=r^2-Nの取りうる幅

r が √N に近い値であるほど f(r) は小さくなり、その分、素因数分解できる可能性が大きくなることが見込まれます。よって、f(r) = r2 - N の形よりも f(r) がなるべく小さくなるような多項式がないかを検討します。そのような多項式を仮に

F(r) = ar2 - 2br + c

とおきます。両辺に a を掛けると

aF(r)=( ar )2 - 2abr + ac
=( ar - b )2 - b2 + ac

となって、b2 - ac = N ならば

aF(r) = ( ar - b )2 - N

なので f(r) = r2 - N と同じ形となります。aF(r) が小さな素因数 p で割り切れるとき、( ar - b )2 ≡ N ( mod p ) であり、因数基地としては ( N / p ) = 1 ( 平方剰余 ) であるものだけを考えればいいので、二次ふるい法と全く同じものが利用できます。

a > 0 のとき、F(r) の最小値は r = b / a のときで -N / a です。r の幅を M としたとき最大値は r = ±M / 2 + b / a のときで aM2 / 4 - N / a なので、どちらも最小化するためには両者の絶対値が等しくなるようにすればよく、

aM2 / 4 - N / a = N / a

より

a ≅ 2( 2N )1/2 / M

となります。

N = b2 - ac ≡ b2 ( mod a )

より上記合同式の解を b とします。この合同式は、前章で紹介した「レーマーのアルゴリズム」で解くことができますが ( *1-1 )、N は a を法として平方剰余であることが前提となることに注意して下さい。最後に

c = ( b2 - N ) / a

で c を決めます。

F(r) の絶対値の最大は | F( b / a ) | で求められます。その値は

| F( b / a ) | = N / a = M√N / 2√2

で f(r) の約 1 / 3 程度にすることができます。この手法を「複数多項式二次ふるい法 ( Multiple Polynomial Quadratic Sieve ; MPQS )」といいます。

図 1-2. F(r) = ar2 - 2br + c の取りうる幅
F(r)=ar^2-2br+cの取りうる幅

実装にあたって、いくつか注意しなければならない点があります。複数多項式二次ふるい法では、aF(r) が 因数基地の素数 p で割り切れるなら

N ≡ ( ar - b )2 ( mod p )

が成り立つのでした。つまり、N ≡ t2 ( mod p ) を満たす t に対して ar - b ≡ ±t ( mod p ) であることになり、この合同式を解く必要があります。t の値は N ≡ t2 ( mod p ) の解、すなわち前回のサンプル・プログラムの中では関数 QC_Solver で求めた変数 sqrt[i] の値になります。まずは、この一次合同式を求めるサンプル・プログラムを以下に示します。

/*
  CalcGCD : ユークリッドの互除法による一次方程式の解の計算

  ax + by = GCD( a, b ) の解を x, y に求める

  a > b を前提としていることに注意

  xIsPos, yIsPos : x, y の符号が正なら true を返す

  戻り値 : a, b の最大公約数 GCD( a, b )
*/
template< class T >
T CalcGCD( T a, T b, T* x, bool* xIsPos, T* y, bool* yIsPos )
{
  *x = 1;
  *y = 0;
  *xIsPos = *yIsPos = true;
  T x1 = 0;
  T y1 = 1;
  bool xIsPos1 = true;
  bool yIsPos1 = true;

  T x2, y2;
  bool xIsPos2, yIsPos2;
  do {
    T q = a / b;
    T r = a % b;

    if ( ( *xIsPos && xIsPos1 ) || ( ( ! *xIsPos ) && ( ! xIsPos1 ) ) ) {
      if ( *x > q * x1 )
        x2 = *x - q * x1;
      else
        x2 = q * x1 - *x;
      xIsPos2 = ( *xIsPos ) ? *x >= q * x1 : *x <= q * x1;
    } else {
      x2 = *x + q * x1;
      xIsPos2 = *xIsPos;
    }

    if ( ( *yIsPos && yIsPos1 ) || ( ( ! *yIsPos ) && ( ! yIsPos1 ) ) ) {
      if ( *y > q * y1 )
        y2 = *y - q * y1;
      else
        y2 = q * y1 - *y;
      yIsPos2 = ( *yIsPos ) ? *y >= q * y1 : *y <= q * y1;
    } else {
      y2 = *y + q * y1;
      yIsPos2 = *yIsPos;
    }

    *x = x1;
    *y = y1;
    *xIsPos = xIsPos1;
    *yIsPos = yIsPos1;
    x1 = x2;
    y1 = y2;
    xIsPos1 = xIsPos2;
    yIsPos1 = yIsPos2;

    a = b;
    b = r;
  } while ( b != 0 );

  return( a );
}

/*
  mod1st : 一次合同式 ax - b = t ( mod p ) の解 x を ans に求める
*/
template< class T >
bool mod1st( T t, T a, T b, T p, T* ans )
{
  if ( a == p ) return( false );

  T x, y;
  bool xIsPos, yIsPos;

  t += b;
  if ( a < p ) {
    CalcGCD( p, a, &x, &xIsPos, &y, &yIsPos );
    if ( ! yIsPos ) y = ( ( y / p ) + 1 ) * p - y;
    *ans = ( t * y ) % p;
  } else {
    CalcGCD( a, p, &x, &xIsPos, &y, &yIsPos );
    if ( ! xIsPos ) x = ( ( x / p ) + 1 ) * p - x;
    *ans = ( t * x ) % p;
  }

  return( true );
}

一次合同式は

ax - b ≡ t ( mod p )

の形なので、p と a が互いに素なら

x ≡ ( t + b ) / a ( mod p )

で求めることができます。ここで除算が入るので、一次方程式

au + pv = 1

から a の逆数を u に求め、

x ≡ ( t + b ) x u ( mod p )

より目的の値 x を得ます ( *1-2 )。一次方程式の計算には「ユークリッドの互除法」を利用します。因数基地の素数 p の対数 log p の和を計算するときに上記合同式を解き、その結果を r1 としたとき、r の初期値を r0 として、r1 - r0 が log p を加算する添字の初期値となります。そこから p ずつ添字を進めながら加算していくのは二次ふるい法の場合と同様です。

例えば 3x - 2 = 2 ( mod 5 ) のときに aF(r) が 5 で割り切れるとしたとき、x ≡ 3 ( mod 5 ) と求められるので、r0 = 9 のとき、x ≡ 3 ≡ 13 ( mod 5 ) より 13 - 9 = 4 が求める添字の初期値となります。

表 1-1. 添字の求め方の例
添字r3r - 2 ( mod 5 )
090
1103
2111
3124
4132
5140

もう一つの注意点は、ふるいにより残った候補を試し割りするときの対象についてです。平方数となる値は aF(r) ですが、a は素数であり、因数基地に含まれているという保証はありません。そのため、試し割りするときは F(r) を対象に行います。しかし、aF(r) には少なくとも一つ a が含まれることになるため ( 因数基地に a が含まれない場合、F(r) に a が含まれると試し割りの結果、候補から外されるので、aF(r) は必ず一つの a を持つことになります )、aF(r) の積が平方数にならない場合が発生してしまいます。そこで、

p ≅ [ 2( 2N )1/2 / M ]1/2

となるような、N を法として平方剰余な素数 p を求め、a = p2 とします。試し割りの対象は F(r) としたとき、p が因数基地に含まれれば試し割りで指数が得られ、含まれなければ p を持つ候補は除外されます。いずれの場合も、aF(r) = p2F(r) は指数が 2 増える形となるため結果に影響しなくなります。

このとき問題となるのが b を求める操作です。N ≡ b2 ( mod a ) ではなく N ≡ b2 ( mod p2 ) を求めなければならないので、レーマーのアルゴリズムをそのまま利用することはできなくなります。幸い、N ≡ b2 ( mod p ) をレーマーのアルゴリズムで解いた結果を t としたとき、u ≡ ( N - t2 ) / 2t ( mod p2 ) を求めれば t + u が解となります ( 補足 1 )。


複数多項式二次ふるい法のサンプル・プログラムを以下に示します。

/*
  FindFactor : F(r) = ar^2 - 2br + c を因数基地 base で素因数分解する ( MPQS用 )

  coef : 各素因数の指数が奇数なら true とする二値を登録する配列へのポインタ
  ar_b : | ar - b | を登録する変数へのポインタ
  frBuf : a・F(r) を登録する変数へのポインタ

  coef は登録する並び順が逆であることに注意

  戻り値 : 本登録を行ったら(因数基地の素因数で完全に割り切れたら) true を返す
*/
bool FindFactor( BigNum::Unsigned r, BigNum::Unsigned a, BigNum::Unsigned b, BigNum::Unsigned c,
                 const std::vector< BigNum::Digit >& base, std::vector< std::vector< bool > >* coef,
                 BigNum::Unsigned* ar_b, BigNum::Unsigned* frBuf )
{
  // coef 用バッファ
  std::vector< bool > cBuff( base.size() + 2, false );
  // 登録する並び順は逆にする
  std::vector< bool >::reverse_iterator ci = cBuff.rbegin();

  // F(r) = ar^2 - 2br + c の計算
  // 同時に -1 の指数を登録する
  *ar_b = a * r;
  BigNum::Unsigned fr( *ar_b * r + c );
  BigNum::Unsigned br( b * r * 2 );
  if ( fr < br ) {
    fr = br - fr;
    *ci = true;
  } else {
    fr = fr - br;
  }
  ++ci;

  // ar - b と a・F(r) の計算
  if ( *ar_b > b )
    *ar_b -= b;
  else
    *ar_b = b - *ar_b;
  *frBuf = fr * a;

  unsigned int exp = 0; // 指数カウント用バッファ

  // 2 の指数の登録
  while ( ( fr % 2 ) == 0 ) {
    ++exp;
    fr /= 2;
  }
  if ( ( exp % 2 ) > 0 )
    *ci = true;
  ++ci;

  // 3 以上の素因数の指数の登録
  for ( std::vector< BigNum::Digit >::const_iterator b = base.begin() ; b != base.end() ; ++b ) {
    exp = 0;
    while ( ( fr % ( *b ) ) == 0 ) {
      ++exp;
      fr /= *b;
    }
    if ( ( exp % 2 ) > 0 )
      *ci = true;
    ++ci;
  }

  // 完全に割り切れた場合だけ本登録を行う
  if ( fr == 1 ) {
    coef->push_back( cBuff );
    return( true );
  } else {
    return( false );
  }
}

/*
  MPQS : 複数多項式二次篩 (MPQS) による n の素因数分解

  n : 素因数分解する対象の数
  u, v : 分解した数を返す変数へのポインタ
  a, b, c : 多項式の係数
  baseMax : 因数基地(Factor Base)の最大素因数
  rMax : 探索を行う最大幅
  t : ターゲット値計算時の定数
  rem : n に掛け合わされた数 ( 8 を法として 1 に合同にするため )

  戻り値 : 素因数分解できたら true を返す
*/
bool MPQS( BigNum::Unsigned n, BigNum::Unsigned* u, BigNum::Unsigned* v,
           BigNum::Unsigned a, BigNum::Unsigned b, BigNum::Unsigned c,
           BigNum::Digit baseMax, BigNum::Digit rMax, double t, BigNum::Unsigned rem )
{
  // r の最小値
  BigNum::Unsigned r0( b / a - rMax / 2 );

  // 因数基地を求める
  std::vector< BigNum::Digit > base;
  FindPrime( n, baseMax, &base );
  if ( base[0] == 2 ) base.erase( base.begin() ); // 2 は除く

  // 因数基地を法とする n の平方根を求める
  std::vector< BigNum::Digit > sqrt( base.size() );
  for ( size_t i = 0 ; i < base.size() ; ++i ) {
    sqrt[i] = QC_Solver( n, BigNum::Unsigned( base[i] ) )[0];
  }

  // 因数基地の素数の対数を割り切れる対象に加算する
  std::vector< double > logSum( rMax, double() );

  // 2 は別扱い ( 偶数なら 8 で割り切れる )
  BigNum::Unsigned ar_b( a * r0 ); // | a・r0 - b | の値
  if ( ar_b > b )
    ar_b -= b;
  else
    ar_b = b - ar_b;
  for ( std::vector< double >::size_type i = ( ( ar_b % 2 ) == 0 ) ? 1 : 0 ; i < logSum.size() ; i += 2 )
    logSum[i] += std::log( 8 );

  BigNum::Unsigned r1;
  for ( std::vector< BigNum::Digit >::size_type i = 0 ; i < sqrt.size() ; ++i ) {
    // a * r1 - b = ±sqrt ( mod base ) になる r1 を計算
    // r0 との差異から log を加える index が求まる
    if ( ! mod1st( BigNum::Unsigned( sqrt[i] ), a, b, BigNum::Unsigned( base[i] ), &r1 ) )
      continue;
    if ( r1 < r0 )
      r1 += ( ( r0 - r1 ) / base[i] + 1 ) * base[i];
    r1 = ( r1 - r0 ) % base[i];
    for ( size_t j = r1[0] ; j < logSum.size() ; j += base[i] )
      logSum[j] += std::log( base[i] );

    if ( ! mod1st( BigNum::Unsigned( base[i] - sqrt[i] ), a, b, BigNum::Unsigned( base[i] ), &r1 ) )
      continue;;
    if ( r1 < r0 )
      r1 += ( ( r0 - r1 ) / base[i] + 1 ) * base[i];
    r1 = ( r1 - r0 ) % base[i];
    for ( size_t j = r1[0] ; j < logSum.size() ; j += base[i] )
      logSum[j] += std::log( base[i] );
  }

  // ターゲット値の計算
  double target = Log( n ) / 2.0 + std::log( rMax ) - t * log( base.back() );

  // ガウスの消去法用の係数行列を求める
  std::vector< std::vector< bool > > coef;
  std::vector< BigNum::Unsigned > r;
  std::vector< BigNum::Unsigned > fr;
  for ( size_t i = 0 ; i < logSum.size() ; ++i ) {
    if ( logSum[i] <= target ) continue;
    BigNum::Unsigned ar_b, frBuf;
    if ( FindFactor( r0 + i, a, b, c, base, &coef, &ar_b, &frBuf ) ) {
      r.push_back( ar_b );
      fr.push_back( frBuf );
    }
  }

  // 解行列は単位行列で初期化
  std::vector< std::vector< bool > > ans;
  for ( size_t i = 0 ; i < coef.size() ; ++i ) {
    ans.push_back( std::vector< bool >( coef.size(), false ) );
    ans.back()[i] = true;
  }

  // ガウスの消去法を使い、平方数となる組み合わせを探す
  while ( FindFactorFromSquare( &coef, &ans, r, fr, base, n, u, v ) ) {
    if ( *u % rem == 0 )
      *u /= rem;
    else
      *v /= rem;
    if ( *u != 1 && *v != 1 ) return( true ); // 割り算の後で 1 になっていないか再確認
  }

  return( false );
}

/*
  MPQS : 複数多項式二次篩 (MPQS) による n の素因数分解

  多項式の係数を自動的に求める

  n : 素因数分解する対象の数
  u, v : 分解した数を返す変数へのポインタ
  baseMax : 因数基地(Factor Base)の最大素因数
  rMax : 探索を行う最大幅
  t : ターゲット値計算時の定数

  戻り値 : 素因数分解できたら true を返す
*/
bool MPQS( BigNum::Unsigned n, BigNum::Unsigned* u, BigNum::Unsigned* v,
           BigNum::Digit baseMax, BigNum::Digit rMax, double t )
{
  // r^2 - n が 8 を因数として持つようにするため
  // n は法を 8 として 1 に合同になるようにする
  BigNum::Unsigned rem = n % 8;
  n *= rem;

  // 多項式の係数を計算する
  // a = p^2 となるようにする ( 但し p は ( n / p ) = 1 を満たす素数 )
  BigNum::Unsigned p( MathLib::SquareRoot( MathLib::SquareRoot( n * 2 ) * 2 / rMax ) );
  while ( ( ! PrimalityTest::MillerRabinTest( p ) ) || ( MathLib::Jacobi( n, p ) < 1 ) )
    ++p;
  BigNum::Unsigned a( p * p );

  // b^2 = n ( mod p^2 ) を満たすように b を定める
  // t^2 = n ( mod p ) としたとき、t + kp が求める解となる
  // 但し、kp = ( n - t^2 ) / 2t ( mod p^2 )
  BigNum::Unsigned b( QC_Solver( n, p ) ); // b にいったん t を求めておく
  BigNum::Unsigned k( n - b * b ); // k にいったん n - t^2 を代入

  // 2t で割り算するために 2t の逆数 div を求める
  // (2t)x + (p^2)y = 1 を満たす x が求める div となる
  // p^2 > 2t を前提としていることに注意
  BigNum::Unsigned x, y;
  bool xIsPos, yIsPos;
  CalcGCD( a, b * 2, &x, &xIsPos, &y, &yIsPos );
  BigNum::Unsigned div( y );
  bool divIsPos = yIsPos;
  if ( ! divIsPos ) div = ( ( div / a ) + 1 ) * a - div;
  b += k * div;
  b %= a;

  // b に ( rMax / 2 ) x a だけゲタ履きする ( r の初期値をゼロにするため )
  b += ( rMax / 2 ) * a;
  // 最後に c = ( b^2 - n ) / a を求める
  BigNum::Unsigned c( ( b * b - n ) / a );

  return( MPQS( n, u, v, a, b, c, baseMax, rMax, t, rem ) );
}

前述した注意点以外は、前章で紹介した二次ふるい法と基本的に変化はありません。FindFactor は因数基地で試し割りをする関数で、これだけは二次ふるい法のサブ・ルーチンを流用できないため新たに作成していますが、ガウスの消去法を行う部分については前回のプログラムがそのまま使えます。

メイン・プログラムは MPQS ですが、係数を任意とするものと、自動的に計算するものの二種類があり、後者が前者を呼び出す形となっています。係数の求め方については先に説明したとおりですが、r の初期値を正数とするため b の値は ( rMax / 2 ) * a だけゲタ履きさせていることに注意して下さい。


*1-1)数値演算法 (11) 素因数分解 -1-」の「4) リュカ数列 ( Lucas Sequence )」を参照

*1-2)数値演算法 (11) 素因数分解 -1-」の「補足 2) 合同式の除算」を参照


2) ペル方程式 ( Pell's Equation )

m を平方数ではない任意の正数としたとき、次の形をした方程式を「ペル方程式 ( Pell's Equation )」といいます。

x2 - my2 = 1

上式は

( x + y√m )( x - y√m ) = 1

と変形することができます。両辺を n 乗しても右辺は 1 のままであり、

( x ± y√m )n=Σk{0→n}( nCkxn-k(±y√m)k )
=Σk{0→n/2}( nC2kxn-2ky2kmk ) ± √mΣk{0→(n-1)/2}( nC2k+1xn-2k-1y2k+1mk )
S1 ± S2√m

より

( x + y√m )n( x - y√m )n=( S1 + S2√m )( S1 - S2√m )
=S12 - S22m

なので、整数解 ( x, y ) が得られたとすれば ( S1, S2 ) は新たな整数解となります。ペル方程式の解は、このようにべき乗することでのみ生成され、どんな整数 m に対しても必ず解を持つことが証明されています。

生成される解は常に大きくなるため、ペル方程式は必ず最小解を持ちます。この最小解を求める方法に「バスカラ-ブラウンカーのアルゴリズム ( Bhaskara-Brouncker Algorithm )」があります。これは次の一連の漸化式から成り立っています。s は ⌊ √m ⌋ ( √m を超えない最大の整数 ) とします。

B0 = 0 ; C0 = 1 ; P0 = 1 ; Q0 = 0

A1 = s ; B1 = s ; C1 = m - s2 ; P1 = s ; Q1 = 1

An+1 = ⌊ ( s + Bn ) / Cn

Bn+1 = An+1Cn - Bn

Cn+1 = Cn-1 + An+1( Bn - Bn+1 )

Pn+1 = Pn-1 + An+1Pn

Qn+1 = Qn-1 + An+1Qn

例えば m = 13 のとき、s = 3 なので、各パラメータは以下のように推移します。

nAnBnCnPnQn
0-0110
133431
211341
312372
4114113
5131185
663411933
711313738
812325671
9114393109
10131649180

先ほど示した漸化式において、次の等式が成り立つことが証明できます ( 補足 2 )。

Pn2 - mQn2 = ( -1 )nCn

m = 13 の場合、n = 5 のときに Cn = 1 となっており、このときの Pn = 18, Qn = 5 を使って実際に計算してみると

Pn2 - mQn2 = 182 - 13 x 52 = 324 - 325 = -1

となり成り立つことがわかります。また、次の n = 10, Pn = 649, Qn = 180 では

Pn2 - mQn2 = 6492 - 13 x 1802 = 421201 - 421200 = 1

となって、これはペル方程式 x2 - 13y2 = 1 の最小解となります。

以前、無理数がユークリッドの互除法を使って連分数の形に変換できることを示しました ( *2-1 )。√13 について実際に試してみると、

√13=3x1+( -3 + √13 )( 3.605551 = 3 x 1 + 0.605551 )
1=1x( -3 + √13 )+( 4 - √13 )( 1 = 1 x 0.605551 + 0.394449 )
-3 + √13=1x( 4 - √13 )+( -7 + 2√13 )( 0.605551 = 1 x 0.394449 + 0.211103 )
4 - √13=1x( -7 + 2√13 )+( 11 - 3√13 )( 0.394449 = 1 x 0.211103 + 0.183346 )
-7 + 2√13=1x( 11 - 3√13 )+( -18 + 5√13 )( 0.211103 = 1 x 0.183346 + 0.027756 )
11 - 3√13=6x( -18 + 5√13 )+( 119 - 33√13 )( 0.183346 = 6 x 0.027756 + 0.016808 )
-18 + 5√13=1x( 119 - 33√13 )+( -137 + 38√13 )( 0.027756 = 1 x 0.016808 + 0.010948 )
119 - 33√13=1x( -137 + 38√13 )+( 256 - 71√13 )( 0.016808 = 1 x 0.010948 + 0.005859 )
-137 + 38√13=1x( 256 - 71√13 )+( -393 + 109√13 )( 0.010948 = 1 x 0.005859 + 0.005089 )
256 - 71√13=1x( -393 + 109√13 )+( 649 - 180√13 )( 0.005859 = 1 x 0.005089 + 0.000770 )

となるので、

√13 = 3 + 1 / [ 1 + 1 / [ 1 + 1 / [ 1 + 1 / [ 1 + 1 / [ 6 + ...

と変換することができます。ここで、ユークリッドの互除法の右辺において、符号を別として An, Pn, Qn の列が現れていることに注意して下さい。連分数を見ると、

√m = A1 + 1 / [ A2 + 1 / [ A3 + 1 / [ A4 + 1 / [ A5 + 1 / [ A6 + ...

が成り立つことが予想できます ( 補足 2 )。

また、Pn / Qn を計算してみると、

3 / 1 = 3

4 / 1 = 4

7 / 2 = 3.5

11 / 3 = 3.666666666...

18 / 5 = 3.6

119 / 33 = 3.606060606...

137 / 38 = 3.605263157...

256 / 71 = 3.605633802...

393 / 109 = 3.605504587...

649 / 180 = 3.605555555...

となり、

√13 = 3.605551275...

に近づく様子が確認できます。このことから Pn / Qn は平方根の近似値として利用できることがわかります。実際、

Pn / Qn = A1 + 1 / [ A2 + 1 / [ A3 + ... + 1 / [ An-1 + 1 / An ] ... ] ]

が成り立ちます ( 補足 2 )。


バスカラ-ブラウンカーのアルゴリズムのサンプル・プログラムを以下に示します。

/*
  BhaskaraBrouckner : バスカラ-ブラウンカーのアルゴリズム

  次の漸化式を計算する。s は n の平方根を超えない最大の整数。

  An+1 = INT( ( s + Bn ) / Cn )
  Bn+1 = An+1 x Cn - Bn
  Cn+1 = Cn-1 + An+1 x ( Bn - Bn+1 )
  Pn+1 = Pn-1 + An+1 x Pn
  Qn+1 = Qn-1 + An+1 x Qn

  初期値は次の通り。

  B0 = 0 C0 = 1 P0 = 1 Q0 = 0
  A1 = s B1 = s C1 = n - s x s P1 = s Q1 = 1

  An は √m の連分数表現を、Pn / Qn は √m の近似値をそれぞれ与える。
  Cn が 1 のときの Pn, Qn は、ペル方程式 x^2 - my^2 = ±1 の解となる。

  s m の平方根を超えない最大の整数
  a,b,c,p,q 漸化式の各パラメータ An, Bn, Cn, Pn, Qn
  c_1, p_1, q_1 漸化式の一つ前のパラメータ Cn-1, Pn-1, Qn-1
*/
template< class T >
void BhaskaraBrouckner( const T& s, T* a, T* b, T* c_1, T* c, T* p_1, T* p, T* q_1, T* q )
{
  *a = ( s + *b ) / *c;
  T b1 = *a * *c - *b;

  // 符号なし整数を利用することを想定し、正数のみで計算
  // 各パラメータが必ず正であることは保証されている
  T buff( *c );
  if ( *b < b1 )
    *c = *c_1 - *a * ( b1 - *b );
  else
    *c = *c_1 + *a * ( *b - b1 );
  *b = b1;
  *c_1 = buff;

  buff = *p;
  *p = *p_1 + *a * *p;
  *p_1 = buff;

  buff = *q;
  *q = *q_1 + *a * *q;
  *q_1 = buff;
}

サンプル・プログラムは、前回の結果 c_1, p_1, q_1 を使って更新するといった処理を行うようになっています。従って、n 番目の値が必要な場合は関数を n 回呼び出して得ることになります。処理自体は漸化式の計算そのままですが、各パラメータは全て常に正であることを前提にしています。これは実際に成り立ち、さらに

0 < An < 2√m

0 < Bn < √m

0 < Cn < 2√m

になります ( 補足 3 )。この結果から、An, Bn, Cn は有限の値であり、その組 ( An, Bn, Cn ) も有限なのでやがて必ず一致する組が現れます。各パラメータは漸化式で求められることから、一致した組が一度現れれば ( An, Bn, Cn ) の組は周期性を持つようになります ( 補足 4 )。


*2-1)数値演算法 (5) 有理数と無理数の演算」の「3) 無理数と連分数」を参照


3) 連分数法 ( Continued Fraction Factorization Method ; CFRAC )

バスカラ-ブラウンカーのアルゴリズムから以下の結果が得られました。

Pn2 - mQn2 = ( -1 )nCn

これは次の合同式が成り立つことを意味しています。

Pn2 ≡ ( -1 )nCn ( mod m )

これにより二次ふるい法と同様の手法を使うことができます。すなわち、( -1 )nCn を因数基地で素因数分解して成功したらそれを保存しておきます。素因数分解可能な ( -1 )nCn の個数が因数基地の大きさを超えたら、ガウスの消去法を使い、( -1 )nCn の積が平方数となる組み合わせを探します。この手法を「連分数法 ( Continued Fraction Factorization Method ; CFRAC )」といいます。

Pn2 - mQn2 = ( -1 )nCn より、Cn が素数 p で割り切れるならば

Pn2 ≡ mQn2 ( mod p )

となります。すなわち ( mQn2 / p ) = 1 であり、

( mQn2 / p ) = ( m / p )( Qn / p )( Qn / p ) = ( m / p )

より m は法 p の平方剰余になります。従って、因数基地としても二次ふるい法と全く同じものが利用できます。しかし、二次ふるい法が因数基地で試し割りする前に候補を絞り込むことができたのに対し、連分数法は試し割りを毎回行う必要があります。ではどこが二次ふるい法と比べて有利かというと、素因数分解の対象となる値の大きさです。二次ふるい法では、素因数分解対象の数を N、ふるいを掛ける幅を M としたとき、その大きさは ±M√N 程度となるのでした。複数多項式二次ふるい法を使ったとしても、その 1 / 3 程度までしか小さくなりません。連分数法の場合、0 < Cn < 2√N という制限のため二次ふるい法と比較してかなり小さくすることができます。


連分数法のサンプル・プログラムを以下に示します。

/*
  BhaskaraBrouckner : バスカラ-ブラウンカーのアルゴリズム(CFRACに使用する)

  次の漸化式を計算する。s は n の平方根を超えない最大の整数。

  Ai+1 = INT( ( s + Bi ) / Ci )
  Bi+1 = Ai+1 x Ci - Bi
  Ci+1 = Ci-1 + Ai+1 x ( Bi - Bi+1 )
  Pi+1 = Pi-1 + Ai+1 x Pi

  ( Qi は不要なため計算は省略 )

  初期値は次の通り。

  B0 = 0 C0 = 1 P0 = 1
  A1 = s B1 = s C1 = n - s x s P1 = s
*/
template< class T >
void BhaskaraBrouckner( const T& n, const T& s, T* a, T* b, T* c_1, T* c, T* p_1, T* p )
{
  *a = ( s + *b ) / *c;
  T b1 = *a * *c - *b;

  // 符号なし整数を利用することを想定し、正数のみで計算
  // 各パラメータが必ず正であることは保証されている
  T buff( *c );
  if ( *b < b1 )
    *c = *c_1 - *a * ( b1 - *b );
  else
    *c = *c_1 + *a * ( *b - b1 );
  *b = b1;
  *c_1 = buff;

  buff = *p;
  *p = ( *p_1 + *a * *p ) % n;
  *p_1 = buff;
}

/*
  FindFactor : (-1)^i・Ci を因数基地 base で素因数分解する ( CFRAC 用 )

  c : Ci
  cIsNeg : Ci が負数なら true
  coef : 各素因数の指数が奇数なら true とする二値を登録する配列へのポインタ

  exp と coef では登録する並び順が逆であることに注意

  戻り値 : 本登録を行ったら(因数基地の素因数で完全に割り切れたら) true を返す
*/
bool FindFactor( BigNum::Unsigned c, bool cIsNeg, const std::vector< BigNum::Digit >& base, std::vector< std::vector< bool > >* coef )
{
  std::vector< bool > cf( base.size() + 1, false ); // coef 用バッファ
  // 登録する並び順は逆にする
  std::vector< bool >::reverse_iterator cfi = cf.rbegin();

  // -1 の指数を登録する
  if ( cIsNeg ) {
    *cfi = true;
  }
  ++cfi;

  // 素因数の指数の登録
  for ( std::vector< BigNum::Digit >::const_iterator b = base.begin() ; b != base.end() ; ++b ) {
    unsigned int exp = 0;
    while ( ( c % ( *b ) ) == 0 ) {
      ++exp;
      c /= *b;
    }
    if ( ( exp % 2 ) > 0 )
      *cfi = true;
    ++cfi;
  }

  // 完全に割り切れた場合だけ本登録を行う
  if ( c == 1 ) {
    coef->push_back( cf );
    return( true );
  } else {
    return( false );
  }
}

/*
  CFRAC : 連分数アルゴリズム (CFRAC) による m の素因数分解

  m : 素因数分解する対象の数
  u, v : 分解した数を返す変数へのポインタ

  戻り値 : 素因数分解できたら true を返す
*/
bool CFRAC( BigNum::Unsigned m, BigNum::Unsigned* u, BigNum::Unsigned* v, BigNum::Digit baseMax )
{
  // バスカラ-ブラウンカーのアルゴリズム用パラメータ
  BigNum::Unsigned s( MathLib::SquareRoot( m ) );
  BigNum::Unsigned a( s );
  BigNum::Unsigned b( s );
  BigNum::Unsigned c_1( 1 ); BigNum::Unsigned c( m - s * s );
  BigNum::Unsigned p_1( 1 ); BigNum::Unsigned p( s );

  // m が平方数ならすぐに結果を返す
  if ( m == s * s ) {
    *u = *v = s;
    return( true );
  }
  // 因数基地を求める
  std::vector< BigNum::Digit > base;
  FindPrime( m, baseMax, &base );

  std::set< BigNum::Unsigned > cPool;    // Ci 用バッファ(既出かどうかの判定)
  std::vector< BigNum::Unsigned > cCand; // Ci 用候補
  std::vector< BigNum::Unsigned > pCand; // Pi 用候補
  std::vector< std::vector< bool > > coef;         // 係数行列(指数が奇数ならtrue)
  bool cIsNeg = true; // Ci が負数なら true にする
  for ( BigNum::Unsigned i( 0 ) ; i < 2 * s ; ++i ) {
    std::set< BigNum::Unsigned, bool >::const_iterator cIt = cPool.find( c );
    if ( cIt == cPool.end() ) {
      cPool.insert( c );
      if ( FindFactor( c, cIsNeg, base, &coef ) ) {
        cCand.push_back( c );
        pCand.push_back( p );
      }
    }
    BhaskaraBrouckner( m, s, &a, &b, &c_1, &c, &p_1, &p );
    cIsNeg = ! cIsNeg;

    // 係数行列の取得は、因数基地の大きさの 2 倍まで取れたら中断する
    if ( coef.size() >= base.size() * 2 ) break;
  }

  // 解行列は単位行列で初期化
  std::vector< std::vector< bool > > ans;
  for ( size_t i = 0 ; i < coef.size() ; ++i ) {
    ans.push_back( std::vector< bool >( coef.size(), false ) );
    ans.back()[i] = true;
  }

  // ガウスの消去法を使い、平方数となる組み合わせを探す
  while ( FindFactorFromSquare( &coef, &ans, pCand, cCand, base, m, u, v ) ) {
    if ( *u != 1 && *v != 1 ) return( true ); // 割り算の後で 1 になっていないか再確認
  }

  return( false );
}

関数 BhaskaraBrouckner は前節で紹介したものと処理の内容は同じですが、連分数法では Qn は使わないため計算を省略している部分が異なります。FindFactor が因数基地で試し割りをする部分であり、CFRAC がメイン・ルーチンとなります。CFRAC 内の FindPrime, FindFactorFromSquare は二次ふるい法と全く同じ関数が利用できます。

Cn の試し割りをする回数は 2√m までに制限しています。これは、0 < Cn < 2√m という制限があるため、おおよそ 2√m 回程度ループすればほぼ全ての Cn が得られるだろうという見積もりからの値です。しかし、一回の周期の中で Cn が毎回異なる値であるという保証はありません ( 前節での m = 13 の例を参照 )。また、係数行列が因数基地の大きさの 2 倍程度になれば素因数が得られるだろうと試算して、そうなったらループを抜けるようにもしています。これにより処理は高速になりますが、素因数が得られない可能性は若干減少することになります。


4) 楕円曲線 ( Elliptic Curve )

以下の式で定義される方程式を「楕円曲線 ( Elliptic Curve )」といいます。

y2 = x3 + ax + b ( 但し 4a3 + 27b2 ≠ 0 )

x が決まれば y2 = (定数) となることから明らかなように、曲線は x 軸について対称になります。y = x3 + ax + b の形ならば定義域は実数全体になりますが、左辺は ( 実数に限定すれば ) ゼロ以上である必要があるので定義域は限定されます。例えば y2 = x3 - x = x( x + 1 )( x - 1 ) の場合、x = 0, ±1 のとき y = 0 で、x ≥ 1, -1 ≤ x ≤ 0 が定義域となります。

4a3 + 27b2 ≠ 0 は楕円曲線が重解を持たないための条件です。例えば、

a = -3s2, b = 2s3

のとき 4a3 + 27b2 = 0 で、

y2=x3 - 3s2x + 2s3
=( x - s )3 + 3sx2 - 6s2x + 3s3
=( x - s )3 + 3s( x - s )2
=( x - s )2( x + 2s )

となって、y = 0 のとき重解が発生します。重解が発生したとき、楕円曲線は尖点 ( カスプ ; Cusp ) や孤立点が発生したり、自己交差するようになります ( 補足 5 )。これ以降は、4a3 + 27b2 = 0 の場合は除外することとします。

図 4-1. 様々な楕円曲線
y2 = x3 - xy2 = x3 + x
y^2=x^3-xy^2=x^3+x
y2 = x3y2 = x3 - 3x + 2
y^2=x^3y^2=x^3-3x+2

楕円曲線上の二点 ( x1, y1 ), ( x2, y2 ) を決め、その傾きを λ で表すと、

λ = ( y1 - y2 ) / ( x1 - x2 )

となります。但し、分母はゼロになれないので x1 ≠ x2 という条件が付きます。x1 = x2 のときは、さらに二つの場合に分けることができて、y1 = -y2 のとき、二点を結んだ直線は x 軸に垂直になります。y1 = y2 のときは接線を意味するので ( 但し y1 = y2 = 0 の場合は x 軸に垂直になるので除外します )、楕円曲線の方程式の両辺を x で微分すれば

2y(dy/dx) = 3x2 + a

より

f'(x) ≡ (dy/dx) = ( 3x2 + a ) / 2y

となって、

λ = f'(x1) = ( 3x12 + a ) / 2y1

という結果が得られます。

三次方程式には ( 重解を無視すれば ) 解が 3 つあるので、第 3 の点 ( x3, y3 ) も存在しているはずです。曲線と接するときは点の個数が減りますが、このときは接点上に座標の等しい 2 点があるとみなします。このとき、( x3, y3 ) と ( x1, y1 ), ( x2, y2 ) それぞれに対して

λ = ( y1 - y3 ) / ( x1 - x3 ) = ( x12 + x1x2 + x32 + a ) / ( y1 + y3 )
λ = ( y2 - y3 ) / ( x2 - x3 ) = ( x22 + x2x3 + x32 + a ) / ( y2 + y3 )

が成り立つので

( y1 + y3 )λ = x12 + x1x3 + x32 + a
( y2 + y3 )λ = x22 + x2x3 + x32 + a

となり、辺々引いて

( y1 - y2=x12 + x1x3 - x22 - x2x3
=( x1 + x2 )( x1 - x2 ) + x3( x1 - x2 )
=( x1 + x2 + x3 )( x1 - x2 )

x1 ≠ x2 であれば、両辺を x1 - x2 で割ると

λ2 = x1 + x2 + x3

より

x3 = λ2 - x1 - x2

となり、y3 は λ = ( y1 - y3 ) / ( x1 - x3 ) より

y3 = λ( x3 - x1 ) + y1

となります。( x1, y1 ) = ( x2, y2 ) のときは上式は使えないので、接線の方程式を

y = λ( x - x1 ) + y1

但し λ = ( 3x12 + a ) / 2y1

として y2 = x3 + ax + b に代入し、

[ λ( x - x1 ) + y1 ]2 = x3 + ax + b

を直接解きます。

(左辺)=λ2( x - x1 )2 + 2λy1( x - x1 ) + y12
=λ2( x - x1 )2 + ( 3x12 + a )( x - x1 ) + ( x13 + ax1 + b )

なので、

λ2( x - x1 )2 + ( 3x12 + a )( x - x1 ) + x13 + ax1 + b - ( x3 + ax + b )
=λ2( x - x1 )2 + ( 3x12 + a )( x - x1 ) - ( x3 - x13 ) - a( x - x1 )
=λ2( x - x1 )2 + ( 3x12 + a )( x - x1 ) - ( x - x1 )( x2 + x1x + x12 ) - a( x - x1 )
=( x - x1 )[ λ2( x - x1 ) + ( 3x12 + a ) - ( x2 + x1x + x12 ) - a ]
=( x - x1 )[ λ2( x - x1 ) - ( x2 + x1x - 2x12 ) ]
=( x - x1 )[ λ2( x - x1 ) - ( x + 2x1 )( x - x1 ) ]
=( x - x1 )2( λ2 - x - 2x1 ) = 0

となって、x = λ - 2x1 より x1 ≠ x2 と同じ結果が得られます。

二点が有理点 ( 有理数からなる点 ) であれば、そこから求められる傾きは有理数であり、もう一つの点も有理点となります。そこで、楕円曲線上の二点 ( x1, y1 ), ( x2, y2 ) を通る直線と楕円曲線とのもう一つの交点 ( x3, y3 ) に対し、以下のような二項演算 ⊕ を定義します。

( x1, y1 ) ⊕ ( x2, y2 ) = ( x3, -y3 )

y3 に "-" が付いていることに注意して下さい。すなわち、演算の結果は求めた交点に対し x 軸について対称な点になりますが、楕円曲線が x 軸に対して対称なので、やはり楕円曲線上の点となります。しかし、まだ ( x1, y1 ) = ( x2, -y2 ) と y1 = y2 = 0 の場合が残っています。このときも含めて上記の演算が成り立つようにするため、ここで無限遠点 ∞ を定義します。無限遠点はその名の通り無限遠にある点とみなされ、これを使って以下の式が成り立つとします。

( x1, y1 ) ⊕ ( x2, -y2 ) = ∞

二点 ( x1, y1 ), ( x2, -y2 ) は無限遠上で楕円曲線と交わり、その対称点もまた無限遠上にあります ( 実際には楕円曲線は無限遠上でつながっていると見なされるので、無限遠点はただ一つと定義されるようです )。また、ある点 ( x1, y1 ) と無限点を結んだ線、すなわち ( x1, y1 ) を通る x 軸に垂直な直線は ( x1, -y1 ) と交わるので

( x1, y1 ) ⊕ ∞ = ( x1, y1 )

となります。よって、∞ は通常の加法におけるゼロ ( これをゼロ元といいます ) と同等のものと考えることができます。また、( x1, y1 ) と ( x2, -y2 ) の演算結果がゼロ元になるということは、ちょうど整数の正負関係と同等であることを意味します ( これを逆元といいます )。以後、この二項演算で定義された楕円曲線 y2 = x3 + ax + b 上の点の集合を E( a, b ) で表します。


上記の二項演算の定義は、さらにある数 n を法としても成り立ちます。まず、二点 ( x1, y1 ), ( x2, y2 ) について、

x1 ≡ x2, y1 ≡ -y2 ( mod n ) のとき

( x1, y1 ) ⊕ ( x2, y2 ) = ∞

とします。また、x1 ≠ x2 ( mod n ) かつ x1 - x2 と n が互いに素なとき

λ ≡ ( y1 - y2 ) / ( x1 - x2 ) ( mod n )

x1 ≡ x2 ( mod n ) かつ y1 + y2 と n が互いに素なとき

λ ≡ ( 3x12 + a ) / 2y1 ( mod n )

とします。このとき、新たな点 ( x3, y3 ) を

x3 ≡ λ2 - x1 - x2 ( mod n )
y3 ≡ λ( x3 - x1 ) + y1 ( mod n )

として二項演算

( x1, y1 ) ⊕ ( x2, y2 ) ≡ ( x3, -y3 ) ( mod n )

を定義し、その集合を E( a, b ) / n で表します。n が奇素数 p であれば、この演算は必ず定義することができます。但し、ここで a3 + 27b2 ≡ 0 ( mod p ) の場合は除くことに注意して下さい。

奇素数 p を法としたとき、集合となる点 ( x, y ) は有限となります。ある点 ( x0, y0 ) があった場合、点の個数は有限なので

( x0, y0 ) ⊕ ( x0, y0 ) = ( x1, y1 )
( x1, y1 ) ⊕ ( x0, y0 ) = ( x2, y2 )
 :
( xn, yn ) ⊕ ( x0, y0 ) = ∞
∞ ⊕ ( x0, y0 ) = ( x0, y0 )

となり巡回して元に戻るはずです。( x0, y0 ) ⊕ ( x0, y0 ) ⊕ ... ( x0, y0 ) と n 回演算を行ったとき、この演算を ( x0, y0 )n で表すことにします。また、( x0, y0 )n = ∞ となるときの最小の正整数を | ( x0, y0 ) | とします。

y2 = x3 + x + 2 ( mod 5 ) を例に取ると、点の一つは ( 1, 2 ) になります。

( 1, 2 ) ⊕ ( 1, 2 )

λ ≡ ( 3 x 12 + 1 ) / ( 2 x 2 ) ≡ 1 ( mod 5 )
x ≡ 12 - 1 - 1 ≡ -1 ≡ 4 ( mod 5 )
y ≡ 1 x ( 4 - 1 ) + 2 ≡ 5 ≡ 0 ( mod 5 )
よって ( x, y ) ≡ ( 4, 0 ) ( mod 5 )

( 4, 0 ) ⊕ ( 1, 2 )

λ ≡ ( 0 - 2 ) / ( 4 - 1 ) ≡ -2 / 3 ≡ 1 ( mod 5 )
x ≡ 12 - 4 - 1 ≡ -4 ≡ 1 ( mod 5 )
y ≡ 1 x ( 1 - 4 ) + 0 ≡ -3 ≡ 2 ( mod 5 )
よって ( x, y ) ≡ ( 1, -2 ) ≡ ( 1, 3 ) ( mod 5 )

( 1, 3 ) ⊕ ( 1, 2 ) = ∞

∞ ⊕ ( 1, 2 ) = ( 1, 2 )

となって、確かに巡回して元に戻る様子が確認できます。また、このとき | ( 1, 2 ) | = 4 となります。

上記の例において、点は ∞, ( 1, 2 ), ( 1, 3 ), ( 4, 0 ) の 4 つありました。以下、この点の個数を | E( a, b ) / p | で表します。このとき、果たして | E( a, b ) / p | と | ( x0, y0 ) | の間の関係はどのようになるのでしょうか。例の場合、他に点は見つからず | E( 1, 2 ) / 5 | = | ( 1, 2 ) | となります。もしそうならないような場合があれば、巡回するパターンが複数あることになりますが、結論から言うと、K を 1 以上の整数としたとき | E( a, b ) / p | = K・| ( x0, y0 ) | が必ず成り立ちます ( 補足 6 )。


5) 楕円曲線法 ( Elliptic Curve Method ; ECM )

ポラードの p - 1 法 ( Pollard's p - 1 Algorithm )」は、ある数 k と奇素数 p に対して k! が p - 1 で割り切れると仮定すると、p を法とする数 c の k! 乗がフェルマーの小定理から 1 に合同になる、すなわち

ck! ≡ 1 ( mod p )

が成り立つことを利用した手法でした。ck! - 1 ( mod N ) は高速に計算することが可能で、上記の条件を満たせばこの値が p で割り切れることになります。同じようなことを楕円曲線においても考えてみます。ある数 k と m ≡ | ( x0, y0 ) | に対して k! が m で割り切れると仮定すると、

( x0, y0 )k! = ∞

が成り立ちます。∞ の場合に p で割り切れるものを利用して、その値と素因数分解対象との最大公約数から素因数を見つけ出すわけです。しかし、∞ は無限遠点という、扱うには少々面倒な値です。そこで、有理点であった ( x, y ) に対し、整数 X, Y, Z を使って

( x, y ) = ( X / Z, Y / Z )

とすることで座標系を ( x, y ) → ( X, Y, Z ) に変換します。このとき、楕円曲線の方程式は

( Y / Z )2 = ( X / Z )3 + a( X / Z ) + b より

Y2Z = X3 + aXZ2 + bZ3

と表されます。( x, y ) の変換式から、任意のゼロでない整数 c に対して ( X, Y, Z ) と ( cX, cY, cZ ) が同じ有理数解 ( x, y ) となることに注意して下さい。すなわち、

( X, Y, Z ) = ( cX, cY, cZ ) ( 但し c ≠ 0 )

が成り立ちます。座標系を ( X, Y, Z ) に変換するもう一つの理由は ∞ が扱いやすくなるという点です。無限遠点を Z → 0 の極限における値と考えれば、∞ = ( X, Y, 0 ) ( mod p ) であり Z は p で割り切れます。従って、このときの Z と素因数分解対象の数 N との最大公約数を計算することによって素因数が見つかる可能性があります。ポラードの p - 1 法と同様に、p は未知数なので、楕円曲線の方程式のパラメータや点の初期値を変えながらいろいろと試す形となります。

ck! が高速に計算できるのと同様に ( x, y )k! も高速に計算することができます。y の値は最終的には不要となるため、以下では x だけ算出します。( u, v ) = ( x, y ) ⊕ ( x, y ) のとき、

λ = ( 3x2 + a ) / 2y

より

u=λ2 -2x
=[ ( 3x2 + a ) / 2y ]2 -2x
=( 3x2 + a )2 / 4( x3 + ax + b ) -2x
=[ ( 9x4 + 6ax2 + a2 ) - 8x( x3 + ax + b ) ] / 4( x3 + ax + b )
=( x4 - 2ax2 - 8bx + a2 ) / 4( x3 + ax + b )
=[ ( x2 - a )2 - 8bx ] / 4( x3 + ax + b )

であり、この公式を使って ( x, y )2i を計算していくことができます。( e, f ) = ( x, y )i, ( g, h ) = ( x, y )i+1, ( u, v ) = ( x, y )2i+1 としたとき、( u, v ) = ( e, f ) ⊕ ( g, h ) なので

λ = ( f - h ) / ( e - g )
u = λ2 - e - g = ( f - h )2 / ( e - g )2 - ( e + g )

より

u( e - g )2=( f - h )2 - ( e + g )( e - g )2
=( f2 - 2fh + h2 ) - ( e3 - e2g - eg2 + g3 )
=( e3 + ae + b ) + ( g3 + ag + b ) - 2fh - ( e3 - e2g - eg2 + g3 )
=-2fh + 2b + a( e + g ) + eg( e + g )
=-2fh + 2b + ( a + eg )( e + g )

また、( g, h ) = ( x, y ) ⊕ ( e, f ) より ( x, y ) = ( g, h ) ⊕ ( e, -f ) なので

λ = ( h + f ) / ( g - e )
x = λ2 - g - e = ( h + f )2 / ( g - e )2 - ( g + e )

より

x( g - e )2=( h + f )2 - ( g + e )( g - e )2
=( h2 + 2hf + f2 ) - ( g3 - g2e - ge2 + e3 )
=( g3 + ag + b ) + ( e3 + ae + b ) + 2hf - ( g3 - g2e - ge2 + e3 )
=2hf + 2b + a( g + e ) + ge( g + e )
=2hf + 2b + ( a + ge )( g + e )

となります。二式を辺々掛け合わせて、

u( e - g )2・x( g - e )2=[ -2fh + 2b + ( a + eg )( e + g ) ][ 2hf + 2b + ( a + ge )( g + e ) ]
=[ 2b + ( a + eg )( e + g ) ]2 - ( 2fh )2
=[ 2b + ( a + eg )( e + g ) ]2 - 4( e3 + ae + b )( g3 + ag + b )
=4b2 + 4b( a + eg )( e + g ) + ( a + eg )2( e + g )2 - 4{ eg( e2 + a )( g2 + a ) + b[ e( e2 + a ) + g( g2 + a ) ] + b2 }

となって b2 を含む項は消えるので、b を含む項と含まない項にそれぞれ分けて計算します。

4b( a + eg )( e + g ) - 4b[ e( e2 + a ) + g( g2 + a ) ]
=4b[ ( ae + ag + e2g + eg2 ) - ( e3 + ae + g3 + ag ) ]
=-4b[ e2( e - g ) - g2( e - g ) ]
=-4b( e + g )( e - g )2
( a + eg )2( e + g )2 - 4eg( e2 + a )( g2 + a )
=( e + g )2( a2 + 2aeg + e2g2 ) - 4eg[ e2g2 + a( e2 + g2 ) + a2 ]
=( e2g2 + a2 )[ ( e + g )2 - 4eg ] + 2aeg[ ( e + g )2 - 2( e2 + g2 ) ]
=( e2g2 + a2 )( e - g )2 - 2aeg( e - g )2
=( e2g2 - 2aeg + a2 )( e - g )2
=( eg - a )2( e - g )2

よって、

u( e - g )2・x( g - e )2=( eg - a )2( e - g )2 - 4b( e + g )( e - g )2
ux( e - g )4=[ ( eg - a )2 - 4b( e + g ) ]( e - g )2
u=[ ( eg - a )2 - 4b( e + g ) ] / x( e - g )2

となります。これで ( x, y )i と ( x, y )i+1 から ( x, y )2i と ( x, y )2i+1 が計算できるようになりました。

( xi, yi ) = ( x, y )i とし、xi = Xi / Zi を先ほど求めた x2i の計算式に代入すると

x2i = X2i / Z2i={ [ ( Xi / Zi )2 - a ]2 - 8b( Xi / Zi ) } / 4[ ( Xi / Zi )3 + a( Xi / Zi ) + b ]
=[ ( Xi2 - aZi2 )2 - 8bXiZi3 ] / 4( Xi3Zi + aXiZi3 + bZi4 )

であり、( X, Y, Z ) とそれを定数倍した ( cX, cY, cZ ) が等しいことから

X2i = ( Xi2 - aZi2 )2 - 8bXiZi3
Z2i = 4Zi( Xi3 + aXiZi2 + bZi3 )

x2i+1 についても同様に、

x2i+1 = X2i+1 / Z2i+1={ [ ( Xi / Zi )( Xi+1 / Zi+1 ) - a ]2 - 4b[ ( Xi / Zi ) + ( Xi+1 / Zi+1 ) ] } / ( X / Z )[ ( Xi / Zi ) - ( Xi+1 / Zi+1 ) ]2
=Z[ ( XiXi+1 - aZiZi+1 )2 - 4bZiZi+1( XiZi+1 + Xi+1Zi ) ] / X( XiZi+1 - Xi+1Zi )2

より

X2i+1 = Z[ ( XiXi+1 - aZiZi+1 )2 - 4bZiZi+1( XiZi+1 + Xi+1Zi ) ]
Z2i+1 = X( XiZi+1 - Xi+1Zi )2

となります。かなり複雑な式ですが、一度コーディングしてしまえば問題はありません。また、全て乗加算とべき乗だけなので、法を決めて合同式で解くこともできます。


まずは Xi と Zi を計算するサンプル・プログラムを以下に示します。

/*
  Diff : ( u1 - u2 ) % p を解く

  u1 < u2 の場合、( u1 + Np - u2 ) % p として負数にならないようにする
*/
template< class U >
  U Diff( U u1, const U& u2, const U& p )
{
  if ( u1 < u2 ) {
    u1 += ( ( u2 - u1 ) / p + 1 ) * p;
  }

  return( ( u1 - u2 ) % p );
}

/*
  ECCoord_xSub_2i : x_2i = ( x^2 - az^2 )^2 - 8bxz^3 ( mod n ) を計算する
*/
template< class T >
T ECCoord_xSub_2i( const T& n, const T& a, const T& b, const T& x, const T& z )
{
  T t( Diff( x * x, a * z * z, n ) );

  return( Diff( t * t, 8 * b * x * z * z * z, n ) );
}

/*
  ECCoord_zSub_2i : z_2i = 4z( x^3 + axz^2 + bz^3 ) ( mod n ) を計算する
*/
template< class T >
T ECCoord_zSub_2i( const T& n, const T& a, const T& b, const T& x, const T& z )
{
  T t( ( x * x * x + a * x * z * z + b * z * z * z ) % n );

  return( ( t * 4 * z ) % n );
}

/*
  ECCoord_xSub_2i_p1 : x_2i+1 = z0[ ( x1x2 - az1z2 )^2 - 4bz1z2( x1z2 + x2z1 ) ] ( mod n ) を計算する
*/
template< class T >
T ECCoord_xSub_2i_p1( const T& n, const T& a, const T& b, const T& x1, const T& z1, const T& x2, const T& z2, const T& z0 )
{
  T t1( Diff( x1 * x2, a * z1 * z2, n ) );
  T t2( ( b * z1 * z2 * ( x1 * z2 + x2 * z1 ) ) % n );

  return( ( z0 * Diff( t1 * t1, 4 * t2, n ) ) % n );
}

/*
  ECCoord_zSub_2i_p1 : z_2i+1 = x0( x1z2 - x2z1 )^2 ( mod n ) を計算する
*/
template< class T >
T ECCoord_zSub_2i_p1( const T& n, const T& a, const T& b, const T& x1, const T& z1, const T& x2, const T& z2, const T& x0 )
{
  T t( Diff( x2 * z1, x1 * z2, n ) );

  return( ( x0 * t * t ) % n );
}

/*
  EllipticCurveCoord : 楕円曲線 y^2 = x^3 + ax + b について、( x, y, z )^k ( mod n ) の x と z を計算する
*/
template< class T >
void EllipticCurveCoord( T k, const T& n, const T& a, const T& b, T* x, T* z )
{
  T x1( *x ); // x_i
  T z1( *z ); // z_i
  T x2( ECCoord_xSub_2i( n, a, b, *x, *z ) ); // x_2i
  T z2( ECCoord_zSub_2i( n, a, b, *x, *z ) ); // z_2i

  // k をビット列で表す(最上位ビットを除いていることに注意)
  std::vector< bool > bit;
  while ( k > 1 ) {
    bit.push_back( ( k & 1 ) == 1 );
    k /= 2;
  }

  for ( std::vector< bool >::reverse_iterator i = bit.rbegin() ; i != bit.rend() ; ++i ) {
    T x3( ECCoord_xSub_2i_p1( n, a, b, x1, z1, x2, z2, *z ) ); // x_2i+1
    T z3( ECCoord_zSub_2i_p1( n, a, b, x1, z1, x2, z2, *x ) ); // z_2i+1
    // べき乗しない場合 2i, 2i+1 を採用
    if ( ! *i ) {
      T t( ECCoord_xSub_2i( n, a, b, x1, z1 ) );
      z1 = ECCoord_zSub_2i( n, a, b, x1, z1 ); // z_i ← z_2i
      x1 = t;                                  // x_i ← x_2i
      x2 = x3; // x_2i ← x_2i+1
      z2 = z3; // z_2i ← z_2i+1
    // べき乗する場合 2i+1, 2i+2 を採用
    } else {
      T t( ECCoord_xSub_2i( n, a, b, x2, z2 ) );
      z2 = ECCoord_zSub_2i( n, a, b, x2, z2 ); // z_i+1 ← z_2i+2
      x2 = t;                                  // x_i+1 ← x_2i+2
      x1 = x3; // x_i ← x_2i+1
      z1 = z3; // z_i ← z_2i+1
    }
  }

  *x = x1;
  *z = z1;
}

処理の仕方は、前章「数値演算法 (11) 素因数分解 -1-」の「4) リュカ数列 ( Lucas Sequence )」で紹介した「レーマーのアルゴリズム」のサンプル・プログラムとよく似ています。べき数 k をビット列で表し、ビットが立っていたら ( Xi, Zi ) = ( X2i+1, Z2i+1 ), ( Xi+1, Zi+1 ) = ( X2i+2, Z2i+2 ) とし、そうでなければ ( Xi, Zi ) = ( X2i, Z2i ), ( Xi+1, Zi+1 ) = ( X2i+1, Z2i+1 ) とする操作を上位側から順番に行います。( X2i, Z2i ), ( X2i+1, Z2i+1 ) の計算用に各サブルーチンを用意していますが、内容は前述の計算をそのまま実行しているだけの単純なものです。但し、n を法として計算し、結果は正数であることを前提としているので、減算処理で負数とならないよう専用の関数 Diff を用意しています。


座標の計算ができてしまえば、素因数分解への対応は比較的簡単にできます。以下にサンプル・プログラムを示します。

/*
  ECM : 楕円曲線法(ECM) による n の素因数分解

  x, y にはあらかじめ初期値が保持されているものとする

  x, y : 素因数分解した結果を返すための変数
  a : 楕円曲線の方程式 y^2 = x^3 + ax + b の係数 a
  maxCnt : 最大試行回数
*/
template< class T >
bool ECM( const T& n, T* x, T* y, const T& a, T maxCnt )
{
  T b( Diff( Diff( *y * *y, *x * *x * *x, n ), a * *x, n ) ); // b = y^2 - x^3 - ax

  // 判別式のチェック
  // 最大公約数を持つのなら解が得られたことになる
  T g( MathLib::GCD( 4 * a * a * a + 27 * b * b, n ) );
  if ( g != 1 ) {
    *x = n / g;
    *y = n / *x;
    return( true );
  }

  T z( 1 ); // Z
  T k( 2 ); // べき数

  while ( k <= maxCnt ) {
    // 座標の k 乗を 10 回計算してから最大公約数がないか確認する
    for ( unsigned int i = 0 ; i < 10 ; ++i ) {
      EllipticCurveCoord( k, n, a, b, x, &z );
      ++k;
    }
    g = MathLib::GCD( z, n );
    if ( g != 1 ) {
      *x = n / g;
      *y = n / *x;
      return( true );
    }
  }

  return( false );
}

引数として、素因数分解対象の数 n、結果を返すための x, y、楕円曲線の方程式 y2 = x3 + ax + b のパラメータ a、試行回数の最大値 maxCnt があります。xy にはあらかじめ座標の初期値を入力しておきます。通常ならば y2 = x3 + ax + b の関係式を満たす ( x, y ) の組を渡す必要がありますが、その代わりにパラメータ b をプログラム内で決定させることで任意の値が利用できるようになります。

判別式 4a3 + 27b2 と n が互いに素でない場合、n のある素因数 p について

4a3 + 27b2 ≡ 0 ( mod p )

となる可能性があります。そのため、この二数に最大公約数がないことを最初に確認しておきます。もしあれば、それが求めたい素因数の一つということになるので、その後の処理は不要となります。


ポラードの p - 1 法 ( Pollard's p - 1 Algorithm )」は、ck! ( mod p ) の中の値 c を変化させながら試行することができますが、小さな値について試して失敗した場合はより大きな値を選ぶしかありませんでした。しかし、楕円曲線法の場合は a, x, y を任意に選ぶことができるので、選択の幅が大きくなります。これを利用して、パラメータを様々に変更しながら試行錯誤するようなプログラムを以下に示します。

/*
  ECM : 楕円曲線法(ECM) による n の素因数分解

  x, y にはあらかじめ初期値が保持されているものとする
  楕円曲線の方程式 y^2 = x^3 + ax + b の係数 a を平方数に更新しながら
  maxCnt1 回試行する。1 回の試行における最大処理回数は maxCnt2 とする

  n : 素因数分解対象の数
  x, y : 素因数分解した結果を返すための変数
  a : 楕円曲線の方程式 y^2 = x^3 + ax + b の係数 a
  maxCnt1 : 最大試行回数
  maxCnt2 : 一回の試行内での最大処理回数
*/
template< class T >
bool ECM( const T& n, T* x, T* y, T a, const T& maxCnt1, const T& maxCnt2 )
{
  for ( T i( 0 ) ; i < maxCnt1 ; ++i ) {
    if ( ECM( n, x, y, a, maxCnt2 ) )
      return( true );
    a = ( a * a ) % n;
  }

  return( false );
}

a は n を法とする平方数で更新するようにしています。a を一つに固定して処理するよりこの方が高速化できるそうです。


6) 性能評価

最後に、この章で紹介した「複数多項式二次ふるい法」「連分数法」「楕円曲線法」についての性能評価をしてみたいと思います。各アルゴリズムについての計測方法は前回と同様で、ランダムに作成した素数を二つ掛け合わせた値を 10 個作成して行い、その処理時間 ( 秒 ) の中央値を求めています。桁数によって処理に失敗する場合もあるため、このときは成功時と失敗時のそれぞれに分けて中央値を算出し、併せて成功した回数も記載します。表に記載された N は素因数分解したときの素数の桁数を表します。従って、N = 3 のときに得られる合成数は、100 x 100 = 10000 ( 5 桁 ) から 999 x 999 < 1000000 ( 6 桁 ) の間となります。

まずは「複数多項式二次ふるい法」を前章での「二次ふるい法」の性能評価と同じ条件 ( 素因数の最大値 = 1000, 5000, 10000 ; F(r) の幅 500000 ; T = 1.5, 2.0 ) で処理した結果を以下に示します。

表 6-1. 複数多項式二次ふるい法 性能評価結果 ( T = 1.5 ) ( 単位 : 秒 )
素因数の最大値
1000500010000
処理時間
(成功)
処理時間
(失敗)
成功回数処理時間
(成功)
処理時間
(失敗)
成功回数処理時間
(成功)
処理時間
(失敗)
成功回数
MPQSQSMPQSQSMPQSQS
N70.5950.55891016.6-101062.315.0910
100.4030.133883.42-101012.9-1010
110.3780.140222.48-10108.66-1010
12-0.131003.860.395848.911.2489
13-0.137005.790.3545019.61.0464
 
図 6-1. 複数多項式二次ふるい法 性能評価結果グラフ ( T = 1.5 )
MPQS法性能評価結果グラフ(T=1.5)
 
表 6-2. 複数多項式二次ふるい法 性能評価結果 ( T = 2.0 ) ( 単位 : 秒 )
素因数の最大値
1000500010000
処理時間
(成功)
処理時間
(失敗)
成功回数処理時間
(成功)
処理時間
(失敗)
成功回数処理時間
(成功)
処理時間
(失敗)
成功回数
MPQSQSMPQSQSMPQSQS
N72.07-101025.2-101074.9-1010
100.4300.214787.19-101028.5-1010
110.3420.189125.80-101022.9-1010
120.4720.165107.771.107615.9-109
13-0.138006.211.122140.93.3235
 
図 6-2. 複数多項式二次ふるい法 性能評価結果グラフ ( T = 2.0 )
MPQS法性能評価結果グラフ(T=2.0)

上表の中で、成功回数については前回の「二次ふるい法」の結果を併せて載せています。QS が二次ふるい法、MPQS が今回の複数多項式二次ふるい法の成功回数を表します。また、グラフには QSMPQS の両方の処理時間を表示し、実線が MPQS、点線が QS をそれぞれ示しています。

T = 1.5, 2.0 のどちらについても、QSMPQS には処理時間・成功回数のどちらも差はないようです。そこで、F(r) の幅を 500000 から 100000 に減らした場合に差が生じるのかを調べてみました。その結果は以下のようになります。なお、T = 1.5 としています。

表 6-3. 複数多項式二次ふるい法 性能評価結果 ( F(r) の幅 100000, T = 1.5 ) ( 単位 : 秒 )
素因数の最大値
1000500010000
処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数
QSN70.245-102.69-1011.1-10
100.1230.082211.41-102.94-10
11-0.071201.160.38852.411.067
12-0.07550-0.2800-0.7810
13-0.07240-0.2460-0.5420
MPQS70.213-103.80-109.30-10
100.3040.084441.230.41784.421.229
11-0.081302.200.37765.480.9527
12-0.082201.630.30019.150.7643
13-0.083102.670.27013.600.5811
 
図 6-3. 複数多項式二次ふるい法 性能評価結果グラフ ( F(r) の幅 100000, T = 1.5 )
MPQS法性能評価結果グラフ(F(r)の幅100000,T=1.5)

こちらは成功回数において MPQS の方がほんの少しだけ多くなっているようにみえますが、誤差の範囲内かもしれません。この結果を見る限り、F(r) の幅が小さくなったことによる効果はそれほど大きくないように思います。

MPQS の一番の利点は多項式の係数をいろいろと変えて処理することができるところです。今回のサンプル・プログラムでは、計算によって係数を一つだけ決めて処理するようにしていますが、特に一つに決めなければならない理由はなく、いくつかのパターンについて処理させることも可能です。さらに、それぞれの係数について並列で処理できるのであれば、より効率よく素因数を見つけることも可能になります。そのようにして初めて、複数多項式二次ふるい法は真の能力を発揮することができるのだと言えます。


続いては「連分数法」です。因数基地の素因数の最大値を 200, 400, 1000, 10000 として処理時間と成功回数を計測しました。

表 6-4. 連分数法 性能評価結果 ( 単位 : 秒 )
素因数の最大値
200400100010000
処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数
N30.003850.025390.007450.0066080.02700.031890.1060.0860759
50.007480.016280.02820.053980.07190.066390.731-10
70.06330.088480.0793-100.171-107.72-10
100.927-100.513-100.834-1011.0-10
112.00-102.17-101.64-1012.3-10
129.28-103.70-104.13-1017.8-10
1320.6-1013.2-1011.4-1031.6-10
14114.254-1043.7-1022.9597-1046.2-10
15126-1039.8396-1096.6-10
 
図 6-4. 連分数法 性能評価結果グラフ
連分数法性能評価結果グラフ

N = 3, 5, 7 のときに失敗する場合がありました。逆に大きな桁数では二次ふるい法よりも成功率が高いという結果になっています。N を 15 まで大きくしてみましたが、成功率は 100% です。しかし、処理時間は桁数が上がるごとに増えていきます。因数基地の素因数の最大値は、30 桁 ( N = 15 ) 程度の数を素因数分解する場合は今回試した中で 1000 くらいがちょうど適しているようです。これより少ないと Cn が素因数分解できる確率は減り、逆に大きいと試し割りに時間がかかるようになります。

連分数法は、試し割りを行う数の大きさを抑えることができる代わりに、二次ふるい法のように候補をあらかじめ絞り込むというようなことはできません。複数多項式二次ふるい法のようにパラメータを振り分けて並列処理するようなこともできないため、処理時間としては二次ふるい法の方が速いとされています。


最後に「楕円曲線法」の結果を以下に示します。まずは、楕円曲線の方程式 y2 = x3 + ax + b のパラメータ a を 2, 7, N / 2 に固定したとき ( N は素因数分解対象の数 ) の計測結果です。x, y の初期値は、それぞれ N / 3, 2N / 3 としています。また、最大処理回数は 100000 に制限しました。

表 6-5. 楕円曲線法 性能評価結果 - パラメータを固定した場合 - ( 単位 : 秒 )
a
27N / 2
処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数
N70.66252.990.32748.770.11262.49
105.1758.625.9661.740.94762.55
114.7564.2219.465.0415.767.72
1210.665.049.3266.1211.768.64
1325.368.0416.268.3470.473.01
 
図 6-5. 楕円曲線法 性能評価結果グラフ - パラメータを固定した場合 -
ECM性能評価結果グラフ(a固定)

a に対する差は特には見られません。7 桁でも失敗することがあることから、パラメータを固定して処理するのはあまりよくないようです。失敗したときの処理時間が最長となり、かつ最大処理回数に依存することは処理の内容からみて明らかです。

a を平方数で更新するサンプル・プログラムを使った計測結果は以下のようになりました。あるパラメータでの処理回数の最大を 1000 回とし、最大 100 回までパラメータを更新しながら処理していきます。従って、延べの最大処理回数は a を固定した場合と同程度となります。

表 6-6. 楕円曲線法 性能評価結果 - パラメータを更新 1 - ( 単位 : 秒 )
a
27N / 2
処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数
N70.361-100.263-100.718-10
101.75-106.59-103.75-10
119.14-104.90-1010.7-10
1210.039.9615.639.285.0238.29
1320.440.2720.340.9511.840.96
 
図 6-6. 楕円曲線法 性能評価結果グラフ - パラメータを更新 1 -
ECM性能評価結果グラフ(a 100回更新)

パラメータをいろいろ切り替えながら処理することによる効果はあるようです。一つのパラメータでの処理をどの程度で打ち切るかは難しいところですが、成功したときと失敗したときの処理時間の差を見る限り、早めに打ち切ってパラメータの更新回数の方を増やした方が有利なようにみえます。そこで、あるパラメータでの処理回数の最大を 100 回とし、最大 1000 回までパラメータを更新しながら処理したところ、次のような結果になりました。

表 6-7. 楕円曲線法 性能評価結果 - パラメータを更新 2 - ( 単位 : 秒 )
a
27N / 2
処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数
N70.197-100.134-100.194-10
102.4923.393.0823.093.0523.09
1115.924.445.3024.3616.424.55
129.5324.949.6225.223.6725.12
13-29.30-26.60-26.40
 
図 6-7. 楕円曲線法 性能評価結果グラフ - パラメータを更新 2 -
ECM性能評価結果グラフ(a 1000回更新)

残念ながら、失敗したときの時間は短くなるものの、成功回数が減少してしまうという結果になりました。逆に、あるパラメータでの処理回数の最大を 10000 回とし、最大 10 回までパラメータを更新しながら処理させると次のようになります。

表 6-8. 楕円曲線法 性能評価結果 - パラメータを更新 3 - ( 単位 : 秒 )
a
27N / 2
処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数処理時間(成功)処理時間(失敗)成功回数
N70.562-100.239-101.11-10
108.70-105.9148.793.52-10
1122.0-1018.0-1021.653.08
1219.354.7812.753.9722.5-10
1332.454.5521.256.8625.657.23
 
図 6-8. 楕円曲線法 性能評価結果グラフ - パラメータを更新 3 -
CM性能評価結果グラフ(a 10回更新)

成功回数は最初の結果とほとんど変わらず、処理時間だけ長くなってしまうという結果になりました。最初のパラメータが最適だったようです。


補足 1) 合同式 x2 ≡ a ( mod pn ) の解法

p を奇素数とします。合同式 x2 ≡ a ( mod pn ) が二つの解 x ≡ ±t ( mod pn ) を持つと仮定すると、ある整数 m を使って

a - t2 = mpn

と表すことができます。さらに合同式 x2 ≡ a ( mod pn+1 ) が二つの解 x ≡ ±s ( mod pn+1 ) を持つとき、ある整数 m' を使って

a - s2 = m'pn+1 = (m'p)pn

と表され、

x ≡ ±s ( mod pn )

でもあるので、整数 k を用いて

s = t + kpn

と表すことができます。これを a - s2 = m'pn+1 に代入すると、

a - ( t + kpn )2 = a - t2 - 2ktpn + kp2n = m'pn+1

より

a - t2 ≡ 2ktpn ( mod pn+1 )

を満たす必要があります。よって、x2 ≡ a ( mod pn ) の解 t が得られたら、t と p は互いに素であることに注意して

kpn ≡ ( a - t2 ) / 2t ( mod pn+1 )

を満たすような kpn を求めれば、s = t + kpn が x2 ≡ a ( mod pn+1 ) の解となります。p = 1 のとき、x2 ≡ a ( mod p ) の解法はすでにわかっているので、解 t が得られたら

kp ≡ ( a - t2 ) / 2t ( mod p2 )

より kp を計算して s = t + kp を求めれば、それが x2 ≡ a ( mod p2 ) の解です。


補足 2) バスカラ-ブラウンカーのアルゴリズムの証明 (1)

バスカラ-ブラウンカーのアルゴリズムの漸化式

B0 = 0 ; C0 = 1 ; P0 = 1 ; Q0 = 0
A1 = s ; B1 = s ; C1 = m - s2 ; P1 = s ; Q1 = 1
An+1 = ⌊ ( s + Bn ) / Cn---(1.1)
Bn+1 = An+1Cn - Bn---(1.2)
Cn+1 = Cn-1 + An+1( Bn - Bn+1 )---(1.3)
Pn+1 = Pn-1 + An+1Pn---(1.4)
Qn+1 = Qn-1 + An+1Qn---(1.5)

を使い、まずは

√m = A1 + 1 / [ A2 + 1 / [ A3 + 1 / [ A4 + 1 / [ A5 + 1 / [ A6 + ... --- (1.6)

を証明します。

En+1 ≡ ( √m + Bn ) / Cn --- (1.7)

とすると、

En - 1 / En+1=( √m + Bn-1 ) / Cn-1 - Cn / ( √m + Bn )
=[ ( √m + Bn-1 )( √m + Bn ) - Cn-1Cn ] / Cn-1( √m + Bn )
=[ ( m - Cn-1Cn ) + ( Bn-1 + Bn )√m + Bn-1Bn ] / Cn-1( √m + Bn )

と変形できます。ここで、

Bn+12 + CnCn+1=( An+1Cn - Bn )2 + Cn[ Cn-1 + An+1( Bn - Bn+1 ) ]
=An+12Cn2 - 2An+1BnCn + Bn2 + Cn( Cn-1 + An+1Bn - An+1Bn+1 )
=An+12Cn2 - An+1BnCn - An+1Bn+1Cn + Bn2 + CnCn-1
=An+1Cn( An+1Cn - Bn ) - An+1Bn+1Cn + Bn2 + CnCn-1
=An+1Bn+1Cn - An+1Bn+1Cn + Bn2 + CnCn-1
=Bn2 + CnCn-1

であり、m = 0 のとき

B12 + C1C0 = s2 + ( m - s2 ) = m

なので、帰納法により

Bn2 + CnCn-1 = m --- (1.8)

となります。(1.8) を先ほどの式に代入すると

En - 1 / En+1=[ Bn2 + ( Bn-1 + Bn )√m + Bn-1Bn ] / Cn-1( √m + Bn )
=( Bn + Bn-1 )( √m + Bn ) / Cn-1( √m + Bn )
=( Bn + Bn-1 ) / Cn-1

となって、(1.2) より

An+1 = ( Bn+1 + Bn ) / Cn

と比較すれば

En - 1 / En+1 = An より

En = An + 1 / En+1 --- (1.9)

という結果が得られます。従って、

√m = E1=A1 + 1 / E2
=A1 + 1 / [ A2 + 1 / E3 ]
=...
=A1 + 1 / [ A2 + 1 / [ ... + 1 / [ An-1 + 1 / En ] ... ] ]
=...
=A1 + 1 / [ A2 + 1 / [ A3 + 1 / [ A4 + ...

となり (1.6) が証明できました。

次に、この結果を利用して

Pn / Qn = A1 + 1 / [ A2 + 1 / [ A3 + ... + 1 / [ An-1 + 1 / An ] ... ] ] --- (1.10)

が成り立つことを帰納法を使って証明します。n = 1 のときは、P1 = s = A1, Q1 = 1 より

P1 / Q1 = A1

n = 2 のときは

P2 = P0 + A2P1 = A1A2 + 1

Q2 = Q0 + A2Q1 = A2

より

P2 / Q2 = A1 + 1 / A2

でどちらも成り立っています。そこで、

Pn / Qn = A1 + 1 / [ A2 + 1 / [ A3 + ... + 1 / [ An-1 + 1 / An ] ... ] ]

が成り立っていると仮定します。(1.4) (1.5) から

Pn / Qn = ( Pn-2 + AnPn-1 ) / ( Qn-2 + AnQn-1 )

でもあるので、

A1 + 1 / [ A2 + 1 / [ A3 + ... + 1 / [ An-1 + 1 / [ An + 1 / An+1 ] ] ... ] ]
=[ Pn-2 + ( An + 1 / An+1 )Pn-1 ] / [ Qn-2 + ( An + 1 / An+1 )Qn-1 ]
=[ ( Pn-2 + AnPn-1 ) + Pn-1 / An+1 ] / [ ( Qn-2 + AnQn-1 ) + Qn-1 / An+1 ]
=( Pn + Pn-1 / An+1 ) / ( Qn + Qn-1 / An+1 )
=( PnAn+1 + Pn-1 ) / ( QnAn+1 + Qn-1 )
=Pn+1 / Qn+1

となり、(1.10) が成り立つことが証明されました。同様のやり方で、

√m=A1 + 1 / [ A2 + 1 / [ A3 + ... + 1 / [ An-1 + 1 / [ An + 1 / En+1 ] ] ... ] ]
=[ Pn-2 + ( An + 1 / En+1 )Pn-1 ] / [ Qn-2 + ( An + 1 / En+1 )Qn-1 ]
=( PnEn+1 + Pn-1 ) / ( QnEn+1 + Qn-1 )

と変形できるので、

√m( QnEn+1 + Qn-1 ) = PnEn+1 + Pn-1

に En+1 の定義式 (1.7) を代入して

√m[ Qn( √m + Bn ) / Cn + Qn-1 ] = Pn( √m + Bn ) / Cn + Pn-1

より

mQn + ( BnQn + CnQn-1 )√m = ( BnPn + CnPn-1 ) + Pn√m

両辺の係数を比較して

mQn = BnPn + CnPn-1

Pn = BnQn + CnQn-1

よって、

Pn2 - mQn2=Pn( BnQn + CnQn-1 ) - Qn( BnPn + CnPn-1 )
=Cn( PnQn-1 - Pn-1Qn )

という結果が得られます。PnQn-1 - Pn-1Qn は n = 1 のとき

P1Q0 - P0Q1 = -1

です。PnQn-1 - Pn-1Qn = ( -1 )n が成り立つと仮定すると、

Pn+1Qn - PnQn+1=( Pn-1 + An+1Pn )Qn - Pn( Qn-1 + An+1Qn )
=Pn-1Qn - PnQn-1
=-( -1 )n = ( -1 )n+1

となって、帰納法により PnQn-1 - Pn-1Qn = ( -1 )n が成り立つことが確認できて、

Pn2 - mQn2 = ( -1 )nCn --- (1.11)

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

補足 3) バスカラ-ブラウンカーのアルゴリズムの証明 (2)

An, Bn, Cn は全て有限の値を取り、その範囲は

0 < An < 2√m
0 < Bn < √m
0 < Cn < 2√m

となります。これを以下で証明します。

An+1 は Bn と Cn のみで決まり、Bn と Cn がともに正であれば必ず正の整数になることに注意して下さい。B1 = ⌊√m⌋, C1 = m - ( ⌊√m⌋ )2 なので、少なくとも A2 までは正の整数となりますが、m が平方数の場合は C1 = 0 となってしまうため、漸化式として成り立たなくなります。よって、m は平方数ではないとします。

まずは、(1.1) (1.7) の式

An+1 = ⌊ ( ⌊√m⌋ + Bn ) / Cn ⌋ --- (1.1)
En+1 = ( √m + Bn ) / Cn --- (1.7)

において

0 < En+1 - An+1 < 1

となることの証明です。最初に

An+1 = ⌊En+1

を証明します。An+1 ≤ ⌊En+1⌋ なのは明らかで、どちらも整数となるので、An+1 + 1 > ⌊En+1⌋ であることが示されれば十分です。

An+1 + 1=⌊ ( ⌊√m⌋ + Bn ) / Cn ⌋ + 1
=⌊ ( ⌊√m⌋ + Bn ) / Cn + 1 ⌋
=⌊ ( ⌊√m⌋ + Bn + Cn ) / Cn

ここで Cn > 0 ならば ⌊√m⌋ + Cn > √m なので

An+1 + 1 > ⌊ ( √m + Bn ) / Cn ⌋ = ⌊En+1

となって、An+1 = ⌊En+1⌋ が示されました。よって、En+1 < An+1 + 1 であり、An+1 < En+1 は明らかなので、

0 < En+1 - An+1 < 1 --- (2.1)

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

次に、少なくとも 0 ≤ k ≤ n においては

0 < Bk+1 < √m
0 < Ck < 2√m

がともに成り立つと仮定します。n = 0 なら

0 < B1 = ⌊√m⌋ < √m
0 < C0 = 1 < 2√m

となるので実際に成り立っています。Bn+1, Cn までがそれぞれ成り立っていると仮定していることに注意して下さい。

補足 1 の (1.9) より

An = En - 1 / En+1

であることがわかっているので、(1.7) と併せて

0 < En+1 - An+1 = 1 / En+2 = Cn+1 / ( √m + Bn+1 )

となり、Bn+1 > 0 と仮定しているので Cn+1 > 0 が成り立ちます。

上式は、n を一つ増やして

1 / En+3 = Cn+2( √m - Bn+2 ) / ( m - Bn+22 )

と変形することができます。補足 1 にて (1.8)

Bn+12 + CnCn+1 = m

が成り立つことを証明していたので これを代入することで

1 / En+3 = ( √m - Bn+2 ) / Cn+1

となります。よって、

0 < Cn+1 / En+3 = √m - Bn+2

より Bn+2 < √m となります。さらに

0 < En+1 - An+1 = 1 / En+2 < 1

より En+2 > 1 となることを利用すると、

1 < En+2 = ( √m + Bn+1 ) / Cn+1 < 2√m / Cn+1

となって、Cn+1 < 2√m が証明されます。

Bn の漸化式 (1.2) から

Bn+2 = An+2Cn+1 - Bn+1

について、Bn+2 ≤ 0 になると仮定すると

An+2Cn+1 ≤ Bn+1

仮定から Bn+1 < √m だったので、

An+2Cn+1 ≤ Bn+1 < √m

となります。すでに Cn+1 > 0 が成り立つことは証明されていて、仮定から Bn+1 > 0 でもあるので An+2 > 0 となることに注意すると、

Cn+1 ≤ An+2Cn+1 ≤ Bn+1 < √m

より Cn+1 < √m となります ( これは Bn+2 ≤ 0 になると仮定したときに導き出された関係式であることに注意して下さい )。

1 > En+2 - An+2 = 1 / En+3 = ( √m - Bn+2 ) / Cn+1

となりますが、Bn+2 ≤ 0 を仮定しているので

( √m - Bn+2 ) / Cn+1 ≥ √m / Cn+1

であり、Cn+1 < √m の関係式から

√m / Cn+1 > 1

です。従って、 1 > 1 となって矛盾することから、Bn+2 > 0 が示されました。

最後に、

An < En = ( √m + Bn+1 ) / Cn+1 < 2√m / Cn+1 ≤ 2√m

であり、これで全ての関係式が証明されたことになります。


補足 4) ペル方程式が解を持つことの証明

バスカラ・ブラウンカーのアルゴリズムがペル方程式 Pn2 - mQn2 = (-1)nCn を満たすことは証明しましたが、まだその中に Pn2 - mQn2 = 1 を満たす解が存在するかどうかはわかっていません。これを証明したのはイタリアの数学者・天文学者「ジョセフ・ルイ・ラグランジュ ( Joseph-Louis Lagrange )」でした。その中では「鳩の巣原理 ( Pigeonhole Principle )」という証明法が使われています。

まず、バスカラ・ブラウンカーのアルゴリズムの Pn, Qn に対する漸化式

Pn = Pn-1 + An+1Pn
Qn = Qn-1 + An+1Qn

に着目すると、An+1 > 0 だったので Pn > Pn-1, Qn > Qn-1 が必ず成り立ち、無数の解の組 ( Pn, Qn ) が生成されることに着目します。これに対し、0 < Cn < 2√m よりペル方程式の右辺は有限な数しか存在しません。右辺の取りうる数 -⌊2√m⌋ から ⌊2√m⌋ を鳩の巣に見立て、この中に鳩として解 ( Pn, Qn ) を入れることを考えると、鳩は無数にいるので、鳩の巣の少なくともひとつは鳩である解 ( Pn, Qn ) が無数に存在することになります。そのような鳩の巣の一つに着目します。

次に、この無数の解を対応する Cn の剰余で分類します。例えば、Cn = 5 だったとき、( Pn, Qn ) = ( 12, 6 ) という解は ( 2, 1 ) として分類されます。ここでも先ほどと同様な「鳩の巣原理」の考え方により、剰余による分類は有限な数しかないので、少なくとも一つの Cn に対して ( Pn, Qn ) は無数にあります。そのような解の組から二つを抽出したとき、常に

Pi ≡ Pj ( mod Cn )
Qi ≡ Qj ( mod Cn )

を満たします。

抽出した二つの解を使い、整数 u, v を次のように定義します。

u + v√m = ( Pi + Qi√m )( Pj - Qj√m )
u - v√m = ( Pi - Qi√m )( Pj + Qj√m )

このとき、

u + v√m=( Pi + Qi√m )( Pj - Qj√m )
=( PiPj - mQiQj ) + ( PjQi - PiQj )√m
u - v√m=( Pi - Qi√m )( Pj + Qj√m )
=( PiPj - mQiQj ) - ( PjQi - PiQj )√m

となるので、

u = PiPj - mQiQj
v = PjQi - PiQj

として一致することに注意して下さい。辺々掛け合わせると

u2 - mv2=( u + v√m )( u - v√m )
=( Pi + Qi√m )( Pj - Qj√m )( Pi - Qi√m )( Pj + Qj√m )
=( Pi2 - mQi2 )( Pj2 - mQj2 ) = Cn2

となるので、( u, v ) もペル方程式の解の一つです。ここで、ゼロでない整数 k, l を使って

Pj = Pi + kCn
Qj = Qi + lCn

と表せることから

u=Pi( Pi + kCn ) - mQi( Qi + lCn )
=( Pi2 - mQi2 ) + ( Pik - mQil )Cn
Pi2 - mQi2 ( mod Cn )
0 ( mod Cn )
v=( Pi + kCn )Qi - Pi( Qi + lCn )
=( PiQi - PiQi ) + ( kQi - lPi )Cn
0 ( mod Cn )

となって、u, v は Cn で割り切れます。a = u / Cn, b = v / Cn とすれば

a2 - mb2 = ( u2 - mv2 ) / Cn2 = 1

より目的のペル方程式が得られます。


補足 5) 楕円曲線の判別式

楕円曲線 E : y2 = x3 + ax + b に対し、

Δ(E) = 4a3 + 27b2

を E の「判別式 ( Discriminant )」といい、Δ( E ) = 0 のとき y = 0 の解に重解が含まれるようになります。一般に、n 次の多項式 Fn : f(x) の判別式 Δ( Fn ) は、f(x) = 0 の解を αk ( k = 1, 2, ... n ) としたとき

Δ( Fn )=( α1 - α2 )2( α1 - α3 )2...( α1 - αn )2( α2 - α3 )2...( αn-1 - αn )2
=Πk{1→n ; j<k}( ( αj - αk )2 )

と表されます。二次式なら

Δ( F2 )=( α1 - α2 )2
=( α1 + α2 )2 - 4α1α2

となりますが、二次式を f(x) = ax2 + bx + c ( a ≠ 0 ) としたときに

ax2 + bx + c=a( x - α1 )( x - α2 )
=ax2 - a( α1 + α2 )x + aα1α2

となるので、係数を比較して

α1 + α2 = -b / a
α1α2 = c / a

が成り立ちます。よって、

Δ( F2 )=( -b / a )2 - 4( c / a )
=( b2 - 4ac ) / a2

となって、b2 - 4ac = 0 のとき二次方程式は重解を持つことになります。

三次式の場合、

Δ( F3 ) = ( α1 - α2 )2( α2 - α3 )2( α3 - α1 )2

となり、これを係数に置き換えるのは少々やっかいです。まず、三次式は

f(x) = x3 + ax + b

の形に限定します。解と係数の関係は、

( x - α1 )( x - α2 )( x - α3 )
=[ x2 - ( α1 + α2 )x + α1α2 ]( x - α3 )
=x3 - ( α1 + α2 + α3 )x2 + ( α1α2 + α2α3 + α3α1 )x - α1α2α3
=x3 + ax + b

より

α1 + α2 + α3 = 0
α1α2 + α2α3 + α3α1 = a
α1α2α3 = -b

となり、これを利用すると

( α1 - α2 )( α2 - α3 )
=α1α2 - α3α1 - α22 + α2α3
=α2( α3 + α1 ) - α3α1 - α22
=22 - ( a - α1α2 - α2α3 ) - α22
=-2α22 - [ a - α2( α3 + α1 ) ]
=-2α22 - ( a + α22 )
=-( 3α22 + a )

と変形できます。同様に、

( α2 - α3 )( α3 - α1 ) = -( 3α32 + a )
( α3 - α1 )( α1 - α2 ) = -( 3α12 + a )

となるので、Δ( F3 ) は

Δ( F3 )=( α1 - α2 )2( α2 - α3 )2( α3 - α1 )2
=-( 3α12 + a )( 3α22 + a )( 3α32 + a )
=-[ 9α12α22 + 3a( α12 + α22 ) + a2 ]( 3α32 + a )
=-[ 27α12α22α32 + 9a( α12α22 + α22α32 + α32α12 ) + 3a2( α12 + α22 + α32 ) + a3 ]
=-[ 27b2 + 9a( α12α22 + α22α32 + α32α12 ) + 3a2( α12 + α22 + α32 ) + a3 ]

とすることができます。

α12α22 + α22α32 + α32α12
=( α1α2 + α2α3 + α3α1 )2 - 2α12α2α3 - 2α1α22α3 - 2α1α2α32
=a2 - 2α1α2α3( α1 + α2 + α3 ) = a2
α12 + α22 + α32
=( α1 + α2 + α3 )2 - 2( α1α2 + α2α3 + α3α1 ) = -2a

より

Δ( F3 )=-[ 27b2 + 9a・a2 - 3a2・2a + a3 ]
=-( 4a3 + 27b2 )

となって、係数に置き換えることができました。


補足 6) | E( a, b ) / p | = | ( x, y ) | の証明 ( ラグランジュの定理 )

E( a, b ) / p には、p を法とする全ての整数解 ( x, y ) が含まれており、p が法となるのでその個数 | E( a, b ) / p | は有限となるのでした。この中からある整数解 ( x0, y0 ) を選び、( x1, y1 ) = ( x0, y0 ) ⊕ ( x0, y0 ) から始めて ( x2, y2 ) = ( x1, y1 ) ⊕ ( x0, y0 ), ( x3, y3 ) = ( x2, y2 ) ⊕ ( x0, y0 ) ... と次々と ( x0, y0 ) との演算結果を求めていくと、| E( a, b ) / p | は有限なのでやがて演算結果は ∞ となり、次の演算で ( x0, y0 ) に戻ります。すなわち ( x0, y0 )n = ∞ となる n が存在し、n の最小数を | ( x0, y0 ) | とするのでした。ここで、k ≠ m ならば ( x0, y0 )k ≠ ( x0, y0 )m となることに注意して下さい。もし、( x0, y0 )k = ( x0, y0 )m ( 但し n > k > m > 0 ) が成り立つと仮定すると、両辺に ( x0, y0 )n-k を演算すると

( x0, y0 )k ⊕ ( x0, y0 )n-k = ∞
( x0, y0 )m ⊕ ( x0, y0 )n-k = ( x0, y0 )n-(k-m)

が等しいことになって、n が ( x0, y0 )n = ∞ となる最小数であるという仮定に矛盾します。

| E( a, b ) / p | = | ( x0, y0 ) | のとき、全ての点が演算結果の中に含まれていることになるので、これで証明は終わったことになります。従って、| E( a, b ) / p | > | ( x0, y0 ) | と仮定します。このとき、演算結果に含まれなかった点 ( x'0, y'0 ) が存在することになります。この点に ( x0, y0 ), ( x0, y0 )2, ... ( x0, y0 )n で演算を行います。もし、この中のある点 ( x'0, y'0 ) ⊕ ( x0, y0 )k と ( x'0, y'0 ) ⊕ ( x0, y0 )m が等しいとすれば ( x0, y0 )k と ( x0, y0 )m が等しいことになりますが、先に示した通り、このときは k = m になります。従って、演算した結果は全て異なる点であることになります。また、最初の集合と、新たに作成された集合の間で同じ点が存在する、すなわち ( x'0, y'0 ) ⊕ ( x0, y0 )k と ( x0, y0 )m が等しいと仮定すると、( x0, y0 )n-k を演算することで

( x'0, y'0 ) ⊕ ( x0, y0 )k ⊕ ( x0, y0 )n-k = ( x'0, y'0 )
( x0, y0 )m ⊕ ( x0, y0 )n-k = ( x0, y0 )m+n-k

となって、( x'0, y'0 ) が最初の集合に含まれることとなり仮定に矛盾します。これによって、ちょうど n 個の新たな集合が生成されたことになります。これで点が残っていなければ証明は完了で、残っているのなら同じ操作を繰り返せば、やはり新しい n 個の点が生成できます。この新しい点はやはり重複がなく、前の二つの集合にある点とも異なるものになります。このように、全ての点が見つかるまでこの操作は続けることができます。また、n 個より少ない点が残るということはありません。もしそうなれば、( x0, y0 ), ( x0, y0 )2, ... ( x0, y0 )n によってちょうど n 個の点が生成されるという事実と矛盾するからです。

このようにして、E( a, b ) / p の点は全て n 個ずつの集合に分けることができます。これを「類別 ( Classification )」といいます。このようにしてちょうど K 個の集合に類別できたとすると、

| E( a, b ) / p | = K・| ( x, y ) |

が成り立ちます。これを一般化した定理は「ラグランジュの定理 ( Lagrange's Theorem )」といい、「群論 ( Group Theory )」における重要な定理の一つとなっています。


<参考文献>
  1. 「素因数分解と素数判定」 David M. Bressoud 著 (SiB access)
  2. Wikipedia

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