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