fork download
copy 
  1. ; generating large random primes
  2.  
  3. (define rand
  4.   (let ((seed 0.3141592654))
  5.     (lambda args
  6.       (set! seed
  7.         (if (pair? args)
  8.             (sin (car args))
  9.             (let ((x (* seed 147.0)))
  10.               (- x (floor x)))))
  11.       seed)))
  12.  
  13. (define (randint first past)
  14.   (inexact->exact (floor (+ (* (- past first) (rand)) first))))
  15.  
  16. (define (primes n) ; list of primes not exceeding n
  17.   (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
  18.     (let loop ((i 0) (p 3) (ps (list 2)))
  19.       (cond ((< n (* p p))
  20.               (do ((i i (+ i 1)) (p p (+ p 2))
  21.                    (ps ps (if (vector-ref bits i) (cons p ps) ps)))
  22.                   ((= i len) (reverse ps))))
  23.             ((vector-ref bits i)
  24.               (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
  25.                   ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
  26.                 (vector-set! bits j #f)))
  27.             (else (loop (+ i 1) (+ p 2) ps))))))
  28.  
  29. (define prime? ; strong pseudoprime to prime bases less than 100
  30.   (let* ((ps (primes 100)) (p100 (apply * ps)))
  31.     (lambda (n)
  32.       (define (expm b e m)
  33.         (let loop ((b b) (e e) (x 1))
  34.           (if (zero? e) x
  35.             (loop (modulo (* b b) m) (quotient e 2)
  36.                   (if (odd? e) (modulo (* b x) m) x)))))
  37.       (define (spsp? n a) ; #t if n is a strong pseudoprime base a
  38.         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  39.             ((odd? d) (if (= (expm a d n) 1) #t
  40.               (do ((r 0 (+ r 1)))
  41.                   ((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
  42.                     (< r s)))))))
  43.       (if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
  44.         (do ((ps ps (cdr ps)))
  45.             ((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))
  46.  
  47. (define sievers (primes (expt 2 16)))
  48.  
  49. (define (rand-prime lo hi)
  50.   (define (rand-odd lo hi)
  51.     (let ((n (randint lo hi)))
  52.       (if (odd? n) n (+ n 1))))
  53.   (let outer-loop ((base (rand-odd lo (- hi 100000))))
  54.     (let ((sieve (make-vector 50000 #t)))
  55.       (do ((ps (cdr sievers) (cdr ps))) ((null? ps))
  56.         (let ((p (car ps)))
  57.           (do ((i (modulo (* -1/2 (+ base p)) p) (+ i p)))
  58.                ((<= 50000 i)))))
  59.       (let inner-loop ((i 0))
  60.         (if (<= 50000 i) (outer-loop (rand-odd lo (- hi 100000)))
  61.           (if (and (vector-ref sieve i) (prime? (+ base i i))) (+ base i i)
  62.             (inner-loop (+ i 1))))))))
  63.  
  64. (display (rand-prime #e1e49 #e1e50)) (newline)
Success #stdin #stdout 1.01s 44472KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
26327081241999932390469193988834384819314000658501
https://ideone.com/jc57fP
language:
Scheme (guile 2.2.4)
created:
6 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!