;lehman's factoring algorithm
(use-modules (srfi srfi-60))
(define (isqrt n)
(if (not (and (positive? n) (integer? n)))
(error 'isqrt "must be positive integer")
(let loop ((x n))
(let ((y (quotient (+ x (quotient n x)) 2)))
(if (< y x) (loop y) x)))))
(define (iroot k n)
(let ((k-1 (- k 1)))
(let loop ((u n) (s (+ n 1)))
(if (<= s u) s
(loop (quotient (+ (* k-1 u) (quotient n (expt u k-1))) k) u)))))
(define (square? n) ; guile needs (use-modules (srfi srfi-60))
(if (negative? n) #f (if (zero? n) 0
(let ((m (modulo n 128)))
(if (positive? (bitwise-and (* m #x8bc40d7d) (* m #xa1e2f5d1) #x14020a)) #f
(let ((large-mod (modulo n 3989930175))) ; (* 63 25 11 17 19 23 31)
(and (let ((m (modulo large-mod 63)))
(zero? (bitwise-and (* m #x3d491df7) (* m #xc824a9f9) #x10f14008)))
(let ((m (modulo large-mod 25)))
(zero? (bitwise-and (* m #x1929fc1b) (* m #x4c9ea3b2) #x51001005)))
(let ((m (* #xd10d829a (modulo large-mod 31))))
(zero? (bitwise-and m (+ m #x672a5354) #x21025115)))
(let ((m (modulo large-mod 23)))
(zero? (bitwise-and (* m #x7bd28629) (* m #xe7180889) #xf8300)))
(let ((m (modulo large-mod 19)))
(zero? (bitwise-and (* m #x1b8bead3) (* m #x4d75a124) #x4280082b)))
(let ((m (modulo large-mod 17)))
(zero? (bitwise-and (* m #x6736f323) (* m #x9b1d499) #xc0000300)))
(let ((m (modulo large-mod 11)))
(zero? (bitwise-and (* m #xabf1a3a7) (* m #x2612bf93) #x45854000)))
(let ((root (isqrt n))) (if (= (* root root) n) root #f)))))))))
(define (wheel n) ; factoring by 2,3,5-wheel
(let ((wheel '#(1 2 2 4 2 4 2 4 6 2 6)))
(let loop ((n n) (f 2) (fs (list)) (w 0))
(if (< n (* f f)) (reverse (cons n fs))
(if (zero? (modulo n f))
(loop (/ n f) f (cons f fs) w)
(loop n (+ f (vector-ref wheel w)) fs
(if (= w 10) 3 (+ w 1))))))))
(define (fermat n) ; fermat's factoring algorithm
(call-with-current-continuation (lambda (return)
(do ((a (inexact->exact (ceiling (sqrt n))) (+ a 1)))
((< (quotient (+ n 9) 6) a) n)
(cond ((square? (- (* a a) n)) =>
(lambda (b) (return (- a b)))))))))
(define (factors n) ; lehman's factoring algorithm
(define (wheel n limit)
(let ((wheel '#(1 2 2 4 2 4 2 4 6 2 6)))
(let loop ((n n) (f 2) (fs (list)) (w 0))
(if (< limit f) (values n fs)
(if (< n (* f f)) (values 1 (cons n fs))
(if (zero? (modulo n f)) (loop (/ n f) f (cons f fs) w)
(loop n (+ f (vector-ref wheel w)) fs (if (= w 10) 3 (+ w 1)))))))))
(define (lehman n)
(call-with-current-continuation (lambda (return)
(do ((k 1 (+ k 1))) ((< (expt n 1/3) k))
(let* ((4kn (* 4 k n)) (sqrt-4kn (sqrt 4kn)))
(do ((a (inexact->exact (ceiling sqrt-4kn)) (+ a 1)))
((< (floor (+ sqrt-4kn (/ (expt n 1/6) (* 4 (sqrt k))))) a))
(let ((b (square? (- (* a a) 4kn))))
(if b (return (gcd (+ a b) n)))))))
(return n))))
(let ((limit (iroot 3 n)))
(call-with-values (lambda () (wheel n limit))
(lambda (n fs)
(if (= n 1) (reverse fs)
(let* ((p (lehman n)) (q (/ n p)))
(if (= q 1) (reverse (cons p fs))
(reverse (cons (max p q) (cons (min p q) fs))))))))))
(display (wheel 13290059)) (newline)
(display (fermat 13290059)) (newline)
(display (factors 13290059)) (newline)