[ create a new paste ] login | about

Link: http://codepad.org/54jgdwKo    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Sep 30:
; aks primality prover

(define (split n xs)
  (let loop ((n n) (xs xs) (zs '()))
    (if (or (zero? n) (null? xs))
        (values (reverse zs) xs)
        (loop (- n 1) (cdr xs) (cons (car xs) zs)))))

(define (make-list n x)
  (let loop ((n n) (xs '()))
    (if (zero? n) xs
      (loop (- n 1) (cons x xs)))))

(define (square x) (* x x))

(define (log2 n) (/ (log n) (log 2)))

(define (expm b e m)
  (define (m* x y) (modulo (* x y) m))
  (cond ((zero? e) 1)
        ((even? e) (expm (m* b b) (/ e 2) m))
        (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))

(define (ilog b n)
  (let loop1 ((lo 0) (b^lo 1) (hi 1) (b^hi b))
    (if (< b^hi n) (loop1 hi b^hi (* hi 2) (* b^hi b^hi))
      (let loop2 ((lo lo) (b^lo b^lo) (hi hi) (b^hi b^hi))
        (if (<= (- hi lo) 1) (if (= b^hi n) hi lo)
          (let* ((mid (quotient (+ lo hi) 2))
                 (b^mid (* b^lo (expt b (- mid lo)))))
            (cond ((< n b^mid) (loop2 lo b^lo mid b^mid))
                  ((< b^mid n) (loop2 mid b^mid hi b^hi))
                  (else mid))))))))

(define (td-prime? n)
  (if (even? n) (= n 2)
    (let loop ((d 3))
      (cond ((< n (* d d)) #t)
            ((zero? (modulo n d)) #f)
            (else (loop (+ d 2)))))))

(define (uniq-factors n)
  (define (uniq-cons x xs)
    (if (null? xs) (list x)
      (if (= x (car xs)) xs
        (cons x xs))))
  (let twos ((n n) (fs '()))
    (if (even? n)
        (twos (/ n 2) (uniq-cons 2 fs))
        (if (= n 1) fs
          (let odds ((n n) (d 3) (fs fs))
            (cond ((< n (* d d))
                    (reverse (uniq-cons n fs)))
                  ((zero? (modulo n d))
                    (odds (/ n d) d (uniq-cons d fs)))
                  (else (odds n (+ d 2) fs))))))))

(define (prime-power? n)
  (if (even? n) (if (= (expt 2 (ilog 2 n)) n) 2 #f)
    (let loop ((a 2) (n n))
      (let* ((b (expm a n n)) (p (gcd (- b a) n)))
        (cond ((or (= p 1) (= a b)) #f)
              ((prime? p) (if (= (expt p (ilog p n)) n) p #f))
              (else (loop (+ a 1) p)))))))

; The order of n modulo r,
; written as ord_r(n),
; is the smallest number k
; such that n^k = 1 (mod r).

(define (ord r n)
  (do ((k 2 (+ k 1)))
      ((= (expm n k r) 1) k)))

; Find the smallest r such
; that ord_r(n) > 4(log2 n)^2.

(define (compute-r n)
  (let ((target (* 4 (square (log2 n)))))
    (do ((r 3 (+ r 1)))
        ((and (td-prime? r)
              (not (= r n))
              (< target (ord r n)))

(define (phi n)
  (let loop ((fs (uniq-factors n)) (t n))
    (if (null? fs) t
      (loop (cdr fs) (* t (- 1 (/ (car fs))))))))

(define (poly-mult-mod xs ys r n)
  (define (times x) (lambda (y) (* x y)))
  (define (plus xs ys)
    (let loop ((xs xs) (ys ys) (zs (list)))
      (cond ((null? xs) (reverse (append (reverse ys) zs)))
            ((null? ys) (reverse (append (reverse xs) zs)))
            (else (loop (cdr xs) (cdr ys)
                        (cons (+ (car xs) (car ys)) zs))))))
  (define (mod-poly xs)
    (let-values (((hs ts) (split r (reverse xs))))
      (reverse (plus hs ts))))
  (define (mod-n x) (modulo x n))
  (let ((xs (reverse xs)))
    (let loop ((xs (cdr xs)) (zs (map (times (car xs)) ys)))
      (if (null? xs) (map mod-n (mod-poly zs))
        (loop (cdr xs) (plus (cons 0 zs)
                             (map (times (car xs)) ys)))))))

(define (poly-power-mod bs e r n)
  (let loop ((bs bs) (e e) (rs (list 1)))
    (if (zero? e) rs
      (loop (poly-mult-mod bs bs r n)
            (quotient e 2)
            (if (even? e) rs
              (poly-mult-mod rs bs r n))))))

(define (binomial-test? a r n)
  (not (equal? (poly-power-mod (list 1 a) n r n)
               (append (list 1)
                       (make-list (- r 1) 0)
                       (list a)))))

(define (aks-prime? n)
  (if (prime-power? n) #f
    (let* ((r (compute-r n))
           (phi-r (phi r))
           (sqrt-phi-r (sqrt phi-r))
           (log2-n (log2 n))
           (sqrt-phi-r-log2-n (* sqrt-phi-r log2-n)))
      (let loop ((i 2))
        (if (< i sqrt-phi-r-log2-n)
            (if (and (not (= n i)) (zero? (modulo n i))) #f
              (loop (+ i 1)))
            (let loop ((a 1))
              (if (< a sqrt-phi-r-log2-n)
                  (if (binomial-test? a r n) (loop (+ a 1)) #f)

(display (compute-r 89)) (newline)
(display (aks-prime? 89)) (newline)


Create a new paste based on this one