Namareba食べたい

備忘録てきなもの。だらだら書いていきます。

G00913

http://ksnctf.sweetduet.info/problem/20

 

問題文に

FLAG_Q20_{first 10-digit prime found in consecutive digits of π}

とある

要するにπの中に現れる最初の10桁の素数を見つけよ

 

とのことらしい。

 

10桁かぁエラトステネスは計算量てきに無理だろうな

フェルマールールか?でも実装出来る自信が・・・てかそれ以前に・・・

 

順番に試すのもアリダナ

と、思いつつとりあえずggると答えが・・・

もうこれでいったれ!

http://answers.yahoo.com/question/index?qid=20090308170117AA3AdKK

 

すこし問題はちがうけど良い記事があった↓

http://kosugi.typepad.com/blog/2007/01/%E3%81%84%E3%81%BE%E3%81%94%E3%82%8D-first-10-digit-prime-found-in-consecutive-digits-of-e-com.html

 

引用ここから

いまごろ { first 10-digit prime found in consecutive digits of e }.com

1/21 の NHK スペシャルを見ていたら,ハイウェイ沿いに掲げられた Google の謎の求人広告の話題が挙がっていた.

{ first 10-digit prime found in consecutive digits of e }.com

あれー,確か昔に聞いたことがあったけど,当時はハナから「どうせ解けないだろうなー」と諦めていたような気がする.しかし,いま改めて向き合ってみると,e が分かっていれば大したことないと分かる.

10 桁というと 32bit で収まりきらないのでチト気を付ける必要があるが,gauche多倍長整数を扱えるのでそちらに任せる.素数判定は,まぁその都度行っても問題ないだろう.ダメだったらそれはそれで.

年末は SICP を読んでいたので,素数判定はフェルマーテストで行ってみる.合格した場合は念のためシラミ潰す.

(define (square x) (* x x))

(define (power base exp m)
  (cond *1 m))
        (else        (modulo (* base (power base (- exp 1) m)) m))))

(define (fermat-test x)
  (if (= 2 x) #t
      (= 2 (power 2 x x))))

(define (prime? x)
  (if (fermat-test x)
    (if (< x 2)
      #f
      (let loop *2 #f)
              (else (loop (if (< 2 n) (+ 2 n) 3) h)))))
    #f))

末尾再帰の if がどうも美しくないような気がするけどソコは触れない約束で.

で,問題の e だけど,ググれば何万桁オーダーで載っているサイトも見つかる.これをテキストファイルにして頭から読み込みながら計算してゆけばいずれ答えは出る.やや bruteforcely だが,他に (数学的にエレガントな?) 解法は思いつかない.

(define low (expt 10 9))
(define high (expt 10 10))

(display
 (let loop *3
   (if (and (>= x low) (prime? x))
       x
     (loop (remainder (+ (- (char->integer (read-char)) 48) (* 10 x)) high)))))

とりあえずコレを q.scm という名前で保存して

cat e.txt | gosh q.scm

と実行すると 7427466391 が求められる.小数点以下 99 桁目だ.意外と浅い場所にあったが,ググってみると "Google" というキーワードと共に沢山のページが引っかかるので,おそらくコレで正解なのだろう (http://7427466391.com 自体はもう存在しないようだ).

しかし,これでいいのか? 既知の e を用いて求めるのはアリなのか!?

e 自体は lim(n → ∞)(1 + 1/n)^n で求められる (ググった).しかし n が大きくなると結果が 2.71 にも届かないうちに NaN ってしまうので無理がある.ここはひとつ,マクローリン展開された Σ(n → ∞)(1/n!) を使ってジワジワ攻めてゆくのが良いはず (ググった).というかマクローリン展開なつかしす(w

gauche の多倍長有理数を使えば,メモリの許す限りの精度で e を算出できる.

(define precision 120)

(define e
  (let loop *4
    (if (> (- precision) (/ (log d) (log 10)))
      x
      (loop (+ 1 n) (/ d n) (+ x (/ d n))))))

(let loop *5
  (if (> n precision)
    #t
    (let *6
  (let *7 (z (+ x d)))
    (cond *8 (floor (* 1000 y)))
           (let *9 x)))
          (else (loop (+ n 1) d z x)))))

1,000 を掛けているのは,実は適当.頭 4 桁くらい一致してればまぁ一桁目は確定でいいかな~なんて.あはは.2 桁目までしか見ないと 960 桁目で実際の値 (e.txt) と食い違ったので,余裕を持ってみた次第 (こんなんでいいのか?

それはさて置き,上記ソースを e.scm という名前で保存して

gosh e.scm | gosh q.scm

とすると,見事 7427466391 が求められる.いやっほう,エレガント!! (微妙).

本当は,e.scm の方をコルーチンにして「呼ばれるたびに次の e の桁を返す関数」みたいにすると格好いいんだろうけど,call/cc が未だに理解できていないのでまた後日.

 

引用ここまで 

 

休日にでもしっかり理解し直したいなぁ

 

*1:zero? exp) 1)
        ((even? exp) (modulo (square (power base (/ exp 2) m

*2:n 2) (h (sqrt x)))
        (cond ((> n h) #t)
              ((zero? (remainder x n

*3:x 0

*4:n 1) (d 1) (x 1

*5:n 0) (e e

*6:x (floor e)))
      (display x)
      (loop (+ 1 n) (* 10 (- e x))))))

さて,これで求められた結果を先ほどの q.scm に食わせて…って,おい,何か違うだろう.e.txt を作る手間が増えただけじゃないか.そうじゃなくて,もっと,こう,なんつーか,エレガンスが必要じゃないか? これだと予め決められた桁数の範囲内でしか試行できない.できれば,「何桁目で答えが出るのか分からない状態で試行を続ける」というアプローチでいきたい.

n = 0 から出発して,前回よりも非常に小さい値 (階乗の逆数) を加算し続けてゆくので e はどんどん収束してゆくはず.例えば,n = 3 において 0.16 を加算しても 1 の位への桁上がりは起きないので一桁目の "2" は確定する.また,n = 5 において 0.04 を加算しても 0.1 の位への桁上がりは起きないので二桁目の "7" は確定する.

ne
  0  1.00000000000000
 1  2.00000000000000
 2  2.50000000000000
 3  2.66666666666667
 4  2.70833333333333
 5  2.71666666666667
 6  2.71805555555556
 7  2.71825396825397
 8  2.71827876984127
  9  2.71828152557319
10  2.71828180114638
11  2.71828182619849
12  2.71828182828617

 

こうやって,確定した値を次々と取り出してゆき,10 桁に達したら素数判定を行い,10 桁目を捨てる.この作業を素数が見つかるまで延々と繰り返す.コードはあまり綺麗ではないが,こんな感じ.

(let loop ((n 1) (d 1) (x 0) (y 0

*7:d (/ d n

*8:= (floor (* 1000 z

*9:fz (floor z)))
             (display fz)
             (flush)
             (loop (+ n 1) (* 10 d) (* 10 (- z fz