Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 3,229
Threads: 42
Joined: Jul 2006
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 (x1) / 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.
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
Posts: 2,761
Threads: 100
Joined: Jul 2005
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*(x1)/2
sum of the house numbers to right: n*(n+1)/2  x*(x+1)/2
x*(x1)/2 = n*(n+1)/2  x*(x+1)/2
x*(x1)/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.
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
Posts: 3,229
Threads: 42
Joined: Jul 2006
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 HP16c program in 17 steps that works like the above albeit glacially slowly:
00143,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
01043, 6, 4 F? 4
011 22 1 GTO 1
012 43 21 RTN
01343,22, 1 LBL 1
014 33 Rv
015 1 1
016 40 +
017 22 A GTO A
I thought it could be done.
 Pauli
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 3,229
Threads: 42
Joined: Jul 2006
For the continued fraction solution to this problem have a look here.
I'll leave a program to calculate this as an exercise :)
 Pauli
Posts: 2,761
Threads: 100
Joined: Jul 2005
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^{2}  y ^{2}  y = 0 is
x_{n+1} = 3*x_{n} + 2*y_{n} + 1
y_{n+1} = 4*x_{n} + 4*y_{n} + 1
Here is a 42S program to display the solutions:
00 { 44Byte 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/cgisys/cgiwrap/hpmuseum/archv018.cgi?read=145904
(I knew I'd seen this before...)
Gerson.
Edited: 13 June 2011, 3:14 a.m.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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 HP42S program:
00 { 70Byte 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 { 63Byte 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*nk*d*d)=1
then
WriteLn(n:10:0,' /',d:9:0,' = ',n/d:19:16);
until Abs((n/d)Sqrt(k))<1e16;
Writeln;
Writeln(' ':14,'Sqrt(',k:3:0,' ) = ',r:19:16)
end.
Edited: 13 June 2011, 9:31 p.m.
Posts: 67
Threads: 2
Joined: Jun 2011
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 { 44Byte 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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 1,545
Threads: 168
Joined: Jul 2005
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.
Posts: 67
Threads: 2
Joined: Jun 2011
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 { 39Byte 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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
Posts: 2,761
Threads: 100
Joined: Jul 2005
Quote:
CLSTK
NOT
Missed that one. It has the advantage to save yet one more byte in the HP42S program.
Gerson.
Posts: 67
Threads: 2
Joined: Jun 2011
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 minichallenge. "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 minichallenge 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 { 33Byte 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
Posts: 2,761
Threads: 100
Joined: Jul 2005
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^{2} + m = 2n^{2}. Are you sure the specified range is correct? Here is a bruteforce oneliner:
1 FOR I%=1 TO 1000: R=SQR(I%*(I%+1)/2): IF (RINT(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.
Posts: 67
Threads: 2
Joined: Jun 2011
Hi Gerson
Quote:
1 FOR I%=1 TO 1000: R=SQR(I%*(I%+1)/2): IF (RINT(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
Posts: 2,761
Threads: 100
Joined: Jul 2005
Oops! Forgot to check the results. It should have been
1 FOR I%=1 TO 1700: R=SQR(I%*(I%+1)/2): IF (RINT(R))<>0 THEN NEXT ELSE PRINT R*R: NEXT
run
1
36
1225
41616
1413721
Gerson.
Posts: 67
Threads: 2
Joined: Jun 2011
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 { 30Byte 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
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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,...
Posts: 3,229
Threads: 42
Joined: Jul 2006
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 67
Threads: 2
Joined: Jun 2011
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
