Efficiently checking whether a number is a perfect power

9

1

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?

JungHwan Min

Posted 2016-07-23T20:28:58.170

Reputation: 4 484

Have you seen this?

– J. M.'s ennui – 2016-07-23T20:43:34.640

@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 – 2016-07-24T03:18:46.243

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 – 2016-07-24T10:02:55.823

Answers

6

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

mikado

Posted 2016-07-23T20:28:58.170

Reputation: 11 100

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 – 2016-07-24T18:25:03.343

@JHM quite correct – mikado – 2016-07-24T18:50:51.500

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 – 2016-07-24T21:56:53.960

5

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

Coolwater

Posted 2016-07-23T20:28:58.170

Reputation: 17 733