Thanks for the link!
The trick is to transform the equation we'd already found
m2 + m - 2n2 = 0
into an easier to solve Pell's equation:
m2 + m = 2n2
m2 + m + 3m2 + 3m + 1 = 2n2 + 3m2 + 3m + 1
4m2 + 4m + 1 = 2n2 + 3(m2 + m) + 1
(2m + 1)2 = 2n2 + 6n2 + 1
(2m + 1)2 - 8n2 = 1
(2m + 1)2 - 2(2n)2 = 1
Let m + 1 = x and 2n = y:
x2 - 2*y2 = 1
Now the problem is ready to be quickly solved by the following HP-42S program:
00 { 70-Byte Prgm } 27 1/X
01>LBL "R" 28 FP
02 2 29 STO 03
03 SQRT 30 1/X
04 STO 01 31 IP
05 IP 32 STO 04
06 STO 02 33 RCL 07
07 LAST X 34 STO 05
08 FP 35 RCL 02
09 STO 03 36 X^2
10 1/X 37 RCL 06
11 IP 38 X^2
12 STO 04 39 STO+ X
13 1 40 -
14 STO 05 41 1
15 STO 06 42 X=Y?
16>LBL 00 43 XEQ 01
17 RCL 02 44 GTO 00
18 STO 07 45>LBL 01
19 RCL* 04 46 RCL 02
20 RCL+ 05 47 2
21 STO 02 48 BASE/
22 RCL 01 49 RCL 06
23 / 50 2
24 IP 51 /
25 STO 06 52 STOP
26 RCL 03 53.END.
keystrokes display
XEQ "R" 1
1
R/S 8
6
R/S 49
35
R/S 288
204
R/S 1681
1189
R/S 9800
6930
R/S 57121
40391
This will solve Pell's equation for any positve integer N, such that N is not a perfect square, as long as machine accuracy allows it:
Generic Pell's Equation Solver ( x2 - N*y2 = 1 )
00 { 63-Byte Prgm } 27 1/X
01>LBL "P" 28 FP
02>STO 00 29 STO 03
03 SQRT 30 1/X
04 STO 01 31 IP
05 IP 32 STO 04
06 STO 02 33 RCL 07
07 LAST X 34 STO 05
08 FP 35 RCL 02
09 STO 03 36 X^2
10 1/X 37 RCL 06
11 IP 38 X^2
12 STO 04 39 RCL* 00
13 1 40 -
14 STO 05 41 1
15 STO 06 42 X=Y?
16>LBL 00 43 XEQ 01
17 RCL 02 44 GTO 00
18 STO 07 45>LBL 01
19 RCL* 04 46 RCL 02
20 RCL+ 05 47 RCL 06
21 STO 02 48 STOP
22 RCL 01 49 .END.
23 /
24 IP
25 STO 06
26 RCL 03
keystrokes display
2 XEQ "P" 3
2
R/S 17
12
R/S 99
70
R/S 577
408
R/S 3363
2378
R/S 19601
13860
R/S 114243
80782
I'll leave optimization and porting to wp34s as an exercise :-)
Gerson.
P.S.: Be warned this is my own algorithm for generating the approximants to sqrt(N), which means there is quite room for improvement.
TurboBDC code
Program Pell;
var b, c, d, k, n, n0, r, s, t, x: Real;
function Sqrt(x: Real): real;
var s, t: Real;
begin
if x<>0 then
begin
s:=x/2;
repeat
t:=s;
s:=(s+x/s)/2
until s=t;
Sqrt:=s
end
else
Sqrt:=0
end;
begin
ClrScr;
ReadLn(k);
ClrScr;
r:=Sqrt(k);
n0:=1;
d:=1;
n:=Int(r);
b:=Frac(r);
c:=Int(1/b);
repeat
t:=n;
n:=c*n+n0;
d:=Int(n/r+0.5);
b:=Frac(1/b);
c:=Int(1/b);
n0:=t;
if (n*n-k*d*d)=1
then
WriteLn(n:10:0,' /',d:9:0,' = ',n/d:19:16);
until Abs((n/d)-Sqrt(k))<1e-16;
Writeln;
Writeln(' ':14,'Sqrt(',k:3:0,' ) = ',r:19:16)
end.
Edited: 13 June 2011, 9:31 p.m.