Last added floating point term of harmonic series 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?

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

