Project Euler Problem 27

問題27: 二次式 n^2 + an + n = p (|a| < 1000, |b| < 1000, p = 素数) において、n に 0から連続する整数を代入したときもっとも大きな n まで式が成り立つ a b の組を求め、その積を答えよ。[=>問題文]
 効率的な素数判定方法を考えていたら2日経ってしまいました。最終的に、処理時間は 1.89秒になりました。これは速いなのか遅いほうなのか、これから他の人の解答を見るのが楽しみです。

問題27の解答

;; Emacs Lisp
(require 'cl)
(require 'math-lib)

(defun problem027 (limit)
  (let* ((b-list (primes-sieve limit))
         (prime-p (make-prime-number-p b-list))
         (pairs (let ((lis) (b))
                  (dolist (b (reverse b-list) lis)
                    (do ((a (1- limit) (1- a)) (end (* -1 limit)))
                        ((= a end))
                      (push (list a b) lis))))))
    (do ((n 1 (1+ n)))
        ((>= 1 (length pairs)) (apply #'* (car pairs)))
      (let ((lis) (n2 (* n n)) (continue t))
        (while (and (consp pairs) continue)
          (let ((pair (pop pairs)))
            (destructuring-bind (a b) pair
              (if (>= n b)
                  (setq continue nil)
                (let ((x (+ n2 (* a n) b)))
                  (if (and (> x 1) (funcall prime-p x))
                      (push pair lis)))))))
        (setq pairs lis)))))

(problem027 1000)

 n^2+an+b=p に n=0 を代入すると b=p となり、bは素数で無ければならないことがわかります。また、n=b を代入すると左辺は b^2+ab+b となりこれは明らかに b で割り切れます。したがって 連続する nにおいて、式が成り立つ n の最大値は高々 b - 1 であることがわかります。
 problem027 では、はじめに a, b の組み合わせをすべて作成 (pairs) しています。但しbは1000未満の素数だけです。メインのループでは n の値を増やしながら、成り立たないものを除去ししてゆき リスト pairs の要素が 1 つになったところで、処理を完了としています。
 なお素数判定関数 prime-p は「素数判定関数作成関数」によって取得しています。

素数判定関数作成関数

(defun make-prime-number-p (primes)
  (lexical-let*
      ((primes (if (> 2 (length primes)) [2 3] (vconcat primes)))
       (capacity (length primes))
       (index (1- capacity))
       (expand-capacity
        (lambda ()
          (setq primes (vconcat primes (make-vector capacity 0)))
          (setq capacity (* 2 capacity))))
       (iterate)
       (iterate-2
        (lambda ()
          (when (= index (1- capacity))
            (funcall expand-capacity))
          (let ((n (aref primes index)) (r 0))
            (while (zerop r)
              (incf n 2)
              (setq r 1)
              (do ((i 0 (1+ i)) (p 0) (end (truncate (sqrt n))))
                  ((or (>= p end) (zerop r)))
                (setq p (aref primes i))
                (setq r (mod n p))))
            (aset primes (incf index) n))))
       (iterate-1
        (lambda ()
          (if (< index (1- capacity))
              (aref primes (incf index))
            (progn
              (setq iterate iterate-2)
              (funcall iterate-2)))))
       (binary-search
        (lambda (x)
          (let ((i 0) (j index) (exists))
            (while (and (<= i j) (not exists))
              (let* ((k (/ (+ i j) 2)) (m (aref primes k)))
                (cond ((< m x) (setq i (1+ k)))
                      ((> m x) (setq j (1- k)))
                      (t (setq exists t)))))
            exists)))
       (prime-number-p
        (lambda (x)
          (let ((last-prime (aref primes index)))
            (cond
             ((< x last-prime) (funcall binary-search x))
             ((> x last-prime)
              (do ((p (funcall iterate) (funcall iterate)))
                  ((>= p x) (= p x))))
             (t t))))))

    (setq iterate iterate-1)
    (lambda (x) (if (evenp x) 
                    nil
                  (funcall prime-number-p x)))))

 最大値のわからない素数判定に使える関数が欲しくて考えた末できたのがこれです。はじめに 素数を格納した vector primes を作成し、これに収まる整数 x についてはバイナリ探索で該当する素数があるかを判定します。primes よりもおおきな素数が渡されたときは、primes を拡張して対応できるようにしています。一度拡張するとそれよりさらに大きい数がくるまではバイナリ探索ですむので、大量のランダムな数を素数判定するには有利なんじゃないかと。(でも条件によるか・・・小さな素因数を持つ合成数の判定はむしろ遅くなるな、すごくダメな気がしてきた。普通に試し割りの方式に修正するかも)。
 無駄に複雑に見えるかもしれませんが、実はもう少し発展させて便利なものをつくろうと画策しています。うまくいくかな。