; 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