OVH Cloud OVH Cloud

Calculs du point de vue de la norme

57 réponses
Avatar
jul
Salut,

Je programme une application qui fait un fort usage des calculs avec
des valeurs d=E9cimales (=E0 virgule) avec les types float ou double. Il
ne me semble pas que ces calculs sont pris directement en charge par
le(s) processeur(s). Ma question est donc la suivante : la norme
impose-t-elle l'identit=E9 des r=E9sultats de calculs entre diff=E9rents
compilateurs ? Si oui, au moyen de quel(s) crit=E8re(s) ?

Je pr=E9cise pour bien me faire comprendre. Il facile d'exiger qu'un
compilateur calcule de mani=E8re correcte sur les entiers (1+1 doit
toujours =EAtre =E9gal =E0 2), ce qui est en soi une garantie
(ext=E9rieure) de l'identit=E9 des r=E9sultats. Mais lorsque le calcul a
lieu sur des nombres d=E9cimaux, le r=E9sultat est souvent arrondi. Il
n'y a pas alors de "r=E9sultat juste" (ex. 1/3 n'a pas d'expression
d=E9cimale finie juste, m=EAme si certaines sont meilleures que
d'autres). Comment sait-on que la mani=E8re de calculer et d'arrondir
sera la m=EAme sur chaque compilateur ?

Merci pour vos =E9claircissements.
Jul.

10 réponses

2 3 4 5 6
Avatar
James Kanze
Arnaud Meurgues wrote:
kanze wrote:
Mais ce n'était qu'une remarque parenthètique, parce que ça ne change
rien au problème. Grosso modo, si je comprends ce que tu veux dire,
c'est que si j'ajoute deux nombres de 6 chiffres, il faut que
j'utilise l'arithmétique à 12 chiffres, et que je garde le résultat
sur 12 chiffres (minimum, évidemment). Que la précision nécessaire
est la somme des précisions en entrée.



Là, je ne suis plus certain de suivre. La précision d'un flottant,
c'est la précision de sa mantisse, il me semble (qui correspond au
nombre de chiffres significatifs). Mais si l'on additionne deux
flottants avec un même nombre de chiffres significatifs, mais des
exposants très éloignés, on fait bien plus que doubler le nombre de
chiffres significatifs nécessaires pour conserver la précision voulue.


Du coup, je ne comprends pas pourquoi il faudrait particulièrement 12
chiffres significatifs pour additionner deux nombres ayant 6 chiffres
significatifs.


Moi non plus, mais c'est ce qu'il a dit. (À mon avis, si j'ajoute 10^-20
à 10^20, même si les deux n'ont que 6 chiffres de précision, il me faut
bien 40 chiffres au moins pour ne pas avoir une perte de précision.)

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


Avatar
Sylvain
James Kanze wrote on 20/05/2006 02:48:

Du coup, je ne comprends pas pourquoi il faudrait particulièrement 12
chiffres significatifs pour additionner deux nombres ayant 6 chiffres
significatifs.


Moi non plus, mais c'est ce qu'il a dit. (À mon avis, si j'ajoute 10^-20
à 10^20, même si les deux n'ont que 6 chiffres de précision, il me faut
bien 40 chiffres au moins pour ne pas avoir une perte de précision.)



là vous faites exprès de ne pas comprendre non ?

j'avais répondu 2h47 avant ce post sur la raison _particulière_ pour
laquelle cela était appliquable ici.

s'il faut redétailler:
- l'exemple qui prétends caractériser des différences de précisions par
des sommes algorithmiquement différentes est biaisé car les valeurs
choisies sont telles que l'on peut toujours réaliser le calcul sans
perte de précision (autre que les arrondis normaux),
- si l'exemple avait sommé 1 millions de nombres dans [0..10^3[ ou
[0..10^6[, il aurait été plus pertinent,
- sommer 1 milliards de valeurs dans [0..1[ aurait également montré des
pertes significatives de précision - d'où le SumPart (oui illustré dans
mon post sur 10^6 éléments, et alors?); j'ai bien noté vos réponses (et
je ne chercherais pas à vous contredire, vous devez avoir nécessairement
raison, et c'est la merde que vous avez noté).

c'est également parce que l'exemple est biaisé, qu'un calcul en double
donne directement le "bon résultat" qlq que soit l'algo et que les
arrangements listés sont possibles.

pour répondre à ce post, en math. appl. comme en physique comme en
info., 10^20 + 10^-20 font 10^20 point barre.

si l'eventuel sous entendu est comment sommer de 10^14 à 10^20 valeurs
proches de 10^-20 (moins de valeurs donnerait une somme négligeable
devant 10^20 pour un résultat float) et qlq valeurs proches de 10^20;
vous aurez compris même sans cette réponse qu'il convient de sommer
entre elles les valeurs en 10^-20 (et pas pas un SumUp) avec d'y ajouter
les valeurs en 10^20.

notons qu'il est aussi possible (facilement) d'écrire une classe
réalisant les opérations de base (+-*/) en précision arbitraire.

Sylvain.


Avatar
Sylvain
jul wrote on 17/05/2006 12:28:

C'est moins le problème d'identité dans le code (utiliser = > après deux calculs identiques) que celui de garantir la même sortie
entre deux exécutables compilés à partir du même code source,
sachant qu'une faible différence dans le résultat d'un calcul peut
entraîner de grosses différences dans ladite sortie.


il existe des problèmes numériques pour lesquels tu pourras le garantir
et d'autres pour lesquels ce sera difficile voire impossible.

préambule: la garantie du "même résutat" entre différents portages du
même code peut revétir un aspect psychologique: le client "croit" le
résultat (quel qu'il soit) s'il est constant, et se prend à douter des
compétences du fournisseur si les résultats changent entre 2 versions du
soft (Wintel vs gnu/linux) ou 2 runs sur 2 architectures différentes
(Intel vs AMD ou Sparc).

dans un tel cas, l'objectif du programme est ""simplement"" de minimiser
les écarts entre architectures/compilateurs et de fournir (d'afficher,
imprimer, tracer) que la partie sure et commune des résultats; ce qui
signifie fournir des résultats avec une précision moindre que celle
utilisée pour les calculs - j'ai indiqué dans ce fil qu'il est possible
(pour certains calculs) d'obtenir des résultats identiques à 10^-8, pour
autant fournir des résultats en 10^-4 peut être suffisant.

si la précision à fournir en sortie doit elle maximale (proche ou égale
à la précision du type et/ou du FPU), vous devrez coder vos expressions
de manière la plus explicite et minimiser (voire supprimer) toute
optimisation (réorganisation) de codes; le résultat sera "moins juste"
(que la valeur théorique) mais plus constant entre environnements.
(ainsi (a + b/c) / (a - b/c) exécuté entièrement via des registres
internes est proche du résultat théorique mais dépend du FPU; alors que
e=b/c; r=(a+e)/(a-e); est "moins précis" (plus d'écart par rapport à la
valeur théorique) mais dépendra moins du FPU et du compilateur.

finalement, la constance des résultats varie fortement selon le type de
calcul; pour de l'algèbre linéaire, obtenir les mêmes résultats avec
différentes config. reste possible; pour des calculs trigonométriques
et/ou log/exp, cela sera plus difficile car ces fonctions sont tabulées
et tous les chips ne contiennent pas les mêmes tables; ici la précision
rendue devra être plus limitée.

Sylvain.

Avatar
Jean-Marc Bourguet
Sylvain writes:

pour des calculs trigonométriques et/ou log/exp, cela sera
plus difficile car ces fonctions sont tabulées et tous les
chips ne contiennent pas les mêmes tables; ici la
précision rendue devra être plus limitée.


Il est vrai que les processeurs qui en disposent (il me
semble que c'est plutôt une spécificité des x86 mais je n'ai
pas été vérifier) et les implémentations soft ne
garantissent en général pas la précision des calculs
trigonométriques (être précis sur le dernier bit peut
demander beaucoup d'efforts, voir
http://perso.ens-lyon.fr/jean-michel.muller/Intro-to-TMD.htm)
et qu'elle est parfois nettement inférieure à ce qui est
représentable.

Mais à ma connaissance, les tables ne sont pas la manière
préférée d'implémenter ces fonctions. Deux techniques sont
couremment utilisées, CORDIC
(http://orochoir.club.fr/Maths/cordic.htm) et la réduction
dans un intervalle par des propriétés mathématiques suivie
de l'utilisation d'un polynome d'approximation.

A+

--
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
Jean-Marc Bourguet
Sylvain writes:

Jean-Marc Bourguet wrote on 17/05/2006 15:19:
Sylvain writes:

tes routines sont intéressantes ! de mon point de vue, elles sont toutes
inadaptées et maladroites (ce n'est pas une provocation!).
En quoi sont-elles inadaptees a l'objectif poursuivi: analyser

l'erreur introduite par differentes manieres de faire la somme
(definie comme etant la difference entre le resultat exact arrondi a
la precision d'un float et le resultat calculer)?


pour cet objectif, elles sont très bien... si utilisées correctement [1]

je me plaçais dans un cadre de traitement numérique et commentais le main()
de James qui ne contient que des SumXXX<float>; là je maintiens que un
SumUp<double> est suffisant.


Ah, un marteau est inadapté quand on veut visser.
D'accord. Quand je vois du code marqué "bench", je regarde
ce que c'est sensé mettre en évidence avant de le critiqué.

Je n'ai pas vu à quoi le [1] faisait référence.

dans tes benchs tu accumules 1.000.000 de valeurs
flottantes de type float, or un float à une précision
de 10^-6, en sommer 10^6 va obligeatoirement faire
perdre toute précision de l'élément ajouté par rapport
à l'accu courant - ie va générer un débordement de
précision et revient à faire un inf + epsilon en
espérant trouver l'epsilon dans le nouveau résultat.


Je ne comprends pas.


tu ne comprends pas quoi ?


Le raisonnement qui te fait penser que la somme d'un million
de nombres en flottant simple précision fera perdre
*nécessairement* tous les chiffres significatifs. J'attends
un seul jeu de données pour lequel c'est le cas avec
sumCompensated, pour ne pas parler de sumPriority.

que le nombre de chiffres significatifs est limité ? que
(float)(500000. + 0.512345) vaut 500000. et non
500000.512345 ?


Je comprends l'idée mais l'exemple est mal choisi.

500000 a besoin de 19 bits, les floats IEEE ont 24 bits de
précision, donc il y a 5 bits de plus, et 500000.512345 sera
approximé avec 500000.5. Si je ne me trompe pas le flottant
qui suit est 500000.53125.

Dans le cadre donne ci-dessus, on peut prouver des
bornes superieures pour l'erreur introduite pour deux
des algorithmes (sumPriority, sumCompensated) qui sont
de l'ordre de quelques ULP (si j'ai bonne memoire, moins
d'un ULP pour sumCompensated, autrement dit, une des
valeurs qui encadrent le resultat exact).



Note si on veut me citer, la borne dépend du nombre
d'éléments sommés. Et pour ceux qui ne l'aurait pas
remarque, sumCompensated est un moyen de doubler la
précision (donc sumCompensated<float> est légèrement moins
bon que sumCompensated<double> parce qu'il a 48 bits de
précision tandis que les doubles en ont 53, le nombre de
bits affectés à l'exposant ne doublant pas et le signe
n'étant naturellement pas doublé).

le problème est que pour la valeur 0.512345 lorsqu'elle
est ajouté à 500000. c'est même son UFP qui se volatilise


L'important, ce n'est pas tant qu'une valeur soit en
pratique totalement ignorée, c'est est-ce que ça affecte le
résultat final.


Pour ceux que le problème de la sommation intéresse:

http://www.ti3.tu-harburg.de/paper/rump/RuOgOi06.pdf

A+

--
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
Sylvain
Jean-Marc Bourguet wrote on 20/05/2006 21:30:

Mais à ma connaissance, les tables ne sont pas la manière
préférée d'implémenter ces fonctions.


pour du hard (surtout les "petits" CPU) c'est la manière préférée; mais
concernant, par exemple, un Itanium avec toute sa quincaillerie, je n'ai
aucune info précise sur sa façon de faire.


Deux techniques sont couremment utilisées, CORDIC
(http://orochoir.club.fr/Maths/cordic.htm) et la réduction
dans un intervalle par des propriétés mathématiques suivie
de l'utilisation d'un polynome d'approximation.


les développements de Cordic et suite de Taylor sont les plus courantes,
le developpement par produit d'Euler existe aussi mais il converge
moins bien et est au final trop couteux.

en formulant ma "mise en garde" précédente contre les résultats de
fonctions trigo. j'avais en tête les fonctions inline de math.h; selon
le domaine (dont traitement signal), on utilise parfois des librairies
soft fournissant des implémentations bcp plus rapide pour ces fonctions
transcendantes (telle DCT), utiliser une telle librairie contribuerait à
maitriser la reproductivité (cross-environement) des résultats.

Sylvain.

Avatar
Jean-Marc Bourguet
Sylvain writes:

Jean-Marc Bourguet wrote on 20/05/2006 21:30:
Mais à ma connaissance, les tables ne sont pas la manière
préférée d'implémenter ces fonctions.


pour du hard (surtout les "petits" CPU) c'est la manière
préférée;


A partir du moment où un CPU a une unité de calcul flottant,
ce n'est plus un petit CPU, avec ou sans guillements à
petit. En particulier si on a eu la place pour mettre des
tables de fonctions trigonométriques.

Les architectures ayant les fonctions trigonométriques
cablées sont rares. J'écrivais dans mon message précédant
que je ne voyais que le x86 mais que je n'avais rien
vérifié. J'ai été voir rapidement ni POWER ni Sparc V9 ni
MIPS-64 n'en ont. Je cherche donc autre chose, mais ça
risque d'être un machin très anecdotique.

Pour du hard des grandes tables ce n'est jamais la méthode
préférée. On peut tabuler un peu puis partir avec des
itérations de rafinement, mais je ne connais pas de tels
algo pour les fonctions trigonométriques. Je continue donc
à penser que CORDIC -- qui utilise une petite table -- ou un
polynome bien choisi -- qui va utiliser a peu près autant de
coefficiants qu'il n'y a d'élément dans la table pour CORDIC
à précision égale -- sont les méthodes que j'emploierais.

mais concernant, par exemple, un Itanium avec toute sa
quincaillerie, je n'ai aucune info précise sur sa façon de
faire.


Je n'ai jamais vu d'info publiée la-dessus, mais je n'ai pas
réellement cherché.

Deux techniques sont couremment utilisées, CORDIC
(http://orochoir.club.fr/Maths/cordic.htm) et la réduction
dans un intervalle par des propriétés mathématiques suivie
de l'utilisation d'un polynome d'approximation.


les développements de Cordic


Voir CORDIC comme un développement est original. J'ai
plutôt tendance à le voir comme une réduction d'intervalle
où on fini par n'avoir plus qu'un point dans l'intervalle.

et suite de Taylor sont les plus courantes, le
developpement par produit d'Euler existe aussi mais il
converge moins bien et est au final trop couteux.


Je ne suis pas un expert en la matière, mais utiliser un
développement de Taylor me semblerait naïf comme approche.
Il doit avoir moyen de calculer des polynomes de moindre
degré avec la même erreur maximale en acceptant qu'il y ait
plus de points où l'évaluation donne cette erreur.

A+

--
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
Jean-Marc Bourguet
Jean-Marc Bourguet writes:

Dans le cadre donne ci-dessus, on peut prouver des
bornes superieures pour l'erreur introduite pour deux
des algorithmes (sumPriority, sumCompensated) qui sont
de l'ordre de quelques ULP (si j'ai bonne memoire, moins
d'un ULP pour sumCompensated, autrement dit, une des
valeurs qui encadrent le resultat exact).



Note si on veut me citer, la borne dépend du nombre
d'éléments sommés. Et pour ceux qui ne l'aurait pas
remarque, sumCompensated est un moyen de doubler la
précision (donc sumCompensated<float> est légèrement moins
bon que sumCompensated<double> parce qu'il a 48 bits de


sumUp<double> et pas sumCompensated<double> naturellement

précision tandis que les doubles en ont 53, le nombre de
bits affectés à l'exposant ne doublant pas et le signe
n'étant naturellement pas doublé).


Je complete, ca fait un peu plus que cela. Il est possible de batir
des cas ou sumCompensated<float> a une meilleure precision que
sumUp<double>.

A+

--
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
Jean-Marc Bourguet
Jean-Marc Bourguet writes:

mais concernant, par exemple, un Itanium avec toute sa
quincaillerie, je n'ai aucune info précise sur sa façon de
faire.


Je n'ai jamais vu d'info publiée la-dessus, mais je n'ai pas
réellement cherché.


New Algorithms for Improved Transcendental Functions on IA-64
Shane Story, Ping Tak Peter Tang
Trouve quelque part sur le site de l'IEEE.

A regarder rapidement ils utilisent une decomposition polynomiale par
morceau.

D'autres papiers sur le site de l'IEEE font reference a CORDIC et a
des polynomes.

Pour le moment, je n'ai que les x86 et les ia64 qui ont des fonctions
trigonometriques cablees.

A+


--
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
Sylvain
Jean-Marc Bourguet wrote on 22/05/2006 15:51:

Pour le moment, je n'ai que les x86 et les ia64 qui ont des fonctions
trigonometriques cablees.



oui, et ?! tu fais collec' ? si j'en trouve je te les mets de coté, promis.

la question (liée à "comment garantir des résutats identiques") ne
passe-t-elle pas plutôt par des commentaires sur les librairies
mathématiques disponibles sur le marché ?

Sylvain.

2 3 4 5 6