sqrt (1.0 - pow (1.0.2)) renvoie -nan

Question

J'ai trouvé un problème de virgule flottante intéressant. Je dois calculer plusieurs racines carrées dans mon code, et l'expression est comme ceci:

sqrt(1.0 - pow(pos,2)) 

où pos va -1,0 à 1,0 dans une boucle. Le -1.0 est bon pour pow, mais quand pos = 1.0, j'obtiens un -nan. Faire des tests, en utilisant gcc 4.4.5 et 12.0 cpi, la sortie de

1.0 - pow(pos,2) = -1.33226763e-15 

et

1.0 - pow(1.0,2) = 0 

ou

poss = 1.0 
1.0 - pow(poss,2) = 0 

Si bien le premier va donner des problèmes, être négatif. Quelqu'un sait pourquoi pow renvoie un nombre inférieur à 0? Le plein code incriminé est ci-dessous:

int main() { 
    double n_max = 10; 
    double a = -1.0; 
    double b = 1.0; 
    int divisions = int(5 * n_max); 
    assert (!(b == a)); 

    double interval = b - a; 
    double delta_theta = interval/divisions; 
    double delta_thetaover2 = delta_theta/2.0; 
    double pos = a; 
    //for (int i = 0; i < divisions - 1; i++) { 
    for (int i = 0; i < divisions+1; i++) { 

    cout<<sqrt(1.0 - pow(pos, 2)) <<setw(20)<<pos<<endl; 

    if(isnan(sqrt(1.0 - pow(pos, 2)))){ 
     cout<<"Danger Will Robinson!"<<endl; 
     cout<< sqrt(1.0 - pow(pos,2))<<endl; 
     cout<<"pos "<<setprecision(9)<<pos<<endl; 
     cout<<"pow(pos,2) "<<setprecision(9)<<pow(pos, 2)<<endl; 
     cout<<"delta_theta "<<delta_theta<<endl; 
     cout<<"1 - pow "<< 1.0 - pow(pos,2)<<endl; 
     double poss = 1.0; 
     cout<<"1- poss "<<1.0 - pow(poss,2)<<endl; 


    } 

    pos += delta_theta; 

} 

return 0; 
} 

Answer 1

Le problème est que les calculs à virgule flottante ne sont pas exacts, et que 1 - 1^2 peut donner de petits résultats négatifs, ce qui donne un sqrt calcul invalide.

Tenir compte de capsulage votre résultat:

double x = 1. - pow(pos, 2.); 
result = sqrt(x < 0 ? 0 : x); 

ou

result = sqrt(abs(x) < 1e-12 ? 0 : x); 

Answer 2

Lorsque vous maintenez incrémenter pos dans une boucle, des erreurs d'arrondi accumulent et dans votre cas, la valeur finale> 1.0. Au lieu de cela, calculez pos par multiplication sur chaque tour pour obtenir seulement une quantité minimale d'erreur d'arrondi.

Answer 3

setprecision(9) va provoquer l'arrondissement. Utilisez un débogueur pour voir quelle est la valeur réelle. Bref, au moins définir la précision au-delà de la taille possible du type que vous utilisez.

Answer 4

Vous ont presque toujours des erreurs d'arrondi lors du calcul avec double, parce que le type double a seulement 15 chiffres décimaux significatifs (52 bits) et beaucoup de nombres décimaux ne sont pas convertibles en nombres à virgule flottante binaire sans arrondi. La norme IEEE contient beaucoup d'efforts pour garder ces erreurs faibles, mais par principe, elle ne peut pas toujours réussir. Pour une introduction complète, voir this document

Dans votre cas, vous devez calculer pos sur chaque boucle et arrondir à 14 chiffres ou moins. Cela devrait vous donner un 0 propre pour le sqrt.

Vous pouvez calco pos dans la boucle comme

pos = round(a + interval * i/divisions, 14); 

avec tour défini comme

double round(double r, int digits) 
{ 
    double multiplier = pow(digits,10); 
    return floor(r*multiplier + 0.5)/multiplier; 
}