#### Fermat-Pell Equations

This code requires the Multi-Precision Integer Arithmetic package.

```%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% x^2-Dy^2=1 using the English method %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% english(D, J, XYs) is true if the elements [X,Y] of XYs are the first J
%   solutions of x^2-Dy^2=1, where D > 1 and non-square.
% e.g. english(`15`, 3, XYs) gives XYs=[[`4`,`1`],[`31`,`8`],[`244`,`63`]]
english(D0, J, XYs):-
string_bigint(D0, D),
less_than([1], D),
square_root(D, SqrtD),
multiply(SqrtD, SqrtD, D1),
less_than(D1, D),
english_1(J, D, SqrtD, [1], [0], [1], [], XYs0),
pairs_strings(XYs0, XYs).

% Finds the list of solutions.
% e.g. english_1(2, [15], [3], [1], [0], [1], [], XYs) gives
%   XYs=[[[4],[1]],[[31],[8]]]
english_1(0, _, _, _, _, _, XYs, XYs):-!.
english_1(J, D, SqrtD, P0, Q0, K0, XYs0, XYs):-
english_3(SqrtD, P0, Q0, K0, R0, Q1), % Q1 = (P0 + Q0*R0) div K0
multiply(P0, R0, P0R0),
multiply(Q0, D, Q0D),
divide(P0R0Q0D, K0, P1, [0]),         % P1 = (P0*R0 + Q0*D) div K0
multiply(R0, R0, R0R0),
subtract(D, R0R0, DR0R0),
divide(DR0R0, K0, K1, [0]),           % K1 = (D - R0*R0) div K0
english_3(SqrtD, P1, Q1, K1, R1, Q2), % Q2 = (P1 + Q1*R1) div K1
multiply(P1, R1, P1R1),
multiply(Q1, D, Q1D),
divide(P1R1Q1D, K1, P2, [0]),         % P2 = (P1*R1 + Q1*D) div K1
multiply(R1, R1, R1R1),
subtract(D, R1R1, DR1R1),
divide(DR1R1, K1, K2, [0]),           % K2 = (D - R1*R1) div K1
english_2(J, D, SqrtD, P2, Q2, K2, XYs0, XYs).

% If K=1 then a solution [P,Q] has been found.
english_2(J, D, SqrtD, P, Q, [1], XYs0, [[P,Q]|XYs]):-!, % Solution found
J1 is J - 1,
english_1(J1, D, SqrtD, P, Q, [1], XYs0, XYs).
english_2(J, D, SqrtD, P, Q, K, XYs0, XYs):-             % Solution not found
english_1(J, D, SqrtD, P, Q, K, XYs0, XYs).

% Finds the largest integer R such that P+QR is divisible by K, and R is less
%   than sqrt(D), and also returns (P+QR) div K.
% e.g. english_3([3], [119], [33], [4], R, Q1) gives R=[1], Q1=[38]
english_3(R, P, Q, K, R, Q1):-
multiply(Q, R, QR),
divide(PQR, K, Q1, [0]), !.   % P+QR is divisible by K
english_3(I, P, Q, K, R, Q1):-  % P+QR is not divisible by K
subtract(I, [1], I1),
english_3(I1, P, Q, K, R, Q1).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% x^2-Dy^2=1 using continued fractions %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% plus_one(D, J, XYs) is true if the elements [X,Y] of XYs are the first J
%   solutions of x^2-Dy^2=1, where D > 1 and non-square.
% e.g. plus_one(`7`, 4, XYs) gives
%   XYs=[[`8`,`3`],[`127`,`48`],[`2024`,`765`],[`32257`,`12192`]]
% e.g. plus_one(`661`, 1, XYs) gives
%   XYs=[[`16421658242965910275055840472270471049`,
%         `638728478116949861246791167518480580`]]
plus_one(D0, J, ABs):-
J > 0,
string_bigint(D0, D),
plus_one_1(D, J, XYs),
pairs_strings(XYs, ABs).

% This is similar to plus_one/3 but uses bigints, not strings.
plus_one_1(D, J, [[P,Q]|XYs]):-
J > 0,
cont_frac(D, Period, As),
even_odd(Period, As, As1),
convergents(As1, PQs),
last(PQs, [P,Q]),
plus_one_2(J, D, P, Q, P, Q, XYs).

plus_one_2(1, _, _, _, _, _, []):-!.
plus_one_2(J, D, P, Q, X, Y, [[X1,Y1]|XYs]):-
multiply(P, X, PX),
multiply(Q, D, QD),
multiply(QD, Y, QDY),
multiply(P, Y, PY),
multiply(Q, X, QX),
J1 is J - 1,
plus_one_2(J1, D, P, Q, X1, Y1, XYs).

% Converts a list of pairs of bigints [X,Y] to a list of pairs of
%   strings [A,B]
pairs_strings([], []).
pairs_strings([[X,Y]|XYs], [[A,B]|ABs]):-
bigint_string(X, A),
bigint_string(Y, B),
pairs_strings(XYs, ABs).

%%%%%%%%%%%%%%%
% x^2-Dy^2=-1 %
%%%%%%%%%%%%%%%

% minus_one(D, J, XYs) is true if the elements [X,Y] of XYs are the first J
%   solutions of x^2-Dy^2=-1, where D > 1 and non-square.
% e.g. minus_one(`5`, 4, XYs) gives
%   XYs=[[`2`,`1`],[`38`,`17`],[`682`,`305`],[`12238`,`5473`]]
% e.g. minus_one(`3001`, 1, XYs) gives
%   XYs=[[`4758959166078874909404390921598306320`,
%         `86871832084583873601621533570831549`]]
minus_one(D0, J, ABs):-
J > 0,
string_bigint(D0, D),
minus_one_1(D, J, XYs),
pairs_strings(XYs, ABs).

% This is similar to minus_one/3 but uses bigints, not strings.
minus_one_1(D, J, [[P,Q]|XYs]):-
J > 0,
cont_frac(D, Period, As),
Period mod 2 =:= 1,
convergents(As, PQs),
last(PQs, [P,Q]),
minus_one_2(J, D, P, Q, P, Q, XYs).

minus_one_2(1,_, _, _, _, _, []):-!.
minus_one_2(J, D, P, Q, X, Y, [[X2,Y2]|XYs]):-
multiply(P, X, PX),
multiply(Q, D, QD),
multiply(QD, Y, QDY),
multiply(P, Y, PY),
multiply(Q, X, QX),
multiply(P, X1, PX1),
multiply(QD, Y1, QDY1),
multiply(P, Y1, PY1),
multiply(Q, X1, QX1),
J1 is J - 1,
minus_one_2(J1, D, P, Q, X2, Y2, XYs).

% cont_frac(D, Period, As) is true if As is the continued fraction
%   corresponding to sqrt(D) up to the term preceding the term which is equal
%   to 2*floor(sqrt(D)), at which point cycling starts with a period of Period.
% e.g. cont_frac([97], Period, As) gives
%   Period=11, As=[[9],[1],[5],[1],[1],[1],[1],[1],[1],[5],[1]]
cont_frac(D, Period, As):-
less_than([1], D),
square_root(D, SqrtD),
multiply(SqrtD, SqrtD, D1),
less_than(D1, D),
cont_frac(D, SqrtD, TwoSqrtD, [0], [1], -1, Period, [[0]|As]).

% cont_frac(D, SqrtD, TwoSqrtD, [0], [1], -1, Period, As) is true if As is the
%   continued fraction corresponding to sqrt(D) up to the term preceding the
%   term which is equal to TwoSqrtD=2*SqrtD where SqrtD=floor(sqrt(D)), at
%   which point cycling with a period equal to Period starts.
% e.g. cont_frac([97], [9], [18], [0], [1], -1, Period, As) gives
%   Period=11, As=[[0],[9],[1],[5],[1],[1],[1],[1],[1],[1],[5],[1]]
cont_frac(D, SqrtD, TwoSqrtD, U0, V0, Period0, Period, [A|As]):-
multiply(U0, U0, T1),
subtract(D, T1, T2),
divide(T2, V0, V1, [0]),   % V1=(D-U0*U0) div V0
divide(T3, V1, A, _),      % A=(SqrtD+U0) div V1
less_than(A, TwoSqrtD), !, % A<2*SqrtD
multiply(A, V1, T4),
subtract(T4, U0, U1),      % U1=A*V1-U0
Period1 is Period0 + 1,
cont_frac(D, SqrtD, TwoSqrtD, U1, V1, Period1, Period, As).
cont_frac(_, _, _, _, _, Period, Period, []).

% convergents(As, PQs) is true if the elements [P,Q] of PQs are successive
%   convergents of the continued fraction, As.
% e.g. convergents([[4],[1],[2],[4],[2],[1]], PQs) gives
%   PQs=[[[4],[1]],[[5],[1]],[[14],[3]],[[61],[13]],[[136],[29]],[[197],[42]]]
convergents([A], [[A,[1]]]):-!.
convergents([A0,A1|As], [[A0,[1]],[P1,A1]|PQs]):-
multiply(A0, A1, T1),
convergents_1(As, A0, P1, [1], A1, PQs).

convergents_1([], _, _, _, _, []).
convergents_1([A|As], P0, P1, Q0, Q1, [[P,Q]|PQs]):-
multiply(A, P1, T1),
multiply(A, Q1, T2),
convergents_1(As, P1, P, Q1, Q, PQs).

% If Period is even, returns [A1,A2,A3,...]
% If Period is odd, returns [A1,A2,A3,...,2*A1,A2,A3,...]
% e.g. even_odd(5, [[3],[1],[1],[1],[1]], As) gives
%   As=[[3],[1],[1],[1],[1],[6],[1],[1],[1],[1]]
even_odd(Period, As, As):-
Period mod 2 =:= 0, !.
even_odd(_, [A|As0], As):-
append([A|As0], [TwoA|As0], As).

%%%%%%%%%%%%%%
% x^2-Dy^2=N %
%%%%%%%%%%%%%%

% plus_n(D, N, J, XYs) is true if the elements [X,Y] of XYs are solutions of
%   x^2-Dy^2=N, where D > 1, N > 0 and N < sqrt(D).  Solutions will be obtained
%   from the first J solutions of x^2-Dy^2=1 and all fundamental solutions of
%   x^2-Dy^2=N.
% e.g. plus_n(`157`, `12`, 0, XYs) gives
%   XYs=[[`13`,`1`],[`10663`,`851`],[`579160`,`46222`]]
% e.g. plus_n(`157`, `12`, 1, XYs) gives
%   XYs=[[`13`,`1`],[`10663`,`851`],[`579160`,`46222`],
%        [`483790960`,`38610722`],
%        [`26277068347`,`2097138361`],
%        [`21950079635497`,`1751807067011`],
%        [`1192216867392577`,`95149264530709`],
%        [`995897062658343427`,`79481238398745359`],
%        [`54092071464191542720`,`4317017278925659678`]]
plus_n(D0, `1`, J, XYs):-!,
plus_one(D0, J, XYs).
plus_n(D0, N0, J, XYs):-
J >= 0,
string_bigint(D0, D),
string_bigint(N0, N),
plus_n_1(D, N, J, XYs0),
insertion_sort(XYs0, XYs1),
pairs_strings(XYs1, XYs).

% This is similar to plus_n/4 but uses bigints, not strings.
plus_n_1(D, N, 0, XYs):-!,       % J=0
plus_n_fundamental(D, N, XYs). % Fundamental solutions of x^2-Dy^2=N
plus_n_1(D, N, J, XYs):-
plus_one_1(D, J, RSs),         % J solutions of x^2-Dy^2=1
plus_n_fundamental(D, N, PQs), % Fundamental solutions of x^2-Dy^2=N
findall(XY, one_solution(PQs, RSs, D, XY), XYs0),
append(PQs, XYs0, XYs).

% plus_n_fundamental(D, N, XYs) is true if the elements [X,Y] of XYs are the
%   fundamental solutions of x^2-Dy^2=N, where N < sqrt(D).
% e.g. plus_n_fundamental([88], [9], XYs) gives
%   XYs=[[[19],[2]],[[47],[5]],[[3],[0]]]
plus_n_fundamental(D, N, XYs):-
multiply(N, N, N2),
less_than(N2, D),
cont_frac(D, _, As),
convergents(As, XYs0),
plus_n_fundamental_1(XYs0, D, N, XYs1),           % Basic solutions
square_root(N, SqrtN),
plus_n_fundamental_2([2], SqrtN, D, N, [], XYs2), % Additional solutions
append(XYs1, XYs2, XYs).

% Selects those convergents which will give a solution of x^2-Dy^2=N.
plus_n_fundamental_1([], _, _, []).
plus_n_fundamental_1([[X,Y]|XYs0], D, N, [[X,Y]|XYs]):-
multiply(X, X, X2),
multiply(Y, Y, Y2),
multiply(D, Y2, DY2),
subtract(X2, DY2, N), !, % X*X-D*Y*Y=N
plus_n_fundamental_1(XYs0, D, N, XYs).
plus_n_fundamental_1([_|XYs0], D, N, XYs):-
plus_n_fundamental_1(XYs0, D, N, XYs).

% Finds additional solutions of x^2-Dy^2=N when N is not square-free.
% If N is a perfect square, then [sqrt(N),0] is a fundamental solution of
%   x^2-Dy^2=N.
% If R^2 is a factor of N, and [X,Y] is the smallest solution of
%   x^2-Dy^2=N/R^2, then [RX,RY] is a fundamental solution of x^2-Dy^2=N.
plus_n_fundamental_2(R, SqrtN, _, _, XYs, XYs):-
less_than(SqrtN, R), !.            % All values of R have been tested
plus_n_fundamental_2(SqrtN, SqrtN, _, N, XYs, [[SqrtN,[0]]|XYs]):-
multiply(SqrtN, SqrtN, N), !.      % N is a perfect square
plus_n_fundamental_2(R, SqrtN, D, N, XYs0, XYs):-
divide(N, R, NR1, [0]),
divide(NR1, R, NR2, [0]),          % R^2 divides N
plus_n_1(D, NR2, 0, [[X,Y]|_]), !, % [X,Y] is a solution of x^2-Dy^2=N/R^2
multiply(X, R, RX),
multiply(Y, R, RY),                % [RX,RY] is a solution of x^2-Dy^2=N
plus_n_fundamental_2(R1, SqrtN, D, N, [[RX,RY]|XYs0], XYs).
plus_n_fundamental_2(R, SqrtN, D, N, XYs0, XYs):-
plus_n_fundamental_2(R1, SqrtN, D, N, XYs0, XYs).

%%%%%%%%%%%%%%%
% x^2-Dy^2=-N %
%%%%%%%%%%%%%%%

% minus_n(D, N, J, XYs) is true if the elements [X,Y] of XYs are solutions of
%   x^2-Dy^2=-N, where D > 1, N > 0 and N < sqrt(D).  Solutions will be
%   obtained from the first J solutions of x^2-Dy^2=1 and all fundamental
%   solutions of x^2-Dy^2=-N.
% e.g. minus_n(`157`, `12`, 0, XYs) gives
%   XYs=[[`50`,`4`],[`2719`,`217`],[`2271269`,`181267`]]
% e.g. minus_n(`157`, `12`, 1, XYs) gives
%   XYs=[[`50`,`4`],[`2719`,`217`],[`2271269`,`181267`],
%        [`123363799`,`9845503`],
%        [`103049745749`,`8224265053`],
%        [`5597138921710`,`446700316396`],
%        [`4675470012106610`,`373143129538396`],
%        [`253947789893540611`,`20267240045357413`],
%        [`212130749816239256561`,`16929876922062299863`]]
minus_n(D0, `1`, J, XYs):-!,
minus_one(D0, J, XYs).
minus_n(D0, N0, J, ABs):-
J >= 0,
string_bigint(D0, D),
string_bigint(N0, N),
minus_n_1(D, N, J, XYs0),
insertion_sort(XYs0, XYs),
pairs_strings(XYs, ABs).

% This is similar to minus_n/4 but uses bigints, not strings.
minus_n_1(D, N, 0, XYs):-!,       % J=0
minus_n_fundamental(D, N, XYs). % Fundamental solutions of x^2-Dy^2=-N
minus_n_1(D, N, J, XYs):-
plus_one_1(D, J, RSs),          % J solutions of x^2-Dy^2=1
minus_n_fundamental(D, N, PQs), % Fundamental solutions of x^2-Dy^2=-N
findall(XY, one_solution(PQs, RSs, D, XY), XYs0),
append(PQs, XYs0, XYs).

% minus_n_fundamental(D, N, XYs) is true if the elements [X,Y] of XYs are the
%   fundamental solutions of x^2-Dy^2=-N, where N < sqrt(D).
% e.g. minus_n_fundamental([55], [6], XYs) gives
%   XYs=[[[7],[1]], [[37],[5]]]
minus_n_fundamental(D, N, XYs):-
multiply(N, N, N2),
less_than(N2, D),
cont_frac(D, _, As),
convergents(As, XYs0),
minus_n_fundamental_1(XYs0, D, N, XYs1),           % Basic solutions
square_root(N, SqrtN),
minus_n_fundamental_2([2], SqrtN, D, N, [], XYs2), % Additional solutions
append(XYs1, XYs2, XYs).

% Selects those convergents which will give a solution of x^2-Dy^2=-N.
minus_n_fundamental_1([], _, _, []).
minus_n_fundamental_1([[X,Y]|XYs0], D, N, [[X,Y]|XYs]):-
multiply(X, X, X2),
multiply(Y, Y, Y2),
multiply(D, Y2, DY2),
subtract(DY2, X2, N), !, % X*X-D*Y*Y=-N
minus_n_fundamental_1(XYs0, D, N, XYs).
minus_n_fundamental_1([_|XYs0], D, N, XYs):-
minus_n_fundamental_1(XYs0, D, N, XYs).

% Finds additional solutions of x^2-Dy^2=-N when N is not square-free.
% If N is a perfect square and [X,Y] is a solution of x^2-Dy^2=-1, then [SX,SY]
%   is a fundamental solution of x^2-Dy^2=-N, where S=sqrt(N).
% If R^2 is a factor of N, and [X,Y] is the smallest solution of
%   x^2-Dy^2=-N/R^2, then [RX,RY] is a fundamental solution of x^2-Dy^2=-N.
minus_n_fundamental_2(R, SqrtN, _, _, XYs, XYs):-
less_than(SqrtN, R), !.             % All values of R have been tested
minus_n_fundamental_2(SqrtN, SqrtN, D, N, XYs, [[SX,SY]|XYs]):-
multiply(SqrtN, SqrtN, N),          % N is a perfect square
minus_one_1(D, 1, [[X,Y]]), !,      % [X,Y] is a solution of x^2-Dy^2=-1
multiply(X, SqrtN, SX),
multiply(Y, SqrtN, SY).             % [SX,SY] is a solution of x^2-Dy^2=-N
minus_n_fundamental_2(R, SqrtN, D, N, XYs0, XYs):-
divide(N, R, NR1, [0]),
divide(NR1, R, NR2, [0]),           % R^2 divides N
minus_n_1(D, NR2, 0, [[X,Y]|_]), !, % [X,Y] is a solution of x^2-Dy^2=-N/R^2
multiply(X, R, RX),
multiply(Y, R, RY),                 % [RX,RY] is a solution of x^2-Dy^2=-N
minus_n_fundamental_2(R1, SqrtN, D, N, [[RX,RY]|XYs0], XYs).
minus_n_fundamental_2(R, SqrtN, D, N, XYs0, XYs):-
minus_n_fundamental_2(R1, SqrtN, D, N, XYs0, XYs).

% Finds one solution of x^2-Dy^2=N or x^2-Dy^2=-N
member([P,Q], PQs),
multiply(P, R, PR),   % PR=P*R
multiply(D, Q, DQ),
multiply(DQ, S, DQS), % DQS=D*Q*S
multiply(P, S, PS),   % PS=P*S
multiply(Q, R, QR),   % QR=Q*R
one_solution_1(PR, DQS, PS, QR, X, Y),
less_than([0], Y).    % Ignore trivial solutions (Y=0)

one_solution_1(PR, DQS, PS, QR, X, Y):-
absolute_difference(PR, DQS, X),
absolute_difference(PS, QR, Y).
one_solution_1(PR, DQS, PS, QR, X, Y):-

% insertion_sort(Xss, Zss) sorts the list of lists of bigints Xss, using the
%   first bigint of each list to decide the ordering.
insertion_sort(Xss, Yss):-insertion_sort_1(Xss, [], Yss).

insertion_sort_1([], Yss, Yss).
insertion_sort_1([Xs|Xss], Yss0, Yss):-
insert(Yss0, Xs, Yss1),
insertion_sort_1(Xss, Yss1, Yss).

% Duplicate elements are not inserted.
insert([[X|_]|Yss], [X|Xs], Zss):-!,
insert(Yss, [X|Xs], Zss).
insert([[Y|Ys]|Yss], [X|Xs], [[Y|Ys]|Zss]):-
less_than(Y, X), !,
insert(Yss, [X|Xs], Zss).
insert(Yss, Xs, [Xs|Yss]). % Y>X or Ys=[]

%%%%%%%%%%%%%%%
% Cx^2-Dy^2=N %
%%%%%%%%%%%%%%%

% c_plus_n(C, D, N, J, XYs) is true if the elements [X,Y] of XYs are solutions
%   of Cx^2-Dy^2=N, where C > 0, D > 0, N > 0, CN < sqrt(CD), and CN is
%   square-free or a perfect square.
%   Solutions will be obtained from the first J solutions of x^2-CDy^2=1 and
%   all fundamental solutions of x^2-CDy^2=CN.
% e.g. c_plus_n(`3`, `25`, `2`, 3, XYs) gives
%   XYs=[[`3`,`1`],[`153`,`53`],[`7953`,`2755`],[`413403`,`143207`]]
% e.g. c_plus_n(`4`, `73`, `3`, 1, XYs) gives
%   XYs=[[`47`,`11`],[`18203`,`4261`],[`214419203`,`50191739`],
%        [`83051151047`,`19440803989`]]
c_plus_n(C0, D0, N0, J, ABs):-
string_bigint(C0, C),
string_bigint(D0, D),
string_bigint(N0, N),
multiply(C, D, CD),
multiply(C, N, CN),
plus_n_1(CD, CN, J, XYs0),
valid_solutions(XYs0, C, XYs1),
insertion_sort(XYs1, XYs),
pairs_strings(XYs, ABs).

%%%%%%%%%%%%%%%%
% Cx^2-Dy^2=-N %
%%%%%%%%%%%%%%%%

% c_minus_n(C, D, N, J, XYs) is true if the elements [X,Y] of XYs are solutions
%   of Cx^2-Dy^2=-N, where C > 0, D > 0, N > 0, CN < sqrt(CD), and CN is
%   square-free or a perfect square.
%   Solutions will be obtained from the first J solutions of x^2-CDy^2=1 and
%   all fundamental solutions of x^2-CDy^2=-CN.
% e.g. c_minus_n(`3`, `29`, `2`, 3, XYs) gives
%   XYs=[[`3`,`1`],[`171`,`55`],[`9573`,`3079`],[`535917`,`172369`]]
c_minus_n(C0, D0, N0, J, ABs):-
string_bigint(C0, C),
string_bigint(D0, D),
string_bigint(N0, N),
multiply(C, D, CD),
multiply(C, N, CN),
minus_n_1(CD, CN, J, XYs0),
valid_solutions(XYs0, C, XYs1),
insertion_sort(XYs1, XYs),
pairs_strings(XYs, ABs).

valid_solutions([], _, []).
valid_solutions([[CX,Y]|XYs0], C, [[X,Y]|XYs]):-
divide(CX, C, X, [0]), !, % 0=CX mod C and X=CX div C
valid_solutions(XYs0, C, XYs).
valid_solutions([_|XYs0], C, XYs):-
valid_solutions(XYs0, C, XYs).

% last(Xs, X) is true if X is the last element of the list Xs.
last([X], X):-!.
last([_|Xs], X):-last(Xs, X).
```