fork download
copy 
  1. ; prime -- a modest library
  2.  
  3. ; primes n -- list of primes not greater than n in ascending order
  4. (define (primes n) ; assumes n is an integer greater than one
  5.   (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
  6.     (let loop ((i 0) (p 3) (ps (list 2))) ; sieve of eratosthenes
  7.       (cond ((< n (* p p))
  8.               (do ((i i (+ i 1)) (p p (+ p 2))
  9.                    (ps ps (if (vector-ref bits i) (cons p ps) ps)))
  10.                   ((= i len) (reverse ps))))
  11.             ((vector-ref bits i)
  12.               (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
  13.                   ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
  14.                 (vector-set! bits j #f)))
  15.             (else (loop (+ i 1) (+ p 2) ps))))))
  16.  
  17. ; prime? n -- #f if provably composite, else #t if probably prime
  18. (define prime? ; strong pseudoprime to prime bases less than 100
  19.   (let ((ps (primes 100))) ; assumes n an integer greater than one
  20.     (lambda (n)
  21.       (define (expm b e m)
  22.         (define (times p q) (modulo (* p q) m))
  23.           (let loop ((b b) (e e) (x 1))
  24.             (if (zero? e) x
  25.               (loop (times b b) (quotient e 2)
  26.                     (if (odd? e) (times b x) x)))))
  27.       (define (spsp? n a) ; #t if n is a strong pseudoprime base a
  28.         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  29.             ((odd? d)
  30.               (if (= (expm a d n) 1) #t
  31.                 (do ((r 0 (+ r 1)))
  32.                     ((or (= (expm a (* d (expt 2 r)) n) (- n 1))
  33.                          (= r s))
  34.                       (< r s)))))))
  35.       (if (member n ps) #t
  36.         (do ((ps ps (cdr ps)))
  37.             ((or (null? ps) (not (spsp? n (car ps))))
  38.               (null? ps)))))))
  39.  
  40. ; factors n -- list of prime factors of n in ascending order
  41. (define (factors n) ; assumes n is an integer, may be negative
  42.   (if (<= -1 n 1) (list n) ; pollard's rho algorithm
  43.     (if (< n 0) (cons -1 (factors (- n)))
  44.       (let fact ((n n) (c 1) (fs (list)))
  45.         (define (f x) (modulo (+ (* x x) c) n))
  46.         (if (even? n) (fact (/ n 2) c (cons 2 fs))
  47.           (if (= n 1) fs
  48.             (if (prime? n) (sort (cons n fs) <)
  49.               (let loop ((t 2) (h 2) (d 1))
  50.                 (cond ((= d 1) (let ((t (f t)) (h (f (f h))))
  51.                                  (loop t h (gcd (- t h) n))))
  52.                       ((= d n) (fact n (+ c 1) fs)) ; cyclic
  53.                       ((prime? d)
  54.                         (fact (/ n d) (+ c 1) (cons d fs)))
  55.                       (else (fact n (+ c 1) fs)))))))))))
  56.  
  57. (display (primes 50)) (newline)
  58.  
  59. (do ((n 1 (+ n 1)))
  60.     ((< 98 n))
  61.   (display n)
  62.   (display ": ")
  63.   (for-each
  64.     (lambda (x)
  65.       (display " ")
  66.       (display x))
  67.     (factors (- (expt 2 n) 1)))
  68.   (newline))
Time limit exceeded #stdin #stdout 15s 5336KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/07Pjb
language:
Scheme (guile 2.2.4)
created:
14 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!