Réduire l'erreur d'arrondi de la somme des carrés des flotteurs

Question

J'essaie de calculer la somme des carrés d'un tableau de flottants.

Comment faire pour réduire l'erreur d'arrondi? J'essaie d'additionner environ 5.000.000 flottants dans la boucle interne de mon programme actuel.

test.cpp:

#include <iostream> 
#include <stdint.h> 
template <typename Sum, typename Element> 
Sum sum(const size_t st, const size_t en) { 
    Sum s = 0; 
    for (size_t i = st; i < en; ++ i) { 
     s += Element(i)*Element(i); 
    } 
    return s; 
} 
int main() { 
    size_t size = 100000; 
    std::cout << "double, float: " 
       << sum<double, float>(0,size) << "\n"; 
    std::cout << "int, int: " 
       << sum<int, int>(0,size) << "\n"; 
} 

Sortie:

double, float: 3.33328e+14 
int, int: 216474736 

Answer 1

Si le format des flotteurs est connu, comme IEEE, puis un tableau indexé par l'exposant d'un flotteur peut être utilisé pour stocker des sommes partielles, puis additionnées pour produire la somme totale. Pendant la mise à jour du tableau, seuls les flottants ayant le même exposant sont additionnés et stockés dans le tableau à l'emplacement approprié. La somme finale va du plus petit au plus grand. Pour C++, le tableau et les fonctions peuvent être membres d'une classe.

Exemple pour les flotteurs où le tableau est passé en tant que paramètre aux fonctions:

/* clear array */ 
void clearsum(float asum[256]) 
{ 
size_t i; 
    for(i = 0; i < 256; i++) 
     asum[i] = 0.f; 
} 

/* add a number into array */ 
void addtosum(float f, float asum[256]) 
{ 
size_t i; 
    while(1){ 
     /* i = exponent of f */ 
     i = ((size_t)((*(unsigned int *)&f)>>23))&0xff; 
     if(i == 0xff){   /* max exponent, could be overflow */ 
      asum[i] += f; 
      return; 
     } 
     if(asum[i] == 0.f){  /* if empty slot store f */ 
      asum[i] = f; 
      return; 
     } 
     f += asum[i];   /* else add slot to f, clear slot */ 
     asum[i] = 0.f;   /* and continue until empty slot */ 
    } 
} 

/* return sum from array */ 
float returnsum(float asum[256]) 
{ 
float sum = 0.f; 
size_t i; 
    for(i = 0; i < 256; i++) 
     sum += asum[i]; 
    return sum; 
} 

Answer 2

Lorsque vous utilisez le type int pour Element, il y a un trop-plein sur la place pour chaque i après std::sqrt(std::numeric_limits<int>::max()), ce qui pourrait être 46341 sur votre système. Il y a également un débordement sur la somme lorsqu'elle atteint std::numeric_limits<int>::max().

Vous pouvez utiliser le type long ou long long au lieu de int pour augmenter ce nombre.

Il est aussi une bonne idée de stocker ou de jeter la première float dans un double ou long double avant de multiplier pour réduire l'erreur sur la virgule flottante opération carré ainsi. L'arrondi sur les dernières étapes d'un ensemble de calculs donne toujours un meilleur résultat que l'arrondi sur les premières étapes, car vous évitez de propager (et d'augmenter) les erreurs de représentation sur vos calculs internes.

Si vous voulez vraiment la précision, et ne veulent pas de réinventer la roue en utilisant des techniques compliquées, vous pouvez utiliser une bibliothèque multi-précision comme GNU Multi-Precision Library ou Boost Multiprecision: https://en.wikipedia.org/wiki/List_of_arbitrary-precision_arithmetic_software

Ils sont plus précis que le type long double de votre système

Answer 3

Si tout ce que vous voulez faire est d'ajouter des carrés de valeurs consécutives, utilisez la formule n*(n+1)*(2n+1)/6 pour calculer la somme des carrés de toutes les valeurs 1-n.

Cela élimine la plupart des effets d'arrondi tant que vous utilisez un type qui peut représenter le résultat. Par exemple;

template<typename Sum> Sum sumsq(size_t n) 
{ 
    // calculates sum of squares from 1 to x 
    // assumes size_t can be promoted to a Sum 

    Sum temp(n);  // force promotion to type Sum 
    return temp * (temp + 1)* (2*temp + 1)/6; 
} 

template<typename Sum> Sum alternate_sum(size_t st, size_t en) 
{ 
     Sum result = sumsq(en - 1); 
     if (st > 0) result -= sumsq(st-1); 
     return result; 
} 

int main() 
{ 
    size_t size = 100000; 
    std::cout << "double, float: " 
       << alternate_sum<double>(0,size) << "\n"; 
    std::cout << "int, int: " 
      << alternate_sum<long long>(0,size) << "\n"; 
} 

Notez que, pour size égal à 100000, en utilisant int pour maintenir le résultat donne un comportement non défini (dépassement d'un type intégral signé).

Les -1 s du alternate_sum() reflètent vos boucles étant de la forme for (size_t i = st; i < en; ++ i)

Vous pouvez éliminer l'utilisation du type size_t comme une caractéristique fixe, mais je vais laisser cela comme un exercice. BTW: puisque vous dites que ce code est dans une boucle interne, il est à noter que cette formule sera significativement plus rapide que les boucles que vous avez utilisées.

Answer 4

Un flotteur a 24 bits significatifs, tandis qu'un double en a 53. Donc, vous avez 29 bits de garde, ce qui est environ 100 fois plus que 5000000.

Donc, les erreurs d'arrondi ne se poseront que pour les valeurs qui sont 100 fois plus petites que les plus grandes.

Notez également que dans l'architecture Intel, les registres à virgule flottante contiennent en réalité des nombres de précision étendus sur 80 bits, dont 63 sont significatifs.

Ensuite, seuls les nombres inférieurs à 100 000 fois les plus grands entraîneront une troncature.

Si vous vous inquiétez vraiment?