[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Jan 24:
; primality checking, revisited

(define (expm b e m)
  (define (m* x y) (modulo (* x y) m))
  (cond ((zero? e) 1)
        ((even? e) (expm (m* b b) (/ e 2) m))
        (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))

(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 (isqrt n)
  (let loop ((x n) (y (quotient (+ n 1) 2)))
    (if (<= 0 (- y x) 1) x
      (loop y (quotient (+ y (quotient n y)) 2)))))

(define (square? n)
  (let ((n2 (isqrt n)))
    (= n (* n2 n2))))

(define (jacobi a n)
  (if (not (and (integer? a) (integer? n) (positive? n) (odd? n)))
      (error 'jacobi "modulus must be positive odd integer")
      (let jacobi ((a a) (n n))
        (cond ((= a 0) 0)
              ((= a 1) 1)
              ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1)))
              ((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n)))
              ((< n a) (jacobi (modulo a n) n))
              ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a)))
              (else (jacobi n a))))))

(define legendre jacobi)

(define (inverse x n)
  (let loop ((x (modulo x n)) (a 1))
    (cond ((zero? x) (error 'inverse "division by zero"))
          ((= x 1) a)
          (else (let ((q (- (quotient n x))))
                  (loop (+ n (* q x)) (modulo (* q a) n)))))))

(define (miller? n a)
  (let loop ((r 0) (s (- n 1)))
    (if (even? s) (loop (+ r 1) (/ s 2))
      (if (= (expm a s n) 1) #t
        (let loop ((r r) (s s))
          (cond ((zero? r) #f)
                ((= (expm a s n) (- n 1)) #t)
                (else (loop (- r 1) (* s 2)))))))))

(define (chain m f g x0 x1)
  (let loop ((ms (digits m 2)) (u x0) (v x1))
    (cond ((null? ms) (values u v))
          ((zero? (car ms)) (loop (cdr ms) (f u) (g u v)))
          (else (loop (cdr ms) (g u v) (f v))))))

(define (lucas? n)
  (let loop ((a 11) (b 7))
    (let ((d (- (* a a) (* 4 b))))
      (cond ((square? d) (loop (+ a 2) (+ b 1)))
            ((not (= (gcd n (* 2 a b d)) 1))
              (loop (+ a 2) (+ b 2)))
            (else (let* ((x1 (modulo (- (* a a (inverse b n)) 2) n))
                         (m (quotient (- n (legendre d n)) 2))
                         (f (lambda (u) (modulo (- (* u u) 2) n)))
                         (g (lambda (u v) (modulo (- (* u v) x1) n))))
                    (let-values (((xm xm1) (chain m f g 2 x1)))
                      (zero? (modulo (- (* x1 xm) (* 2 xm1)) n)))))))))

(define (prime? n)
  (cond ((or (not (integer? n)) (< n 2))
          (error 'prime? "must be integer greater than one"))
        ((even? n) (= n 2)) ((zero? (modulo n 3)) (= n 3))
        (else (and (miller? n 2) (miller? n 3) (lucas? n)))))

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

(do ((k 2 (+ k 1))) ((< 1000 k))
  (if (prime? (mersenne k))
      (begin (display k)
             (display #\tab)
             (display (mersenne k))
             (newline))))


Output:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
2	3
3	7
5	31
7	127
13	8191
17	131071
19	524287
31	2147483647
61	2305843009213693951
89	618970019642690137449562111
107	162259276829213363391578010288127
127	170141183460469231731687303715884105727
521	6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151
607	531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127

Timeout


Create a new paste based on this one


Comments: