Thanks for the link!

The trick is to transform the equation we'd already found

m^{2} + m - 2n^{2} = 0

into an easier to solve Pell's equation:

m^{2} + m = 2n^{2}

m^{2} + m + 3m^{2} + 3m + 1 = 2n^{2} + 3m^{2} + 3m + 1

4m^{2} + 4m + 1 = 2n^{2} + 3(m^{2} + m) + 1

(2m + 1)^{2} = 2n^{2} + 6n^{2} + 1

(2m + 1)^{2} - 8n^{2} = 1

(2m + 1)^{2} - 2(2n)^{2} = 1
Let m + 1 = x and 2n = y:

x^{2} - 2*y^{2} = 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 ( x^{2} - N*y^{2} = 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. *