sqrt, carrés parfaits et erreurs en virgule flottante

Dans la fonction sqrt de la plupart des langues (bien que C et Haskell m’intéressent principalement), existe-t-il des garanties quant à la restitution exacte de la racine carrée d’un carré parfait? Par exemple, si je fais sqrt(81.0) == 9.0 , est-ce sûr ou risque-t-il que sqrt retourne 8.999999998 ou 9.00000003?

Si la précision numérique n’est pas garantie, quelle serait la meilleure méthode pour vérifier qu’un nombre est un carré parfait? Prenez la racine carrée, prenez le sol et le plafond et assurez-vous qu’ils correspondent bien au nombre initial?

Je vous remercie!

Dans IEEE 754 à virgule flottante, si la valeur double précision x est le carré d’un nombre représentable non négatif y (c.-à-d. Y * y == x et que le calcul de y * y n’entraîne aucun arrondi, dépassement ou dépassement), alors sqrt (x) retournera y.

Tout cela parce que sqrt doit être correctement arrondi par la norme IEEE 754. En d’autres termes, sqrt (x), pour tout x, sera le double le plus proche de la racine carrée réelle de x. Ce carré qui fonctionne pour les carrés parfaits est un simple corollaire de ce fait.

Si vous voulez vérifier si un nombre à virgule flottante est un carré parfait, voici le code le plus simple auquel je puisse penser:

 int issquare(double d) { if (signbit(d)) return false; feclearexcept(FE_INEXACT); double dd = sqrt(d); asm volatile("" : "+x"(dd)); return !fetestexcept(FE_INEXACT); } 

J’ai besoin du bloc vide asm volatile qui dépend de dd car sinon votre compilateur pourrait être intelligent et “optimiser” le calcul de dd .

J’ai utilisé quelques fonctions étranges de fenv.h , à savoir feclearexcept et fetestexcept . C’est probablement une bonne idée de regarder leurs pages de man .

Une autre stratégie que vous pouvez peut-être utiliser est de calculer la racine carrée, de vérifier si elle a défini des bits dans les 26 bits les plus bas de la mantisse et de vous plaindre si tel est le cas. J’essaie cette approche ci-dessous.

Et j’avais besoin de vérifier si d est égal à zéro, sinon il peut renvoyer true pour -0.0 .

EDIT : Eric Postpischil a suggéré que bidouiller avec la mantisse serait peut-être mieux. Étant donné que issquare ci-dessus ne fonctionne pas dans un autre compilateur populaire, clang , j’ai tendance à être d’accord. Je pense que le code suivant fonctionne:

 int _issquare2(double d) { if (signbit(d)) return 0; int foo; double s = sqrt(d); double a = frexp(s, &foo); frexp(d, &foo); if (foo & 1) { return (a + 33554432.0) - 33554432.0 == a && s*s == d; } else { return (a + 67108864.0) - 67108864.0 == a; } } 

Ajouter et soustraire 67108864.0 de a a pour effet d’effacer les 26 bits inférieurs de la mantisse. Nous aurons a retour exactement quand ces éléments étaient clairs au départ.

Selon cet article , qui traite de la démonstration de la correction de la racine carrée à virgule flottante IEEE:

La norme IEEE-754 pour l’arithmétique binary en virgule flottante [1] requirejs que le résultat d’une opération de division ou de racine carrée soit calculé comme si il était en précision infinie, puis arrondi à l’un des deux nombres en virgule flottante les plus proches de la précision spécifiée. qui entourent le résultat infiniment précis

Puisqu’un carré parfait pouvant être représenté exactement en virgule flottante est un entier et que sa racine carrée est un entier pouvant être représenté avec précision, la racine carrée d’un carré parfait doit toujours être parfaitement correcte.

Bien entendu, rien ne garantit que votre code sera exécuté avec une bibliothèque à virgule flottante IEEE conforme.

@tmyklebu a parfaitement répondu à la question. En complément, voyons une alternative peut-être moins efficace pour tester un carré parfait de fractions sans directive asm.

Supposons que nous ayons un sqrt conforme à la norme IEEE 754 qui arrondit correctement le résultat.
Supposons que les valeurs exceptionnelles (Inf / Nan) et les zéros (+/-) soient déjà gérées.
Décomposons sqrt(x) en I*2^mI est un entier impair.
Et où I couvre n bits: 1+2^(n-1) < = I < 2^n .

Si n > 1+floor(p/2)p est la précision de la virgule flottante (par exemple, p = 53 et n> 27 en double précision)
Alors 2^(2n-2) < I^2 < 2^2n .
Comme I est impair, I^2 est aussi et s'étend donc sur> p bits.
Ainsi, I n'est pas la racine carrée exacte de tout point flottant pouvant être représenté avec cette précision.

Mais étant donné I^2<2^p , pourrions-nous dire que x était un carré parfait?
La réponse est évidemment non. Une expansion de taylor donnerait

 sqrt(I^2+e)=I*(1+e/2I - e^2/4I^2 + O(e^3/I^3)) 

Ainsi, pour e=ulp(I^2) jusqu’à sqrt(ulp(I^2)) la racine carrée est correctement arrondie à rsqrt(I^2+e)=I ... (arrondi au plus proche pair ou tronqué ou mode étage).

Nous devrions donc affirmer que sqrt(x)*sqrt(x) == x .
Mais test ci - dessus ne suffit pas, par exemple, en supposant IEEE 754 double précision, sqrt(1.0e200)*sqrt(1.0e200)=1.0e200 , où 1.0e200 est exactement 99999999999999996973312221251036165947450327545502362648241750950346848435554075534196338404706251868027512415973882408182135734368278484639385041047239877871023591066789981811181813306167128854888448 dont le premier facteur premier est 2^613 , à peine un carré parfait de n'importe quelle fraction ...

Nous pouvons donc combiner les deux tests:

 #include  bool is_perfect_square(double x) { return sqrt(x)*sqrt(x) == x && squared_significand_fits_in_precision(sqrt(x)); } bool squared_significand_fits_in_precision(double x) { double scaled=scalb( x , DBL_MANT_DIG/2-ilogb(x)); return scaled == floor(scaled) && (scalb(scaled,-1)==floor(scalb(scaled,-1)) /* scaled is even */ || scaled < scalb( sqrt((double) FLT_RADIX) , DBL_MANT_DIG/2 + 1)); } 

EDIT: Si nous voulons nous limiter au cas des entiers, nous pouvons également vérifier cet floor(sqrt(x))==sqrt(x) ou utiliser des hacks de bits sales dans squared_significand_fits_in_precision ...

Au lieu de sqrt(81.0) == 9.0 , essayez 9.0*9.0 == 81.0 . Cela fonctionnera toujours tant que le carré est dans les limites de la magnitude en virgule flottante.

Edit: Je ne savais probablement pas bien ce que je voulais dire par “magnitude à virgule flottante”. Ce que je veux dire est de garder le nombre dans la plage des valeurs entières pouvant être conservées sans perte de précision, inférieur à 2 ** 53 pour un double IEEE. Je m’attendais également à ce qu’il y ait une opération distincte pour s’assurer que la racine carrée était un entier.

 double root = floor(sqrt(x) + 0.5); /* rounded result to nearest integer */ if (root*root == x && x < 9007199254740992.0) /* it's a perfect square */