Statistics

# Statistics & Probability for the HP-41

Overview

1°)  One Variable ( Grouped Data )

1-a)  Moments of a Distribution  ( Program#1 )
1-b)  Moments of a Distribution  ( Program#2 )
1-c)  Arithmetic, Geometric and Harmonic Means

2°)  Two Variables ( Grouped Data )

2-a)  Standard Deviation, Covariance, Coefficient of Correlation
2-b)  Linear Regression

3°)  Probability Functions

3-a)  Chi2-Probability Function
3-b)  Gaussian Probability Function
3-c)  F-Distribution Function

3-c1)  Program#1
3-c2)  Inverse FDF function

3-d)  Student's t-Distribution & UTPT Functions

3-d1)  Program#1
3-d2)  UTPT function
3-d3)  Inverse UTPT function

3-e)  Poisson distribution
3-f)  Binomial Distribution
3-g) Negative Binomial Distribution
3-h) Pascal Distribution
3-i)  Geometric Distribution
3-j)  Hypergeometric Distribution

4°)  Classic Problems

4-a)  Binomial Coefficients ( Focal Program & M-Code Routine )
4-b)  Probability of no repetition in a sample
4-c)  Hypergeometric Distribution ( Urn Problem )
4-d)  Multivariate Hypergeometric Distribution
4-c)  Permutations with repetitions

-See also "Quantiles for the HP-41" & "Least-Squares Approximation for the HP-41"

-The latest addition is paragraph 3-c2)

1°)  One Variable ( Grouped Data )

a) Moments of a Distribution ( Program#1 )

-This program calculates the arithmetic mean, the standard deviation, the skewness and the kurtosis
of a given set of p points  x1 , x2 , ............ ,  xp
with respective frequencies  n1 , n2 , ............ ,  np

Data Registers:    R00 = Sum ni                                         when the routine stops,  R05 = m = arithmetic mean

R01 = Sum ni xi    R03 = Sum ni xi3
R02 = Sum ni xi2   R04 = Sum ni xi4

Flag:   F01                  CF 01  if you want to use the sample standard deviation
Subroutines: /            SF 01  if you want to use the population standard deviation

 01  LBL A 02  ST+ 00 03  RCL Y 04  STO T 05  * 06  ST+ 01 07  * 08  ST+ 02 09  * 10  ST+ 03 11  * 12  ST+ 04 13  RDN 14  RTN 15  LBL "STAT" 16  RCL 01 17  RCL 00 18  / 19  ENTER^ 20  STO 05 21  6 22  * 23  RCL 02  24  * 25  RCL 03 26  4 27  * 28  - 29  * 30  RCL 04 31  + 32  RCL 00 33  / 34  R^ 35  X^2 36  X^2 37  3 38  * 39  - 40  RCL 03  41  R^ 42  RCL 02        43  * 44  3 45  * 46  - 47  RCL 00 48  / 49  RCL 05 50  3 51  Y^X 52  ST+ X 53  + 54  RCL 05 55  X^2 56  RCL 00  57  * 58  RCL 02        59  X<>Y 60  - 61  RCL 00 62  FC? 01 63  DSE X 64  / 65  ST/ Z 66  ST/ Z 67  ST/ Y 68  SQRT 69  ST/ Y 70  RCL 05        71  END

( 95 bytes / SIZE 006 )

 STACK INPUTS#i OUTPUTS#i final OUTPUTS T / / µ4 Z / / µ3 Y xi / s X ni xi m L / / V = s2

m = arithmetic mean                                               V = s2 = variance
s = sample standard deviation if CF 01                        µ3 = skewness
s = population standard deviation if SF 01                   µ4 = kurtosis

Example:    Given the following data

 xi 2 7 9 12 15 19 ni 3 5 8 14 7 2

-Clear registers R00 thru R04  ( for instance execute CLS  if your statistics registers are R00 thru R05 )
-Then,  xi  ENTER^   ni S+  ( [A]  in user mode )  for all the points

2  ENTER^      7  ENTER^      9  ENTER^     12  ENTER^    15  ENTER^     19  ENTER^
3  S+               5  S+                8  S+              14  S+               7   S+               2   S+                in user mode

-Finally, simply press  R/S  ( or XEQ "STAT" )  it yields:

CF 01                                                               SF 01

m = 10.8718                                                       m = 10.8718
RDN    s  =  4.0012                                           RDN    s  =  3.9496
RDN   µ3 = -0.3406 = skewness                        RDN   µ3 = -0.3542 = skewness
RDN   µ4 =  3.0605 =  kurtosis                          RDN   µ4 =  3.2238 =  kurtosis

-To correct wrong inputs, key in  xi  ENTER^   ni  CHS  S+  ( [A] in user mode )  for the wrong ( xi , ni )
or add   LBL a   CHS   just before  LBL A  and key in  xi  ENTER^   niS-  ( [a]  in user mode )

-For a normal distribution,  µ3 = 0  &  µ4 =  3

b) Moments of a Distribution ( Program#2 )

-This routine computes the centered moment  µk =  [ (1/n) Sum ni.( xi - m )k ]/sk   for a given k-value
where  m is the arithmetic mean and s is the standard deviation.
µk is a dimensionless quantity.

Data Registers:       •  R00 = p = number of points ,                  ( Registers R00 thru R2p  are to be initialized before executing "RCMK" )

•  R01 = x1     •  R03 = x2      .......    •  R2p-1 = xp
•  R02 = n1     •  R04 = n2      .......    •  R2p   = np

Flag:   F01                  CF 01  if you want to use the sample standard deviation
Subroutines: /            SF 01  if you want to use the population standard deviation

 01  LBL "RCMK" 02  CLA 03  STO O 04  RCL 00 05  ST+ X 06  STO M 07  ENTER^ 08  CLX 09  LBL 01 10  RCL IND Y 11  ST+ N 12  DSE Z 13  RCL IND Z 14  * 15  + 16  DSE Y 17  GTO 01 18  RCL N 19  / 20  X<>Y 21  LBL 02 22  RCL IND M  23  DSE M 24  RCL IND M 25  R^ 26  ST- Y 27  RDN 28  X^2 29  X<>Y 30  ST* Y 31  X<>Y 32  ST+ P 33  X<> L 34  X<>Y 35  X<> O           36  Y^X 37  LASTX 38  X<> O           39  * 40  + 41  DSE M 42  GTO 02 43  X<>Y 44  RCL P 45  RCL N  46  ST/ T 47  FC? 01 48  DSE X 49  / 50  SQRT 51  ENTER^ 52  X<> O           53  Y^X 54  ST/ Z 55  X<> O  56  X<> Z 57  CLA 58  END

( 97 bytes / SIZE 2p+1 )

 STACK INPUTS OUTPUTS T / mk = Sum ni ( xi - m )k/n Z / s Y / m X k µk = mk/sk L / k

m = arithmetic mean
s = sample standard deviation if CF 01
s = population standard deviation if SF 01

Example:

 xi 2 7 9 12 15 19 ni 3 5 8 14 7 2

-Store these 12 numbers into

R01         R03         ..............................            R11
R02         R04         ..............................            R12         respectively

-There are  p = 6  data points whence   6  STO 00

CF 01    4   XEQ "RCMK"  >>>>  µ4 =  3.0605   RDN   m = 10.8718   RDN   s =  4.0012   RDN   m4 = 784.4261       (  µ4 = kurtosis )
SF 01    4           R/S            >>>>  µ4 =  3.2238   RDN   m = 10.8718   RDN   s =  3.9496   RDN   m4 = 784.4261

CF 01    5           R/S           >>>>  µ5 =  -2.2513   RDN   m = 10.8718   RDN   s =  4.0012   RDN   m5 = -2308.7606    ( in 10 seconds )
SF 01    5           R/S           >>>>  µ5 =  -2.4024   RDN   m = 10.8718   RDN   s =  3.9496   RDN   m5 = -2308.7606

c)  Arithmetic, Geometric and Harmonic Means

-The arithmetic mean m is defined by   m = ( n1 x1 + n2 x2 + ...... + np xp ) / ( n1 + n2 + .....  + np )
-The geometric mean m' is defined by  m' = [ ( x1^n1 ).( x2^n2 ) ...... ( xp^np ) ] 1/n                                        with   n =  n1 + n2 + .....  + np
-The harmonic mean m" is defined by   m" = ( n1 + n2 + .....  + np ) / ( n1/x1 + n2/x2 + ...... + np/xp )

Data Registers:       •  R00 = p = number of points                   ( Registers R00 thru R2p  are to be initialized before executing "MEANS" )

•  R01 = x1     •  R03 = x2      .......    •  R2p-1 = xp
•  R02 = n1     •  R04 = n2      .......    •  R2p   = np
Flags: /
Subroutines: /

 01  LBL "MEANS" 02  RCL 00 03  ST+ X 04  1 05  CLA 06  LBL 01 07  RCL IND Y 08  STO O 09  ST+ P 10  DSE Z 11  RCL IND Z     12  / 13  ST+ N 14  CLX 15  RCL O  16  LASTX 17  * 18  ST+ M  19  X<> L 20  RCL O             21  Y^X 22  * 23  DSE Y 24  GTO 01           25  RCL P  26  ST/ M 27  ST/ N 28  1/X 29  Y^X 30  RCL N             31  1/X 32  X<>Y 33  RCL M  34  CLA 35  END

( 63 bytes / SIZE 2p+1 )

 STACK INPUTS OUTPUTS Z / m" Y / m' X / m

m = arithmetic mean
m' = geometric mean
m" = harmonic mean

Example:

 xi 2 7 9 12 15 19 ni 3 5 8 14 7 2

-Store these 12 numbers into

R01         R03         ..............................            R11
R02         R04         ..............................            R12         respectively

-There are 6 points whence   6  STO 00

XEQ "MEANS"  >>>>   m  = 10.8718        ( execution time = 7 seconds )
RDN    m' =  9.8020
RDN    m" =  8.0549

2°)  Two Variables ( Grouped Data )

a) Standard Deviation, Covariance, Coefficient of Correlation

-This program calculates the arithmetic mean, the standard deviations and the covariance
of a given set of p points  ( x1 , y1 )  , ( x2 , y2 )  ............ ,  ( xp , yp )
with respective frequencies  n1 , n2 , ............ ,  np

Data Registers:    R00 = Sum ni xi     R02 = Sum ni yi      R04 = Sum ni xi yi    and when the routine stops,  R06 = mx = arithmetic mean
R01 = Sum ni xi2    R03 = Sum ni yi2    R05 = n = Sum ni                                                  R07 = my = arithmetic mean

Flag:   F01                  CF 01  if you want to use the sample standard deviation & covariance
Subroutines: /            SF 01  if you want to use the population standard deviation & covariance

-If your statistics registers are R00 thru R05, lines 01 to 22 may be replaced by

LBL A
X<> Z
LBL 01
S+                      ( Sigma+ )
LASTX
DSE Z
GTO 01
RTN

 01  LBL A 02  ST+ 05 03  STO 06 04  X<>Y 05  * 06  STO 07 07  ST+ 02 08  LASTX          09  *  10  ST+ 03  11  X<>Y 12  ENTER^ 13  ST* 07  14  RCL 06  15  *  16  ST+ 00 17  *  18  ST+ 01 19  RCL 07 20  ST+ 04 21  RCL 06 22  RTN 23  LBL "STAT2" 24  RCL 04 25  RCL 02 26  STO 07 27  RCL 00 28  RCL 05 29  ST/ 07 30  / 31  STO 06          32  * 33  - 34  RCL 01 35  RCL 06  36  X^2 37  RCL 05 38  * 39  - 40  RCL 07 41  X^2 42  RCL 05 43  * 44  RCL 03          45  X<>Y 46  - 47  RCL 05  48  FC? 01 49  DSE X 50  ST/ T 51  ST/ Z 52  / 53  SQRT  54  ST/ T 55  X<>Y 56  SQRT             57  ST/ T 58  END

( 83 bytes / SIZE 008 )

 STACK INPUTS#i OUTPUTS#i final OUTPUTS T / / r Z xi / cov(x,y) Y yi / sy X ni ni sx L / / Vx = sx2

s = sample standard deviation if CF 01                             cov(x,y) = sample covariance if CF 01                   r = coefficient of correlation
s = population standard deviation if SF 01                        cov(x,y) = population covariance if SF 01

V = s2 = variance       ,     R06 = mx = arithmetic mean   ,    R07 = my = arithmetic mean

Example:    With the following data

 xi 2 5 10 16 21 yi 4 7 12 19 24 ni 3 4 7 5 2

-Clear registers R00 thru R05  ( for instance execute CLS  if your statistics registers are R00 thru R05 )
-Then,  xi  ENTER^   yi  ENTER^    ni   S+  ( [A]  in user mode )  for all the points

2  ENTER^      5  ENTER^      10  ENTER^     16  ENTER^    21  ENTER^
4  ENTER^      7  ENTER^      12  ENTER^     19  ENTER^    24  ENTER^
3  S+                4  S+                7   S+                5   S+              2   S+            in user mode

-Finally, simply press  R/S  ( or XEQ "STAT2" )  it yields:

if CF 01                         sx  =  5.9622                                       if SF 01                               sx = 5.8185
RDN            sy =  6.3808                                                              RDN            sy  =  6.2270
RDN   cov(x,y) = 38.0143                                                              RDN   cov(x,y) = 36.2041
RDN            r   =  0.9992                                                               RDN            r   =  0.9992

-And we have also

R06 =  mx = 10.3810
R07 =  my = 12.7143

-To correct wrong inputs, key in  xi  ENTER^  yi  ENTER^   ni  CHS  S+  ( [A] in user mode )  for the wrong ( xi , yi , ni )
or add   LBL a   CHS   just before  LBL A  and key in  xi  ENTER^  yi  ENTER^  ni  S-  ( [a]  in user mode )

b) Linear Regression

Data Registers:    R00 = Sum ni xi     R02 = Sum ni yi      R04 = Sum ni xi yi        R06 = mx = arithmetic mean
R01 = Sum ni xi2    R03 = Sum ni yi2    R05 = n = Sum ni         R07 = my = arithmetic mean

Flag:   F01 is used by "STAT2"
Subroutine:  "STAT2"

 01  LBL "LR+"  02  XEQ "STAT2"  03  R^  04  R^  05  LASTX  06  /  07  RCL 06  08  RCL Y  09  *  10  RCL 07  11  X<>Y  12  -  13  X<>Y  14  END

( 29 bytes / SIZE 008 )

 STACK INPUTS OUTPUTS Z / r Y / p X / m

where   y = m x + p  is the least-squares line and  r = linear correlation coefficient

Example:    With the same data as in paragraph 2-a)

XEQ "LR+"  gives   m = 1.0694
RDN   p =  1.6130
RDN   r  =  0.9992

Notes:

-Though "STAT2" produces different outputs if flag F01 is set or clear, "LR+" yields the same results ( with perhaps a small difference in the last decimal )
-See also "LR" & "LSL" listed in "Least-Squares Approximation for the HP-41"

3°)  Probability Functions

a)  Chi2-Probability Function

-This program computes P( X2 | m ) = 2 -m/2 [ 1/Gamma(m/2) ] §0X^2   t -1+m/2 e -t/2 dt

-It uses the series exansion:    P( X2 | m ) = (X2/2)m/2 e -(X^2)/2 [ 1/Gamma(m/2) ]  [ 1 + Sumk=1,2,...  X2k/[ (m+2)(m+4) ..... (m+2k) ]

Data Registers: /
Flags: /
Subroutine:  "GAM+" or "GAM" or "1/G"   ( cf "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )

 01  LBL "PXN" 02  STO N 03  X<>Y 04  STO O 05  2 06  ST+ Y 07  / 08  XEQ "GAM+" 09  RCL N 10  2 11  / 12  E^X 13  ST* Y 14  X<> L 15  RCL O            16  2 17  / 18  Y^X 19  X<>Y 20  / 21  STO M 22  1 23  RCL O            24  1 25  ENTER^  26  LBL 01 27  X<> T 28  RCL N 29  * 30  R^ 31  2 32  + 33  STO T 34  / 35  STO T  36  X<>Y 37  ST+ Y 38  X#Y? 39  GTO 01          40  RCL M 41  * 42  RCL N 43  SIGN 44  X<> O            45  X<>Y 46  CLA 47  END

( 78 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y m m X X2 P( X2 | m ) L / X2

where m is a positive integer ( m > 0 )

Example:

7   ENTER^
12  XEQ "PXN"  >>>>   P (12 | 7 ) = 0.899441131  ( in 14 seconds )

b)  Gaussian Probability Function

"GPF"  calculates  P(x) = (1/sqrt(2.pi)) § -infinityx  exp (-t2/2).dt    via the error function if x >= 0 and via the complementary error function if x < 0

Formulae:

P(x) = [ 1 + erf ( x/sqrt(2) ) ] / 2      if x >= 0
P(x) = (1/2) erfc ( x/sqrt(2) )            if x < 0

Q(x) = 1 - P(x)  verifies  Q(x) = P(-x)

Data Registers: /
Flag:  F10
Subroutine:  "ERF"  ( cf "Error Function for the HP-41" )

 01  LBL "GPF"  02  CF 10  03  X<0?  04  SF 10  05  2  06  SQRT  07  /  08  XEQ "ERF"  09  FS? 10  10  X<>Y  11  .5  12  ST* Y  13  FS?C 10  14  CLX  15  +  16  END

( 34 bytes / SIZE 000 )

 STACK INPUT OUTPUT X x P(x)

Examples:

3.14  XEQ "GPF"  >>>>  P(3.14) = 0.999155261 = Q(-3.14)                  ( in 8 seconds )
-1         R/S          >>>>    P(-1)   = 0.158655254 =   Q(1)                       ( in 7 seconds )

3-c)  F-Distribution Function

3-c1)  Program#1

-This program calculates    Q( x | n1 , n2 ) =  (n1/n2)(n1)/2  Gamma[(n1+n2)/2] / [ Gamma(n1/2) Gamma(n2/2) ] §xinfinity  t -1+(n1)/2 ( 1 + t.n1/n2 ) -(n1+n2)/2 dt

via the incomplete beta function:   Q( x | n1 , n2 ) = Iy( n2/2 , n1/2 )   with  y = n2/(n2+n1.x)

Data Registers:   R00 thru  R05
Flag:   F10
Subroutines:   "BETA"  "IBETA"  ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )

 01  LBL "FDF"      02  RCL Z 03  STO 01 04  * 05  STO Z 06  RCL Y 07  STO 05 08  + 09  ST/ Z 10  / 11  CF 10             12  X>Y? 13  SF 10 14  X>Y? 15  X<>Y 16  STO 00  17  2 18  ST/ 01 19  ST/ 05 20  RCL 01          21  RCL 05 22  XEQ "BETA" 23  X<> 05 24  RCL 01 25  FS? 10 26  X<>Y 27  RCL 00 28  XEQ "IBETA" 29  FS? 10 30  CHS 31  RCL 05 32  FS?C 10          33  ST+ Y 34  / 35  END

( 67 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS Z n1 / Y n2 / X x >= 0 Q( x | n1 , n2 )

Examples:

7    ENTER^
6    ENTER^
4.3   XEQ "FDF"  >>>>  Q ( 4.3 |  7 , 6 ) = 0.04764079963  ( in 24 seconds )

7    ENTER^
9    ENTER^
0.4      R/S       >>>>       Q ( 0.4 | 7 , 9 ) =  0.879873501  ( in 27 seconds )

Notes:

-We have used the relation   Iy( a , b ) = 1 - I1-y( b , a )  to reduce execution time and roundoff errors when y > 0.5
-  Q( 0 | n1 , n2 ) =  1

3-c2)  Inverse FDF function

-"IFDF" employs Halley's iterative method to solve  Q( y | n1 , n2 ) = x   i-e   to find  y = IFDF ( x | n1 , n2 )

tk+1 = tk - 2 f(tk) f '(tk) / [ 2 ( f '(tk) )2 - f(tk) f ''(tk) ]      where   f = FDF

Data Registers:   R00 thru R11: temp
Flags: /
Subroutines:   "BETA"         ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )
"FDF"           ( cf paragraph above )

-Line 21,  1  is the 1st guess whch could probably be replaced by a better estimation
-Lines 59-60 are probably superfluous

 01  LBL "IFDF"  02  STO 06  03  RDN  04  STO 08  05  X<>Y  06  STO 07  07  2  08  ST/ Z  09  /  10  XEQ "BETA"  11  RCL 07  12  RCL 08  13  /  14  RCL 07  15  2  16  /  17  Y^X 18  X<>Y  19  /  20  STO 09  21  1  22  STO 10  23  LBL 01  24  VIEW 10  25  RCL 07  26  RCL 08  27  RCL 10  28  XEQ "FDF"    29  RCL 06  30  -  31  STO 11  32  RCL 10  33  RCL 07  34  2 35  /  36  1  37  -  38  STO 03  39  Y^X  40  RCL 07  41  RCL 10  42  *  43  RCL 08  44  ST+ Y  45  /  46  STO 04          47  RCL 07  48  RCL 08  49  +  50  2  51  / 52  STO 05  53  Y^X  54  /  55  RCL 09  56  *  57  RCL 03  58  RCL 10  59  X=0?  60  SIGN  61  /  62  RCL 07  63  RCL 08          64  /  65  RCL 05  66  *  67  RCL 04  68  / 69  -  70  RCL 11  71  STO T  72  *  73  2  74  /  75  +  76  /  77  ST+ 10  78  ABS  79   E-4  80  X

( 110 bytes / SIZE 012 )

 STACK INPUTS OUTPUTS Z n1 / Y n2 / X x ifdf ( x | n1 , n2 )

Where  y = ifdf ( x | n1 , n2 )  is the solution of   FDF ( y | n1 , n2 ) = x

Examples:

7     ENTER^
6     ENTER^
0.05  XEQ "IFDF"  >>>>   IFDF ( 0.05 | 7 , 6 ) = 4.206658486                     ---Execution time = 72s---

3     ENTER^
5     ENTER^
0.001 XEQ "IFDF"  >>>>  IFDF ( 0.001 | 3 , 5 ) = 33.20246321                    ---Execution time = 93s---

Notes:

-The execution times suppose that "IBETA" calls the M-Code routine HGF+
-Line 79  E-4  is a small number to control the precision of the result.
-Since Halley's method is of order 3, almost all the decimals are already correct !

3-d)  Student's t-Distribution & UTPT Functions

3-d1)  Program#1

- "STD"  computes  A( t | n ) = n -1/2 / Beta(0.5 ; n/2 )  §-t+t  ( 1 + x2/n ) -(n+1)/2  dx

-This routine also computes   n -1/2 / Beta(0.5 ; n/2 )  §t+infinity  ( 1 + x2/n ) -(n+1)/2  dx = [ 1 - A ( t | n ) ] / 2 = UTPT(t,n)  on the HP-48

Formula:   A( t | n ) = 1 - Iy( n/2 , 1/2 )   with   y = n/(n+t2)

Data Registers:   R00 thru  R05
Flag:  F10
Subroutines:   "GAM+" ( or "GAM" )  "IBETA"  ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )

 01  LBL "STD"     02  X^2 03  STO Z 04  RCL Y 05  STO 00 06  + 07  ST/ Z 08  / 09  CF 10 10  X>Y? 11  SF 10 12  X>Y? 13  X<>Y 14  X<> 00 15  .5 16  * 17  STO 05 18  LASTX 19  + 20  XEQ "GAM+" 21  STO 04 22  RCL 05 23  XEQ "GAM+" 24  PI 25  SQRT 26  * 27  RCL 04 28  / 29  X<> 05 30  .5 31  FS? 10 32  X<>Y 33  RCL 00 34  XEQ "IBETA" 35  ENTER^ 36  CHS 37  RCL 05 38  + 39  FC?C 10 40  X<>Y 41  LASTX           42  ST/ Z 43  ST+ X 44  / 45  X<>Y 46  END

( 83 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS Y n UTPT(t,n) X t >= 0 A( t | n )

Examples:

4   ENTER^
1   XEQ "STD"   >>>>            A( 1 | 4 ) = 0.626099033        ( in 21 seconds )
RDN     UTPT( 1 , 4 ) = 0.186950483

4   ENTER^
10    R/S         >>>>            A( 10 | 4 ) = 0.999437997          ( in 16 seconds )
RDN     UTPT( 10 , 4 ) = 0.0002810018113

-We have used the symmetry relation   Iy( a , b ) = 1 - I1-y( b , a )  to reduce execution time and roundoff errors when y > 0.5
- A( 0 | n ) = 0

3-d2)  UTPT function

-There is however a possible loss of significant digits with "STD" if y > 0.5
-To overcome that difficulty, the following program employs a quadratic transformation to calculate the hypergeometric function F. It yields:

UTPT (t | n)  = { 1/[ n.B(n/2,1/2) ] }  yn/2  F ( 1 , n ; 1+n/2 ; y / [ 2 + 2 (1-y)1/2 ] )     where  B = beta function  and  y = n/(n+t2)

Data Registers:   R00 thru  R06
Flags:  /
Subroutines:   "1/G+" ( or "GAM+" )    ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )
"HGF"                                           ( cf  "Hypergeometric functions for the HP-41" )

 01  LBL "UTPT"  02  STO 04  03  X<>Y  04  STO 02  05  2  06  /  07  STO 03  08  XEQ "1/G+"  09  STO 06  10  RCL 03 11  .5  12  +  13  XEQ "1/G+"  14  PI  15  SQRT  16  *  17  RCL 06  18  /  19  STO 06  20  1 21  STO 01  22  ST+ 03  23  RCL 02   24  RCL 02  25  RCL 04         26  X^2  27  STO T  28  +  29  ST/ Z  30  / 31  STO 05  32  X<>Y  33  SQRT  34  1  35  +  36  ST+ X  37  /  38  XEQ "HGF"  39  RCL 05  40  RCL 02 41  2  42  /  43  Y^X  44  *  45  RCL 02   46  RCL 06         47  *  48  /  49  END

( 77 bytes / SIZE 007 )

 STACK INPUTS OUTPUTS Y n / X t >= 0 UTPT(t,n)

Examples:

•  4    ENTER^
2.8  XEQ "UTPT"  >>>>   UTPT ( 2.8 | 4 ) = 0.02440577530           ---Execution time = 10s---

•  7    ENTER^
10     R/S        >>>>   UTPT ( 10 | 4 ) = 1.069710145 E-5            ---Execution time = 8s---

Note:

-If you to compute  A(t | n) , simply add   ENTER^   ST+ X  1  -  CHS  X<>Y   after line 48.

3-d3)  Inverse UTPT function

-IUTPT uses the same expression to calculte UTPT(t,n)  and find a solution of the equation  UTPT(t,n) - A = 0  where  A and n are given
-It employs Halley's method which is even better than Newton's method:

tk+1 = tk - 2 f(tk) f '(tk) / [ 2 ( f '(tk) )2 - f(tk) f ''(tk) ]     with   f =  UTPT

Data Registers:   R00 thru  R09
Flags:  /
Subroutines:   "1/G+" ( or "GAM+" )    ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )
"HGF"                                           ( cf  "Hypergeometric functions for the HP-41" )

 01  LBL "IUTPT"  02  STO 06  03  RDN  04  STO 02  05  1  06  STO 01  07  STO 03  08  +  09  2  10  /  11  STO 08  12  X<>Y  13  STO 04  14  RCL 02  15  2  16  /  17  ST+ 03 18  STO 05  19  XEQ "1/G+"  20  STO 07  21  RCL 08  22  XEQ "1/G+"   23  PI  24  SQRT  25  *  26   RCL 07  27  /  28  ST* 04  29  LBL 01  30  VIEW 06  31  RCL 02  32  RCL 02  33  RCL 06  34  X^2 35  STO T  36  +  37  STO 09  38  ST/ Z  39  /  40  STO 07  41  X<>Y  42  SQRT  43  1  44  +  45  ST+ X  46  /  47  XEQ "HGF"   48  RCL 07  49  RCL 05  50  Y^X  51  * 52  RCL 02  53  /  54  RCL 04  55  -  56  RCL 07  57  RCL 08  58  Y^X  59  RCL 02  60  SQRT  61  /  62  RCL 02   63  RCL 06  64  *  65  RCL 08          66  *  67  RCL 09  68  X^2 69  /  70  R^  71  *  72  RCL 07   73  /  74  -  75  /  76  ST+ 06  77  ABS  78   E-8  79  X

( 117 bytes / SIZE 010 )

 STACK INPUTS OUTPUTS Z A Y n / X t0 IUTPT(A,n)

where  t0  is an initial guess  and  t = IUTPT(A,n) is  the solution of the equation  UTPT(t,n) = A

Example1:          A = 0.0001 ,  n = 4  ,  with the initial approximation    t0 = 10

0.0001   ENTER^
4        ENTER^
10       XEQ "IUTPT"  >>>>   t = 13.03367172                   ---Execution time = 22s---

Example2:          A = 0.025 ,  n = 21  ,  with the initial approximation    t0 = 1

0.025  ENTER^
21     ENTER^
1      XEQ "IUTPT"   >>>>   t = 2.079613847                   ---Execution time = 63s---

Example3:          A = 0.001 ,  n = 21  ,  with the initial approximation    t0 = 1

0.001  ENTER^
21     ENTER^
1      XEQ "IUTPT"   >>>>   t = 3.527153671                   ---Execution time = 66s---

3-e)  Poisson Distribution

-The function is defined by   P(X=k) = exp(-L) Lk / k!        where  L > 0  and  k is a non-negative integer

Data Registers: /
Flags: /
Subroutines: /

 01  LBL "POI"  02  RCL Y  03  X<>Y  04  Y^X  05  LASTX  06  FACT  07  /  08  X<>Y  09  E^X  10  /  11  END

( 20 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y L / X k P( X=k )

Example:           L = 5  ,  k = 6

5   ENTER^
6   XEQ "POI"   >>>>   P(X=6) = 0.146222808

-The cumulative distribution is given by   P( X <= k ) = exp(-L) SUMj=0,1,...,k  Lj / j!

Data Register:  R00 = k , k-1 , ..... , 0
Flags: /
Subroutines: /

 01  LBL "POI+"  02  X=0?  03  GTO 02  04  STO 00  05  CLX  06  LBL 01  07  RCL Y  08  RCL 00  09  Y^X  10  LASTX  11  FACT  12  /  13  +  14  DSE 00  15  GTO 01  16  LBL 02  17  1  18  +  19  RCL Y  20  E^X  21  /  22  END

( 36 bytes / SIZE 001 )

 STACK INPUTS OUTPUTS Y L / X k P( X<=k )

Example:           L = 5  ,  k = 6

5   ENTER^
6   XEQ "POI+"   >>>>   P(X<=6) = 0.762183463

3-f)  Binomial Distribution

-This routine computes  P(X=k) = Cnk pk (1-p)n-k    where   0 <= k <= n  are integers and  0 < p < 1

Data Registers:  R00 thru R02: temp
Flags: /
Subroutine:  "CNK"  or an M-Code CNK  cf §4-a) below

 01  LBL "BNP"  02  STO 00  03  X<>Y  04  STO 01  05  X<>Y  06  Y^X  07  X<>Y  08  STO 02  09  RCL 00  10  XEQ "CNK"  11  *  12  1  13  RCL 01  14  -  15  RCL 02  16  RCL 00  17  -  18  Y^X  19  *  20  END

( 32 bytes / SIZE 003 )

 STACK INPUTS OUTPUTS Z n / Y p / X k P( X=k )

Example:        n = 100 ,  p = 1/6

100   ENTER^
6     1/X
20    XEQ "BNP"   >>>>   P(X=20) = 0.067861995

-The cumulative distribution could be calculated by calling "BNP"  (k+1) times.
-The following routine is faster, at least for large k-values:

Data Registers:  R00 thru R03: temp
Flags: /
Subroutines: /

 01  LBL "BNP+" 02  STO 00 03  SIGN 04  X<>Y 05  STO 01 06  - 07  ST/ 01 08  X<>Y 09  STO 02 10  Y^X 11  SIGN 12  RCL 00         13  X=0? 14  GTO 02  15  CLX 16  STO 03 17  LASTX 18  LBL 01 19  LASTX  20  RCL 01         21  * 22  RCL 02 23  RCL 03  24  - 25  * 26  RCL 03  27  1 28  + 29  STO 03         30  / 31  + 32  DSE 00 33  GTO 01         34  SIGN 35  LBL 02  36  LASTX 37  END

( 50 bytes / SIZE 004 )

 STACK INPUTS OUTPUTS Z n / Y p / X k P( X<=k )

Example:        n = 100 ,  p = 1/6

100   ENTER^
6     1/X
15    XEQ "BNP+"   >>>>   P(X<=15) = 0.387657550

3-g)  Negative Binomial Distribution

-Let  r a positive integer, k a non-negative integer and a real p  0 < p < 1

P(X=k) = Ck+r-1r-1 pr (1-p)k

Data Registers:  R00-R01: temp
Flags: /
Subroutine:  "CNK"  or an M-Code CNK  cf §4-a) below

 01  LBL "JRP"  02  STO 00  03  SIGN  04  X<>Y  05  STO 01  06  -  07  RCL 00  08  Y^X  09  RCL 01  10  RCL Z  11  Y^X  12  *  13  RCL 00  14  R^  15  1  16  -  17  ST+ Y  18  XEQ "CNK"  19  *  20  END

( 34 bytes / SIZE 002 )

 STACK INPUTS OUTPUTS Z r / Y p / X k P( X=k )

Example:       r = 5  ,  p = 1/7

5   ENTER^
7   1/X
10  XEQ "JRP"  >>>>   P(X=10) = 0.012748996

"JRP+"  computes the cumulative negative binomial distribution

Data Registers:  R00 thru R03 temp
Flags: /
Subroutines: /

 01  LBL "JRP+" 02  STO 00       03  RDN 04  STO 01 05  X<>Y 06  STO 02 07  Y^X 08  ABS 09  RCL 00       10  X=0? 11  GTO 02 12  CLX 13  STO 03 14  X<>Y 15  LBL 01 16  LASTX 17  RCL 02  18  RCL 03       19  + 20  * 21  1 22  ST+ 03 23  RCL 01  24  - 25  * 26  RCL 03       27  / 28  + 29  DSE 00 30  GTO 01       31  SIGN 32  LBL 02  33  LASTX 34  END

( 47 bytes / SIZE 004 )

 STACK INPUTS OUTPUTS Z r / Y p / X k P( X<=k )

Example:       r = 5  ,  p = 1/7

5   ENTER^
7   1/X
10  XEQ "JRP+"  >>>>   P(X<=10) = 0.051581690

3-h)  Pascal Distribution

-The function is defined as   P(X=k) = Ck-1r-1 pr (1-p)k-1    where   r and k are positive integers ( r <= k )  and  0 < p < 1

Data Registers:  R00-R01-R02: temp
Flags: /
Subroutine:  "CNK"  or an M-Code CNK  cf §4-a) below

 01  LBL "PRP"   02  STO 00 03  1 04  - 05  RDN 06  STO 01 07  RDN 08  STO 02       09  SIGN 10  - 11  XEQ "CNK" 12  RCL 01 13  RCL 02 14  Y^X 15  * 16  1 17  RCL 01        18  - 19  RCL 00 20  RCL 02 21  - 22  Y^X              23  * 24  END

( 36 bytes / SIZE 003 )

 STACK INPUTS OUTPUTS Z r / Y p / X k P( X=k )

Example:       r = 10  ,  p = 0.6

10   ENTER^
0.6  ENTER^
15   XEQ "PRP"  >>>>  P(X=15) = 0.123958563

-Now,  "PRP+"  calculates  P(X<=k)

Data Registers:  R00 thru R03: temp
Flags: /
Subroutines: /

 01  LBL "PRP+" 02  RDN 03  STO 01       04  X<>Y 05  STO 02 06  STO 03 07  R^ 08  XY 11  - 12  STO 00       13  RCL 01 14  RCL 02 15  Y^X 16  ABS 17  ISG 00 18  LBL 01 19  DSE 00  20  X<0? 21  RTN 22  LASTX 23  1 24  RCL 01 25  - 26  * 27  RCL 03  28  ST* Y 29  1 30  + 31  STO 03  32  RCL 02 33  - 34  / 35  + 36  GTO 01       37  LBL 02 38  CLX 39  END

( 53 bytes / SIZE 004 )

 STACK INPUTS OUTPUTS Z r / Y p / X k P( X<=k )

Example:       r = 10  ,  p = 0.6

10   ENTER^
0.6  ENTER^
15   XEQ "PRP"  >>>>  P(X<=15) = 0.403215551

Note:      P(X=k)  is in L-register

3-i)  Geometric Distribution

-This probability function is defined by                                    P(X=k)  =  p (1-p)k-1                  where   0 < p < 1  and  k is a positive integer
-The cumulative distribution may easily be computed since     P(X<=k) = 1 - (1-p)k

-So a simple program can return both results.

Data Registers: /
Flags: /
Subroutines: /

 01  LBL "GP"  02  1  03  RCL Z  04  -  05  ENTER^  06  X<> Z  07  DSE X  08  INT  09  Y^X  10  ST* Z  11  *  12  1  13  X<>Y  14  -  15  X<>Y  16  END

( 27 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y p P( X<=k ) X k P( X=k )

Example:    p = 0.6

0.6    ENTER^
4      XEQ "GP"   >>>>    P(X=4)  = 0.0384
RDN    P(X<=4) = 0.9744

3-j)  Hypergeometric Distribution

-Here,   P(X=k)  =  Cmk Cn-ma-k / Cna    where   n , m , a , k  are integers with  m <= n ,   k <= a <= m  ,  a - k <= n - m    ( cf § 4-c) for a concrete example )

Data Register:  R00: temp
Flags: /
Subroutine:    "CNK"  ( cf 4-a) )

 01  LBL "HGD"  02  RDN  03  STO 00  04  X<>Y  05  R^  06  ST- 00  07  XEQ "CNK"  08  RDN  09  XEQ "CNK"  10  ST/ Z  11  X<> T  12  X<>Y  13  -  14  RCL 00  15  XEQ "CNK"  16  *  17  END

( 40 bytes / SIZE 001 )

 STACK INPUTS OUTPUTS T n / Z m / Y a / X k P(X=k)

Example:     n = 30  ,  m = 12  ,  a = 10

30  ENTER^
12  ENTER^
10  ENTER^
5   XEQ "HGD"  >>>>   P(X=5) = 0.225856303

- "HGD+" below computes  P(X<=k)

Data Registers:  R00 thru R05: temp
Flags: /
Subroutine:  "CNK"  or an M-Code CNK  cf §4-a)

 01  LBL 'HGD+" 02  STO 00 03  RDN 04  STO 01 05  RDN 06  STO 02 07  X<>Y 08  STO 03 09  - 10  RCL 01 11  + 12  STO 05 13  X<0? 14  CLX 15  STO 04 16  RCL 00 17  - 18  X>0? 19  GTO 02 20  CHS 21  STO 00 22  RCL 02 23  RCL 04 24  XEQ "CNK" 25  RCL 03 26  RCL 02 27  - 28  RCL 01 29  RCL 04 30  - 31  XEQ "CNK" 32  * 33  RCL 03 34  RCL 01 35  XEQ "CNK" 36  / 37  ABS 38  ISG 00 39  LBL 01 40  DSE 00 41  X<0? 42  RTN 43  LASTX 44  RCL 01         45  RCL 04  46  - 47  * 48  RCL 02 49  RCL 04 50  - 51  * 52  RCL 04 53  1 54  + 55  STO 04 56  ST/ Y 57  RCL 05         58  - 59  / 60  + 61  GTO 01 62  LBL 02 63  CLX 64  END

( 90 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS T n / Z m / Y a / X k P(X<=k)

Example1:       n = 30  ,  m = 12  ,  a = 10

30  ENTER^
12  ENTER^
10  ENTER^
5   XEQ "HGD+"  >>>>   P(X<=5) = 0.881728367

Example2:       n = 100  ,  m = 84  ,  a = 41

100   ENTER^
84    ENTER^
41    ENTER^
37    XEQ "HGD+"   >>>>  P(X<=37) =  0.958543068

4°)  Classic Problems

4-a)  Binomial Coefficients

-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  Cnk =  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 Cnk L / k

Example:      n = 100  ,  k = 84

100  ENTER^
84   XEQ "CNK"  >>>>   C10084 = 1.345860628 1018    ( in 5 seconds )

-One could of course use the buil-in FACT function, but there will be an OUT OF RANGE if n > 69

-The following M-Code routine performs the same calculations.
-Moreover, 0 is returned if n or k or n-k is negative or fractional.
-However, it does not check for alpha data or overflow.
-Synthetic register Q is used.

08B   "K"
00E   "N"
003   "C"
2A0   SETDEC
3C8   CLRKEY
128   WRIT 4(L)
0E8   WRIT 3(X)
04E   C=0 ALL
0A8  WRIT 2(Y)
2FE  ?C#0 MS
0BF  JC+23d=+17h
084   CLRF 5                 C
0ED  ?NCXQ                 =
064   193B                   frc(C)
2EE   ?C#0 ALL
097   JC+18d=+12h
2FE  ?C#0 MS
07F   JC+15d=+0Fh
084   CLRF 5                 C
0ED  ?NCXQ                 =
064   193B                   frc(C)
2EE   ?C#0 ALL
057   JC+10d=+0Ah
2BE  C=-C-1 MS         C=-C
070   N=C ALL
10E   A=C ALL
01D  ?NCXQ               C=
060   1807                   A+C
2FE  ?C#0 MS
01B   JNC+03
3A5   ?NCGO          ?NCgoto
052    14E9                  RDN
10E    A=C ALL
0F0    C<>N ALL
01D  ?NCXQ               C=
060   1807                   A+C
2FE  ?C#0 MS
01F   JC+03
070   N=C ALL
04E   C=0 ALL
268   WRIT 9(Q)
35C  PT=12
050   LD@PT- 1
0A8  WRIT 2(Y)
10E   A=C ALL
0B0   C=N ALL
36E   ?A#C ALL
36B   JNC-19d=-13h
1BE   A=A-1 MS
01D  ?NCXQ               C=
060   1807                   A+C
10E   A=C ALL
135   ?NCXQ               C=
060   184D                  A*C
0A8  WRIT 2(Y)
00E   A=0 ALL             A
35C  PT=12                  =
162   A=A+1@PT        1
01D  ?NCXQ               C=
060   1807                   A+C
268   WRIT 9(Q)
10E   A=C ALL
0AE  A<>C ALL
261   ?NCXQ               C=
060   1898                   A/C
3CC  ?KEY
360    ?C RTN                                        the routine stops here if you press a key ( useful if, for example,  n = 10000000 and k = 1000000 )
31B    JNC-29d=-1Dh

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

Example:      n = 400  ,  k = 100

400  ENTER^
100   XEQ "CNK"  >>>>   C400100 = 2.241854802 1096    ( in 12 seconds )

-The last decimals should be 792 instead of 802

-With, for instance, n = 4 and k = 0.5  "CNK" returns 0.
-This may be useful when "CNK" is used in the cumulative distributions above.

4-b)  Probability of no repetition in a sample

-The typical exercise is:  "In a room containing m persons, what is the probability P that nobody has the same birthday?"  ( neglecting the leap years )

-The answer is  P = ( 1 - 1/n ).( 1 - 2/n )........( 1 - (m-1)/n )           with  n = 365

Data Registers: /
Flags: /
Subroutines: /

 01  LBL "NRS"  02  X<>Y  03  ENTER^  04  SIGN  05  LBL 01  06  LASTX  07  *  08  X<>Y  09  ST/ Y  10  X<>Y  11  DSE L  12  ""   13  DSE Z  14  GTO 01  15  END

( 27 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y n n X m probability

Example:    n = 365 ,  m = 24

365   ENTER^
24    XEQ "NRS"  >>>>  probability = 0.4617   ( in 8 seconds )

4-c)  Hypergeometric Distribution ( Urn Problem )

"An urn contains m black marbles and (n-m) white marbles. You take simultaneously ( or successively without replacement )  "a"  marbles.
What is the probability P that you have drawn exactly k black marbles ( and (a-k) white ones ) ?"

-The answer is:    P = Cmk Cn-ma-k / Cna

-So we can use the program "HGD" listed in § 3-j above.

Example:    The urn contains 60 marbles ( 15 black ones and 45 white ones )  You take 10 marbles. Probability that you have exactly 3 black marbles?

n = 60 , m = 15 , a = 10 , k = 3

60  ENTER^
15  ENTER^
10  ENTER^
3   XEQ "HGD"  >>>>   P = 0.273864227

4-d)  Multivariate Hypergeometric Distribution

-The previous exercise may be generalized if the urn contains    n1 marbles of color 1  ,   n2 marbles of color 2  ,   .............................  ,  nk marbles of color k
-You take simultaneously - or successively without replacement -  m marbles.
-What is the probability P that you have drawn exactly    p1 marbles of color 1  ,   p2 marbles of color 2  ,   .............................  ,  pk marbles of color k
( of course  m =  p1 +  p2 + ............  +  pk )

-The program below computes the answer   P = Cn1p1 Cn2p2 ......  Cnkpk   / Cn1+n2+....+nkp1+p2+....+pk

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

•  R01 = n1     •  R03 = n2      .......    •  R2k-1 = nk
•  R02 = p1     •  R04 = p2      .......    •  R2k   = pk
Flags: /
Subroutine:   "CNK"  ( cf 4-a) )

 01  LBL "MHGD" 02  RCL 00 03  ST+ X 04  STO N 05  1 06  LBL 01 07  RCL IND Y 08  DSE Z 09  RCL IND Z 10  X<>Y 11  XEQ "CNK"  12  * 13  DSE Y 14  GTO 01 15  0 16  ENTER^ 17  LBL 02 18  RCL IND N 19  + 20  DSE N 21  X<>Y 22  RCL IND N  23  + 24  X<>Y 25  DSE N 26  GTO 02 27  XEQ "CNK"  28  / 29  END

( 58 bytes / SIZE 2k+1 )

 STACK INPUT OUTPUT X / probability

Example:        An urn contains 12 marbles ( 3 blue ones , 5 red ones , 4 green ones )  You take simultaneously 7 marbles
Calculate the probability P that you have exactly  2 blue marbles , 3 red marbles and 2 green marbles ?

k = 3     n1 = 3     n2 =  5     n3 = 4     are to be stored into      R00     R01     R03     R05
p1 = 2     p2 = 3      p3 = 2     are to be stored into                  R02     R04     R06      respectively

XEQ "MHGD"  >>>>  probability = 0.227272727 = 5/22    ( in 6 seconds )

4-c)  Permutations with repetitions

-This program calculates  P = n! / ( n1! n2! ....... nk! )   where  n = n1 + n2 + ...... + nk

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

•  R01 = n1   •  R02 = n2   .....................   •  Rkk = nk
Flags: /
Subroutines: /

 01  LBL "PERM" 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

( 43 bytes / SIZE k+1 )

 STACK INPUTS OUTPUTS Y / n X / P

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

( k = 4 )      n1 = 12  ,  n2 = 20  ,  n3 = 41  ,  n4 = 84   ( n = 12 + 20 + 41 + 84 = 157 )

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

XEQ "PERM"  >>>>  P = 9.078363541 1074  ( in 21 seconds )

Example2:     How many anagrams can you write with  the 16-letter "word"   AAABBCCCCEEEEELL

( k = 5 )     n1 = 3 A's  ,  n2 = 2 B's  ,  n3 = 4 C's  ,  n4 = 5 E's  ,  n5 = 2 L's    (  n = 16 )

5  STO 00      3  STO 01     2  STO 02      4  STO 03     5  STO 04     2  STO 05

XEQ "PERM"   yields   302,702,400  words ( including AAABBCCCCEEEEELL itself )  So, this "word" has 302,702,399  ( meaningless ) anagrams

Hint:   Store the largest ni into the last register, this will reduce execution time.

References:

[1]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4