Perfect Powers

/* perfect_power(+N, ?X, ?Y) is true if N=X^Y, where N>=2, X>=2, Y>=2 and    */
/*   X is minimal.                                                           */ 
/* e.g. perfect_power(387420489, 3, 18).                                     */
perfect_power(N, X, Y):-
  prime_power(N, Z, P),
  perfect_power_1(Z, P, X, Y).
  
perfect_power_1(Z, P, X, Y):-
  perfect_power(Z, X, Q), !,
  Y is P * Q.
perfect_power_1(X, Y, X, Y).

/* prime_power(+N, ?X, ?P) is true if N=X^P, where N>=2, X>=2 and P>=2, and  */
/*   P is prime.                                                             */ 
/* e.g. prime_power(387420489, 19683, 2).                                    */
prime_power(N, X, P):-
  N >= 2,
  log2(N, LogN),
  Ps=[2,3,5,7,11,13,17,19,23,29], % Valid for N < 2^31
  prime_power_1(Ps, N, LogN, X, P).

prime_power_1([P|_], N, LogN, X, P):-
  P =< LogN,
  Temp is LogN // P + 1,
  power(2, Temp, H),
  binary_search(2, H, N, P, X), !.
prime_power_1([P|Ps], N, LogN, X, Q):-
  P =< LogN,
  prime_power_1(Ps, N, LogN, X, Q).

/* log2(+N, ?LogN) is true if LogN=trunc(log2(N)), that is, one less than    */
/*   the number of significant binary digits in N.                           */
/* e.g. log2(123456789,26).                                                  */
log2(N, LogN):-
  N > 0,
  log2(-1, N, LogN).

log2(I, 0, I):-!.
log2(I, N, LogN):-
  I1 is I + 1,
  N1 is N // 2,
  log2(I1, N1, LogN).

/* binary_search(+Low, +High, +N, +P, ?X) is true if N=X^P for some value of */
/*   X between Low and High inclusive.                                       */
/* e.g. binary_search(2, 4, 1162261467, 19, 3).                              */
binary_search(Low, High, N, P, X):-
  Low1 is Low - 1,
  binary_search_1(Low1, High, High, N, P, X).

binary_search_1(Low, High, Highest, N, P, X):-
  Low + 1 < High,
  Mid is (Low + High) // 2,
  aln(P * ln(Mid)) =< 2147483647, !, % Avoid integer overflow
  binary_search_2(Mid, P, N, Low, Low1, High, High1),
  binary_search_1(Low1, High1, Highest, N, P, X).
binary_search_1(_, High, Highest, N, P, High):-
  High < Highest,
  power(High, P, N).

binary_search_2(Mid, P, N, _, Mid, High, High):-
  power(Mid, P, Temp),
  Temp < N, !.
binary_search_2(Mid, _, _, Low, Low, _, Mid).

/* power(+X, +N, ?Y) is true if Y is the result of raising X to the power N. */
/* e.g. power(3, 19, 1162261467).                                            */
power(X, N, Y):-
  N >= 0,
  power_1(N, X, 1, Y).

power_1(0, _, Y, Y):-!.
power_1(1, X, Y0, Y):-!,
  Y is Y0 * X.
power_1(N, X, Y0, Y):-
  N mod 2 =:= 1, !,
  N1 is N // 2,
  X1 is X * X,
  Y1 is Y0 * X,
  power_1(N1, X1, Y1, Y).
power_1(N, X, Y0, Y):-
  N1 is N // 2,
  X1 is X * X,
  power_1(N1, X1, Y0, Y).

LPA Index     Home Page