hp41programs

Apery

Apery Numbers for the HP-41


Overview
 

 1°)  Apéry's Constant
 2°)  Apéry Numbers

   a)  3 Focal Programs
   b)  2 M-Code Routines
 

1°)  Apéry's Constant
 

-This program calculates the decimals of Zeta(3) = SUMk=1,2,......  1/k3
-In 1979, the French mathematician Roger Apéry  ( University of Caen in Normandy ) proved that zeta(3) is irrational.
-Since then, Zeta(3) is called Apéry's Constant.

-The following routine employs the series:  Zeta(3) = (5/2) SUMk=1,2,...... (-1)k+1 (k!)2 / k3 / (2k)!
-More explicitly, it uses

        Zeta(3) = u1 - u2 + u3 - .......   with   u1 = 5/4  and   uk+1 / uk = k3 / ( 4k+2) / ( k+1)2
 
 

Data Registers:   R00 = k

                              R01-R02- .............. -Rnn = the digits of zeta(3) by groups of 6 digits  ( except R01 which contains a 7-digit number )
                              Rn+1- .................... -R2n = the successive uj

Flag:  F01 is cleared at the end
Subroutines: /

-Line 59 is a synthetic TEXT0 which may be replaced by another NOP like  STO X   LBL 00  ....
-Synthetic registers M N O Q are used.
 
 

  01  LBL "APERY"
  02  CLRG
  03  STO O          
  04   E6
  05  STO N 
  06  SQRT
  07  /
  08  ST+ O
  09  ISG X
  10  ST+ O
  11  STO M
  12  125 E4
  13  STO 01
  14  STO IND O
  15  SIGN
  16  STO 00
  17  CF 01
  18  LBL 01
  19  FC?C 01
  20  SF 01
  21  RCL O
  22  RCL 00
  23  LBL 02
  24  ST* IND Y
  25  ISG Y
  26  GTO 02
  27  RCL 00
  28  ISG X
  29  CLX
  30  STO Q          
  31  RCL O          
  32  XEQ 06
  33  RCL O
  34  RCL 00
  35  LBL 03
  36  ST* IND Y
  37  ISG Y
  38  GTO 03
  39  RCL O
  40  XEQ 06
  41  RCL O
  42  RCL 00
  43  LBL 04
  44  ST* IND Y
  45  ISG Y
  46  GTO 04
  47  RCL 00
  48  ISG 00
  49  CLX
  50  RCL 00
  51  +
  52  ST+ X
  53  STO Q
  54  RCL O          
  55  XEQ 06
  56  RCL IND O
  57  X=0?
  58  ISG O
  59  ""
  60  X=0?
  61  ISG M
  62  X<0?
  63  GTO 07
  64  RCL O
  65  RCL M
  66  LBL 05
  67  RCL IND Y
  68  FS? 01
  69  CHS
  70  ST+ IND Y
  71  ISG Y
  72  RDN
  73  ISG Y
  74  GTO 05
  75  GTO 01
  76  LBL 06
  77  RCL IND X
  78  RCL Q          
  79  MOD
  80  ST- IND Y
  81  LASTX
  82  ST/ IND Z
  83  CLX
  84  RCL N
  85  *
  86  ISG Y
  87  X<0?
  88  RTN
  89  ST+ IND Y
  90  X<>Y
  91  GTO 06
  92  LBL 07
  93  CF 01
  94  RCL M
  95  INT
  96  DSE X
  97  STO M
  98  LBL 08
  99  RCL IND X
100  RCL N          
101  MOD
102  ST- IND Y
103  X<> IND Y
104  LASTX
105  /
106  DSE Y
107  ST+ IND Y
108  CLX
109  SIGN
110  X<>Y
111  X#Y?
112  GTO 08
113  RCL M
114   E3
115  /
116  ISG X
117  CLA
118  END

 
    ( 204 bytes / SIZE 2n+1 )
 
 

      STACK         INPUT      OUTPUT
           X          n > 1         1.nnn

  where the 6n decimals ( 6n+1 digits ) of zeta(3) are stored in R01 thru Rnn

Example:

    7   XEQ "APERY"  >>>>   1.007                                  ---Execution time = 9mn40s---

-We find in registers R01 thru R07:

    R01 = 1202056       R04 = 399738     R07 = 764985
    R02 =  903159        R05 = 161511
    R03 =  594285        R06 = 449990

-Whence   Zeta(3) = 1.202056,903159,594285,399738,161511,449990,764985..........................

-In fact, the last decimal ( 5 ) should be a "6"

Notes:

-Add zeros on the left if there are less than 6 digits in a register ( read for instance 001234 instead of 1234 )
-In the example above, R00 = 63 so the last computed term was u63

-Since there are at most 319 registers, only 159 registers may be used to store the digits of zeta(3)
-So, only 954 decimals may be computed, the last one(s) beeing doubtful
-However, since we use an alternating series, the roundoff-errors should remain small.

-I've used this program to compute 600 decimals of zeta(3)    100  XEQ "APERY" >>>>  1.100  ( about 3 minutes with V41 in turbo mode )
-ALL the decimals are correct !
-R00 = 985 so the last uk was u985

-With n = 159 , the last k-value should be about 1572.
-So this does not exceed the limit ( 9999 ) to perform correctly the multiplications and divisions i-e without exceeding E10.
 

2°)  Apéry's Numbers
 

     a) 3 Focal Programs
 

-Apéry's Numbers are defined by  An = SUMk=0,1,...,n  (n+k)!2 / k!4  / (n-k)!2                   n = 0 , 1 , 2 , ..................

-The fisrt few are  1  5  73  1445  33001  ....     Sloane's  A005259

-These numbers may also be computed by the formula  An4F3 ( -n , -n , n+1 , n+1 ; 1 , 1 , 1 ; 1 )   where  4F3 is a generalized hypergeometric function.
-The following program uses this method
 

Data Registers:   R00 is unused if you employ the M-Code routine "HGF+"  R01 to R07: temp
Flags: /
Subroutine:   The M-Code routine  "HGF+"  or the focal program  "GHGF"  ( cf "Hypergeometric funtions for the HP-41" )

-Line 17 may be replaced by  XEQ "GHGF"
 
 

 01  LBL "APNB"
 02  CHS
 03  STO 01
 04  STO 02
 05  1
 06  STO 05
 07  STO 06
 08  STO 07
 09  X<>Y
 10  -
 11  STO 03
 12  STO 04
 13  4
 14  PI
 15  INT
 16  1
 17  HGF+
 18  END

 
( 28 bytes / SIZE 008 )
 
 

      STACK        INPUT      OUTPUT
           X             n             An

  where n is a non-negative integer

Examples:

   41   XEQ "APNB"  >>>>   A41 = 4.944386782 E59                      ---Execution time = 22s---
  100  XEQ "APNB"  >>>>  A100 = 2.824655679 E149                    ---Execution time = 48s---
  329  XEQ "APNB"  >>>>  A329 = 1.990511251 E499                    ---Execution time = 147s---

Notes:

-In the last 2 examples, the exponent appears in an ambiguous display:  press  ENTER^    LOG   INT   to get the correct result.
-If n > 329 the results are wrong.

-The following variant employs the formula  An = SUMk=0,1,...,n  (n+k)!2 / k!4  / (n-k)!2
 

Data Registers:   R00 = n ,  R01 thru R03: temp
Flags: /
Subroutines: /
 
 

 01  LBL "APNB"
 02  STO 00        
 03  SIGN
 04  STO 01
 05  STO 02
 06  STO 03
 07  FIX 0
 08  LBL 01
 09  SIGN
 10  RCL 00        
 11  ST+ Y
 12  RCL 02
 13  ST- Z
 14  +
 15  *
 16  RCL 02        
 17  X^2
 18  /
 19  X^2
 20  RCL 03
 21  *
 22  RND
 23  STO 03        
 24  ST+ 01
 25  RCL 00
 26  ISG 02
 27  CLX
 28  RCL 02
 29  X<=Y?
 30  GTO 01
 31  RCL 01        
 32  FIX 9
 33  END

 
    ( 49 bytes / SIZE 004 )
 
 

      STACK        INPUT      OUTPUT
           X             n             An

  where n is a non-negative integer

Example:

   67   XEQ "APNB"  >>>>   A67 = 1.529563779 E99                      ---Execution time = 51s---

Notes:

-Line 22  RND  may be replaced by the M-Code routine  RND0.  In this case, delete lines 32 & 07

-There is an OUT OF RANGE if  n > 67
-The routine hereunder overcomes this limitation:  it returns the mantissa in X and the exponent in Y.
 

Data Registers:   R00 = n ,  R01 thru R05: temp
Flag:  F10 is cleared
Subroutines: /
 
 

 01  LBL "APNB"
 02  STO 00         
 03  SIGN
 04  STO 01
 05  STO 02
 06  STO 03
 07  CLX
 08  STO 04
 09   E80
 10  STO 05
 11  FIX 0
 12  SF 10
 13  LBL 01
 14  SIGN
 15  RCL 00         
 16  ST+ Y
 17  RCL 02
 18  ST- Z
 19  +
 20  *
 21  RCL 02
 22  X^2
 23  /
 24  X^2
 25  RCL 03
 26  *
 27  FS? 10
 28  RND
 29  STO 03         
 30  ST+ 01
 31  RCL 05
 32  X>Y?
 33  GTO 00
 34  ST/ 01
 35  ST/ 03
 36  LOG
 37  ST+ 04
 38  CF 10
 39  LBL 00
 40  RCL 00         
 41  ISG 02
 42  CLX
 43  RCL 02
 44  X<=Y?
 45  GTO 01
 46  RCL 01
 47  ENTER^
 48  LOG
 49  INT
 50  ST+ 04
 51  10^X
 52  /
 53  RCL 04         
 54  X<>Y
 55  FIX 9
 56  CF 10
 57  END

 
     ( 84 bytes / SIZE 006 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             /       exponent
           X             n       mantissa

  where n is a non-negative integer

Examples:

   100   XEQ "APNB"  >>>>  2.824655679  X<>Y  149  whence   A100 = 2.824655679 E149              ---Execution time = 88s---

-Likewise:

     A329 = 1.990511241 E499                      ---Execution time = 4m49s---
     A330 = 6.731192145 E500                      ---Execution time = 4m50s---
    A1000 = 8.811881599 E1525                    ---Execution time = 14m33s---
 

     b) 2 M-Code Routines
 

-This routine takes a non-negative integer n  in X-register and returns An in X-register.
-There is no check for alpha data.
-The alpha "register" is cleared
 

082  "B"
00E  "N
010  "P"
001  "A"
2A0  SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
04E  C=0 ALL
1E8  O=C
268  Q=C
168  M=C
35C C=
050  1
1A8  N=C
228  P=C
178  C=M                                                    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
1BE  A=A-1 MS
168  M=C
001  C=
060  AB+1
138  C=L
025  C=
060  AB+C
2EE  ?C#0 ALL
153  JNC+42d
070  N=C ALL
138  C=L
10E  A=C ALL
178  C=M
01D  C=
060  A+C
0B0  C=N ALL
13D  C=
060  AB*C
178  C=M
269  C=
060  AB/C
178  C=M
269  C=
060  AB/C
0AE  A<>C ALL
10E  A=C ALL
158  M=C ALL
0CE  C=B ALL
149  C=
060  AB*CM
1F8  C=O
158  M=C ALL
1B8  C=N
149  C=
060  AB*CM
0AE A<>C ALL
10E  A=C ALL
1E8  O=C
0CE  C=B ALL
1A8  N=C
278  C=Q
158  M=C ALL
238  C=P
031  C=
060  AB+CM
0AE  A<>C ALL
268  Q=C
0CE  C=B ALL
228  P=C
23B  JNC-57d                            goto Loop
1A0  A=B=C=0
278  C=Q
158  M=C ALL
238  C=P
031  C=
060  AB+CM
0E8  X=C
345  ?NCGO
042  CLA

( 83 words )
 
 

      STACK        INPUTS      OUTPUTS
           X             n             An
           L             /             n

  where n < 330  is a non-negative integer

Examples:

    67    XEQ "APNB"  >>>>   A67  = 1.529563780  E99                      ---Execution time = 19s---
   100   XEQ "APNB"  >>>>   A100 = 2.824655678 E149                    ---Execution time = 29s---
   329   XEQ "APNB"  >>>>   A329 = 1.990511240 E499                    ---Execution time = 98s---

Notes:

-In the last 2 examples, the exponent is displayed in an ambiguous form ( E-51 & E-01 respectively )
-Press  ENTER^  LOG   INT  to get the correct value.

-For n > 329, the result may be wrong.
-There will be an infinte loop if n is negative or fractional, press any key to stop the loop.

-The following version of "APNB"  returns the mantissa in X-register and the exponent in Y-register.
 

082  "B"
00E  "N
010  "P"
001  "A"
2A0  SETDEC
3C8  CLRKEY
0F8   C=X
128   L=C
04E  C=0 ALL
0A8  Y=C
168   M=C
1E8  O=C
228  P=C
35C  C=
050   1
0E8   X=C
1A8  N=C
178  C=M                                     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
1BE  A=A-1 MS
168  M=C
001  C=
060  AB+1
138  C=L
025  C=
060  AB+C
2EE  ?C#0 ALL
1DB  JNC+59d
070  N=C ALL
138  C=L
10E  A=C ALL
178  C=M
01D  C=
060  A+C
0B0  C=N ALL
13D  C=
060  AB*C
178  C=M
269  C=
060  AB/C
178  C=M
269  C=
060  AB/C
0AE  A<>C ALL
10E  A=C ALL
158  M=C ALL
0CE  C=B ALL
149  C=
060  AB*CM
1F8  C=O
158  M=C ALL
1B8  C=N
149  C=
060  AB*CM
0AE A<>C ALL
10E  A=C ALL
1E8  O=C
0CE  C=B ALL
1A8  N=C
0B8  C=Y
158  M=C ALL
0F8  C=X
031  C=
060  AB+CM
0CE  C=B
0E8   X=C
0AE  A<>C ALL
0A8  Y=C
21C  PT=2
2E2  ?C#0@PT
22B  JNC-59d                                     goto Loop
10E  A=C ALL
04E  C=0 ALL
130  LDI S&X
100  100
1C6  A=A-C S&X
0AE  A<>C ALL
0A8  Y=C
238   C=P
20E   C=A+C ALL
228   P=C
1F8  C=O
0AE  A<>C ALL
246  C=A-C S&X
1E8  O=C
38B  JNC-15d                                     goto Loop
1A0  A=B=C=0
0B8  C=Y
158  M=C ALL
0F8  C=X
031  C=
060  AB+CM
046  C=0 S&X
0E8  X=C
01A  A=0 M
238   C=P
14E  A=A+C ALL
34E  ?A#0 ALL
3A0  ?NC RTN
130   LDI S&X
013  013
3EE  LSHFA ALL
266  C=C-1 S&X
35E  ?A#0 MS
3EB  JNC-03
38E   RSHFA ALL
0BA  A<>C M
0A8  Y=C
345  ?NCGO
042  CLA

( 116 words )
 
 

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

  where n is a non-negative integer

Examples:

   100   XEQ "APNB"  >>>>  2.824655678  X<>Y  149  whence   A100 = 2.824655678 E149              ---Execution time = 29s---

-Likewise:

     A329 = 1.990511240 E499                      ---Execution time = 98s---
     A330 = 6.731192131 E500                      ---Execution time = 99s---
    A1000 = 8.811881564 E1525                    ---Execution time = 305s---

Notes:

-Registers Z & T are saved but Y-register is lost.
-The alpha register is cleared.
-There will be an infinte loop if n is negative or fractional, press any key to stop the loop.
 

References:

[1]  http://mathworld.wolfram.com/AperysConstant.html
[2]  http://mathworld.wolfram.com/AperyNumber.html