Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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
2  
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

5 Answers 5

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.

enter image description here

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.

enter image description here

share|improve this answer
    
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 –  rm -rf 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 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.

share|improve this answer
    
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
1  
Ignore what I said about EvenQ! Need coffee! –  Sal Mangano Jan 21 '12 at 15:30
up vote 5 down vote accepted

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] &;
share|improve this answer
    
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... –  rm -rf Feb 20 '12 at 14:21
    
@R.M you would do better to ask Sasha. –  Mr.Wizard Feb 20 '12 at 15:22
2  
@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
    
@KennyColnago would you please post your complete function as an answer? –  Mr.Wizard Oct 19 '12 at 0:24
    
Speaking of the code you say you got from Sasha: have a look at the definition of NumberTheory`PowersRepresentationsDump`ProbablePerfectSquareQ[] sometime... –  J. M. Oct 20 '12 at 15:46

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?).

share|improve this answer
    
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

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]
share|improve this answer
1  
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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