# Fast integer square-root

I'm looking for the highest-performance method of calculating integer square roots in Mathematica of very big arbitrary-precision numbers.

As an example testcase, I use:

n = 10^1000000 - 3 ^ 2095903


On my Core i7 linux machine, calculating the integer (floor) square root using the straightforward method takes 2.68 seconds:

In[1]:= n = 10^1000000 - 3^2095903;
In[2]:= First@Timing@Floor[Sqrt[n]]
Out[2]= 2.68984


However, the same machine can calculate the integer square root much faster using the GMP library. Here's an example python program doing just that:

from gmpy2 import mpz, isqrt
from time import time

n = mpz(10**1000000 - 3**2095903)

t0 = time()
s = isqrt(n)
t1 = time()

print(t1 - t0)


On the same machine, this is approximately 160 times faster:

\$ python3 benchmark.py
0.01685619354248047


The frustrating thing is that Mathematica already depends on the GMP library for its arbitrary-precision integers, but I cannot seem to make Mathematica use the GMP implementation of the integer square root.

Instead I suspect that its Floor[Sqrt[...]] function falls back on a generic algorithm that must calculate its argument to sufficient precision using arbitrary-precision floating-point approximation.

Is there any way to speed up the calculation of integer square roots, preferably approaching GMP's raw performance?

• this is more a question of how the programming languages store the numbers with data structures. python could calculate using big int down to a regular int representation of the number while mathematica just stores the symbolics and does no actual calculations of what kind of a number 10^.... actually is before you start that clock. – mathreadler Feb 2 '18 at 22:23
• Mathematica will evaluate the 1,000,000 digit integer before assigning it to the variable 'n', as can be established using FullForm[n] or Head[n]. This particular aspect of Mathematica's behavior is pretty predictable. – reddish Feb 2 '18 at 23:29
• But what does evaluate mean? It usually means parse and interpret so that we can store into some data structure. There are quite many ways you can do that symbolically and we might need to help mathematica to tell which kind of answer and precision we would like. – mathreadler Feb 3 '18 at 7:50
• I don't understand the point you are trying to make. – reddish Feb 3 '18 at 21:26
• Unfortunately NumberTheoryIntegerSqrt[] is even slower than the naive approach for the OP's example. – J. M.'s ennui Mar 4 '18 at 20:57

Sqrt is using exact methods in an effort to pull out "small" squares. This is going to take time. A direct approach, as already noted, would do the square root numerically. For purposes of Floor extraction, it suffices to use as precision half the digit size.

floorSqrt[n_Integer] :=
Module[{prec = RealExponent[1.001*n]/2},
Floor[Sqrt[SetPrecision[n, prec]]]]


Test:

n = 10^1000000 - 3^2095903;
AbsoluteTiming[sn1 = Floor[Sqrt[n]];]
AbsoluteTiming[sn2 = floorSqrt[n];]
sn1 === sn2

(* Out[410]= {2.774283, Null}

Out[411]= {0.016654, Null}

Out[412]= True *)


For variety, here is a top-level implementation of an integer-based method. My guess is it is similar to Zmmermann's, but he may well have had some extra efficiencies. The idea is to split the number into an upper and lower part, using a power of 4 for the split size so that we can shift back by a factor of 2. That is, write a = 4^n*b+c. Recursively compute the integer sqrt of b, multiply by 2^n, and use the usual Taylor approximation to get a correction that estimates sqrt(a). Last step is to iteratively repair that estimate. We use integer multiplication, squaring, and the integer Quotient function. The recursion has to bottom out so we provide some base size below which we use the numeric approximation method noted already.

iSqrt[a_, baselen_] := Catch[Module[
{quarterscale = Ceiling[RealExponent[a, 2.]/4], aUpper, aLower,
sqrt, diff},
If[quarterscale < baselen, Throw[Floor[Sqrt[N[a, 2*baselen + 4]]]]];
aUpper = BitShiftRight[a, 2*quarterscale];
aLower = a - BitShiftLeft[aUpper, 2*quarterscale];
sqrt = BitShiftLeft[iSqrt[aUpper, baselen], quarterscale];
sqrt += Quotient[aLower, (2*sqrt)];
diff = a - sqrt^2;
While[Abs[diff] >= 2 sqrt + 1 || diff < 0,
sqrt += Quotient[diff, (2*sqrt)];
diff = a - sqrt^2;
];
sqrt
]]


This is not as fast as the numerical approximation method but it is not bad either. Below is on my slow laptop.

n = 10^1000000 - 3^2095903;
AbsoluteTiming[sn1 = Floor[Sqrt[n]];]
AbsoluteTiming[sn2 = floorSqrt[n];]
AbsoluteTiming[sn3 = iSqrt[n, 10];]
sn1 == sn2 == sn3

(* Out[642]= {7.20206, Null}

Out[643]= {0.0479314, Null}

Out[644]= {0.128822, Null}

Out[645]= True *)


So it's around 2.5 times slower than the numeric code. An internal implementation might do better in terms of avoiding some amount of overhead, and might possibly make up a good chunk from that factor of 2.5.

• I am impressed by the performance of this approach -- it's really in the same ballpark as GMP's integer square root performance (0.022417 seconds on my system for this vs 0.016856 seconds for GMP). The fact that the result passes through an evaluation that depends on precision is slightly unnerving, though. – reddish Feb 2 '18 at 19:21
• It seems from here that the isqrt() in Python also uses finite precision. Check the use of hbc (which is slightly larger than half the bit length) in the main loop for isqrt_fast_python – Daniel Lichtblau Feb 2 '18 at 20:19
• Nice find, curious! I made a C implementation of the same calculation using just GMP functions, and it also executes in 16.8 ms. It could be that GMP also does something similar. – reddish Feb 2 '18 at 21:39
• It seems GMP uses an integer method due to Paul Zimmermann: gmplib.org/manual/Square-Root-Algorithm.html – Daniel Lichtblau Feb 2 '18 at 21:57
• It's quite remarkable that those methods show identical performance. – reddish Feb 3 '18 at 21:30

Not directly addressing the question about GMP library, but you can get a Sqrt much faster starting with an extended precision float.

n = 10^1000000 - 3^2095903;
(g = Floor[Sqrt[n]]) // AbsoluteTiming // First


6.22126

(h = Floor[Sqrt[N[n, Log[10, n]]]]) // AbsoluteTiming // First


0.086576

g == h


True

using this you maybe want a validity check, (h + 1)^2 > n > h^2` takes very little time.

• I should have read this more carefully, I failed to notice it uses very similar precision control to what I posted. – Daniel Lichtblau Feb 2 '18 at 18:32