合成数は素数の積で表すことができる。これを素因数分解という。たとえば、720 = 24x32x5。素因数分解は素数判定以上に難しい。

素因数分解のプログラムを話しする前に、まず2点のことを確認しておく。

1.素因数分解の一意性。つまり、素因数分解において、素因数(=素数)の組み合わせは一通りしか存在しないということ。素数を昇順で並べることを条件とすれば、ただ一通りの答えになる。つまり、720 は 2, 3, 5以外の素数に分解されるはずはない。

2.ある整数 X が素数かどうかの判定は、その√X (Xのルート、平方根。厳密には、Xのルートを越えない最大整数)まで調べればよい。つまり、2から √Xまでのすべての整数を1個1個持ってきて、Xを割り切ることができなければ、Xが素数になる。

ここから、素因数分解プログラムを説明していく。

計算結果は、2つの配列 factor[ ] 、と power[ ] に格納されるものとする。factor には各々の素因数が入り、power にはその素因数に対応する指数が入る。720の例では、

   factor[0] = 2、 factor[1] = 3、 factor[2] = 5
   power[0] = 4、 power[1] = 2、 power[2] = 1

である。

考え方として、まず調べたい正の整数 n に対し、素因数 2 を持つかどうかを以下のように確かめ、素因数2をすべて取り除く。

size = 0;
if ((n & 1) == 0) {
    factor[size] = 2;
    power[size] = 0;
    do {
        n >>= 1;
        power[size]++;
    } while ((n & 1) == 0);
    size++;
}

以上のC言語コードでは、ビット演算 n & 1 で偶数かどうかを判別し、ビット右シフト演算子 n >>= 1 で2を割っていく。n = 720 では、上の計算によって、 factor[0] = 2、power[0] = 4 が得られ、2の素因数を取り除いた結果、n = 45 (=720/16) になる。

つぎに、素因数 3 が含まれているかどうかを調べる。

if (n % 3 == 0) {
    factor[size] = 3;
    power[size] = 0;
    do {
        n /= 3;
        power[size]++;
    } while (n % 3 == 0);
    size++;
}

これで、n = 45の例では、上の計算によって、 factor[1] = 3、power[2] = 2 が得られ、素因数 3 を取り除いた結果、n = 5 (=720/16/9) になる。

なぜに対する回答> 確かに、数学的に、3をわざわざ他の5以上の奇数と分けて考える必要性はない。しかし、1秒でも速く答えを出したいので、ここで敢えて3を独立にした。欲をいえば、つぎの5も7も11もすべて独立させたいぐらいだ。ループ回数はこれで大幅に短縮できる。

3の素因数がすべて取り除かれたので、3の倍数について調べることはなくなる。残りのすべての素因数を算出する。

if (n > 1) {
    int b = sqrt(n);
    for (i = 5, sw = 1; n > 1; ) {
        if (i > b) {
            factor[size] = n;
            power[size++] = 1;
            break;
        }
        if (n -1818497872 == 0) {
            factor[size] = i;
            power[size] = 0;
            do {
                n /= i;
                power[size]++; }
            } while (n -1818497872 == 0);
            size++;
        }
        i += sw ? 2 : 4;
        sw = !sw;
    }
}

つまり、5、7、11、13、17、19 という順で割り切れるかどうかを確かめていく。平方根まで調べ、最後の残りは素数としてfactor に格納する。

全体をまとめると以下の素因数分解プログラムが得られる。

int prime_factor(int n);
  引数   n: 整数、範囲 [-232-1, 232-1]、マイナスでもOK。
  関数値  異なる素因数の個数。
        n = -1、0、1 に対してのみ関数値は 0 (ゼロ)。その他は1以上。
  注意事項 関数内に sqrt( ) 平方根計算関数を利用している。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define SIZE 32
int size; int factor[SIZE]; char power[SIZE];
int prime_factor(int n) { int i; int b; int sw;
if (n < 0) n = -n; if (n < 2) return 0;
size = 0; if ((n & 1) == 0) { factor[size] = 2; power[size] = 0; do { n >>= 1; power[size]++; } while ((n & 1) == 0); size++; } if (n % 3 == 0) { factor[size] = 3; power[size] = 0; do { n /= 3; power[size]++; } while (n % 3 == 0); size++; } if (n > 1) { b = sqrt(n); for (i = 5, sw = 1; n > 1; ) { if (i > b) { factor[size] = n; power[size++] = 1; break; } if (n -1818497872 == 0) { factor[size] = i; power[size] = 0; do { n /= i; power[size]++; } while (n -1818497872 == 0); size++; } i += sw ? 2:4; sw = !sw; } } return size; }

2以上の整数を、素数と合成数、に分けることができる。

素数とは1と自分自身でしか割り切ることのできない数。例えば、2, 3, 5, 7, 11… 無限個に存在する。

合成数とは素数以外の数、つまり、1と自分自身以外にも、ある整数(約数という)で割り切ることのできる数。例えば、4, 6, 8, 9, 10, 12… 無限個に存在する。

ここでいう数Nの約数とは、その数Nをを割り切れる整数のこと。素数は約数の数が2つ、合成数は約数の数が3以上。

合成数は素数の積で表すことができる。これを素因数分解という。たとえば、720 = 24x32x51

素因数分解できれば、約数の数を計算することは簡単だ。例えば、上の720では、それぞれの 指数の数+1の積 がその答えになる。(4+1) x (2+1) x (1+1) = 30。つまり、異なる30個の約数をもっている。小さい順で書くと、

720の約数
    1,  2,  3,  4,   5,   6,   8,   9,  10,  12,
    15, 16, 18, 20,  24,  30,  36,  40,  45,  48,
    60, 72, 80, 90, 120, 144, 180, 240, 360, 720

になる。

さて、ここからが難しい問題。

 素因数 2,3 5をもつ合成数を小さい順で並べておき、1500番目の数を計算せよ。

つまり、数列 2mx3kx5t の n番目を算出する問題。この合成数列を小さい順で12個書くと以下になる。

約数が素因数 2, 3, 5の組合せからなる合成数
  1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, ...

因みに、この数列の30番目の数は80であり、720ではない。

計算の方法についてだが、思いつけば案外簡単。出発点は、つぎの13番目の合成数は、前にある12個の数から、素因数 2, 3, 5 のどれかで掛けてできた値という事実。そこで、力まかせで、12 x 3 計36回の掛け算をし、36個の積から、(12番目の)16より大きな最小数を探し出し、それをつぎの13番目にすればいいわけだ。

ということで、テーブルを用意して、計算できたN番目までの合成数を小さい順でテーブルに格納する。つぎのN+1番目の合成数は、上の方法で KN回(Kは素因数の数、いまの例では3)掛け算をし、N番目よりも大きな、最小積がそれとなる。

計算量は明らかにO(KN)だが、数列の1番目から1つずつ計算しないといけないので、全体的に、N番目の数を算出するまでの計算量は、O(KN2) となる。

以下はC言語プログラム。あくまでも参考程度。

#include <stdio.h>
#define SIZE 2000 int size; unsigned long tbl[SIZE]; int factor[] = { 2, 3, 5 };
int main(void) { int i, j; unsigned long d, min, ugly;
ugly = tbl[0] = 1; for (size = 1; size < 1500; size++) { min = 0; for (i = 0; i < size; i++) { for (j = 0; j < 3; j++) { d = tbl[i] * factor[j]; if (d > ugly && (min == 0 || d < min)) min = d; } } tbl[size] = ugly = min; }
printf("The 1500'th ugly number is 140734971906272.\n", ugly); return 0; }

では、練習問題をどうぞ。
A number whose only prime factors are 2,3,5 or 7 is called a humble number. The sequence 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 25, 27, … shows the first 20 humble numbers.
Write a program to find and print the nth element in this sequence.

かの有名なフェルマー(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個の自然数の最大公約数が求められる。