labunix's blog

labunixのラボUnix

「primes」コマンドが無いので、安全素数を「openssl/factor」で。

■例のシェルスクリプト大喜利の1問目のお題の安全素数について。
 投稿しても採用されないと思いますので。。。

 シェルスクリプト大喜利 第十二回 ~冬季シェリンピック(なんだその痛い題名は!)
 http://pub.usp-lab.com/shogiri/shogiri12.html

■jpman projectにはあるようだけど、
 「factor」コマンドのバイナリには「primes」の文字列が無い。

 http://linuxjm.sourceforge.jp/html/bsd-games/man6/primes.6.html

$ strings /usr/bin/factor  | grep -A 1 prime 
Print the prime factors of each specified integer NUMBER.  If none
are specified on the command line, read them from standard input.

■「Plan9のツールを入れる」のとは別の話の気がする。。。

$ apt-file search bin/primes
9base: /usr/lib/plan9/bin/primes

$ apt-cache search ^9base\$
9base - Plan 9 userland tools

■「openssl」と「factor」や、「clisp」で遊んだことがあるので、
 安全素数も簡単に出来そう。

 bashで素数×素数の生成(openssl/factor)
 http://d.hatena.ne.jp/labunix/20130227

 clispで素数と素数の数を得る(1億まで)
 http://d.hatena.ne.jp/labunix/20111207

■こうかな。opensslとfactorを仲良く。
 awkの「print」文だとたまに小数を出す。

$ for n in `seq 1 100`;do \
    echo $((0x`openssl prime "$n" | awk '($3!="not"){print $1}'`)) ; \
  done | grep -v "^0\$" | \
  for next in `xargs`;do \
    echo $next | awk '{print ($1-1)/2}' | xargs factor 2>/dev/null | \
    awk '(NF==2){print ($0*2)+1}'; \
  done
5
7
11
23
47
59
83

■awkを「printf」に変えて横並び。

$ for n in `seq 1 100`;do \
    echo $((0x`openssl prime "$n" | awk '($3!="not"){print $1}'`)) ; \
  done | grep -v "^0\$" | \
  for next in `xargs`;do \
    echo $next | awk '{printf "%d",($1-1)/2}' | xargs factor | \
      awk '(NF==2){print ($0*2)+1}'; \
  done | xargs echo;
5 7 11 23 47 59 83

■答え合わせ。合っているようだ。

 http://ja.wikipedia.org/wiki/%E5%AE%89%E5%85%A8%E7%B4%A0%E6%95%B0

■せっかくなので、wikiにある500まで。
 時間も3秒ほど。大してかからない。

$ env LANG=C date; \
  for n in `seq 1 500`;do \
    echo $((0x`openssl prime "$n" | awk '($3!="not"){print $1}'`)) ; \
  done | grep -v "^0\$" | \
  for next in `xargs`;do \
    echo $next | awk '{printf "%d",($1-1)/2}' | xargs factor | \
      awk '(NF==2){print ($0*2)+1}';done | xargs echo; \
  env LANG=C date;
Sun Feb  9 18:47:44 JST 2014
5 7 11 23 47 59 83 107 167 179 227 263 347 359 383 467 479
Sun Feb  9 18:47:47 JST 2014

■以下には2903まであるので、試してみましょう。

 A005385 Safe primes p: (p-1)/2 is also prime. 
 http://oeis.org/A005385

■3000まで約17秒。
 columnだと縦になるのか。。。「2963 2999」はおまけ。

$ env LANG=C date; \
  for n in `seq 1 3000`;do \
    echo $((0x`openssl prime "$n" | awk '($3!="not"){print $1}'`)) ; \
  done | grep -v "^0\$" | \
  for next in `xargs`;do \
    echo $next | awk '{printf "%d",($1-1)/2}' | xargs factor | \
      awk '(NF==2){print ($0*2)+1}';done | column -c 80; \
  env LANG=C date;
Sun Feb  9 18:55:27 JST 2014
5	59	227	467	719	1019	1367	1823	2099	2819
7	83	263	479	839	1187	1439	1907	2207	2879
11	107	347	503	863	1283	1487	2027	2447	2903
23	167	359	563	887	1307	1523	2039	2459	2963
47	179	383	587	983	1319	1619	2063	2579	2999
Sun Feb  9 18:55:44 JST 2014