; generating large random primes

(define rand
  (let ((seed 0.3141592654))
    (lambda args
      (set! seed
        (if (pair? args)
            (sin (car args))
            (let ((x (* seed 147.0)))
              (- x (floor x)))))
      seed)))
 
(define (randint first past)
  (inexact->exact (floor (+ (* (- past first) (rand)) first))))

(define (primes n) ; list of primes not exceeding n
  (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
    (let loop ((i 0) (p 3) (ps (list 2)))
      (cond ((< n (* p p))
              (do ((i i (+ i 1)) (p p (+ p 2))
                   (ps ps (if (vector-ref bits i) (cons p ps) ps)))
                  ((= i len) (reverse ps))))
            ((vector-ref bits i)
              (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
                  ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
                (vector-set! bits j #f)))
            (else (loop (+ i 1) (+ p 2) ps))))))

(define prime? ; strong pseudoprime to prime bases less than 100
  (let* ((ps (primes 100)) (p100 (apply * ps)))
    (lambda (n)
      (define (expm b e m)
        (let loop ((b b) (e e) (x 1))
          (if (zero? e) x
            (loop (modulo (* b b) m) (quotient e 2)
                  (if (odd? e) (modulo (* b x) m) x)))))
      (define (spsp? n a) ; #t if n is a strong pseudoprime base a
        (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
            ((odd? d) (if (= (expm a d n) 1) #t
              (do ((r 0 (+ r 1)))
                  ((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
                    (< r s)))))))
      (if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
        (do ((ps ps (cdr ps)))
            ((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))

(define sievers (primes (expt 2 16)))

(define (rand-prime lo hi)
  (define (rand-odd lo hi)
    (let ((n (randint lo hi)))
      (if (odd? n) n (+ n 1))))
  (let outer-loop ((base (rand-odd lo (- hi 100000))))
    (let ((sieve (make-vector 50000 #t)))
      (do ((ps (cdr sievers) (cdr ps))) ((null? ps))
        (let ((p (car ps)))
          (do ((i (modulo (* -1/2 (+ base p)) p) (+ i p)))
               ((<= 50000 i)))))
      (let inner-loop ((i 0))
        (if (<= 50000 i) (outer-loop (rand-odd lo (- hi 100000)))
          (if (and (vector-ref sieve i) (prime? (+ base i i))) (+ base i i)
            (inner-loop (+ i 1))))))))

(display (rand-prime #e1e49 #e1e50)) (newline)