fork(2) download
copy 
  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
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
(3119 4261)
3119
(3119 4261)
https://ideone.com/b3AN9T
language:
Scheme (guile 2.2.4)
created:
6 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!