[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Aug 1:
; 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))


Output:
1831829999999999999999999999999
48873824999999999999999999999999
82912829999999999999999999999999
91351259999999999999999999999999
97747649999999999999999999999999
97747649999999999999999999999999
119339219999999999999999999999999
119339219999999999999999999999999
120585464999999999999999999999999
132417284999999999999999999999999
154609454999999999999999999999999
155705549999999999999999999999999
163258094999999999999999999999999
166291124999999999999999999999999
168468299999999999999999999999999
168468299999999999999999999999999
169729559999999999999999999999999
177747569999999999999999999999999
184849664999999999999999999999999
189309119999999999999999999999999
195495299999999999999999999999999
207086879999999999999999999999999
208888679999999999999999999999999
209309099999999999999999999999999
209309099999999999999999999999999
210194984999999999999999999999999
217176959999999999999999999999999
231366134999999999999999999999999
236411174999999999999999999999999
240164924999999999999999999999999
251411159999999999999999999999999
285795509999999999999999999999999
295254959999999999999999999999999
306561254999999999999999999999999
323963639999999999999999999999999
326035709999999999999999999999999
329459129999999999999999999999999
347266919999999999999999999999999
367266899999999999999999999999999
380029649999999999999999999999999
384504119999999999999999999999999
396290894999999999999999999999999
420389969999999999999999999999999
458933474999999999999999999999999
460525064999999999999999999999999
479864384999999999999999999999999
480329849999999999999999999999999
487897409999999999999999999999999
507957449999999999999999999999999
507957449999999999999999999999999
511696184999999999999999999999999
513873359999999999999999999999999
514023509999999999999999999999999
531320789999999999999999999999999
549488939999999999999999999999999
561275714999999999999999999999999
571696124999999999999999999999999
579894314999999999999999999999999
581711129999999999999999999999999
591065474999999999999999999999999
603317714999999999999999999999999
635629994999999999999999999999999
647581934999999999999999999999999
658918259999999999999999999999999
659233574999999999999999999999999
676065389999999999999999999999999
686876189999999999999999999999999
694533839999999999999999999999999
694533839999999999999999999999999
718677959999999999999999999999999
732446714999999999999999999999999
734488754999999999999999999999999
760059299999999999999999999999999
760059299999999999999999999999999
766440674999999999999999999999999
779188409999999999999999999999999
787581794999999999999999999999999
807611804999999999999999999999999
832942109999999999999999999999999
837431594999999999999999999999999
850269419999999999999999999999999
852221369999999999999999999999999
853002149999999999999999999999999
874068194999999999999999999999999
877056179999999999999999999999999
877491614999999999999999999999999
886650764999999999999999999999999
887176289999999999999999999999999
892671779999999999999999999999999
899443544999999999999999999999999
902671769999999999999999999999999
902671769999999999999999999999999
907041134999999999999999999999999
908797889999999999999999999999999
911260349999999999999999999999999
921050129999999999999999999999999
921560639999999999999999999999999
931560629999999999999999999999999
933347414999999999999999999999999
947626679999999999999999999999999
951935984999999999999999999999999
959728769999999999999999999999999
959728769999999999999999999999999
964713749999999999999999999999999
971095124999999999999999999999999
990329339999999999999999999999999
999668669999999999999999999999999
999668669999999999999999999999999
1025854829999999999999999999999999
1030914884999999999999999999999999
1039308269999999999999999999999999
1047971924999999999999999999999999
1049143094999999999999999999999999
1052866814999999999999999999999999
1062671609999999999999999999999999
1077146069999999999999999999999999
1078602524999999999999999999999999
1083947864999999999999999999999999
1084668584999999999999999999999999
1089398309999999999999999999999999
1104533429999999999999999999999999
1126350224999999999999999999999999
1132040909999999999999999999999999
1132040909999999999999999999999999
1138302164999999999999999999999999
1159788629999999999999999999999999
1177686509999999999999999999999999
1180884704999999999999999999999999
1181064884999999999999999999999999
1188977789999999999999999999999999
1189908719999999999999999999999999
1204142939999999999999999999999999
1230839609999999999999999999999999
1230839609999999999999999999999999
1239368129999999999999999999999999
1246154909999999999999999999999999
1277401124999999999999999999999999
1285389104999999999999999999999999
1296320024999999999999999999999999
1318467149999999999999999999999999
1323542219999999999999999999999999
1325719394999999999999999999999999
1332521189999999999999999999999999
1340058719999999999999999999999999
1340058719999999999999999999999999
1342386044999999999999999999999999
1342821479999999999999999999999999
1344353009999999999999999999999999
1346560214999999999999999999999999
1347025679999999999999999999999999
1347671324999999999999999999999999
1352130779999999999999999999999999
1395328934999999999999999999999999
1404908504999999999999999999999999
1425058634999999999999999999999999
1427566139999999999999999999999999
1438151714999999999999999999999999
1439848409999999999999999999999999
1447340894999999999999999999999999
1447400954999999999999999999999999
1470358889999999999999999999999999
1470358889999999999999999999999999
1473121649999999999999999999999999
1485914429999999999999999999999999
1485914429999999999999999999999999
1491154664999999999999999999999999
1492671179999999999999999999999999
1492671179999999999999999999999999
1493151659999999999999999999999999
1493181689999999999999999999999999
1502370869999999999999999999999999
1538346809999999999999999999999999
1558376819999999999999999999999999
1615223609999999999999999999999999
1616454839999999999999999999999999
1653752099999999999999999999999999
1677265589999999999999999999999999
1706004299999999999999999999999999
1719728009999999999999999999999999
1798887089999999999999999999999999
1814082269999999999999999999999999
1866694829999999999999999999999999
2013421409999999999999999999999999
2074172099999999999999999999999999
2090478389999999999999999999999999
2096003909999999999999999999999999
2105733629999999999999999999999999
2126874749999999999999999999999999
2264081819999999999999999999999999
2332249919999999999999999999999999
2355373019999999999999999999999999
2361769409999999999999999999999999
2408285879999999999999999999999999
2438676239999999999999999999999999
2443150709999999999999999999999999
2459456999999999999999999999999999
2467715249999999999999999999999999
2469306839999999999999999999999999
2554802249999999999999999999999999
2595582989999999999999999999999999
2597354759999999999999999999999999
2629276649999999999999999999999999
2791768979999999999999999999999999
2809817009999999999999999999999999
2838315479999999999999999999999999
2866273409999999999999999999999999
2894681789999999999999999999999999
2900807909999999999999999999999999
2916273359999999999999999999999999
2940717779999999999999999999999999


Create a new paste based on this one


Comments: