fork download 
  1. ; omega less than a thousand
  2.  
  3. (define sort #f)
  4. (define merge #f)
  5. (let ()
  6.   (define dosort
  7.     (lambda (pred? ls n)
  8.       (if (= n 1)
  9.           (list (car ls))
  10.           (let ((i (quotient n 2)))
  11.             (domerge pred?
  12.                      (dosort pred? ls i)
  13.                      (dosort pred? (list-tail ls i) (- n i)))))))
  14.   (define domerge
  15.     (lambda (pred? l1 l2)
  16.       (cond
  17.         ((null? l1) l2)
  18.         ((null? l2) l1)
  19.         ((pred? (car l2) (car l1))
  20.          (cons (car l2) (domerge pred? l1 (cdr l2))))
  21.         (else (cons (car l1) (domerge pred? (cdr l1) l2))))))
  22.   (set! sort
  23.     (lambda (pred? l)
  24.       (if (null? l) l (dosort pred? l (length l)))))
  25.   (set! merge
  26.     (lambda (pred? l1 l2)
  27.       (domerge pred? l1 l2))))
  28.  
  29. (define (unique eql? xs)
  30.   (cond ((null? xs) '())
  31.         ((null? (cdr xs)) xs)
  32.         ((eql? (car xs) (cadr xs))
  33.           (unique eql? (cdr xs)))
  34.         (else (cons (car xs) (unique eql? (cdr xs))))))
  35.  
  36. ; primes n -- list of primes not greater than n in ascending order
  37. (define (primes n) ; assumes n is an integer greater than one
  38.   (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
  39.     (let loop ((i 0) (p 3) (ps (list 2))) ; sieve of eratosthenes
  40.       (cond ((< n (* p p))
  41.               (do ((i i (+ i 1)) (p p (+ p 2))
  42.                    (ps ps (if (vector-ref bits i) (cons p ps) ps)))
  43.                   ((= i len) (reverse ps))))
  44.             ((vector-ref bits i)
  45.               (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
  46.                   ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
  47.                 (vector-set! bits j #f)))
  48.             (else (loop (+ i 1) (+ p 2) ps))))))
  49.  
  50. ; prime? n -- #f if provably composite, else #t if probably prime
  51. (define prime? ; strong pseudoprime to prime bases less than 100
  52.   (let ((ps (primes 100))) ; assumes n an integer greater than one
  53.     (lambda (n)
  54.       (define (expm b e m)
  55.         (define (times p q) (modulo (* p q) m))
  56.           (let loop ((b b) (e e) (x 1))
  57.             (if (zero? e) x
  58.               (loop (times b b) (quotient e 2)
  59.                     (if (odd? e) (times b x) x)))))
  60.       (define (spsp? n a) ; #t if n is a strong pseudoprime base a
  61.         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  62.             ((odd? d)
  63.               (if (= (expm a d n) 1) #t
  64.                 (do ((r 0 (+ r 1)))
  65.                     ((or (= (expm a (* d (expt 2 r)) n) (- n 1))
  66.                          (= r s))
  67.                       (< r s)))))))
  68.       (if (member n ps) #t
  69.         (do ((ps ps (cdr ps)))
  70.             ((or (null? ps) (not (spsp? n (car ps))))
  71.               (null? ps)))))))
  72.  
  73. ; factors n -- list of prime factors of n in ascending order
  74. (define (factors n) ; assumes n is an integer, may be negative
  75.   (if (<= -1 n 1) (list n) ; pollard's rho algorithm
  76.     (if (< n 0) (cons -1 (factors (- n)))
  77.       (let fact ((n n) (c 1) (fs (list)))
  78.         (define (f x) (modulo (+ (* x x) c) n))
  79.         (if (even? n) (fact (/ n 2) c (cons 2 fs))
  80.           (if (= n 1) fs
  81.             (if (prime? n) (sort < (cons n fs))
  82.               (let loop ((t 2) (h 2) (d 1))
  83.                 (cond ((= d 1) (let ((t (f t)) (h (f (f h))))
  84.                                  (loop t h (gcd (- t h) n))))
  85.                       ((= d n) (fact n (+ c 1) fs)) ; cyclic
  86.                       ((prime? d)
  87.                         (fact (/ n d) (+ c 1) (cons d fs)))
  88.                       (else (fact n (+ c 1) fs)))))))))))
  89.  
  90. (define (omega n) ; omega of a single number
  91.   (length (unique = (factors n))))
  92.  
  93. (define (omegas n) ; omegas less than n
  94.   (let* ((sieve (make-vector (+ n 1) 0)))
  95.     (do ((ps (primes (quotient n 2)) (cdr ps))) ((null? ps))
  96.       (do ((j (car ps) (+ j (car ps)))) ((<= n j))
  97.         (vector-set! sieve j (+ (vector-ref sieve j) 1))))
  98.     (do ((i 2 (+ i 1))) ((= i n))
  99.       (display i) (display #\tab)
  100.       (display (max (vector-ref sieve i) 1))
  101.       (display #\tab) (display (omega i))
  102.       (newline))))
  103.  
  104. (omegas 100)
Success #stdin #stdout 0.04s 4176KB
comments (?)
stdin 
Standard input is empty
stdout 
2	1	1
3	1	1
4	1	1
5	1	1
6	2	2
7	1	1
8	1	1
9	1	1
10	2	2
11	1	1
12	2	2
13	1	1
14	2	2
15	2	2
16	1	1
17	1	1
18	2	2
19	1	1
20	2	2
21	2	2
22	2	2
23	1	1
24	2	2
25	1	1
26	2	2
27	1	1
28	2	2
29	1	1
30	3	3
31	1	1
32	1	1
33	2	2
34	2	2
35	2	2
36	2	2
37	1	1
38	2	2
39	2	2
40	2	2
41	1	1
42	3	3
43	1	1
44	2	2
45	2	2
46	2	2
47	1	1
48	2	2
49	1	1
50	2	2
51	2	2
52	2	2
53	1	1
54	2	2
55	2	2
56	2	2
57	2	2
58	2	2
59	1	1
60	3	3
61	1	1
62	2	2
63	2	2
64	1	1
65	2	2
66	3	3
67	1	1
68	2	2
69	2	2
70	3	3
71	1	1
72	2	2
73	1	1
74	2	2
75	2	2
76	2	2
77	2	2
78	3	3
79	1	1
80	2	2
81	1	1
82	2	2
83	1	1
84	3	3
85	2	2
86	2	2
87	2	2
88	2	2
89	1	1
90	3	3
91	2	2
92	2	2
93	2	2
94	2	2
95	2	2
96	2	2
97	1	1
98	2	2
99	2	2
https://ideone.com/x1Sxw
language:
Scheme (guile 2.2.4)
created:
12 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!