Re: HP-41 Programming Challenge



#22

Just in case this gets buried in the board: I have posted my version of the prime finder program here.


Please post your results or any improvements you might think of.

Thanks,

Axel


#23

Hi, Axel;

I made minor changes to reduce proram size (bytes used) but they are not significant so far. As I am able to show an actual improved "version", I'm posting.

For now, congrats! I like seeing others programs. And yours work pretty fine!

Best regards.

Luiz (Brazil)

#24

here's my entry. it misses out "3".

LBL "MOO"
5
LBL 00
STO 00
9
+
6
/
STO 02
RCL 00
SQRT
3
STO 03
/
INT
1
+
3
*
STO 01
XEQ 04
RCL 00
2
+
STO 00
9
+
6
/
STO 02
1
STO 03
RCL 00
SQRT
ENTER^
INT
STO 01
RDN
XEQ 02
RCL 00
4
+
GTO 00
LBL 01
RCL 01
RCL 01
*
RCL 00
-
SQRT
LBL 02
FRC
X=0?
GTO 03
LBL 04
RCL 03
ST+ 01
RCL 02
RCL 01
X<=Y?
GTO 01
RCL 00
VIEW X
LBL 03
END


#25

Hi Hugh,

you got a bug in your program. It finds 35 as a primnumber. I think it does not really check for 5 beeing a divider of the number thats beeing tested.
I will look deeper into it now to find where it is.

Regards,
Harry


#26

aaawwl. and i thought i had it nailed.... :-)

after more testing, i am not happy with this submission of mine. it appears to start quite fast and after a minute or so, is way ahead of the division programs listed here. however, its long term performance suffers because the algorithm can have worst cases later on.

i submitted another program below, which (although a bit long) is complementary to the division programs. its a lot slower at first, but gets better (against division for bigger numbers). say you start at 100,000. its happy.

the ideal agorithm would be some mix of trial division and the "tortoise".

regards,

#27

You only have forty-nine program lines.

tm


#28

Hi, Trent;

I always like to challenge myself first and then I go for the fight!

I tried once to write a program for the HP25 and I did not succeed. It was supposed to compute nPr and nCr (permutations and combinations given n and r). I succeeded writing it for the HP33, that is essentially an HP25 with subroutine and a few extra functionality. I also wrote it for the HP55 and HP38.

The key its the fact that both HP55 and HP38 have [n!], primarily needed to compute both nPr and nCr. Both HP25 and HP33 do not have [n!] and it is necessary to implement it OR add some math simplification based on numbers entered. I simply coded the original nCr and nPr expresions so they fitted in the HP33 memory. The only restriction is that neither n nor r might be greater than 69. Would it be possible to write a program for the HP25 to compute both nPr and nCr? In my program for the HP33C/E, you simply do this:

key in:  press:   display
[f]CLEAR[PRGM]
n [ENTER] n
r [R/S] nCr
[x<>y] nPr
Pressing [f]CLEAR[PRGM] before [R/S] or before inputting data ensures that the calculator is "pointing" to the first line of the program. I will add the listing below, but should you prefer trying from scratch, you should simply not to look at it |^). I accept challenges if I think I have a chance to succeed, but I do not feel comfortable challenging... I cannot explain why! I'd rather showing what I achieved and ask for new, different, clever solutions.

Thanks.

Luiz (Brazil)

Program Listing (HP33E/C)

01 STO 1
02 Rv
03 GSB 14
04 RCL 0
05 RCL 1
06 -
07 +
08 GSB 10
09 RCL 1
10 GSB 14
11 +
12 ÷
13 g RTN
14 g INT
15 ENTER
16 STO 0
17 1
18 -
19 g x=0
20 g RTN
21 ×
22 f LASTx
23 GTO 17


Edited: 29 July 2003, 1:23 a.m.


#29

Hi Luiz,

Actually, by using the n! function, you restrict yourself to n <= 69 on the older calcs, or n<= 253 on then newer ones. By programming a loop instead that keeps multiplying and dividing as needed, you will need more time (and space, of course), but you will get an extended range of input values that the program can compute for.

I have no access to my calcs at the moment (they're at home, of course, and I'm at work right now), but I'll try to find my old HP41 program I wrote for these functions. Our math teacher got a thrill out of presenting problems that could not be solved with the calculator, and that's why I wrote a program to circumvent these limits.

I could, however, write a HP49 program if I find the time to - I have that one at work with me. But it would be fun to see if it fits into the HP25...

Cheers,
Victor


#30

Hi, Viktor;

thank you for the additional information.

In fact, I was aware of the fact that the [n!] limits for a single factorial are not the same when computing nCr or nPr. I tried to mention this fact when I wrote "add some math simplification based on numbers entered", but I simply missed it. Thank you.

I'm not sure if both HP11C and HP15C were the first HP pocket calculators to offer Py,x and Cy,x as standard resources, but I know the algorithm used in both works things out in a way these limits are surpassed. What I want to know is if it's worth detecting overflow before it happens and mostly if it is possible to do it. I did not take spare time to think of it, but it seems to me it's not possible to detect such overflow occurrence in other way than using calculator's overflow detection after each successive multiplication to accomplish simplified terms for permutation and combination. (did I write a concise, understandable text? I hope so...:)

If I understand it correctly, your program for the HP41 extends [FACT] limits so it computes nCr and nPr with n or c >69 if final result is <= 9.99999999E99. Is it correct?

I'm curious to see it.

About developing a program for the HP49G: I think using RPL would be very interesting, mostly to handle internal factorial limits. What makes me lazy about this is that all RPL models, since the HP28C, already have COMB an PERM as standard resources...

Best regards, Victor.

Luiz (Brazil)


#31

Hi Luiz,

I have not found a way to detect overflow in advance, so the program has to go for it and see if it overflows or not.

My HP41 program circumvents the limit of FACT by explicitly multiplying and dividing things out. The implementation is not difficult. Consider

             n!
nPr = --------
(n-r)!
For n=7 and r=3, this yields
           7x6x5x4x3x2x1
7P3 = ---------------
4x3x2x1
leaving only 7*6*5=210 after cancelling out the corresponding factors in enumerator and denominator. Similarly,
              n!
nCr = ----------
r!(n-r)!
so
          7x6x5x4x3x2x1     7x6x5
7C3 = --------------- = ------- = 35
3x2x1x4x3x2x1 3x2x1
In the second case, it can be debated whether doing the multiplications (more accurate, but higher risk of overflow) or the divisions (less accurate, lower risk of overflow) first is better. I think I have used an intermediate way doing both alternatingly, so 7C3 is evaluated as 7:3x4:2*5 (=35, at least on the HP49G I have available right now).

I will try to find my program, but be warned: Chances are not that high, it's a long time since. If I don't find it, let's simply write new programs. It's fun, after all! Let's try to fit it into the HP25!

And let's not forget to thank Harry for starting this challenge, and Axel for restarting the thread.

Cheers, Victor


#32

Well, here's mine.
Some pitfalls:
. If you calculate COMB(8,3) as 8/3*7/2*6, the answer is not correct (it's 56.00000001)
. If you calculate COMB(335,167) you get an overflow, yet
the answer is < 1e100
. Calculate C(n,p) as C(n,min(p,n-p)) to avoid lengthy calculations (eg C(1000,999) = C(1000,1))
. C(n,0) = 1

This does it all, and the contents of reg Z are preserved
(and end up in Y afterwards). Answers larger than 1e10 may
be slightly inaccurate due to roundoff errors. Not much I can do about that.

01*LBL"COMB"
02 RCL Y
03 X<>Y
04 ST- Y
05 X<Y?
06 X<>Y
07 RDN
08 1 E-3
09 *
10 X<> L
11 STO Z
12 ISG L
13 X=0?
14 GTO 00
15*LBL 02
16 RDN
17 ST* Y
18 LASTX
19 INT
20 ST/ Z
21 DSE Y
22 ISG L
23 GTO 02
24*LBL 00
25 RDN
26 CLX
27 1 E3
28 *
29 END


#33

I meant of course, that my program *avoided* all those pitfalls..
Werner

#34

six bytes shorter, and still keeping the Z-reg
(and even the T-reg in case of C(n,0) or C(n,n)):

01*LBL"COMB"
02 ST- Y
03 X>Y?
04 X<>Y
05 +
06 1 E-3
07 ST* L
08 X<>Y
09 ISG L
10 X=0?
11 GTO 00
12*LBL 02
13 ST* Y
14 LASTX
15 INT
16 ST/ Z
17 RDN
18 DSE X
19 ISG L
20 GTO 02
21*LBL 00
22 RDN
23 1 E3
24 *
25 END

#35

Luiz, Victor,

Here is a program for the HP-25 that calculates either C(n,r) or P(n,r) for all cases where C(n,r) and P(n,r) are less than 1e100. The program fits comfortably in the 49 steps of program memory available on the HP-25. I haven't tested these routines on an actual HP-25, but I have tested them on the simulator found on this museums HP-25 page.

The C(n,r) routine is based on the algorithm described by Victor. It makes use of the identity C(n,r) = C(n,n-r), selecting the minimum value of r and n-r for the calculations. In the loop body, the division is done first, to avoid overflows in cases where the intermediate calculations exceed the dynamic range of the HP-25. Since doing division first can result in non-integer results, the final result is rounded to the nearest integer. This routine does work for the case C(355,167), returning a value of 3.0443588e99. Incidentially, the HP-11C also works for this case, returning a value of 3.044358e99. There is a possibility that the final rounding does not produce the correct result, but all the cases I've tested it with seem to work fine so far.

The P(n,r) routine uses the identity P(n,r) = n x (n-1) x (n-2) x ... x (n-r+1). No divisions required here, so the result is always an integer. There is no limitation on the values of n and r, as long as P(n,r) is less than 1e100. For example, P(5000,27) is calculated to be 6.9446225e99.

To use the C(n,r) routine, position the program counter at step 0 if it's not there already, key in n, hit enter, key in r and hit R/S.

To use the P(n,r) routine, position the program counter at step 30 (GTO 30), key in n, hit enter, key in r and hit R/S.

Enjoy,

Eamonn.

[pre]
01 sto 0 ; X=r, Y=n. Store r in R0
02 x<>y
03 sto 1 ; X=n, Y=r. Store n in R1
04 x<>y ; X=r, Y=n
05 - ; X=n-r
06 rcl 0 ; X=r, Y=n-r
07 x<y ; Want to use the minimum of r and (n-r). If r < n-r ...
08 gto 11 ; then skip ahead 2 steps
09 x<>y ; X=n-r, Y=r
10 sto 0 ; store n-r in R0. This will be used as the loop counter
11 1 ; X=1
12 sto 2 ; R2 = result
13 Rv ; Let's preserve the Z register while we're at it
14 rcl 0 ; X = loop counter
15 x=0 ; have we reached the end of the loop yet?
16 gto 24 ; if so, then exit loop (pc + 8)
17 sto/ 2 ; update the result
18 rcl 1 ; recall next numerator factor
19 sto* 2 ; update the result
20 1
21 sto- 0 ; decrement loop count
22 sto- 1 ; decrement numerator factor
23 gto 14 ; loop back (pc-9)
24 rcl 2 ; place the result in X
25 .
26 5
27 + ; round result to nearest integer
28 int ; by adding 0.5 and taking the integer part
29 gto 00 ; done
30 x<>y ; P(n,r) routine. X=n, Y=r
31 sto 0 ; Store n in R0
32 x<>y ; X=r, Y=n
33 1 ; X=1, Y=r, Z=n
34 - ; X=r-1, Y=n
35 - ; X=n-r+1
36 rcl 0 ; X=n, y=n-r+1
37 1
38 - ; subtract 1 from X
39 x<y ; are we done yet?
40 gto 43 ; if so, then display the result
41 sto* 0 ; update result
42 gto 37 ; loop
43 rcl 0 ; display the result
[\pre]


#36

Sorry about the last post - used the wrong formatting code for the listing. This should be easier to read.
---------------------------------------------------

Luiz, Victor,

Here is a program for the HP-25 that calculates either C(n,r) or P(n,r) for all cases where C(n,r) and P(n,r) are less than 1e100. The program fits comfortably in the 49 steps of program memory available on the HP-25. I haven't tested these routines on an actual HP-25, but I have tested them on the simulator found on this museums HP-25 page.

The C(n,r) routine is based on the algorithm described by Victor. It makes use of the identity C(n,r) = C(n,n-r), selecting the minimum value of r and n-r for the calculations. In the loop body, the division is done first, to avoid overflows in cases where the intermediate calculations exceed the dynamic range of the HP-25. Since doing division first can result in non-integer results, the final result is rounded to the nearest integer. This routine does work for the case C(355,167), returning a value of 3.0443588e99. Incidentially, the HP-11C also works for this case, returning a value of 3.044358e99. There is a possibility that the final rounding does not produce the correct result, but all the cases I've tested it with seem to work fine so far.

The P(n,r) routine uses the identity P(n,r) = n x (n-1) x (n-2) x ... x (n-r+1). No divisions required here, so the result is always an integer. There is no limitation on the values of n and r, as long as P(n,r) is less than 1e100. For example, P(5000,27) is calculated to be 6.9446225e99.

To use the C(n,r) routine, position the program counter at step 0 if it's not there already, key in n, hit enter, key in r and hit R/S.

To use the P(n,r) routine, position the program counter at step 30 (GTO 30), key in n, hit enter, key in r and hit R/S.

Enjoy,

Eamonn.

01 sto 0  ; X=r, Y=n. Store r in R0 
02 x<>y
03 sto 1 ; X=n, Y=r. Store n in R1
04 x<>y ; X=r, Y=n
05 - ; X=n-r
06 rcl 0 ; X=r, Y=n-r
07 x<y ; Want to use the minimum of r and (n-r). If r < n-r ...
08 gto 11 ; then skip ahead 2 steps
09 x<>y ; X=n-r, Y=r
10 sto 0 ; store n-r in R0. This will be used as the loop counter
11 1 ; X=1
12 sto 2 ; R2 = result
13 Rv ; Let's preserve the Z register while we're at it
14 rcl 0 ; X = loop counter
15 x=0 ; have we reached the end of the loop yet?
16 gto 24 ; if so, then exit loop (pc + 8)
17 sto/ 2 ; update the result
18 rcl 1 ; recall next numerator factor
19 sto* 2 ; update the result
20 1
21 sto- 0 ; decrement loop count
22 sto- 1 ; decrement numerator factor
23 gto 14 ; loop back (pc-9)
24 rcl 2 ; place the result in X
25 .
26 5
27 + ; round result to nearest integer
28 int ; by adding 0.5 and taking the integer part
29 gto 00 ; done
30 x<>y ; P(n,r) routine. X=n, Y=r
31 sto 0 ; Store n in R0
32 x<>y ; X=r, Y=n
33 1 ; X=1, Y=r, Z=n
34 - ; X=r-1, Y=n
35 - ; X=n-r+1
36 rcl 0 ; X=n, y=n-r+1
37 1
38 - ; subtract 1 from X
39 x<y ; are we done yet?
40 gto 43 ; if so, then display the result
41 sto* 0 ; update result
42 gto 37 ; loop
43 rcl 0 ; display the result


#37

try C(90,7). It should be 7471375560, but your program returns (probably) 7471375562, if you would have done the
calculations as 90/7*89/6*88/5*87/4*86/3*85/2*84, it would have returned 7471375566. That's why I calculate it as
90/1*89/2*88/3 etc guaranteeing integers all along.
BTW, there's no need to do the divisions first in your case.
Either the result will overflow, or none of the intermediate calculations will.
Werner


#38

Werner,

I don't have a HP-25 at hand, but I calculated C(90,7) on a HP-11C, using the algorithm implemented in my original program. It gives a result of 7471375562, as you expected, which is incorrect.

When I checked my program on the HP-25 simulator I didn't come across any failure cases. This seems to be because of the 17 digits of percision that the simulator uses. On the simulator, the program I wrote does return the correct value of C(90,7) = 7471375560.

I had a look at the programs that you wrote for the HP-41 and now I see how you very elegantly avoided the overflow that can occur with the intermediate calculations.

Here is a version of your HP-41 program for the HP-25. It will not overflow for C(n,r) < 1e100 and should return the correct value for C(n,r) < 1e10. It does not preserve the Z register, as your program does, but then again the HP-25 instruction set isn't nearly as powerful as the HP-41 instruction set. It requires the use of several registers, but uses your method to calculate the minimum of r and n-r to save a couple of steps. There's still room to have both this program and the other program that calculates P(n,r) in the 49 steps of HP-25 program memory.

Regards,

Eamonn.

01  -      ; X=n-r
02 lastx ; X=r, Y=n-r
03 x>=y ;
04 x<>y ; place minimum of r and n-r in X
05 + ; X=n
06 sto 1 ; Store n in R1
07 lastx ; X=min(r,n-r)
08 sto0 ; store in R0. Use as loop counter
09 1 ;
10 sto 2 ; R2 is divisor. Initial value is 1
11 eex
12 3
13 chs
14 sto 3 ; R3 is result so far. Initial value is 1e-3
15 rcl 0 ; X = loop counter
16 x=0 ; have we reached the end of the loop yet?
17 gto 27 ; if so, then exit (pc + 10)
18 rcl 1 ; recall next numerator factor
19 sto* 3 ; update the result
20 rcl 2 ; recall next denominator factor
21 sto/ 3 ; update the result
22 1
23 sto- 0 ; decrement loop count
24 sto- 1 ; decrement numerator factor
25 sto+ 2 ; increment denominator factor
26 gto 15 ; loop (pc-11)
27 rcl 3 ; place the result on the stack
28 eex
29 3
30 * ; multiply result by 1e3
31 gto 00 ; done
#39

I'm not a statistical power user but my 32SII returns 1,6712 e 105 for C(355,167) and my 11C flashes with 9,999999 99 as a result.


#40

Oops I mis-typed in my original message.

That should be C(335, 167) = 3.0443588e99, not C(355, 167). C(355,167) is indeed 1.671163e105.

Rgds.,

Eamonn.

#41

Hi all;

I read the posts in this thread again and I can tell you, if you allow me so: if I am the teacher and these are your answers for the task, I'd never feel more proud for them! I'm impressed.

I think that dealing with "restrictions" is teasing, and dealing with resourceful equipment may be confusing. Using different resources and combining their power so you can achieve results is sometimes a matter of "taste" and specific knowledge.

On the other hand, using all available resources that are not enough to accomplish the task without an specific "treatment" and find a clever "way out" with them is way beyond "taste" and specific knowledge.

I'm amazed.

Best regards.

Luiz (Brazil)

#42

heres another idea. also misses 2 & 3. a bit of a
tortoise, but should have better linearity.

LBL "FOO"
0
STO 04
5
LBL 00
STO 01
0
STO 06
2
STO 08
RCL 01
1
-
STO 07
LBL 05
RCL 07
RCL 07
2
/
INT
STO 09
2
*
x<>y?
GTO 06
1
ST+ 06
RCL 09
STO 07
GTO 05
LBL 06
RCL 06
x=0?
GTO 03
1
ST- 06
LBL 02
RCL 08
STO 02
1
STO 03
RCL 07
STO 09
LBL 08
RCL 09
RCL 09
2
/
INT
STO 09
2
*
x<>y?
XEQ 09
RCL 09
x=0?
GTO 10
RCL 02
RCL 02
*
RCL 01
MOD
STO 02
GTO 08
LBL 10
RCL 03
1
x=y?
GTO 07
CHS
RCL 01
+
x=y?
GTO 07
RCL 06
STO 09
LBL 11
RCL 09
x=0?
GTO 03
RCL 03
STO 02
XEQ 09
RCL 01
1
-
x=y?
GTO 07
1
ST- 09
GTO 11
LBL 07
2
0
4
6
RCL 01
X<Y?
GTO 04
RCL 08
1
+
STO 08
4
x<>y?
GTO 02
LBL 04
RCL 01
VIEW X
LBL 03
RCL 04
4
MOD
2
+
STO 04
RCL 01
+
GTO 00
LBL 09
RCL 02
RCL 03
*
RCL 01
MOD
STO 03
END

Possibly Related Threads...
Thread Author Replies Views Last Post
  HP-41(CL): The easiest way to transfer FOCAL programs from a Linux PC to the HP-41 Geir Isene 13 800 12-05-2013, 02:40 AM
Last Post: Hans Brueggemann
  HP 41 Advanced Programming Tips Michael Fehlhammer 0 154 11-28-2013, 06:11 PM
Last Post: Michael Fehlhammer
  HHC 2013 RPN Programming Challenge - Final Thought Jeff O. 7 381 10-17-2013, 11:02 AM
Last Post: Jeff O.
  Programming CMT EPROM modules for HP-41 Paul Berger (Canada) 4 310 01-22-2013, 06:04 PM
Last Post: Paul Berger (Canada)
  HP-41 CMT EPROM & ZEPROM Programming Dan Grelinger 3 229 11-26-2012, 02:42 AM
Last Post: Diego Diaz
  HHC 2012 RPN Programming Challenge Conundrum Jeff O. 15 543 10-08-2012, 03:34 PM
Last Post: Gerson W. Barbosa
  Mini-challenge: HHC2012 RPL programming contest with larger input David Hayden 14 621 10-05-2012, 10:36 PM
Last Post: David Hayden
  Weekend programming challenge: Euler's Totient function Allen 36 1,152 06-03-2012, 10:39 PM
Last Post: Paul Dale
  41-MCODE: a weekend challenge Ángel Martin 3 195 03-19-2012, 06:49 AM
Last Post: Mike (Stgt)
  A Sunday Programming Challenge Allen 61 1,120 03-17-2012, 06:48 PM
Last Post: Allen

Forum Jump: