[ create a new paste ] login | about

Link: http://codepad.org/J9oc7Q3B    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Feb 5:
; n+1 primality prover

(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))))))))))

(define (selfridge n)
  (let loop ((d-abs 5) (sign 1))
    (let ((d (* d-abs sign)))
      (cond ((< 1 (gcd d n)) (values d 0 0))
            ((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4)))
            (else (loop (+ d-abs 2) (- sign)))))))

(define (lucas p q m n) ; right-to-left
  (define (even e o) (if (even? n) e o))
  (define (mod n) (if (zero? m) n (modulo n m)))
  (let ((d (- (* p p) (* 4 q))))
    (let loop ((un 1) (vn p) (qn q) (n (quotient n 2))
               (u (even 0 1)) (v (even 2 p)) (k (even 1 q)))
      (if (zero? n) (values u v k)
        (let ((u2 (mod (* un vn))) (v2 (mod (- (* vn vn) (* 2 qn))))
              (q2 (mod (* qn qn))) (n2 (quotient n 2)))
          (if (even? n) (loop u2 v2 q2 n2 u v k)
            (let* ((uu (+ (* u v2) (* u2 v)))
                   (vv (+ (* v v2) (* d u u2)))
                   (uu (if (and (positive? m) (odd? uu)) (+ uu m) uu))
                   (vv (if (and (positive? m) (odd? vv)) (+ vv m) vv))
                   (uu (mod (/ uu 2))) (vv (mod (/ vv 2))))
              (loop u2 v2 q2 n2 uu vv (* k q2)))))))))

(define (u p q m n) (let-values (((u v k) (lucas p q m n))) u))

(define (facts n)
  (let loop ((f 2) (fs (list)) (fp 1) (r n))
    (cond ((< n (* fp fp)) (reverse fs))
          ((< r (* f f)) (reverse (cons r fs)))
          ((zero? (modulo r f))
            (loop f (if (member f fs) fs (cons f fs))
                  (* fp f) (/ r f)))
          (else (loop (+ f 1) fs fp r)))))

(define (nplus1-prime? n)
  (let-values (((d p q) (selfridge n)))
    (if (< 1 (gcd d n)) #f
      (if (positive? (u p q n (+ n 1))) #f
        (let f-loop ((fs (facts (+ n 1))))
          (if (null? fs) #t
            (let d-loop ((p p) (q q))
              (if (= (gcd (u p q n (/ (+ n 1) (car fs))) n) 1)
                  (f-loop (cdr fs))
                  (d-loop (+ p 2) (+ p q 1))))))))))

(display (nplus1-prime? (+ #e1e12 39))) (newline)
(display (nplus1-prime? (+ #e1e12 61))) (newline)


Output:
1
2
#t
#t


Create a new paste based on this one


Comments: