main(S) :- rho(5, S).

mcd(A, 0, D) :- D is abs(A).
mcd(A, B, D) :- R is A mod B, mcd(B, R, D).

even(X) :- 0 is X mod 2.

pot(_, 0, N, P) :- P is 1 mod N.
pot(X, E, N, P) :-
    F is E div 2,
    pot(X, F, N, Y),
    (even(E) -> P is (Y * Y) mod N
              ; P is (X * ((Y * Y) mod N)) mod N).

g(X, E, N, Y) :- pot(X, E, N, P), Y is (1 + P) mod N.

rho(F, D) :-
    E is 2 ^ (F + 2),
    N is 1 + 2 ^ (2 ^ F),
    S is 3,
    g(S, E, N, T),
    g(T, E, N, L),
    mcd(T - L, N, C),
    rho_1(C, E, N, T, L, D).

rho_1(1, E, N, T, L, D) :-
    g(T, E, N, U),
    g(L, E, N, W),
    g(W, E, N, V),
    mcd(U - V, N, C),
    rho_1(C, E, N, U, V, D).

rho_1(D, _, _, _, _, D).