Project Euler: Problem 7

クリスマスも年末年始も…素数を数えるんだ。

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6^(th) prime is 13.

What is the 10001^(st) prime number?

http://projecteuler.net/index.php?section=problems&id=7

10001番目の素数を求めよとのこと。

例によって Gauche で書いてみた。

#!/usr/bin/env gosh

(define (prime-fermat? n)
  (case n
    [(1) #f]
    [(2) n]
    [else
      (if (even? n)
        #f
        (= (remainder (expt 2 (- n 1)) n) 1))])) ; (gcd n 2) ==> 1 for all odd number except 1

(define (prime-inexact? n) ; fast but inexact
  (case n
    [(1) #f]
    [(2) n]
    [else
      (if (even? n)
        #f
        (let loop [(lis '(3 5 7 11 13 17 19 23 29 31))] ; magic of first 10 prime
          (if (pair? lis)
            (if (= n (car lis))
              n
              (if (= (remainder n (car lis)) 0)
                #f
                (loop (cdr lis))))
            n)))]))

(define (prime-exact? n) ; slow but exact
  (case n
    [(1) #f]
    [(2) n]
    [else
      (if (even? n)
        #f
        (let loop [(i 3) ; minumum odd prime
                   (m (floor->exact (sqrt n)))]
          (if (<= i m)
            (if (= (remainder n i) 0)
              #f
              (loop (+ i 2) m))
            n)))]))

(define (prime? n) ; try faster test first
  (and (prime-inexact? n)
       (prime-fermat?  n)
       (prime-exact?   n)))

(define (main args)
  (define plis '())

  (let loop [(n   3)     ; try tests from 3 and other odd numbers
             (l   1)     ; (length lis) ==> 1
             (lis '(2))] ; calculated prime before n
    (if (< l 10001)
      (if (prime? n)
        (loop (+ n 2) (+ l 1) (cons n lis))
        (loop (+ n 2) l lis))
      (set! plis lis)))
  (print (car plis))
  0)

これまで使っていた素数判定の関数をちょっくら作り直してみた。いい加減、フェルマー法じゃなくて Miller-Rabin 等に書き換えようかとも思ったけどむしろ遅くなりそうな気もするし、当面は現状維持。超適当な判定方法 *1 から順番に詳細な判定をするようにしたことでかなりの高速化をすることができた。

とりあえず、最後は完全に割ってみて素数かどうか判定してる *2 けど、求める素数の範囲がある程度まで絞れるのならばカーマイケル数を埋め込みで弾いてもよかったかも知れない。

*1:prime-inexact?: 小さい素数で適当に割ってみて合成数かどうか判定する

*2:prime-exact?