[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Sep 25:
(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 (iroot k n)
  (let loop ((u n))
    (let* ((s u)
           (t (+ (* (- k 1) s)
                 (quotient n (expt s (- k 1)))))
           (u (quotient t k)))
      (if (< u s) (loop u) s))))

(define (digits n . args)
  (let ((b (if (null? args) 10 (car args))))
    (let loop ((n n) (d '()))
      (if (zero? n) d
          (loop (quotient n b)
                (cons (modulo n b) d))))))

(define square?
  (let ((q11 (make-vector 11 #f))
        (q63 (make-vector 63 #f))
        (q64 (make-vector 64 #f))
        (q65 (make-vector 65 #f)))
    (do ((k 0 (+ k 1))) ((< 5 k))
      (vector-set! q11 (modulo (* k k) 11) #t))
    (do ((k 0 (+ k 1))) ((< 31 k))
      (vector-set! q63 (modulo (* k k) 63) #t))
    (do ((k 0 (+ k 1))) ((< 31 k))
      (vector-set! q64 (modulo (* k k) 64) #t))
    (do ((k 0 (+ k 1))) ((< 32 k))
      (vector-set! q65 (modulo (* k k) 65) #t))
    (lambda (n)
      (if (not (integer? n)) (error 'square? "must be integer")
        (if (< n 1) #f
          (if (not (vector-ref q64 (modulo n 64))) #f
            (let ((r (modulo n 45045)))
              (if (not (vector-ref q63 (modulo r 63))) #f
                (if (not (vector-ref q65 (modulo r 65))) #f
                  (if (not (vector-ref q11 (modulo r 11))) #f
                    (let ((q (iroot 2 n)))
                      (if (= (* q q) n) q #f))))))))))))

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

(define (pseudoprime? n)
  (and (spsp? n 2) (spsp? n 3) (spsp? n 5)))

(define (td n limit)
  (let twos ((n n) (fs (list)))
    (cond ((= n 1) fs)
          ((even? n) (twos (/ n 2) (cons 2 fs)))
          (else (let odds ((n n) (f 3) (fs fs))
                  (cond ((< limit f) fs)
                        ((< n (* f f)) (reverse (cons n fs)))
                        ((zero? (modulo n f))
                          (odds (/ n f) f (cons f fs)))
                        (else (odds n (+ f 2) fs))))))))

(define (rho n limit)
  (let loop ((t 2) (h 2) (d 1) (c 1) (limit limit))
    (define (f x) (modulo (+ (* x x) c) n))
    (cond ((zero? limit) #f)
          ((= d 1) (let ((t (f t)) (h (f (f h))))
                     (loop t h (gcd (- t h) n) c (- limit 1))))
          ((= d n) (loop 2 2 1 (+ c 1) (- limit 1)))
          ((pocklington? d) d)
          (else (rho d (- limit 1))))))

(define (pocklington? n)
  (if (not (pseudoprime? n)) #f
    (if (< n #e1e12) (= (car (td n #e1e6)) n)
      (let* ((n-1 (- n 1))
             (root2 (iroot 2 n)) (root3 (iroot 3 n))
             (fs (td n-1 1000)) (f (apply * fs)))
        (let f-loop ((f f) (fs fs))
          (if (< f root3)
              (let ((p (rho (/ n-1 f) #e1e6)))
                (if p (f-loop (* f p) (cons p fs)) 'unknown))
              (let a-loop ((a 2))
                (if (not (= (expm a n-1 n) 1)) #f
                  (let loop ((fs fs))
                    (let ((g (gcd (- (expm a (/ n-1 (car fs)) n) 1) n)))
                      (cond ((< 1 g n) #f)
                            ((= g n) (a-loop (+ a 1)))
                            ((pair? (cdr fs)) (loop (cdr fs)))
                            ((< root2 f) #t)
                            (else (let* ((ds (digits n f))
                                         (c2 (car ds)) (c1 (cadr ds))
                                         (s (- (* c1 c1) (* 4 c2))))
                                    (not (square? s)))))))))))))))

(define (mersenne n) (- (expt 2 n) 1))

(display (pocklington? (mersenne 89))) (newline)
(display (pocklington? (mersenne 521))) (newline)
(display (pseudoprime? 3825123056546413051)) (newline)
(display (pocklington? 3825123056546413051)) (newline)


Output:
1
2
3
4
#t
#t
#t
#f


Create a new paste based on this one


Comments: