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: /
 
 

01  LBL "CNK" 02  ST- Y 03  SIGN 04  STO M         05  CLX

06  LASTX  07  X>Y? 08  X<>Y 09  X=0? 10  GTO 02

11  LBL 01         12  X<>Y 13  ISG X 14  CLX 15  ST* M

16  X<>Y 17  ST/ M  18  DSE X 19  GTO 01         20  LBL 02

21  X<> M         22  X<>Y 23  RDN 24  END

 
   ( 41 bytes / SIZE 000 )
 
 

STACK	INPUTS	OUTPUTS
T	T	n
Z	Z	T
Y	n	Z
X	k	Ckn
L	/	k

 
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 )
 
 

STACK	INPUTS	OUTPUTS
T	T	n
Z	Z	T
Y	n	Z
X	k	Ckn
L	/	k

  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: /
 
 

01  LBL "LCNK"  02  X<>Y  03  STO 00     04  RCL Y  05  -  06  X>Y?  07  X<>Y  08  STO 01         09  X<>Y  10  1  11  +  12  STO 02  13  STO 05  14   E6  15  STO 04  16  5.005  17  STO 03  18  LBL 01  19  RCL 03  20  RCL 00

21  RCL 02  22  ENTER^  23  SIGN  24  +  25  X>Y?  26  GTO 00   27  STO 02         28  LBL 02         29  ST* IND Z  30  ISG Z  31  GTO 02  32  R^  33  LBL 03  34  RCL IND X  35  RCL 04  36  X>Y?  37  GTO 03  38  MOD  39  ST- IND Y  40  X<> IND Y

41  LASTX  42  /  43  ISG Y  44  GTO 04   45  STO IND Y  46  X<>Y  47  INT  48   E3  49  /  50  5  51  +  52  STO 03   53  GTO 01         54  LBL 04  55  ST+ IND Y  56  ENTER^  57  DSE Z  58  LBL 03  59  X<> Z  60  ISG X

61  GTO 03  62  GTO 01  63  LBL 00  64  RCL 03   65  FRC  66   E3  67  *  68  .004  69  +  70  STO 03         71  LBL 05  72  RCL IND X  73  RCL 01  74  MOD  75  ST- IND Y  76  LASTX  77  ST/ IND Z  78  CLX  79  RCL 04  80  *

81  DSE Y  82  X<0?  83  GTO 06   84  ST+ IND Y  85  X<>Y  86  GTO 05  87  LBL 06  88  RCL IND 03  89  X=0?  90  DSE 03  91  RCL 03         92  DSE 01  93  GTO 05  94   E-3  95  +  96  STO 03   97  END

 
     ( 148 bytes / SIZE??? )
 
 

STACK	INPUT	OUTPUT
Y	n	/
X	k	eee.005

  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: /
 
 

01  LBL "MLN" 02  RCL 00         03  1 04  - 05  RCL IND 00

06  LBL 01 07  RCL IND Y  08  X=0? 09  GTO 03         10  LBL 02

11  X<>Y 12  ISG X           13  CLX 14  ST* L 15  X<>Y

16  ST/ L 17  DSE X  18  GTO 02         19  LBL 03 20  RDN

21  DSE Y  22  GTO 01         23  LASTX 24  END

 
   ( 42 bytes / SIZE k+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	n
X	/	P

  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"
 
 

01  LBL "MLN"  02  CLA  03  RCL 00          04  1  05  -  06  SF 25  07  RCL IND 00  08  LBL 01   09  RCL IND Y  10  LBL 02   11  X<>Y  12  ISG X  13  CLX

14  ST* L  15  FS? 25  16  GTO 02          17   E99  18  ST/ L  19  CLX  20  99  21  ST+ M  22  RDN  23  ST* L  24  SF 25  25  LBL 02   26  X<>Y

27  ST/ L  28  DSE X  29  GTO 02          30  RDN  31  DSE Y  32  GTO 01   33  LASTX  34  ENTER^  35  ENTER^  36  LOG  37  INT  38  X<> M  39  STO N

40  SF 24  41  10^X  42  *  43  X<>Y  44  RCL M          45  ST+ N  46  10^X  47  /  48  RCL N   49  X<> Z  50  CF 24  51  CF 25  52  FIX 9

53  "  "  54  ARCL Y   55  >" E"  56  FIX 0  57  CF 29  58  ARCL Z         59  FIX 9  60  SF 29  61  AVIEW  62  END

 
     ( 110 bytes / SIZE kkk+1 )
 
 

STACK	INPUTS	OUTPUTS
T	/	n
Z	/	exponent
Y	/	mantissa
X	/	P

  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 )
 
 

01  LBL "MLN"  02  RCL 00  03  1  04  STO Z  05  -  06  RCL IND 00  07  LBL 01  08  RCL IND Y  09  ST+ Y  10  CNK  11  ST* Z  12  X<> T  13  DSE Y  14  GTO 01  15  RCL Z  16  END

 
( 34 bytes / SIZE kkk+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	n
X	/	P
L	/	n1

  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 )
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	Z	0
Y	Y	n
X	k	P
L	L	k

    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 )
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	Z	Z
Y	/	exp
X	k	mantissa
L	/	n

    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.
 
 

01  LBL "LMLN"   02  " N=?"   03  PROMPT   04  STO 00   05  " N1=?"   06  PROMPT   07  STO 01   08  1   09  +   10  STO 02          11  STO 05    12   E6   13  STO 04   14  5.005   15  STO 03   16  LBL 01   17  RCL 03   18  RCL 00   19  RCL 02   20  ENTER^   21  SIGN   22  +   23  X>Y?   24  GTO 00   25  STO 02   26  LBL 02   27  ST* IND Z   28  ISG Z

29  GTO 02   30  R^   31  LBL 03   32  RCL IND X    33  RCL 04   34  X>Y?   35  GTO 03    36  MOD   37  ST- IND Y   38  X<> IND Y   39  LASTX   40  /   41  ISG Y   42  GTO 04   43  STO IND Y   44  X<>Y   45  INT   46  RCL 04    47  SQRT   48  /   49  5   50  +   51  STO 03          52  GTO 01   53  LBL 04   54  ST+ IND Y   55  SIGN   56  -

57  STO Z   58  LBL 03   59  X<> Z   60  ISG X   61  GTO 03    62  GTO 01    63  LBL 00   64  RCL 03   65  FRC   66  RCL 04   67  SQRT   68  *   69  .004   70  +   71  STO 03    72  SF 10   73  1   74  STO 02          75  RCL 01   76  ST- 00   77  LBL 10   78  ISG 02   79  LBL 12   80  FIX 0   81  CF 29   82  " N"   83  ARCL 02   84  >"=?"

85  FIX 4   86  SF 29   87  PROMPT    88  STO 01   89  RCL 00   90  X<>Y   91  -   92  X<0?   93  GTO 12    94  STO 00          95  1   96  X<>Y   97  X#Y?   98  X=0?   99  CF 10 100  X<>Y 101  RCL 01 102  X#Y? 103  GTO 00 104  FS? 10 105  GTO 10 106  GTO 07 107  LBL 00 108  RCL 03 109  LBL 05 110  RCL IND X 111  RCL 01 112  MOD

113  ST- IND Y  114  LASTX 115  ST/ IND Z 116  CLX 117  RCL 04  118  * 119  DSE Y 120  X<0? 121  GTO 06        122  ST+ IND Y 123  X<>Y 124  GTO 05 125  LBL 06 126  RCL IND 03 127  X=0? 128  DSE 03 129  RCL 03 130  DSE 01 131  GTO 05 132  FS? 10 133  GTO 10 134  LBL 07 135  RCL 03 136   E-3 137  + 138  STO 03 139  END

 
   ( 215 bytes / SIZE ??? )
 
 

STACK	INPUT	OUTPUT
X	/	eee.005

  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 )