Binomial & Multinomial Coefficients for the HP-41

Overview
 

 1°)  Binomial Coefficients

  1-a)  Focal Program
  1-b)  M-Code Routine
  1-c)  Long Integers

 2°)  Multinomial Coefficients

  2-a)  Focal Programs#1&2
  2-b)  Focal Program#3
  2-c)  M-Code Routine#1
  2-d)  M-Code Routine#2
  2-e)  Long Integers
 

-The program given in paragraph 1-a) and the 1st program in paragraph 2-a) are also listed in "Statistics and Probability for the HP-41"
  which also contains an M-Code routine "CNK"
 

1°)  Binomial Coefficients
 

     1-a)  Focal Program
 

-This routine is only a slight modification of a program by Keith Jarret listed in "Synthetic Programming Made Easy"
-I've simply added lines 09-10-20 to avoid a data error if k = 0 or k = n

-"CNK" computes  C_n^k =  n! / [ k! (n-k)! ]
 

Data Registers: /
Flags: /
Subroutines: /
 
 

 
   ( 41 bytes / SIZE 000 )
 
 

 
Example:      n = 100  ,  k = 84

    100  ENTER^
     84   XEQ "CNK"  >>>>   C₈₄¹⁰⁰ = 1.345860628 10¹⁸    ( in 5 seconds )
 

     1-b)  M-Code Routine
 

-This routine uses the 13-digit routines to get a better precision.
 

08B  "K"
00E  "N"
003  "C"
2A0  SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
070   N=C ALL
2BE  C=-C
10E  A=C
0B8  C=Y
01D  C=
060   A+C
0E8   X=C
0B0  C=N ALL
2BE  C=-C
025   C=
061   AB+C
2FE  ?C<0
023   JNC+04
0F8  C=X
0F0  C<>N ALL
0E8  X=C
1A0  A=B=C=0
35C  C=
050   1
0A8  Y=C
0EE   B<>C ALL
089   AB
064   STO Q+
0B0   C=N ALL
3CC  ?KEY                      Loop               Press a key to stop an infinite loop or if you think the routine is too slow
0E7   JC+28d
2EE   ?C#0 ALL
0D3   JNC+26d
0F8   C=X
02E   B=0 ALL
0FA  B<>C M
0AE  A<>C ALL
001   C=
060   AB+1
0E8   X=C
0B0   C=N ALL
269   C=
060   AB/C
0D1  RCL
064   Q+
149   C=
060   AB*CM
0A8  Y=C
089   AB
064   STO Q+
0B0   C=N ALL
02E   B=0 ALL
0FA  B<>C M
0AE  A<>C ALL
009   C=
060   AB-1
070   N=C ALL
323   JNC-28d
3A5  ?NC GO
052   RDN

( 62 words )
 
 

  with  0 <= k <= n

Examples:

    400   ENTER^
    100   XEQ "CNK"  >>>>   C₁₀₀⁴⁰⁰ = 2.241854791 E96                ---Execution time = 13s---

    1234  ENTER^
     567   XEQ "CNK"   >>>>   C₅₆₇¹²³⁴ = 1.166825443 E-32                ---Execution time = 75s---

-This last result seems meaningless !
-In fact it is 1.166825443 E368            ( an HP-48 gives  1.16682544286 E368 )

-Is it possible to get the correct exponent ?
-Simply press  LOG  the integer part is 368. And for the mantissa,  FIX 9   CLA   ARCL X   ANUM   ENTER^   LOG   INT   10^X   /

-This method may be employed provided the exponent does not exceed 499.

-Such numbers are difficult to handle though their behaviour is quite logical:  0  +  will produce  OUT OF RANGE
-But  .2  Y^X  computes correctly the 5th root of  C₅₆₇¹²³⁴  i-e  4.105832027 E73

-Even if the exponent exceeds 499 , the result is still useful:

   2345  ENTER^
    456   XEQ "CNK"   >>>>   C₄₅₆²³⁴⁵ = 1.030717292 E-00                ---Execution time = 61s---

-Press  ENTER^  LOG  ( which gives -499.9868604 )  1000  +  INT  to get the correct exponent i-e  500

  thus,  C₄₅₆²³⁴⁵ = 1.030717292 E500

-This will work until  E999  and even after...

Notes:

-Like the focal program "CNK" , this M-Code routine calculates    C_pⁿ   where  p = inf ( k , n-k )  for faster results since   C_kⁿ =  C_n-kⁿ
-If the inputs are fractional or negative, there will be an infinite loop: press any key to stop the loop.
-Of course, it would be better to test the inputs at the beginning of the routine and return 0 or DATA ERROR in this case.
-I let you perform this improvement...
 

     1-c)  Long Integers
 

-"LCNK" employs several registers to compute and store the binomial oefficients by groups of 6 digits
-It first computes n(n-1) ..... (n-k+1)
-Then, the result is divided by k(k-1).... 1 = k!
 

Data Registers:            R00 = n  ,  R01-R02: temp   ,  R03 = ppp.005 = control number of the result  ,  R04 = 10^6

                                      R05-R06-R07 ............... -Rpp =  the digits of C_n^k  which must be read from  Rpp to R05
Flags: /
Subroutines: /
 
 

 
     ( 148 bytes / SIZE??? )
 
 

  where   0 < k < n < 10000  and   eee.005  is the control number of the result

Example:      n = 100  ,  k = 41

    100  ENTER^
     41   XEQ "CNK"  >>>>   9.005                            ---Execution time = 6m29s---

    R09 = 20116 ,  R08 = 440213 ,  R07 = 369968 ,  R06 = 50635 ,  R05 = 175200

Thus,  C₄₁¹⁰⁰ =   20116,440213,369968,050635,175200

Notes:

-Add zeros on the left if there are less than 6 digits ( like in register R06 in the example above )

-Though the final result is in registers R09 thru R05, more registers were used to calculate n(n-1) ..... (n-k+1)
-In this example, SIZE 018 was necessary.

-Obviously, this program is not fast: even with V41 in turbo mode, about 2 minutes are required to get C₅₀₀¹⁰⁰⁰
-The control number is  54.005
-The first digits are  270288,240945, ......  and the last digits are  .... 799821,216320  ( approximately 2.70288240945 E299 )
-Moreover, SIZE 244 is needed.
 

2°)  Multinomial Coefficients
 

     2-a)  Focal Programs#1&2
 

-This program calculates  P =  ( n₁ , n₂ , ....... n_k ) ! = n ! / ( n₁! n₂! ....... n_k! )   where  n = n₁ + n₂ + ...... + n_k
 

Data Registers:     •  R00 = k                                          ( Registers R00 thru Rkk are to be initialized before executing "MLN" )

                                •  R01 = n₁   •  R02 = n₂   .....................   •  Rkk = n_k
Flags: /
Subroutines: /
 
 

 
   ( 42 bytes / SIZE k+1 )
 
 

  where  P = ( n₁ , n₂ , ....... n_k ) ! = n ! / ( n₁! n₂! ....... n_k! )   with  n = n₁ + n₂ + ...... + n_k

Example:     Calculate  157! / ( 12! 20! 41! 84! )

   ( k = 4 )      n₁ = 12  ,  n₂ = 20  ,  n₃ = 41  ,  n₄ = 84   ( n = 12 + 20 + 41 + 84 = 157 )

   4  STO 00   12  STO 01   20  STO 02   41  STO 03   84  STO 04

   XEQ "MLN"  >>>>  P = 9.078363541 10⁷⁴                                   ---Execution time = 21s---

Notes:

-Store the largest n_i into the last register, this will reduce execution time.
-Lines 08-09-19 are only useful if there is an i for which n_i = 0 , otherwise, they may be deleted.
-"MLN" may also be used to compute the binomial coefficients ( with k = 2 )

-If you get OUT OF RANGE line 14, you can press  E99  ST/ L  RDN  R/S
-This may be done several times - say q - and when the program stops, add  q x 99 to the exponent of the result in X.

-Of course, your HP-41 can do the job for you, like in the variant hereunder:
 

Data Registers:     •  R00 = k > 1                                        ( Registers R00 thru Rkk are to be initialized before executing "MLN" )

                                •  R01 = n₁   •  R02 = n₂   .....................   •  Rkk = n_k

Flags:  F24-F25-F29
Subroutines: /

-Line 55 ,   >" E"  means "append space E"
 
 

 
     ( 110 bytes / SIZE kkk+1 )
 
 

  where  P = ( n₁ , n₂ , ....... n_k ) ! = n ! / ( n₁! n₂! ....... n_k! )   with  n = n₁ + n₂ + ...... + n_k

Example1:     Calculate ( 12 , 20 , 41 , 84 ) !  = 157! / ( 12! 20! 41! 84! )

    n₁ = 12  ,  n₂ = 20  ,  n₃ = 41  ,  n₄ = 84   ( n = 12 + 20 + 41 + 84 = 157 )

   4  STO 00   12  STO 01   20  STO 02   41  STO 03   84  STO 04

   XEQ "MLN"  >>>>  the HP-41 displays:  9.078363541 E74                                   ---Execution time = 29s---

                  and                  X = 9.078363541 E74
                              RDN   Y = 9.078363541 = mantissa
                              RDN   Z =  74  = exponent
                              RDN   T = 157 = n  which is saved in T-register.

Example2:     Calculate  ( 76 , 107 , 112 , 184 ) !  = 479! / ( 76! 107! 112! 184! )

     4   STO 00
    76  STO 01   107  STO 02   112   STO 03   184   STO 04

     XEQ "MLN"   >>>>   the HP-41 displays  3.849801810 E273                                   ---Execution time = 107s---

                  and                  X = 9.999999999 E99           ( out of range: the result can't be displayed by normal means )
                              RDN   Y = 3.849801810 = mantissa
                              RDN   Z =  273  = exponent
                              RDN   T =  479 = n  which is saved in T-register.

Notes:

-If you want to use this program with some n_i = 0 , add  LBL 03 after line 29  and  X=0?  GTO 03  after line 09.
-On the other hand, if you don't want to get the result in the alpha register, delete lines 50 to 61.
 

     2-b)  Focal Program#3
 

-The multinomial coefficient  ( n₁ , n₂ , ....... n_k ) !   may also be computed by a product of binomial coefficients:

   ( n₁ , n₂ , ....... n_k ) !  =  C_n(k-1)^n(k)+n(k-1) + C_n(k-2)^{n(k)+n(k-1)+n(k-2)} + ................... + C_n(1)ⁿ    where  n = n₁ + n₂ + ...... + n_k
 

Data Registers:     •  R00 = k > 1                                        ( Registers R00 thru Rkk are to be initialized before executing "MLN" )

                                •  R01 = n₁   •  R02 = n₂   .....................   •  Rkk = n_k

Flags:  /
Subroutine:   "CNK"  ( M-Code routine or focal program )
 
 

 
( 34 bytes / SIZE kkk+1 )
 
 

  where  P = ( n₁ , n₂ , ....... n_k ) ! = n ! / ( n₁! n₂! ....... n_k! )   and   n = n₁ + n₂ + ...... + n_k

Example:     Calculate  ( 16 , 24 , 41 , 48 ) !

   4  STO 00   16  STO 01   24  STO 02   41  STO 03   48  STO 04

   XEQ "MLN"  >>>>  P = 9.227558922 10⁶⁹                                   ---Execution time = 11s---
 

     2-c)  M-Code Routine#1
 

-Like "CNK",  "MLN"  uses the 13-digit routines to get a better precision.
-It takes k in X-register and gives ( n₁ , n₂ , ....... n_k ) ! in X-register.
 

Data Registers:        R00  unused                                          ( Registers R01 thru Rkk are to be initialized before executing "MLN" )

                                •  R01 = n₁   •  R02 = n₂   .....................   •  Rkk = n_k

Flags:  /
Subroutines: /
 

08E  "N"
00C  "L"
00D  "M"
3C8  CLRKEY
0F8   C=X
128   L=C
38D  ?NCXQ
008   C->S&X                        changes k into a hexadecimal number in S&X field of register C
106   A=C S&X
378   C=c
03C   RCR 3
146   A=A+C S&X
1BC  RCR 11
0A6  A<>C S&X
070   N=C ALL                       CPU register N now contains 2 absolute addresses in hexadecimal:  R00 , Rkk
106   A=C S&X
130   LDI S&X
200   200h
306   ?A<C S&X                     checks that  Rkk  does exist
381   ?NC GO
00A   NONEXISTENT
1A0   A=B=C=0
35C   C=
050    1
0E8   X=C
0EE   B<>C ALL
089   AB
064   STO Q+
0B0  C=N ALL
270   RAMSLCT
038   READATA
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
0AE  A<>C ALL
0A8   Y=C                                                   Y-register now contains n(k)
0B0   C=N ALL                 LOOP1
260   SETHEX
266   C=C-1 S&X
070   N=C ALL
106   A=C S&X
03C  RCR 3
0A6  A<>C S&X
306   ?A<C S&X
3A0  ?NC RTN
270   RAMSLCT
038   READATA
10E  A=C ALL
046   C=0 S&X
270   RAMSLCT
0AE  A<>C ALL
2A0  SETDEC
068   Z=C                         LOOP2
3CC  ?KEY
360   ?C RTN                                              Press a key to stop the routine if you think it's too slow or to stop an infinite loop
2EE   ?C#0 ALL
363    JNC-20d                GOTO LOOP1
0B8   C=Y
02E   B=0 ALL
0FA  B<>C M
0AE  A<>C ALL
001   C=
060   AB+1
0A8   Y=C
078   C=Z
269   C=
060   AB/C
0D1  RCL
064   Q+
149   C=
060   AB*CM
0E8   X=C
089   AB
064   STO Q+
078   C=Z
02E   C
0FA  ->
0AE  AB
009   C=
060   AB-1
323   JNC-28d                       GOTO LOOP2

( 81 words )
 
 

    where  P = ( n₁ , n₂ , ....... n_k ) ! = n ! / ( n₁! n₂! ....... n_k! )   with  n = n₁ + n₂ + ...... + n_k

Examples:

   •   Calculate  ( 16 , 24 , 41 , 48 ) !

  16  STO 01   24  STO 02   41  STO 03   48  STO 04

   4   XEQ "MLN"  >>>>  P = 9.227558919 E69                                   ---Execution time = 10s---
 

   •   Calculate  ( 76 , 107 , 112 , 184 ) !

   76  STO 01   107  STO 02   112   STO 03   184   STO 04

   4    XEQ "MLN"   >>>>   P =  3.849801799 E273                                   ---Execution time = 38s---

-The exponent is actually displayed  E-27 , press  ENTER^  LOG to get the correct exponent
-The mantissa may be viewed by an M-code routine like "VM"  or press  FIX 9   CLA   ARCL X

Notes:

-If k = 0 or 1 , this routine returns 1.
-If k is fractional or negative, 38D  008  changes k in a non-negative hexadecimal integer.
-However, the contents of registers R01 to Rkk are not checked for alpha data and if these registers contain negative or fractional numbers,
 this will produce an infinte loop: simply press any key to stop the routine.

-You can use  LOG  like with "CNK" to use the results if the exponent is > 99 .
 

     2-d)  M-Code Routine#2
 

-This variant also takes k in X-register and returns the mantissa of ( n₁ , n₂ , ....... n_k ) !  in X-register and its exponent inY-register
  which is of course more explicit than the routine above, at least if the exponent is > 99
 
 

Data Registers:        R00  unused                                          ( Registers R01 thru Rkk are to be initialized before executing "MLN" )

                                •  R01 = n₁   •  R02 = n₂   .....................   •  Rkk = n_k

Flags:  /
Subroutines: /
 

08E  "N"
00C  "L"
00D  "M"
3C8  CLRKEY
0F8   C=X
38D  ?NCXQ
008   C->S&X                        changes k into a hexadecimal number in S&X field of register C
106   A=C S&X
378   C=c
03C   RCR 3
146   A=A+C S&X
1BC  RCR 11
0A6  A<>C S&X
070   N=C ALL                       CPU register N now contains 2 absolute addresses in hexadecimal:  R00 , Rkk
106   A=C S&X
130   LDI S&X
200   200h
306   ?A<C S&X                     checks that  Rkk  does exist
381   ?NC GO
00A   NONEXISTENT
04E   C=0 ALL
0A8   Y=C
35C   C=
050    1
0E8   X=C
0B0  C=N ALL
270   RAMSLCT
038   READATA
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
0AE  A<>C ALL
128   L=C                                                   L-register now contains n(k)
0B0   C=N ALL                 LOOP1
260   SETHEX
266   C=C-1 S&X
070   N=C ALL
106   A=C S&X
03C  RCR 3
0A6  A<>C S&X
306   ?A<C S&X
14B  JNC+41
270   RAMSLCT
038   READATA
10E  A=C ALL
046   C=0 S&X
270   RAMSLCT
0AE  A<>C ALL
2A0  SETDEC
268   Q=C                         LOOP2
3CC  ?KEY
360   ?C RTN                                              Press a key to stop the routine if you think it's too slow or to stop an infinite loop
2EE   ?C#0 ALL
363    JNC-20d                GOTO LOOP1
138   C=L
02E   B=0 ALL
0FA  B<>C M
0AE  A<>C ALL
001   C=
060   AB+1
128   L=C
278   C=Q
269   C=
060   AB/C
04E  C=0 ALL
158  M=C ALL
0F8  C=X
149   C=
060   AB*CM
0EE  B<>C ALL
0E8   X=C
01A  A=0 M
0B8  C=Y
20E  C=A+C ALL
0A8  Y=C
278   C=Q
02E   C
0FA  ->
0AE  AB
009   C=
060   AB-1
303   JNC-32d                       GOTO LOOP2
2A0  SETDEC
1A0  A=B=C=0
0F8  C=X
0EE  B<>C
025  C=
060  AB+C
0E8  X=C
0B8  C=Y
2EE  ?C#0 ALL
3A0  ?NC RTN
10E  A=C ALL
130   LDI S&X
013   013
3EE  LSHFA
266  C=C-1 S&X
35E  ?A#0 MS
3EB  JNC-03
38E  RSHFA
0BA  A<>C M
0A8  Y=C
3E0   RTN

( 103 words )
 
 

    with  P =  ( n₁ , n₂ , ....... n_k ) ! = n ! / ( n₁! n₂! ....... n_k! )  = mantissa  x 10^exp    and    n = n₁ + n₂ + ...... + n_k

Examples:

  •   Calculate  ( 76 , 107 , 112 , 184 ) !

    76  STO 01   107  STO 02   112   STO 03   184   STO 04

     4   XEQ "MLN"   >>>>   P =  3.849801799   X<>Y   273   whence  P =  3.849801799  E273                 ---Execution time = 38s---

  •   Calculate  ( 123 , 456 , 789 , 1234 ) !

    123  STO 01   456  STO 02   789   STO 03   1234   STO 04

     4   XEQ "MLN"   >>>>   P =  5.596708693   X<>Y   1311  thus  P = 5.596708693  E1311                    ---Execution time = 182s---

Notes:

-If k = 0 or 1 , this routine returns 1.
-If k is fractional or negative, 38D  008  changes k in a non-negative hexadecimal integer.

-However, the contents of registers R01 to Rkk are not checked for alpha data and if these registers contain negative or fractional numbers,
 this will produce an infinte loop: simply press any key to stop the routine.

-Of course, "MLN" may also be used to compute binomial coefficients with  k = 2

>>>>  Store the largest n_i into the last register to reduce the execution time.
 

     2-e)  Long Integers
 

-If we want the exact result, we must use several registers to store all the digits - here by groups of 6 digits.

-"LMNL" is quite similar to "LCNK" listed in paragraph 1-c)
 

Data Registers:            R00-R01-R02: temp   ,  R03 = ppp.005 = control number of the result  ,  R04 = 10^6

                                      R05-R06-R07 ............... -Rpp =  the digits of  ( n₁ , n₂ , ....... n_k ) !  which must be read from  Rpp to R05
Flags:  F10-F29
Subroutines: /

-Line 84 ,  >  ( before "=?" )  denotes the append character.
 
 

 
   ( 215 bytes / SIZE ??? )
 
 

  with   0 < n_i < n < 10000  and   eee.005  is the control number of the result

Example:    Compute  P = ( 16 , 24 , 41 , 48 ) ! = 129! / ( 16! 24! 41! 48! )

  XEQ "LMLN"  >>>>  "N=?"
          129           R/S     "N1=?"
           48            R/S     "N2=?"                       ( 12m50s )
           41            R/S     "N3=?"                        ( 7m28s )
           24            R/S     "N4=?"                        ( 3m12s )
           16            R/S     16.005                         ( 1m45s )

    and we get in registers R16 thru R05

    P = 9227,558919,428890,741185,817619,801447,846479,197604,492597,531243,374784,200000

-Approximately:    9.227558919 E69

Notes:

-Key in the numbers n_i in decreasing order to save execution time & data registers.
-In this example, SIZE  must be at least 032.
-Add zeros on the left if there less than 6 digits ( for instance read  001234 if a data register contains 1234 )