中国式剰余定理のC言語によるプログラム例

中国式剰余定理の C 言語による素朴な実装の例。ここでいう「素朴」とは、Gauss の方法としてよく知られているアルゴリズムをそのままプログラムにしたという意味です。この記事では述べませんが、中国式剰余定理の実装においては、さまざまな効率化・高速化の方法があります。

プログラム例

#include <stdio.h>

/* product of an array */
int prod_array( const int a[], int size )
{
  int i, res;

  res = 1;
  for ( i = 0; i < size; i++ ) { 
    res *= a[i]; 
  }

  return res;
}

/* least non-negative residue of a modulo 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;
}

/* a^{-1} mod m */
int inv_mod( int a, int m )
{
  int b, x, u, q, abs_m, tmp; 

  abs_m = ( m < 0 ) ? -m : m;
  b = m; x = 1; u = 0; 
  while ( b > 0 ) {
    q = a / b; 
    tmp = u; u = x - q * u; x = tmp;
    tmp = b; b = a - q * b; a = tmp;
  }

  return ( x < 0 ) ? abs_m + x : x;
}

/* Chinese Remainder Theorem (CRT) */
void chrem( 
  const int a[], const int m[], 
  int size, int *x, int *M )
{
   int M_i, t_i, i;

   /* Gauss's algorithm */
   *M = prod_array( m, size );
   *x = 0;
   for( i = 0; i < size; i++ ) {
     M_i = *M / m[i];
     t_i = inv_mod( M_i, m[i] );
     *x = mod( *x + a[i] * M_i * t_i, *M );
   }
}

/* 
  連立合同式
  x = 1 mod 3
  x = 2 mod 5
  x = 3 mod 7
  の整数解 x mod M を計算する
  → (x, M) = (52, 105)
*/
int main(void)
{
  const int a[] = { 1, 2, 3 }; 
  const int m[] = { 3, 5, 7 };
  const int size = 3;
  int x, M;

  chrem( a, m, size, &x, &M );
  printf( "(x, M) = (%d, %d)\n", x, M );

  return 0;
}

参考文献

  • 高木貞治: 初等整数論講義 第2版, 共立出版, 1971

関連記事

法mに関する逆元を計算するCプログラムの例

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

プロフィール

よしいず

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

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

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

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

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

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