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

[Gerald H]

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

x ^ 2 + D * y ^ 2 = 4 * p

given 0 < D < 4*p & p prime < 2.5*10^11 or returns 0 if there is no solution..

eg For input

{4563 , 4567}

the programme returns

{4 , 2} in Ans

& indeed

4 ^ 2 + 4563 * 2 ^ 2 = 4 * 4567

The sub-programme SQRTMODP is here

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

Ans►L1:
Ans(1)►Y:
L1(2)►Z:
IF Ans≥2.5E11
THEN
MSGBOX " π ≥ 2.5E11":
ELSE
IF Z==2
THEN
√(8-Y):
IFTE(FRAC(Ans),0,{Ans,1}):
ELSE
{-Y,Z}:
RUN SQRTMODP:
IF Ans
THEN
ABS((Ans-Y) MOD 2*Z-Ans)►B:
2*Z►A:
INT(2*√Z)►L:
WHILE B>L
REPEAT
A MOD B►R:
B►A:
R►B
END:
√((4*Z-B2)/Y)►C:
IFTE(FRAC(Ans),0,{B,C}):
END:
END:
END: