fork download
  1. ; http://programmingpraxis.com/contents/themes/#Prime Numbers
  2.  
  3. (define (expm b e m)
  4. (define (times p q) (modulo (* p q) m))
  5. (let loop ((b b) (e e) (x 1))
  6. (if (zero? e) x
  7. (loop (times b b) (quotient e 2)
  8. (if (odd? e) (times b x) x)))))
  9.  
  10. (define (primes n) ; assumes n is an integer greater than one
  11. (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
  12. (let loop ((i 0) (p 3) (ps (list 2)))
  13. (cond ((< n (* p p))
  14. (do ((i i (+ i 1)) (p p (+ p 2))
  15. (ps ps (if (vector-ref bits i) (cons p ps) ps)))
  16. ((= i len) (reverse ps))))
  17. ((vector-ref bits i)
  18. (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
  19. ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
  20. (vector-set! bits j #f)))
  21. (else (loop (+ i 1) (+ p 2) ps))))))
  22.  
  23. (define prime? (let ((ps (primes 100)))
  24. (lambda (n) ; assumes n is an integer greater than one
  25. (define (spsp? n a)
  26. (define (f d s)
  27. (do ((r 0 (+ r 1)))
  28. ((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
  29. (< r s))))
  30. (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  31. ((odd? d) (if (= (expm a d n) 1) #t (f d s)))))
  32. (if (member n ps) #t
  33. (do ((ps ps (cdr ps)))
  34. ((or (null? ps) (not (spsp? n (car ps)))) (null? ps)))))))
  35.  
  36. (define (factors n)
  37. (if (<= -1 n 1) (list n) (if (< n 0) (cons -1 (factors (- n)))
  38. (let fact ((n n) (c 1) (fs (list)))
  39. (define (f x) (modulo (+ (* x x) c) n))
  40. (if (even? n) (fact (/ n 2) c (cons 2 fs))
  41. (if (= n 1) fs (if (prime? n) (sort (cons n fs) <)
  42. (let loop ((t 2) (h 2) (d 1))
  43. (cond ((= d 1) (let ((t (f t)) (h (f (f h))))
  44. (loop t h (gcd (- t h) n))))
  45. ((= d n) (fact n (+ c 1) fs)) ; cyclic
  46. ((prime? d) (fact (/ n d) (+ c 1) (cons d fs)))
  47. (else (fact n (+ c 1) fs)))))))))))
  48.  
  49. (display (factors 41748850938502584251))
Success #stdin #stdout 0.15s 4664KB
stdin
Standard input is empty
stdout
(6461335039 6461335109)