OVH Cloud OVH Cloud

Gros probleme étrange de precision pour les long double

25 réponses
Avatar
Cédric BILLOT
Voila, je suis sur un projet qui demande énormément de précision, je travail
donc en long double.
Je calcule la dérivé d'une fonction en un point grace a une methode simple
(f(x+h)-f(x-h))/2h) mais voila,
pour 80% des cas tout va bien; et parfois, je me retrouvais avec 0, alors
que la reponse du calcul fait sous Maple donnais un résultat.
J'ai épuré mon code autour de l'erreur jusqu'a avoir le code suivant.

long double h = 0.00000000000000000000001;
long double a = 10132.2;
if(a == a + h) cout<<"ERREUR DE CALCUL."<<endl;

et la!!!! c'est le drame!!!!

en effet a==a+h pourquoi?
Je ne pense pas pourtant arriver a la limite de précision du long double....

Quelqu'un a-t-il la solution a mon problème???

Merci d'avance,
@+Fab

10 réponses

1 2 3
Avatar
Horst Kraemer
On Thu, 1 Apr 2004 21:34:31 +0200, "Cédric BILLOT"
wrote:

Oui en effet; merci pour vos explication ;-)
Par contre comment je pourais m'en sortir pour avoir la meilleur
précision???
-Utilisé un h plutôt gros???
-Utilisé un h calculé en fonction de je ne sais pas trop koi??
En gros j'ai besoin d'une précision assez importante, et ma fonction qui
fait des dérivés partiels est utilisé plein de fois pour remplire une
matrice de taille 500*5 donc je ne peut pas me permettre non plus de perdre
trop de temps....
Est-ce qu'il existe un type de précision superieur au long double ???
Vous avez parlé du type Extended, je l'avais utilisé dans ma jeunesse, mais
je ne sais plus le faire (y-a-t-il un .h a inclure...).

Pour info j'utilise VC++ 6 et aussi Borland C++ 6 Builder (d'ailleur, c'est
surtout sous ce dernier que j'aimerais résoudre ce problème).
Je suis sur machine x86 sous Win2k.


Tu m'as mal compris. Il s'agit d'un problème numérique général, c.a.d.
du problème qu'une certaine méthode *mathématique* correcte ne peut
pas être employé pour un problème *numérique*.

Je te donne un exemple qui te convainvra peut-être. Tu sais que la
dérivé de la fonction exp(x) est exp(x). Maintenant on essaye de
calculer la dérivée de exp pour x=1 (qui est égale à exp(1)) par la
formule (f(1+h)-f(1-h))/(2h) pour de petits h. La loi mathématique dit
que pour h->0 ce quotient converge vers exp(1). La loi numérique des
flottants dit que le résultat s'éloigne de exp(1) si h tend vers 0:





Voila, je suis sur un projet qui demande énormément de précision, je
travail


donc en long double.
Je calcule la dérivé d'une fonction en un point grace a une methode
simple


(f(x+h)-f(x-h))/2h) mais voila,
pour 80% des cas tout va bien; et parfois, je me retrouvais avec 0,
alors


que la reponse du calcul fait sous Maple donnais un résultat.
J'ai épuré mon code autour de l'erreur jusqu'a avoir le code suivant.

long double h = 0.00000000000000000000001;
long double a = 10132.2;
if(a == a + h) cout<<"ERREUR DE CALCUL."<<endl;

et la!!!! c'est le drame!!!!

en effet a==a+h pourquoi?


Pour garantir que a+h > a il faut que h/a soit supérieur à
LDBL_EPSILON.



Je ne pense pas pourtant arriver a la limite de précision du long
double....



Chez moi LDBL_EPSILON vaut 1E-19 et ton h/a vaut environ 1E-27.

Le problème général est que dans un calcul numérique avec des
flottants le quotient (f(a+h)-f(a-h))/2h ne converge pas vers f'(a)
pour h->0 parce que la soustraction de deux valeurs trés proches
détruit la précision du résultat au lieu de l'augmenter.

--
Horst







Avatar
Horst Kraemer
On Thu, 1 Apr 2004 21:34:31 +0200, "Cédric BILLOT"
wrote:

Oui en effet; merci pour vos explication ;-)
Par contre comment je pourais m'en sortir pour avoir la meilleur
précision???
-Utilisé un h plutôt gros???
-Utilisé un h calculé en fonction de je ne sais pas trop koi??
En gros j'ai besoin d'une précision assez importante, et ma fonction qui
fait des dérivés partiels est utilisé plein de fois pour remplire une
matrice de taille 500*5 donc je ne peut pas me permettre non plus de perdre
trop de temps....
Est-ce qu'il existe un type de précision superieur au long double ???
Vous avez parlé du type Extended, je l'avais utilisé dans ma jeunesse, mais
je ne sais plus le faire (y-a-t-il un .h a inclure...).

Pour info j'utilise VC++ 6 et aussi Borland C++ 6 Builder (d'ailleur, c'est
surtout sous ce dernier que j'aimerais résoudre ce problème).
Je suis sur machine x86 sous Win2k.


Tu m'as mal compris. Il s'agit d'un problème numérique général, c.a.d.
du problème qu'une certaine méthode *mathématique* correcte ne peut
pas être employé pour un problème *numérique*.

Je te donne un exemple qui te convainvra peut-être. Tu sais que la
dérivé de la fonction exp(x) est exp(x). Maintenant on essaye de
calculer la dérivée de exp pour x=1 (qui est égale à exp(1)) par la
formule (f(1+h)-f(1-h))/(2h) pour de petits h. La loi mathématique dit
que pour h->0 ce quotient converge vers exp(1). La loi numérique des
flottants dit que le résultat s'éloigne de exp(1) pour de petits h si
h tend vers 0:

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

int main()
{
double d1 = 1e-7;
double d2 = 1e-9;
double d3 = 1e-11;
double d4 = 1e-13;
double d5 = 1e-15;

printf("%20.15enn", exp(1));
printf("%20.15en", (exp(1+d1)-exp(1-d1))/(2*d1));
printf("%20.15en", (exp(1+d2)-exp(1-d2))/(2*d2));
printf("%20.15en", (exp(1+d3)-exp(1-d3))/(2*d3));
printf("%20.15en", (exp(1+d4)-exp(1-d4))/(2*d4));
printf("%20.15en", (exp(1+d5)-exp(1-d5))/(2*d5));

return 0;
}


Résultat

2.718281828459045e+00 ( C'est le "bon" résultat )

2.718281829036965e+00 ( h = 1E-7 )
2.718281824783640e+00 ( h = 1E-9 )
2.718291021962249e+00 ( h = 1E-11 )
2.717753322736827e+00 ( h = 1E-13 )
2.840067590825779e+00 ( h = 1E-15 )


La méthode de calculer la dérivée par cette formule pour de petits h
est tout simplement *numériquement fausse*.

--
Horst

Avatar
Jean-Marc Bourguet
"Cédric BILLOT" writes:

Oui en effet; merci pour vos explication ;-)
Par contre comment je pourais m'en sortir pour avoir la meilleur
précision???


Te renseigner sur l'analyse numerique. Desole, mais je ne crois pas
que des renseignements pris sur usenet puisse etre suffisant si tu as
la moindre exigences quant au resultat.

--
Jean-Marc
FAQ de fclc++: http://www.cmla.ens-cachan.fr/~dosreis/C++/FAQ
C++ FAQ Lite en VF: http://www.ifrance.com/jlecomte/c++/c++-faq-lite/index.html
Site de usenet-fr: http://www.usenet-fr.news.eu.org

Avatar
kanze
Jean-Marc Bourguet wrote in message
news:...
"Cédric BILLOT" writes:

long double h = 0.00000000000000000000001;
long double a = 10132.2;
if(a == a + h) cout<<"ERREUR DE CALCUL."<<endl;
et la!!!! c'est le drame!!!!
en effet a==a+h pourquoi?
Je ne pense pas pourtant arriver a la limite de précision du
long double....


Si 17-18 chiffres significatifs et il t'en faudrait 28...
A bon, pourquoi???

La Mantisse est codé sur combien de bit pour le long double?


ca depend, 52 64 ou plus (si les long double sont sur 128 bits, je ne
connais ca que sur des stations Unix).


Ou des gros IBM ? Théoriquement, au moins, il peut y avoir moins
aussi ; la norme n'exige que 10 chiffres décimaux.

Ceci dit, je crois qu'il y a un problème plus profond. Parce que même
avec 100 bits ou plus de mantisse, il va y avoir une perte de précision
notable, et un résultat sans signification. Si l'algorithme fait qu'on
ajoute des valeurs avec des magnitudes très différentes, il faut changer
l'algorithme.

--
James Kanze GABI Software mailto:
Conseils en informatique orientée objet/ http://www.gabi-soft.fr
Beratung in objektorientierter Datenverarbeitung
9 place Sémard, 78210 St.-Cyr-l'École, France, +33 (0)1 30 23 00 34




Avatar
Jean-Marc Bourguet
writes:

Jean-Marc Bourguet wrote in message
La Mantisse est codé sur combien de bit pour le long double?


ca depend, 52 64 ou plus (si les long double sont sur 128 bits, je ne
connais ca que sur des stations Unix).


Ou des gros IBM ?


Possible, j'ai jamais utilise.

--
Jean-Marc
FAQ de fclc++: http://www.cmla.ens-cachan.fr/~dosreis/C++/FAQ
C++ FAQ Lite en VF: http://www.ifrance.com/jlecomte/c++/c++-faq-lite/index.html
Site de usenet-fr: http://www.usenet-fr.news.eu.org



Avatar
Marco
Cédric BILLOT wrote:
Oui en effet; merci pour vos explication ;-)
Par contre comment je pourais m'en sortir pour avoir la meilleur
précision???
-Utilisé un h plutôt gros???
-Utilisé un h calculé en fonction de je ne sais pas trop koi??
En gros j'ai besoin d'une précision assez importante, et ma fonction qui
fait des dérivés partiels est utilisé plein de fois pour remplire une
matrice de taille 500*5 donc je ne peut pas me permettre non plus de pe rdre
trop de temps....
Est-ce qu'il existe un type de précision superieur au long double ???
Vous avez parlé du type Extended, je l'avais utilisé dans ma jeunes se, mais
je ne sais plus le faire (y-a-t-il un .h a inclure...).


si tu veux utiliser un h plus gros, utilise des expressions "d'ordre
supérieur pour évaluer la dérivée", qui prennent les valeurs de u (x),
u(x+h), u(x-h) mais aussi u(x+2*h), u(x-2*h). Tu peux également utilise r
des points non régulièrement distribués et avoir une précision en core
meilleur (points de Gauss-Lobatto :-) ). Moi, j'aborderais la question
sous un angle éléments finis tu vois, ou tu dois résoudre v = u'.
L'avatange des éléments finis, c'est que la montée en ordre se fait de
manière immédiate, ce qui n'est pas le cas des différences finies. La je
rejoins une autre personne, tu dois lire des bouquins d'analyse
numérique, qui traitent du calcul de la dérivée. Si tu as une expre ssion
explicite de ta fonction, dérive la, ça coute pas cher et c'est ce qu i
marche le mieux !
500*5, c'est de la blague, petite matrice, puis si tu arretais windaube
tu t'en rendrais compte...

Avatar
Cédric BILLOT
Merci beaucoup pour toute ces precisions .
Je peux en effet calculer formellement la derive et l evaluer en la valeur
souhaite , mais alors mon programme
n est plus autonome.
Je vais faire des recherche dans des livres d analyse numerique pour voir s
il y a d autre solution ....
Merci beaucoup
@+Ced
"Marco" a écrit dans le message de
news:c4jait$m8v$
Cédric BILLOT wrote:
Oui en effet; merci pour vos explication ;-)
Par contre comment je pourais m'en sortir pour avoir la meilleur
précision???
-Utilisé un h plutôt gros???
-Utilisé un h calculé en fonction de je ne sais pas trop koi??
En gros j'ai besoin d'une précision assez importante, et ma fonction qui
fait des dérivés partiels est utilisé plein de fois pour remplire une
matrice de taille 500*5 donc je ne peut pas me permettre non plus de
perdre

trop de temps....
Est-ce qu'il existe un type de précision superieur au long double ???
Vous avez parlé du type Extended, je l'avais utilisé dans ma jeunesse,
mais

je ne sais plus le faire (y-a-t-il un .h a inclure...).


si tu veux utiliser un h plus gros, utilise des expressions "d'ordre
supérieur pour évaluer la dérivée", qui prennent les valeurs de u(x),
u(x+h), u(x-h) mais aussi u(x+2*h), u(x-2*h). Tu peux également utiliser
des points non régulièrement distribués et avoir une précision encore
meilleur (points de Gauss-Lobatto :-) ). Moi, j'aborderais la question
sous un angle éléments finis tu vois, ou tu dois résoudre v = u'.
L'avatange des éléments finis, c'est que la montée en ordre se fait de
manière immédiate, ce qui n'est pas le cas des différences finies. La je
rejoins une autre personne, tu dois lire des bouquins d'analyse
numérique, qui traitent du calcul de la dérivée. Si tu as une expression
explicite de ta fonction, dérive la, ça coute pas cher et c'est ce qui
marche le mieux !
500*5, c'est de la blague, petite matrice, puis si tu arretais windaube
tu t'en rendrais compte...

Avatar
Nicolas
Je vais faire des recherche dans des livres d analyse numerique...
Je pense que c'etais la première chose à faire!! personne n'évalue de

dérivée
comme tu l'as fait à cause de l'instabilité numérique.
Tu peux essayer entre autre approximation les polynome de
chebyshev.
http://www.library.cornell.edu/nr/bookcpdf/c5-8.pdf

Et n'oublie pas que tu as le Numerical Receipes in C qui t'attends avec ses
algorithmes
performants et son code sources diponible... voir "derivative"

http://www.library.cornell.edu/nr/ (le bouquin en ligne en PDF)
http://www.nr.com/ (le site du bouquin )

Bon travail.

--
Nicolas

Avatar
Nicolas
ptain que de fautes, honte a moi...
Avatar
kanze
Pierre Maurette <maurette.pierreATfree.fr@> wrote in message
news:...
"Cédric BILLOT" typa:

long double h = 0.00000000000000000000001;
long double a = 10132.2;
if(a == a + h) cout<<"ERREUR DE CALCUL."<<endl;
et la!!!! c'est le drame!!!!
en effet a==a+h pourquoi?
Je ne pense pas pourtant arriver a la limite de précision du long
double....


Si 17-18 chiffres significatifs et il t'en faudrait 28...
A bon, pourquoi???

La Mantisse est codé sur combien de bit pour le long double?


Pour info, formats des "réels" de la FPU Intel:
Single Double Extended
bits total 32 64 80
bits exposant 8 11 15
bits mantisse 23 52 63
^^


Juste un détail, mais pour l'extended, c'est bien 64.

précision (bits) 24 53 64

La mantisse est positive, l'exposant "biasé", c.a.d en gros signé
centré sur 0.


C-à-d en gros, non signé, mais on soustrait une valeur constante pour en
obtenir la valeur réele.

Si vous faites un contrôle du total des bits, n'oubliez pas le bit du
cygne. Et il n'y a pas d'erreur, l'extended a un format un peu
particulier.


Je dirais plutôt que ce sont les float et les double qui ont un format
un peu particulier:-). C'est eux où on ne représente pas le bit de poids
fort de la mantisse (parce que dans un nombre normalisé, il serait
toujours un).

Mais ce que tu documentes, c'est le hardware. Rien n'oblige un
compilateur de l'utiliser, et en fait, avec VC++ 6.0, j'ai
sizeof(double) == sizeof(long double). Les deux vaut 8, c-à-d donc que
le compilateur utilise le type double de la FPU et pour double et pour
long double.

--
James Kanze GABI Software mailto:
Conseils en informatique orientée objet/ http://www.gabi-soft.fr
Beratung in objektorientierter Datenverarbeitung
9 place Sémard, 78210 St.-Cyr-l'École, France, +33 (0)1 30 23 00 34




1 2 3