Last added floating point term of harmonic series

Question

In a system using IEEE 754 single precision floating point numbers, if we start calculating the sum $\sum_{i=0}^n 1/i$ , because of the precision the sum will not go to infinity but after a term N, any term added won't change the sum, which in our case must be ~15.4 . What's the way of finding this term N?

Answer 1

We have $\sum_{i=0}^n 1/i\approx \log n + 0.577$, which is easily close enough for what we want. The mantissa is $24$ bits, so you need $\frac 1N \lt 2^{-24+E}$, where $E$ is the exponent for the current value of the sum. $E(n)=\lfloor \log_2 \log (n + 0.577)\rfloor$ This gives $\frac 1n \lt 2^{-24}(\log n+0.577)$ I let Alpha do the numerics, getting $n=1154191$ If you need the exact answer, you should worry about the rounding on $\frac 1N$ when it is added into the sum. The LSB of the sum will be $2^{-20}$, so if we just need $\frac 1N$ less than this, we need $N=2^{20}+1$

Answer 2

A rough estimation: $\epsilon$ is about $5\times10^{-8}$. If $\oplus$ is machine arithmetik addition we have $$1 \oplus \epsilon= 1$$ and for arbitrary $n$ $$ n \oplus n\delta = n$$ for all $$\delta \lt \delta_x$$ with $\delta_x \in (\epsilon, 2\epsilon)$. For 32 bit arithmetic we have $$\epsilon=2^{−24} ≈ 5.96\cdot 10 {-8}$$ We have approximately $$\sum_{i=1}^{n}{\frac{1}{i}}\approx \ln(n)+0.6$$

So if $$5 \cdot 10^{-8}\cdot\bigoplus_{i=1}^{n}{\frac{1}{i}} \approx\frac{1}{n}$$ the sum won't increase anymore. $\bigoplus_{i=1}^{n}{\frac{1}{i}}$ is the sum $\sum_{i=1}^{n}{\frac{1}{i}}$ calculated using machine arithmetic. But if $\bigoplus_{i=1}^{n}{\frac{1}{i}}$ approximately $\sum_{i=1}^{n}{\frac{1}{i}}$ we can substitue it by $\ln(n)+0.6$ and get $$5 \cdot 10^{-8}\cdot \ln(n)\approx\frac{1}{n}$$ and further $$n\ln{(n)} \approx 10^7$$ The solution of this equation is $$n\approx10^{6}$$ and $$\sum_{i=1}^{n}{\frac{1}{i}} \approx \ln(n)+0.6 \approx14.4$$

The folowing Python 3.5 program will estimate eps, estimate $n$ and try to find $n$ by calculating the harmonic series. Note that some of the calculation may be not in float but may use a higher precision. I use python arrays to simulate 32-bit arithmetic according to https://stackoverflow.com/a/2232650/754550 .

import array
a=array.array('f')
a.append(1.0)
a.append(1.0)
a.append(1.0)
print('bytes of single float:',a.itemsize)
print('estimate eps:')
n=0
while True:
    n+=1
    a[1]=a[0]
    a[1]+=1.0/2**n
    if (a[0]==a[1]):
        print('eps ~ 2 **',-n)
        break

print('')
estimated_sum=14.4
print('find n for estimated sum',estimated_sum)
eps=1.0/2**24
a[0]=estimated_sum
for i in range(2):
    a[1]=a[0]
    delta=a[0]*(eps/2**i)
    a[1]+=delta
    print('n =',int(1/delta),'estimated_sum==estimated_sum+1/n (',a[0]==a[1],')')

print('') 
print('harmonic series:')
print('calculate n such that h(n)=h(n)+1/n')
n=0
a[1]=0
while True:
    n+=1
    a[2]=1.0/n
    # remebemr the (n-1)th partial sum of the harmonic series
    # calculate the n-th partial sum of the  harmonic series
    # terminate if the partial sum does not change anymore
    a[0]=a[1]
    a[1]+=a[2]
    if (a[0]==a[1]): 
        print('n =',n)
        print('h(n) = ',a[0])
        break    

This prints out the following:

bytes of single float: 4
estimate eps:
eps ~ 2 ** -24

find n for estimated sum 14.4
n = 1165084 estimated_sum==estimated_sum+1/n ( False )
n = 2330168 estimated_sum==estimated_sum+1/n ( True )

harmonic series:
calculate n such that h(n)=h(n)+1/n
n = 2097152
h(n) =  15.403682708740234

Answer 3

The IEEE singles has $1$ bit for sign, $8$ bits for exponent and $24$ bits for significand precision ($23$ bits explicitly stored).

In IEEE, when you carry out some computation, you are suppose to compute the result in higher precision internally and then round the number to nearest representable bit. Adding a number greater than half the LSB of the significand precision will cause the result rounded up, less than half the LSB will be rounded down. If the number is exactly half the LSB, the default is round to even. For our purposes, we don't know whether it will round up or down for this particular case.

When $N \approx 2^{21}$, $H_n \approx 15.1333$. This takes away $3$ bits from the significant precision. The LSB corresponds to $2^{-23+3} = 2^{-20}$. When $N = 2^{21}$ or $N = 2^{21}+1$, the $\frac{1}{N}$ will be small than half of the LSB and the sum terminates.

This matches what @lhf mentioned in comment.

By a similar argument, the IEEE doubles has $1$ bit for sign, $11$ bit for exponent and $53$ bits for significand precision ($52$ bits explicitly stored).

At $N \approx 2^{48}$, $H_n \approx 33.848$. This takes away $5$ bits from the significand precision. The LSB corresponds to $2^{-52+5} = 2^{-47}$. At either $N = 2^{48}$ or $N = 2^{48} + 1$, $\frac{1}{N}$ will be smaller than half of the LSB and the sum terminates. On my computer, the sum stop updating at $N = 2^{48} + 1$.