Cornacchia-Smith法のC言語による実装例

奇素数 p と、p で割り切れない正整数 d に対して、方程式 p = x^2 + d * y^2 の整数解 (x, y) を計算するアルゴリズムである Cornacchia-Smith 法 の C 言語による実装例。

プログラム例

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

typedef int fp2[2];  /* F_{p^2} の元を表す */

/* a を m で割った最小非負剰余 */
int mod( int a, int m )
{
  if ( m < 0 ) { m = -m; }
  if ( a < 0 ) { 
    a = m - ((-a) % m);
  }
  else { 
    a = a % m; 
  }

  return a;
}

/* Jacobi Symbol */
int jacobi_symbol( int a, int m )
{
  int t, temp;

  a = mod( a, m );  
  t = 1;
  while ( a != 0 ) {
    while ( a % 2 == 0 ) {
      a = a / 2;
      if ( m % 8 == 3 || m % 8 == 5 ) { t = -t; }
    }
    temp = a; a = m; m = temp;
    if ( a % 4 == 3 && m % 4 == 3 ) { t = -t; }
    a = a % m;
  }

  if ( m == 1 ) { return t; }
  return 0;
}

/* F_{p^2} での積:c = ab */
void fp2_mult( fp2 c, const fp2 a, const fp2 b, int omega2, int p )
{
  fp2 temp;  /* c *= b のときに必要 */

  temp[0] = mod( a[0] * b[0] + a[1] * b[1] * omega2, p );
  temp[1] = mod( a[0] * b[1] + a[1] * b[0], p );
  c[0] = temp[0]; c[1] = temp[1];
}

/* F_{p^2} での累乗:r = a^n */
void fp2_pow( fp2 r, const fp2 a, unsigned n, int omega2, int p )
{
  fp2 b; 
  b[0] = mod( a[0], p ); b[1] = mod( a[1], p );  /* b = a */
  r[0] = 1; r[1] = 0;  /* r = 1 */
  while ( n > 0 ) {
    if ( n % 2 == 1 ) { 
      fp2_mult( r, r, b, omega2, p );  /* r *= b */ 
    }
    fp2_mult( b, b, b, omega2, p );    /* b *= b */
    n /= 2;
  }
}

/* 平方根 mod p */
int fp_sqrt( int a, int p )
{
  int t, omega2;
  fp2 r;

  /* a が平方剰余でないときは 0 を返す */
  if ( jacobi_symbol( a, p ) != 1 ) { return 0; }

  do {
    t = rand() % p;
    omega2 = mod( t * t - a, p ); 
  } while ( omega2 == 0 || jacobi_symbol( omega2, p ) != -1 );

  r[0] = t; r[1] = 1;
  fp2_pow( r, r, (p+1)/2, omega2, p );  /* 必ず F_{p} に入る */

  return r[0];  /* r[1] = 0 */
}

/* sqrt(a) の整数部分 */
unsigned isqrt( unsigned a )
{
  unsigned x, y;

  y = a;
  do {
    x = y;
    y = (x + a / x) / 2;
  } while ( y < x );

  return x;
}

/* 平方数かどうか */
int is_sqrt( unsigned a )
{
  unsigned b = isqrt(a); 
  return b * b == a;
}

/* p = x^2 + d * y^2 の整数解 (x, y) を求める */
void cornacchia_smith( int p, int d, int *x, int *y )
{
  int x0, a, b, c, t, temp;

  if ( jacobi_symbol( -d, p ) == -1 ) {
    *x = *y = 0; 
    return;
  }

  a = p; 
  b = fp_sqrt( -d, p ); 
  if ( 2 * b < p ) { b = p - b; }
  c = (int) isqrt( (unsigned) p );
  while ( c < b ) {
    temp = a % b; a = b; b = temp;
  }

  t = p - b * b;  /* 必ず t >= 0 */
  if ( t % d != 0 || !is_sqrt( t / d ) ) { 
    *x = *y = 0; 
    return; 
  }
  *x = b; *y = isqrt( t / d );
}

int main(void)
{
  int p = 593;  /* 奇素数 */
  int d = 2;    /* p と互いに素な正整数 */
  int x, y;

  cornacchia_smith( p, d, &x, &y );

  printf( "p = %d\n", p );
  printf( "d = %d\n", d );
  printf( "p = x^2 + d * y^2 の整数解: (x, y) = (%d, %d)\n", x, y );

  return 0;
}

実行結果:

p = 593
d = 2
p = x^2 + d * y^2 の整数解: (x, y) = (9, 16)

参考文献

  • Crandall (著)、Pomerance (著)、和田秀男 (監訳): 素数全書 計算からのアプローチ、朝倉書店、2010

関連記事

奇素数を法とする平方根を計算するCプログラムの例
Jacobi記号の計算のC言語による実装例
累乗根の整数部分を求めるNewton反復に関する命題

【theme : プログラミング
【genre : コンピュータ

プロフィール

よしいず

Author:よしいず
MATHEMATICS.PDFというウェブサイトを運営しています。

管理の都合上、トラックバックとコメントはオフにしてあります。ブログ経験者なら分かっていただけると思いますが、スパム(アダルトやその他の宣伝)ばかりなのが現実です。

リンクは自由です。当サイトの記事に対する間違いの指摘・意見・感想などを述べた記事からのリンクは歓迎です。ただし、ブログ記事アップ直後はミスが多く、頻繁に修正します。場合によっては削除する可能性もあります。その際、何も断りもなく修正・削除しますがご了承ください。内容を参考にする場合には投稿後一週間ほど様子を見てからにしてください(笑)。

記事の間違いを指摘するときは、その具体的箇所、理由(仕様に反するなど)・根拠(参考にした文献など)、代替案(同じ結果を得るための正しいやり方)も教えてください。そうしないと、(指摘される側および第三者はその時点では無知の状態なので、)どこが間違いなのか分かりませんし、本当に間違っているのかどうかが判断・検証できません。実際、間違いだと指摘されたことが結局は正しかったというケースもありますので。

このブログのタイトル一覧

リンク
月別アーカイブ
カテゴリ
最新記事
検索フォーム
RSSリンクの表示