hp41programs

Factorials

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) = nn (n-1)n-1  ....... 22 11
 

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) =  nn (n-1)n-1  ....... 22 11

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