素数全書 問題1.17

ちょこっと計算にもっともふさわしい問題である。ただし、(10^4)^2=100,000,000=一億以下の素数を求めておく必要がある。いままで使っていた既知の素数で割ってみるテストを繰り返すという schemeプログラムだと時間がかかりすぎてどうにもならなかったので、エラトステネスの篩を作るプログラムを作成した。リストで作るとどうやらメモリが足りなくなるので、ここはvectorを使った。

(define (eratosthenes n)
  (let ((ls (make-vector (+ n 1) 1)))
    (vector-set! ls 0 n)
    (vector-set! ls 1 0)
    (letrec ((check-list (lambda (x lst)
                           (let loop ((k (* x 2)))
                             (cond
                               ((> k n) lst)
                               (else (begin
                                       (vector-set! lst k 0)
                                       (loop (+ k x)))))))))
      (check-list 2 ls)
      (let loop ((k 3))
      (cond
        ((> k n) ls)
        ((= (vector-ref ls k) 0) (loop (+ k 2)))
        (else (begin
                (check-list k ls)
                (loop (+ k 2)))))))))
      

(define (count-vect v)
  (let ((n (vector-ref v 0)))
    (let loop ((k 0)(result 0))
      (cond
        ((> k n) result)
        ((= (vector-ref v k) 1) (loop (+ k 1) (+ result 1)))
        (else (loop (+ k 1) result))))))


(define (test n lst)
  (let loop ((k 1) (result 0))
    (let ((f (abs (+ (* k k) k -1354363))))
      (cond
        ((> k n) result)
        ((= (vector-ref lst f) 1) (loop (+ k 1) (+ result 1)))
        (else (loop (+ k 1) result))))))


(define (test2 n lst)
  (let loop ((k 1) (result 0) (l '()) )
    (let ((f (abs (+ (* k k) k -1354363))))
      (cond
        ((> k n) result)
        ((memq f l) (loop (+ k 1) result l))
        ((= (vector-ref lst f) 1) (loop (+ k 1) (+ result 1) (append (list f) l)))
        
        (else (loop (+ k 1) result l))))))

> (define primes (eratosthenes 100000000))
> (count-vect primes)
5761455
> (test 10000 primes)
5356
> (test2 10000 primes)
5356


一億以下の素数の数は、5761455個であり、1〜10000のxについて問題の多項式の(絶対)値が素数であるのは、5356個である。またtest2の結果より、これらには重複はない。

よって、53.56%の素数生成率である。