fork(1) download
copy 
  1. ; display random factored numbers, easily
  2.  
  3. (define bitwise-and logand)
  4.  
  5. (define (prime? n) ; baillie-wagstaff
  6.   (define ps (list 2 3 5 7 11 13 17 19 23 29 31
  7.     37 41 43 47 53 59 61 67 71 73 79 83 89 97))
  8.   (define (expm b e m)
  9.     (define (m* x y) (modulo (* x y) m))
  10.     (cond ((zero? e) 1)
  11.           ((even? e) (expm (m* b b) (/ e 2) m))
  12.           (else (m* b (expm (m* b b) (/ (- e 1) 2) m))))) 
  13.   (define (square? n)
  14.     (define (isqrt n)
  15.       (if (not (and (positive? n) (integer? n)))
  16.           (error 'isqrt "must be positive integer")
  17.           (let loop ((x n))
  18.             (let ((y (quotient (+ x (quotient n x)) 2)))
  19.               (if (< y x) (loop y) x)))))
  20.     (let ((m (modulo n 128)))
  21.       (if (positive? (bitwise-and (* m #x8bc40d7d) (* m #xa1e2f5d1) #x14020a)) #f
  22.         (let ((large-mod (modulo n 3989930175))) ; (* 63 25 11 17 19 23 31)
  23.           (and (let ((m (modulo large-mod 63)))
  24.                  (zero? (bitwise-and (* m #x3d491df7) (* m #xc824a9f9) #x10f14008)))
  25.                (let ((m (modulo large-mod 25)))
  26.                  (zero? (bitwise-and (* m #x1929fc1b) (* m #x4c9ea3b2) #x51001005)))
  27.                (let ((m (* #xd10d829a (modulo large-mod 31))))
  28.                  (zero? (bitwise-and m (+ m #x672a5354) #x21025115)))
  29.                (let ((m (modulo large-mod 23)))
  30.                  (zero? (bitwise-and (* m #x7bd28629) (* m #xe7180889) #xf8300)))
  31.                (let ((m (modulo large-mod 19)))
  32.                  (zero? (bitwise-and (* m #x1b8bead3) (* m #x4d75a124) #x4280082b)))
  33.                (let ((m (modulo large-mod 17)))
  34.                  (zero? (bitwise-and (* m #x6736f323) (* m #x9b1d499) #xc0000300)))
  35.                (let ((m (modulo large-mod 11)))
  36.                  (zero? (bitwise-and (* m #xabf1a3a7) (* m #x2612bf93) #x45854000)))
  37.                (let ((root (isqrt n))) (if (= (* root root) n) root #f)))))))
  38.   (define (jacobi a m)
  39.     (if (not (integer? a)) (error 'jacobi "must be integer")
  40.       (if (not (and (integer? m) (positive? m) (odd? m)))
  41.           (error 'jacobi "modulus must be odd positive integer")
  42.           (let loop1 ((a (modulo a m)) (m m) (t 1))
  43.             (if (zero? a) (if (= m 1) t 0)
  44.               (let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))
  45.                 (let loop2 ((a a) (t t))
  46.                   (if (even? a) (loop2 (/ a 2) (* t z))
  47.                     (loop1 (modulo m a) a
  48.                            (if (and (= (modulo a 4) 3)
  49.                                     (= (modulo m 4) 3))
  50.                                (- t) t))))))))))
  51.   (define (strong-pseudoprime? n a)
  52.     (let loop ((r 0) (s (- n 1)))
  53.       (if (even? s) (loop (+ r 1) (/ s 2))
  54.         (if (= (expm a s n) 1) #t
  55.           (let loop ((r r) (s s))
  56.             (cond ((zero? r) #f)
  57.                   ((= (expm a s n) (- n 1)) #t)
  58.                   (else (loop (- r 1) (* s 2)))))))))(define (selfridge n)
  59.     (let loop ((d-abs 5) (sign 1))
  60.       (let ((d (* d-abs sign)))
  61.         (cond ((< 1 (gcd d n)) (values d 0 0))
  62.               ((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4)))
  63.               (else (loop (+ d-abs 2) (- sign)))))))
  64.   (define (lucas p q m n) ; right-to-left
  65.     (define (even e o) (if (even? n) e o))
  66.     (define (mod n) (if (zero? m) n (modulo n m)))
  67.     (let ((d (- (* p p) (* 4 q))))
  68.       (let loop ((un 1) (vn p) (qn q) (n (quotient n 2))
  69.                  (u (even 0 1)) (v (even 2 p)) (k (even 1 q)))
  70.         (if (zero? n) (values u v k)
  71.           (let ((u2 (mod (* un vn))) (v2 (mod (- (* vn vn) (* 2 qn))))
  72.                 (q2 (mod (* qn qn))) (n2 (quotient n 2)))
  73.             (if (even? n) (loop u2 v2 q2 n2 u v k)
  74.               (let* ((uu (+ (* u v2) (* u2 v)))
  75.                      (vv (+ (* v v2) (* d u u2)))
  76.                      (uu (if (and (positive? m) (odd? uu)) (+ uu m) uu))
  77.                      (vv (if (and (positive? m) (odd? vv)) (+ vv m) vv))
  78.                      (uu (mod (/ uu 2))) (vv (mod (/ vv 2))))
  79.                 (loop u2 v2 q2 n2 uu vv (* k q2)))))))))
  80.   (define (powers-of-two n)
  81.     (let loop ((s 0) (n n))
  82.       (if (odd? n) (values s n)
  83.         (loop (+ s 1) (/ n 2)))))
  84.   (define (strong-lucas-pseudoprime? n)
  85.     ; assumes odd positive integer not a square
  86.     (call-with-values
  87.       (lambda () (selfridge n))
  88.       (lambda (d p q)
  89.         (if (zero? p) (= n d)
  90.           (call-with-values
  91.             (lambda () (powers-of-two (+ n 1)))
  92.             (lambda (s t)
  93.               (call-with-values
  94.                 (lambda () (lucas p q n t))
  95.                 (lambda (u v k)
  96.                   (if (or (zero? u) (zero? v)) #t
  97.                     (let loop ((r 1) (v v) (k k))
  98.                       (if (= r s) #f
  99.                         (let* ((v (modulo (- (* v v) (* 2 k)) n))
  100.                                (k (modulo (* k k) n)))
  101.                           (if (zero? v) #t (loop (+ r 1) v k))))))))))))))
  102.   (if (not (integer? n)) (error 'prime? "must be integer"))
  103.   (if (or (< n 2) (square? n)) #f
  104.     (let loop ((ps ps))
  105.       (if (pair? ps)
  106.           (if (zero? (modulo n (car ps))) (= n (car ps)) (loop (cdr ps)))
  107.           (and (strong-pseudoprime? n 2)
  108.                (strong-pseudoprime? n 3)
  109.                (strong-lucas-pseudoprime? n))))))
  110.  
  111. ; composite number between lo and hi, with its prime factorization
  112. (define (random-composite lo hi) ; j cryptology (2003) 16: 287-289
  113.   (let loop ((n hi) (ss (list)))
  114.     (if (< 1 n) (let ((t (+ (random n) 1))) (loop t (cons t ss)))
  115.       (let* ((rs (filter prime? ss)) (r (apply * rs)))
  116.         (if (and (< 1 (length rs)) (< lo r hi) (< (random 1.0) (/ r hi)))
  117.             (values r rs)
  118.             (random-composite lo hi))))))
  119.  
  120. (call-with-values
  121.   (lambda () (random-composite #e1e9 #e1e10))
  122.   (lambda (r rs) (display r) (newline) (display rs) (newline)))
Success #stdin #stdout 0.98s 10960KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
7151251193
(79 89 223 4561)
https://ideone.com/106ySG
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!