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(A, D, AD), multiply(B, C, BC), multiply(A, C, AC), multiply(B, D, BD), add(AD, BC, E), absolute_difference(AC, BD, F), cornacchia_3(XYs, [E,F], XY). cornacchia_3([[C,D]|XYs], [A,B], XY):- multiply(A, D, AD), multiply(B, C, BC), multiply(A, C, AC), multiply(B, D, BD), absolute_difference(AD, BC, E), add(AC, BD, F), 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 add(B, [1], B1), 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 add(A, [1], A1), 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).