Zenn
🙃

dc(1) をつかって素数を生成する

2021/03/21に公開
shell
dc
素数
tech

前回の記事以降 dc にハマっている蜜柑です． 今回は dc と shell を使って素数 (らしきもの) を生成しました．

フェルマの小定理

フェルマの小定理とは，素数についての性質の一つを表した定理のひとつです．

この定理は，素数 p と p の倍数ではない適当な整数 a について次式が成り立つというものです．

a^{p-1} \equiv 1 \hspace{-1ex} \mod p

この定理の対偶を利用したものにフェルマテストがあります． これは，n と互いに素な a が次式を満たすなら n が合成数であることを示すものです．

a^{n-1} \not\equiv 1 \hspace{-1ex} \mod n

方針

素数 (らしき数) を次の方法で生成します:

  1. 適当に思いついた数を n，既知の素数を a としてフェルマテストを行います．
  2. フェルマテストの結果 n は合成数でないと判断された場合，a を別の素数に置き換えて再び 1. に戻ります．
  3. 合成数だと判断されたら n を思いつき直します．
  4. 飽きるまで， 2. と 3. を行います．
  5. どうやら n は素数かもしれません．

dc の使いみち

dc は無限精度の計算機を謳っているだけあって，大きな桁の数も難なく計算してくれます． 更に GNU 版の dc は冪剰余の演算子を備えているので，最初に示した式を簡単に計算できます[1]．

例えば，n=32529，a=11 としてフェルマテストを実行するなら次のようにします． 結果は 1 ではないので，今回の n は素数ではないようです．

$ echo "11 32529d1-r|p" | dc
22108

適当に数を思いつく

適当に数を思いつくために /dev/urandom を利用します． 適当さは適当でいいので /dev/random を使うべきではありません．

/dev/urandom から読み取れるものはバイナリデータですが，バイナリデータをそのまま扱うことは苦労が多いのでテキストに変換します． od(1) はどこにでもあるので，これを使います． od を次のように実行すると，ランダムな 16 進数を 64 バイト分出力してくれます． いろいろなオプションがついていますが，それぞれ，オフセットを表示しないこと (-An)，64 バイトだけ表示すること (-N64)，8 バイト区切りの 16 進数で表示すること (-tx8) を指定しています．

$ od -An -N64 -tx8 /dev/urandom
 6cbff7a8b328a334 cd02607bfca7f24a
 87222740f2b9496d 6e3ace39a3cb82df
 7a6342f2bb52161b 0f468f29f9f3ef56
 b31f49269323dc29 b9d63f2ae7ddb6df

しかし，これでは 64 バイトの適当な数というよりも 8 バイトの適当な数が 8 つという感じです． そこで sed(1) を利用して空白文字を取り去ります． 後々 dc に渡すことも考え，大文字に変換することも同時に行ってしまいます． 次に示すコマンドは GNU 版の sed では動きますが，それ以外は -z オプションで怒られるかもしれません．

$ od -An -N64 -tx8 /dev/urandom | \
> sed -ze 's/\(.*\)/\U\1/g; s/[[:space:]]//g'
2F841CA1352FACBCEE4481E502AE1E8286EC562496741EE9B56ECB0085D092733B65E2F685E120F97603CA247A556E39350F64CC55D1B4420EE09FC1CDDB6CA3

これで今度こそ適当な数字を思いつくことができました．

既知の素数を用意する

このセクションはもはやありません。 ここにあったコマンドを実行しても当時の目的は果たされません。(2021-04-25)

かつてここにあったもの

フェルマテストの a として利用する適当な素数を用意する必要があります． 怪しげな素数ですが，ここに過去の自分が用意した素数の一覧があるのでそれを使います． prime.txt という名前でまとめて保存しておきます．

$ for i in `seq 0 9`; do
> curl https://www.kusaremkn.com/prime/$i.txt >>prime.txt;
> done

とりあえず書いてみる

そして次のように適当なスクリプトを書きます．

prime.sh
#!/bin/sh

# フェルマの小定理
# fermat(a, n)
fermat() {
	if [ $# -ne 2 ]; then
		return 1
	fi
	if [ `echo "[1sr]sP0sr16i$1 $2d1-r|1=Plrp" | dc` -eq 0 ]; then
		return 1
	fi
	return 0
}

# 素数 (っぽい) かどうか調べる
# isPrime(n)
isPrime() {
	if [ $# -ne 1 ]; then 
		return 1
	fi
	cat prime.txt | while read line; do
		printf "."
		if ! fermat $line $1; then
			printf "#"
			return 1
		fi
	done
}

# 適当な数を思いつく
# randGen()
randGen() {
	od -An -N64 -tx8 /dev/urandom | \
		sed -ze 's/\(.*\)/\U\1/g; s/[[:space:]]//g'
}

RandNum=`randGen`
while ! isPrime $RandNum; do
	RandNum=`randGen`
done
echo
echo "`echo "16i${RandNum}p" | dc` may be prime number."

実行してみる

randGen() の od のオプション -N64 を適当に変更して，次のコマンドを実行した結果の概要です．

$ time sh prime.sh
オプション 生成された素数 実行時間 (real)
-N4 (32bit) 2732855357 (素数) 1m26.050s
-N8 (64bit) 1309708607\
123597119 (素数)		 1m50.535s
-N16 (128bit) 1915996383\
2120558763\
1910472604\
200883991 (素数)		 3m42.332s
-N32 (256bit) 1103595478\
4997404420\
6585173347\
4602803964\
1301281974\
4438750812\
9675453800\
04961237 (素数)		 21m45.643s
-N64 (512bit) ¯\_(ツ)_/¯ NaN

ああ，あまりにもおそすぎます． 見たところ，素数ではない数は a=2 のフェルマテストで脱落し，かかっている時間の殆どは素数らしき数が本当に素数がどうかを判断するところにありそうです． そこで，a を 2, 3, 5, 7, 11, 13, 17, 19 くらいに絞って実験し直します． ホビー用途の素数なので擬素数が混ざっていても大して問題にはなりません．

prime2.sh
--- prime.sh	2021-03-21 01:46:19.726158589 +0900
+++ prime2.sh	2021-03-21 01:55:52.887732224 +0900
@@ -18,7 +18,7 @@
 	if [ $# -ne 1 ]; then 
 		return 1
 	fi
-	cat prime.txt | while read line; do
+	for line in 2 3 5 7 11 13 17 19; do
 		printf "."
 		if ! fermat $line $1; then
 			printf "#"

改訂版での実行結果は次のとおりです:

オプション 生成された素数 実行時間 (real)
-N4 (32bit) 3085616531 (素数) 0m0.021s
-N8 (64bit) 1476409253\
5285870949 (素数)		 0m0.033s
-N16 (128bit) 2285957700\
8916006975\
9998898705\
825808489 (素数)		 0m1.313s
-N32 (256bit) 4108798460\
1097619771\
1779692954\
2846983472\
4855398896\
5209709038\
9879828337\
8274057 (素数)		 0m0.570s
-N64 (512bit) 9533664173\
9419813679\
0984503783\
6154684221\
2583977960\
1243603755\
1969554340\
2952713495\
8674724972\
5200953556\
5464412058\
3392109147\
8883257305\
1010090813\
9778012640\
5727 (素数)		 0m14.179s
-N128 (1024bit) 3254941820\
9324001573\
1807475136\
9735125857\
8869541575\
2810515522\
7253914761\
3578272604\
2088571921\
9090978348\
3237423812\
2296951841\
7373903126\
2364559528\
3487164146\
5215004822\
6956222569\
3448123410\
2302911507\
3570729925\
4351092010\
4290664385\
3210692795\
0686049633\
4686527330\
6387476759\
6649858689\
3196569838\
4595940637\
2816237742\
75613597 (判別不能)		 2m5.427s

まとめ

フェルマの小定理を利用することで，50 行に満たない小さなプログラムから 150 桁を超える素数と 300 桁程度の素数らしき数を生成することができました． このスクリプトは稀に擬素数を吐くようですが，手元では擬素数を得られませんでした．

素数判定には Wolfram Alpha を利用しましたが，150 桁以上の数を因数分解できることには驚きです． なお，最後の 300 桁超の数は素数かどうか判別できていません． 因数分解してくださる猛者がいらっしゃれば，どうぞその時間を他の有益なことに使ってください．

追記 (03-21)

この記事のタイトル dc(1) をつかって素数を生成する は果たして妥当なのでしょうか…… dc は素数の判定をしてるだけでは……?

そんなことはさておき，2048 bit の素数らしき数ができました． 4 つ同時に走らせて 13m11.442s でした． 乱数の機嫌が良ければ現実的な時間で巨大な素数を生成できそうです．

131822638207803451818222022043343534930349405962170460137775369668243\
918198952377075315111533957942709442817646728644891254720889768060014\
831953685980296248989791817421998664142932806083697183594643459671264\
103538309336726400009142483274239518814459042898661126353825986566197\
871625123634052835406923728126868778315373679588807579690636642899050\
757525570047426707170818378819340026823598914837248888758194882296768\
022436684305342911252486968812434482586534738474186602525743258649195\
024587461692524389968889416687705118513245048530860428736113964395229\
87273359220494812974966157894261748260319895792642965522293861009
脚注

  1. GNU 版ではない dc の場合は Wikipedia の dc の この節 にあるコードを読んでください． できます． ↩︎

Discussion