Quoting from the Wikipedia article on Ramanujan:

"Imagine that you are on a street with houses marked 1 through n. There is a house in between (x) such that the sum of the house numbers to left of it equals the sum of the house numbers to its right. If n is between 50 and 500, what are n and x?"

Other HP RPN calculators can be used, but only on the wp34s this can be done in 18 steps or less (LBL and RTN included) using the method I have used. It is possible Ramanujan's method based on continued fractions allows for an even smaller program, but I am not aware of it. A solution based on his method is much welcome, no matter the program size, in case anyone knows it.

Interesting challenge. I've no idea about the continued fraction solution.

I suspect the 16c can do this in 18 steps too although I might be wrong here. Anyway the 34s can do a bit better than that:

	001: LBL A

	002: 5

	003: 0

	004: 0

	005: STO I

	006: RCL I

	007: x^2

	008: RCL+ I

	009: SR 01

	010: SQRT

	011: FC? C

	012: RTN

	013: DSZ I

	014: BACK 08

This requires the device to be in integer mode with a sufficient word size (19 or 20 minimum I think -- I used 64) and decimal base. Press A and x is returned in the display. n is in register I.

If you don't allow an integer mode start, change steps 2, 3 & 4 to:

	002: 5

	003: S.L 02

	004: BASE 10

Actually, the base doesn't matter here. We just need to be in integer mode. I love the S.L and S.R commands :-)

How did I arrive at this? A little algebra. The sum of the numbers to the left of x is x (x-1) / 2. The sum of the numbers to the right is (n + x + 1) (n - x) / 2. From this you get that n (n+1) / 2 = x^2 at the solution both sides must be integral. The divide by two is safe since n^2 - n is always even. So we need to check the square root for being integer. Fortunately integer mode square root sets carry if there is a remainder.

- Pauli

Edit: return n in register I instead of n+1.

Edited: 12 June 2011, 9:22 p.m.

And an improvement. Replace lines 2 through 4 with MASKR 09 giving a 12 step solution:

	001: LBL A       entry label

	002: MASKR 09    511

	003: STO I       save in the n register

	004: RCL I       beginning of loop: get n

	005: x^2         n^2

	006: RCL+ I      n^2 + n

	007: SR 01       (n^2 + n) / 2 no chance of carry here

	008: SQRT

	009: FC? C       integer answer for the square root?

	010: RTN         yes then we've finished

	011: DSZ I       decrement the loop counter

	012: BACK 08     back to the top of the loop

- Pauli

Very quick solution! I wonder if Ramanujan was as fast :-)

I have not used a good method, it would have taken too long if done by hand (too many quadratic equations to solve). Anyway, it works:

 sum of the house numbers to left: x*(x-1)/2

sum of the house numbers to right: n*(n+1)/2 - x*(x+1)/2

x*(x-1)/2 = n*(n+1)/2 - x*(x+1)/2

x*(x-1)/2 - (n*(n+1)/2 - x*(x+1)/2) = 0

which simplifies to

2*x^2 - n^2 - n = 0 

Now, all is left to do is solving the quadratic equations for n = 50 and above and checking for positive roots.



001 LBL A

002 FILL

003 2

004 ENTER^

005 0

006 R^

007 x^2

008 RCL+ L

009 +/-

010 SLVQ

011 FP

012 x=0?

013 SKIP 03

014 R^

015 INC X

016 GTO A

017 Rv

018 RTN



50 A   -->   x

x<>y   -->   n

Regards,

Gerson.

If we can input the loop starting value before running the program...

	001: LBL A

	002: STO I

	003: RCL I

	004: x^2

	005: RCL+ I

	006: SR 01

	007: SQRT

	008: FC? C

	009: RTN

	010: INC I

	011: BACK 08

50 A to execute. x is in the display. RCL I for n.

- Pauli

And a version that doesn't use a register:

	001: LBL A

	002: FILL

	003: *

	004: +

	005: SR 01

	006: SQRT

	007: FC? C

	008: RTN

	009: Rv

	010: INC X

	011: GTO A

As above, 50 A to run. x in in the display. X<>Y to get n.

And if we had included the fused multiply add instruction steps 3 and 4 could be merged. Sadly it isn't in the firmware.

Here is a HP-16c program in 17 steps that works like the above albeit glacially slowly:

	001-43,22, A  LBL A

	002-      36  ENTER

	003-      36  ENTER

	004-      36  ENTER

	005-      20  *

	006-      40  +

	007-       2  2

	008-      10  /

	009-  43  25  SQRT

	010-43, 6, 4  F? 4

	011-   22  1  GTO 1

	012-  43  21  RTN

	013-43,22, 1  LBL 1

	014-      33  Rv

	015-       1  1

	016-      40  +

	017-   22  A  GTO A

I thought it could be done.

- Pauli

If you are going to try the 16c version, remember to set the word size large enough. And start with a guess for n around 250 or you'll be waiting and waiting and waiting...

- Pauli

Only 100 seconds to find the answer, starting from 250, WS set to 20. Fast enough! My having resorted to SLVQ looks really silly, that's what happens when we stick to the first idea :-)

Cheers,

Gerson.

Edited: 12 June 2011, 11:06 p.m.

I thought about using SLVQ too but figured the setup would be too expensive in steps. Then I remembered that integer square root did the FP test for me...

- Pauli

Quote:

---

If we can input the loop starting value before running the program...

---

Yes, we CAN! Sorry for not having warned you about that.

Gerson.

For the continued fraction solution to this problem have a look here.

I'll leave a program to calculate this as an exercise :-)

- Pauli

For solving, I'll choose Dario Alpern's solver. A program based on continued fractions would be a bit larger but would run significantly faster. The general solution for 2*x² - y ² - y = 0 is

xn+1 = 3*xn + 2*yn + 1

yn+1 = 4*xn + 4*yn + 1

Here is a 42S program to display the solutions:

00 { 44-Byte Prgm }

01>LBL "R"

02 0

03 STO 00

04 STO 01

05>LBL 00

06 3

07 RCL* 00

08 2

09 RCL* 01

10 +

11 1

12 +

13 4

14 RCL* 00

15 3

16 RCL* 01

17 +

18 1

19 +

20 STO 01

21 X<>Y

22 STO 00

23 STOP

24 GTO 00

25.END.

Notice the solutions are successive approximants to sqrt(2), when divided by each other.

They used to solve this kind of problem centuries ago, only it took a little longer :-)

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv018.cgi?read=145904

(I knew I'd seen this before...)

Gerson.

Edited: 13 June 2011, 3:14 a.m.

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.

Hi Gerson,

if you want it really fast, then try this one on your 42S(or on the WS34). This program will stop and show matching "n" and "x" in Stack y and x. Continue with R/S to see the next pair. I've left out the "50" and "500" condition. You'll see why :-)

00 { 44-Byte Prgm }

01>LBL "S"

02 0

03 1

04>LBL 00

05 RCL ST Y

06 RCL+ ST Y

07 X<>Y

08 RCL+ ST Z

09 RCL+ ST Z

10 RCL ST Y

11 X^2

12 2

13 ×

14 RCL ST Y

15 X^2

16 X<Y?

17 X<>Y

18 Rv

19 RCL ST Y

20 RCL× ST T

21 STOP

22 Rv

23 Rv

24 GTO 00

25 .END. 

Wow! Really fast and optimized!

And you can optimize it even further if you replace lines 12 and 13 with

12 STO+ ST X

This will save one byte and one step.

Thank you for your valuable contribution!

Gerson.

I believe the 34s will require the first few steps to be 0 ENTER 1. The 34s stores numbers one digit per line and so 0 1 is just 1.

That will be something to watch out for if converting programs from the 42s / 41c to the 34s.

FYI.

Hi Gerson and Gene,

your comments made me delete some lines. But now the values 0 and 1 must be on the stack. I have also added some comments, because I think the code is rather kryptisch :-)


A and B are auxiliarey "vars". Instead of using the complicated continuos
fraction algorithm the following rules apply:

A_new = A_old + B_old ; B_new = B_old + 2*A_old

x=A*B n= MIN(2A^2;B^2)

put 0 and 1 on the stack and then XEQ "S"

00 { 39-Byte Prgm }	

01>LBL "S"	

02 RCL ST Y			#RCL A_old

03 RCL+ ST Y		#add B_old giving A_new

04 X<>Y				#get B_old in front

05 RCL+ ST Z	 	#add A_old twice

06 RCL+ ST Z	 	#giving B_new

07 RCL ST Y	 		#get A_new to calculate 2*A_new^2

08 X^2	

09 STO+ ST X	

10 RCL ST Y	 		#get B_new to calculate B_new^2

11 X^2	

12 X<Y?				#we need the smaller one

13 X<>Y				#put the larger one in front if necessary

14 Rv				#and drop it, leaving "n" on stack

15 RCL ST Y			#get B_new

16 RCL× ST T		#multiply with A_new giving "x"	 

17 STOP				#show "n" and "x" on stack Y and X

18 Rv				#do away with "n" and "x"

19 Rv				#leave A_new and B_new an stack Y and X

20 GTO "S"			#calculate the next pair, A_new and B_new now becoming A_old and B_old

21 .END.	

Günter,

I'd rather keep the values 0 and 1 inside the program. It is possible to achieve this on the wp34s using two steps. That's a good reason to define 0^0 = 1 :-)

001 LBL B

002 CLSTK

003 y^x

004 RCL Y

005 RCL+ Y

006 X<>Y

007 RCL+ Z

008 RCL+ Z

009 RCL Y

010 X^2

011 STO+ X

012 RCL Y

013 X^2

014 X<? Y

015 X<>Y

016 Rv

017 RCL Y

018 RCL× T

019 STOP

020 Rv

021 Rv

022 BACK 18

Regards,

Gerson.

There is another way, in fact lots of ways.

    iC 00

    iC 01

or

    CLSTK

    NOT

Might make for an interesting challenge all on its on.

Pity I didn't define complex versions of the internal constants -- then it would be possible in a single step :-)

- Pauli

Quote:

---

    CLSTK

    NOT

---

Missed that one. It has the advantage to save yet one more byte in the HP-42S program.

Gerson.

As I'm aiming more for the 42S and here for the Free42, I'd like to leave the input outside of the program, not only would the inclusion add 2 superficial lines it also required an additional local label.

When playing around to get this piece done I didn't really squeeze the last byte out of it, because my main problem was to recall something I dealt with 25 years ago, when I ran into a similar mini-challenge. "Scrooge" didn't remember how rich he was, but knew he had something between 500.000 and a million Golddollars, and they could be layed down either as a square or a triangle. Like in this mini-challenge there was a lot discussion how to optimize the necessary loops. Then I came with my 1 line BASIC program and the discussion was over.

okay, here is my next optimization, with a slightly different approach, and some lines less :-)

put 0 and 1 on the stack, then XEQ "S"

00 { 33-Byte Prgm }

01>LBL "S"

02 RCL+ ST Y

03 STO+ ST Y

04 X<>Y

05 RCL ST Y

06 X^2

07 STO+ ST X

08 RCL ST Y

09 X^2

10 X<Y?

11 X<>Y

12 Rv

13 RCL ST Y

14 RCL× ST T

15 STOP

16 Rv

17 Rv

18 GTO "S"

19 .END.

BTW, I'm not at all an expert in math. The solution is based on simply analyzing the figures found with brute force that time 25 years ago. The hint to the continoues fraction approach now sheds some light on my findings.

Günter

Quote:

---

"Scrooge" didn't remember how rich he was, but knew he had something between 500.000 and a million Golddollars, and they could be layed down either as a square or a triangle.

---

Yes, that's essentially the same problem as it can be described by the same equation, m² + m = 2n². Are you sure the specified range is correct? Here is a brute-force one-liner:

1 FOR I%=1 TO 1000: R=SQR(I%*(I%+1)/2): IF (R-INT(R))<>0 THEN NEXT ELSE PRINT I%*I%: NEXT

run

 1

 64 

 2401

 82944

Ok

Quote:

---

BTW, I'm not at all an expert in math.

---

Neither am I.

Quote:

---

The solution is based on simply analyzing the figures found with brute force that time 25 years ago.

---

Very insightful analysis, I would say. Congrats!

Gerson.

Hi Gerson

Quote:

---

1 FOR I%=1 TO 1000: R=SQR(I%*(I%+1)/2): IF (R-INT(R))<>0 THEN NEXT ELSE PRINT I%*I%: NEXT
run

1

 64 

 2401

 82944

Ok

---

You would need R to be squared. I% is the base of the triangle, thus the results are wrong.

concerning the boundaries, I now believe the upper limit of "Scrooge's" dollars was two millions.

Günter

Oops! Forgot to check the results. It should have been

1 FOR I%=1 TO 1700: R=SQR(I%*(I%+1)/2): IF (R-INT(R))<>0 THEN NEXT ELSE PRINT R*R: NEXT

run

 1

 36

 1225

 41616

 1413721

Gerson.

Hi Gerson

squeezing is fun :-)

so here is the next condensed Free42 approach, that's really short.

clear the stack and then XEQ "T"

R/S will show the next pair of m and n

00 { 30-Byte Prgm }

01>LBL "T"

02 RCL ST X

03 RCL+ ST Z

04 STO+ ST X

05 RCL+ ST Z

06 1

07 +

08 RCL ST X

09 RCL+ ST T

10 RCL+ ST Z

11 STOP

12 GTO "T"

13 .END.

Günter

Really impressive! On the 34s you would need only 11 steps, or 13 if you wanted to press A and R/S for each pair:

001 LBL A

002 CLSTK

003 RCL X

004 RCL+ Z

005 STO+ X

006 RCL+ Z

007 INC X

008 RCL X

009 PSE 10

010 RCL+ T

011 RCL+ Z

012 STOP

013 BACK 10

Gerson.

001 LBL A

002 CLx

003 2

004 SQRT

005 INC X

006 x^2

007 *

008 INC X

009 IP

010 STO

011 BACK 08

This will generate the sequence

1, 6, 35, 204, 1189, 6930, 40931,...

010 STO

????

- Pauli

STOP

Too cold here right now (6 °C), hence the typo :-)

Gerson.

---------------

001 LBL A

002 # c2

003 2

004 SQRT

005 INC X

006 x^2

007 *

008 CEIL

009 STOP

010 BACK 07 

Edited: 27 June 2011, 8:29 p.m.

Hi,

as the 42S doesn't have the INC nor the BACK commands, I think there's little room, if at all, to make that code shorter.

I think I leave this how it is now.

Have fun, Günter