[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/lquSkg9s    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Feb 21:
; 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) 4)))
          (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))


Output:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
d = 5, p = 1, q = -1
loop 0 1 2 1 -1 49
loop 1 1 54 3 1 24
loop 1 3 54 7 1 12
loop 1 21 54 47 1 6
loop 1 17 54 73 1 3
loop 38 77 18 89 1 1
loop 96 63 61 62 1 0
#f

d = 13, p = 1, q = -3
loop 0 1 2 1 -3 50
loop 0 1 2 7 9 25
loop 7 7 3 31 81 12
loop 7 19 3 7 27 6
loop 7 34 3 94 36 3
loop 26 28 58 52 9 1
loop 21 70 32 13 81 0
#f


Create a new paste based on this one


Comments: