かの有名なフェルマー(1665年)は、つぎのフェルマー数 Fn
    Fn = 22n + 1 (n = 0, 1, 2, ...)

は素数であろうと予想した。この数は、n = 5, 6, 7, 8, ... に応じて非常に速くとてつもなく大きくなる。したがってコンピュータなしの時代には、これらの値を計算することはほとんど不可能であった。確かに、

    F0 = 3  F1 = 5  F2 = 17  F3 = 257  F4 = 65537

は素数である。しかし、オイラー(1732年)により

    F5 = 232+1 = 4294967297 = 641 x 6700417

は合成数であることが後にわかった。

本関数ではペパンの判定法を利用して、フェルマー数が素数であるかどうかを調べる。ペパンの判定法によると、n が 1 以上のとき

    3 (Fn-1) / 2 ≡ -1 (mod Fn)

の成り立つことが Fn が素数であるための必要十分条件。

いまのコンピュータ時代では、多くの人たちの努力にもかかわらず、F4より大きな素数は1つも見つかっていない。なお、フェルマー数が合成数の場合、その素因数はすべて

    k 2m+1 (m≧n+2)

の形をしている。

<関数名>
  pepan フェルマー(Fermat)数 Fn=22n+1 が素数であるかどうかを判定する

<形式>
  int pepan(int n);

<引数>
  n フェルマー数における正整数

<関数値>
  素数なら 1、非素数(つまり、合成数)なら0、内部エラーが起きたときは -1。

用例
  pepan(3);

<関数本体>
  pepan.c

素数となるフェルマー数(フェルマー素数)を見つけることはなぜそんなに重要なのかというと、一つには、全く関係がなさそうな幾何学の正多角形の作図問題と深い関係があることをガウスが発見したからである。

ガウスは、n の値が

    ・2 のベキである
    ・フェルマー素数である
    ・これらの2種類の数の積である

ときに、かつそのときにだけ、正 n 角形は定規とコンパスで作図可能であることを証明した。ガウスの墓石に刻まれた正17角形は、ちょうどフェルマー数 F2に対応するものである。したがって作図可能な奇数正多角形はフェルマー素数より、

    正3角形、5角形、15角形、17角形、51角形、85角形、...

などのみとなる。

Mp = 2p - 1

(p は素数)の形の数をメルセンヌ (Mersenne) 数という。

p が2、3、5、7、13、17、19、31、61、89、107、127、521、607、1279、2203、2281、3217、4253、4423、9689、9941、11213、19937、21701、23209、44497、86243、110503、132049、216091、756839、859433 のとき

Mp

は素数であることがわかっている。ちなみに M859433は10進数で約26万桁の数であり、1994年1月にLucasテストによって発見された巨大素数である。

Lucasテストとは、ルーカス (Lucas) が1876年に発見したもので、メルセンヌ数が素数かどうかを簡単に判定できる。

   

L0 = 4、 Li+i =(Li2 - 2) mod (2p - 1)

で、

Lp-1

が 0 なら

Mp

は素数、0 以外なら非素数という。

Lucas テストのおかげで、コンピュータの力による素数の発見競争は、これからも続けられ、記録は塗り替えられていくであろう。

<関数名>
  lucas ---- メルセンヌ(Mersenne)数が素数であるかどうかを判定する

<形式>
  int lucas(int p);

<引数>
  p メルセンヌ数における素数

<関数値>
  素数なら1、非素数(つまり、合成数)なら0、内部エラーが起きたときは -1。

用例
  lucas(3217);

<関数本体>
  lucas.c

ペル方程式(Pell Equation)は整数解を求める不定方程式の1種で、一般論が完成されている数少ない方程式である。

ペル方程式の名はオイラーが誤って命名したものとされている。ペルは実在の数学者であるが、この方程式に対する研究は全くなかった。

ペル方程式

    x2 - M y2 = c  (M は平方数でない正整数、c = +1, -1)

において、自明な解 x = 1、y = 0 以外の x、y が最も小さい正の整数である解を、そのペル方程式の最小解という。方程式の右辺 c = -1 としたときには解がない場合が多いが、しかし -1 に対する最小解が存在する場合には、必ずその最小解のほうが +1 に対する最小解よりも小さい。

例えば、M = 2 の場合の方程式については、x、y のペア (1,1)、(3,2)、(7,5)、 (17,12)、(41,29)、(99,70)、(239,169)、(577,408)... がその解であるが、 -1 に対する最小解は x=1、y=1 であり、+1 に対する最小解は x=2、y=3 である。

<関数名>
  pellEqu ---- ペル方程式 x2 - M y2 = c の最小解を得る。ただし、M は平方数でない正整数、c = 1, -1。

<形式>
  int pellEqu(long *x1, long *y1, long *x2, long *y2, int m);

<引数>
  x1, y1 (出力)ペル方程式 x2 - M y 2 = 1 に対する最小解
  x2, y2 (出力)ペル方程式 x2 - M y 2 = -1 に対する最小解
  m    (入力) ペル方程式の係数 M (1,4,9,16...以外の非平方数)

<関数値>
  c = 1 または -1 の両方に対する最小解が得られたときは 1、
  c = -1 に対する最小解が存在しなかったときは 0。
  c = 1 に対する最小解も c = 1 に対する最小解も得られなかったときは -1。

<注意事項>
  ペル方程式の解は必ず存在するが、限られた整数範囲(long整数で表現できる範囲内)のなかで解が見つからなかったときには、関数値は -1 となる。

用例

<関数本体>
  pellEqu.c

<説明>
  ペル方程式の一般解はつぎの漸化式で求められる。

      xn+1 = x1xn + M y1yn
      yn+1 = x1yn + y1xn

    ただし x1、y1は最小解を表す。したがってペル方程式を解くにはその最小解を求めればよい。
  M の値によって、最小解が意外に大きな数になることがある。例えば、M=61、c=1 に対する最小解は x=1766319049、y=226153980 にもなる。
  最小解はつぎのように漸化的に求められる。

      s0 = 0  t0 = 1
      a0 = M1/2 の整数部
      x-1 = 1  y-1 = 0  x0 = a0  y0 = 1
      sk+1 = aktk - sk  tk+1 = (M - sk+12) / tk
      ak+1 = (a0 + sk+1) / tk+1の整数部
      xk+1 = ak+1 xk + xk-1  yk+1 = ak+1 yk + yk-1

  tk+1の値が 1 になったときの、xk、ykの値がペル方程式の最小値となる。

奇素数 p と p で割れない a に対して、

    x2≡a (mod p)

なる x は存在するときと、存在しないときがある。

たとえば、p = 7 のとき

    12≡1、 22≡4、32≡2、42≡2、 52≡4、62≡1

であるから、a≡1 または 2 または 4 のときのみ解がある。解があるとき、
a を mod p での平方剰余という。a がある数を平方して p で割った余りと
なるからである。解がないときは、平方非剰余という。このとき、

    (a/p)= 1   解があり
    (a/p)=-1   解がなし

と定義される。また、この (a/p) を平方剰余記号という。

<関数名>
  jacobi ---- 平方剰余記号 (a/p) の値を計算する

<形式>
  int jacobi(long a, long p);

<引数>
  a、p それぞれ平方剰余記号の分子、および分母

<関数値>
  平方剰余記号の値 +1 または -1。引数により計算不能の場合は 0。

<注意事項>
  平方剰余記号の分母 p は奇素数でなければならない。また、分子 a は p で割れないものでなければならない。特に奇素数の判定については手間がかかるので、本関数内では奇素数かどうかの判定はしていない。

用例

<関数本体>
  jacobi.c

<説明>
  平方剰余記号の値を計算するのに、以下の性質が利用される。

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

      (2/p) = 1、p≡1 または 7 (mod 8) 時
      (2/p) = -1、p≡3 または 5 (mod 8) 時

      奇数 a と互いに素な p に対して、ガウス定理より
      (a/p) = (p/a)、p≡1 または a≡1 (mod 4) 時
      (a/p) = -(p/a)、p≡3 かつ a≡3 (mod 4) 時

  上の性質から、a の値が奇数か偶数かで異なる処理をする。偶数の場合は、 (2/p)の値をかけておき、a 自身の値を半分にする。奇数の場合は、分母と分子を入れ替え、ユークリッド互除法を適用する。この繰り返しを a が 1 になるまで行う。
  本アルゴリズムでは平方剰余記号の値を log(p) に比例時間で計算される。

すこし難しくなるかな。

合同の知識は、暗号の理論等にも良く使われてて、コンピュータについてさらに詳しく知りたい人には基本知識の1つ。難しく感じたら、初等整数論の教科書等を見るといいだろう。

一次合同式を解くとは、a、b を自然数とするとき、a x ≡ b (mod m) となる x を求めることをいう。

例えば、3x ≡ 2 (mod 7) の解は、 x ≡ 3 (mod 7) で与えられる。

一次合同式が解を持たないこともある。例えば、 3x ≡ 1 (mod 6) は解を持たない。

<関数名>
  congruence ---- 一次合同式 a x≡b (mod m) を満たす x を求める

<形式>
  int congruence(int a, int b, int m);

<引数>
  a、b、m ともに自然数

<関数値>
  合同式の解 x の値 (0≦x<m)。解がない時に、関数値は -1 となる。

用例

<関数本体>
  congruence.c

<説明>
  合同式 a x≡b (mod m) を解くには、もう一つ自明な合同式 m x≡0 (mod m) を利用して、同値な変形を繰り返せばよい。すなわち、x の係数 a と m に対して ユークリッド互除法 を用いて、 x の係数を小さく変形していく。

これからは、整数論関係のプログラム(関数)をのせることにするね。

最初の関数は、最大公約数を求めてくれるもの。ユークリッド互除法ともいう。

2つ以上の自然数に共通する約数(公約数)のうち、最大のものを最大公約数という。 ちなみに、いくつもの自然数に共通する倍数(公倍数)のうち最小のものを最小公倍数という。

最小公約数は 1、最大公倍数は無限大なので、最小公約数や最大公倍数を考えることはふつうはない。

<関数名>
  gcd ---- 最大公約数を算出する

<形式>
  int gcd(int n1, int n2);

<引数>
  n1、n2 2つの自然数

<関数値>
  2つの自然数の最大公約数

用例
  gcd(1996, 2000);

<関数本体>
  gcd.c

<説明>
  最大公約数を求める方法は、2千年以前の古代ギリシャ人により発見され、ユークリッドの著書「Elements」に詳しく書かれており、ユークリッドのアルゴリズムといわれている。
  Lameの定理により、このアルゴリズムは小さい方の自然数の桁数の5倍以内で必ず終了するが判り、計算量は O(logN) になる。
 また、n個数の最大公約数を求めるには、まず最初の2つについてここの方法を適用して最大公約数を求め、つぎに3つ目の数と計算した最大公約数から再び最大公約数を計算する。このように最後の数まで繰] り返せば、n個の自然数の最大公約数が求められる。

2進数字列を、整数値に変換する関数。

たとえば、2進数字列 "101011" → 数値 43。

<関数名>
  binsz2int ---- 2進数字列を整数値に変換する

<形式>
  int binsz2int(char *str);

<引数>
  str 2進数 ASCIZ文字列。先頭に符号 (+またはー) があってもOK。

<関数値>
  入力の2進数字列に対応する整数値(10進数)

<注意事項>
  整数値のオーバーフローに注意。

用例
  binsz2int("-1101");

<関数本体>
  binsz2int.c

<説明>
  2進数字列の先頭から、まずスペースを読み飛ばす。つぎに符号を処理し、2進数字の各桁を変換していく。

10進数字列を、整数値に変換する関数。

たとえば、10進数字列 "12345" → 数値 12345。

<関数名>
  sz2int ---- 10進数字列を整数値に変換する

<形式>
  int sz2int(char *str);

<引数>
  str 10進数 ASCIZ文字列。先頭に符号 (+またはー) があってもOK。

<関数値>
  入力の10進数字列に対応する整数値(10進数)

<注意事項>
  整数値のオーバーフローに注意。

用例
  sz2int("-6119");

<関数本体>
  sz2int.c

<説明>
  10進数字列の先頭から、まずスペースを読み飛ばす。つぎに符号を処理し、10進数字の各桁を変換していく。

整数値を、16進数字列への変換関数。

たとえば、値 180 が 数字列(16進数字列) "B4" に変換される。

<関数名>
  int2hexsz ---- 整数値を16進数字列へ変換する

<形式>
  void int2hexsz(char *str, int num);

<引数>
  str  (出力) 16進数ASCIZ文字列
  num (入力) 変換したい整数値(マイナスの値でもOK)

<関数値>
  なし

用例
  char str[15];
  int2hexsz(str, 12345);

<関数本体>
  int2hexsz.c

<説明>
  再帰法を使用。先行の数字については自分自身を呼出して対応し、最終桁だけを変換して処理する。その際、A-F への変換に注意すること。

整数値を、2進数文字列に変換する関数。

たとえば、値 13 が 数字列(2進数列) "1101" に変換される。

<関数名>
  int2binsz ---- 整数値を2進数字列へ変換する

<形式>
  void int2binsz(char *str, int num);

<引数>
  str  (出力) 2進数ASCIZ文字列
  num (入力) 変換したい整数値(マイナスでもOK)

<関数値>
  なし

<注意事項>
  マイナスの整数値が、マイナス記号が先頭につく2進数字列に変換される。
  例: -100 → "-1100100"

用例
  char str[50];
  int2binsz(str, 12345);

<関数本体>
  int2binsz.c

<説明>
  再帰法を使用。先行の数字については自分自身を呼出して対応し、最終桁だけを変換して処理する。
  なお、ASCIZ文字列とは、最後にヌル(Null) コードで終了する文字列のこと。