Racine cubique et approximation

victor.debouille@gmail.com a écrit

Bonjour, je suis plus ou moins débutant en C, et je rencontre un
problème assez génant avec les racines cubiques.

Que ce soit avec la fonction pow de math.h :
pow(x3, 1./3)
Ou avec la fonction cbrt de cette même librairie
cbrt(x3)

Si mon x3 vaut 64, la racine cubique va valloir 3.99999999999 (et des
poussières) et non 4. Ce qui est assez embetant pour vérifier si ma
racine cubique est entière ou non.

Auriez vous une idée pour vérifier que ma racine cubique est entière?

Merci d'avance !

Tu fais l'approximation
Tu arrondis à l'entier le plus proche
Tu portes au cube
Tu compares a la veleur initiale

Si les deux sont égale, c'est tout bon

victor.debouille@gmail.com writes:

Bonjour, je suis plus ou moins débutant en C, et je rencontre un
problème assez génant avec les racines cubiques.

Que ce soit avec la fonction pow de math.h :
pow(x3, 1./3)
Ou avec la fonction cbrt de cette même librairie
cbrt(x3)

Si mon x3 vaut 64, la racine cubique va valloir 3.99999999999 (et des
poussières) et non 4. Ce qui est assez embetant pour vérifier si ma
racine cubique est entière ou non.

Auriez vous une idée pour vérifier que ma racine cubique est entièr e?

Les calculs en virgule flottante sont des calculs approchés. Pour
tester l'égalité de a et de b on teste en réalité si la différenc e est
petite, au moyen du test

fpclassifiy(a - b) == FP_ZERO

Exemple de programme:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

int main(void)
{
double a = 64.;
double b = 4.;
double c = 4.;

b-= pow(a, 1./3);
c-= cbrt(a);

printf("b = %.64fn", b);
if(fpclassify(b) == FP_ZERO) {
printf("Classificationn");
}
if(b == 0.) {
printf("Égalitén");
}

printf("c = %.64fn", c);
if(fpclassify(c) == FP_ZERO) {
printf("Classificationn");
}
if(c == 0.) {
printf("Égalitén");
}
return EXIT_SUCCESS;
}

chez moi (amd64 FreeBSD 6.2) la sortie est

b = 0.0000000000000004440892098500626161694526672363281250000000000000
c = 0.0000000000000000000000000000000000000000000000000000000000000000
Classification
Égalité

En pratique, la seule réponse à laquelle on sait répondre (presque)
sans calcul est ``la valeur de la racine cubique est entière à epsilon
près'' (presque sans calcul parceque le coût d'une soustraction est
négligeable devant les autres opérations en virgule flottante).

Pour avoir une réponse catégorique, un moyen est effectivement
d'arrondir à l'entier le plus proche et d'élever au cube.

Pour plus de détails, lire le manuel de la bibliothèque math (man math
dans FreeBSD) et en particulier la rubrique fpclassify (man
fplclassifiy dans FreeBSD).
--
Amicalement,
Michaël

En news:86bq8yu23x.fsf@Llea.celt.neu, Michaël Grünewald va escriure:

En pratique, la seule réponse à laquelle on sait répondre (presque)
sans calcul est ``la valeur de la racine cubique est entière à epsilon
près'' (presque sans calcul parceque le coût d'une soustraction est
négligeable devant les autres opérations en virgule flottante).

Et sauf qu'il faut que tu calcules "epsilon" pour l'opération « racine
cubique »...
(OK, ce n'est nécessaire qu'une seule fois)

Pour avoir une réponse catégorique, un moyen est effectivement
d'arrondir à l'entier le plus proche et d'élever au cube.

Et cela ne marche que si tu ne dépasses pas l'espace de représentation en
flottant des nombres entiers, souvent 1<<53, c'est-à-dire environ 1<<18 au
cube. Et 1<<18 c'est 262144, ce n'est pas un nombre très grand, donc il vaut
mieux tenir compte de ce « biais » avant d'être catégorique.



Antoine

"Antoine Leca" <root@localhost.invalid> writes:

En news:86bq8yu23x.fsf@Llea.celt.neu, Michaël Grünewald va escriure:

En pratique, la seule réponse à laquelle on sait répondre (presque)
sans calcul est ``la valeur de la racine cubique est entière à epsil on
près'' (presque sans calcul parceque le coût d'une soustraction est
négligeable devant les autres opérations en virgule flottante).

Et sauf qu'il faut que tu calcules "epsilon" pour l'opération « racine
cubique »...
(OK, ce n'est nécessaire qu'une seule fois)

Je n'ai pas été très précis. Ce que je voulais dire exactement c'est
qu'on a une fonction g (racine cubique approchée) qui calcule la
fonction f (racine cubique) avec une erreur E, le tout vérifiant

| f(x) - g(x) | < E(x)

pour tout x dans le domaine de ces fonctions. Dans les cas favorables
on sait répondre à la question

``pour epsilon > 0, existe-t-il un nombre entier
epsilon-proche de f(x)?''

à l'aide de notre valeur approchée g(x) et de notre marge d'erreur
E(x).

Je suis loin de m'y connaître en calcul numérique, mais il me semble
que dans ce cas particulier on peut s'attendre à une erreur E(x) de la
forme O(Racine cubique de l'ordre de grandeur de x) où l'ordre de
grandeur est le morceau en 10^n dans l'écriture scientifique de x.



Pour en revenir à la question de l'OP, on peut dire que le calcul
numérique n'est pas sans difficulté, et qu'utiliser les nombres
machine pour tester des identités algébriques est délicat.

Pour avoir une réponse catégorique, un moyen est effectivement
d'arrondir à l'entier le plus proche et d'élever au cube.

Et cela ne marche que si tu ne dépasses pas l'espace de représentatio n en
flottant des nombres entiers, souvent 1<<53, c'est-à-dire environ 1<<18 au
cube. Et 1<<18 c'est 262144, ce n'est pas un nombre très grand, donc il vaut
mieux tenir compte de ce « biais » avant d'être catégorique.

Merci pour tes précisions.
--
Amicalement,
Michaël

Pour avoir une réponse catégorique, un moyen est effectivement
d'arrondir à l'entier le plus proche et d'élever au cube.

Et cela ne marche que si tu ne dépasses pas l'espace de représentation en
flottant des nombres entiers, souvent 1<<53, c'est-à-dire environ 1<<18 au
cube. Et 1<<18 c'est 262144, ce n'est pas un nombre très grand, donc il vaut
mieux tenir compte de ce « biais » avant d'être catégorique.

Je ne comprend pas trop la remarque. L'algorithme (arrondi élevé au
cube) ne pose pas de problème ici, non ? C'est juste que si la racine
cubique que l'on obtient est supérieure à ton 2^18, alors c'est que le
nombre entier qu'on a voulu traiter au départ n'a probablement pas pu
être représenté fidèlement sur un flottant ; mais enfin, ce n'est pas
lié à l'algo lui-même, c'est un problème annexe.

michaelgrunewald@yahoo.fr (Michaël Grünewald) writes:

Les calculs en virgule flottante sont des calculs approchés. Pour
tester l'égalité de a et de b on teste en réalité si la différence est
petite, au moyen du test

fpclassifiy(a - b) == FP_ZERO

Ce n'est certainement pas le but du fpclassify provenant de la norme C qui
est -- comme son nom l'indique -- d'indiquer la nature du flottant (NaN,
infini, etc...) Donc

fpclassify(a-b) == FP_ZERO

devrait a priori toujours etre equivalent a

a-b == 0.0

Je peux me gourer en ce qui concerne les cas particuliers (NaN...) mais ca
ne fait certainement pas un test d'egalite approximatif -- et l'interface
n'est pas bonne pour cela, un test d'egalite approximatif prendrait les
deux nombres en parametre (si je suis en train de manipuler deux nombres
proches de 1E-30, il ne sont pas d'office egaux).

A+

--
Jean-Marc
FAQ de fclc: http://www.isty-info.uvsq.fr/~rumeau/fclc
Site de usenet-fr: http://www.usenet-fr.news.eu.org

En news:fjloid$vk3$1@news.ciril.fr, Blaise Potard va escriure:

Pour avoir une réponse catégorique, un moyen est effectivement
d'arrondir à l'entier le plus proche et d'élever au cube.

Et cela ne marche que si tu ne dépasses pas l'espace de
représentation en flottant des nombres entiers, souvent 1<<53,
c'est-à-dire environ 1<<18 au cube. Et 1<<18 c'est 262144, ce n'est
pas un nombre très grand, donc il vaut mieux tenir compte de ce «
biais » avant d'être catégorique.

Je ne comprend pas trop la remarque. L'algorithme (arrondi élevé au
cube) ne pose pas de problème ici, non ? C'est juste que si la racine
cubique que l'on obtient est supérieure à ton 2^18, alors c'est que le
nombre entier qu'on a voulu traiter au départ n'a probablement pas pu
être représenté fidèlement sur un flottant ; mais enfin, ce n'est pas
lié à l'algo lui-même, c'est un problème annexe.

Tu me fais douter (et je n'aurais pas voulu y passer autant de temps).

Soit un nombre entier positif N, supérieur à 1<<DBL_MANTISSA (donc pas
exactement représentable comme double). Appelons n sa représentation en
double. On en extrait la racine cubique avec x=pow(n,1/3.), on arrondit à
l'entier le plus proche (rint, nearbyint, floor(x-.5), ceil(x+.5), comme
vous voulez, pour le moment je m'en moque) ce qui donne y, puis on repasse
au cube par deux multiplications z=y*y*y, est-on sûr de toujours obtenir
exactement la même représentation que n ?

Déjà, je suis quasiment certain que cela dépend des paramètres d'arrondi,
puisque la double multiplication va donner un résultat « inexact ». En
particulier, si les paramètres ne sont pas les mêmes entre l'arrondi de N à
n d'une part, et l'arrondi pour z d'autre part.
Mais comment peuvent-ils être différents ? En dehors d'une utilisation
explicite ou imprévues (genre un appel de bibliothèque au milieu ;) ) de
fesetround, une situation à problème potentiel est le cas où le premier
arrondi est fait par le compilateur, et le second est fait au moment de
l'exécution...

J'ai écrit un petit truc qui montre le problème en caricaturant :

#include <math.h>
#include <fenv.h>
#include <stdio.h>

int teste(int racine) {
fesetround(FE_TONEAREST); /* simule le comportement d'un compilateur */
double n=1.0*racine*racine*racine;
fesetround(FE_DOWNWARD); /* retour au mode exécution (supposé) */
double x=pow(n,1/3.);
double y=floor(x+.5);
double z=y*y*y;

if( n!=z || racine%20000==0 )
printf("n=%6d³ y-x=%8.6e (>>%4.2f) z=%20.0f == n -> %dn",
racine, y-x, y==x?0:-log((y-x)/racine)/log(2.), z, n==z);
return n==z;
}

int main() {
int X;

printf("«Petits» nombresn");
/* 208060 est juste en dessous de la racine cubique de 1<<53 */
for( X=1; X < 208060; ++X) teste(X);

printf("«Grands» nombresn");
for( X00000; X<300030; ++X) teste(X);
return 0;
}



J'en obtiens un sur 4 en erreur dans les « grands » nombres... et aucun
avant 208000.
Au passage, on peut voir aussi l'erreur sur la racine cubique (la question
originale) : pour ma libc, j'ai une erreur par défaut quasi-systématique de
3 chiffres binaires (erreur relative de -1>>50 sur x).


Antoine

Jean-Marc Bourguet <jm@bourguet.org> writes:

michaelgrunewald@yahoo.fr (Michaël Grünewald) writes:

Les calculs en virgule flottante sont des calculs approchés. Pour
tester l'égalité de a et de b on teste en réalité si la différ ence est
petite, au moyen du test

fpclassifiy(a - b) == FP_ZERO

Ce n'est certainement pas le but du fpclassify provenant de la norme C qui
est -- comme son nom l'indique -- d'indiquer la nature du flottant (NaN,
infini, etc...)

Aïe, je croyais pourtant me souvenir qu'il existait un moyen simple de
tester si des nombres sont ``à peu près égaux''.

Je peux me gourer en ce qui concerne les cas particuliers (NaN...) mais ca
ne fait certainement pas un test d'egalite approximatif -- et l'interface
n'est pas bonne pour cela, un test d'egalite approximatif prendrait les
deux nombres en parametre (si je suis en train de manipuler deux nombres
proches de 1E-30, il ne sont pas d'office egaux).

Argument imparable qui éclaire la débilité de ma proposition. Merci!
(Et désolé pour cette contribution décidément pas très pertinente .)
--
Cordialement,
Michaël

Le 11/12/07 22:41, dans <86lk80di42.fsf@Llea.celt.neu>, « Michaël
Grünewald » <michaelgrunewald@yahoo.fr> a écrit :

Aïe, je croyais pourtant me souvenir qu'il existait un moyen simple de
tester si des nombres sont ``à peu près égaux''.

FAQ § 11.5

--
Éric Lévénez -- <http://www.levenez.com/>
Unix is not only an OS, it's a way of life.