A simple wp34s mini-challenge
#1

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.

#2

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.

#3

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

#4

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.

#5

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

#6

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

#7

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

#8

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.

#9

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

#10

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.

#11

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

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


- Pauli

#12

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*x2 - y 2 - 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.

#13

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.

#14

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.




#15

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.

#16

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.

#17

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.
#18

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.

#19

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

#20

Quote:
    CLSTK
NOT

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

Gerson.

#21

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

#22

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, m2 + m = 2n2. 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.

#23

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

#24

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.

#25

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

#26

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.

#27

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,...
#28

010 STO

????


- Pauli

#29

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.

#30

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



Possibly Related Threads…
Thread Author Replies Views Last Post
  Simple Tetris. free for you to improve on cyrille de Brébisson 3 1,979 11-20-2013, 05:43 PM
Last Post: Erwin Ried
  [HP-Prime] Simple Game (Bugs) CompSystems 1 1,464 11-01-2013, 10:18 AM
Last Post: Han
  HPCC Mini Conference 2013 hugh steers 6 2,342 09-13-2013, 04:27 PM
Last Post: BruceH
  Simple Math Question Namir 2 1,450 08-09-2013, 06:13 PM
Last Post: Eddie W. Shore
  Picture from the Mini-HHC (LA) Geir Isene 11 3,278 07-07-2013, 01:06 PM
Last Post: Jim Horn
  HP-42S with Electroluminescent screen and simple I/O port Jose Poyan 8 2,859 03-27-2013, 07:11 PM
Last Post: Jose Poyan
  Simple sample programs for the HP-41CX? Tom Lewis 5 2,140 03-25-2013, 07:11 PM
Last Post: Allen
  Simple? programming question [WP34S] Shawn Gibson 3 1,501 03-15-2013, 11:56 AM
Last Post: Didier Lachieze
  My birthday, so a little commemorative mini-challenge ! Valentin Albillo 43 8,744 03-07-2013, 03:44 AM
Last Post: Walter B
  Good (simple?) calculation for benchmarking? Jedidiah Smith 28 7,553 03-01-2013, 05:13 PM
Last Post: Harald

Forum Jump: