; perrin pseudoprimes (define (make-matrix rows columns . value) (do ((m (make-vector rows)) (i 0 (+ i 1))) ((= i rows) m) (if (null? value) (vector-set! m i (make-vector columns)) (vector-set! m i (make-vector columns (car value)))))) (define (matrix-rows x) (vector-length x)) (define (matrix-cols x) (vector-length (vector-ref x 0))) (define (matrix-ref m i j) (vector-ref (vector-ref m i) j)) (define (matrix-set! m i j x) (vector-set! (vector-ref m i) j x)) (define-syntax for (syntax-rules () ((for (var first past step) body ...) (let ((ge? (if (< first past) >= <=))) (do ((var first (+ var step))) ((ge? var past)) body ...))) ((for (var first past) body ...) (let* ((f first) (p past) (s (if (< first past) 1 -1))) (for (var f p s) body ...))) ((for (var past) body ...) (let* ((p past)) (for (var 0 p) body ...))))) (define (matrix-multiply a b) (let ((ar (matrix-rows a)) (ac (matrix-cols a)) (br (matrix-rows b)) (bc (matrix-cols b))) (if (not (= ac br)) (error 'matrix-multiply "incompatible matrices") (let ((c (make-matrix ar bc 0))) (for (i ar) (for (j bc) (for (k ac) (matrix-set! c i j (+ (matrix-ref c i j) (* (matrix-ref a i k) (matrix-ref b k j))))))) c)))) (define (matrix-power m n) (cond ((= n 1) m) ((even? n) (matrix-power (matrix-multiply m m) (/ n 2))) (else (matrix-multiply m (matrix-power (matrix-multiply m m) (/ (- n 1) 2)))))) (define perrin-matrix '#(#(0 1 0) #(0 0 1) #(1 1 0))) (define perrin-init '#(#(3) #(0) #(2))) (define (perrin n) (if (= n 0) 3 (if (= n 1) 0 (if (= n 2) 2 (vector-ref (vector-ref (matrix-multiply (matrix-power perrin-matrix (- n 2)) perrin-init) 2) 0))))) (display (perrin 200)) (newline) (define (matrix-multiply-mod a b m) (let ((ar (matrix-rows a)) (ac (matrix-cols a)) (br (matrix-rows b)) (bc (matrix-cols b))) (if (not (= ac br)) (error 'matrix-multiply-mod "incompatible matrices") (let ((c (make-matrix ar bc 0))) (for (i ar) (for (j bc) (for (k ac) (matrix-set! c i j (modulo (+ (matrix-ref c i j) (* (matrix-ref a i k) (matrix-ref b k j))) m))))) c)))) (define (matrix-power-mod b e m) (cond ((= e 1) b) ((even? e) (matrix-power-mod (matrix-multiply-mod b b m) (/ e 2) m)) (else (matrix-multiply-mod b (matrix-power-mod (matrix-multiply-mod b b m) (/ (- e 1) 2) m) m)))) (define (perrin-mod n) (if (= n 0) 3 (if (= n 1) 0 (if (= n 2) 2 (vector-ref (vector-ref (matrix-multiply-mod (matrix-power-mod perrin-matrix (- n 2) n) perrin-init n) 2) 0))))) (display (perrin-mod 200)) (newline) (display (modulo (perrin-mod 200) 200)) (newline) (define (perrin-prime? n) (zero? (perrin-mod n))) (display (perrin-prime? 271441)) (newline) (define prime? (let ((seed 3141592654)) (lambda (n) (define (rand) (set! seed (modulo (+ (* 69069 seed) 1234567) 4294967296)) (+ (quotient (* seed (- n 2)) 4294967296) 2)) (define (expm b e m) (define (times x y) (modulo (* x y) m)) (let loop ((b b) (e e) (r 1)) (if (zero? e) r (loop (times b b) (quotient e 2) (if (odd? e) (times b r) r))))) (define (spsp? n a) (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1))) ((odd? d) (let ((t (expm a d n))) (if (or (= t 1) (= t (- n 1))) #t (do ((s (- s 1) (- s 1)) (t (expm t 2 n) (expm t 2 n))) ((or (zero? s) (= t (- n 1))) (positive? s)))))))) (if (not (integer? n)) (error 'prime? "must be integer") (if (< n 2) #f (do ((a (rand) (rand)) (k 25 (- k 1))) ((or (zero? k) (not (spsp? n a))) (zero? k)))))))) (display (prime? 271441)) (newline) (define (perrin-pseudo-primes n) (do ((i 2 (+ i 1))) ((< n i)) (when (and (perrin-prime? i) (not (prime? i))) (display i) (newline)))) ; about ten minutes -- linear ; > (perrin-pseudo-primes 1000000) ; 271441 ; 904631 (define (perrin-pseudo-primes n) (let loop ((p-2 3) (p-1 0) (p 2) (i 3)) (when (< i n) (let ((p+1 (+ p-2 p-1))) (when (and (zero? (modulo p+1 i)) (not (prime? i))) (display i) (newline)) (loop p-1 p p+1 (+ i 1)))))) ; about two minutes -- quadratic ; > (perrin-pseudo-primes 1000000) ; 271441 ; 904631