素数全書練習問題1.5 - Finite Groups Fun 有限群あるいはちょこっと計算ファンの日記

n番目の素数を与える Gandhi formula を証明せよという問題。（GahndiとprimeでWeb検索するとどうしてもprime minister Gahndiが出てくるというオチが...）

ヒントを参考に、 $1/(2^d-1)=\sum_{n=1}^\infty 1/2^{nd}$ から、d桁毎に1が立っている二進小数展開ができ、それが合成数のところでうまく打ち消しあうように $\mu(d)$ を掛けてあるのでエラトステネスの篩ができるのだと気づいたが、きちんと証明しようとするとなかなかうまくいかずかなりの間悩んでいた。

和の部分を変形して,二進小数展開での篩の形にはできるのである。（-1/2はN=1,d=1の例外の項を打ち消すため。）

$\sum_{N=1}^\infty \quad \sum_{d|P_{n-1}! \quad & \quad d|N} \mu(d)/2^N - 1/2$

一週間ぐらい悩んで、やっと、Nと $P_{n-1}!$ が素なとき

$\sum_{d|P_{n-1}! \quad & \quad d|N} \mu(d)=1$

素でないとき

$\sum_{d|P_{n-1}! \quad & \quad d|N} \mu(d)=0$

であることに気付いた。２つ目の式が０になるのは、 $K=GCD(P_{n-1}!,\quad N)$ としたとき、 $\sum_{d|K} \mu(d)=0$ という、メビウス関数の性質からである。

$P_{n-1}!$ と素ということは、 $P_{n-1}$ 以下の素因子を持たないと同値。

つまり、問題の和は小数点以下N桁目が、Nが $P_1,P_2...,P_{n-1}$ で割れるとき0、そうでない時 1 になっているという数である。よって初めて１が立つのは $P_n$ のところであり、 $-log_2$ をとると $P_n-1$ 以上 $P_n$ 未満の数となる。これに１を加えて小数を切り捨てると $P_n$ となって証明完了である。

ちなみに $P_n$ がでるときにできた篩は、 $P_n$ 以上の素数の部分は１になっており、合成数でも $P_n$ 以上の素因子のみを含むとそこで１になっている。つまり、この篩が初めて間違うのは $(P_n)^2$ である。

とまあ、理論的にはこの通りなのだが、本当にそうかと計算したくなって、schemeでプログラムを作ってみた。手習い段階のschemeプログラムなので恥ずかしいのだが、とりあえず公開してみる。

　schemeでは整数演算は有効桁に制限がなく、分数も整数比で形式的に計算されるため、計算精度は考慮しなくていい点が整数論とすごく相性がいいと思う。（最近だとMathematicaを使うのが普通なのかも？）

(define (binary-disp x n)
  (display "0.")
  (let loop((x x) (count 1))
    (cond
      ((> count n) (newline))
      (else
       (let ((x2 (* x 2)))
         (if (>= x2 1) (display "1") (display "0")) 
         (loop (- x2 (truncate x2)) (+ count 1)))))))

(define (next-prime plist n)
  (let ((last-p (list-ref plist (sub1 n))))
    (let loop ((k (+ last-p 2)) (check-list plist))
      (cond
        ((null? check-list) k)
        (else
         (let ((sqrt-k (sqrt k)) (p (car check-list)))
           (cond
             ((> p sqrt-k) k)
             ((zero? (modulo k p)) (loop (+ k 2) plist))
             (else (loop k (cdr check-list))))))))))

(define prime
  (let ((plist '(2 3 5)) (listed-n 3))
    (lambda (n)
      (let loop()
        (cond
          ((zero? n) 1)
          ((<= n listed-n) (list-ref plist (sub1 n)))
          (else 
           (let ((k (next-prime plist listed-n)))
             (set! plist (append plist (list k)))
             (set! listed-n (add1 listed-n))
             (loop))))))))

(define (mu n)
  (let loop ((k 1) (result 1))
    (let ((p (prime k)))
      (cond
        ((> p n) result)
        ((zero? (modulo n p))
         (begin
           (cond
             ((zero? (modulo (quotient n p) p)) 0)
             (else (loop (add1 k) (* -1 result))))))
        (else (loop (add1 k) result))))))

(define (factorial n)
  (let loop ((d 1) (result 1))
    (cond
      ((> d n) result)
      (else (loop (add1 d) (* result d))))))

(define (test n)
  (let ((f (factorial (prime (sub1 n)))))
    (let loop ((d 1) (result (- 1/2)))
      (cond
      ((> d f) result)
;      ((> d f) (- 1 (/ (log result) (log 2))))
      ((zero? (modulo f d)) (loop (add1 d) (+ result (* (mu d) (/ 1 (sub1 (expt 2 d)))))))
      (else (loop (add1 d) result))))))

n=5で試してみる。(binary-disp (test 5) 200)を実行すればよい。結果とその解釈はつぎのようになる。

0.000000000010100010100010000010100000100010100010000010000010100000100010100000100010
           ↑↑  ↑↑  ↑    ↑↑    ↑  ↑↑  ↑    ↑    ↑↑    ↑  ↑↑    ↑  ↑
　　　　　11 13  17 19 23    29 31   37  41 43  47   53     59 61  67 71 73    79  83

00001000000010001010001010001000000010000010001000001010001000001010000...
   ↑      ↑  ↑↑  ↑↑  ↑      ↑    ↑     
   89     97 101 103 107 109 113  *121* 127

$P_{4}=7$ までは消えて、 $P_{5}=11$ で1が立っている。以下、エラトステネスの篩が作られているのだが、 $P_{4}!=7!$ と素になる最少の合成数である11*11=121で、予想通りに間違って１が立っていることが確認できる。