Dr. Goulu

“The Right Way to Calculate Stuff” est une page pour informaticiens dans mon genre : elle contient des petits “snippets” de code utile pour contourner certains pièges tendus par les maths et la géométrie.

Par exemple, si on veut calculer sin(x)/x, il faut faire attention à ce qui se passe pour x=0, sinon on obtiendra une erreur de “division by zero” au lieu du résultat attendu, qui est 1. Un bon programmeur prévoira le coup en écrivant (en C++) :

if (x == 0.) return 1. else return sin(x)/x;

mais ce n’est peut-être pas suffisant. Si x est très petit, la division de deux nombres très petits approximés par nos chers ordinateurs a des chances d’être assez fausse. Pour des raisons expliquées en détail sur la page, il faut s’assurer en fait que x² est encore un nombre numériquement assez grand par rapport à 1, donc un très bon programmeur écrira ceci :

if (1. + x*x == 1.) return 1. else return sin(x)/x;

ce qui parait mathématiquement absurde puisqu’on imagine pas que 1+x² puisse valoir 1 pour x≠0. Mais informatiquement, c’est très possible.

L’article donne ainsi une série de recettes pour calculer des expressions dangereuses, par exemple :

• 1-cos(x) → 2*sin(x/2)²
• acos(1-x) → 2*asin(sqrt(x/2))
• et bien d’autres

Mais les recettes les plus intéressantes (en tout cas pour moi) concernent la géométrie.

Pour calculer l’angle entre 2 vecteurs u et v, l’ “approche pachyderme” comme disait Burgi consiste à évaluer acos(dot(u,v)), dot étant la fonction renvoyant le produit scalaire des 2 vecteurs. mais si u ou v ne sont pas parfaitement normalisés et que dot renvoie un résultat un pouillème plus grand que 1 : crash.

une deuxième approche est: 2*atan2(abs(u-v),abs(u+v)). Pas mal mais couteux en temps de calcul. La meilleure implémentation est

if (dot(u,v) < 0.)

  return M_PI - 2*asin(abs(-v-u)/2)

else

  return 2*asin(abs(v,u)/2);

La partie vraiment intéressante concerne l’interpolation sphérique linéaire. Comme c’est assez spécifique à la 3D, je publie donc la suite de l’article  sur mon blog 3Dmon