Cuando se hace aritmética en coma flotante, se debe considerar que cada adición induce un error relativo de aproximadamente el float $\epsilon$ definido como el menor número positivo tal que $1+\epsilon == 1$ se evalúa como verdadero. Este número puede ser consultado mediante std::numeric_limits<float>::epsilon() en C++, o FLT_EPSILON en C.
Así que para generar una secuencia cuyo valor dependa del orden de la suma, encuentra una secuencia que sepas que es $\approx 1$ y cuyos términos posteriores son menores que $\epsilon$ pero que, sin embargo, puede acumularse hasta algo $>\epsilon$ . Aquí hay un código que hace precisamente eso:
#include <limits>
#include <iostream>
#include <cmath>
int main()
{
float zeta_2_ls = 0;
//Make n strictly > 1/sqrt(eps), or else the two sums will be identical:
unsigned long long n = std::sqrt(100000000/std::numeric_limits<float>::epsilon());
for(unsigned long long ii = 1; ii <= n; ++ii)
{
zeta_2_ls += 1.0f/(ii*ii);
}
float zeta_2_sl = 0;
for(unsigned long long ii = n; ii > 0; --ii)
{
zeta_2_sl += 1.0f/(ii*ii);
}
std::cout << std::setprecision(12);
std::cout << std::fixed;
std::cout << "Zeta(2), summed largest to smallest: " << zeta_2_ls << std::endl;
std::cout << "Zeta(2), summed smallest to largest: " << zeta_2_sl << std::endl;
std::cout << "Zeta(2), using pi^2/6 to float prec: " << M_PI*M_PI/6.0f << std::endl;
}
La salida:
Zeta(2), summed largest to smallest: 1.644725322723
Zeta(2), summed smallest to largest: 1.644933938980
Zeta(2), using pi^2/6 to float prec: 1.644934066848
Se puede calcular que el Distancia ULP de la respuesta exacta es sólo 1 para la serie sumada de pequeño a grande, pero 1751 para la de grande a pequeño.
He aquí otro ejemplo, en el que los términos llegan a cero, la serie converge, pero el valor de la suma depende del número de términos
#include <stdio.h>
#include <math.h>
float saw_tooth(float x, unsigned int terms)
{
float y = 0;
for(unsigned int ii = 1; ii <= terms; ++ii)
{
if(ii & 1)
{
y += sinf(ii*x)/ii;
}
else
{
y -= sinf(ii*x)/ii;
}
}
y *= 2.0/M_PI;
return y;
}
int main()
{
float x = M_PI;
unsigned int terms = 10;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 100000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 1000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 10000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 100000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
x = M_PI/2.0;
terms = 10;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 100000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 1000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 10000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 100000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
printf("Why is the so far from zero? sinf(M_PI) = %g\n", sinf(M_PI));
}
La salida:
Saw tooth(3.14159) with 10 terms = -5.30531e-07
Saw tooth(3.14159) with 100000 terms = -0.0055615
Saw tooth(3.14159) with 1000000 terms = -0.0555674
Saw tooth(3.14159) with 10000000 terms = -0.491704
Saw tooth(3.14159) with 100000000 terms = -0.6344
Saw tooth(1.5708) with 10 terms = 0.531527
Saw tooth(1.5708) with 100000 terms = 0.499998
Saw tooth(1.5708) with 1000000 terms = 0.500003
Saw tooth(1.5708) with 10000000 terms = 0.500004
Saw tooth(1.5708) with 100000000 terms = 0.616448
Why is the so far from zero? sinf(M_PI) = -8.74228e-08
Por supuesto, la reducción de la gama está causando la evaluación de sinf(ii*M_PI) para convertirse en ~ FLT_EPSILON y por lo tanto, los términos son más o menos. . a la deriva. Pero, para ser honesto, no puedo explicar por qué se desvían tanto, dado el efecto mitigador de la división por ii .
