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

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?

  • $\begingroup$ 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. $\endgroup$ Commented Feb 2, 2018 at 22:23
  • $\begingroup$ 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. $\endgroup$
    – reddish
    Commented Feb 2, 2018 at 23:29
  • $\begingroup$ 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. $\endgroup$ Commented Feb 3, 2018 at 7:50
  • $\begingroup$ I don't understand the point you are trying to make. $\endgroup$
    – reddish
    Commented Feb 3, 2018 at 21:26
  • $\begingroup$ Unfortunately NumberTheory`IntegerSqrt[] is even slower than the naive approach for the OP's example. $\endgroup$ Commented Mar 4, 2018 at 20:57

2 Answers 2


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]]]]


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;

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.

  • 1
    $\begingroup$ 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. $\endgroup$
    – reddish
    Commented Feb 2, 2018 at 19:21
  • 1
    $\begingroup$ 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 $\endgroup$ Commented Feb 2, 2018 at 20:19
  • 1
    $\begingroup$ 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. $\endgroup$
    – reddish
    Commented Feb 2, 2018 at 21:39
  • 1
    $\begingroup$ It seems GMP uses an integer method due to Paul Zimmermann: gmplib.org/manual/Square-Root-Algorithm.html $\endgroup$ Commented Feb 2, 2018 at 21:57
  • $\begingroup$ It's quite remarkable that those methods show identical performance. $\endgroup$
    – reddish
    Commented Feb 3, 2018 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


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


g == h


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

  • 2
    $\begingroup$ I should have read this more carefully, I failed to notice it uses very similar precision control to what I posted. $\endgroup$ Commented Feb 2, 2018 at 18:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.