fork download
  1. ;lehman's factoring algorithm
  2.  
  3. (use-modules (srfi srfi-60))
  4.  
  5. (define (isqrt n)
  6. (if (not (and (positive? n) (integer? n)))
  7. (error 'isqrt "must be positive integer")
  8. (let loop ((x n))
  9. (let ((y (quotient (+ x (quotient n x)) 2)))
  10. (if (< y x) (loop y) x)))))
  11.  
  12. (define (iroot k n)
  13. (let ((k-1 (- k 1)))
  14. (let loop ((u n) (s (+ n 1)))
  15. (if (<= s u) s
  16. (loop (quotient (+ (* k-1 u) (quotient n (expt u k-1))) k) u)))))
  17.  
  18. (define (square? n) ; guile needs (use-modules (srfi srfi-60))
  19. (if (negative? n) #f (if (zero? n) 0
  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.  
  39. (define (wheel n) ; factoring by 2,3,5-wheel
  40. (let ((wheel '#(1 2 2 4 2 4 2 4 6 2 6)))
  41. (let loop ((n n) (f 2) (fs (list)) (w 0))
  42. (if (< n (* f f)) (reverse (cons n fs))
  43. (if (zero? (modulo n f))
  44. (loop (/ n f) f (cons f fs) w)
  45. (loop n (+ f (vector-ref wheel w)) fs
  46. (if (= w 10) 3 (+ w 1))))))))
  47.  
  48. (define (fermat n) ; fermat's factoring algorithm
  49. (call-with-current-continuation (lambda (return)
  50. (do ((a (inexact->exact (ceiling (sqrt n))) (+ a 1)))
  51. ((< (quotient (+ n 9) 6) a) n)
  52. (cond ((square? (- (* a a) n)) =>
  53. (lambda (b) (return (- a b)))))))))
  54.  
  55. (define (factors n) ; lehman's factoring algorithm
  56. (define (wheel n limit)
  57. (let ((wheel '#(1 2 2 4 2 4 2 4 6 2 6)))
  58. (let loop ((n n) (f 2) (fs (list)) (w 0))
  59. (if (< limit f) (values n fs)
  60. (if (< n (* f f)) (values 1 (cons n fs))
  61. (if (zero? (modulo n f)) (loop (/ n f) f (cons f fs) w)
  62. (loop n (+ f (vector-ref wheel w)) fs (if (= w 10) 3 (+ w 1)))))))))
  63. (define (lehman n)
  64. (call-with-current-continuation (lambda (return)
  65. (do ((k 1 (+ k 1))) ((< (expt n 1/3) k))
  66. (let* ((4kn (* 4 k n)) (sqrt-4kn (sqrt 4kn)))
  67. (do ((a (inexact->exact (ceiling sqrt-4kn)) (+ a 1)))
  68. ((< (floor (+ sqrt-4kn (/ (expt n 1/6) (* 4 (sqrt k))))) a))
  69. (let ((b (square? (- (* a a) 4kn))))
  70. (if b (return (gcd (+ a b) n)))))))
  71. (return n))))
  72. (let ((limit (iroot 3 n)))
  73. (call-with-values (lambda () (wheel n limit))
  74. (lambda (n fs)
  75. (if (= n 1) (reverse fs)
  76. (let* ((p (lehman n)) (q (/ n p)))
  77. (if (= q 1) (reverse (cons p fs))
  78. (reverse (cons (max p q) (cons (min p q) fs))))))))))
  79.  
  80. (display (wheel 13290059)) (newline)
  81. (display (fermat 13290059)) (newline)
  82. (display (factors 13290059)) (newline)
Success #stdin #stdout 0.04s 43496KB
stdin
Standard input is empty
stdout
(3119 4261)
3119
(3119 4261)