Último término añadido en coma flotante de la serie armónica

Question

En un sistema que utiliza números de punto flotante de precisión simple IEEE 754, si empezamos a calcular la suma $\sum_{i=0}^n 1/i$ Por la precisión, la suma no llegará al infinito, pero después de un término N, cualquier término que se añada no cambiará la suma, que en nuestro caso debe ser ~15,4 . ¿Cuál es la forma de encontrar este término N?

Answer 1

Tenemos $\sum_{i=0}^n 1/i\approx \log n + 0.577$ que es fácilmente lo suficientemente cerca para lo que queremos. La mantisa es $24$ bits, por lo que se necesita $\frac 1N \lt 2^{-24+E}$ , donde $E$ es el exponente del valor actual de la suma. $E(n)=\lfloor \log_2 \log (n + 0.577)\rfloor$ Esto da $\frac 1n \lt 2^{-24}(\log n+0.577)$ Dejo que Alpha hacer los números, obteniendo $n=1154191$ Si necesitas la respuesta exacta, debes preocuparte por el redondeo en $\frac 1N$ cuando se añade a la suma. El LSB de la suma será $2^{-20}$ , por lo que si sólo necesitamos $\frac 1N$ menos que esto, necesitamos $N=2^{20}+1$

Answer 2

Una estimación aproximada: $\epsilon$ se trata de $5\times10^{-8}$ . Si $\oplus$ es la adición aritmética de la máquina tenemos $$1 \oplus \epsilon= 1$$ y para la arbitrariedad $n$ $$ n \oplus n\delta = n$$ para todos $$\delta \lt \delta_x$$ con $\delta_x \in (\epsilon, 2\epsilon)$ . Para la aritmética de 32 bits tenemos $$\epsilon=2^{−24} ≈ 5.96\cdot 10 {-8}$$ Tenemos aproximadamente $$\sum_{i=1}^{n}{\frac{1}{i}}\approx \ln(n)+0.6$$

Así que si $$5 \cdot 10^{-8}\cdot\bigoplus_{i=1}^{n}{\frac{1}{i}} \approx\frac{1}{n}$$ la suma no aumentará más. $\bigoplus_{i=1}^{n}{\frac{1}{i}}$ es la suma $\sum_{i=1}^{n}{\frac{1}{i}}$ calculado con la aritmética de la máquina. Pero si $\bigoplus_{i=1}^{n}{\frac{1}{i}}$ aproximadamente $\sum_{i=1}^{n}{\frac{1}{i}}$ podemos sustituirlo por $\ln(n)+0.6$ y obtener $$5 \cdot 10^{-8}\cdot \ln(n)\approx\frac{1}{n}$$ y además $$n\ln{(n)} \approx 10^7$$ La solución de esta ecuación es $$n\approx10^{6}$$ y $$\sum_{i=1}^{n}{\frac{1}{i}} \approx \ln(n)+0.6 \approx14.4$$

El siguiente programa de Python 3.5 estimará eps, estimará $n$ y tratar de encontrar $n$ mediante el cálculo de la serie armónica. Tenga en cuenta que algunos de los cálculos pueden no estar en float, sino que pueden utilizar una precisión mayor. Utilizo arrays de python para simular la aritmética de 32 bits según 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    

Esto imprime lo siguiente:

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

Los solteros del IEEE tienen $1$ bit por signo, $8$ bits para el exponente y $24$ bits para la precisión del significante ( $23$ bits almacenados explícitamente).

En el IEEE, cuando se lleva a cabo algún cálculo, se supone que se calcula el resultado con una precisión mayor internamente y luego se redondea el número al bit representable más cercano. Si un número es mayor que la mitad del LSB de la precisión del significante, el resultado se redondeará hacia arriba, y si es menor que la mitad del LSB, se redondeará hacia abajo. Si el número es exactamente la mitad del LSB, por defecto se redondea a par. Para nuestros propósitos, no sabemos si se redondeará hacia arriba o hacia abajo para este caso particular.

Cuando $N \approx 2^{21}$ , $H_n \approx 15.1333$ . Esto quita $3$ bits de la precisión significativa. El LSB corresponde a $2^{-23+3} = 2^{-20}$ . Cuando $N = 2^{21}$ o $N = 2^{21}+1$ El $\frac{1}{N}$ será menor que la mitad del LSB y la suma termina.

Esto coincide con lo que mencionó @lhf en el comentario.

Por un argumento similar, los dobles del IEEE tienen $1$ bit por signo, $11$ para el exponente y $53$ bits para la precisión del significante ( $52$ bits almacenados explícitamente).

En $N \approx 2^{48}$ , $H_n \approx 33.848$ . Esto quita $5$ bits de la precisión del significante. El LSB corresponde a $2^{-52+5} = 2^{-47}$ . En cualquiera de los dos casos $N = 2^{48}$ o $N = 2^{48} + 1$ , $\frac{1}{N}$ será menor que la mitad del LSB y la suma termina. En mi ordenador, la suma deja de actualizarse en $N = 2^{48} + 1$ .