#### Cornacchia's Algorithm

This code requires the Multi-Precision Integer Arithmetic package.

A list of the first few numbers which are the sum of the squares of their halves (see squares/3 below), can be found here.

```%%%%%%%%%%%%%%
% x^2+Dy^2=P %
%%%%%%%%%%%%%%

% cornacchia(D, P, X, Y) is true if (X,Y) is a solution of x^2+Dy^2=P where
%   0 < D < P, and P is prime.
% e.g. cornacchia(`3`, `141592653589`, X, Y) gives X=`295139`, Y=`134766`
%   that is, 295139^2 + 3*134766^2 = 141592653589
% This also succeeds for some composite values of P, for example,
%   cornacchia(`4`, `221`, X, Y) gives X=`5`, Y=`7`
cornacchia(D0, P0, X0, Y0):-
string_bigint(D0, D),
string_bigint(P0, P),
cornacchia_1(D, P, X, Y),
bigint_string(X, X0),
bigint_string(Y, Y0).

cornacchia_1(D, P, X, Y):-
less_than([0], D),
less_than(D, P),
subtract(P, D, P_D),
sqrtmod(P_D, P, R),         % R = sqrt(-D) (modulo P)
cornacchia_2(P, R, P, X),
multiply(X, X, XX),
subtract(P, XX, P_XX),
divide(P_XX, D, PXXD, [0]), % 0 = (P-X*X) mod D, PXXD = (P-X*X) div D
square_root(PXXD, Y),       % Y = trunc(sqrt(PXXD))
multiply(Y, Y, PXXD).       % Y*Y = PXXD

cornacchia_2(_, R, P, R):-
multiply(R, R, RR),
\+ less_than(P, RR), !.     % R*R <= P
cornacchia_2(R0, R1, P, X):-
divide(R0, R1, _, R2),      % R2 = R0 mod R1
cornacchia_2(R1, R2, P, X).

%%%%%%%%%%%%%%%%%%%%%%%%
% x^2+y^2=P1*P2*P3*... %
%%%%%%%%%%%%%%%%%%%%%%%%

% cornacchia(Ps, XYs) is true if the elements [X,Y] of the list XYs are
%   solutions of x^2+y^2=N, where N is the product of the primes in the list
%   Ps, each odd prime is equal to 1 mod 4, and X < Y.
% e.g. cornacchia([`101`,`3541`,`27961`], XYs) gives
%   XYs=[[`1`,`100000`],[`48320`,`87551`],[`19801`,`98020`],[`64700`,`76249`]]
%        so, for example, 48320^2+87551^2 = 101*3541*27961 = 10000000001
cornacchia(Ps, XYs):-
cornacchia_1(Ps, ABs),
findall(XY, cornacchia_2(ABs, XY), XYs).

% e.g. cornacchia_1([`73`,`137`], XYs) gives XYs=[[[8],[3]], [[11],[4]]]
cornacchia_1([], []).
cornacchia_1([P0|Ps], [[X,Y]|XYs]):-
string_bigint(P0, P),
cornacchia_1([1], P, X, Y),
cornacchia_1(Ps, XYs).

% e.g. findall(XY, cornacchia_2([[[8],[3]], [[4],[11]]], XY), XYs) gives
%   XYs=[[`1`,`100`], [`65`,`76`]]
cornacchia_2([XY|XYs], [X0,Y0]):-
cornacchia_3(XYs, XY, [X,Y]),
sort_two(X, Y, X1, Y1),
bigint_string(X1, X0),
bigint_string(Y1, Y0).

% E^2 + F^2 = (A^2 + B^2)(C^2 + D^2) = (AD +/- BC)^2 + (AC -/+ BD)^2
cornacchia_3([], XY, XY).
cornacchia_3([[C,D]|XYs], [A,B], XY):-
multiply(B, C, BC),
multiply(A, C, AC),
multiply(B, D, BD),
absolute_difference(AC, BD, F),
cornacchia_3(XYs, [E,F], XY).
cornacchia_3([[C,D]|XYs], [A,B], XY):-
multiply(B, C, BC),
multiply(A, C, AC),
multiply(B, D, BD),
cornacchia_3(XYs, [E,F], XY).

% "Sorts" the two bigints.
sort_two(X, Y, X, Y):-
\+ less_than(Y, X), !.
sort_two(X, Y, Y, X).

%%%%%%%%%%%%%%%%%%%%
% x^2+y^2=x*10^U+y %
%%%%%%%%%%%%%%%%%%%%

% squares(Ps, U, XYs) is true if the product of the primes in the list Ps is
%   equal to 10^(2U)+1, and the elements [X,Y] of the list XYs are such that
%   X^2+Y^2 = X*10^U+Y.
% e.g. squares([`101`,`3541`,`27961`], `5`, XYs) gives XYs=
%   [[`25840`,`43776`],[`74160`,`43776`],[`17650`,`38125`],[`82350`,`38125`]]
%       so, for example, 74160^2+43776^2=7416043776
%   where 10^10+1 = 10000000001=101*3541*27961
% e.g. squares([`101`,`1369778187490592461`,`722817036322379041`], `19`, XYs)
%   gives XYs=[[`135370711001883190`,`1155587244919948276`],
%              [`9864629288998116810`,`1155587244919948276`]]
%   where 10^38+1 = 101*1369778187490592461*722817036322379041
squares(Ps, U0, XYs):-
cornacchia(Ps, ABs),
string_bigint(U0, U),
power([10], U, U10),         % U10 = 10^U
squares_1(ABs, U10, XYs).

squares_1([], _, []).
squares_1([[A0,B0]|ABs], U10, [[X0,Y0],[X1,Y0]|XYs]):-
string_bigint(A0, A),
divide(A, [2], _, [0]),      % A is even
string_bigint(B0, B),
subtract(U10, A, U10A),
divide(U10A, [2], X, [0]),   % X = (10^U-A) / 2
divide(B1, [2], Y, [0]),     % Y = (B+1) / 2
multiply(X, X, XX),
multiply(Y, Y, YY),
add(XX, YY, N),              % N = X^2 + Y^2
bigint_string(N, N0),
bigint_string(X, X0),
bigint_string(Y, Y0),
concat(X0, Y0, N0), !,       % No leading zeroes are required
subtract(U10, X, U10X),
bigint_string(U10X, X1),     % X1 = 10^U-X
squares_1(ABs, U10, XYs).
squares_1([[A0,B0]|ABs], U10, [[X0,Y0],[X1,Y0]|XYs]):-
string_bigint(A0, A),        % A is odd
string_bigint(B0, B),
subtract(U10, B, U10B),
divide(U10B, [2], X, [0]),   % X = (10^U-B) / 2
divide(A1, [2], Y, [0]),     % Y = (A+1) / 2
multiply(X, X, XX),
multiply(Y, Y, YY),
add(XX, YY, N),              % N = X^2 + Y^2
bigint_string(N, N0),
bigint_string(X, X0),
bigint_string(Y, Y0),
concat(X0, Y0, N0), !,       % No leading zeroes are required
subtract(U10, X, U10X),
bigint_string(U10X, X1),     % X1 = 10^U-X
squares_1(ABs, U10, XYs).
squares_1([_|ABs], U10, XYs):- % No solution
squares_1(ABs, U10, XYs).

% concat(A, B, C) is true if C is the result of appending the string B to the
%   string A.
concat(A, B, C):-
string_chars(A, As),
string_chars(B, Bs),
append(As, Bs, Cs),
string_chars(C, Cs).
```