2 posts tagged “gauche”

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 ((zero? exp) 1)
        ((even? exp) (modulo (square (power base (/ exp 2) m)) 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 ((n 2) (h (sqrt x)))
        (cond ((> n h) #t)
              ((zero? (remainder x n)) #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 ((x 0))
   (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 ((n 1) (d 1) (x 1))
    (if (> (- precision) (/ (log d) (log 10)))
      x
      (loop (+ 1 n) (/ d n) (+ x (/ d n))))))

(let loop ((n 0) (e e))
  (if (> n precision)
    #t
    (let ((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" は確定する．

n	e
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))
  (let ((d (/ d n)) (z (+ x d)))
    (cond ((= (floor (* 1000 z)) (floor (* 1000 y)))
           (let ((fz (floor z)))
             (display fz)
             (flush)
             (loop (+ n 1) (* 10 d) (* 10 (- z fz)) 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 が未だに理解できていないのでまた後日．

scheme の処理系 guile, gauche を試してみた．

どちらもインタプリタフロントエンドで readline が効かないようなので何とかする．

~/.guile を作る．

(use-modules (ice-9 readline))
(activate-readline)

Gauche-readline を入れて gosh-rl を起動する．

guile, gauche 両者とも，シンボルが case-sensitive なのは良いと思う．case-insensitive だと

(equal? "a" (symbol->string (string->symbol "a")))

が成り立たないかもしれなくて気持ち悪い．

ところで，guile に何故か nil という自己評価フォームが予め定義されている．

nil => nil
'nil => nil
(eqv? nil 'nil) => #t
(eqv? nil '()) => #f
(eqv? nil #f) => #f

なんだこれ…．何のためにあるんだろう．R4RS には nil という単語は見当たらないし，R5RS にも注釈として

Note: Programmers accustomed to other dialects of Lisp should be aware that Scheme distinguishes both #f and the empty list from the symbol nil.

という一文があるだけで，これでは nil の在り方を定めているようには読めない．てゆーか guile でも

(equal? nil '()) => #f

だし．