[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Mar 4:
; combined n +/- 1 primality prover

; primes n -- list of primes not greater than n in ascending order
(define (primes n) ; assumes n is an integer greater than one
  (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
    (let loop ((i 0) (p 3) (ps (list 2))) ; sieve of eratosthenes
      (cond ((< n (* p p))
              (do ((i i (+ i 1)) (p p (+ p 2))
                   (ps ps (if (vector-ref bits i) (cons p ps) ps)))
                  ((= i len) (reverse ps))))
            ((vector-ref bits i)
              (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
                  ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
                (vector-set! bits j #f)))
            (else (loop (+ i 1) (+ p 2) ps))))))

; prime? n -- #f if provably composite, else #t if probably prime
(define prime? ; strong pseudoprime to prime bases less than 100
  (let ((ps (primes 100))) ; assumes n an integer greater than one
    (lambda (n)
      (define (expm b e m)
        (let loop ((b b) (e e) (x 1))
          (if (zero? e) x
            (loop (modulo (* b b) m) (quotient e 2)
                  (if (odd? e) (modulo (* b x) m) x)))))
      (define (spsp? n a) ; #t if n is a strong pseudoprime base a
        (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
            ((odd? d) (if (= (expm a d n) 1) #t
              (do ((r 0 (+ r 1)))
                  ((or (= (expm a (* d (expt 2 r)) n) (- n 1))
                       (= r s))
                    (< r s)))))))
      (if (< n 2) #f (if (member n ps) #t
        (do ((ps ps (cdr ps)))
            ((or (null? ps) (not (spsp? n (car ps))))
              (null? ps))))))))

; factors n -- list of prime factors of n in ascending order
(define (factors n) ; assumes n is an integer, may be negative
  (if (<= -1 n 1) (list n) ; pollard's rho algorithm
    (if (< n 0) (cons -1 (factors (- n)))
      (let fact ((n n) (c 1) (fs (list)))
        (define (f x) (modulo (+ (* x x) c) n))
        (if (even? n) (fact (/ n 2) c (cons 2 fs))
          (if (= n 1) fs
            (if (prime? n) (sort < (cons n fs))
              (let loop ((t 2) (h 2) (d 1))
                (cond ((= d 1) (let ((t (f t)) (h (f (f h))))
                                 (loop t h (gcd (- t h) n))))
                      ((= d n) (fact n (+ c 1) fs)) ; cyclic
                      ((prime? d)
                        (fact (/ n d) (+ c 1) (cons d fs)))
                      (else (fact n (+ c 1) fs)))))))))))

(define next-prime ; uses 2,3,5,7-wheel
  (let ((gap (vector 1 10 9 8 7 6 5 4 3 2 1 2 1 4
      3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4 3 2
      1 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 6 5 4 3 2
      1 2 1 6 5 4 3 2 1 4 3 2 1 2 1 6 5 4 3 2 1 4
      3 2 1 6 5 4 3 2 1 8 7 6 5 4 3 2 1 4 3 2 1 2
      1 4 3 2 1 2 1 4 3 2 1 8 7 6 5 4 3 2 1 6 5 4
      3 2 1 4 3 2 1 6 5 4 3 2 1 2 1 4 3 2 1 6 5 4
      3 2 1 2 1 6 5 4 3 2 1 6 5 4 3 2 1 4 3 2 1 2
      1 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4 3 2 1 4 3 2
      1 2 1 4 3 2 1 2 1 10 9 8 7 6 5 4 3 2 1 2)))
    (lambda (n)
      (if (< n 2) 2 (if (< n 3) 3 (if (< n 5) 5
        (let loop ((n (+ n (if (even? n) 1 2))))
          (if (prime? n) n (loop (+ n
            (vector-ref gap (modulo n 210))))))))))))

(define (expm b e m) ; modular exponentiation
  (let loop ((b b) (e e) (x 1))
    (if (zero? e) x
      (loop (modulo (* b b) m) (quotient e 2)
            (if (odd? e) (modulo (* b x) m) x)))))


(define (fermat-test? n ps)
  (let loop ((ps ps) (a 1))
    (if (null? ps) #t
      (if (and (= (expm a (- n 1) n) 1)
               (= (gcd (expt a (/ (- n 1) (car ps))) n) 1))
          (loop (cdr ps) 1)
          (loop ps (+ a 1))))))

(define (jacobi a m) ; jacobi symbol
  (let loop1 ((a (modulo a m)) (m m) (t 1))
    (if (zero? a) (if (= m 1) t 0)
      (let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))
        (let loop2 ((a a) (t t))
          (if (even? a) (loop2 (/ a 2) (* t z))
            (loop1 (modulo m a) a
                   (if (and (= (modulo a 4) 3)
                            (= (modulo m 4) 3))
                       (- t) t))))))))

(define (selfridge n) ; initialize lucas sequence
  (let loop ((d-abs 5) (sign 1))
    (let ((d (* d-abs sign)))
      (cond ((< 1 (gcd d n)) (values d 0 0))
            ((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4)))
            (else (loop (+ d-abs 2) (- sign)))))))

(define (lucas p q m n) ; lucas sequences u[n] and v[n] and q^n (mod m)
  (define (even e o) (if (even? n) e o))
  (define (mod n) (if (zero? m) n (modulo n m)))
  (let ((d (- (* p p) (* 4 q))))
    (let loop ((un 1) (vn p) (qn q) (n (quotient n 2))
               (u (even 0 1)) (v (even 2 p)) (k (even 1 q)))
      (if (zero? n) (values u v k)
        (let ((u2 (mod (* un vn))) (v2 (mod (- (* vn vn) (* 2 qn))))
              (q2 (mod (* qn qn))) (n2 (quotient n 2)))
          (if (even? n) (loop u2 v2 q2 n2 u v k)
            (let* ((uu (+ (* u v2) (* u2 v)))
                   (vv (+ (* v v2) (* d u u2)))
                   (uu (if (and (positive? m) (odd? uu)) (+ uu m) uu))
                   (vv (if (and (positive? m) (odd? vv)) (+ vv m) vv))
                   (uu (mod (/ uu 2))) (vv (mod (/ vv 2))))
              (loop u2 v2 q2 n2 uu vv (* k q2)))))))))

(define (lucas-test? n ps)
  (call-with-values
    (lambda () (selfridge n))
    (lambda (d pp qq)
      (let loop ((ps ps) (p pp) (q qq))
        (if (null? ps) #t
          (call-with-values
            (lambda () (lucas p q n (+ n 1)))
            (lambda (un+1 vn+1 qn+1)
              (call-with-values
                (lambda () (lucas p q n (/ (+ n 1) (car ps))))
                (lambda (un+1/p vn+1/p qn+1/p)
                  (if (and (zero? un+1) (= (gcd un+1/p n) 1))
                      (loop (cdr ps) pp qq)
                      (loop ps (+ p 2) (+ p q 1))))))))))))

(define-syntax while
  (syntax-rules ()
    ((while pred? body ...)
      (do () ((not pred?)) body ...))))

(define (bls-prime? n . args)
  (define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs))))
  (define (cycle xs) (set-cdr! (last-pair xs) xs) xs)
  (define wheel (list 2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8
    4 2 4 2 4 8 6 4 6 2 4 6 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10))
  (define (rho n c b)
    (define (f y) (modulo (+ (* y y) c) n))
    (define (g p x y) (modulo (* p (abs (- x y))) n))
    (let stage1 ((x 2) (y (+ 4 c)) (z (+ 4 c)) (j 1) (q 2) (p 1))
      (if (= j b) #f ; timeout
        (if (= x y) (rho n (+ c 1) (- b j)) ; cycle
          (if (= j q) (let ((t (f y))) (stage1 y (f y) z (+ j 1) (* q 2) (g p y t)))
            (if (positive? (modulo j 100)) (stage1 x (f y) z (+ j 1) q (g p x y))
              (let ((d (gcd p n)))
                (if (= d 1) (stage1 x (f y) y (+ j 1) q (g p x y))
                  (if (and (< 1 d n) (bls-prime? d)) d ; factor
                    (let stage2 ((k 1) (z (f z)))
                      (if (= k 100) (rho n (+ c 1) (- b j))
                        (let ((d (gcd (- z x) n)))
                          (if (= d 1) (stage2 (+ k 1) (f z))
                            (if (= d n) (rho n (+ c 1) (- b j))
                              (if (bls-prime? d) d ; factor
                                (rho n (+ c 1) (- b j)))))))))))))))))
  (define (remove-twos n)
    (let loop ((f 1) (r n))
      (if (odd? r) (values f (list 2) r)
        (loop (* f 2) (/ r 2)))))
  (define (join x xs) (if (member x xs) xs (cons x xs)))
  (define (enough? n p f1 f2)
    (< n (* (max (+ (* p f1) 1) (- (* p f2) 1)) (+ (* p p f1 f2 1/2) 1))))
  (if (not (prime? n)) #f ; sanity check
    (let-values (((b1) (if (pair? args) (car args) 20000000))
                 ((b2) (if (and (pair? args) (pair? (cdr args)))
                           (cadr args) 10000000))
                 ((f1 f1s r1) (remove-twos (- n 1)))
                 ((f2 f2s r2) (remove-twos (+ n 1))))
      (let loop ((p 3) (ws (cons 2 (cons 2 (cons 4 (cycle wheel))))))
        (cond ((and (< p b1) (not (enough? n p f1 f2)))
                (case (modulo (+ n 1) p)
                ((0) (while (zero? (modulo r2 p))
                  (set! f2 (* f2 p)) (set! f2s (join p f2s)) (set! r2 (/ r2 p))))
                ((2) (while (zero? (modulo r1 p))
                  (set! f1 (* f1 p)) (set! f1s (join p f1s)) (set! r1 (/ r1 p)))))
                (when (< r1 (* p p))
                  (set! f1 (- n 1)) (set! f1s (join r1 f1s)) (set! r1 1))
                (when (< r2 (* p p))
                  (set! f2 (+ n 1)) (set! f2s (join r2 f2s)) (set! r2 1))
                (loop (+ p (car ws)) (cdr ws)))
        (else (set! p (min p b1))
              (when verbose?
                (display (list 'wheel n p f1 f1s r1 f2 f2s r2)) (newline))
              (let loop ((more1? #t) (more2? #t) (more3? #t))
                (cond ((and (not (enough? n p f1 f2)) (or more1? more2?))
                        (if more1?
                            (let ((f (rho r1 1 b2)))
                              (if (not f) (loop #f #t #t)
                                (begin (set! f1 (* f1 f)) (set! f1s (join f f1s))
                                       (set! r1 (/ r1 f)) (loop #t #t#t))))
                            (let ((f (rho r2 1 b2)))
                              (if (not f) (loop #f #f #t)
                                (begin (set! f2 (* f2 f)) (set! f2s (join f f2s))
                                       (set! r2 (/ r2 f)) (loop #f #t #t))))))
                      (more3?
                        (when verbose?
                          (display (list 'rho n p f1 f1s r1 f2 f2s r2)) (newline))
                        (loop #f #f #f))
                      ((not (enough? n p f1 f2)) #f) ; failure to prove
                      ((< r1 f1)
                        (when verbose? (display "fermat only") (newline))
                        (fermat-test? n f1s)) ; condition 1 only
                      ((< r2 f2)
                        (when verbose? (display "lucas only") (newline))
                        (lucas-test? n f2s)) ; condition 2 only
                      (else
                        (when verbose? (display "fermat and lucas") (newline))
                        (and (fermat-test? n (cons r1 f1s)) ; 1 and 3
                             (lucas-test? n (cons r2 f2s)))))))))))) ; 2&4


(define rand #f)
(define randint #f)
(let ((two31 #x80000000) (a (make-vector 56 -1)) (fptr #f))
  (define (mod-diff x y) (modulo (- x y) two31)) ; generic version
  ; (define (mod-diff x y) (logand (- x y) #x7FFFFFFF)) ; fast version
  (define (flip-cycle)
    (do ((ii 1 (+ ii 1)) (jj 32 (+ jj 1))) ((< 55 jj))
      (vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
    (do ((ii 25 (+ ii 1)) (jj 1 (+ jj 1))) ((< 55 ii))
      (vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
    (set! fptr 54) (vector-ref a 55))
  (define (init-rand seed)
    (let* ((seed (mod-diff seed 0)) (prev seed) (next 1))
      (vector-set! a 55 prev)
      (do ((i 21 (modulo (+ i 21) 55))) ((zero? i))
        (vector-set! a i next) (set! next (mod-diff prev next))
        (set! seed (+ (quotient seed 2) (if (odd? seed) #x40000000 0)))
        (set! next (mod-diff next seed)) (set! prev (vector-ref a i)))
      (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle)))
  (define (next-rand)
    (if (negative? (vector-ref a fptr)) (flip-cycle)
      (let ((next (vector-ref a fptr))) (set! fptr (- fptr 1)) next)))
  (define (unif-rand m)
    (let ((t (- two31 (modulo two31 m))))
      (let loop ((r (next-rand)))
        (if (<= t r) (loop (next-rand)) (modulo r m)))))
  (init-rand 19380110) ; happy birthday donald e knuth
  (set! rand (lambda seed
    (cond ((null? seed) (/ (next-rand) two31))
          ((eq? (car seed) 'get) (cons fptr (vector->list a)))
          ((eq? (car seed) 'set) (set! fptr (caadr seed))
                                 (set! a (list->vector (cdadr seed))))
          (else (/ (init-rand (modulo (numerator
                  (inexact->exact (car seed))) two31)) two31)))))
  (set! randint (lambda args
    (cond ((null? (cdr args))
            (if (< (car args) two31) (unif-rand (car args))
              (floor (* (next-rand) (car args)))))
          ((< (car args) (cadr args))
            (let ((span (- (cadr args) (car args))))
              (+ (car args)
                 (if (< span two31) (unif-rand span)
                   (floor (* (next-rand) span))))))
          (else (let ((span (- (car args) (cadr args))))
                  (- (car args)
                     (if (< span two31) (unif-rand span)
                       (floor (* (next-rand) span))))))))))

(define (rand-prime n . args)
  (let ((b (if (pair? args) (car args) 10)))
    (let loop ((r (randint 1 b)) (n n))
      (if (= n 1) (next-prime r)
        (loop (+ (* r b) (randint b)) (- n 1))))))

(define verbose? #f)

(let ((n (rand-prime 42)))
  (display "Proving primality of ")
  (display n) (display ": ")
  (display (bls-prime? n)) (newline))


Output:
1
Proving primality of 710429394483947895803394585080463740040049: #t


Create a new paste based on this one


Comments: