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 An
= 4F3 ( -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