# Efficiently checking whether a number is a perfect power

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]]];
write = OpenWrite["powers"];


(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?

• Have you seen this? – J. M.'s torpor Jul 23 '16 at 20:43
• @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? – JungHwan Min Jul 24 '16 at 3:18
• 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. – mikado Jul 24 '16 at 10:02

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}}} *)


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}}} *)

• 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? – JungHwan Min Jul 24 '16 at 18:25
• @JHM quite correct – mikado Jul 24 '16 at 18:50
• 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] := ... ) – JungHwan Min Jul 24 '16 at 21:56

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