(DM42) Harmonic function based on Bernoulli - 32 digits

[Gjermund Skailand]

Harmonic function
H(n) = SUM[1/i] , i=1 ..n
Direct calculation used for values of n < 40
H(n) = 1 + 1/2 + 1/3 ... 1/n

Asymptotic expansion used for larger values of n >= 40
H(n) = ln(n) + gamma + 1/2n - SUM [B[k]/(2k n^(2k))]
H(n) = ln(n) + gamma + 1/2n - ( B[1]/(2*n^2) + B[2]/(4*n^4) + B[3]/(6*n^6) ... )

A version with 32 digit accuracy using array of 10 pre-computed Bernoulli numbers
Largest error is less than 1E-32 (for values below 1000), typical relative error is about 1E-33

Should run on Free42,
run in nstack mode, alternatively insert LNSTK after line 02

Note:
/ is division
x is multiply

00 { 121-Byte Prgm }@ DM42
01 LBL "H(n)" @ Harmonic function, H(n) = SUM[1/i] , i=1 ..n
02 FUNC 11 @ H(n) = 1 + 1/2 + 1/3 ... 1/n
03 40
04 X>Y @ uses asymptotic expansion with Bernoulli numbers
05 GTO 01 @ small numbers are calculated directly
06 CLX
07 RCL ST Y @ H(n) =
08 LN @ ln(n)
09 5.772156649015328606065120900824024E-1 @ gamma with 34 DIGITS
10 2
11 RCLx ST T
12 1/X
13 R^N 4
14 X^2
15 1
16 INDEX "Bn"
17 LBL 00
18 RCLx ST Y
19 RCLEL
20 RCLIJ
21 CLX
22 2
23 x
24 /
25 RCL/ ST Y
26 STO- ST T
27 DROP
28 I+
29 FC? 77
30 GTO 00
31 DROPN 2
32 +
33 GTO 03
34 LBL 01 @ direct sum for small values
35 CLX
36 LBL 02
37 RCL ST Y
38 1/X
39 +
40 DSE ST Y
41 GTO 02
42 LBL 03
43 +
41 END

Program to expand 10 Bernoulli numbers into an array,

00 { 124-Byte Prgm }
01 Lbl "B10"
02 10
03 1
04 NEWMAT
05 EDIT
06 LBL 00
07 CLX
08 1
09 6
10 /
11 ->
12 1
13 30
14 /
15 +/-
16 ->
17 1
18 42
19 /
20 ->
21 1
22 30
23 /
24 +/-
25 ->
26 5
27 66
28 /
29 ->
30 691
31 2730
32 /
33 +/-
34 ->
35 7
36 6
37 /
38 ->
39 3617
40 510
41 /
42 +/-
43 ->
44 43867
45 798
46 /
47 ->
48 174611
49 330
50 /
51 +/-
52 EXITALL
53 STO "Bn"
54 END


Best regards
Gjermund Skailand
2022-12-12