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