fork(1) download
copy 
  1. ; programming with prime numbers
  2.  
  3. (define (primes n)
  4.   (if (or (not (integer? n)) (< n 2))
  5.       (error 'primes "must be integer greater than one")
  6.       (let* ((len (quotient (- n 1) 2))
  7.              (bits (make-vector len #t)))
  8.         (let loop ((i 0) (p 3) (ps (list 2)))
  9.           (cond ((< n (* p p))
  10.                   (do ((i i (+ i 1)) (p p (+ p 2))
  11.                        (ps ps (if (vector-ref bits i) (cons p ps) ps)))
  12.                       ((= i len) (reverse ps))))
  13.                 ((vector-ref bits i)
  14.                   (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
  15.                       ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
  16.                     (vector-set! bits j #f)))
  17.                 (else (loop (+ i 1) (+ p 2) ps)))))))
  18.  
  19. (define (td-prime? n . args)
  20.   (if (or (not (integer? n)) (< n 2))
  21.       (error 'td-prime? "must be integer greater than one")
  22.       (let ((limit (if (pair? args) (car args) 1000000)))
  23.         (if (even? n) (= n 2)
  24.           (let loop ((d 3))
  25.             (cond ((< limit d)
  26.                     (error 'td-prime? "limit exceeded"))
  27.                   ((< n (* d d)) #t)
  28.                   ((zero? (modulo n d)) #f)
  29.                   (else (loop (+ d 2)))))))))
  30.  
  31. (define (td-factors n . args)
  32.   (if (or (not (integer? n)) (< n 2))
  33.       (error 'td-factors "must be integer greater than one")
  34.       (let ((limit (if (pair? args) (car args) 1000000)))
  35.         (let twos ((n n) (fs '()))
  36.           (if (even? n)
  37.               (twos (/ n 2) (cons 2 fs))
  38.               (let odds ((n n) (d 3) (fs fs))
  39.                 (cond ((< limit d)
  40.                         (error 'td-factors "limit exceeded"))
  41.                       ((< n (* d d))
  42.                         (reverse (cons n fs)))
  43.                       ((zero? (modulo n d))
  44.                         (odds (/ n d) d (cons d fs)))
  45.                       (else (odds n (+ d 2) fs)))))))))
  46.  
  47. (define prime?
  48.   (let ((seed 3141592654))
  49.     (lambda (n)
  50.       (define (rand)
  51.         (set! seed (modulo (+ (* 69069 seed) 1234567) 4294967296))
  52.         (+ (quotient (* seed (- n 2)) 4294967296) 2))
  53.       (define (expm b e m)
  54.         (define (times x y) (modulo (* x y) m))
  55.         (let loop ((b b) (e e) (r 1))
  56.           (if (zero? e) r
  57.             (loop (times b b) (quotient e 2)
  58.                   (if (odd? e) (times b r) r)))))
  59.       (define (spsp? n a)
  60.         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  61.             ((odd? d)
  62.               (let ((t (expm a d n)))
  63.                 (if (or (= t 1) (= t (- n 1))) #t
  64.                   (do ((s (- s 1) (- s 1))
  65.                        (t (expm t 2 n) (expm t 2 n)))
  66.                       ((or (zero? s) (= t (- n 1)))
  67.                         (positive? s))))))))
  68.       (if (not (integer? n))
  69.           (error 'prime? "must be integer")
  70.           (if (< n 2) #f
  71.             (do ((a (rand) (rand)) (k 25 (- k 1)))
  72.                 ((or (zero? k) (not (spsp? n a)))
  73.                   (zero? k))))))))
  74.  
  75. (define (rho-factors n . args)
  76.   (define (cons< x xs)
  77.     (cond ((null? xs) (list x))
  78.           ((< x (car xs)) (cons x xs))
  79.           (else (cons (car xs) (cons< x (cdr xs))))))
  80.   (define (rho n limit)
  81.     (let loop ((t 2) (h 2) (d 1) (c 1) (limit limit))
  82.       (define (f x) (modulo (+ (* x x) c) n))
  83.       (cond ((zero? limit) (error 'rho-factors "limit exceeded"))
  84.             ((= d 1) (let ((t (f t)) (h (f (f h))))
  85.                        (loop t h (gcd (- t h) n) c (- limit 1))))
  86.             ((= d n) (loop 2 2 1 (+ c 1) (- limit 1)))
  87.             ((prime? d) d)
  88.             (else (rho d (- limit 1))))))
  89.   (if (not (integer? n))
  90.       (error 'rho-factors "must be integer")
  91.       (let ((limit (if (pair? args) (car args) 1000)))
  92.         (cond ((<= -1 n 1) (list n))
  93.               ((negative? n) (cons -1 (rho-factors (- n) limit)))
  94.               ((even? n)
  95.                 (if (= n 2) (list 2)
  96.                   (cons 2 (rho-factors (/ n 2) limit))))
  97.               (else (let loop ((n n) (fs '()))
  98.                       (if (prime? n)
  99.                           (cons< n fs)
  100.                           (let ((f (rho n limit)))
  101.                             (loop (/ n f) (cons< f fs))))))))))
  102.  
  103. (display (primes 100)) (newline)
  104. (display (length (primes 1000000))) (newline)
  105. (display (td-prime? 600851475143)) (newline)
  106. (display (td-factors 600851475143)) (newline)
  107. (display (prime? 600851475143)) (newline)
  108. (display (prime? 2305843009213693951)) (newline)
  109. (display (rho-factors 600851475143)) (newline)
Success #stdin #stdout 0.96s 5356KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
78498
#f
(71 839 1471 6857)
#f
#t
(71 839 1471 6857)
https://ideone.com/QC8bC
language:
Scheme (guile 2.2.4)
created:
11 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!