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:

  http://mathworld.wolfram.com/AperysConstant.html
  http://mathworld.wolfram.com/AperyNumber.html