;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)