codepad

Language:

; sophie germain primes

(define (inverse x m)
  (let loop ((x x) (a 0) (b m) (u 1))
    (if (positive? x)
        (loop (modulo b x) u x (- a (* (quotient b x) u)))
        (if (= b 1) (modulo a m) (error 'inverse "must be coprime")))))

(define (drop n xs)
  (let loop ((n n) (xs xs))
    (if (or (zero? n) (null? xs)) xs
      (loop (- n 1) (cdr xs)))))

(define sort #f)
(define merge #f)
(let ()
  (define dosort
    (lambda (pred? ls n)
      (if (= n 1)
          (list (car ls))
          (let ((i (quotient n 2)))
            (domerge pred?
                     (dosort pred? ls i)
                     (dosort pred? (list-tail ls i) (- n i)))))))
  (define domerge
    (lambda (pred? l1 l2)
      (cond
        ((null? l1) l2)
        ((null? l2) l1)
        ((pred? (car l2) (car l1))
         (cons (car l2) (domerge pred? l1 (cdr l2))))
        (else (cons (car l1) (domerge pred? (cdr l1) l2))))))
  (set! sort
    (lambda (pred? l)
      (if (null? l) l (dosort pred? l (length l)))))
  (set! merge
    (lambda (pred? l1 l2)
      (domerge pred? l1 l2))))

(define (isqrt n)
  (if (not (and (positive? n) (integer? n)))
      (error 'isqrt "must be positive integer")
      (let loop ((x n))
        (let ((y (quotient (+ x (quotient n x)) 2)))
          (if (< y x) (loop y) x)))))

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

(define (primes n)
  (let ((sieve (make-vector n #t)))
    (let loop ((p 2) (ps (list)))
      (cond ((= n p) (reverse ps))
            ((vector-ref sieve p)
              (do ((i p (+ i p))) ((<= n i))
                (vector-set! sieve i #f))
              (loop (+ p 1) (cons p ps)))
            (else (loop (+ p 1) ps))))))

(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 (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 (strong-pseudoprime? 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 (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 (chain n u v u2 v2 d q m)
  (let loop ((u u) (v v) (u2 u2) (v2 v2) (q q) (qkd q) (m m))
    (if (zero? m) (values u v qkd)
      (let* ((u2 (modulo (* u2 v2) n))
             (v2 (modulo (- (* v2 v2) (* 2 q)) n))
             (q (modulo (* q q) n)))
        (if (odd? m)
            (let* ((t1 (* u2 v)) (t2 (* u v2))
                   (t3 (* v2 v)) (t4 (* u2 u d))
                   (u (+ t1 t2)) (v (+ t3 t4))
                   (qkd (modulo (* qkd q) n)))
              (loop (modulo (quotient (if (odd? u) (+ u n) u) 2) n)
                    (modulo (quotient (if (odd? v) (+ v n) v) 2) n)
                    u2 v2 q qkd (quotient m 2)))
            (loop u v u2 v2 q qkd (quotient m 2)))))))

(define (powers-of-two n)
  (let loop ((s 0) (n n))
    (if (odd? n) (values s n)
      (loop (+ s 1) (/ n 2)))))

(define (strong-lucas-pseudoprime? n)
  ; assumes odd positive integer not a square
  (call-with-values
    (lambda () (selfridge n))
    (lambda (d p q)
      (if (zero? p) (= n d)
        (call-with-values
          (lambda () (powers-of-two (+ n 1)))
          (lambda (s t)
            (call-with-values
              (lambda () (chain n 1 p 1 p d q (quotient t 2)))
              (lambda (u v qkd)
                (if (or (zero? u) (zero? v)) #t
                  (let loop ((r 1) (v v) (qkd qkd))
                    (if (= r s) #f
                      (let* ((v (modulo (- (* v v) (* 2 qkd)) n))
                             (qkd (modulo (* qkd qkd) n)))
                        (if (zero? v) #t (loop (+ r 1) v qkd))))))))))))))

(define prime?
  (let ((ps (primes 100)))
    (lambda (n)
      (if (not (integer? n)) (error 'prime? "must be integer"))
      (if (or (< n 2) (square? n)) #f
        (let loop ((ps ps))
          (if (pair? ps)
              (if (zero? (modulo n (car ps))) (= n (car ps)) (loop (cdr ps)))
              (and (strong-pseudoprime? n 2)
                   (strong-pseudoprime? n 3)
                   (strong-lucas-pseudoprime? n))))))))

(define (germain b f m)

(define (p c) (- (* c 3003 (expt 10 b)) 1))
  (define (q c) (- (* c 6006 (expt 10 b)) 1))
  (define (r c) (- (* c 15015 (expt 10 (- b 1))) 1))

(define (p-start p) (inverse (* 3003 (expt 10 b)) p))
  (define (q-start p) (inverse (* 6006 (expt 10 b)) p))
  (define (r-start p) (inverse (* 15015 (expt 10 (- b 1))) p))

(let ((sieve (make-vector (+ m 1) #t)) ; ignore sieve[0]
        (ps (drop 6 (primes f)))) ; sieving primes

; sieve on p
    (do ((ps ps (cdr ps))) ((null? ps))
      (do ((i (p-start (car ps)) (+ i (car ps)))) ((< m i))
        (vector-set! sieve i #f)))

; sieve on q
    (do ((ps ps (cdr ps))) ((null? ps))
      (do ((i (q-start (car ps)) (+ i (car ps)))) ((< m i))
        (vector-set! sieve i #f)))

; sieve on r
    (do ((ps ps (cdr ps))) ((null? ps))
      (do ((i (r-start (car ps)) (+ i (car ps)))) ((< m i))
        (vector-set! sieve i #f)))

; collect and test primes
    (let ((gs (list))) ; sophie-germain primes
      (do ((c 1 (+ c 1))) ((< m c) (sort < gs))
        (when (and (vector-ref sieve c) (prime? (p c)))
          (when (prime? (q c)) (set! gs (cons (p c) gs)))
          (when (prime? (r c)) (set! gs (cons (r c) gs))))))))

(for-each
  (lambda (p) (display p) (newline))
  (germain 25 100000 100000))

Private [?]

Run code