Factorials for the HP-41

Overview
 

 1°)  Integer Arguments

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

 2°)  Real Arguments

  a)  Focal Program
  b)  M-Code Routine
  c)  Log ( n! )

 3°)  Multifactorials

  a)  Focal Program
  b)  M-Code Routine
  c)  Log ( M(n) )
  d)  Long Integers

 4°)  Superfactorials

  a)  Focal Program
  b)  M-Code Routine
  c)  Log ( Sf(n) )
  d)  Long Integers

 5°)  Hyperfactorials

  a)  Focal Program
  b)  3 M-Code Routines
  c)  Log ( H(n) )
  d)  Long Integers
 
 

1°)  Integer Arguments
 

     a)  Focal Program
 

-The HP-41 has of course a built-in function  FACT which works perfectly to get  n!  provided  n < 70  ( otherwise, OUT OF RANGE )
-The internal form

  3B1   ?NCXQ    C=
  060    18EC      fact(C)         may be used if  n < 100

-So, we must write programs if we want to overcome these limitations.

-"FCT" returns n! in register X ( 9.999999999 E99 if n > 69 ), the mantissa in Y-register and the exponent in Z-register
-It also displays the result which is stored in the alpha register.
 

Data Registers: /
Flags:  F24-F25-F29
Subroutines: /

-Line 48 is:   append space E
 
 

01  LBL "FCT"  02  SIGN  03  0  04  X<>Y  05  X<> L  06  X=0?  07  GTO 02         08  SF 25  09  LBL 01  10  ST* L  11  FS? 25

12  GTO 01   13   E99  14  ST/ L  15  CLX  16  99  17  ST+ Z  18  RDN  19  ST* L  20  SF 25  21  LBL 01         22  DSE X

23  GTO 01   24  LBL 02   25  X<> L  26  ENTER^  27  ENTER^  28  LOG  29  INT  30  X<> T  31  STO M         32  SF 24  33  10^X

34  *  35  X<>Y  36  R^  37  ST+ M  38  10^X  39  /  40  0  41  X<> M         42  X<> Z  43  CF 24  44  CF 25

45  FIX 9  46  "  "  47  ARCL Y  48  >" E"  49  FIX 0  50  CF 29  51  ARCL Z        52  FIX 9  53  SF 29  54  AVIEW   55  END

 
     ( 99 bytes / SIZE 000 )
 
 

STACK	INPUTS	OUTPUTS
Z	/	exponent
Y	/	mantissa
X	n	n !

  and  n !  is stored in alpha register and displayed at the end

Example:

  357   XEQ "FCT"   >>>>   8.608917246 E757  is displayed                  ---Execution time = 78s---

     Press the backarrow key   9.999999999 E99
                     RDN                 8.608917246         the mantissa
                     RDN                 757                        the exponent

-The last 3 digits of the mantissa should be 198 instead of 246
 

     b)  M-Code Routine#1
 

-"FCT"  calculates n!  by computing  n(n-1)(n-2).... 1  with a 13-digit routine to get a better precision.

-Synthetic register Q is used.
-There is no check for alpha data.
 

094   "T"
003   "C"
006   "F"
2A0   SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
1A0  A=B=C=0
35C  C=
050   1
0E8   X=C
0EE   B<>C ALL
089   AB
064   STO Q+
138   C=L
070   N=C ALL                            LOOP
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
3A0  ?NC RTN
02E   C
0FA  ->
0AE  AB
0D1  RCL
064  Q+
149  C=
060  AB*CM
0E8  X=C
089  AB
064  STO Q+
0B0  C=N ALL
02E   C
0FA  ->
0AE  AB
009   C=
060   AB-1
35B  JNC-21d

( 37 words )
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	Z	Z
Y	Y	Y
X	n	n !
L	/	n

  where  n is a non-negative integer.

Examples:

   10    XEQ "FCT"  >>>>   362880

  100   XEQ "FCT"  >>>>   9.332621544 E-43                   ---Execution time = 7s---

-The exponent seems wrong but press  ENTER^  LOG  and you get  157.9700037  so the exponent is 157   and  100 ! =  9.332621544 E157

-This method may be used provided the exponent remains smaller than  500, for instance:

  253   XEQ "FCT"  >>>>   5.173460992 E499                     ---Execution time = 17s---

-But let's try   449   XEQ "FCT"  >>>>  3.851930518 E-3                      ---Execution time = 31s---

 Here, press ENTER^  LOG  1000  +  and it yields  997.5856784  thus,   449 ! = 3.851930518 E997

-Likewise, we can find   450 ! = 1.733368733 E1000  and so on ...

Notes:

-Though the HP-41 cannot display correctly the numbers > 9.999999999 E99, we can get the correct exponent by such tricks.

-If the argument is not a non-negative integer, there will be an infinite loop: press any key to stop the loop.
-Of course, it would be preferable to test the input at the beginning of the routine...

Variant:

-The routine will run about 33% faster if we do not use   089   064  ( AB  STO Q+ )  and  0D1  064  ( RCL Q+ )
 

094   "T"
003   "C"
006   "F"
2A0   SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
04E  C=0 ALL
268  Q=C
35C  C=
050   1
0E8   X=C
138   C=L
070   N=C ALL                     LOOP
3CC  ?KEY
360   ?C RTN                                              Press a key to stop the routine if you think it's too slow
2EE   ?C#0 ALL
09B   JNC+19d
278   C=Q
10E  A=C ALL
0F8  C=X
0EE  B<>C ALL
0B0  C=N ALL
13D  C=
060  AB*C
0EE  B<>C ALL
0E8  X=C
0AE  A<>C ALL
268  Q=C
0B0  C=N ALL
02E   C
0FA  ->
0AE  AB
009   C=
060   AB-1
353  JNC-22d                       GOTO LOOP
278  C=Q
10E  A=C ALL
0F8  C=X
0EE  B<>C
04E  C=0 ALL
025  C=
060  AB+C
0E8  X=C
3E0  RTN

( 45 words )
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	Z	Z
Y	Y	Y
X	n	n !
L	/	n

  where  n is a non-negative integer.

Example:

   449   XEQ "FCT"  >>>>  3.851930518 E-3                      ---Execution time = 20s---   ( instead of 31s )

   ENTER^  LOG  1000  +  and it yields  997.5856784  thus,   449 ! = 3.851930518 E997
 

     c)  M-Code Routine#2
 

-It is of course better to get the exponent more explicitly.

-This program takes a non-negative integer n in register X and returns the mantissa of n! in X-register and the exponent in Y-register.
-If n is negative or fractional, you'll get a DATA ERROR message.
-However, there is no check for alpha data.

-Press any key to stop the routine if you wish.
 

094   "T"
003   "C"
006   "F"
2A0   SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
2FE  ?C#0 MS
0B5   ?CGO
0A3   DATAERROR
084   C
0ED  =
064   frc(C)
2EE  ?C#0 ALL
0B5   ?CGO
0A3   DATAERROR
078   C=Z
028   T=C
0B8   C=Y
068   Z=C
04E  C=0 ALL
0A8  Y=C
35C   C=
050    1
0E8   X=C
138   C=L
2EE  ?C#0 ALL
3A0  ?NC RTN
070   N=C ALL                     LOOP
3CC  ?KEY
360   ?C RTN                                              Press a key to stop the routine if you think it's too slow
2EE   ?C#0 ALL
0A3   JNC+20d
1A0  A=B=C=0
0F8  C=X
0EE  B<>C ALL
0B0  C=N ALL
13D  C=
060  AB*C
0EE  B<>C ALL
0E8  X=C
01A  A=0 M
0B8  C=Y
20E  C=A+C ALL
0A8  Y=C
0B0  C=N ALL
02E   C
0FA  ->
0AE  AB
009   C=
060   AB-1
34B  JNC-23d                       GOTO LOOP
1A0  A=B=C=0
0F8  C=X
0EE  B<>C
025  C=                                 Theoretically, the mantissa may be rounded to 10 if the 13-digit mantissa is > 9.999999999499
060  AB+C                            but that seems highly improbable in practice.
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

( 72 words )
 
 

STACK	INPUTS	OUTPUTS
Y	/	exponent
X	n	mantissa
L	/	n

  where n is a non-negative integer

Examples:

    0     XEQ "FCT"   >>>>   1  X<>Y  0  so  0! = 1 E0 = 1

    7     XEQ "FCT"   >>>>   5.040   X<>Y   3          so   7! =  5.040 E3 = 5040

  357   XEQ "FCT"   >>>>   8.608917198  X<>Y   757    whence   357 ! =  8.608917198  E757                ---Execution time = 16s---

  449   XEQ "FCT"  >>>>    3.851930518  X<>Y   997    whence   449 ! =  3.851930518  E997                 ---Execution time = 20s---

 1234  XEQ "FCT"   >>>>   5.108498145  X<>Y  3280   whence  1234 ! = 5.108498145 E3280               ---Execution time = 58s---

 2345  XEQ "FCT"   >>>>   4.452388605  X<>Y  6886   whence  2345 ! = 4.452388605 E6886               ---Execution time = 112s---

Notes:

-Register Z & Y are saved in T & Z respectively.
-This routine is actually much faster than the 1st one !

-Even with 13-digit routines, thousands of multiplications may produce small roundoff-errors:
-In the last example, the exact value is   2345 ! = 4.45238860677855988 ......  E6886   ( I will not write the 6887 digits... )
 

     d)  Long Integers
 

-Except for small arguments, FACT & "FCT" only return approximate results.
-If we want all the digits, several data registers must be used to store them.

-"LFCT" calculates n!  by groups of 6 digits.  The REGMOVE function from the X-functions module is used.
 

Data Registers:            R00 = number of zeros to be added to the result in R01-R02 ................

                                      Rmm , Rm-1 , Rm-2 , ............... , R01 =  the digits of n !
Flags: /
Subroutines: /
 
 

01  LBL "LFCT"  02  CLRG   03   E3  04  /  05  10  06  +  07  STO M         08  9  09  FACT  10  STO 01   11  SIGN

12  STO N  13   E6  14  STO O  15  LBL 01   16  SIGN  17  RCL N   18  RCL O         19  SQRT  20  /  21  +  22  ENTER^

23  CLX  24  LBL 02   25  RCL IND Y  26  RCL M   27  INT  28  *  29  +  30  ENTER^  31  ENTER^  32  RCL O         33  ST/ Z

34  MOD  35  STO IND Z  36  RDN  37  INT  38  ISG Y  39  GTO 02   40  ST+ IND Y  41  X#0?  42  ISG N   43  CLX  44  RCL 01

45  X#0?  46  GTO 03  47  6  48  ST+ 00  49  RCL N   50  RCL O         51  /  52  2.001  53  +  54  REGMOVE  55  DSE N

56  LBL 03   57  ISG M  58  GTO 01   59  RCL M         60  INT  61  DSE X  62  RCL N   63   E-3  64  +  65  CLA  66  END

 
    ( 109 bytes / SIZE ??? )
 
 

STACK	INPUTS	OUTPUTS
Y	/	n
X	n	mmm.001

  where  n is an integer ( 9 < n < 1000 ) and  mmm.001  is the control number of the result

Example:

   100   XEQ "LFCT"   >>>>   23.001                               ---Execution time = 8m23s---

   R23 = 93   R22 = 326215   R21 = 443944   R20 = 152681   R19 = 699238   R18 = 856266   R17 = 700490   R16 = 715968   R15 = 264381
   R14 = 621468   R13 = 592963   R12 = 895217   R11 = 599993   R10 = 229915   R09 = 608941   R08 = 463976   R07 = 156518   R06 = 286253
   R05 = 697920   R04 = 827223   R03 = 758251   R02 = 185210   R01 = 916864

   and   R00 = 24  so we must add  24 zeros at the right of the number above:

 100 ! =  93,326215,443944,152681,699238,856266,700490,715968,264381,621468,592963,895217,599993,229915,608941,463976,156518,
              286253,697920,827223,758251,185210,916864,000000,000000,000000,000000.

Note:

-Add zeros on the left if there are less than 6 digits in a register - except the first one ( for instance read 001234 instead of 1234 )
 

2°)  Real Arguments
 

     a)  Focal Program
 

-The 2 following routines are quite similar to those listed in "Gamma Function for the HP-41"
-They simply use the formula  x ! = Gam(x+1)

-This program employs the continued fraction that s used in the program "GAM3":

       Gam(x) = [ x^(x-1/2) ] sqrt(2.Pi) exp [ -x + ( 1/12 )/( x + ( 1/30 )/( x + ( 53/210)/( x + (195/371)/( x + ...  )))) ]
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "FCT"  02  FRC  03  X#0?  04  GTO 00          05  LASTX  06  FACT  07  RTN  08  LBL 00  09  1  10  LASTX  11  5

12  X<>Y  13  LBL 01   14  1  15  +  16  ST* Z  17  X<Y?  18  GTO 01         19  STO Y  20  .1  21  +  22  210

23  *  24  53  25  X<>Y  26  /  27  RCL Y          28  +  29  30  30  *  31  1/X  32  RCL Y   33  +

34  12  35  *  36  1/X  37  RCL Y          38  -  39  E^X  40  PI  41  ST+ X   42  SQRT  43  R^  44  /

45  *  46  X<>Y  47  R^  48  .5  49  ST+ Y          50  *  51  Y^X  52  ST* Y   53  *  54  END

 
    ( 79 bytes / SIZE 000 )
 
 

STACK	INPUT	OUTPUT
X	x	x !

  where  x is a real number ( integer or not ).

Examples:

   PI    XEQ "FCT"   >>>>   7.188082733
  -1.28    R/S      >>>>    -4.526689379
   69.7    R/S      >>>>      3.343680763 E99
    10      R/S      >>>>      3628800

Notes:

-If n is a negative integer, you'll get DATA ERROR
-If n is a non-negative integer, the built-in function FACT is used and the routine stops at line 07.

-It would be more rigorous to replace line 20 by   371   *   195   X<>Y   /   RCL Y    but this doesn't seem to change the results significantly.
 

     b)  M-Code Routine
 

-The same continued fraction is used again.
-There is no check for alpha data.
 

094  "T"
003  "C"
006  "F"
2A0  SETDEC
0F8  C=X
128  L=C
3C8  CLRKEY
084  C
0ED  =
064  frc(C)
2EE  ?C#0
06F  JC+13d
0F8  C=X
2FE  ?C<0
027  JC+04
3B1  C=
060  fact(C)
02B  JNC+05
04E  C
27A  =
130
099  MAXR
0E8  X=C
3E0  RTN
268  Q=C
138  C=L
070  N=C
10E  A=C
04E  C
35C  =
190  6
0E8  X=C
2BE  C=-C
01D  ?CXQ
061  C=A+C
2FE  ?C<0
05B  JNC+11d
0F8  C=X
10E  A=C
278  C=Q
35C  PT=12
2FE  ?C<0
013  JNC+02
162  A=A+1 @PT
01D  C=
060  A+C
070  N=C ALL
0B0  C=N  ALL
10E  A=C
04E  C
35C
0D0  =
1D0
050  3.71
135  C=
060  A*C
04E  C
35C
050  =
250
150  1.95
269  C=
060  AB/C
239  C=
060  1/AB
0B0  C=N ALL
025  C=
060  AB+C
04E  C
35C
090  =
050
226  21
13D  C=
060  AB*C
239  C=
060  1/AB
04E  C
35C  =
150
0D0  5.3
13D  C=
060  AB*C
0B0  C=N ALL
025  C=
060  AB+C
04E  C
35C  =
0D0
226  30
13D  C=
060  AB*C
239  C=
060  1/AB
0B0  C=N ALL
025  C=
060  AB+C
04E  C
35C
050  =
090
226  12
13D  C=
060  AB*C
239  C=
060  1/AB
0B0  C=N ALL
2BE  C=-C
025  C=
061  AB+C
089  AB
064  STO Q+
04E  C
35C  =
150
266  0.5
10E  A=C
0B0  C=N ALL
01D  C=
060  A+C
0E8  X=C
0B0  C=N ALL
084  C
115  =
06C  LN(C)
0F8  C=X
13D  C=
060  AB*C
0D1  RCL
064  Q+
031  C=
060  AB+CM
044  C
035  =
068  exp(AB)
138  C=L
0E8  X=C
089  AB                          loop
064  STO Q+
0F8  C=X
10E  A=C
0B0  C=N ALL
36E  ?A#C
07B  JNC+15d
3CC  ?KEY                                        Press any key to stop the loop if you want
360  ?C RTN
04E  C
35C  =
050   1
01D  C=
060   A+C
0E8  X=C
0A9  AB
064  <>Q+
0F8  C=X
269  C=
060  AB/C
363  JNC-20d                 goto loop
271  C
064  =
1EE
1EE  2 x PI ( 13 digits )
0EE  B<>C
00E  A=0
305  C=
060  sqrt(AB)
0D1  RCL
064  Q+
149  C=
060  AB*CM
0E8  X=C
3E0  RTN

( 172 words )
 
 

STACK	INPUT	OUTPUT
X	x	x !

  where  x is a real number ( integer or not ).

Examples:

     PI      XEQ "FCT"   >>>>     7.188082733
  -1.28    XEQ "FCT"   >>>>    -4.526689374
   69.7    XEQ "FCT"   >>>>      3.343680766 E99
  -41.7   XEQ "FCT"   >>>>     -3.545755999 E-49                     in less than 5 seconds
    10      XEQ "FCT"   >>>>      3628800
    -3      XEQ "FCT"   >>>>       9.999999999 E99

Notes:

-For all negative integers, FCT returns  9.999999999 E99
-If n is a non-negative integer, the built-in function FACT is used: the results may be exploited provided n < 100

   99  XEQ "FCT"  >>>>  9.332621544 E-45  which is in fact  9.332621544 E155    ( press  ENTER^ LOG  to get the correct exponent )

-If x is not an integer, x must be smaller than about  120.755.....   ( with x = 120.756 the result is wrong )

  120.755  XEQ "FCT"  gives  2.497842660 E-00  which is actually  2.497842660 E200    ( press  ENTER^ LOG  to get the correct exponent )
 

     c)  Log ( n! )
 

-This short program computes Log(n!) where n > -1 and it's especially usefull to get an idea of n! when n is a large number.
-It's almost the same as "LOGAM" listed in "Gamma Function for the HP-41"
-We employ an asymptotic series.
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "LOGFCT"  02  1  03  +  04  LOG  05  STO M             06  16  07  LASTX  08  LBL 01  09  1  10  +

11  LOG  12  ST+ M   13  X<> L  14  X<Y?  15  GTO 01            16  ENTER^  17  X^2  18  SIGN  19  LASTX  20  30

21  *  22  1/X  23  -  24  X<>Y   25  12  26  *  27  /  28  X<>Y               29  -  30  10

31  LN  32  /  33  X<>Y  34  ENTER^           35  LOG  36  *  37  +  38  X<>Y   39  PI  40  *

41  ST+ X  42  SQRT  43  LOG  44  +  45  0  46  X<> M             47  -  48  END

 
   ( 69 bytes / SIZE 000 )
 
 

STACK	INPUT	OUTPUT
X	x	Log (x!)

  where  x > -1.

Examples:

  1000   XEQ "LOGFCT"  >>>>  2567.604644
  1234              R/S            >>>>  3280.708294                              ---Execution time = 3s---
  2345              R/S            >>>>  6886.648594

Notes:

-The results are not very accurate to get n! but the routine is fast !
-In the last example,  FRC  10^X  gives  4.452398215  whence  2345 ! = 4.45239(8215) E6886

-The correct result is  4.45238860677855988 ......  E6886  we have at most 6 significant digits:
-That is logical since  FRC  gave  0.648594  with 6 significant digits too...
 

3°)  Multifactorials
 

     a)  Focal Program
 

-Factorials may be generalized as follows:

  double factorial   n!! = n(n-2)(n-4) .......          all the factors remaining positive.
   triple factorial    n!!! = n(n-3)(n-6) .......

   ..........................................................

   multifactorials  n !....! = n(n-k)(n-2k).....        where k = the number of  "!"s
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "MFCT"  02  1  03  STO T  04  RDN  05  X<=0?  06  GTO 02  07  LBL 01  08  ST* Z  09  RCL Y  10  -  11  X>0?  12  GTO 01  13  LBL 02  14  X<> Z  15  END

 
( 30 bytes / SIZE 000 )
 
 

STACK	INPUTS	OUTPUTS
Y	k	k
X	n	n !!......!

  where n is a non-negative integer and  k = the number of "!"s

Examples:

     2    ENTER^
   100  XEQ "MFCT"  >>>>   100 !! = 3.424322474 E79

    3    ENTER^
  100    R/S        >>>>    100 !!! = 1.745486670 E53

    4    ENTER^
  100    R/S        >>>>    100 !!!! = 1.746406994 E40

Note:

-This program may also be used to compute factorials ( k = 1 )
 

     b)  M-Code Routine
 

-An M-Code routine will give faster and more accurate results.
-This one is very similar to "FCT" listed in §1-c)

-It takes k in Y-register & n in X-register and returns the mantissa of n!....! ( k "!"s ) in X-register and the exponent in Y-register.
-There is no check for alpha data.
-Press any key to stop the routine if it's too slow.
 

094   "T"
003   "C"
006   "F"
00D  "M"
2A0   SETDEC
3C8  CLRKEY
0B8  C=Y
05E  C= | C |
2BE  C=-C
070   N=C ALL                         CPU register N contains  (-k)
0F8   C=X
128   L=C
04E  C=0 ALL
0A8  Y=C
35C   C=
050    1
0E8   X=C
138   C=L
2EE  ?C#0 ALL
3A0  ?NC RTN
268   Q=C                      LOOP
3CC  ?KEY
360   ?C RTN                                              Press a key to stop the routine if you think it's too slow
2EE   ?C#0 ALL
0AB   JNC+21d
2FE   ?C#0 MS
09F   JC+19d
1A0  A=B=C=0
0F8  C=X
0EE  B<>C ALL
278  C=Q
13D  C=
060  AB*C
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
10E  A=C ALL
0B0  C=N ALL
01D  C=
060   A+C
343  JNC-24d                       GOTO LOOP
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

( 65 words )
 
 

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

  where n is a non-negative integer & k is the number of "!"s

Examples:

      5      ENTER^
   1234   XEQ "MFCT"   >>>>   2.641604251   X<>Y    657   so   1234 !!!!! = 2.641604251 E657                           ---Execution time = 12s---

     6       ENTER^
 12345   XEQ "MFCT"   >>>>   7.556280567   X<>Y   7526   so   12345 !!!!!! = 7.556280567 E7526                  ---Execution time = 103s---

Notes:

-Press any key to stop the loop if you wish
-Actually, only | k | is taken into account: k = +5 or k = -5 will give the same output.
-If n or k are not integers, the product n(n-k)(n-2k).....  will be performed until the factor becomes non-positive.

-Of course, "MFCT" may also be used to compute factorials ( with k = 1 ), for instance:

     1      ENTER^
  2345   XEQ "MFCT"  >>>>     4.452388605  X<>Y  6886   whence  2345 ! = 4.452388605 E6886               ---Execution time = 166s---

-The reults are identical to those given by "FCT" listed in paragraph 1-c) ... and the execution times too !
 

     c)  Log ( M(n) )
 

-"LOGMF" calculates  Log (M(n)) = Log n + Log (n-k) + ......
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "LOGMF"  02  STO M  03  0  04  X<>Y  05  LBL 01  06  LOG  07  +  08  RCL M  09  RCL Z  10  -  11  STO M  12  X>0?  13  GTO 01  14  CLX  15  STO M  16  +  17  END

 
( 33 bytes / SIZE 000 )
 
 

STACK	INPUT	OUTPUT
Y	k	k
X	n > 0	Log M(n)

  where  n is a positive integer

Examples:

     2      ENTER^
   100    XEQ "LOGMF"   >>>>   Log ( 100 !! )  =  79.53457468                     ---Execution time = 20s---

   999     ENTER^
123456  XEQ "LOGMF"  >>>>   Log ( 123456 ! ..... ! ) = 578.0564932           ---Execution time = 53s---

Notes:

-In the last example, there are 999  "!"s
-We can write an M-Code routine for a better precision:
 

086   "F"
00D  "M"
00C  "L"
2A0  SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
1A0  A=B=C=0
0E8   X=C
089   AB
064   STO Q+
138   C=L
070   N=C ALL            Loop
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
3A0   ?NC RTN
2FE   ?C#0 MS
360    ?C RTN
088   C
115   =
06C  Log(C)
0D1  RCL
064  Q+
031   C=
060  AB+CM
0E8  X=C
089  AB
064  STO Q+
0B8  C=Y
2BE  C=-C
10E   A=C ALL
0B0  C=N ALL
01D  C=
060  A+C
34B  JNC-23d                    goto Loop

( 36 words )
 
 

STACK	INPUTS	OUTPUTS
Y	k > 0	k
X	n > 0	Log M(n)
L	/	n

  where  n & k positive integers

Examples:

     2      ENTER^
   100    XEQ "LMF"   >>>>   Log ( 100 !! )  =  79.53457466                     ---Execution time = 13s---

   999     ENTER^
123456  XEQ "LMF"  >>>>   Log ( 123456 ! ..... ! ) = 578.0564933           ---Execution time = 35s---

Note:

-With k = 1,  "LMF" computes Log (n!)
 

     d)  Long Integers
 

-Like "LFCT", "LMFCT" computes M(n) by groups of 6 digits
 

Data Registers:            R00 = number of zeros to be added to the result in R01-R02 ................

                                      Rmm , Rm-1 , Rm-2 , ............... , R01 =  the digits of  M(n)
Flags: /
Subroutines: /
 
 

01  LBL "LMFCT"  02  STO M  03  X<>Y  04  STO N  05  1.001  06  STO O  07  CLX  08  STO 00  09  SIGN  10  STO 01            11  LBL 01  12  RCL O  13  RCL M  14  LBL 02  15  ST* IND Y  16  ISG Y  17  GTO 02  18  RCL O

19  LBL 03  20  RCL IND X  21   E6  22  X>Y?  23  GTO 05            24  MOD  25  ST- IND Y  26  X<> IND Y  27  LASTX  28  /  29  ISG Y  30  GTO 04  31  STO IND Y  32  X<>Y  33  INT  34   E3  35  /  36  ISG X

37  STO O  38  GTO 06  39  LBL 04  40  ST+ IND Y  41  ENTER^  42  DSE Z  43  LBL 05            44  X<> Z  45  ISG X  46  GTO 03  47  LBL 06  48  RCL 01  49  X#0?  50  GTO 07  51  6  52  ST+ 00  53  RCL O  54  FRC

55   E3  56  *  57  DSE X  58  STO Y  59  LASTX  60  X^2  61  /  62  2.001  63  +  64  REGMOVE  65  CLX  66   E3  67  /  68  ISG X  69  STO O  70  LBL 07            71  CLX  72  SIGN

73  RCL M  74  RCL N  75  -  76  STO M  77  X>Y?  78  GTO 01           79  RCL O  80  FRC  81   E3  82  ST* Y  83  1/X  84  +  85  CLA  86  END

 
    ( 147 bytes / SIZE ??? )
 
 

STACK	INPUTS	OUTPUTS
Y	k > 0	/
X	n > 0	mmm.001

  where  n & k positive integers ( n < 10000 ) and mmm.001 is the control number of the result

Example:

   1234   ENTER^
   9999   XEQ "LMFCT"  >>>>  6.001   and  R00 = 0  ( no zero must be added to the result )

-We have  R06 = 36    R05 = 139932    R04 = 60227    R03 = 617285    R02 = 723555    R01 = 932975

-Thus,  9999 !!....! ( 1234 "!"s ) = 36,139932,060227,617285,723555,932975   ( approximately  3.613993206 E31 )

Notes:

-If there are less than 6 digits, add zeros on the left ( for instance, read  R04 = 060227 )
-With k = 1 , "LMFCT" will calculate n !
 

4°)  Superfactorials
 

     a)  Focal Program
 

-'SFCT" calculates Sf(n) = n! (n-1)! (n-2)! ....... 2! 1! 0!
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "SFCT"  02  SIGN  03  LASTX  04  LBL 01  05  FACT  06  ST* Y  07  X<> L  08  DSE X  09  GTO 01  10  RDN  11  END

 
( 24 bytes / SIZE 000 )
 
 

STACK	INPUT	OUTPUT
X	n	Sf(n)

  where  n is a non-negative integer & Sf(n) =  n! (n-1)! (n-2)! ....... 2! 1! 0!

Examples:

  0   XEQ "SFCT"   >>>>   Sf(0) =   1
  4          R/S           >>>>   Sf(4) =  288
 16         R/S           >>>>  Sf(16) = 1.890966834 E90

-For n = 17 , there is an OUT OF RANGE

Note:

-Clifford A. Pickover also defines another superfactorial in "Keys to Infinity" page 102 ( reference 6 )  n$ = n!^...^n!^n!  where n!  is repeated n! times.
-They seem very difficult to compute except for very small n-values...
 

     b)  M-Code Routine
 

-We could modify the M-Code routine "FCT" but the following program ( "SMFCT" ) combines "MFCT" & "SFCT"
-So, it can be used to compute factorials ( k = 1 ) , multifactorials ( k > 1 ) and superfactorials ( k = -1 )
 

094   "T"
003   "C"
006   "F"
00D  "M"
013   "S"
2A0   SETDEC
3C8  CLRKEY
104  CLRF 8
0B8  C=Y
2FE  ?C#0 MS
013   JNC+02
108   SETF 8
05E  C= | C |
2BE  C=-C
070   N=C ALL                         CPU register N contains  - | k |
0F8   C=X
128   L=C
04E  C=0 ALL
0A8  Y=C
35C   C=
050    1
0E8   X=C
138   C=L
2EE  ?C#0 ALL
12B  JNC+37d
268   Q=C                      LOOP
3CC  ?KEY
360   ?C RTN                                              Press a key to stop the routine if you think it's too slow
2EE   ?C#0 ALL
0AB   JNC+21d
2FE   ?C#0 MS
09F   JC+19d
1A0  A=B=C=0
0F8  C=X
0EE  B<>C ALL
278  C=Q
13D  C=
060  AB*C
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
10E  A=C ALL
0B0  C=N ALL
01D  C=
060   A+C
343  JNC-24d                       GOTO LOOP
10C  ?FSET 8
053   JNC+10d
138   C=L
02E   C
0FA  ->
0AE  AB
009   C=
060   AB-1
128   L=C
2EE  ?C#0 ALL
2EF   JC-35d                         GOTO LOOP
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

( 81 words )
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	Z	Z
Y	k	exponent
X	n	mantissa

  where n is a non-negative integer & k is the number of "!"s  for multifactorials.   k = -1  for superfactorials

Examples:

     1      ENTER^
  2345   XEQ "SMFCT"  >>>>     4.452388605  X<>Y   6886   whence  2345 ! = 4.452388605 E6886               ---Execution time = 111s---

     6       ENTER^
 12345   XEQ "SMFCT"   >>>>   7.556280567   X<>Y   7526   so   12345 !!!!!! = 7.556280567 E7526             ---Execution time = 103s---

   -1     ENTER^
   41    XEQ "SMFCT"   >>>>      4.885831876   X<>Y   873  whence  Sf(41) = 4.885831876  E873                   ---Execution time = 36s---

   -1     ENTER^
  100   XEQ "SMFCT"   >>>>      2.703176855   X<>Y   6940   thus    Sf(100) = 2.703176855  E6940               ---Execution time = 214s---
 

Note:

-"SMFCT" may also be used with k = -2 , -3 , .... to calculate  n!! (n-1)!! (n-2)!! ...... and   n!!! (n-1)!!! (n-2)!!! ......   though I don't know if it's really usefull.
 

     c)  Log ( Sf(n) )
 

-"LOGSF" calculates  Log (Sf(n)) = Log [ n! (n-1)! (n-2)! ....... 2! 1! ]
-It uses the same formula as "LOGFCT" in paragraph 2-c)
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "LOGSF"  02  CLA  03  STO N   04  LBL 10        05  RCL N            06  1  07  +  08  LOG  09  ST- M  10  16  11  LASTX

12  LBL 01   13  1  14  +  15  LOG  16  ST- M  17  X<> L             18  X<Y?  19  GTO 01   20  ENTER^  21  X^2  22  SIGN

23  LASTX   24  30  25  *  26  1/X  27  -  28  X<>Y             29  12  30  *  31  /  32  X<>Y  33  -

34  10  35  LN  36  /  37  X<>Y  38  ENTER^          39  LOG  40  *  41  +  42  X<>Y  43  PI  44  *

45  ST+ X  46  SQRT  47  LOG  48  +  49  ST+ M  50  DSE N  51  GTO 10           52  X<> M  53  CLA  54  END

 
    ( 79 bytes / SIZE 000 )
 
 

STACK	INPUT	OUTPUT
X	n	Log [ Sf(n) ]

  where  n is a non-negative integer & Sf(n) =  n! (n-1)! (n-2)! ....... 2! 1! 0!

Example:

       10   XEQ "LOGSF"  >>>>   Log (Sf(10)) = 27.82338353             ---Execution time = 65s---

Notes:

-The last 2 decimals are wrong
-A better formua - with more terms in the asymptotic series - would give of course a better accuracy...

-The variant hereunder employs a more elementary method which is preferable if n is not too large:

  Log Sf(n) = 1 Log n + 2 Log (n-1) + 3 Log (n-2) + ...... + n Log 1
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "LOGSF"  02  STO M  03  SIGN  04  0  05  LBL 01  06  RCL M  07  LOG  08  RCL Z  09  *  10  +  11  ISG Y  12  CLX  13  DSE M  14  GTO 01  15  END

 
( 31 bytes / SIZE 000 )
 
 

STACK	INPUT	OUTPUT
X	n > 0	Log Sf(n)

  where  n is a positive integer

Examples:

    10   XEQ "LOGSF"  >>>>   Log (Sf(10)) = 27.82338336             ---Execution time = 5s---

   100  XEQ "LOGSF"  >>>>   Log (Sf(100)) = 6940.431879           ---Execution time = 54s---

Note:

-The following M-Code routine uses the same method:
 

086   "F"
013   "S"
00C  "L"
2A0  SETDEC
3C8  CLRKEY
0F8   C=X
070   N=C ALL
1A0  A=B=C=0
128   L=C
089   AB
064   STO Q+
138   C=L                              Loop
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
02E   C
0FA  ->
0AE  AB
001   C=
060   AB+1
128   L=C
088   C
115   =
06C  Log(C)
0B0  C=N ALL
13D  C=
060   AB*C
0D1  RCL
064  Q+
031   C=
060  AB+CM
0E8  X=C
089  AB
064  STO Q+
0B0  C=N ALL
02E   C
0FA  ->
0AE  AB
009  C=
060  AB-1
070  N=C ALL
2EE  ?C#0 ALL
317  JNC-30d                           goto Loop
3E0  RTN

( 43 words )
 
 

STACK	INPUT	OUTPUT
X	n > 0	Log Sf(n)
L	/	n

  where  n is a positive integer

Examples:

    10    XEQ "LSF"  >>>>   Log (Sf(10)) = 27.82338336             ---Execution time = 3s---

   100   XEQ "LSF"  >>>>   Log (Sf(100)) = 6940.431874           ---Execution time = 31s---

  1000  XEQ "LSF"  >>>>   Log (Sf(1000)) = 1177245.511         ---Execution time = 5m24s---
 

     d)  Long Integers
 

-If we want all the digits, several data registers must be used to store them.

-"LSFCT" calculates Sf(n)  by groups of 6 digits.  The REGMOVE function from the X-functions module is used.
 

Data Registers:            R00 = number of zeros to be added to the result in R01-R02 ................

                                      Rmm , Rm-1 , Rm-2 , ............... , R01 =  the digits of Sf(n)
Flags: /
Subroutines: /

-Line 90 is a three-byte GTO 10
 
 

01  LBL "LSFCT"  02  STO 01   03   E6  04  STO N  05  1.001  06  STO O   07  CLX  08  STO 00   09  SIGN  10  X<> 01  11  LBL 10   12  RCL N  13  SQRT  14  /  15  2  16  +  17  STO M            18  LBL 01  19  RCL O  20  RCL M

21  INT  22  LBL 02  23  ST* IND Y  24  ISG Y  25  GTO 02   26  RCL O   27  LBL 03  28  RCL IND X  29  RCL N  30  X>Y?  31  GTO 05   32  MOD  33  ST- IND Y  34  X<> IND Y  35  LASTX  36  /  37  ISG Y  38  GTO 04           39  STO IND Y  40  X<>Y

41  INT  42  RCL N  43  SQRT  44  /  45  ISG X  46  STO O   47  GTO 06   48  LBL 04  49  ST+ IND Y  50  ENTER^  51  DSE Z  52  LBL 05   53  X<> Z  54  ISG X  55  GTO 03           56  LBL 06  57  RCL 01  58  X#0?  59  GTO 07  60  6

61  ST+ 00  62  RCL O   63  FRC  64  RCL N   65  SQRT  66  *  67  DSE X  68  STO Y  69  RCL N   70  /  71  2.001  72  +  73  REGMOVE  74  CLX  75  RCL N             76  SQRT  77  /  78  ISG X  79  STO O   80  LBL 07

81  ISG M  82  GTO 01   83  CLX  84  SIGN  85  RCL M   86  INT  87  2  88  -  89  X>Y?  90  GTO 10  91  RCL O            92  FRC  93   E3  94  ST* Y  95  1/X  96  +  97  CLA  98  END

 
    ( 164 bytes / SIZE ??? )
 
 

STACK	INPUT	OUTPUT
X	n	mmm.001

  where  n is an integer ( 1 < n < 1000 ) and  mmm.001  is the control number of the result

Examples:

   •  10  XEQ "LSFCT"  >>>>   4.001                               ---Execution time = 87s---

   R04 = 6658    R03 = 606584    R02 = 104736    R01 = 522240  &  R00 = 6 zeros to add on the right

-Whence   Sf(10) = 6658,606584,104736,522240,000000    ( approximately  6.658606584 E27 )

  •  41  XEQ "LSFCT"  >>>>   118.001                            ---Execution time = several hours---   ( V41 in turbo mode is preferable )

   R118 = 4885   R117 = 831876  .......................  R02 = 724851    R01 = 200000        R00 = 168

-You can check that  Sf(41) =

  4885831876911199922535923458596922414263022771295934189443274148270199837389494004711909972672987233752
  5619370340788379771991357781209553324584667109104833395981620779095261652107715556366032754444863043260
  0869714329763066401198372444861903682765000887343163113829267580062367575521410631749338487229360389773
  2239349531719939588145981827516350700200493273745612968699958924658210871599738389102330295346404136128
  8323220093103506346463801828526038711001495475735727877014022882457875749302679561878377839682128010664
  1374901877816660520680981245203287001765885516903031502642011575504076480041176841698560314015530446023
  8022531783028878525746632438552494719100713807005964600719820199559851929303724851200000000000000000000
  0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
  00000000000000000000000000000000000000000000000000.
 

Notes:

-Add zeros on the left if there are not enough digits in a register ( for instance; read  001234  instead of  1234 )
-Though the control numbers could handle n = 999, there are not enough registers to store all the digits of  Sf(999)
 

5°)  Hyperfactorials
 

     a)  Focal Program
 

-'HFCT" calculates H(n) = nⁿ (n-1)^n-1  ....... 2² 1¹
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "HFCT"  02  1  03  LBL 01  04  RCL Y  05  ENTER^  06  Y^X  07  *  08  DSE Y  09  GTO 01  10  END

 
( 22 bytes / SIZE 000 )
 
 

STACK	INPUT	OUTPUT
X	n	H(n)

  where  n is a non-negative integer & H(n) =  nⁿ (n-1)^n-1  ....... 2² 1¹

Examples:

  1   XEQ "HFCT"   >>>>   H(1)  =   1
  4          R/S           >>>>   H(4)   =  27648
 14         R/S           >>>>   H(14) =  1.847398448 E99                   ---Execution time = 8s---

Note:

  n = 15  produces an OUT OF RANGE
 

     b)  3 M-Code Routines
 

-"HFCT" computes the hyperfactorial of a non-negative integer with 13-digit routines.
-Synthetic register P is cleared ( it could be replaced by register Z or another stack-register T or L )

-There is no check for alpha data.
 

094   "T"
003   "C"
006   "F"
008   "H"
2A0   SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
04E   C=0 ALL
0A8  Y=C
35C   C=
050   1
0E8   X=C
138   C=L
070   N=C ALL
2EE  ?C#0 ALL
3A0  ?NC RTN
006   A=0 S&X              the 8 words from this one shift the digits of n into the right part of the mantissa field of register C
39C  PT=0                      then, it's stored into synthetic registers P & Q
1A2  A=A-1@PT
362   ?A#C @PT            Thus, it's simpler and faster to decrement these numbers by 1: we simply have to use  27A  C=C-1 M  to perform k-1 multiplications.
01B   JNC+03
3DA  RSHFC M
3E3   JNC-04
046   C=0 S&X
27A  C=C-1 M
268   Q=C                       Loop1
228   P=C                        Loop2
3CC  ?KEY
360   ?C RTN                                              Press a key to stop the routine if you think it's too slow
1A0  A=B=C=0
0F8  C=X
0EE  B<>C ALL
0B0  C=N ALL
13D  C=
060   AB*C
0EE  B<>C ALL
0E8  X=C
01A  A=0 M
0B8  C=Y
20E  C=A+C ALL
0A8  Y=C
238   C=P
27A  C=C-1 M
37B  JNC-17d                goto loop2
0B0  C=N ALL
02E   C
0FA  ->
0AE  AB
009   C=
060   AB-1
070  N=C ALL
278  C=Q
27A  C=C-1 M
323   JNC-28d                goto  loop1
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

( 75 words )
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	Z	Z
Y	Y	exponent
X	n	mantissa
L	L	n

  where n is a non-negative integer

Examples:

   14    XEQ "HFCT"   gives    1.847398449   X<>Y   99     whence   H(14) =  1.847398449  E99                     ---Execution time = 3s---
  100   XEQ "HFCT"   yields   3.455370850   X<>Y  9014    thus      H(100) = 3.455370850  E9014                ---Execution time = 148s---
  234   XEQ "HFCT"   yields   6.891346366   X<>Y  59196   so       H(234) = 6.891346366  E59196               ---Execution time = 13m08s---

Note:

-The following variant avoids using synthetic register P
-It's slightly slower but it saves 3 words:
 

094   "T"
003   "C"
006   "F"
008   "H"
2A0   SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
04E   C=0 ALL
0A8  Y=C
35C   C=
050   1
0E8   X=C
138   C=L
070   N=C ALL
2EE  ?C#0 ALL
3A0  ?NC RTN
006   A=0 S&X              Loop1
39C  PT=0
1A2  A=A-1@PT
362   ?A#C @PT
01B   JNC+03
3DA  RSHFC M
3E3   JNC-04
27A  C=C-1 M
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
1A0  A=B=C=0
0F8  C=X
0EE  B<>C ALL
0B0  C=N ALL
13D  C=
060   AB*C
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
27A  C=C-1 M
37B  JNC-17d                goto loop2
0B0  C=N ALL
02E   C
0FA  ->
0AE  AB
009   C=
060   AB-1
070   N=C ALL
2EE  ?C#0 ALL
2F7   JNC-34d                goto  loop1
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

( 72 words )
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	Z	Z
Y	Y	exponent
X	n	mantissa
L	L	n

  where n is a non-negative integer

Example:   ( same results as above )

  100   XEQ "HFCT"   yields   3.455370850   X<>Y  9014    thus      H(100) = 3.455370850  E9014                ---Execution time = 149s---

Note:

-The 2 routines listed above calculate  m^m  by performing m-1 multiplications.
-The following one performs much less operations, at least if n is not too small.
-So, we can expect faster and more accurate results.

-Synthetic registers O & P are cleared, except for n = 0.
-There is no check for alpha data.
 

094   "T"
003   "C"
006   "F"
008   "H"
2A0   SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
04E   C=0 ALL
0A8  Y=C
35C   C=
050   1
0E8   X=C
138   C=L
2EE  ?C#0 ALL
3A0  ?NC RTN
268   Q=C                     Loop1
070   N=C ALL
05A  C=0 M
228   P=C
0B0  C=N ALL
046   C=0 S&X
1E8   O=C
0B0   C=N ALL            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
10E   A=C ALL
04E   C
35C  =
090   2
261   C=
060   A/C
070   N=C ALL
084   C
0ED  =
064   frc(C)
2BE  C=-C
10E  A=C ALL
0F0  C<>N ALL
01D  C=
060   A+C
0F0  C<>N ALL
2FA  ?C#0 M
083   JNC+16d
1A0  A=B=C=0
1F8  C=O
0EE  B<>C 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
14E  A=A+C ALL
238  C=P
20E  C=A+C ALL
0A8  Y=C
0B0  C=N ALL
2EE  ?C#0 ALL
08B  JNC+17d
1A0  A=B=C=0
1F8  C=O
0EE  B<>C ALL
158  M=C ALL
1F8  C=O
149  C=
060  AB*CM
0EE  B<>C ALL
1E8  O=C
01A  A=0 M
238  C=P
1EE  C=C+C ALL
20E  C=A+C ALL
228  P=C
25B  JNC-53d                         goto loop2
21B  JNC-61d                         goto loop1
278  C=Q
02E   C
0FA  ->
0AE  AB
009   C=
060   AB-1
2EE  ?C#0 ALL
3C7  JC-08
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
05B  JNC+11d
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
04E   C=0 ALL
1E8   O=C
228   P=C
3E0   RTN

( 109 words )
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	Z	Z
Y	Y	exponent
X	n	mantissa
L	L	n

  where n is a non-negative integer

Examples:

   14    XEQ "HFCT"   gives    1.847398449   X<>Y   99     whence   H(14)  = 1.847398449  E99                     ---Execution time = 4.6s---
  100   XEQ "HFCT"   yields   3.455370850   X<>Y  9014    thus      H(100) = 3.455370853  E9014                ---Execution time = 64s---
  234   XEQ "HFCT"   yields   6.891346366   X<>Y  59196   so       H(234) = 6.891346393  E59196               ---Execution time = 186s---

Notes:

-Though the routine is slightly slower for n = 14, we save more than 1 minute for n = 100 and more than 10 minutes for n = 234.
-For n = 100 & 234, roundoff-errors are also smaller.

-There will be an infinite loop if n is negative or fractional: press any key to stop the loop.
 

     c)  Log ( H(n) )
 

-"LOGHF" simply uses the formula::

   Log ( H(n) ) = n Log n + ( n-1 ) Log ( n-1 ) + ................. + 2 Log 2 + 1 Log 1
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "LOGHF"  02  0  03  LBL 01  04  RCL Y  05  ENTER^  06  LOG  07  *  08  +  09  DSE Y  10  GTO 01  11  END

 
( 24 bytes / SIZE 000 )
 
 

STACK	INPUT	OUTPUT
X	n > 0	Log H(n)

  where  n is a positive integer

Examples:

  4   XEQ "LOGHF"   >>>>     Log H(4)   =   4.441663720
 14           R/S             >>>>     Log H(14)  =  99.26656058                   ---Execution time = 7s---
100          R/S             >>>>    Log H(100) =  9014.538503                   ---Execution time = 47s---
 

Note:

-The following M-Code routine uses the same method:
 

086   "F"
008   "H"
00C  "L"
2A0  SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
1A0  A=B=C=0
0E8   X=C
089   AB
064   STO Q+
138   C=L
070   N=C ALL                             Loop
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
3A0  ?NC RTN
088   C
115   =
06C  Log(C)
0B0  C=N ALL
13D  C=
060   AB*C
0D1  RCL
064  Q+
031   C=
060  AB+CM
0E8  X=C
089  AB
064  STO Q+
0B0  C=N ALL
02E   C
0FA  ->
0AE  AB
009  C=
060  AB-1
343  JNC-24d                goto Loop

( 37 words )
 
 

STACK	INPUT	OUTPUT
X	n	Log H(n)
L	/	n

  where  n is a non-negative integer

Examples:

   100   XEQ "LHF"  >>>>   Log H(100)  =  9014.538495          ---Execution time = 29s---
   234   XEQ "LHF"  >>>>   Log H(234)  =  59196.83830          ---Execution time = 71s---
  1000  XEQ "LHF"  >>>>   Log H(1000) = 1392926.737          ---Execution time = 5m09s---
 

     d)  Long Integers
 

-"LHFCT" calculates H(n)  by groups of 6 digits.  The REGMOVE function from the X-functions module is used.
 

Data Registers:            R00 = number of zeros to be added to the result in R01-R02 ................

                                      Rmm , Rm-1 , Rm-2 , ............... , R01 =  the digits of H(n)
Flags: /
Subroutines: /
 
 

01  LBL "LHFCT"  02   E3  03  /  04  ISG X  05  STO N   06  1.001  07  STO O   08  CLX  09  STO 00  10  SIGN  11  STO 01   12  LBL 10  13  RCL N  14  STO M           15  LBL 01  16  RCL O  17  RCL M

18  INT  19  LBL 02  20  ST* IND Y  21  ISG Y  22  GTO 02  23  RCL O   24  LBL 03   25  RCL IND X  26   E6  27  X>Y?  28  GTO 05           29  MOD  30  ST- IND Y  31  X<> IND Y  32  LASTX  33  /  34  ISG Y

35  GTO 04  36  STO IND Y  37  X<>Y  38  INT  39   E3  40  /  41  ISG X  42  STO O   43  GTO 06           44  LBL 04   45  ST+ IND Y  46  ENTER^  47  DSE Z  48  LBL 05  49  X<> Z  50  ISG X  51  GTO 03

52  LBL 06  53  RCL 01  54  X#0?  55  GTO 07  56  6  57  ST+ 00  58  RCL O            59  FRC  60   E3  61  *  62  DSE X  63  LASTX  64  /  65  STO O  66  LASTX  67  /  68  2.001

69  +  70  REGMOVE  71  ISG O  72  LBL 07  73  ISG M  74  GTO 01   75  ISG N  76  GTO 10  77  RCL O            78  FRC  79   E3  80  ST* Y  81  1/X  82  +  83  CLA  84  END

 
    ( 146 bytes / SIZE ??? )
 
 

STACK	INPUT	OUTPUT
X	n	mmm.001

  where  n is an integer ( 0 < n < 1000 ) and  mmm.001  is the control number of the result

Example:

   •  10  XEQ "LHFCT"  >>>>   6.001                                              ---Execution time = 144s---

   R06 = 215                 R03 = 562091             and    R00 = 12   ( 12 zeros must be added on the right )
   R05 = 779412           R02 = 680268
   R04 = 229418           R01 = 288000

-Whence  H(10) = 215,779412,229418,562091,680268,288000,000000,000000.

-Approximately  H(10) = 2.157794122 E44
 

References:

[1]   N.J.A. Sloane's On-Line Encyclopedia of Integer Sequences: www.research.att.com/~njas/sequences
[2]   http://mathworld.wolfram.com/Multifactorial.html
[3]   http://oeis.org/A000178
[4]   http://mathworld.wolfram.com/Hyperfactorial.html
[5]   http://oeis.org/A002109
[6]   Clifford A. Pickover - "Keys to Infinity" - John Wiley & Sons -  ISBN  0-471-11857-5