Analyse du bug Debian 153548 dans la glibc

Note: ce bug a été corrigé dans le CVS de la glibc le 2003-01-20.

Alors que je testais la bibliothèque GNU C (glibc) avec des pires cas pour l'arrondi correct, j'ai trouvé cinq valeurs (complètement) incorrectes sur ma machine PowerPC pour la fonction cosinus, parmi 1576 pires cas. Habituellement, quand il y a un bug, le nombre de valeurs fausses est bien plus grand. Alors, pourquoi seulement cinq?

D'abord, rappelons les valeurs incorrectes. Comme je l'ai dit dans mon rapport de bug, le problème se produit seulement dans un petit intervalle (aux alentours de 0,80 à 0,83):

Error: cos (line 33)
  [arg.]  x = 0.80190127184058835
  [corr]  y = 0.69534156199418473
  [libm]  z = 0.71867942238767868
Error: cos (line 34)
  [arg.]  x = 0.8107741500588278
  [corr]  y = 0.68893751943804493
  [libm]  z = 0.72482073253360613
Error: cos (line 35)
  [arg.]  x = 0.81808070960270574
  [corr]  y = 0.68362323105565981
  [libm]  z = 0.72983510326718326
Error: cos (line 36)
  [arg.]  x = 0.81909979079003903
  [corr]  y = 0.6828791149804222
  [libm]  z = 0.73053139174408888
Error: cos (line 37)
  [arg.]  x = 0.82926610298753134
  [corr]  y = 0.67541714390565444
  [libm]  z = 0.73743588312363018

Comme Lukas Geyer l'a remarqué (voir mon rapport de bug), c'est le sinus qui est calculé au lieu du cosinus.

J'ai alors regardé dans le CVS de la glibc et trouvé pour /cvs/glibc/libc/sysdeps/ieee754/dbl-64/sincos32.c,v:

revision 1.12
date: 2003/01/20 05:25:30;  author: roland;  state: Exp;  lines: +1 -1
2003-01-20  Segher Boessenkool  <...>

        * sysdeps/ieee754/dbl-64/sincos32.c (__mpcos): Really
        compute the cosine, not the sine, even if x > 0.8.

et le patch correspondant:

Index: sincos32.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/dbl-64/sincos32.c,v
retrieving revision 1.11
retrieving revision 1.12
diff -d -u -r1.11 -r1.12
--- sincos32.c  26 Aug 2002 22:40:37 -0000      1.11
+++ sincos32.c  20 Jan 2003 05:25:30 -0000      1.12
@@ -214,7 +214,7 @@
   __add(&a,&b,&c,p);
   if (x>0.8)
   { __sub(&hp,&c,&b,p);
-    __c32(&b,&a,&c,p);
+    __c32(&b,&c,&a,p);
   }
   else __c32(&c,&a,&b,p);     /* a = cos(x+dx)     */
   __mp_dbl(&a,&y,p);

Noter que ce bug ne peut se produire que sur des cas difficiles à arrondir puisqu'il se trouve dans la fonction __mpcos, qui calcule le cosinus en multiprécision (cette fonction n'est pas appelée dans les cas faciles, car une faible précision est suffisante pour décider de l'arrondi).

La borne inférieure 0,8 est maintenant évidente, mais qu'en est-il de la borne supérieure? Celle-ci peut être trouvée dans le fichier s_sin.c, où la fonction __cos est définie:

  else if (k < 0x3feb6000 ) {/* 2^-27 < |x| < 0.855469 */
[...]
    return (res==res+1.020*cor)? res : cslow2(x);

}    /*   else  if (k < 0x3feb6000)    */

  else if (k <  0x400368fd ) {/* 0.855469  <|x|<2.426265  */;

Ainsi, la borne supérieure est 0,855469, et avec un cas difficile à arrondir, la fonction cslow2 appelle __mpcos(x,0). Si l'argument est supérieur à 0,855469, alors __cos appelle __mpcos1(x), qui fait une réduction d'argument et appelle __c32 directement: la fonction erronée n'est pas appelée.