Cinderellaで数学：素数
エラトステネスのふるい
ｎ以下の素数をすべて見つける方法に、「エラトステネスの篩（ふるい）」というものがあります。何も難しいことではなく､ある素数の倍数を消していくだけのことです。たとえば、100以下の素数でしたら、２から始めて、３，５・・・の倍数を消していけばよいのです。難しいのはそれを抽象的に考えた場合ですが、高校までの数学ではそこまでは扱いません。
さて、エラトステネスの篩で素数をリストアップするプログラムを作ります。「リスト」処理ができるCindyScriptならではのエレガントな方法も紹介します。
まずは、基本形。2から始めて$${\sqrt{n}}$$ まで、倍数を消すことを「繰り返す」方法です。なお，n までやらなくても，$${\sqrt{n}}$$までやればいいことが「知られて」いますのでそれを使います。その証明は略します。
はじめに、２からｎまでの自然数を用意します。それを「リスト」にします。
nlist = 2 .. n
で２からｎまでの自然数のリストが作れます。ピリオドが2個です。具体的に使う場合はｎに目的の数を代入します。n=20 としてprintln()で表示してみましょう。
エラトステネスのふるいでは、このリストの中から「倍数」を消していきます。その手順を具体的に考えてみましょう。
(1) ２の倍数を消す
(2) ３の倍数を消す
(3) ４は、(1)で消えているので、次の５の倍数を消す
(4) つぎの６は(1)で消えているので、次の７の倍数を消す
これを続けて、20までやればおわり。実際には20はすでに消えていますが。
では具体的に考えていきましょう。
(1)２の倍数を消す
(1-1) リストからk番目の要素をひとつ取り出す。
(1-2) それが２の倍数かどうか調べる。
(1-3) ２の倍数であればそれをリストから消す。
これを繰り返します。「繰り返す」ということは、(1-1)のｋを 1,2,3 ・・ と変えていくことになります。ところがここで問題が生じます。
具体的に考えてみましょう。
２，３，４，５，６，７・・
で、２番目の要素３を取り出して２の倍数かどうか調べます。２の倍数ではないので消しません。したがって、リストはそのままです。
次に３番目の要素４を取り出して２の倍数かどうか調べます。２の倍数ですので消します。すると、リストは次のようになります。
２，３，５，６，７・・
すると、「次に取り出す」のは４番目の要素ではなく、やっぱり３番目の要素です。
ややこしくなりそうですね。対策を考えましょう。
＜その１＞
ｋ番目を消したら次もｋ番目、そうでなければ次は（ｋ＋１）番目とする。「もし〜ならば〜」なので、条件判断の if を使えばよい。
＜その２＞
直接消してしまわないで、もうひとつリストを用意して、そちらに消したかどうかの情報として，消したら０，残っていれば１を書いていく。
では、やってみましょう。まず＜その１＞から。
n = 20;
nlist = 2 .. n;
a = 2; // ２の倍数を消す
k = 1; // 取り出すリストの要素番号
while(k <= length(nlist),
b = nlist_k;
if(mod(b, a) == 0,
nlist = remove(nlist, b);
,
k = k + 1;
);
);
println(nlist); 結構ややこしいですが、動作を追ってみましょう。while のところです。
始めは k は 1 です。nlist_k で1番目の要素を取り出して b とします。1番目は2なので，b=2 です。すると，mod(b,a) は割り切れて0 になるので，nlist から remove で 2 を削除します。このとき，kの値は変わりませんので，繰り返すときにはやはり1番目の要素を取り出します。しかし，さきほど2を削除したので，リストの1番目は３です。
すると，こんどは2で割り切れないので，リストは変更せず，k = k + 1 で，kは2になります。2番目の要素は4です。4は2で割り切れるのでリストから削除します。このとき nlist は[3,5,6,7,8,9,・・・,20] になっています。
kの値は変わらないので，2番目の要素 5 をbに代入して・・・ と続きます。
これで，２の倍数が削除されて奇数だけが残ります。
次に，＜その２＞です。
直接消してしまわないで、もうひとつリストを用意して、そちらに消したかどうかの情報として，消したら０，残っていれば１を書いていきます。
n = 20;
nlist = 2 .. n;
flag = apply(2 .. n, # = 1); // 消したかどうかのリスト はじめすべて１
a = 2;
repeat(n - 1, k, // nlist の要素の数は n-1 個
b = nlist_k;
if(mod(b, a) == 0, // a(今は2)で割り切れたら
flag_(i + k) = 0; // flag を0にする
);
);
println(nlist);
println(flag);動作結果です。２の倍数のところが0になっています。
このあと
(2) ３の倍数を消す
(3) ４は、(1)で消えているので、次の５の倍数を消す
(4) つぎの６は(1)で消えているので、次の７の倍数を消す
をやっていきます。
＜その１＞ですと，nlist に残っているのは奇数だけなので，1番目の要素 3 を取り出して，３の倍数を削除していけばいいですね。
＜その２＞ですと，nlist はそのままですが，残っている数は flag を参照すればわかるので，次に残っている数（2番目の要素である３）で同じようにすればいいですね。
やってみてください。
リスト処理の関数を使って素数のリストを作る
さきほどやったエラトステネスの篩では、while または repeat の繰り返しを使って、素数だけを残していきます。しかし、CindyScriptにはリストを処理する関数が用意されており、これを使うとプログラムが驚くほど簡単になります。
約数のリストを作る
まず、ある数の約数のリストを作ることを考えます。あるｎに対し、1からnまでの数について、nを割りきる数を選び出します。リストの中からある条件をみたすものを選び出すには select という関数を使います。次の１行で、nの約数のリストが yaku に作られます。
yaku = select(1 .. n, mod(n, #) == 0)
select(list, expr) という関数は、list のリストから、expr の条件を満たすものを選び出してリストにします。今の場合、1 .. n が対象となるリストになります。expr の条件は mod(n, #) == 0 ですが、#はリストの各要素を指します。select()関数は、リストの最初から順に要素を#へ代入し、条件を満たすかどうかを判定していきます。その結果、1からnまでの自然数のうち、nの約数が選ばれて、リストyakuができます。
次に，これを関数として定義します。
yakusuu(n):= select(1 .. n, mod(n, #) == 0);
yakusuu(20) とすると20の約数のリストができます。
素数の判定
1とn以外に約数を持たない数nが素数です。nの約数のリストを作って、その要素の個数が2(1とnだけ)ならばnは素数であると判定できるわけです。たとえば、191が素数かどうかを判定してみましょう。
これを応用すれば100までの素数をリストアップすることができます。使うのは、やはりselect()です。つまり、2から100までの自然数のリストを作っておき、そこから、素数であるものをselectすればいいのです。プログラムは次の3行だけです。
yakusuu(n):= select(1 .. n, mod(n, #) == 0);
primelist = select(2 .. 100, length(yakusuu(#)) == 2);
println(primelist);
このように、リストを使うと繰り返し処理が簡単にできます。
素因数分解
与えられた数を素因数分解します。2から始めて、順に素数で割っていきます。割り切れたら、その商をさらに同じ素数で割ります。割り切れなかったら次の素数で割ります。これをｎ以下の素数のうち一番大きな数まで繰り返しますが，割り切れて商が1になったらそこで終わります。
(1) n以下の素数のリスト primelist を作る
(2) primelist の1番目の要素２から始めて上の手続きを行なう。そのとき、ある素数で何回割り切れたか記録する。
(3) 商が１になったら終わり、結果を表示する。CindyScriptではTeX形式での数式表現ができますので、指数を使って結果を表示します。
次のスクリプトがその例です。
yakusuu(n):= select(1 .. n, mod(n, #) == 0);
prime(n):= select(1 .. n, length(yakusuu(#)) == 2);
n = 360; // 目的の数
primelist = prime(n);
pn = length(primelist);
kosu = apply(1 .. pn, # = 0); // 割り切れた素数の個数を記録するリスト
m = n; // m は最初がn そのあとは素数で割った商
i = 1;
while(m > 1, // mが1になったら終わる
p = primelist_i; // 素数のリストから一つ取り出して
if(mod(m, p) == 0, // 割り切れたら
m = div(m, p); // mは商にして
kosu_i = kosu_i + 1; // kosu の i 番目に割り切れた数をカウント
,
i = i+1; // 割り切れなかったら次の素数
);
);
disp = ""; // 結果の表示をするTexの式を作る
repeat(pn, i,
if(kosu_i != 0,
if(length(disp) > 0,disp = disp + "・"); // かけ算の印
disp = disp + "$" + text(primelist_i);
if(kosu_i > 1, disp = disp + "^" + text(kosu_i));
disp=disp + "$";
);
);
drawtext([0, 1], n + "=" + disp, size -> 20); // 描画面に結果を表示
最後の行を，
println(n + "=" + disp);
にすれば，コンソールに
360=$2^3$・$3^2$・$5$
と表示されます。