codepad

Language:

; shanks' square form factorization algorithm

(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 square?
  (let ((q11 (make-vector 11 #f))
        (q63 (make-vector 63 #f))
        (q64 (make-vector 64 #f))
        (q65 (make-vector 65 #f)))
    (do ((k 0 (+ k 1))) ((< 5 k))
      (vector-set! q11 (modulo (* k k) 11) #t))
    (do ((k 0 (+ k 1))) ((< 31 k))
      (vector-set! q63 (modulo (* k k) 63) #t))
    (do ((k 0 (+ k 1))) ((< 31 k))
      (vector-set! q64 (modulo (* k k) 64) #t))
    (do ((k 0 (+ k 1))) ((< 32 k))
      (vector-set! q65 (modulo (* k k) 65) #t))
    (lambda (n)
      (if (not (integer? n)) (error 'square? "must be integer")
        (if (< n 1) #f
          (if (not (vector-ref q64 (modulo n 64))) #f
            (let ((r (modulo n 45045)))
              (if (not (vector-ref q63 (modulo r 63))) #f
                (if (not (vector-ref q65 (modulo r 65))) #f
                  (if (not (vector-ref q11 (modulo r 11))) #f
                    (let ((q (isqrt n)))
                      (if (= (* q q) n) q #f))))))))))))

(define (rho f)
  (let* ((a (car f)) (b (cadr f)) (c (caddr f))
         (d (- (* b b) (* 4 a c))) (s (isqrt d)) (l (abs c))
         (r (if (< l (quotient s 2))
                (- (* 2 l (quotient (+ b s) (* 2 l))) b)
                (- (* 2 l) b))))
    (list c r (quotient (- (* r r) d) (* 4 c)))))

(define (reduce f)
  (let* ((a (car f)) (b (cadr f)) (c (caddr f))
         (d (- (* b b) (* 4 a c))) (s (isqrt d)))
    (if (< (abs (- s (* 2 (abs a)))) b s) f (reduce (rho f)))))

(define (squfof n)
  ; assumes n is odd, composite, and not a perfect square
  (define (enqueue f m q l)
    (let* ((x (abs (caddr f))) (g (quotient x (gcd x m))))
      (if (and (<= g l) (not (member g q))) (cons g q) q)))
  (define (inv-sqrt f c)
    (list (- c) (cadr f) (* -1 c (car f))))
  (let* ((n4 (= (modulo n 4) 1)) (m (if n4 1 2))
         (d (if n4 n (* 4 n))) (s (isqrt d))
         (b (if n4 (+ (* 2 (quotient (- s 1) 2)) 1) (* 2 (quotient s 2))))
         (delta (quotient (- (* b b) d) 4))
         (f (list 1 b delta)) (l (isqrt s)) (bound (* 4 l)))
   (let loop ((i 2) (f (rho f)) (q (enqueue f m '() l)))
      (cond ((or (< bound i) (= (caddr f) 1)) #f)
            ((square? (caddr f)) => (lambda (c)
              (if (member c q)
                  (loop (+ i 1) (rho f) (enqueue f m q l))
                  (let pool ((g (reduce (inv-sqrt f c))))
                    (let ((new-g (rho g)))
                      (if (= (cadr g) (cadr new-g))
                          (let ((c (abs (caddr g))))
                            (if (odd? c) c (/ c 2)))
                          (pool new-g)))))))
            (else (loop (+ i 1) (rho f) (enqueue f m q l)))))))

(display (squfof 633003781))

Private [?]

Run code