数値演算法

(6) 素数判定法

RSA暗号」には公開鍵として二つの素数が必要となることを以前説明しました。サンプル・プログラムの中では非常に小さな数を使って公開鍵を作成したので、素数であるか判別するためには実際に素因数分解できるか試していたのに対し、実際の処理においては数百桁もの数値を扱うため素因数分解を使って素数であるかを判断することは不可能になります ( そもそも RSA 暗号は、素因数分解が困難であることを前提にした暗号です )。実は、素数であるかどうかを判定することは、素因数分解を行わなくても可能です。素数でないことが分かったとしても、それを素因数に分解することはできないというのも不思議な気がしますが、例えば以前紹介した「フェルマーの小定理」を使えば、合成数であることは素因数分解をしなくても判断可能です。
この章では、素数であるかどうかを判定する方法について紹介したいと思います。

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

1) エラトステネスのふるい

素数であるかどうかを判定する最も簡単な方法は、判定する数より小さな値全てを使って除算を行い、その剰余を調べることです。ある数 N に対して 2 から N - 1 までの数で割ってみて、一つでも剰余が 0 になる数が見つかったら N は合成数であると判定することができます。

上記方法のサンプル・プログラムを以下に示します。

/*
  PowerRoot : ニュートン - ラフソン法を用いてべき乗根以上の整数最小値を求める

  c : べき乗根を求める数値
  e : 指数
*/
BigNum::Unsigned PowerRoot( BigNum::Unsigned c, const BigNum::Unsigned& e )
{
  c <<= 64; // 64bit ゲタ履き
  BigNum::Unsigned ans( c );

  while ( true ) {
    BigNum::Unsigned fx( ans );
    fx.pow( e );
    fx -= c;
    BigNum::Unsigned dfx( ans );
    dfx.pow( e - 1 );
    dfx *= e;
    if ( fx < dfx ) break;

    ans -= fx / dfx;
  }
  // ゲタ履き分を元に戻す
  ans >>= 32;
  c >>= 64;

  BigNum::Unsigned n( ans );
  if ( n.pow( e ) != c )
    ans += 1;

  return( ans );
}

/*
  IsPrime : 試し割りによる数 prime の素数判別
*/
bool IsPrime( const BigNum::Unsigned& prime )
{
  if ( prime < 2 ) return( false );
  if ( prime == 2 ) return( true );
  if ( ( prime & 1 ) == 0 ) return( false );

  BigNum::Unsigned maxNum( PowerRoot( prime, 2 ) );
  for ( BigNum::Unsigned div( 3 ) ; div < maxNum ; div += 2 )
    if ( prime % div == 0 )
      return( false );

  return( true );
}

試し割りを行う必要のある数は、判定対象の数の平方根以下の範囲で充分です。なぜなら、平方根を越える素因数が存在する場合、平方根より小さい素因数も必ず存在するからです。平方根の算出には、前章で紹介した「ニュートン - ラフソン法」を利用しています。また、計算には、以前作成した「多倍長整数クラス ( BigNum::Unsigned )」を利用しています。
試し割りは確実に素数判定を行うことができますが、桁数の大きな数に対しては、判定するまでにかなりの時間がかかります。

ある範囲の値の中から全ての素数を得るための方法として「エラトステネスのふるい ( Sieve of Eratosthenes )」があります。「エラトステネス ( Eratosthenes )」は古代ギリシャの学者で、夏至の日の南中高度が場所によって異なることを利用して、地球の大きさを概算ながら求めたことでも有名です。

エラトステネスのふるいのアルゴリズムは次のようになります。

  1. 整数を 2 から昇順で配列に並べる。また、素数リストとして空の配列を用意する。
  2. 配列の先頭にある数を素数リストに追加する。
  3. 追加した数の倍数を配列から消去する。
  4. 配列にある要素の最大値が素数リストの最大値の平方より小さければ、配列に残った要素は全て素数なので、
    素数リストに追加して処理を終了する。そうでなければ 2 に戻る。
素数リスト
配列        2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
  ↓
素数リスト  2
配列           3     5     7     9    11    13    15    17    19
  ↓
素数リスト  2  3
配列                 5     7          11    13          17    19
  ↓
素数リスト  2  3     5
配列                       7          11    13          17    19
  ↓
素数リスト  2  3     5     7          11    13          17    19
配列

素数リストにある素数の倍数は全て配列から除外されています。従って、もし配列の中に合成数が存在するとしたら、それは素数リストの最大値より大きな素因数を持っていることになります。つまり、それは素数リストの最大値の平方より大きいことになり、もし、配列の要素の最大値が素数リストの最大値の平方より小さければ、配列の中には合成数は残っていないということになります。アルゴリズムの 4 番めの内容が成り立つのはそのためです。

エラトステネスのふるいを使った素数表作成用のサンプル・プログラムを以下に示します。

/*
  Eratosthenes : エラトステネスのふるい

  maxNum : 素数表の最大値
  result : 素数表へのポインタ
*/
template< class T >
void Eratosthenes( T maxNum, vector< T >* result )
{
  vector< bool > isPrime( maxNum, true ); // 素数なら true (最初は全数値を true に初期化)
  result->clear();
  T num; // 素数探索開始位置 + 1

  for ( num = 2 ; num <= maxNum ; ++num ) {
    // 素数リストの最大値の平方が数列の最大値を越えたら、数列には素数しか残っていない
    if ( result->size() > 0 ) {
      T i( result->back() * result->back() );
      if ( i > maxNum ) break;
    }
    // 素数を探索して素数リストへ追加
    if ( ! isPrime[num - 1] ) continue;
    result->push_back( num );
    // 追加した素数の倍数を全て除去
    for ( T i = num * num ; i <= maxNum ; i += num )
      isPrime[i - 1] = false;
    // 数列に残った要素の最大値を探索
    while ( ! isPrime[maxNum - 1] )
      if ( --maxNum < 2 ) break;
  }
  // 数列に残った素数を素数リストへコピー
  while ( num <= maxNum ) {
    if ( isPrime[num - 1] )
      result->push_back( num );
    ++num;
  }
}

整数を並べた配列を実際に用意する代わりに、ここでは素数か合成数かを判定する bool値の配列を使い、値を除去する代わりに falseにしています。常にインデックスが数値を示すようにしておきたい他に、vectorに対して要素の追加・除去を行った場合の効率の悪さを避けるための処置になります。
素数の倍数を除去する処理は、対象の素数の平方から始めています。それより小さい値は、より小さな素数による処理によってすでに除去されていることが分かっているためです。

ところで、前の章でゼータ関数を紹介しました。ゼータ関数は、次のように表される関数です。

ζ(s) = Σk{1→∞}( 1 / ks ) = 1 + 1/2s + 1/3s + 1/4s + ...

この式の両辺に 2s を掛けてみます。

2sζ(s)=2s + (2/2)s + (2/3)s + (2/4)s + (2/5)s + (2/6)s + ...
=( 1 + (1/2)s + (1/3)s + ... ) + 2s( 1 + (1/3)s + (1/5)s + ... )

分母が 2 の倍数のべき乗である項は 2s で約分して前半に集め、それ以外は 2s で括って後半に集めています。前半部は ζ(s) になるので、上式は次のように変形することができます。

( 1 - 1/2s )ζ(s) = 1 + (1/3)s + (1/5)s + (1/7)s + (1/9)s + (1/11)s + ...

左辺には、分母が全て奇数のべき乗である項だけが残ります。これはちょうど、エラトステネスのふるいで、最初の素数である 2 を素数リストに追加して、配列から 2 の倍数を除外したこととよく似ています。エラトステネスのふるいの処理に倣い、次の素数である 3 の s 乗を掛けてみましょう。

3s( 1 - 1/2s )ζ(s)=3s + 1 + (3/5)s + (3/7)s + (3/9)s + (3/11)s + ...
=( 1 + (1/3)s + (1/5)s + ... ) + 3s( 1 + (1/5)s + (1/7)s + (1/11)s + ... )

ここでは、分母が 3 の倍数のべき乗である項を 3s で約分して前半に集め、それ以外は 3s で括って後半に集めています。すると今度は、前半部が ( 1 - 1/2s )ζ(s) と等しくなるので、上式は次のように変形することができます。

( 1 - 1/2s )( 1 - 1/3s )ζ(s) = 1 + (1/5)s + (1/7)s + (1/11)s + (1/13)s + (1/17)s + ...

次は 5、その次は 7 というように、先頭に表れた分母 ( これは常に素数になります ) を使って同じような処理を繰り返すと、左辺には p を素数として ( 1 - 1/ps )が次々と表れ、右辺からは逆に、分母が p の倍数のべき乗である項が次々と消去されていきます。この操作を無限に繰り返したとき、右辺は 1 に収束していくため、次の等式を導き出すことができます。

ζ(s) = Σk{1→∞}( 1 / ks ) = Π{p∈素数}( 1 - 1/ps )-1

これは「オイラーの積の公式」と呼ばれています。s = 1 のとき左辺は調和級数となり、発散します。任意の素数に対して ( 1 - 1/ps )-1 > 1 なので、その積となる右辺は素数が無限にある場合、同様に発散します。つまり、この等式から、素数が無限にあることを示すことができます。


2) フェルマーテスト

フェルマーの小定理をもう一度以下に示します。

フェルマーの小定理

p を任意の素数、a を p と互いに素な整数としたとき

ap-1 ≡ 1 ( mod p )

が成り立つ

残念ながら、この命題の逆は成り立ちません。ある数 p に対し、p > a かつ p と互いに素である任意の整数 a に対して、ap-1 ≡ 1 ( mod p ) が常に成り立つとしても、p が素数であると断定することはできません。しかし、フェルマーの小定理の対偶として、一つでも ap-1 ≡ 1 ( mod p ) が成り立たない数 a が見つかれば、p が合成数であると断定することはできます。よって、これを利用して素数判定を行うことが可能です。この判定法は「フェルマーテスト ( Fermat Primality Test )」と呼ばれます。

以下に、フェルマーテストのサンプル・プログラムを示します。

/*
  ModularPower : 繰り返し自乗法を使った法 n のべき乗計算( a の k 乗を n で割った余りを求める )
*/
template< class T >
T ModularPower( const T& a, T k, const T& n )
{
  if ( a == 0 || n == 0 ) return( 0 ); // 底や法が 0 の場合は 0 を返す
  if ( k == 0 ) return( 1 % n );       // 指数が 0 の場合は法 n における 1

  T mod( a % n ); // n を法とした a
  T ans( ( ( k & 1 ) > 0 ) ? mod : 1 ); // 求める値

  for ( k >>= 1 ; k > 0 ; k >>= 1 ) {
    mod = ( mod * mod ) % n;
    if ( ( k & 1 ) > 0 )
      ans = ( ans * mod ) % n;
  }

  return( ans );
}

/*
  CalcGCD : ユークリッドの互除法を使った a, b の最大公約数の計算
*/
template< class T >
T CalcGCD( T a, T b )
{
  if ( a < b ) std::swap( a, b );
  if ( b == 0 ) return( a );

  do {
    T r = a % b;
    a = b;
    b = r;
  } while ( b != 0 );

  return( a );
}

/*
  FermatTest : フェルマーテスト

  p : 判定対象の数

  戻り値 : 素数なら true を返す
*/
bool FermatTest( const BigNum::Unsigned& p )
{
  if ( p == 1 ) return( false );
  if ( p == 2 ) return( true );
  if ( ( p & 1 ) == 0 ) return( false );

  BigNum::Unsigned a( 2 ); // a = 2
  BigNum::Unsigned gap( p ); // テスト回数は最大100回程度
  gap /= 100;
  if ( gap == 0 ) gap += 1;
  BigNum::Unsigned pm( p - 1 ); // p - 1

  while ( a < p ) {
    BigNum::Unsigned g( CalcGCD( a, p ) );
    if ( g != 1 ) return( false );

    BigNum::Unsigned pow = ModularPower( a, pm, p );
    if ( pow != 1 ) return( false );
    a += gap;
  }

  return( true );
}

フェルマーテストでは最初に、任意の数 a と判定対象 p が互いに素であるかをチェックしています。チェックには、ユークリッドの互助法が利用されています ( CalcGCD )。a は 1 < a < p の範囲に限定しているため、ここで互いに素でないことが判明すれば、p が 1 と p 以外に約数を持つことになり、合成数であると判断することができます。
互いに素であることが分かったら、ap-1 ( mod p ) を求めます。この計算には繰り返し自乗法が利用されています ( ModularPower )。結果が 1 でなければ、p は合成数であると判断できます。そうでなければ他の数を a として処理を続けます。これを何度か繰り返し、合成数であると判断されなかった場合は素数である可能性が「高い」として true を返します。

全ての a に対して常に ap-1 ≡ 1 ( mod p )が成り立つ合成数は、「R.D.カーマイケル ( Robert Daniel Carmichael )」によって 1910 年に指摘されたため、「カーマイケル数 (Carmichael Numbers )」と呼ばれます。合成数 n がカーマイケル数である必要十分条件は以下のようになります。

  1. n は奇数である。
  2. n を割る素数 p に対し、n は p2 では割り切れない。
  3. n - 1 は p - 1 で割り切れる。

これはカーマイケル数に対する「コルセルトの判定法 ( Korselt's Criterion )」と呼ばれています。最小のカーマイケル数は 561 = 3 x 11 x 17 であり、各素因数は相異なるため、コルセルトの判定法の 2 番めを満たしています。また、各素因数から 1 を引いた値で 561 - 1 = 560 を割ってみると、

560 / ( 3 - 1 ) = 280
560 / ( 11 - 1 ) = 56
560 / ( 17 - 1 ) = 35

となるため、3 番めを満たしていることにもなります。コルセルトの判定法を使えば、フェルマーテストで素数と判定されてしまったカーマイケル数を見つけ出すことが可能です。しかし、そのためには素因数に分解する必要があるため、巨大な数に対しては結局利用することができないということになります。


3) Solovay-Strassen 素数判定法

まず、二次合同式 x2 ( mod p ) の解としてどのような値があるかを検討してみます。例として、x2 ( mod 5 ) を求めてみた場合、x は 0 から 4 までの値を考えればよく、

02 ≡ 0 (mod 5)
12 ≡ 1 (mod 5)
22 ≡ 4 (mod 5)
32 = 9 ≡ 4 (mod 5)
42 = 16 ≡ 1 (mod 5)

なので、解は 0, 1, 4 になります。同様に、p を素数としていくつかの値について解いてみると、結果は次のようになります。

表 3-1. 二次合同式 x2 ( mod p ) の解
p
 5 71113
x00000
11111
24444
34299
41253
54312
61310
7510
8912
943
1019
114
121

ここで、各列の上側と下側は、順序を逆にしただけで全く同じパターンを形成していることがわかります。これは、( p - x )2 = p2 - 2px + x2 ≡ x2 ( mod p ) が成り立つことから容易に理解できると思います。また、上半分だけ着目したときに、各値は相異なっていることが予想でき、これは他の素数を p にしたときも実際に成り立ちます。従って、1 から p - 1 の間に、解は必ず ( p - 1 ) / 2 個存在することになります。

0 でない数に対し、p を法として平方数に合同な数を「p を法とした平方剰余」といいます。逆に、合同でない数は「p を法とした平方非剰余」と呼びます。例えば、1, 4 は 5 を法とした平方剰余、2, 3 は 5 を法とした平方非剰余になります。

平方剰余と平方非剰余を表す便利な記号として「ルジャンドル記号 ( Legendre Symbol )」があります。p を法とした a のルジャンドル記号は (a/p) のように表されますが、ここでは "( a / p )" のように表すことにします。( a / p ) は、a が p を法として平方剰余である場合に 1、平方非剰余である場合に -1 と定義されます。よって、( 1 / 5 ) = ( 4 / 5 ) = 1, ( 2 / 5 ) = ( 3 / 5 ) = -1 になります。

ルジャンドル記号を使った定理に「オイラーの基準 ( Euler's Criterion )」というものがあります。

オイラーの基準

p を奇素数、a を p と互いに素な任意の数としたとき、次の合同式が成り立つ。

a(p-1)/2 ≡ ( a / p ) ( mod p )

p が奇素数であれば、オイラーの基準は任意の a に対して成り立ちます。しかし、p が合成数であるときは合同式が成り立たない a が必ず存在するため、これを利用して素数判定を行なうことが可能になります。

オイラーの基準の左辺は、繰り返し自乗法によって計算可能であるのに対し、右辺のルジャンドル記号は、1 から ( p - 1 ) / 2 までの値を平方して剰余を計算し、a と一致するものがないかを調べる必要があります。しかし、平方剰余に関する有名な法則を利用することで、ルジャンドル記号は簡単に計算することができます。それは、「平方剰余の相互法則」と呼ばれています。

平方剰余の相互法則

1) p を奇素数とする。-1 は p を法として、p ≡ 1 ( mod 4 ) のとき平方剰余、p ≡ 3 ( mod 4 ) のとき平方非剰余である。

( -1 / p ) =  1,  p ≡ 1 ( mod 4 )
           = -1,  p ≡ 3 ( mod 4 )

2) p を奇素数とする。2 は p を法として、p ≡ 1 または 7 ( mod 8 ) のとき平方剰余、p ≡ 3 または 5 ( mod 8 ) のとき平方非剰余である。

( 2 / p ) =  1,  p ≡ 1 or 7 ( mod 8 )
          = -1,  p ≡ 3 or 5 ( mod 8 )

3) p, q を相異なる奇素数とする。( q / p ) は、p ≡ 1 ( mod 4 ) または q ≡ 1 ( mod 4 ) のとき ( p / q ) に、p ≡ q ≡ 3 ( mod 4 ) のとき -( p / q ) に等しい。

( q / p ) =  ( p / q ),  p ≡ 1 ( mod 4 ) or q ≡ 1 ( mod 4 )
          = -( p / q ),  p ≡ q ≡ 3 ( mod 4 )

また、a, b 二つの値の積は、二つの値が、どちらも平方剰余または平方非剰余である場合は「平方剰余」に、片側が平方剰余でもう一方が平方非剰余である場合は「平方非剰余」になります。これを ( a / p ) を使って表すと、

( a / p )( b / p ) = ( ab / p )

が成り立つことになります。これを「平方剰余の積法則」といいます。これらの法則を利用すると、例えば次のようなルジャンドル記号も「簡単に」計算することができます。

(192329/450691) = ((89 x 2161)/450691) = (89/450691)(2161/450691)
                = (450691/89)(450691/2161) = (84/89)(1203/2161) = (2/89)(2/89)(3/89)(7/89)(3/2161)(401/2161)
                = (89/3)(89/7)(2161/3)(2161/401) = (2/3)(5/7)(1/3)(156/401) = (-1)(5/7)(2/401)(2/401)(3/401)(13/401)
                = (-1)(7/5)(401/3)(401/13) = (-1)(2/5)(2/3)(11/13) = (-1)(-1)(-1)(11/13)
                = (-1)(13/11) = (-1)(2/11)
                = (-1)(-1) = 1

a と p が同じものどうしの積は、そのルジャンドル記号が正か負かに関わらず必ず 1 になります。また、常に ( 1 / p ) = 1 が成り立ちます ( 1 より大きな任意の数 p に対して 12 ≡ 1 ( mod p ) なので )。a を「素因数分解」して積の形に分解してから分子と分母を反転し、分子を分母で割った余りを計算することを繰り替えしていくと、やがて分子は 2 以下になります。そうなったら、平方剰余の相互法則を利用して正負の判定を行なえばいいわけです。

先程「簡単に」計算することができると書きましたが、上記処理には「素因数分解」が必要となります。もし、a や p が巨大な桁数を持っていたら、この方法は利用できなくなります。しかし実は、平方剰余の相互法則は、a や p が素数でなくても成り立ちます ( 但し、2) の場合を除き、分子の値は奇数である必要があります )。この一般化されたルジャンドル記号は「ヤコビ記号 ( Jacobi Symbol )」と呼ばれています。上記の例を、素因数分解を行なわずに処理してみましょう。

(192329/450691) = (450691/192329) = (66033/192329)
                = (192329/66033) = (60263/66033)
                = (66033/60263) = (5770/60263) = (2/60263)(2885/60263)
                = (60263/2885) = (2563/2885)
                = (2885/2563) = (322/2563) = (2/2563)(161/2563)
                = (-1)(2563/161) = (-1)(148/161) = (-1)(2/161)(2/161)(37/161)
                = (-1)(161/37) = (-1)(13/37)
                = (-1)(37/13) = (-1)(11/13)
                = (-1)(13/11) = (-1)(2/11)
                = (-1)(-1) = 1

ヤコビ記号を求めるサンプル・プログラムを以下に示します。

/*
  Jacobi : ヤコビ記号を求める

  a : 判定対象の数
  p : 法

  戻り値 : 1 ... 平方剰余 ; -1 ... 平方非剰余 ; 0 ... GCD(a,p) != 1
*/
template< class T >
int Jacobi( T a, T p )
{
  if ( ( p == 1 ) || ( ( p & 1 ) == 0 ) ) return( 0 ); // p が 3 以上の奇数でなければ 0 を返す
  if ( CalcGCD( a, p ) != 1 ) return( 0 );

  // (1/p) = 1
  if ( a == 1 ) return( 1 );

  // (2/p) =  1, ( p ≡ 1 or 7 (mod 8) )
  //       = -1, ( p ≡ 3 or 5 (mod 8) )
  if ( a == 2 ) {
    T modP = p % 8;
    if ( modP == 1 || modP == 7 )
      return( 1 );
    else
      return( -1 );
  }

  // a から因数 2 を除外する
  unsigned int evenCnt = 0;
  while ( ( a & 1 ) == 0 ) {
    a >>= 1;
    ++evenCnt;
  }
  int ans = ( ( evenCnt & 1 ) != 0 ) ? Jacobi( T( 2 ), p ) : 1;

  // (a/p) =  (p/a), ( p ≡ 1 or a ≡ 1 (mod 4) )
  //       = -(p/a), ( p ≡ a ≡ 3 (mod 4) )
  if ( p % 4 == 3 && a % 4 == 3 )
    ans *= -1;
  if ( a == 1 ) return( ans );

  p %= a;

  return( ans * Jacobi( p, a ) );
}

ヤコビ記号を求めることができるようになれば、Solovay-Strassen 素数判定法の処理は簡単に実装できます。

/*
  SolovayStrassenTest : Solovay-Strassenテスト

  p : 判定対象の数

  戻り値 : 素数なら true を返す
*/
template< class T >
bool SolovayStrassenTest( const T& p )
{
  if ( p == 1 ) return( false );
  if ( p == 2 ) return( true );
  if ( ( p & 1 ) == 0 ) return( false );

  T a( 2 );   // a = 2
  T gap( p ); // テスト回数は最大100回程度
  gap /= 100;
  if ( gap == 0 ) gap += 1;
  T half( p - 1 ); // ( p - 1 ) / 2
  half >>= 1;

  while ( a < p ) {
    T pow = ModularPower( a, half, p ); // a^( ( p - 1 ) / 2 ) (mod p)
    int jacobi = Jacobi( a, p );        // (a/p)
    if ( jacobi == 0 ) return( false );
    if ( ! ( ( pow == 1 && jacobi == 1 ) || ( pow == p - 1 && jacobi == -1 ) ) )
      return( false );
    a += gap;
  }

  return( true );
}

Solovay-Strassen 素数判定法では、フェルマーテストでのカーマイケル数のような、素数のフリをする合成数は存在しません。判定する数を n としたとき、n が合成数であれば、1 から n - 1 のうち少なくとも半分は合成数であることを示す数になります。よって、k 回のテストに対し、合成数に対して誤った判定を返す確率は ( 1 / 2 )k であり、100 回のテストでは 8E-31 の確率となります。

Solovay-Strassen素数判定法は、アメリカの集合論学者である「ロバート・ソロベイ ( Robert Martin Solovay )」と、ドイツの数学者「フォルカー・ストラッセン ( Volker Strassen )」によって発表された素数判定法です。今では後述するミラー・ラビン素数判定法に取って代わったものの、RSA 暗号が実現可能であることを示したことで、歴史的には非常に重要な意味を持ったアルゴリズムです。


4) ミラー・ラビン素数判定法

Solovay-Strassen 素数判定法のように、カーマイケル数のような素数のフリをする数が存在しない判定法として、他に「ミラー・ラビン素数判定法 ( Miller-Rabin Primality Test )」があります。ミラー・ラビン素数判定法では、素数の持つ次の性質が利用されています。

pを奇素数とし、

p - 1 = 2kq ( q は奇数 )

と表す。また、a を p で割り切れない任意の数とする。このとき、次の性質が成り立つ。
  1. aq ≡ 1 ( mod p )
  2. K = 2i ( 但し i は 0 から k - 1 までの整数 ) としたとき、aKq ≡ -1 ( mod p ) が成り立つ K が必ず一つ存在する。

この性質を利用して、次のように処理を行います。

  1. 判定対象の数を p、奇数を q として p - 1 = 2kq の形に分解して q と k を求める。
  2. 2 から p - 1 の間の任意の数 a を決める。
  3. aq ≡ 1 ( mod p ) が成り立てば素数と判定する。
  4. i = 0, 1, 2, ... k-1 の中で、K = 2i として aKq ≡ -1 ( mod p ) が一つでも成り立てば素数と判定する。
  5. 3 と 4 をどちらも満たさなかった場合、合成数と判定して終了する。素数と判定されたら 2 に戻って処理を繰り返す。

最小のカーマイケル数 561 をミラー・ラビン素数判定法でチェックしてみると、次のようになります。

561 - 1 = 560 = 24 x 35 より、a = 2 としてチェックすると、

235 ≡ 263 ( mod 561 )
22x35 ≡ 2632 ≡ 166 ( mod 561 )
24x35 ≡ 1662 ≡ 67 ( mod 561 )
28x35 ≡ 672 ≡ 1 ( mod 561 )

aq ( mod p ) = 235 ( mod 561 )は 1 でも -1 でもなく、a2q ( mod p ), a4q ( mod p ), a8q ( mod p )も -1 ではないため、a = 2によって、561 が合成数であることが示されました。

ミラー・ラビン素数判定法のサンプル・プログラムを以下に示します。

/*
  MillerRabinTest : ミラー・ラビンテスト

  p : 判定対象の数

  戻り値 : 素数なら true を返す
*/
template< class T >
bool MillerRabinTest( const T& p )
{
  if ( p == 1 ) return( false );
  if ( p == 2 ) return( true );
  if ( ( p & 1 ) == 0 ) return( false );

  T q( p - 1 );
  T k;
  while ( ( q & 1 ) == 0 ) {
    q >>= 1;
    k += 1;
  }

  T a( 2 ); // a = 2
  T gap( p ); // テスト回数は最大100回程度
  gap /= 100;
  if ( gap == 0 ) gap += 1;
  T pm( p - 1 ); // p - 1

  while ( a < p ) {
    T pow( ModularPower( a, q, p ) ); // a^q ( mod p )
    a += gap;
    if ( pow == 1 ) continue;  // a^q ≡ 1 ( mod p )なら素数と判定
    if ( pow == pm ) continue; // a^q ≡ -1 ( mod p )なら素数と判定
    T i( 1 );
    while ( i < k ) {
      pow = ModularPower( pow, T( 2 ), p );
      if ( pow == pm ) break; // a^(2^i)q ≡ -1 ( mod p )なら素数と判定
      i += 1;
    }
    if ( i == k ) return( false );
  }

  return( true );
}

ミラー・ラビン素数判定法は最初、「ゲイリー・ミラー ( Gary L.Miller )」によって決定的素数判定法として発表されました。ミラーは、リーマン予想が真とした場合、最大でも K( logep )2 ( K は定数 ; これを、ランダウの記号を用いて O( ( logep )2 ) と表します ) までの間に p が合成数であることを示す数 a が必ず存在することを示しました。1990 年に、「Eric Bach」によって定数 K は 2 までに小さくなっています。よって、2 から 2( logep )2 以下の整数 ( 但し、この値が p - 1 より大きい場合は p - 1 ) までの範囲の a に対して判定処理をおこなうことで、p が素数か合成数かを確実に判定することができます。

「マイケル・ラビン ( Michael Oser Rabin )」は、合成数に対してミラー・ラビン素数判定法を行ったとき、a の中の少なくとも 75% は合成数であることを示すことのできる数になることを示しました。k 回のテストを実施した場合、合成数に対して誤った結果を返す確率は ( 1 / 4 )k であり、例えば 100 回のテストでは 6E-61 という非常に小さな値になります。これは、Solovay-Strassen素数判定法よりも誤った判定を行なう確率が低いことになります。
決定的素数判定法の場合に比べてパフォーマンスは圧倒的によく、ミラーの証明が、未だ証明されていないリーマン予想を利用しているため、実際には確率的素数判定法の方が主に利用されています。


5) AKS 素数判定法

入力されるデータのサイズに対して、最大処理時間を入力サイズの多項式で表現できるものを「多項式時間のアルゴリズム ( Polynomial Time Algorithm )」と言います。例えばバブルソートは、要素の比較・交換処理回数が入力サイズ n に対して 2n( n - 1 ) に比例するため、O 記法を用いると O( n2 ) と表され、多項式時間のアルゴリズムになります。
また、入力されたデータに対し、常に同じ経路で計算を行い、常に同じ結果を返すようなアルゴリズムを「決定的アルゴリズム ( Deterministic Algorithm )」といいます。コンピュータを使って実現されたアルゴリズムは、決定的アルゴリズムになります。
判定問題が多項式時間の決定的アルゴリズムであったとき、その問題は「クラス P」に属するといいます。「AKS 素数判定法」は、2002 年 8 月 6 日「PRIMES is in P」と題された論文で発表され、論文を発表したインド工科大学の「Manindra Agarwal」教授と、2 人の学生「Neeraj Kayal」「Nitin Saxena」の頭文字を取って名付けられました。論文のタイトルにある P は「クラス P」のことを指しており、AKS 素数判定法が、クラス P に属する世界初の素数判定アルゴリズムであることを示しています。

AKS 素数判定法においても、フェルマーの小定理が利用されています。まず、フェルマーの小定理の対偶を考えます。

フェルマーの小定理の対偶

a, n を互いに素な整数としたとき

an ≡ a ( mod n )

成り立たなければ n は合成数である。

このままでは、合成数が素数と判定されてしまう場合があるため、次のように改良します。

a, n を互いに素な整数としたとき、n が素数である必要十分条件は

( X - a )n ≡ Xn - a ( mod n )

である。

ここで両辺は、X を変数とする n 次の多項式になります。つまり、次数が等しい両辺の項に対して n を法として合同な数を求めたとき、1 次から n - 1 次の項の係数までが 0 に等しくなれば、n は素数であると判断できます。

しかし、これをまともに処理した場合、最大で n 回の判定処理を行わなければならないため、合同式をさらに変形します。

( X - a )n ≡ Xn - a ( mod Xr - 1, n )

ここで r は、次の条件を満たす素数である必要があります。

q = P( r - 1 ) のとき ( ここで P(x)は x の最大素因数を表す )、r = kq + 1 として
  1. q ≥ 4√r log2 n
  2. nk ≡ 1 ( mod r ) が成り立たないこと

( mod Xr - 1, n ) とは、多項式を Xr - 1 で割ったときの剰余式について、各係数に対し n を法とする合同な数を評価することを意味します。これで、評価すべき多項式の次数は最大でも r - 1 次となり、r が十分に小さければより短い時間で処理を行うことができるようになります。しかし、n が合成数であったとき、上記方法では誤って素数と判定してしまう可能性があります。そこで、a として 1 から 2√r log2 n までの値を使って判定することで誤りを回避することができます。このことを証明した論文が「PRIMES is in P」です。

上記合同式を利用して素数判定を行った場合、その結果は次のようになります。

  1. n が素数の時、全ての合同式が成り立つ
  2. n が素数のべき乗 ( n' )k ( k > 1 ) である場合、合同式が成り立つ場合と成り立たない場合がある
  3. n が r の倍数である場合、全ての合同式が成り立つ
  4. n がその他の合成数である場合、少なくとも一つの合同式は成り立たない

2 と 3 の場合は、合同式を利用する前にチェックすることになります。

AKS 素数判定法を実現する上で最も難しそうなのが、多項式による剰余計算の部分です。しかし、除数が Xr - 1 と特殊な形をしているため、実際には簡単に実装することができます。

まず、多項式 ( X - a )n は、二項定理を使って次のように展開することができます。

( X - a )n = Σk{0→n}( C(n,k) Xn-kak )

但し、C(n,k) は二項係数で、C(n,k) = n! / k!(n-k)!

ここで、C(n,k)ak = ck とすると、

( X - a )n = Σk{0→n}( ckXn-k ) = c0Xn + c1Xn-1 + c2Xn-2 + ... + cn

になります。多項式の場合も、普通の除算と同じように筆算を行うやり方で剰余を求めることができるので、まずは商の最上位を c0Xn-r として、c0Xn-r( Xr - 1 ) を上式から減算します。すると、n 次の項が消えて、次のようになります。

c1Xn-1 + c2Xn-2 + ... + (cr + c0)Xn-r + ... + cn

同じように、商の次の項を c1Xn-r-1 として、c1Xn-r-1( Xr - 1 ) を上式から減算します。

c2Xn-2 + ... + (cr + c0)Xn-r + (cr+1 + c1)Xn-r-1 + ... + cn

これを、n - r + 1 次の項を消去するところまで繰り返すと、次のようになります。

( cr + c0 )Xn-r + ( cr+1 + c1 )Xn-r-1 + ... + ( c2r-1 + cr-1 )Xn-2r+1 + c2rXn-2r + ... + cn

n - r > r - 1 であれば、まだ除算ができることになるので、さらに ( cr + c0 )Xn-2r( Xr - 1 ) を上式から減算して

( cr+1 + c1 )Xn-r-1 + ... + ( c2r-1 + cr-1 )Xn-2r+1 + ( c2r + cr + c0 )Xn-2r + ... + cn

この操作を、最大次数が r - 1 以下になるまで繰り返します。最終的に、K を 0 以上の整数として n - Kr 次の項が r - 1 次になり、その係数は cKr + c(K-1)r + ... + c0 になります。また、それより次数の一つ小さい項の係数は、cKr+1 + c(K-1)r+1 + ... + c1 になり、以下、ci の各添字番号が一つずつ増えていきます。
少し分かりづらいと思いますので、( X - 1 )5 を X4 - 1 で割った余りを算出してみましょう。

( X - 1 )5 = c0X5 + c1X4 + c2X3 + c3X2 + c4X + c5

但し、c0 = 1, c1 = -5, c2 = 10, c3 = -10, c4 = 5, c5 = -1

c0X( X4 - 1 ) を引くと、c1X4 + c2X3 + c3X2 + ( c4 + c0 )X + c5

c1( X4 - 1 ) を引くと、c2X3 + c3X2 + ( c4 + c0 )X + ( c5 + c1 )

ci に数値を代入して、10X3 - 10X2 + 6X - 6

上記計算を見ると、c0 を持つ項の次数は n を r で割ったときの剰余になります。よって、係数を求めるためには、n を r で割ったときの剰余を次数とする項を開始点として、最も次数の高い項の係数から順番に加算する処理を実行すればよいことになります。0 次まで処理したら r - 1 次へ移動して、同じ処理を続けます。何回ループしても、最後は必ず 0 次の項が処理されることになります。

多項式どうしの除算処理による剰余を求めるサンプル・プログラムを以下に示します。

/*
  Unsigned2Buff : BigNum::Unsigned 型から BigNum::UnsignedBuff 型への変換

  下二桁分を抽出して BigNum::UnsignedBuff 型へ変換する

  n : 対象の多倍長整数
*/
BigNum::UnsignedBuff Unsigned2Buff( BigNum::Unsigned n )
{
  if ( n.size() < 2 )
    return( n[0] );
  else
    return( ( n[1] << ( sizeof( BigNum::Digit ) * CHAR_BIT ) ) + n[0] );
}

/*
  PolyMod : 多項式 ( x - a )^n を x^r - 1 で割ったときの余りの係数を求める

  ( x - a )^n [ = ( c[0]x^n + c[1]x^n-1 + ... + c[n]x^0 ) ] mod ( x^r - 1 )
                   但し c[i] = C(n,i)a^i

  n : 多項式の次数
  a : 多項式の第二項
  r : 除数となる多項式の次数
  res : 余りの係数へのポインタ
*/
void PolyMod( const BigNum::Unsigned& n, const BigNum::Unsigned& a, BigNum::UnsignedBuff r, vector< BigNum::Signed >* res )
{
  if ( n < r ) return;

  res->assign( r, 0 );
  BigNum::Unsigned coef( 1 ); // C(n,i)の初期値 = C(n,0) = 1
  BigNum::Unsigned coefA( n );  // C(n,i)に掛ける値 = coefA/coefB (初期値 = n/1)
  BigNum::Unsigned coefB( 1 );
  BigNum::Signed powA( 1 );  // a^iの初期値 = a^0 = 1
  BigNum::UnsignedBuff i = Unsigned2Buff( n % r ); // 加算処理を行う開始点となる次数

  for ( BigNum::Unsigned num ; num <= n ; num += 1 ) {
    //(*res)[i] += BigNum::Signed( coef ) * powA % BigNum::Signed( n );
    (*res)[i] += BigNum::Signed( coef ) * powA;
    ( coef *= coefA ) /= coefB;
    //( ( powA *= BigNum::Signed( a ) ) %= BigNum::Signed( n ) ).revSign();
    ( powA *= BigNum::Signed( a ) ).revSign();
    coefA -= 1;
    coefB += 1;
    i = ( i > 0 ) ? i - 1 : r - 1;
  }
}

/*
  PolyModTest : PolyMod用テストルーチン

  n : 多項式の次数
  a : 多項式の第二項
  r : 除数となる多項式の次数
*/
void PolyModTest( const BigNum::Unsigned& n, const BigNum::Unsigned& a, BigNum::UnsignedBuff r )
{
  vector< BigNum::Signed > res;
  cout << "( x - " << a << " )^" << n << "/ ( x^" << r << " - 1 ) = ";
  PolyMod( n, a, r, &res );
  for ( BigNum::UnsignedBuff i = res.size() ; i > 0 ; --i )
    cout << " " << ( ( res[i-1] < 0 ) ? "" : "+" ) << res[i-1] << "x^" << i - 1;
  cout << endl;
}

サンプル・プログラムの中で、PolyMod が多項式 ( x - a )n を xr - 1 で割ったときの剰余を計算する関数になります。n を r で割った剰余を計算し、それを添字として扱う必要があるので、n % r の値を BigNum::UnsignedBuff 型に変換するため関数 Unsigned2Buff を用意しています。BigNum::Unsigned のメンバ関数 operator[] は多倍長整数を配列とみなして指定した添字の桁の数値を取得するためのもので、これを使って最下位から二桁分を抽出します。

二項係数を表す変数 coef は、初期値を 1 としています。これは C( n, 0 ) = n! / 0!( n - 0 )! = 1 が任意の整数 n について成り立つためです。また、

C( n, k + 1 )=n! / ( k + 1 )![ n - ( k + 1 ) ]!
=n( n - 1 )( n - 2 )...[ n - ( k - 1 ) ]( n - k ) / ( k + 1 )!
=n( n - 1 )( n - 2 )...[ n - ( k - 1 ) ]( n - k ) / ( k + 1 )k!
=[ ( n - k ) / ( k + 1 ) ]{ n( n - 1 )( n - 2 )...[ n - ( k - 1 ) ] / k! }
=[ ( n - k ) / ( k + 1 ) ][ n! / k!( n - k )! ]
=[ ( n - k ) / ( k + 1 ) ]C( n, k )

より、C( n, 1 ) = nC( n, 0 ) = n, C( n, 2 ) = [ ( n - 1 ) / 2 ]C( n, 1 ) = n( n - 1 ) / 2 ... のように、前に求めた値から簡単に計算することができます。

上記計算を行った場合、例えば n = 100, a = 4, r = 3 の時、結果は次のようになります。

( x - 4 )^100/ ( x^3 - 1 ) =
 -303696310102118487255424664408553239076929401524563149137089022624x^2
 -545636784969068263337519668712869806067000247202107914064841861375x^1
 +849333095071186751108321853853434376180390778492292335904038406000x^0

パラメータが大きくなると係数の桁数も多くなるため、これがなるべく小さくなるように合同式を活用することを考慮すると、a ≡ b ( mod n ) が成り立つとき、任意の整数 K に対して

a + K ≡ b + K ( mod n )

a - K ≡ b - K ( mod n )

aK ≡ bK ( mod n )

が成り立つため、最終的に剰余を評価すればいいのであれば、サンプル・プログラムの中で (*res)[i] に係数を加算するのではなく、n で除算した剰余を加算しても問題ないことになります。同様に、powA ( = ai ) を掛ける代わりに、n による剰余を掛けても合同式は成り立ちます ( コメントアウトされた部分に切り替えることで実現可能です )。
しかし、二項係数には除算が含まれるため、剰余を計算すると式が成り立たなくなります。これは、

C( n, 1 ) = n ≡ 0 ( mod n )

に対し、C( n, 2 ) = n( n - 1 ) / 2 は n が偶数の場合 n で割り切れるとは限らないことからも容易に理解できると思います。

ここまでできたら、あとは r を求める処理とメインルーチンを作成するだけです。

/*
  PowerTest : 整数のべき乗で表せる数かどうかをチェックする

  n : 対象の自然数
*/
template< class T >
bool PowerTest( const T& n )
{
  if ( n < 4 ) return( false );

  T m( 2 );
  T k( 0 );

  // n < 2^k となる k を求める
  while ( m < n ) {
    m <<= 1;
    if ( m == n ) return( true );
    k += 1;
  }

  // ( n^(1/i) )^i = n となる i が存在するか
  for ( T i( k ) ; i > 0 ; i -= 1 )
    if ( PowerRoot( n, i ).pow( i ) == n ) return( true );

  return( false );
}

/*
  GetMaxPrimeFactor : num の最大の素因数を求める

  primeList : 素数表
*/
size_t GetMaxPrimeFactor( size_t num, const vector< size_t >& primeList )
{
  for ( vector< size_t >::const_reverse_iterator cri = primeList.rbegin() ;
        cri != primeList.rend() ; ++cri ) {
    if ( ( num % *cri ) == 0 )
      return( *cri );
  }

  return( 1 );
}

/*
  log2n : log2n以上の最小値を求める
*/
size_t log2n( size_t n )
{
  size_t l = 0;
  for ( size_t sqP = 1 ; sqP < n ; sqP <<= 1 )
    l += 1;

  return( l );
}

/*
  AKS_r : AKS素数判定用のパラメータ rを求める

  p : 素数判定を行う数

  戻り値 : r ( p が合成数と判明した場合は 0 を返す)
*/
size_t AKS_r( size_t p )
{
  size_t r;
  vector< size_t > primeList;
  Eratosthenes( p - 1, &primeList );
  for ( vector< size_t >::const_iterator cit = primeList.begin() ;
        cit != primeList.end() ; ++cit ) {
    r = *cit;
    if ( CalcGCD( p, r ) != 1 ) // pは合成数
      return( 0 );
    // r - 1 = qk ( q は r - 1 の最大の素因数 )
    size_t q = GetMaxPrimeFactor( r - 1, primeList );
    size_t k = ( r - 1 ) / q;
    size_t sqR = std::sqrt( r ); // r の平方根以上の最大整数
    if ( sqR * sqR < r ) ++sqR;
    if ( q < 4 * sqR * log2n( p ) )
      continue;
    if ( ModularPower( p, k, r ) != 1 )
      return( r );
  }

  return( p );
}

/*
  AKSTest : AKS素数判定法

  p : 判定対象の数

  戻り値 : 素数なら true を返す
*/
bool AKSTest( size_t p )
{
  // べき乗根ならば合成数
  if ( PowerTest( p ) ) return( false );

  size_t r = AKS_r( p );
  if ( r == 0 ) return( false );
  //if ( p <= r ) return( true );

  vector< BigNum::Signed > coef;
  size_t sqR = std::sqrt( r );
  if ( sqR * sqR < r ) ++sqR;
  for ( size_t a = 1 ; a <= 2 * sqR * log2n( p ) ; a += 1 ) {
    PolyMod( p, a, r, &coef );
    for ( size_t i = 1 ; i < coef.size() - 1 ; ++i ) {
      if ( ( coef[i] % BigNum::Signed( p ) ) != 0 ) {
        cout << "( x - " << a << " )^" << p << " % x^" << r << " - 1 -> ";
        cout << coef[i] << "x^" << p - i << endl;
        cout << coef[i] << " % " << p << " = " << coef[i] % p << endl;
        return( false );
      }
    }
  }

  return( true );
}

除数となる多項式の次数 r は関数 AKS_rで求めています。判定対象の数 p より小さな値の中から先に示した条件を満たすものを探し、最小の値を返します。ループ内において、最初に pr の最大公約数を求め、その値が 1 以外であった場合は p は合成数ということになるので、それを示す戻り値 0 を返しています。その後は前述した内容の通りに以下のような処理を行っています。

  1. r が素数でなければ無視
  2. r - 1 = qk を満たす最大の素因数 q を求める ( GetMaxPrimeFactor )
  3. q < 4√r log2 p ならば無視
  4. nk ≡ 1 ( mod r ) でなければ r を返す

メインルーチンは AKSTest になります。この中でコメントアウトされた部分を復元すると、prになった場合は多項式を求めずに処理を終了することができます。小さな値に対しては後の処理が不要となるため、通常は復元した方が処理は早くなります。上記では、動作確認をする目的でコメントアウトとなっています。

さて、せっかく苦労して作成したサンプル・プログラムですが、残念なことに、今まで紹介した判定方法に比べると処理時間は非常に長くなります。そもそも、r - 1 の最大素因数を求めるために、p - 1 までの素数をエラトステネスのふるいを使って全て求めているので、このサンプル・プログラムに関して言えばわざわざ AKS 素数判定法を利用する意味はなくなります。多項式時間のアルゴリズムが処理の早いアルゴリズムであるというわけではないということで、三桁の値の処理でもしばらくは結果が返ってきません。
しかし、フェルマーテストや Solovay-Strassen 素数判定法、ラビン・ミラー素数判定法は、素数である確率が「非常に高い」ということしかわかりません。このような素数は「確率的素数 ( Probable Prime )」と呼ばれています。素数を求めるときにはこれで充分信頼できますが、例えばこれらの判定法で合成数であることが見抜けないような「疑似素数」が悪用された場合に、それらを見つけることはできなくなります。実用的ではないものの、把握できる範囲内の時間で処理が可能であり、確実に素数判定ができるアルゴリズムという意味では、画期的な発見だったわけです。


素数判定法として他に「楕円曲線素数判定法 ( Elliptic Curve Primality Proof; ECPP )」があります。これは内容を勉強中の段階なので、また別の機会で紹介したいと思います。


<参考文献>
  1. 「はじめての数論」 Joseph H.Silverman著 (ピアソン・エデュケーション)
  2. 「素数に憑かれた人たち リーマン予想への挑戦」 John Derbyshire著 (日経BP社)
  3. 多項式時間素数判定アルゴリズムについて
  4. University of Miami's Department of Computer Science Solovay-Strassen Primality Test (Home page: Burton Rosenberg)
  5. Wikipedia

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