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