The Museum of HP Calculators

HP Forum Archive 13

[ Return to Index | Top of Index ]

Re: HP-41 Programming Challenge
Message #1 Posted by Axel Poqué on 28 July 2003, 1:15 p.m.

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

      
Re: HP-41 Programming Challenge
Message #2 Posted by Vieira, Luiz C. (Brazil) on 28 July 2003, 6:47 p.m.,
in response to message #1 by Axel Poqué

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)

      
Re: HP-41 Programming Challenge
Message #3 Posted by hugh steers on 28 July 2003, 7:39 p.m.,
in response to message #1 by Axel Poqué

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

            
Re: HP-41 Programming Challenge
Message #4 Posted by Harry (Germany) on 29 July 2003, 6:42 a.m.,
in response to message #3 by hugh steers

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

                  
Re: HP-41 Programming Challenge
Message #5 Posted by hugh on 31 July 2003, 6:06 a.m.,
in response to message #4 by Harry (Germany)

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,

      
What About an HP-25 Programming Challenge?
Message #6 Posted by Trent Moseley on 28 July 2003, 10:30 p.m.,
in response to message #1 by Axel Poqué

You only have forty-nine program lines.

tm

            
Sometimes I challenge myself... (edited)
Message #7 Posted by Vieira, Luiz C. (Brazil) on 29 July 2003, 1:15 a.m.,
in response to message #6 by Trent Moseley

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.

                  
Re: Sometimes I challenge myself... (edited)
Message #8 Posted by Victor Koechli on 29 July 2003, 2:06 a.m.,
in response to message #7 by Vieira, Luiz C. (Brazil)

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

                        
Allow me some questions?
Message #9 Posted by Vieira, Luiz C. (Brazil) on 29 July 2003, 2:34 a.m.,
in response to message #8 by Victor Koechli

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)

                              
Re: Allow me some questions?
Message #10 Posted by Victor Koechli on 29 July 2003, 5:43 a.m.,
in response to message #9 by Vieira, Luiz C. (Brazil)

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

                                    
Re: Allow me some questions?
Message #11 Posted by Werner Huysegoms on 29 July 2003, 8:00 a.m.,
in response to message #10 by Victor Koechli

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

                                          
Re: Allow me some questions?
Message #12 Posted by Werner Huysegoms on 30 July 2003, 3:04 a.m.,
in response to message #11 by Werner Huysegoms

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

                                          
Re: Allow me some questions?
Message #13 Posted by Werner Huysegoms on 30 July 2003, 7:42 a.m.,
in response to message #11 by Werner Huysegoms

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

                                    
Re: Allow me some questions?
Message #14 Posted by Eamonn on 29 July 2003, 3:05 p.m.,
in response to message #10 by Victor Koechli

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]

                                          
Re: Allow me some questions?
Message #15 Posted by Eamonn on 29 July 2003, 3:12 p.m.,
in response to message #14 by Eamonn

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 

                                                
Re: Allow me some questions?
Message #16 Posted by Werner Huysegoms on 30 July 2003, 3:00 a.m.,
in response to message #15 by Eamonn

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

                                                      
Updated HP-25 C(n,r) program
Message #17 Posted by Eamonn on 30 July 2003, 8:50 p.m.,
in response to message #16 by Werner Huysegoms

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
                                          
Re: Allow me some questions?
Message #18 Posted by Julián Miranda (Spain) on 29 July 2003, 5:50 p.m.,
in response to message #14 by Eamonn

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.

                                                
Re: Allow me some questions?
Message #19 Posted by Eamonn on 29 July 2003, 6:16 p.m.,
in response to message #18 by Julián Miranda (Spain)

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.

                  
Gentlemen, that's amazing!
Message #20 Posted by Vieira, Luiz C. (Brazil) on 31 July 2003, 1:48 a.m.,
in response to message #7 by Vieira, Luiz C. (Brazil)

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)

      
Re: HP-41 Programming Challenge
Message #21 Posted by hugh on 29 July 2003, 8:08 a.m.,
in response to message #1 by Axel Poqué

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


[ Return to Index | Top of Index ]

Go back to the main exhibit hall