isqrt is a mess

[stevengj]

The isqrt function is poorly documented, but a reasonable definition of isqrt(n) would be the largest integer m such that m*m <= n. This coincides with the definition:

isqrt(x::Integer) = oftype(x, trunc(sqrt(x)))

that we are using, more or less.

However, for Int128 values, the sqrt function converts to a Float64, which discards too many bits for the result to be accurate. e.g.

z = int128(typemax(Int64))
isqrt(z*z) - z

gives 1 rather than 0 (and there are other numbers where the difference is larger). A simple fix would be to use BigFloat with sufficient precision for isqrt(x::Int128), although I'm not sure if there is a more efficient solution (or whether we care).

[StefanKarpinski]

This feels like one of those cases where you can use Float64 to approximate and then patch up the answer quickly if it's not quite correct.

[StefanKarpinski]

The computation n>>((128-leading_zeros(n))>>1) gets you into the range [sqrt(n),sqrt(2n)].

[stevengj]

I wonder if you can use Newton's method, operating on Int128, to get within ±n for some small n, and then just exhaustively search.

[stevengj]

The following seems to work, though it may be suboptimal:

function isqrt2(x::Int128)
    s = convert(Int128, trunc(sqrt(x)))
    s = div(s + div(x),s), 2) # Newton step
    while s*s > x
        s -= 1
    end
    while (s+1)*(s+1) <= x
        s += 1
    end
    return s
end

[StefanKarpinski]

I can't prove that it's correct, but a single iteration of Newton's method with integer ops works in all cases I've tried:

function isqrt(x::Int128)
    s = convert(Int128,trunc(sqrt(x)))
    (s + div(x,s)) >> 1
end

[stevengj]

You're right, the "fixup" loops in my version never seem to be needed. I wonder why it's not off by one occasionally?

[StefanKarpinski]

I think you're always close enough from the floating-point approximation that a single Newton step gets you there.

[StefanKarpinski]

Ok, this is incredibly hand-wavy at best, but Newton's method for sqrt converges quadratically, so every iteration should double the number of accurate digits. Since the floating-point computation gives 52 correct bits, that means that one more iteration should give 104 correct bits. Since the largest possible sqrt for a 128-bit integer only requires 64 bits, this suffices. Of course, that's some pretty vigorous handwaving.

[RauliRuohonen]

I was also about to handwave the same way :) Anyhow, loops or assert would be nice anyway to ensure that wrong results are never returned no matter what (<- simple paranoia).

BTW, the "(s+1)(s+1) <= x" check won't work near typemax(Int128), because (s+1)(s+1) wraps around to negatives. "(s+1)*(s+1)-x <= 0" is better.

[stevengj]

I know an exact Newton iteration would suffice. I was worried more about the round off errors in the integer divisions, but upon reflection I think those can be bounded to show that the final answer can't be off by more than one.

[stevengj]

Or rather, that the final answer is off by less than 1.

[StefanKarpinski]

Btw, isqrt is wrong for Int64 as well:

julia> isqrt(9223372030926249000)^2
9223372030926249001

The same trick fixes it:

julia> function Base.isqrt(x::Int64)
           s = convert(Int64,trunc(sqrt(x)))
           (s + div(x,s)) >> 1
       end
isqrt (generic function with 4 methods)

julia> isqrt(9223372030926249000)^2
9223372024852248004

julia> isqrt(9223372030926249000)
3037000498

[stevengj]

I would think we can do something cheaper for Int64...

[StefanKarpinski]

Probably, although if it's more than adding the result of a comparison, it's probably not much cheaper.

[GunnarFarneback]

Here's an off by one example.

julia> a = int128(typemax(Int64))
9223372036854775807

julia> isqrt(a*a-1)
9223372036854775807

In general x*x-y for large x and small y are susceptible to having isqrt off by one by the current implementation.

[StefanKarpinski]

Yup, good find. That's definitely a problem. Looks like 128-bit at least may need two Newton iterations.

[stevengj]

You don't need another Newton iteration to correct an off-by-one error.

[StefanKarpinski]

Perhaps not, but I'm finding all sorts of problems with these now that I'm looking harder.

[jiahao]

Perhaps you can get a more systematic sense of what's needed by taking a random sample and looking at how many Newton iterations it takes for the solution to change by less than (say) 0.1.

[GunnarFarneback]

Throwing Newton iterations at this problem does not solve anything since they won't converge for x^2-1 numbers. The off by one solution x is moved to x-1, but another Newton iteration takes you back to x.

To make matters worse isqrt is currently rather ugly in the small end:

julia> isqrt(3)                                                                 
2

julia> isqrt(0)
ERROR: integer division error
 in isqrt at intfuncs.jl:317

[StefanKarpinski]

Yes, the former isqrt behavior was arguably better, even though it gave numbers larger than desired sometimes. What we have now is rather a mess. This is why I didn't want to just "throw Newton iterations" at it, as you put it.

[GunnarFarneback]

It's not that far off. Obviously 0 needs special casing and negative values could do with a more specific domain error than the one thrown from the sqrt call. Otherwise it's just an off by one correction that's needed.

Some of the handwaving above can be tightened up by the following observation:

If y1 = sqrt(x)*(1+e), one exact Newton iteration will give
y2 = 0.5 * (sqrt(x)*(1+e)+sqrt(x)/(1+e))
Using the identity 1/(1+e) = 1-e+e^2/(1+e) gives
y2 = sqrt(x)*0.5*(1+e+1-e+e^2/(1+e)) = sqrt(x)*(1+e^2/(2*(1+e)))
Assuming e < 1 we have
sqrt(x) <= y2 <= sqrt(x) + e^2/(2*(1+e))*sqrt(x)

Even quite crude bounds on the relative error e from the trunc(sqrt(x)) computation in double precision are sufficient to prove that sqrt(x) <= y2 <= sqrt(x) + 1 in exact mathematics and taking the integer divisions into account that sqrt(x) - 1 <= y2 <= sqrt(x) + 1.

[GunnarFarneback]

The assumption should be e>-1, obviously, so that 1+e doesn't switch sign.