Define L = (a+b)/2 and R = (a-b)/2, so that a * b = L² - R² for all a, b.
If there exist positive integers a, b such that a * b = n with |a-b| <= k n^(1/4) [constant k]
Then there exist L, R (either both integers or both half-integers) with L² - R² = n and |R| <= (1/2) k n^(1/4)
From which one can deduce:
R² <= (k²/4) √n
L² - R² >= L² - (k²/4) √n
n >= L² - (k²/4) √n
L² <= n + (k²/4) √n
L <= √(n + (k²/4) √n)
L <= √n + k²/8
If k <= 1, then L must be ceil(√n) or ceil(√n) - 1/2
If k <= 2√2, then L must be ceil(√n), ceil(√n) - 1/2, or ceil(√n) + 1/2
(my original estimate of 2 steps was due a slight miscalculation)
Adding epsilon would re-introduce an exponential component to the search. In fact, epsilon = 3/4 gives you total power. But this is only relevant if the reduction (https://cstheory.stackexchange.com/a/4785) can be made precise enough to bound L, U to within a fixed small epsilon. It could even be the case that you can bound to an arbitrarily small epsilon, but only at the expense of exploding other parameters. (B, T, k).
The funny thing is, it's often the case that already √(m^2 - 4N) > N^(1/4), so we can answer "No" to whether such a factorization exists, without even trying a single attempt.
From the continued fraction of √n, you get a series of best rational approximations a/b to √n. Each such approximation gives you:
a/b ~ √n
a² ~ b² n
The first approximation is a = floor(√n), b = 1
These approximations satisfy |a² - b² n| < 2 √n.
(this can be used for generating relations as input to Dixon's method)
I did make one mistake. I forgot to mention the half-integer case for x and y, which requires some additional logic. I implemented it fully python:
https://gdl.space/obiyukevit.py
N is 256-bit a sample number.