(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

[Werner]

Hi Gjermund!
Doing the same with a Horner-like scheme, independent of stack mode (I think).
I also stored gamma in a variable..

00 { 75-Byte Prgm }
01▸LBL "Hn"
02 FUNC 11
03 40
04 X>Y?
05 GTO 00
06 R↓
07 X^2
08 0
09 INDEX "B10"
10 I-
11▸LBL 01
12 RCLIJ
13 STO+ ST X
14 ×
15 RCLEL
16 X<>Y
17 ÷
18 +
19 RCL÷ ST Y @ in 4STK mode, R^ / would work as well
20 I-
21 FC? 77
22 GTO 01
23 X<>Y
24 SQRT @ n
25 LN
26 LASTX
27 STO+ ST X
28 1/X
29 +
30 RCL+ "GAMMA"
31 X<>Y
32 -
33 RTN
34▸LBL 00
35 CLX
36▸LBL 02
37 RCL ST Y
38 1/X
39 +
40 DSE ST Y
41 GTO 02
42 END

Cheers, Werner

[Gjermund Skailand]

Hi Werner,
as always, a shorter and faster version from you.
I did test the accuracy by expanding all part sums on the stack ( I love the large stack ), and summing in reverse order, from smallest to largest. But apart from the first three items I found that it did not matter.
DM42 is my new favourite calculator. Wold anyone happen to know if it is possible to permanently store variables in any other way than storing the entire state?
br Gjermund

[Werner]

Hi Gjermund.
Sadly, copy/paste works only for programs, and whatever is in the X-register.
You can't export variables. Would be nice to have though.
Cheers, Werner