CodeFR

Articles

Documentation

Algorithmes tri

Graphique

Langages

Systeme et Assembleur

Post mortem de projets

Programmes Hello World

Recherche

Pin ups

Juin 2007

Dernières modifications

1)	HomePage
2)	6502
3)	Byte pair encoding
4)	Documentation
5)	DataCompression
6)	Comb Sort
7)	Shell Sort
8)	Quick Sort
9)	Arbre rouge-noir
10)	Bogo Sort
11)	Gnome Sort
12)	polycubic

Mégaphone

hixcod, 12:10, Jan 09, 2010: Happy New year !

MooZ, 09:41, Jan 05, 2010: Bonne année! (erm...)

MooZ, 13:51, Jan 05, 2009: Bonne année!

MooZ, 07:19, Jan 04, 2008: Diantre! Je me suis fait doublé moi aussi. Bonne année quand même (et il faut que je trouve une pinup).

Axioplase, 15:39, Jan 03, 2008: Merde, Froggy m'a devancé sur 2008 !

Froggy, 22:21, Jan 01, 2008: Super mega cool

MooZ, 12:34, Oct 27, 2007: Yeah! Seven!

Seven, 16:41, Oct 18, 2007: Prout !

MooZ, 09:47, Mars 05, 2007: Feed the wiki comme dirait l'autre!

hixcod, 19:26, Mars 04, 2007: et moi donc

En lire plus…

similairecommentaire

Polynôme cubique

La manière la plus répandue de calculer les racines réelles d'une équation du troisième degré est d'utiliser la technique dite trigonométrique (parce qu'il y a des cosinus qui se battent en duel).

p(x) = ax**3 + bx**2 + cx + d = 0

Comme nous sommes d'immondes coders fainéants, je vais vous faire grâce de l'explication et surtout de la démonstration. Si vous voulez vraiment savoir d'où sort tout ça, je vous conseille d'aller lire l'article correspondant dans mathworld (cache).

#include <math.h>

int cubicSolver (double a, double b, double c, double d, double*x0, double* x1, double* x3) {
    double d, q, r, s, t;
    /* Tout d'abord on 'normalise' le polynome */
    b /= a;
    c /= a;
    d /= a;

    /* On substitue x = y - a / 3 afin d'éliminer le terme au carré pour avoir une équation */
    /* de la forme: x^3 + px + q = 0 */
    q = ((3.0 * c) - (b * b)) / 9.0;
    r = ((9.0 * c * b) - (27.0 * d) - (2.0 * b * b * b)) / 54.0; 

    /* Formule de Cardano */
    d = (q * q * q) + (r * r); 

    if(d > 0.0f) {
        /* On a une seule solution */
        d = sqrt(d);

        /* cbrt(x) calcule la racine cubique réelle de x */
        s = cbrt(r + d);
        t = cbrt(r - d);
      
        *x0 = *x1 = *x2 = -(b / 3.0) + s + t;   
      
        return 1;
    }

    if (d < 0.0f) {
        /* 3 racines réelles disctinctes */
        double theta = acos(r / sqrt(-(q * q * q)));
        s = 2.0 * sqrt(-q);
        *x0 = (s * cos(theta / 3.0)) - (b / 3.0);
        *x1 = (s * cos((theta + (2.0 * M_PI)) / 3.0)) - (b / 3.0);
        *x2 = (s * cos((theta - (4.0 * M_PI)) / 3.0)) - (b / 3.0);
 
        return 3;
    }

    /* 2 racines réelles */
    s = cbrt(r);
    *x0 = -(b / 3.0) + (2.0 * s);   
    *x1 = *x2 = -(b / 3.0) - s;

    return 2;
}

Attaquons maintenant un peu d'optimisation. Plutôt que de les balancer par importance, je vais suivre l'algo et les ajouter au fur et à mesure.

La première chose optimisable est la 'normalisation', en effet, on divise 3 fois par a. Une division prenant bien plus de temps qu'une multiplication, on peut remplacer

 b /= a;
 c /= a;
 d /= a;

par

 double inv_a = 1.0 / a;
 b *= inv_a;
 c *= inv_a;
 d *= inv_a;

On peut ensuite remarquer que l'expression ''b / 3.0"" revient assez souvent :

si une seule racine : 1 fois ;
3 fois s'il y a 3 racines ;
sinon 2 fois.

mais cette expression se retrouve aussi au carré dans la ligne

  q = ((3.0 * c) - (b * b)) / 9.0;

en utilisant une variable temporaire, on peut donc encore réduire le nombre de division. on obtient à la place de "q = ....' :

  double b_3 = b / 3.0;
  q = (c / 3.0) + b_3 * b_3;

et on remplace tous les '(b / 3.0)' par 'b_3'

Il reste à présent une dernière optimisation qui ne nuit pas trop à la lisibilité de l'algo. Elle se trouve dans le cas où il y a 3 racines.

  double theta = acos(r / sqrt(-(q * q * q)));
  s = 2.0 * sqrt(-q);
  *x0 = (s * cos(theta / 3.0)) - (b / 3.0);
  *x1 = (s * cos((theta + (2.0 * M_PI)) / 3.0)) - (b / 3.0);
  *x2 = (s * cos((theta - (4.0 * M_PI)) / 3.0)) - (b / 3.0);

se réécrit (et en prenant en compte l'optimisation précédente):

  double theta = acos(r / sqrt( -(q*q*q))) / 3.0;
  s = 2.0 * sqrt(-q);
  *x0 = (s * cos(theta)) - b_3;
  *x1 = (s * cos( theta + (2.0 * M_PI) / 3.0) - b_3;
  *x2 = (s * cos( theta - (4.0 * M_PI) / 3.0) - b_3;

Créé par: MooZ dernière modification: Mardi 03 of Mars, 2009 [11:52:46 UTC] par MooZ

Connexion