[Project Euler] Problem 27 「二次式素数」

オイラーは以下の二次式を考案している:

n2 + n + 41.

この式は, n を0から39までの連続する整数としたときに40個の素数を生成する. しかし, n = 40 のとき 40² + 40 + 41 = 40(40 + 1) + 41 となり41で割り切れる. また, n = 41 のときは 41² + 41 + 41 であり明らかに41で割り切れる.

計算機を用いて, 二次式 n² - 79n + 1601 という式が発見できた. これは n = 0 から 79 の連続する整数で80個の素数を生成する. 係数の積は, -79 × 1601 で -126479である.

さて, |a| < 1000, |b| < 1000 として以下の二次式を考える (ここで |a| は絶対値): 例えば |11| = 11, |-4| = 4である。

n² + an + b

n = 0 から始めて連続する整数で素数を生成したときに最長の長さとなる上の二次式の, 係数 a, b の積を答えよ.

素数の表があれば解けそうです.

1. 素数の表を作る.
2. 係数a, bから二次式を作る手続きを用意する.
3. 二次式から素数のリストを作る手続きを用意する.
4. aとbを与えられた範囲で変化させ, 素数のリストを作る.
5. 素数のリストが一番長いものを求める.

手続きは次のようになります.

(require srfi/1)

;; エラトステネスのふるい
(define (sieve pv n)
  (define (loop i)
    (when (< i (vector-length pv))
      (vector-set! pv i #f)
      (loop (+ i n))))
  (loop (+ n n)))

(define (find-prime pv start)
  (cond ((<= (vector-length pv) start) #f)
        ((vector-ref pv start) start)
        (else (find-prime pv (+ start 1)))))

(define (eratosthenes pv)
  (define (loop n)
    (when (and n (< n (sqrt (vector-length pv))))
      (sieve pv n)
      (loop (find-prime pv (+ n 1)))))
  (loop 2))
  
(define (prime-vector n)
  (let ((pv (make-vector (+ n 1) #t)))
    (eratosthenes pv)
    (vector-set! pv 0 #f)
    (vector-set! pv 1 #f)
    pv))

(define (prime-vect-to-list pv)
  (define (loop ans n)
    (if (< n (vector-length pv))
        (let ((p (find-prime pv n)))
          (if p 
              (loop (cons p ans) (+ p 1)) 
              (reverse ans)))
        (reverse ans)))
  (loop () 2))

;; 1000万以下の素数表を作る
(define prime-vect (prime-vector 10000000))

(define (prime? n) (vector-ref prime-vect (abs n)))

(define (quadric a b)
  (lambda (n) (+ (* n n)
                 (* a n)
                 b)))

(define (make-prime-list a b)
  (define (loop rslt n)
    (let ((v ((quadric a b) n)))
      (if (prime? v)
          (loop (cons v rslt) (+ n 1))
          (list a b (reverse rslt)))))
  (loop () 0))

(define ans
  (apply append
         (map (lambda (a)
                (map (lambda (b) (make-prime-list a b))
                     (iota 1999 -999)))
              (iota 1999 -999))))
 
(define sorted 
  (sort ans 
        (lambda (x y) (> (length (list-ref x 2)) 
                         (length (list-ref y 2))))))

計算してみます. 少し時間がかかります.

ようこそ DrRacket, バージョン 5.3.3 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (make-prime-list 1 41)
(1
 41
 (41 43 47 53 61 71 83 97 113 131 151 173 197 223 251 281
  313 347 383 421 461 503 547 593 641 691 743 797 853 911
  971 1033 1097 1163 1231 1301 1373 1447 1523 1601))
> (car sorted)
(-61
 971
 (971 911 853 797 743 691 641 593 547 503 461 421 383 347
  313 281 251 223 197 173 151 131 113 97 83 71 61 53 47 43
  41 41 43 47 53 61 71 83 97 113 131 151 173 197 223 251 281
  313 347 383 421 461 503 547 593 641 691 743 797 853 911 971
  1033 1097 1163 1231 1301 1373 1447 1523 1601))
> (* -61 971)
-59231
>