# Fastest square number test

What is the fastest possible square number test in Mathematica 7, both for machine size and big integers?

I presume in version 8 the fastest will be a dedicated C LibraryLink function.

-
Even in version 7, you have MathLink. I don't know how much the MathLink overhead is, but on the local machine MathLink should be very fast as it uses memory mapped files for communication. – Szabolcs Jan 21 '12 at 14:45
@Szabolcs; you're right; I was curious what would be best in native Mathematica, but if someone cares to create a robust MathLink executable and share it here I certainly won't turn it down! – Mr.Wizard Jan 21 '12 at 15:01
The overhead of messaging involved in MatheLink would kill it. Possibly Mathematica 8's ability to load a DLL would help if you used C. But consider if your application really warrants this compared to what I offered below. The speed up you may get is probably not worth the effort. – Sal Mangano Jan 21 '12 at 16:04

More info as requested by @Mr.Wizard. For $n$ below the $\approx 2*10^9$ limit, Compile gives the fastest solutions. For larger $n$, Sasha used JacobiSymbol with four primes 13, 19, 17, and 23, before resorting to the expensive IntegerQ[Sqrt[n]]. The number of ambiguous cases where JacobiSymbol[n,p]=0 decreases as the size of the prime $p$ increases. So using larger $p$ helps filter out more candidates before Sqrt must be called. Similarly, using more primes filters more candidates. However, the computation of JacobiSymbol slows as the number of and size of $p$ increases (no free lunch). As a rough balance, I used SquareQ08.

SquareQ08[n_] :=
JacobiSymbol[n, 541] =!= -1 && JacobiSymbol[n, 547] =!= -1 &&
JacobiSymbol[n, 557] =!= -1 && JacobiSymbol[n, 563] =!= -1 &&
JacobiSymbol[n, 569] =!= -1 && JacobiSymbol[n, 647] =!= -1 &&
JacobiSymbol[n, 653] =!= -1 && JacobiSymbol[n, 659] =!= -1 &&
IntegerQ[Sqrt[n]]
SetAttributes[SquareQ08, Listable]

-
Would it be advantageous to select the number of primes to use based on the size of the integer? – s0rce Oct 24 '12 at 20:31
@s0rce My rough tests showed that about 50% of candidate $n$ were filtered out for each JacobiSymbol test, and that this percentage was independent of the size of $n$. Your mileage may vary... – KennyColnago Oct 24 '12 at 23:55
@KennyColnago I just found this thread, and tried it myself. I tried isSq2 as well as SquareQ08 on sets of 100 random integers in various ranges (10^8 to 10^9, 10^9 to 10^10, 10^10 to 10^11, and 10^11 to 10^12) with RepeatedTiming. In each case, the results were remarkably consistent - isSq2 was about 7.4 x 10^{-5} seconds and SquareQ08 was about 4.5 x 10^{-4} seconds. This is on Mathematica 10.3; I wonder what version you were using and if the performance of isQq2 was somehow improved. – rogerl Dec 24 '15 at 17:09
@rogerl I don't recall which version was used. On version 10.3, I find isSq2 about 6 to 10 times faster than SquareQ08, which is written for integers of any size. I do not understand how isSq2, which is compiled, works with integers larger than 2 billion. Can you explain? – KennyColnago Dec 25 '15 at 1:09
@KennyColnago Well, my machine is a 64-bit machine, so my reading of the tea leaves is that Compile will work with 64-bit integers (or about $10^18$). Is that right? And, just to be clear, your result regarding isQq1 and SquareQ08 that you give in your comment above is at variance with your original answer, right? Do you know what's going on, or am I just confused? – rogerl Dec 25 '15 at 1:44

I voted for all three previous answer because they all taught me something. However they, being Compile solutions, are not helpful with big integers.

At least on my system, Sal Mangano's code appears reducible to this without loss of speed:

isSq2 = Compile[n, Floor@# == # & @ Sqrt @ n];


For big integers between about 2*10^9 and 2*10^11 I am currently using this code from Sasha:

SquareQ =
JacobiSymbol[#, 13] =!= -1 &&
JacobiSymbol[#, 19] =!= -1 &&
JacobiSymbol[#, 17] =!= -1 &&
JacobiSymbol[#, 23] =!= -1 &&
IntegerQ@Sqrt@# &;


For integers larger than that I am using code (modified) from Daniel Lichtblau:

SquareQ2 = # == Round@# & @ Sqrt @ N[#, Log[10, #] + $MachinePrecision] &;  - Could you add references or a short explanation as to how JacobiSymbol comes into play here. I don't know the underlying workings of this and I'd like to know... – R. M. Feb 20 '12 at 14:21 @R.M you would do better to ask Sasha. – Mr.Wizard Feb 20 '12 at 15:22 @rm -rf JacobiSymbol[n,p]=1 for any odd prime$p$when$n$is square. JacobiSymbol[n,p]=1 for some odd primes$p$when$n$is not square. JacobiSymbol[n,p]=0 for any odd prime$p$that divides$n$, where$n$may or may not be square. JacobiSymbol[n,p]=-1 otherwise. Hence, JacobiSymbol[n,p]=-1 means$n$is not square. Evaluating JacobiSymbol[n,p] is much faster than Sqrt[n]. Use JacobiSymbol to filter candidate$n$before passing to Sqrt[n]. I tested with 8$p$near 541 on$n>10^{11}$and found it faster than SquareQ by @Sasha and SquareQ2 by @Daniel Lichtblau. – KennyColnago Oct 18 '12 at 16:41 Also, this is all academic. No serious application would use this method. Mathematica is missing integer remainder square root function$(N-\lfloor \sqrt{N} \rfloor^2)$which would solve this problem. The computation of that function is more efficient than multiplication of two numbers of the size of$N\$. – mhp Nov 6 '14 at 8:31
Not really. It took me several weeks to develop the algorithm, just to find out that someone had already beat me to it. See: hal.inria.fr/inria-00072854/PDF/RR-3805.pdf – mhp Nov 6 '14 at 8:42

I don't think there are any built-in functions for this but the following is probably fast enough for most purposes.

isSq = Compile[{{n, _Integer}}, With[{test = Sqrt[n]},
Floor[test] == test]];


Does 1 million integers in under a second.

Timing[Table[isSq[i], {i, 1, 1000000}]][[1]]
(*
0.76195
*)


This is under 2 orders of magnitude faster than the un-compiled equivalent, by the way.

-
The Sqrt operation is comparatively expensive; I believe one should filter out numbers that are destined to fail with a less expensive test before progressing to Sqrt. – Mr.Wizard Jan 21 '12 at 15:08
Adding a test for EvenQ gives a slight speed up but not much. If the numbers are very large it probably helps more. Testing for ending in 0,1,4,6,9, or 25 is possible but I have not tried it and seem like overkill. – Sal Mangano Jan 21 '12 at 15:23
Ignore what I said about EvenQ! Need coffee! – Sal Mangano Jan 21 '12 at 15:30

## Update

Sorry for my ignorance not taking into account that the question specifically asked for a Mathematica 7 solution. I updated the complete post.

## Mathematica 7

In Mathematica 7 we don't have the option the compile code into a C-library which includes the thread parallelization which can be turned on when using RuntimeAttributes->Listable and Parallelization->True. Therefore, acl's solution will not run in Mathematica 7 because the RuntimeAttributes option for Compile was introduced in version 8.

This leaves the possibility to not compile the used function and make it a normal Mathematica function where you can set the attribute Listable. I tried this, but it was horrible slow.

After a bit of research I found a nice solution which uses some number-properties in base 16. Since (at least in V7) it seems somewhat hard to return lists of True|False, I use 0 and 1 where 0 means no square.

fPat = Compile[{{numbers, _Integer, 1}},
With[{l = Length[numbers]},
Module[{n = 0, i = 0, h = 0, test = 0.0, result = Table[0, {l}]},
For[i = 1, i <= l, ++i,
n = numbers[[i]];
h = BitAnd[15, n];
If[h > 9, Continue[]];
If[h != 2 && h != 3 && h != 5 && h != 6 && h != 7 && h != 8,
test = Sqrt[n];
result[[i]] = test == Floor[test];
];
];
result
]
]
];


Comparing this with the almost one-liner of Sal gives

data = Table[i, {i, 1, 10^6}];

fSal = Compile[{{n, _Integer}},
With[{test = Sqrt[n]}, Floor[test] == test]];

BarChart[{Timing[fSal /@ data][[1]], Timing[fPat[data]][[1]]
}, ChartLabels -> {"Sal Mangano", "Patrick V7"},
ChartStyle -> {Gray, Green}]


I leave it to you to decide whether such a C-like programming style is worth the small speed up.

## Mathematica 8

The fastest way (using Mathematica only) I know is to compile a C-library and process all data in parallel. Since most computers these days have at least 2 cores, this gives a boost. In Mathematica 8 the compilation to a C-library does not copy the data when it is called.

To make the computation parallel you have to use the Parallization option and the compiled function must be Listable. If you are sure of your input-data, you can additionally switch off most of the data-checks by using RuntimeOptions set to "Speed".

Update I include here the parallelized version of the Mathematica 7 code above:

data = Table[i, {i, 1, 10^6}];

fSal = Compile[{{n, _Integer}},
With[{test = Sqrt[n]}, Floor[test] == test]];
fAcl = Compile[{{n, _Integer}},
With[{test = Sqrt[n]}, Floor[test] == test],
RuntimeAttributes -> {Listable}];
fPat = Compile[{{n, _Integer}},
With[{test = Sqrt[n]}, Floor[test] == test],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"];

fPat2 = Compile[{{numbers, _Integer, 1}},
With[{l = Length[numbers]},
Module[{n = 0, i = 0, h = 0, test = 0.0, result = Table[0, {l}]},
For[i = 1, i <= l, ++i,
n = numbers[[i]];
h = BitAnd[15, n];
If[h > 9, Continue[]];
If[h != 2 && h != 3 && h != 5 && h != 6 && h != 7 && h != 8,
test = Sqrt[n];
result[[i]] = test == Floor[test];
];
];
result
]
], CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"
];

BarChart[{Timing[fSal /@ data][[1]], Timing[fAcl[data]][[1]],
Timing[fPat[data]][[1]],
Timing[fPat2[data]][[1]]},
ChartLabels -> {"Sal Mangano", "acl", "Patrick",
"Patrick V7 parallel"},
ChartStyle -> {Gray, Gray, Darker[Green], Green}]


The results here come from my MacBook in battery-save mode which has 2 Intel cores. The disadvantage is that you need a C-compiler installed on your system which is most likely not true for the majority of Mathematica users.

-
Here's some background info on the method employed in fPat and fPat2 for others. I was almost done implementing it, but I see you beat me to it :) I believe that will be the fastest of all approaches – R. M. Jan 21 '12 at 21:14
@R.M, sorry ;-) I was ashamed that I didn't read the question carefully and I felt that I should correct this. – halirutan Jan 21 '12 at 21:52
btw, really nice is that in a compiled function result[[i]] = test == Floor[test]; works as expected for integers and gives a speed-up compared to its if-then-else counterpart. – halirutan Jan 21 '12 at 21:53
Very nice, +1. Mind you, @Mr.W will probably not consider this mathematica code, though: stackoverflow.com/a/8891064/559318 [see comments] :) – acl Jan 21 '12 at 23:40
At least I didn't use Goto ;-) – halirutan Jan 21 '12 at 23:50

I am not sure how to speed up each comparison (as in, I spent half an hour trying different things and didn't manage to), but making the compiled function listable speeds things up quite a bit.

If isSq is the direct implementation that Sal gave, simply make it listable and compare:

isSqL = Compile[
{{n, _Integer}}, With[{test = Sqrt[n]}, Floor[test] == test],
RuntimeAttributes -> {Listable}
];


and then compare:

Timing[Table[isSq[i], {i, 1, 10^6}]]; // Timing
isSq /@ Range[1, 10^6]; // Timing
isSqL[Range[1, 10^6]]; // Timing
(*
{0.697799, Null}
{0.545856, Null}
{0.150171, Null}
*)


ie, it's 3-4 times faster.

What makes you say Sqrt is expensive? (ie, compared to what?).

-
True, but this does not speed up the sq test per se it is just a faster way to apply the function across a list than Table is. Good idea though! – Sal Mangano Jan 21 '12 at 16:25
@Sal That's true, but I can't get the test faster. Looking at CompilePrint[isSq], it's hard to see anything more costly thatn Sqrt`. – acl Jan 21 '12 at 16:34
Perhaps, you'd like to use AbsoluteTiming for these tests, especially if stuff is running in parallel. See the "How can I improve the speed of eigenvalue decompositions for large matrices?" thread. – user21 Jan 21 '12 at 17:49
@ruebenko that is a good point (in this case, though, they give the same answer, so I'll leave this the way it is--but your point about running them in parallel stands) – acl Jan 21 '12 at 20:17