効率的に約数の個数を求めるアルゴリズムを考えているが、まだ四苦八苦している状態。つまり、1~500万までの整数について、それぞれの約数の数を一気に求めたい。

個々の整数なら、素因数分解して、因数の指数の積で約数の数が分かるのだが、1個1個やっているのでは、遅すぎて話にならない。

それよりも多少高速なプログラムは以下の通り。それでも数秒かかってしまう。

#define MAX   5000000
int c[MAX+10];  /* 約数の数を記録する */
void calc(void)
{
    int i, j;
    for (i = 1; i <= MAX; i++) c[i] = 1;
    for (j = 2; j <= MAX; j++) {
        for (i = j; i <= MAX; i+=j) c[i]++;
    }
}

プログラム自体は約数の定義そのものなので、面白みがないけど、圧倒的に分かりやすい。

さて、約数の数について、規則性を見るために、500万までの整数のうち、約数の個数が 1, 2, 3, …100となる整数の個数をリストアップしてみた。

約数の数が1の整数は1だけ。約数の数が2の整数は素数、つまり500万までの整数のうち、素数となる整数は348513個もある。だが、約数の数が8の整数が一番多いことに不思議さを感じている。ほかにも面白いことが沢山潜んでいるようだ。

(1,1) (2,348513) (3,331) (4,978982) (5,15) (6,189382) (7,6) (8,1118592) (9,630) (10,34670) (11,2) (12,429792) (13,2) (14,8631) (15,213) (16,689491) (17,1) (18,38545) (19,1) (20,71017) (21,89) (22,675) (23,1) (24,362994) (25,13) (26,201) (27,386) (28,16023) (29,0) (30,11399) (31,0) (32,248237) (33,21) (34,20) (35,11) (36,61776) (37,0) (38,7) (39,11) (40,50345) (41,0) (42,2722) (43,0) (44,984) (45,209) (46,0) (47,0) (48,146533) (49,2) (50,703) (51,3) (52,243) (53,0) (54,3921) (55,4) (56,9754) (57,1) (58,0) (59,0) (60,15024) (61,0) (62,0) (63,70) (64,50410) (65,2) (66,208) (67,0) (68,11) (69,0) (70,292) (71,0) (72,31359) (73,0) (74,0) (75,33) (76,1) (77,2) (78,59) (79,0) (80,15341) (81,84) (82,0) (83,0) (84,2909) (85,0) (86,0) (87,0) (88,400) (89,0) (90,1291) (91,1) (92,0) (93,0) (94,0) (95,0) (96,27891) (97,0) (98,27) (99,10) (100,700)

約数の数が与えられて、もとの整数を高速に計算するプログラムも考えたい。

最後に、約数の数が1, 2, 3, … 100となる、それぞれの最小整数をリストアップしておく。(ご注意、あくまでも500万までの整数のうち。拡大していけば、0の部分はなくなる。)

(1,1) (2,2) (3,4) (4,6) (5,16) (6,12) (7,64) (8,24) (9,36) (10,48) (11,1024) (12,60) (13,4096) (14,192) (15,144) (16,120) (17,65536) (18,180) (19,262144) (20,240) (21,576) (22,3072) (23,4194304) (24,360) (25,1296) (26,12288) (27,900) (28,960) (29,0) (30,720) (31,0) (32,840) (33,9216) (34,196608) (35,5184) (36,1260) (37,0) (38,786432) (39,36864) (40,1680) (41,0) (42,2880) (43,0) (44,15360) (45,3600) (46,0) (47,0) (48,2520) (49,46656) (50,6480) (51,589824) (52,61440) (53,0) (54,6300) (55,82944) (56,6720) (57,2359296) (58,0) (59,0) (60,5040) (61,0) (62,0) (63,14400) (64,7560) (65,331776) (66,46080) (67,0) (68,983040) (69,0) (70,25920) (71,0) (72,10080) (73,0) (74,0) (75,32400) (76,3932160) (77,746496) (78,184320)
(79,0) (80,15120) (81,44100) (82,0) (83,0) (84,20160) (85,0) (86,0) (87,0) (88,107520) (89,0) (90,25200) (91,2985984) (92,0) (93,0) (94,0) (95,0) (96,27720) (97,0) (98,233280) (99,230400) (100,45360)

1以外は、すべて偶数のようだ。

指定された約数の数となる、最小整数を求める (Smallest number with exactly n divisors) プログラム。それも面白そう。

約数の数が素数の場合は簡単に分かるが、素数以外の場合はどうなっているんだろう。

p11176.jpg

上図は整数の2進表示において、連続した1の並びの最大値の個数を表すもの。たとえば、3桁の2進数(つまり、整数0から7まで)は以下の通りだが、連続した1の並びの最大値は左から 0,1,1,2,1,1,2,3になっている。2進数101では、ビット1は2箇所にあるが、連続した1の並びの長さはそれぞれ1なので、全体としても1だ。

   000 001 010 011 100 101 110 111

そこで、同じ最大値のものをまとめて、最長値が0のものが1つ(000)、最長値が1のものが4つ(001,010,100,101)、最長値が2のものが2つ(011,110)、最長値が3のものが1つ(
111)になり、つまり 1 4 2 1という並び(上から3段目)が得られる。

ビット数をn、連続した1の並びの最長値をkとすると、整数それぞれにnとkはきまった対応関係になる。それを関数T(n,k)で表すと、つぎのような漸化式が成り立つ。

p11176-1.jpg

ただ、値がすぐに大きくなっていくので、多倍長整数か、実数を利用すること。

問題 UVa #11176 – Winning Streakを解くのに、こういう知識が勉強になる。

オイラーの関数φ(n)の値が与えられた時の、nを求める問題。つまり、

  φ(n)=x のxに対し、n=φ-1(x)

を考える。

整数の素因数分解と同様、逆関数の計算はとても難しい。だから問題として出されているかもしれない。

1. 元問題 — UVa #11073 – Euler’s Totient Function
  xの値の範囲は 1<=x<=1000000000。
nの値が複数ある場合、昇順で表示する必要もある。

2. 実例
  x=2に対して、n=3,4,6
  x=4に対して、n=5,8,10,12
  x=583925760に対して、書ききれないが、異なるnは19639個存在する。

3. 考え方
 xの値が小さければ、例えば、100万までなら、φ(n)の値をすべてテーブルに計算登録し、xでlookupする方法も考えられるが、10億もの巨大xではこの方法は無理だろう。

幸い、元問題の入力値xの個数はそれほど多くないらしく、素直に個々にxの逆関数を計算していくことにした。

自分では以下の漸近式を使った。

p11073.jpg

つまり、2から約32000(10億の平方根)までのの素数をひとつひとつ持ってきて、xを(p-1)*pkで割り切れるかどうかを確かめる。割り切れるのであれば、x/((p-1)*pk)を新たなxとして、その逆関数を再呼び出して分解していく。計算できたnをひとつひとつテーブルに登録し、最後にソートをかけ画面出力する。

ちなみに、本方法では、途中で打ち切ることがあっても、nの値がダブって算出されることはない。

なお、x=1に対しては、n=1,2と結果表示し、他の奇数のxに対しては、解がなしと表示する。

4. 実例で説明
 基本的な考え方は上記の通りだが、細かいところについては以下の実例をよく見て、理解すること。

x=4について考える。

(1-1) まず最初の素数p=2を持ってくる。k=0とし、(p-1)*pkを見ていく。なお、k=1については(1-2)の所、k=2については(1-3)の所で調べる。
  p=2, k=0では、(p-1)*pk=1, pk+1=2, x/((p-1)*pk))=4になり、φ-1(4)→2*φ-1(4)になる。
 (1-1-1) 素数2は上で既に使ったので、ここでは次の素数3で調べる。p=3,k=0のもとでは、(p-1)*pk=2, pk+1=3, 4/((p-1)*(pk)=2になり、2*φ-1(4)→6*φ-1(2)になる。
  (1-1-1-1) つぎの素数5で調べる番になっているが、2を(p-1)*pk=4*5kで割り切ることはありえないので、一旦ここで打ち切り。
 (1-1-2) つぎの素数5で調べる。p=5,k=0のもとでは、(p-1)*pk=4, pk+1=5, 4/((p-1)*(pk)=1になり、2*φ-1(4)→10*φ-1(1)になる。φ-1(1)は1であるので、最初のnが見つかり、n=10。
(1-2) (1-1)の続きをやる。つまり、p=2, k=1。
  p=2, k=1では、(p-1)*pk=2, pk+1=4, x/((p-1)*pk))=2になり、
φ-1(4)→4*φ-1(3)になる。
 (1-2-1) 次の素数3で調べる。p=3,k=0のもとでは、(p-1)*pk=2, pk+1=3, 4/((p-1)*(pk)=2になり、
4*φ-1(2)→12*φ-1(1)になり、2番目のn=12がみつかる。
(1-3) (1-2)の続きをやる。つまり、p=2, k=2。
  p=2, k=2では、(p-1)*pk=4, pk+1=8, x/((p-1)*pk))=1になり、
φ-1(4)→8*φ-1(1)になり、3番目のn=8がみつかる。

(2-1) 2番目の素数p=3を持って来て調べる。
  p=3, k=0では、(p-1)*pk=2, pk+1=3, x/((p-1)*pk))=2になり、
φ-1(4)→3*φ-1(2)になるが、つぎの素数5ではそれ以上nを見つけることはないので、打ち切る。

(3-1) 3番目の素数p=5を持って来て調べる。
  p=5, k=0では、(p-1)*pk=4, pk+1=5, x/((p-1)*pk))=1になり、
φ-1(4)→5*φ-1(1)になり、4番目のn=5がみつかる。

ということで、なかなか理解するのは大変かもしれないが、個々の素数をつぎつぎと持ってきてテストし、つぎの深いステップでは、同じ素数は使ってはいけないところが味噌。また、同じステップのところでは、kの値も0からインクリメントして調べる。いわば、3次元的に調べるイメージになる。

p11073-1.jpg

5. 結果
あまりよい計算方法がないらしく、素数をひとつひとつで確かめていき、効率の面を心配していた。結果的にランク1位の実行速度になり、皆も相当苦労していたのだなと実感した。

まだまだ研究余地の残っている問題だと思う。

オイラーの関数φ(n)とは、1からnまでの整数のうちnと互いに素なものの個数のこと。

たとえば、1,2,3,4,5,6のうち6と互いに素なのは、1,5の二つなので、φ(6)=2。素数pについては、1からpまでの整数のうち互いに素でないのはp自身だけなので、φ(p)=p-1となることは直ちに理解できよう。なお、φ関数は英語では、totient functionと呼ばれている。

1からのφ関数の値は以下の通り。

 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18, 8, 12

φ(n)の計算に重要なもの(性質)をいくつかピックアップしておく。

totient.jpg

整数の素因数分解と同様、個々のφ(n)の計算と、多数のnのφ(n)の計算とを分けて、その計算法を考えよう(数学ではそんなことを考える必要性は全くないが)。

個々のφ(n)の計算では、nの素因数分解を利用する方法が最も効率的だろう。nに含まれている素因数pkが分かれば、(n/Pk)*(pk-1)で計算していけばよし。

多数のnをまとめてそのφ(n)を計算するなら、以下のプログラムを利用しよう。φ(n)の性質(互いに素)を利用した方法。

#define MAX 1000000
int phi[MAX+10];
char prime[MAX+10];
void euler_phi(void)
{
    register int i, j;
    for (i = 0; i < MAX; i++) phi[i] = i;
    phi[2] = 1;
#if 0
    memset(prime, 0, sizeof(prime));
#endif
    for (j = 4; j < MAX; j += 2) prime[j] = 1, phi[j] -= phi[j] >> 1;
    for (i = 3; i < MAX; i++) {
        if (prime[i]) continue;
        phi[i] -= phi[i] / i;
        for (j = i+i; j < MAX; j += i) prime[j] = 1, phi[j] -= phi[j] / i;
    }
#if 0
    for (i = 1; i <= 100; i++) printf("(-452953008,-2082935264) ", i, phi[i]);
    printf("\n");
#endif
}

これでいっぺんに、MAX(100万)までの値が計算される。

問題意識は UVa #11480 Jimmy’s Balls による。

つまり、自然数Nを3つの自然数 x, y, z に分割する時の分割数(組合せ数)。ただし、x < y < z。自然数とは1以上の整数のこと。

たとえば、N=10に対して、(x, y, z)は、(1,2,7), (1,3,6), (1,4,5), (2,3,5) の4通りしか選べないので、分割数は4になる。

p11480.png

割とよく遭遇しそうな問題だが、一般解は以下の数式で表される。

  分割数 f = 3p(p-1)+1  if m = 0

      = 3p(p-1)+mp  if m > 0

ただし、p=N/6 (6で割った時の商)、m=N%6(6で割った時の余り)。

いくつかの検証データ(N,f)を書いておく。

 (5,0) (6,1) (7,1) (8,2) (9,3) (10,4) (11,5) (12,7) (13,8) (14,10) (15,12) (100,784) (500,20584)

一般化して、Nをn個に分割する時の数式ってあるのだろうか。

2以上の整数nについて、その素因数の和sf(n)の計算法を考える。

例えば、2 = 2 なので、sf(2) = 2。
6 = 2*3 なので、sf(6) = 2+3 = 5。
8 = 2*2*2 なので、sf(8) = 2+2+2 = 6。

2から30までの整数についての sf(n)の値は以下の通り。
2,3,4,5,5,7,6,6,7,11,7,13,9,8,8,17,8,19,9,10,13,23,9,10,15,9,11,29,10.

#include <stdio.h>
#include <string.h>
#define MAX 500000
#define ROOT_MAX 707
int sf[MAX+1];
void sieve(void)
{
    register int i, j, k;
    /* memset(sf, 0, sizeof(sf)); */
    for (j = 2; j <= MAX; j <<= 1) {
        for (i = j; i <= MAX; i += j) sf[i] += 2;
    }
    for (k = 3; k <= MAX; k += 2) {
        if (sf[k] == 0) {
            for (j = k; j <= MAX; j *= k) {
                for (i = j; i <= MAX; i += j) sf[i] += k;
                if (k > ROOT_MAX) break;
            }
        }
    }
}

計算できる最大値 MAX は適宜変えること。ROOT_MAXはMAXの平方根。関数の実行が終了した時点で、配列sf に2以上、MAXまでの素因数の和が格納される。

2ヶ月もこのブログを更新していない.忙しかった.問題のほうはさすがにここに来て解くペースがだいぶ落ちてきた.1日1問でも厳しい日が続いている.でも諦めたくない.

さて,上の計算を見てみよう.2の3乗の4乗の5乗の下一桁はなんだ,との問いかけだけど,答えられる? ちなみに,数の大きさは対数を取って見積もったら,489桁であることが分かる.

べき乗のべき乗のべき乗の...で表した巨大整数の剰余計算は今回のテーマ.この種の一般化した問題はどういうわけか,正面から答える書物は少ない.誰でも思いつく問題なのに,計算の仕方を言及する書物はあまりみかけない.だからプログラミングコンテスト問題として出されたのかもしれない.

問題の正確な仕様は ここ を参照すること.つまり,与えられた整数 a1, a2, a3, … , aN および 法 m に対し,次の式を計算すること.

   a1^a2^a3^…^aN mod m

なお,各変数の取る範囲は 1 <= ai <= 1000, 1 <= N <= 10, 2 <= m <= 10000.

では,まず重要な定理を紹介しよう.

<離散対数定理 (Discrete logarithm theorem)>
 gとmとが互いに素であれば,

    gx ≡ gy (mod m)

が成り立つための必要十分条件は x ≡ y (mod φ(m)) である.なお,φ(m) はオイラー関数 (mと素であり,しかも m より小さい自然数の数を表す.) </離散対数定理>

離散対数定理に従えば,g と m が互いに素という条件を満足さえすれば,べき乗の高さは1段下げることができる.例えば,

    275 (mod 5)

を計算するときに,値を 2x と設定すれば,2と法5は互いに素なので,

    2x ≡ 275 (mod 5)
    x ≡ 75 (mod 4)  ← φ(5)=4

というふうになり,2段あるべき乗のところが1段に下がり,電卓でもつぎのように値を計算することができる.つまり,
    75 = 16807
    4で割った時のあまりは 3. ということで,x ≡ 3.

さらに,x ≡ 3 をもとの式のところに戻すと, 275 ≡ 2x ≡ 23 ≡ 8 ≡ 3 (mod 5).つまり,計算結果は 3 である.

このように,再帰手法を使えば,べき乗の高さは10段であろうと,繰り返し処理し最終的にべき乗の高さが1段までに下がり,計算できる訳だ.

さて,べき乗の底 (つまり g )と法 m とが互いに素でないときに,どう計算すれば良いのか.整数論関係の書物はこのことについてあまり言及しない.めんどくさいから書かないか,一般論としての議論は難しいか,理由は分からない.しかし,プログラミング問題は容赦なく,われわれの嫌がる所を聞いてくる.人を困らせるのがプログラミング問題の目的だから?

法が合成数の場合,問題解決に良く役立つのは中国人剰余定理 (Chinese remainder theorem) というもの.ここでもこの定理を使ってみよう.

<中国人剰余定理 (Chinese remainder theorem)>
 連立1次合同式

    x ≡ b1 (mod m1)
    x ≡ b2 (mod m2)
    …
    x ≡ bn (mod mn)

の場合、 m1, m2, … mn互いに素であれば、

    m = m1m2…mn

を法として、ただ一つの解がある。 また、

    (m/mi) xi ≡ 1 (mod mi)

の解 xi を求めることができ,

    x ≡ b1 (m/m1) x1 + b2 (m/m2) x2 + … + bn (m/mn) xn (mod m)

とすれば、この x はもとの合同式すべてを満足する解となる。</中国人剰余定理>

この中国人剰余定理に従い,問題の法 m をいかに分割するかをつぎのように考える.

法 m を 互いに素な m1 と m2 との積に分け,m1 を底 g とも互いに素であるようにする.つまり,

    m = m1 m2
    m1 ⇔ m2 互いに素.このことを整数論では (m1, m2) = 1 と表記されることもある.
    m1 ⇔ g 互いに素,つまり, (m1, g) = 1

とする.それで,

    gx (mod m1) については, g と m1 は互いに素なので,離散対数定理により計算する.
    gx (mod m2) については,以下のように考える.

底 g と法 m2 との最大公約数 gcd(g, m2) を d (d > 1 であることは明らか) とすると,gx は dx で割り切れる.もし,dx は m2 で割り切れるなら,gx も m2 で割り切れることになる.つまり,

    gx ≡ 0 (mod m2)

となる.では,dx は m2 で割り切れるだろうか.

m2の大きさは本問題の仕様から分かるように,5000以下であり,素因数 d を含むならば,d の数は最大でも 12 である.なぜなら,

    m2 < m/2 = 5000 < 8192 = 213 <= d13

という不等式が成立するからである.一方,gx のほうは, x の値がほとんどの場合 12 以上であろう.勿論,このことはプログラムの中できちんと確認しなければいけないが,x の値が 12 以上と分かれば,gx は 法 m2 で割り切れることになる.

これらの議論をつぎのようにまとめる.

 剰余計算 gx (mod m) に対し,
  1) x の値が12以下であれば,直接,剰余計算を行い,終了.
  2) g と m とが互いに素であるかを確かめる,素であれば,離散対数定理により計算し終了.
  3) 互いに素でなければ,g と m との共有素因子をすべて m から取り除き,残った値を m1 とする.
  4) gx (mod m1) を離散対数定理により計算し,値を b1 とする.
  5) 中国人剰余定理により,gx ≡ b1(m/m1)k (mod m) で最終結果を算出する.ただし, k は (m/m1)k ≡ 1 (mod m1) より得る. 

上の計算は勿論べき乗の各高さに繰り返し適応すべきである.

では,冒頭の計算を例としてここで取り上げる.

2345 (mod 10) では, 2 と 10 が互いに素でないので,法 10 を 5 と 2 に分割する.

 ⇒ 2345 (mod 5) を離散対数定理により計算する.つまり, 345 (mod 4) を計算する.

  ⇒ 345 (mod 4) では 3 と 法 4 は互いに素であるので,続いて,また離散対数定理を適応する.つまり,45 (mod 2) を計算する.なお,法 2 はオイラー関数 φ(4)=2 とのことから.

   ⇒ 45 (mod 2) では,5 が12より小さいので,直接剰余計算を行い,答は勿論 0 である. 2進数表現では,偶数の下1桁は常に 0 なので.

  ⇒ 値を戻し,345≡30≡1 (mod 4) となり,

 ⇒ 2345≡21≡2 (mod 5) となる.

さらに,中国人剰余定理により,もとの計算 2345≡2*(10/5)*3≡12≡2 (mod 10) となる.

つまり,下一桁は 2 である.なお,3 は (10/5)k≡1 (mod 5) の解 k = 3 による.

以下は本問題を解くC言語プログラム.きちんと整理してなくて,見苦しいところもあろうが,整数論の多くの知識が詰みこまれていて,勉強にはなるだろう.

C Source Solution for UVa P10692 (Huge Mod) Online Judge

少し説明すると,

剰余計算関数 int mod(int a, int p, int m) は,ap (mod m) を直接高速計算する.

最大公約数計算関数 int gcd(int a, int b) は整数 a と b との最大公約数を算出する.

素因数分解関数 int prime_factor(int n) は整数 n (n > 1) のすべての異なる素因数を算出し,外部配列 factor[10] に格納する.異なる素因数の数は関数のリターン値に入っている.なお,高速化するため,素因数 2 と 3 は関数の中では,特別に処理している.

オイラー関数の計算 int phi(int n) (n > 1) は,上の素因数分解関数を利用して,オイラー関数 φ(n) の値を算出する.

合同式を解く関数 int congruence(int a, int m) は,a*x ≡1 (mod m) を解き, x の値をリターンする.

本問題の最重要関数 int huge_mod(int a0, int *a, int size, int m) は a0^a[0]^a[1]^a[2]^…^a[size-1]  (mod m) を計算する関数.高速化するための工夫箇所もあるが,上述したアルゴリズムが含まれている.

メイン関数 int main(void) では入出力を司る部分.工夫として,1以降の入力をすべて打ち切るところにある.1のべき乗はつねに1だからである.

プログラミングコンテストの会場では,これらのソースリスト(整数論の基礎知識に関する部分)が手元になければ,数時間で本問題を解くことは大変困難であろう.

ちなみに,本プログラムの計算時間は 0.002秒とUVa判定プログラムが 報告 している.ゼロ秒には到達できていないが,遅い計算法ではなかろう.

与えられた奇数 n が素数であるかどうかを高速に判定することはなかなか難しいが,比較的高速な判定法に,Miller-Rabin テストと呼ばれるものがある.

次のように,ブール関数 Suspect(b, n) をまず定義する.b は基底を表す整数,n は奇数.

Function Suspect(b, n), where b is an integer that is called the base, and n is an odd integer. It returns one of the boolean values TRUE or FALSE.

  • Let t be the highest power of 2 so that 2t devides n-1 and u be the biggest odd integer that devides n-1.
    This means we can write n-1 = 2tu.

  • Let x0 be bu mod n.
  • For all i from 1 to t, let xi be (xi-1)2 mod n.
  • If, for any i from 1 to t, xi = 1 and xi-1 <> 1 and xi-1 <> n-1, then return FALSE.
  • else, if xt <> 1, then return FALSE.
  • else return TRUE.

関数の値がFALSE であれば,n は素数でないことは直ちに言えるが,一方,関数の値がTRUE の場合には,極まれに n は素数でないこともある.

しかし,異なる基底で複数回計算し,それぞれの関数値がすべてTRUEであれば,素数であることが保証されることも知られている.

具体的に,1000000(百万)までの 奇数 n については,Suspect(2, n) と Suspect(3, n) とで計算した関数の値が両方とも TRUE であれば,n が素数であることが保証されている.

さらに,すべての32ビット 奇数 n (つまり, n < 232) までに拡大したければ,Suspect(2, n), Suspect(7, n) および Suspect(61, n) の3つの異なる基底で計算した関数値が すべてTRUE であれば,n が素数であることが保証されていることを利用すればよい.

以下はC言語プログラム.

unsigned long long calc_pow(int a, unsigned k, unsigned n)
{
    unsigned long long p;
    unsigned bit;
p = 1; for (bit = 0x80000000U; bit > 0; bit >>= 1) { if (p > 1) p = (p*p) ; if (k & bit) p = (p*a) ; } return p; }
int suspect(int b, unsigned n) { unsigned i; unsigned t, u; unsigned long long x0, x1;
u = n-1; t = 0; while ((u & 1) == 0) { t++; u >>= 1; } x0 = calc_pow(b, u, n); for (i = 1; i <= t; i++) { x1 = (x0*x0) ; if (x1 == 1 && x0 != 1 && x0 != n-1) return 0; x0 = x1; } return x1 == 1; }
int prime_test(unsigned n) { if (n <= 1) return 0; /* 0, 1 */ if (n == 2) return 1; /* 2 */ if (n == 3) return 1; /* 3 */ if ((n & 1) == 0) return 0; /* other even integers */
if (n <= 1000000) { if (!suspect(2, n)) return 0; if (!suspect(3, n)) return 0; } else { if (!suspect(2, n)) return 0; if (!suspect(7, n)) return 0; if (!suspect(61, n)) return 0; } return 1; }

上記のプログラムでは,
  calc_pow 関数は ak mod n を計算する.
  最後の prime_test 関数 は任意の32ビット内の整数について,それが素数であるかどうかを判定する.

繰り返しになるが,Suspect関数を1つだけ使って素数判定すると間違った結果を出すこともあるが,異なる基底を組み合わせて複数回計算するば,素数判定は正しくできる.

計算時間もかなり高速だと思われます.

<メモ>
本記事は オンラインプログラムコンテスト問題 10956 Prime Suspect による.

階乗 N! を素因数分解して得られる素因数の数を P(N) とすると、P(N) に対応する最小の N を求めて欲しいのが ACM #10856 Recover Factorial という問題。

例: 4! = 4*3*2*1 = 2*2*3*2 ⇒ P(4) = 4
   6! = 6*5*4*3*2*1 = 3*2*5*2*2*3*2 ⇒ P(6) = 7
   10! = 10*9*8*7*6*5*4*3*2*1 = 2*5*3*3*2*2*2*7*2*3*5*2*2*3*2
     ⇒ P(10) = 15

N と P(N) との対応関係は最初のところでは以下の通り。
   N  0 1 2 3 4 5 6 7  8  9  10 …
  P(N) 0 0 1 2 4 5 7 8 11 13 15 …

注意して欲しいのは、P(N) の値は連続しておらず、3, 6, 9, 10, 12 等にはならないこと。また、問題の入力データは、0 <= P(N) <= 10000001 の範囲内であること。

さて、試行錯誤の末、つぎの解き方にたどり着いた。

 (1) 1 ~ 2703663 までの N を素因数分解し、それぞれの素因数の数を 配列 F[1..N] に格納する。(素数のエラトステネスのふるいを活用。)
 (2) (1)で求まった F配列 を累計して、P(N) に関するテーブルを作成する。漸近式: P[N] = P[N-1] + F[N]、P[0] = 0。
 (3) バイナリサーチで、入力値として指定された P(N) を (2)のテーブルから探し出す。

(1)のところの値 2703663 は、最大のP[N] = 10000001時のN である。また、素数のエラトステネスのふるいを活用するところは、自分なりの工夫を多少したつもりだが、全体的に素直な解き方だろう。しかも、計算結果は満足のいく 0.455 秒、現時点ではベスト1を記録

参考のため、C言語プログラムを載せておく。エラトステネスのふるい(の変形)と、バイナリサーチの部分は、多少ためになるかもしれないと思ってるので。
p10856.c

<問題の詳細>
ACM問題集にある問題: 10168 Summation of Four Primes (4素数の和) を例にプログラムを考えてみたい。

問題の内容は、与えられた整数 N (N<=10000000)に対し、Nを4つの素数の和で表せというもの。もしNは4つの素数の和で表現できなければ、"Impossible."で答えろという。 また、4つの素数について、条件は一切なく、好きな組み合わせで良いという。なお、Nの個数に関し、その上限はとくに問題文に記述されていないが、数十個、数百個のNを素早く計算できないといけないだろう。

<自動判定の話>
つくったプログラムが正しいかどうかは、ここに貼り付けて送信すれば、その場で判定される。フォームに問題ナンバー (Problem #)(この問題では 10168)、登録者ID (ID with suffix)(サーバに登録すればその場で発行してくれる。登録料は勿論必要なし。誰でも参加できる。)、プログラムに使用した言語(C/C++/Pascal/Java の4つに限定されているが、C/C++が無難。Perl 等のスクリプト言語は残念ながら採用見込無) を記入または選択して、フォームの下半分の Source code: のところにソースプログラムを貼り付けるだけでOK。なお、コメント (Comment)欄の記入は必要なし。

自動判定の原理は簡単。送信したプログラムがその場でコンパイルされ、実行される。出てきた答えと用意された答えとの照合が行われ、ひとつでも答えが合わないと Wrong Answer(略してWA) を食らう。無事パスしたら、Accepted (略してAC)をもらう。自分の成績簿にその結果がその場で反映される。プログラム問題の多さ、自動判定レスポンスの良さ等と相まって、世界中(とくにコンピュータサイエンス関係の大学生)に大人気。多くのひとが日夜、多くの難問にチャレンジしている。残念ながら、日本では認知度がそれほど高くないようで、上位100以内にランキングされる日本人はとても少ない。

ちなみに、ほかの多くのプログラムの自動判定結果を眺めるのも楽しい。実行に使用したCPU時間、メモリ量、言語、登録者、問題ナンバーがこれで判る。

<実行時間の厳しさ>
さて、大変重要なことだが、許される実行時間は僅か10秒まで。10秒過ぎると実行が無条件に打ち切られ、Time Limit Exceeded (略してTLE) を食らう。しかも、サーバのCPUはペンティアムIII。どちらかというと遅いPC。ということで、10秒以内に数十個か数百個のNを計算し、正解を出さないといけない。数分間、数十分、数日かけてやっと答えの出るプログラムはACされない。プロの世界はやはり厳しいものだ。

<休憩>

<問題の整理>
問題にもどって、どう解答するかをまず見通しを立てよう。つぎの事実を知ると心強い。

ゴールドバッハの予想: 4以上の偶数は2つの素数の和で表すことができる。証明はされていないが、コンピュータの計算可能な範囲内では数学者達が確認済みなので、信じていいだろう。

どういうことかというと、4 = 2+2、6=3+3、8=3+5 のように、等号の左に4以上の任意の偶数が与えられると、等号を満たす2つの素数は必ず存在するということだ。(2, 3, 5 は素数)

ちなみに、大数学者オイラー宛のゴールドバッハからの手紙には、正確にいうと、つぎの等価予想が書いてあった。「5よりも大きい自然数は3つの素数の和であらわすことができる。」

これらの事実をまとめると、問題解決の戦略は以下のようにできる。

N が奇数ならば、なるべく大きな素数(Pとする)+2+残りの偶数1つ(Cとする)、に分解、
N が偶数ならば、なるべく大きな素数(P)+3+残りの偶数1つ(C) に分解することにしよう。

つまり、C=N-P-(2か3) で、N と P から C を算出。

2以外の素数はすべて奇数なので、N を以上のように奇数と偶数に分けたわけだ。また、C はゴールドバッハの予想より、2つの素数の和にさらに分解できるので、これで、N に対し、合計4つの素数が得られるわけだ。

P をなるべく大きな素数にするのは、C の値を小さくしたいから。C の値が5000以内なら、すべての組合せを事前に計算しておくこともできる。

5000という間隔で見積もると、N の最大値=5000*2000 になるので、約 2000個のP と 約2500対(偶数だけでいいので、5000の半分で済む)の和がCになる素数ペアを用意しておけば、ほぼゼロ秒(計算量の世界では O(1) という)で計算できる。

なお、2000個のP から正しいものを1つ選び出すのは難しくないだろう。2000個のPを昇順で並べておき、N/5000 番目のPを使えばよいのだから。

<ルックアップ・テーブルの多用>
上の文章に、事前に用意する云々の言葉が数回出ていたが、10秒という難関を突破するのに、とくにここのACM問題解答プログラムでは多く用いられている技法のひとつ。つまり、結果を事前に、手元のコンピュータを使って計算しておき、自動判定のためのプログラムソースに計算結果をテーブルとして埋め込むのだ。

(補足説明) ただ、もうひとつの厳しい判定条件にも気をつけないといけない。つまり、ソースプログラムの 長さは最大40KB まで。なんでもかんでもテーブルにしてはいけないという制限条件。

手元のコンピュータでは数分でも、数日でも構わないので、弾き出された結果を、テーブルに加工しておき、サーバに送り込む。入力データに合う結果をテーブルから探し出せばいいので、インチキといえばインチキ。だが、その方法でしか10秒クリアできない問題も結構あるという。

もうすこしわかり易く書くと、手元のコンピュータで、2000個のP、及び、Cに対応する2500対の素数ペアを算出しておく。それらをテーブルに書き込み、判定用「インチキ」プログラムを特別につくっておく。サーバに送り込まれた判定用「インチキ」プログラムでは、入力データの N に対し、N/5000 番目のPを使う。つぎに、C=N-P-(2か3)という式でCを算出し、テーブルから、Cに対応する素数2つを選び出す。それで、4つの素数が揃ったことになる。それを答えとして出力すればOK。プログラムの計算時間は限りなくゼロに近い。ベストタイムが狙えるかも。

ここでいうテーブルの実体は勿論配列。例えば今回の偶数を2つの素数の和で表す場合は、配列のインデックスはその偶数/2の値、配列の値は2つの素数。

(補足説明) ただ、よく考えたら、値が5000以内の素数判定なら、ノーマルな方法でも余裕で1秒以内でクリアできるので、Cの計算についてはテーブルなしでも全然OK。しかし、Pについてはやはりテーブルに用意したほうが無難だろう。

最後に境界条件をまとめよう。

最小の素数が2なので、N は7以下(ゼロやマイナスの整数も含む)ならば “Impossible.”と答えないといけない。

8以上のN はすべて4つの素数の和で表すことができるが、
   N=8=2+2+2+2
   N=9=2+3+2+2
   N=10=3+3+2+2
では、Cに偶素数ペア 2+2 を使うので、特殊ケースと考えていいかもしれない。11以上のNなら、Cは奇素数ペアの和で表現できる。

(補足説明) 偶素数:2 ひとつだけ。奇素数:3,5,7...、無数にある。偶素数+奇素数=奇数になるので、偶数になるには、偶素数同士、つまり、2+2 の一通りしかない。

<残りはプログラミングのみ>
プログラムの見通しはこれで立てたので、残りの作業はプログラミングのみだ。アイディア的に面白いのは実はここまで。プログラミングはいわば確認作業のようなもの。

テーブルのデータが多いので、面白くないと思うが、参考のため、ソースプログラムを掲載しておく。

<<素数テーブルの生成プログラム>> p10168-pre.c

<<解答プログラム>> p10168.c
 最後の偶数については、強引に2つの素数に分解している。

<<解答プログラム2(高速版)>> p10168-2.c
テーブルに改良したバージョン。最後の偶数Cに対応する素数ペアは片方判れば良いので、片方だけをテーブルにしてある。ベストタイムを狙ったつもりだが、結果は 0.016秒、4位になっている。これ以上速くするには、printf 関数も自作しないといけないかもしれない。

<<解答プログラム3(高速版2)>> p10168-3.c
printf 関数を使わないバージョン。送信して確認した結果、実行時間は 0.010秒、3位に浮上。これ以上速くすることは私に無理。

<ちょっとした遊び>
問題文の最後に、さりげなく変わった文章が載っている。

「You can fool some people all the time, all the people some of the time but you cannot fool all the people all the time.」 

そうかもしれないね。つぎのようなことも言えるかも。

「We can love some people all the time, all the people some of the time but we cannot love all the people all the time.」