Prime calculator sum(1.0002^n/n, 1, 20000)

[Albert Chan]

Note that z = 0.0002 is a small number, perfect for taylor series approximation.

\(f = \Large {(1+z)^k \over k} \displaystyle \normalsize
= {1\over k} + \sum_{j=1}^k \binom{k}{j}{z^j \over k}
= {1\over k} + \sum_{j=1}^k \binom{k-1}{j-1}{z^j \over j}
\)

\(\small \displaystyle\sum_{k=j}^n \binom{k-1}{j-1}
= \sum_{k=j}^n \left[\binom{k}{j} - \binom{k-1}{j}\right]
= \binom{n}{j} - \binom{j-1}{j}
= \binom{n}{j} \)       // middle terms telescopingly cancelled

Flip the order of sum, \(\large \sum _{k=1}^n \sum _{j=1}^k
≡ \sum _{j=1}^n \sum _{k=j}^n\)

\(\displaystyle\sum_{k=1}^n f
= H(n) + \sum_{j=1}^n \displaystyle\sum_{k=j}^n \binom{k-1}{j-1} {z^j \over j}
= H(n) + \sum_{j=1}^n \binom{n}{j}{z^j\over j} \)

H(n), we use Harmonic Number asymptotic formula:

lua> n, z = 20000, 0.0002
lua> euler_gamma = 0.5772156649015329
lua> Hn = (1-1/(6*n))/(2*n) + euler_gamma + log(n)
lua> Hn
10.480728217229327

Sum the other terms, until converged.
With small z, we may not need n terms for convergence.

lua> j, t, s = 1, 1, -1
lua> repeat s=s+t/j; j=j+1; t=t*(n-j+1)/j*z; until s==s+t/j
lua> s = n*z * (s+1)
lua> Hn + s       -- = sum(1.0002^k/k, k = 1 .. 20000)
28.14397383505517

lua> j                -- much less than 20000 terms
30

[Joe Horn]

(12-28-2020 04:21 PM)Didier Lachieze Wrote:  

(12-28-2020 03:02 PM)Joe Horn Wrote:  Anybody know how to tell the DM42 not to update the stack display while running? Thanks in advance!

You can try with RefLCD set to 0.

Awesome, thanks! Doing that makes the DM42 run the above program in 40 seconds (powered by USB).

[Gil]

Just note that, due to roundings, the last digits of the
given result 28.14397383505516924725756564071173
are — not surprisingly — incorrect after checking with Wolfram Alpha
(458 [55], instead of 173).

Gil

[Albert Chan]

(12-29-2020 01:57 AM)Gil Wrote:  Just note that, due to roundings, the last digits of the
given result 28.14397383505516924725756564071173
are — not surprisingly — incorrect after checking with Wolfram Alpha
(458 [55], instead of 173).

The errors mostly comes from implementation of Y^X.
For integer X, Free42 were patched to repeated squaring: Y^X = (Y^(X//2))^2 * Y^(X%2)

Example, for 7 squarings:
1.0002 128 Y^X               → 1.025927868161809572716800732564789
1.0002 LN 128 × E^X       → 1.025927868161809572716800732564790

Extending to 14 squarings:
1.0002 16384 Y^X           → 26.48218820206050462500040752260032
1.0002 LN 16384 × E^X   → 26.48218820206050462500040752260470 (last digit should be 1)

Because all terms are using the same base, almost all errors fall on the same side.
And, it get worse and worse, with 1 ULP error for 7 squarings to 471-32 = 439 ULP for 14.

I replaced Y^X by E^(X*LN(Y)), and the sum ends in 447, much closer to true result.

[Claudio L.]

(12-27-2020 07:03 PM)John Keith Wrote:  

(12-27-2020 03:04 AM)Gil Wrote:  On the HP50G, I get the wished answer by writing
directly sum(1.0002^n/n, 1, 20000).

For the HP 50, the correct syntax is '\GS(n=1, 20000, 1.0002^n/n)' where \GS is the Sigma character (right-shift S).

On a physical HP 50g, I get 28.1439738407 in 354.6 seconds, or just under 6 min.

Oh, are we doing benchmarks now? Rewriting the expression in RPL as a FOR NEXT loop, a physical 50g with newRPL gives the answer in 6.4 seconds for 32 digits and 5.0 seconds with 16 digits precision.
As a curiosity, I ran the same on a physical Prime G1 with newRPL and it took 11.9 seconds with 32 digits which makes no sense... unless the clock is stuck at 100 MHz, will need to investigate this and report again, sorry.

[Gil]

Differences of times might be substantial — on the same "device/model" —
according to how you calculate:
- with integer limits or with no-integer (i. e. limits with dot);
- with build-in summation (Greek S) or the loop For... Next.

With Emu48 Android Application
on my phone Samsung galaxy S6,
I got the following results:

7.7" with points after 0., 1. and 2000.:
\<< 0. 1. 20000.
FOR i 1.0002 i ^ i / +
NEXT
\>>

22.8" without points after 0, 1 and 20000:
\<< 0 1 20000
FOR i 1.0002 i ^ i / +
NEXT
\>>

9.1" without points, i. e. limits like integers,
but mandatory with -105 SF!
\<< '\GS(n=1,20000,1.0002^n/n)' -105. SF EVAL
\>>

8.9" with points already included after 1. and 20000.:
\<< '\GS(n=1.,20000.,1.0002^n/n)' EVAL
\>>

Conclusion:
The self-made program is the fastest... providing you do not
forget to put dots to the limits (no-integers).

Regards, [/align]
Gil