; lucas pseudoprimality tester

; based on trn.c from http://www.trnicely.net/misc/bpsw.html

; jacobi a m -- jacobi symbol
(define (jacobi a m)
  (if (not (integer? a)) (error 'jacobi "must be integer")
    (if (not (and (integer? m) (positive? m) (odd? m)))
        (error 'jacobi "modulus must be odd positive integer")
        (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))))))))))

; standard lucas primality test -- assume n odd integer > 2, not a square
(define (lucas-prime? n)
  (let loop ((d 5))
    (if (or (< 1 (gcd d n)) (< -1 (jacobi d n)))
        (loop (* (+ (abs d) 2) (if (positive? d) -1 1)))
        (let ((p 1) (q (/ (- 1 d) 2)))
          (display "d = ") (display d) (display ", p = ")
          (display p) (display ", q = ") (display q) (newline)
          (let loop ((u 0) (u2 1) (v 2) (v2 p) (q q) (bits (/ (+ n 1) 2)))
            (display "loop ") (display u) (display " ") (display u2)
            (display " ") (display v) (display " ") (display v2)
            (display " ") (display q) (display " ") (display bits) (newline)
            (if (zero? bits) (zero? u) ; #t if u==0 else #f
              (let ((u2 (modulo (* u2 v2) n))
                    (v2 (modulo (- (* v2 v2) (+ q q)) n))
                    (q (modulo (* q q) n)))
                (if (odd? bits)
                    (let* ((u (+ (* u2 v) (* u v2)))
                           (u (/ (if (odd? u) (+ u n) u) 2))
                           (u (modulo u n))
                           (v (+ (* v2 v) (* d u2 u)))
                           (v (/ (if (odd? v) (+ v n) v) 2))
                           (v (modulo v n)))
                      (loop u u2 v v2 q (quotient bits 2)))
                    (loop u u2 v v2 q (quotient bits 2))))))))))

; two tests: 97 is prime, 99 is composite, but both are reported composite
(display (lucas-prime? 97)) (newline) (newline)
(display (lucas-prime? 99))

d = 5, p = 1, q = -2
loop 0 1 2 1 -2 49
loop 1 1 56 5 4 24
loop 1 5 56 17 16 12
loop 1 85 56 63 62 6
loop 1 20 56 62 61 3
loop 12 76 87 36 35 1
loop 78 20 1 62 61 0
#f

d = 13, p = 1, q = -6
loop 0 1 2 1 -6 50
loop 0 1 2 13 36 25
loop 13 13 57 97 9 12
loop 13 73 57 85 81 6
loop 13 67 57 34 27 3
loop 14 1 16 13 36 1
loop 90 13 65 97 9 0
#f