As some have seen me talking in IRC, I have my own rendition of Fermat factoring. I compared it to what I thought was a typical approach and my method asymptotically replaces a square root operation with a subtract and an increment. In complexity terms:

Square root is on the level of multiplication and division, so O(M(n))
Subtraction is O(n)
Increment, if my calculations are correct, is O(1) (in the case of the algorithm described below, invert bit 1, if that returns 0, invert bit 2, ... to increment by 2)

So while Fermat factoring isn't as efficient as the general number field sieve, the method I have is more efficient than a typical approach to Fermat factoring.

---

To describe the idea, notice that a²-b²=(a-b)(a+b). This means that, given our number c, if we can find a difference of two squares so that they equal c, then we have a pair of factors for c. For example, 10²-5²=75, so 75=(10-5)(10+5)=5*15. If c>1 and is a positive odd integer, it will always have a trivial solution (the difference between consecutive squares is 1,3,5,7,..., so this is "trivial"). Any other solutions are non-trivial and imply that c is composite. If only the trivial solution exists, c is prime.
The next few steps toward solving this formula are the typical approaches and I do the same thing. We notice that a>sqrt(c).
So then with our algorithm, we can initialise a=cieling(sqrt(c). That is the smallest that a can be, so now we enter the naive algorithm, checking if the difference a²-c is a perfect square:
ceiling(sqrt(c))→a
While notperfectsqaure(a*a-c)
a+1→a
EndWhile
sqrt(a*a-c)→b
Return {a-b,a+b}
However, the next step away from the naive approach is to notice that every time you add 1 to a, the difference increases by 2a+1. So you can save some computations by keeping track of the difference and updating at each iteration:
ceiling(sqrt(c))→a
a*a-c→b2
While notperfectsqaure(b2)
b2+2a+1→b2
a+1→a
EndWhile
sqrt(b2)→b
Return {a-b,a+b}
However, that notperfectsquare() routine will likely require computing the square root at least once (twice if done naively). My approach takes the optimisation another step further by keeping track of how far away b2 is from a perfect square. With this, note that if b2 is the next lowest perfect square below a²-c, then the one above is b2+2*b2+1, so if the difference between b and b2 exceeds 2*b2+1, we need to increment b and adjust the difference:
ceiling(sqrt(c))→a
a*a-c→diff
floor(sqrt(diff))→b
diff-b*b→diff
While diff>0
  diff+a+a+1→diff
  a+1→a
  While diff>=(2b+1)
    diff-b-b-1→diff
    b+1→b
  EndWhile
EndWhile
Return {a-b,a+b}
As it is, the inner while loop converges to iterating 1 time pretty quickly (bounded the square root of the difference between the initial guess for a and c). We can say then that the original square root required at each iteration is replaced by 1 subtract and 1 increment.

Rewriting this for some more efficiency, we can do:
ceiling(sqrt(c))→a
a*a-c→diff
floor(sqrt(diff))→b
diff-b*b→diff
2a+1→a
2b+1→b
While diff>0
  diff+a→diff
  a+2→a
  While diff>=b
    diff-b→diff
    b+2→b
  EndWhile
EndWhile
floor(a/2)→a
floor(b/2)→b
Return {a-b,a+b}
At the chip level, those +2 can also be increments requiring O(1) time instead of additions requiring O(n) time. This means that each iteration of the algorithm requires about:
- 2 add/sub operations
- 2 increments
- 3 compares (that can be changed to simply checking for overflow, so these can be trivially removed)

And the maximum number of iterations are based on how far the actual a value is from the original estimate. If we know C is an odd composite (as, for example, RSAs) then this is a pretty crummy upper bound of c/3-sqrt(c) in the worst case (c=3*p).