素数鉛筆

京都大学生協限定素数価格577円で販売している素数ものさしは皆様よくご存知のことと思う。それでは素数えんぴつ（A Prime Pencil）はどうだろう。このあたりで販売しているようだ。本日twitterでWolfram Japanさんが宣伝していたのでさっそく試してみた。いや，購入したのではなくて，「Juliaでプログラミング教育をはじめよう（仮題）」の練習問題を作ったのだ。本家のほうはバリバリのMathematica関数プログラミングなので，スマートすぎておじぃさんには手に負えない。

注意点：今日も型の整合性のところで難渋した。適当な演算では自動型変換は行われないので，明示的に型を指定する必要があった。10^mが  Int128 (-2^{127} 〜 2^{127} - 1)を必要とする場面で，10^m → Int128(10)^mとすべきところに嵌まった。また，zeros()とかones() で配列を準備するところで，型指定が必要なのに気付かなかった。p = zeros(Int128,N)とか q = ones(Int,N) などとしなければならない。

結果は以下の通りで出力は省略するが，23桁には３つの素数（ 96686312646216567629137， 57686312646216567629137，95918918997653319693967）があり，最大の24桁では１つの素数（357686312646216567629137）が得られた。これ以上はありません。これを鉛筆に焼いて左から削っていけばよい。そうそう，素数鉛筆に印刷されている切り捨て可能素数の定義を説明するのを忘れていた。与えられた素数の上位桁（または下位桁）から１桁づつとりのぞいた数が再び素数であり，最後に１桁の素数が残るものである。

P. S. @ksasaoさんが右から削る場合にも言及していたので，こちらも加えた。いずれも終了判定はしていないので人力でがんばってほしい。


using Primes

#左から削っていく場合の切り捨て可能素数
function pencil1(N)
  p=zeros(Int128,N)
  q=ones(Int,N)
  r=Int128(10)
  p[1]=3
  p[2]=7
  nb = 0
  nt = 2
  while max(nb,nt) < N-1
    nb=nb+1
    for j in 1:9
      k= (r^q[nb])*j+p[nb]
      if isprime(k)
        nt=nt+1
        p[nt]=k
        q[nt]=q[nb]+1
      end
    end
  end
  println(q)
  println(p)
end

#右から削っていく場合の切り捨て可能素数
function pencil2(N)
  p=zeros(Int128,N)
  q=ones(Int,N)
  r=Int128(10)
  p[1]=2
  p[2]=3
  p[3]=5
  p[4]=7
  nb = 0
  nt = 4
  while max(nb,nt) < N-1
    nb=nb+1
    for j in 1:9
      k= 10*p[nb]+j
      if isprime(k)
        nt=nt+1
        p[nt]=k
        q[nt]=q[nb]+1
      end
    end
  end
  println((q,p))

end