素数判定

2/9を終えたところ. nが341550071728321未満のときまでは素数判定がO(log^3 n)か. Miller-Rabinテストはやっぱり偉いんだなー.
ところで本業の方をschemeで書くという暴挙はどうだろう. LLLとか実装すんの. 大変そう.

;Refs:http://www.ice.nuie.nagoya-u.ac.jp/~h003149b/lang/miller.html
(define (decomposition-k&q n) ;n=q 2^k
  (let loop [(k 0) (q n)]
    (if (even? q)
        (loop (+ k 1) (/ q 2))
      (values k q))))

(define (expt-mod a m n) ;a^m mod n
  (let loop [(a a) (m m) (c 1)]
    (cond [(= m 0) c]
          [(even? m) (loop (remainder (* a a) n) (/ m 2) c)]
          [else (loop a (- m 1) (remainder (* a c) n))])))

(define (Fermat-test n a) ;a^{n-1} mod n ?= 1
  (= 1 (expt-mod a (- n 1) n)))

;; Let n = 1 + 2^k q (t:odd)
;; Check a sequence, a^q, a^{2q}, ..., a^{2^{k-1} q}=a^{(n-1)/4}, a^{2^k q}=a^{(n-1)/2}
(define (Miller-Rabin-test n a) ;See textbooks
  (receive
   (k q)
   (decomposition-k&q (- n 1))
   (let loop [(i 0) (r (expt-mod a q n))]
     (and (< i k)
          (cond [(= r 1) (= i 0)]
                [(= r (- n 1)) #t]
                [else (loop (+ i 1) (remainder (* r r) n))])))))

(define (is-prime? n) ;Thresholds are given by Jaeschke
  (cond ((= n 1) #f)
        ((= n 2) #t)
        ((< n 4759123141)
         (and (Miller-Rabin-test n 2)
              (if (> n 7) (Miller-Rabin-test n 7) #t)
              (if (> n 61) (Miller-Rabin-test n 61) #t)))
        ((< n 341550071728321) ;about 2^{48.2790..}
         (and (Miller-Rabin-test n 2)
              (Miller-Rabin-test n 3)
              (Miller-Rabin-test n 5)
              (Miller-Rabin-test n 7)
              (Miller-Rabin-test n 11)
              (Miller-Rabin-test n 13)
              (Miller-Rabin-test n 17)))))