[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on May 27:
; mersenne primes

(define (range . args)
  (case (length args)
    ((1) (range 0 (car args) (if (negative? (car args)) -1 1)))
    ((2) (range (car args) (cadr args) (if (< (car args) (cadr args)) 1 -1)))
    ((3) (let ((le? (if (negative? (caddr args)) >= <=)))
           (let loop ((x(car args)) (xs '()))
             (if (le? (cadr args) x)
                 (reverse xs)
                 (loop (+ x (caddr args)) (cons x xs))))))
    (else (error 'range "unrecognized arguments"))))

(define (filter pred? xs)
  (let loop ((xs xs) (ys '()))
    (cond ((null? xs) (reverse ys))
          ((pred? (car xs))
            (loop (cdr xs) (cons (car xs) ys)))
          (else (loop (cdr xs) ys)))))

(define (prime? n)
  (letrec (
    (expm (lambda (b e m)
      (let ((times (lambda (x y) (modulo (* x y) m))))
        (cond ((zero? e) 1) ((even? e) (expm (times b b) (/ e 2) m))
        (else (times b (expm (times b b) (/ (- e 1) 2) m)))))))
    (digits (lambda (n)
      (let loop ((n n) (ds '()))
        (if (zero? n) ds (loop (quotient n 2) (cons (modulo n 2) ds))))))
    (isqrt (lambda (n)
      (let loop ((x n) (y (quotient (+ n 1) 2)))
        (if (<= 0 (- y x) 1) x (loop y (quotient (+ y (quotient n y)) 2))))))
    (square? (lambda (n) (let ((n2 (isqrt n))) (= n (* n2 n2)))))
    (jacobi (lambda (a n)
      (let loop ((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) (* (loop 2 n) (loop (/ a 2) n)))
              ((< n a) (loop (modulo a n) n))
              ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (loop n a)))
              (else (loop n a))))))
    (inverse (lambda (x m)
      (let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w m))
        (if (zero? w) (modulo a m)
          (let ((q (quotient g w)))
            (loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w))))))))
    (strong-pseudo-prime? (lambda (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))))))))))
    (chain (lambda (m f g u v)
      (let loop ((ms (digits m)) (u u) (v v))
        (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)))))))
    (lucas-pseudo-prime? (lambda (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 (jacobi 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)))))))))))
    (if (not (integer? n)) (error 'prime? "must be integer")
      (if (< n 2) #f (if (even? n) (= n 2) (if (zero? (modulo n 3)) (= n 3)
        (and (strong-pseudo-prime? n 2)
             (strong-pseudo-prime? n 3)
             (lucas-pseudo-prime? n))))))))

(define (mersenne-prime? p)
  (let ((m (- (expt 2 p) 1)))
    (if (not (prime? m)) #f
      (let loop ((i 3) (s 4))
        (if (< p i) (zero? s)
          (loop (+ i 1)
            (modulo (- (* s s) 2) m)))))))

(display (cons 2 (filter mersenne-prime? (range 256))))


Output:
1
(2 3 5 7 13 17 19 31 61 89 107 127)


Create a new paste based on this one


Comments: