Fermat factorization method limit

I am trying to implement Fermat's factorization (Algorithm C in The Art of Computer Programming, Vol. 2). Unfortunately in my edition (ISBN 81-7758-335-2), this algorithm is printed incorrectly. what should be the condition on factor-inner loop below? I am running the loop till y <= n [passed in as limit].

(if (< limit y) 0 (factor-inner x (+ y 2) (- r y) limit))

Is there anyway to avoid th开发者_C百科is condition altogether, as it will double the speed of loop?

(define (factor n) 
  (let ((square-root (inexact->exact (floor (sqrt n))))) 
    (factor-inner (+ (* 2 square-root) 1)
                  1 
                  (- (* square-root square-root) n)
                  n)))

(define (factor-inner x y r limit)
  (if (= r 0)
    (/ (- x y) 2)
    (begin 
      (display x) (display " ") (display y) (display " ") (display r) (newline)
      ;;(sleep-current-thread 1)
      (if (< r 0)
        (factor-inner (+ x 2) y (+  r x) limit)
        (if (< limit y)
          0
          (factor-inner x (+ y 2) (- r y) limit))))))

The (< limit y) check is not necessary because, worst-case, the algorithm will eventually find this solution:

x = N + 2

y = N

It will then return 1.

Looking through Algorithm C, it looks like the issue is with the recursion step, which effectively skips step C4 whenever r < 0, because x is not incremented and r is only decremented by y.

Using the notation of a, b and r from the 1998 edition of Vol. 2 (ISBN 0-201-89684-2), a Scheme version would be as follows:

(define (factor n) 
  (let ((x (inexact->exact (floor (sqrt n)))))
    (factor-inner (+ (* x 2) 1)
                  1
                  (- (* x x) n))))

(define (factor-inner a b r)
  (cond ((= r 0) (/ (- a b) 2))
        ((< 0 r) (factor-inner a       (+ b 2) (- r b)))
        (else    (factor-inner (+ a 2) (+ b 2) (- r (- a b))))))

EDIT to add: Basically, we are doing a trick that repeatedly checks whether

r <- ((a - b) / 2)*((a + b - 2)/2) - N

is 0, and we're doing it by simply tracking how r changes when we increment a or b. If we were to set b to b+2 in the expression for r above, it's equivalent to reducing r by the old value of b, which is why both are done in parallel in step C4 of the algorithm. I encourage you to expand out the algebraic expression above and convince yourself that this is true.

As long as r > 0, you want to keep decreasing it to find the right value of b, so you keep repeating step C4. However, if you overshoot, and r < 0, you need to increase it. You do this by increasing a, because increasing a by 2 is equivalent to decreasing r by the old value of a, as in step C3. You will always have a > b, so increasing r by a in step C3 automatically makes r positive again, so you just proceed directly on to step C4.

It's also easy to prove that a > b. We start with a manifestly greater than b, and if we ever increase b to the point where b = a - 2, we have

N = (a - (a - 2))/2 * ((a + (a - 2) - 2)/2 = 1 * (a - 2)

This means that N is prime, as the largest factor it has that is less than sqrt(N) is 1, and the algorithm has terminated.