Challenge: optimizing isqrt

[Akira Li]

Steven D'Aprano <steve+comp.lang.python at pearwood.info> writes:

> There is an algorithm for calculating the integer square root of any
> positive integer using only integer operations:
>
> def isqrt(n):
>     if n < 0: raise ValueError
>     if n == 0:
>         return 0
>     bits = n.bit_length()
>     a, b = divmod(bits, 2)
>     x = 2**(a+b)
>     while True:
>         y = (x + n//x)//2
>         if y >= x:
>             return x
>         x = y
>
> This returns the integer part of the square root of n, that is, the greatest
> whole number less than or equal to the square root of n:
>
> py> isqrt(35)
> 5
> py> isqrt(36)
> 6
>
>
> That makes it equivalent to int(math.sqrt(n)), which also happens to be
> much, much faster, at least for small values of n. However, for large
> values of n, using floating point intermediate calculations fail:
>
> py> import math
> py> int(math.sqrt(2**3000))
> Traceback (most recent call last):
>   File "<stdin>", line 1, in <module>
> OverflowError: long int too large to convert to float
>
> Another problem is that, above a certain size, the resolution of floats is
> larger than 1, so you can't convert every int into a float without loss:
>
> py> float(2**90-1) == 2**90-1
> False
>
> which means that using math.sqrt is not correct:
>
> py> isqrt(2**90-1)
> 35184372088831
> py> int(math.sqrt(2**90-1))  # Off by one.
> 35184372088832
>
>
> So, the challenge is to identify when it is safe to optimise isqrt(n) as
> int(math.sqrt(n)):
>
> Q1: What is the largest value of M, such that 
>
> all(isqrt(i) == int(math.sqrt(n)) for n in range(M))
>
> returns True?
>
> I have done a brute force test, and in nine hours confirmed that M is at
> least 7627926244. That took nine hours, and I expect that a brute force
> test of every int representable as a float would take *months* of
> processing time.
>
> Q2: For values above M, is there a way of identifying which values of n are
> okay to use the optimized version?
>
> Q3: What is the largest value of n beyond which you can never use the float
> optimization?
>
>
> You can assume that Python floats are IEEE-754 C doubles, and that
> math.sqrt() is correctly rounded.

Where do you want to use your optimized isqrt(i)? 

There could be specialized algorithms that work only in narrow specific
circumstances e.g., the inverse square root (1/sqrt(x)) implementation
from Quake III Arena that has 0x5f3759df constant in it (only of
historical interest now).

If you want to work with very large (thousands, millions of digits)
integers then gmp library might be faster then the default Python
integer implementation.


--
Akira