9
$\begingroup$

Goal

The goal is to efficiently check whether a number is a perfect power.

Attempts

It is possible to check whether a number is a perfect power using FactorInteger, as seen in (42251), but that method is very slow, especially for large numbers.

Since the upper bound of my list of interest is around 10^18, I realized that if a number in my list is a perfect power, its power would be at most 59(= Floor[Log2[10^18]]).

Thus I wrote the following code:

powerq[in_] := Evaluate[Or @@ Array[FractionalPart[in^(1/Prime[#])] == 0 &, PrimePi[59]]];
read = OpenRead["numbers"];
write = OpenWrite["powers"];
Do[If[powerq[#], Write[write, #], , Break[]]&@Read[read], Infinity]

(I am using Read and Write because my list has approximately 2 billion numbers (~40 GB) and cannot be stored in RAM; I used FractionalPart instead of IntegerQ because IntegerQ does not work symbolically)

This was slow, so I decided to use Compile

powerqcomp = Compile[{{in, _Integer}},
 Evaluate[
  Or @@ Array[FractionalPart[in^(1/Prime[#])] == 0 &, 
    PrimePi[59]]], CompilationTarget -> "C"]

The problem was that the Compiled function uses machine real, so it sometimes returns fallacious answer, such as:

powerqcomp[2^15]
(* False *)

So, I used a different approach.

If a number $n$ can be expressed as $a^b$ ($a$ and $b$ are positive integers and $b > 1$), then $b \leq \log_{2}{n}$. Then, I could check whether there exists $a$ such that $n = a^b$ for at least one $b$.

Since my list consists of large numbers, trying all $a$ from $1$ to $n^{(1/b)}$ is out of the question. Perhaps I could use BinarySearch, but then I would need to generate Range[n^(1/b)]; I would need to make my own function that performs a binary search.

This is the code I made that uses binary search. It does work, but it is still too slow (and gives CompiledFunction::cfn errors sometimes).

powerqcomp2 = 
 Compile[{{in, _Integer}}, 
  Module[{up = 10^9, down = 1, 
    blist = Select[{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 
       43, 47, 53, 59}, # <= Log[2, in] &]},
   Catch[Do[
     While[True,
      Which[#^b == in, Throw[True],
         #^b > in, up = #,
         True, down = #
         ] &[Ceiling[(up + down)/2]];
      If[up - down == 1, Break[]]
      ];
     up = 10^9;
     down = 1,
     {b, blist}]; False]]]

Question

How could one create code that efficiently checks whether a number is a perfect power?

$\endgroup$
  • $\begingroup$ Have you seen this? $\endgroup$ – J. M. is away Jul 23 '16 at 20:43
  • $\begingroup$ @J.M. I read the article, but it seems like Mathematica is not designed to perform operations in there. E.g. separating the mantissa and the exponent from a floating-point number. MantissaExponent calculates the mantissa and exponent, which is slower than simple extraction. Plus, bitwise operations are not permitted on floating-point reals in Mma... How could that algorithm be implemented in Mma? $\endgroup$ – JungHwan Min Jul 24 '16 at 3:18
  • $\begingroup$ As there are only a million cubes in the range of interest (and less of higher powers) perhaps you should build a list of all these and only check for exact squares. $\endgroup$ – mikado Jul 24 '16 at 10:02
6
$\begingroup$

I offer the following as a fast way of testing cubic and higher powers

primes = Select[Range[59], PrimeQ]

Get a list of all the relevant powers up to a specified limit

list[nmax_] := Sort[Flatten[Table[Range[2, Floor[nmax^(1/p)]]^p, {p, Drop[primes, 1]}]]];

For example

list[1000]

(* {8, 27, 32, 64, 125, 128, 216, 243, 343, 512, 729, 1000} *)

Define an AssociationMap and Lookup function

 asc = AssociationMap[True&, list[10^19]];
 exactpower1 = Lookup[asc, #, False] &;

Test it on a million random integers in under a second

AbsoluteTiming[Tally[exactpower1 /@ RandomInteger[{10^18, 2*10^18}, {10^6}]]]
(*  {0.678293, {{False, 1000000}}} *)

ADDED

In combination with a check for higher order powers, the following can rapidly test for possible squares. It uses the logic described in https://gmplib.org/manual/Perfect-Square-Algorithm.html to make tests modulo various numbers before falling back to a full test involving a square root.

exactsquarefull[n_] := IntegerQ[Sqrt[n]]

 maketest[n_] := Module[{list, asc},
  list = Union[Mod[Range[n]^2, n]];
  asc = AssociationMap[True &, list];
  possiblesquare[n] = Lookup[asc, Mod[#, n], False] &]

 maketest /@ {256, 9, 5, 7, 13, 17, 97};

 exactsquare[n_] := 
 possiblesquare[256][n] && possiblesquare[9][n] && 
  possiblesquare[5][n] && possiblesquare[7][n] && 
  possiblesquare[13][n] && possiblesquare[17][n] && 
  possiblesquare[97][n] && exactsquarefull[n]

This matches the expected result on all numbers up to 10^6

 {exactsquare[#] == exactsquarefull[#]} & /@ Range[1000000] // Union

 (* {{True}} *)

and checks 10^6 large numbers in under 3 seconds

 AbsoluteTiming[Tally[exactsquare /@ RandomInteger[{10^18, 2*10^18}, {1000000}]]]
 (* {2.76558, {{False, 1000000}}} *)
$\endgroup$
  • $\begingroup$ Thanks! Works like a charm. However, I think there's a typo in your first definition of asc. Shouldn't it have True & not True? $\endgroup$ – JungHwan Min Jul 24 '16 at 18:25
  • $\begingroup$ @JHM quite correct $\endgroup$ – mikado Jul 24 '16 at 18:50
  • $\begingroup$ Just to note, the functions exactpower1 and exactsquare both return False when a non-Integer argument is used. I fixed it by modifying them so that they only take integers as arguments (i.e. exactpower1[n_Integer] := ... and exactsquare[n_Integer] := ... ) $\endgroup$ – JungHwan Min Jul 24 '16 at 21:56
5
$\begingroup$

There is this way:

SetAttributes[test, Listable]
test[n_] := FirstPosition[Reverse[#[[2 ;; Ceiling[Length[#]/2]]] &[
                  Divisors[n]]], _?(IntegerQ[Log[#, n]] &), 0] =!= 0
Pick[#, test[#]] &[Range[1000]]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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