HP Forums
(38G) x^2+D*y^2=p Diophantine Solution - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (38G) x^2+D*y^2=p Diophantine Solution (/thread-3450.html)



(38G) x^2+D*y^2=p Diophantine Solution - Gerald H - 03-21-2015 07:43 AM

The programme CORNACCHIA finds the unique integer solution { x, y } of

x ^ 2 + D * y ^ 2 = p

given D < p & p prime or returns 0 if there is no solution.

The sub-programme SQRTMODP is here

http://www.hpmuseum.org/forum/thread-3448.html

eg For input

{ 23 694283029607 }

the programme returns

{ 829512 16409 } in Ans

& indeed

829512 ^ 2 + 23 * 16409 ^ 2 = 694283029607

Ans►L1:
Ans(1)►Y:
L1(2)►Z:
{-Y,Z}:
RUN SQRTMODP:
IF Ans
THEN
MAX(Z-Ans,Ans)►B:
Z►A:
INT(√Z)►L:
WHILE B>L
REPEAT
A MOD B►R:
B►A:
R►B
END:
√((Z-B^2)/Y)►C:
IF FRAC(Ans)
THEN
0:
ELSE
{B,C}:
END:
ELSE
0:
END: