; 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