fork download
copy 
  1. ; generating large random primes
  2.  
  3. (define (drop-while pred? xs)
  4.   (let loop ((xs xs))
  5.     (if (or (null? xs) (not (pred? (car xs)))) xs
  6.       (loop (cdr xs)))))
  7.  
  8. (define (fortune xs)
  9.   (let loop ((n 1) (x #f) (xs xs))
  10.     (cond ((null? xs) x)
  11.           ((< (rand) (/ n))
  12.             (loop (+ n 1) (car xs) (cdr xs)))
  13.           (else (loop (+ n 1) x (cdr xs))))))
  14.  
  15. (define (expm b e m)
  16.   (define (m* x y) (modulo (* x y) m))
  17.   (cond ((zero? e) 1)
  18.         ((even? e) (expm (m* b b) (/ e 2) m))
  19.         (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
  20.  
  21. (define rand
  22.   (let ((seed 0.3141592654))
  23.     (lambda args
  24.       (set! seed
  25.         (if (pair? args)
  26.             (sin (car args))
  27.             (let ((x (* seed 147.0)))
  28.               (- x (floor x)))))
  29.       seed)))
  30.  
  31. (define (randint first past)
  32.   (inexact->exact (floor (+ (* (- past first) (rand)) first))))
  33.  
  34. (define (primes n) ; list of primes not exceeding n
  35.   (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
  36.     (let loop ((i 0) (p 3) (ps (list 2)))
  37.       (cond ((< n (* p p))
  38.               (do ((i i (+ i 1)) (p p (+ p 2))
  39.                    (ps ps (if (vector-ref bits i) (cons p ps) ps)))
  40.                   ((= i len) (reverse ps))))
  41.             ((vector-ref bits i)
  42.               (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
  43.                   ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
  44.                 (vector-set! bits j #f)))
  45.             (else (loop (+ i 1) (+ p 2) ps))))))
  46.  
  47. (define prime? ; strong pseudoprime to prime bases less than 100
  48.   (let* ((ps (primes 100)) (p100 (apply * ps)))
  49.     (lambda (n)
  50.       (define (expm b e m)
  51.         (let loop ((b b) (e e) (x 1))
  52.           (if (zero? e) x
  53.             (loop (modulo (* b b) m) (quotient e 2)
  54.                   (if (odd? e) (modulo (* b x) m) x)))))
  55.       (define (spsp? n a) ; #t if n is a strong pseudoprime base a
  56.         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  57.             ((odd? d) (if (= (expm a d n) 1) #t
  58.               (do ((r 0 (+ r 1)))
  59.                   ((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
  60.                     (< r s)))))))
  61.       (if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
  62.         (do ((ps ps (cdr ps)))
  63.             ((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))
  64.  
  65. (define (pocklington p k a)
  66.   (and (odd? p) (prime? p)
  67.        (< 1 k (* 2 (+ p 1)))
  68.        (positive? (modulo k p))
  69.        (let ((n (+ (* 2 k p) 1)))
  70.          (and (< 1 a n)
  71.               (= (expm a (/ (- n 1) 2) n) (- n 1))
  72.               (= (gcd (+ (expm a k n) 1) n) 1)
  73.               n)))) ; return n if prime, else #f
  74.  
  75. (define (ratchet p lo hi verbose?)
  76.   (define (pock p k a n)
  77.     (and (= (expm a (/ (- n 1) 2) n) (- n 1))
  78.          (= (gcd (+ (expm a k n) 1) n) 1)))
  79.   (let loop ((k (randint lo hi)))
  80.     (let ((n (+ (* 2 k p) 1)))
  81.       (cond ((pock p k 2 n) (when verbose? (display (list p k 2)) (newline)) n)
  82.             ((pock p k 3 n) (when verbose? (display (list p k 3)) (newline)) n)
  83.             (else (loop (randint lo hi)))))))
  84.  
  85. (define (rand-prime lo hi . args)
  86.   (let ((verbose? (if (pair? args) (car args) #f)))
  87.     (if (< hi 10000) (fortune (drop-while (lambda (x) (< x lo)) (primes hi)))
  88.       (let loop ((p (fortune (cdr (primes 10000)))))
  89.         (if (<= lo p) p
  90.           (let ((k-lo (ceiling (/ (- lo 1) 2 p)))
  91.                 (k-hi (floor (/ (- hi 1) 2 p))))
  92.             (loop (ratchet p (if (< (* 2 p) k-lo) 1 k-lo)
  93.                            (min (* 2 p) k-hi) verbose?))))))))
  94.  
  95. (display (pocklington 2333 2001 2)) (newline)
  96. (display (pocklington 2333 2017 2)) (newline)
  97.  
  98. (display (rand-prime #e1e49 #e1e50 #t)) (newline)
Success #stdin #stdout 0.88s 44616KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
9336667
#f
(607 1022 2)
(1240709 147645 2)
(366368960611 227126727978 2)
(166424366512554387349117 236422002865808466771968 3)
(78692764113142983989384966585012763035290304513 618 2)
97264256443844728210879818699075775111618816378069
https://ideone.com/eVR8RA
language:
Scheme (guile 2.2.4)
created:
7 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!