Faster PowersRepresentation using IntegerPartitions

Question

Earlier I posted a question about taking fast integer square roots of known integer perfect squares.

The reason this came up is I was trying to find a faster way of mimicking the PowersRepresentations[] function in Mathematica.

Using IntegerPartitions[] I can speed up this by quite a bit - but then I needed to take the square roots of the output from IntegerPartitions - which ends up costing me about a factor of two or three in speed. Though still a lot faster than using PowersRepresentations[] directly.

For example:

n = 8174;
parts = 6;

Length@PowersRepresentations[n, parts, 2] // AbsoluteTiming

(* {1.1861, 31934} *)

Length@IntegerPartitions[n, {parts}, 
   Range[0, IntegerPart[Sqrt[n]]]^2] // AbsoluteTiming

(* {0.021446, 31934} *)

Length@Round[
   Sqrt[N[
     IntegerPartitions[n, {parts}, 
      Range[0, IntegerPart[Sqrt[n]]]^2]]]] // AbsoluteTiming

(* {0.039502, 31934} *)

This might be a silly wish - but is there a way of cutting out the middle man so to speak. I'd like the output of IntegerPartitions[] to be the integers themselves and not their squares. This seems unlikely but does anyone know a workaround or am I destined to have to take the square roots after the fact? Having said this I would think a workaround would be no faster anyway. I'll delete the question if it is asking the impossible.

(Oh, and one thing that puzzles me - why does PowersRepresentations[] run so much faster on a repeat run even if I clear the system cache?)

Answer 1

This is an extended comment, with pictures.

PowersRepresentations is implemented in top-level code, so you may take inspiration from its existing code and modify it to output what you want. In particular, PowersRepresentations[.. , .., 2] calls the internal NumberTheory`PowersRepresentationsDump`squareRepresentation, which is also implemented in top-level code. Perhaps you may be able to use that code as inspiration and modify it to output what you want.

(* The single call below is needed to force the loading of definitions *)
PowersRepresentations[10, 2, 2]; 

(* remove the ReadProtected attribute *)
ClearAttributes[PowersRepresentations, ReadProtected]

(* retrieve the definitions of the function, if any *)
??PowersRepresentations

Rooting through the above you will notice the following, where I've stripped the NumberTheory`PowersRepresentationsDump` for readability:

PowersRepresentations[n_Integer?Positive, d_Integer?Positive, 2] := 
 Block[{res}, 
  res = NumberTheory`PowersRepresentationsDump`squareRepresentation[d, n]; 
  res /; FreeQ[res, $Failed]
]

This shows that PowersRepresentations[n, d, 2], which seems to be the application of interest to you, reduces to a call to NumberTheory`PowersRepresentationsDump`squareRepresentation[d, n] (notice the switch in the order of the arguments).

Repeating the same steps, you can dig further into the definition of e.g. squareRepresentation, TaxiCabNumbers, and other functions that seem to be the actual implementations of the underlying algorithms.