(38G) x^2+D*y^2=p Diophantine Solution

[Gerald H]

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: