root/misc/dl.c

/*
  Renvoie la valeur approchée de fonctions mathématiques :
  Exponentiel, cosinus, sinus hyperbolique, arctan, ...
  grâce aux développement limités de ces fonctions.
  ---
  Si quelqu'un trouve comment calculer Tan, TanH, ArcCosH, qu'il
  m'envoie un email.
  ---
  ---
  Auteur: Victor Stinner aka haypo -- http://www.haypocalc.com
  Date de creation: 28 avril 2001
  Code source distribue sous licence GNU GPL
*/

#include <stdio.h>   // printf
#include <math.h>    // fabsl
#include <float.h>   // DBL_DIG
#include <stdbool.h> // type boolean: bool

//=============== DEFINE =======================================================


/* typedef (long double) NombreReel;
#define Abs(a) fabsl(a) */

typedef long NombreEntier;  // Utilise long pour les nombres entiers
typedef double NombreReel;  // Utiise double pour les nombres rééls

// Valeur absolue d'un nombre réel
NombreReel Abs (NombreReel x)
{
  return fabs(x);
}

#define Precision 12  // MAX = 20 en long double
#define PuissancePrecision 1e-12 // =1e-(Precision), MAX = -20 en long double
#define AfficheExemple(A,B) printf (A); printf ("(%.3g) = %.12g\n",x,B(x));

NombreReel Racine (NombreReel x)
{
  // Racine(0) = 0 (évite de tourne en rond avec cette algorithme)
  if (x==0) return 0;

  // Racine(0) = 1 (évite de tourne en rond avec cette algorithme)
  if (x==1) return 1;

  // Racine(nombre négatif) = i * Racine(abs(nombre négatif)) : ERREUR !
  // (ne supporte pas les nombres complexes)
  if (x<0) {
    printf ("ERREUR Racine : Résultat complexe ... ");
    return -1;
  }

  // Avec la méthode de Newton (x(n+1) = x(n) - f(x)/f'(x)) et la fonction
  // f(x)= x^2 - b (b étant le nombre dont on cherche la racine)
  // -> f'(x) = 2x, on trouve la suite qui tend vers la racine de b :
  // x(n+1) = x(n) - (x(n)^2 - b)/2x(n) = x(n) - x(n)/2 +b/2x(n)
  //        = (x(n) + b/x(n))/2
  NombreReel Racine,SauveX;

  // Le nombre dont on cherche la racine
  Racine = x;

  // Cherche une valeur proche de la racine (x=racine du nombre)
  // Pour 2<x, la racine est inférieure à x/2
  if (2<x) x /= 2;

  SauveX = x+1;
  while (PuissancePrecision<Abs(SauveX-x)) {
    SauveX = x;
    x = (x + Racine/x)/2;
  }

  return x;
}

//=============== COSINUS-SINUS-TANGENTE =======================================
NombreReel Cos (NombreReel x)
{
  NombreReel
    ValCosinus,
    Xpuissance,Factoriel,i,Quotient;

  // COS(x) = 1 -x^3/3! +x^5/5! -x^7/7 +...

  ValCosinus = 1;     // Premier terme de la suite
  x = x*x;            // Eleve x au carre (pour acc‚l‚rer le calcul de
                       // Xpuisscance)
  Xpuissance = x;     // = x^2, x^4, x^6, x^8, ...
  Factoriel = 2;      // = 2!,  4!,  6!,  8!,  ...
  i = 2;              // Iteration (utilis‚e pour le calcul du Factoriel)
  Quotient = Xpuissance/Factoriel;    // Calcule le quotient

  do {
    ValCosinus = ValCosinus -Quotient;                 // Soustrait un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    i = i +1; Factoriel = Factoriel*i; // Passe au factoriel suivant
    i = i +1; Factoriel = Factoriel*i; // (= *i*(i+1) ) et i=i+2

    ValCosinus = ValCosinus +Xpuissance/Factoriel;        // Ajoute un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    i = i +1; Factoriel = Factoriel*i; // Passe au factoriel suivant
    i = i +1; Factoriel = Factoriel*i; // (= *i*(i+1) ) et i=i+2
    Quotient = Xpuissance/Factoriel;    // Calcule le quotient
  } while (PuissancePrecision<Abs(Quotient));

  return ValCosinus;
}

NombreReel Sin (NombreReel x)
{
  NombreReel
    ValSin,
    Xpuissance,Factoriel,i,Quotient;

  // SIN(x) = x -x^3/3! +x^5/5! -x^7/7 +...

  ValSin = x;         // Premier terme de la suite
  Xpuissance = x*x*x; // = x, x^3, x^5, x^7, ...
  x = x*x;            // Eleve x au carre (pour acc‚l‚rer le calcul de
                       // Xpuisscance)
  Factoriel = 6;      // = 1!,3!,  5!,  7!,  ...
  i = 3;              // Iteration (utilis‚e pour le calcul du Factoriel)
  Quotient = Xpuissance/Factoriel;    // Calcule le quotient

  do {
    ValSin = ValSin -Quotient;  // Soustrait un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    i = i +1; Factoriel = Factoriel*i; // Passe au factoriel suivant
    i = i +1; Factoriel = Factoriel*i; // (= *i*(i+1) ) et i=i+2

    ValSin = ValSin +Xpuissance/Factoriel;        // Ajoute un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    i = i +1; Factoriel = Factoriel*i; // Passe au factoriel suivant
    i = i +1; Factoriel = Factoriel*i; // (= *i*(i+1) ) et i=i+2
    Quotient = Xpuissance/Factoriel;    // Calcule le quotient
  } while (PuissancePrecision<Abs(Quotient));
  return ValSin;
}

NombreReel Tan (NombreReel x)
{
  NombreReel
    ValSin,ValCos,
    Xpuissance,Factoriel,i,Quotient;

  // TAN(x) = SIN(x)/COS(x)
  //           x -x^3/3 +x^5/5 -x^7/7 +x^9/9 - ...
  //        = -------------------------------------
  //           1 -x^2/2 +x^4/4 -x^6/6 +x^8/8 - ...

  ValCos = 1;         // Premier terme de la suite de cosinus
  ValSin = x;         // Premier terme de la suite de sinus
  Xpuissance = x*x;   // = x^2, x^3, x^4, x^5, ...
  Factoriel = 2;      // = 2!, 3!, 4!, 5!,  ...
  i = 2;              // Iteration : 2,3,4,5,6,...
  Quotient = Xpuissance/Factoriel;    // Calcule le quotient

  do {
    ValCos = ValCos -Quotient;          // Soustrait un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    i = i +1; Factoriel = Factoriel*i; // Passe au factoriel suivant

    ValSin = ValSin -(Xpuissance/Factoriel);           // Soustrait un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    i = i +1; Factoriel = Factoriel*i; // Passe au factoriel suivant

    ValCos = ValCos +(Xpuissance/Factoriel);              // Ajoute un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    i = i +1; Factoriel = Factoriel*i; // Passe au factoriel suivant

    ValSin = ValSin +(Xpuissance/Factoriel);              // Ajoute un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    i = i +1; Factoriel = Factoriel*i; // Passe au factoriel suivant

    Quotient = Xpuissance/Factoriel;    // Calcule le quotient
  } while (PuissancePrecision<Abs(Quotient));

  return ValSin/ValCos;
}


//=============== EXPONENTIEL ==================================================
NombreReel Exp (NombreReel x)
{
  NombreReel
    ValExp,
    Xpuissance,q,f,i;

  // EXP(x) = x^0/0! +x^1/1! +x^2/2! +x^3/3! +... +x^i/i!

  ValExp = 1; // Premier terme de la suite
  Xpuissance = x;
  f = 1;
  i = 1;

  q = Xpuissance/f;
  do {
    ValExp = ValExp +q;
    Xpuissance = Xpuissance*x;
    i = i +1;
    f = f*i;
    q = Xpuissance/f;
  } while (PuissancePrecision<Abs(q));

  return ValExp;
}

NombreReel Ln (NombreReel x)
{
  NombreReel
    ValLn,
    Xpuissance,q,f,i;
  bool
    ResultatOppose;

  // LN(1+x) =    (0<x<=1)
  // x
  // -x^2/2
  // +x^3/3
  // -x^4/4
  // +x^5/5
  // -...

  // ln(1)=0 !
  if (x==1) return 0;

  // x<0 ?
  if (x<=0) {
    printf ("La fonction LN n'est pas definie sur ]-inf,0] !\n");
    return 0;
  }

  // 1<x ? Calcule -ln(1/x) dans ce cas là
  if (1<x) {
    x = 1/x;
    ResultatOppose = true;
  } else ResultatOppose = false;

  // On calcule Ln(1+x), donc on fait -1
  x--;

  // Initialise les variables
  ValLn = 0;
  Xpuissance = x;  // x,x^2,x^3,x^4,...
  i = 1;

  q = Xpuissance/i;
  do {
    ValLn += q;
    Xpuissance *= x;
    i++;

    ValLn -= Xpuissance/i;
    Xpuissance *= x;
    i++;

    q = Xpuissance/i;
  } while (PuissancePrecision<Abs(q));

  if (ResultatOppose)
    return (-ValLn);
  else
    return ValLn;
//  return ResultatOppose? -ValLn:ValLn;
}

//=============== ARCSIN, ARCCOS, ARCTAN =======================================
NombreReel ArcSin (NombreReel x)
{
  NombreReel
    ValArcSin,k,xp,q;
  NombreEntier si1,si2;

  // ArcSin (x) =
  //  x
  // +x^3/3* (1/2)
  // +x^5/5* (1*3)/(2*4)
  // +x^7/7* (1*3*5)/(2*4*6)
  // +...

  ValArcSin = x;
  k = 1;
  xp = x;
  si1 = 1;
  si2 = 2;
  x = x*x;      // Eleve x au carr‚
  q = 1;

  do {
    k = k*si1/si2;
    si1 += 2;
    si2 += 2;
    xp = xp*x;
    q = xp*k/si1;
    ValArcSin = ValArcSin +q;
  } while (PuissancePrecision<Abs(q));

  return ValArcSin;
}

NombreReel ArcCos (NombreReel x)
{
  // ArcCos(x) = pi/2 -ArcSin(x)

  return M_PI/2 -ArcSin(x);
}

NombreReel ArcTan (NombreReel x)
{
  NombreReel
    ValArcTan,
    Xpuissance,Denominateur,Quotient;

  // ARCTAN(x) = x -x^3/3 +x^5/5 -x^7/7 +...
  // --> POUR 0 <= x < 1

  if (1<=x) {
    printf (" <-OVERFLOW-> ");
    return 0;
  }

  // Initialise les variables
  ValArcTan    = 0;  // Premier terme de la suite
  Xpuissance   = x;  // = x, x^3, x^5, x^7, ...
  Denominateur = 1;  // = 1, 3, 5, 7, ...
  Quotient = Xpuissance/Denominateur;   // Quotient des deux
  x = x*x;   // Eleve x au carre (pour optimiser le calcul de Xpuissance

  // S'arrête quand on a atteind une précision suffisante
  do {
    ValArcTan = ValArcTan +Quotient;// Ajoute un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    Denominateur = Denominateur+2;  // Passe au d‚nominateur suivant (+ 2)

    ValArcTan = ValArcTan -Xpuissance/Denominateur; // Soustrait un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    Denominateur = Denominateur+2;  // Passe au d‚nominateur suivant (+ 2)
    Quotient = Xpuissance/Denominateur; // Calcule le quotient
  } while (PuissancePrecision<Abs(Quotient));

  // Renvoi la valeur de l'arctan(x)
  return ValArcTan;
}
//=============== COSH, SINH, TANH =============================================
NombreReel CosH (NombreReel x)
{
  NombreReel
    ValCosH,
    Xpuissance,Factoriel,i,Quotient;

  // COSH(x) = 1/2*(exp(x)+exp(-x))
  //         = 1 +x^2/2! +x^4/4! +...

  ValCosH = 1;        // Premier terme de la suite
  x = x*x;            // Eleve x au carre (pour acc‚l‚rer le calcul de
                       // Xpuisscance)
  Xpuissance = x;     // = x^2, x^4, x^6, x^8, ...
  Factoriel = 2;      // = 2!,  4!,  6!,  8!,  ...
  i = 3;              // Iteration = 3, 5, 7, 9, ...
  Quotient = Xpuissance/Factoriel;    // Calcule le quotient

  do {
    ValCosH = ValCosH +Quotient;    // Ajoute un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    Factoriel = Factoriel*i*(i+1);
    i = i +2;                       // Passe au factoriel suivant

    Quotient = Xpuissance/Factoriel;    // Calcule le quotient
  } while (PuissancePrecision<Abs(Quotient));

  return ValCosH;
}

NombreReel SinH (NombreReel x)
{
  NombreReel
    ValSinH,
    Xpuissance,Factoriel,i,Quotient;

  // SINH(x) = 1/2*(exp(x)-exp(-x))
  //         = x +x^3/3! +x^5/5! +...

  ValSinH = x;        // Premier terme de la suite
  Xpuissance = x*x*x; // = x^3, x^5, x^7, x^9, ...
  x = x*x;            // Eleve x au carre (pour acc‚l‚rer le calcul de
                       // Xpuisscance)
  Factoriel = 6;      // = 3!,  5!,  7!,  9!,  ...
  i = 4;              // Iteration = 3, 5, 7, 9, ...
  Quotient = Xpuissance/Factoriel;    // Calcule le quotient

  do {
    ValSinH = ValSinH +Quotient;    // Ajoute un terme
    Xpuissance = Xpuissance*x;      // Passe 
 la puissance suivante (* x^2)
    Factoriel = Factoriel*i*(i+1);
    i = i +2;                       // Passe au factoriel suivant

    Quotient = Xpuissance/Factoriel;    // Calcule le quotient
  } while (PuissancePrecision<Abs(Quotient));

  return ValSinH;
}

NombreReel TanH (NombreReel x)
{
  // TanH(x) = SinH(x)/CosH(x)
  // Désolé, j'ai pas trouvé mieux pour l'instant
  return SinH(x)/CosH(x);
}


//=============== ARCCOSH, ARCSINH, ARCTANH ====================================
NombreReel ArcCosH (NombreReel x)
{
  // ArcCosH(x) = Ln(x+Racine(x^2-1)), tout simplement !
  x += Racine(x*x-1);
  return Ln(x);
}

NombreReel ArcSinH (NombreReel x)
{
  NombreReel
    ValArcSinH,Rapport,XPuissance,Total;
  NombreEntier
    Impair,Pair;

  // ArcSinH (x) =
  //  x
  // -x^3/3* (1/2)
  // +x^5/5* (1*3)/(2*4)
  // -x^7/7* (1*3*5)/(2*4*6)
  // +...

  XPuissance = x*x*x;
  ValArcSinH = x-XPuissance/6;
  Rapport = (NombreReel)(1*3)/(2*4); // = (1*3)/(2*4), (1*3*5)/(2*4*6), (1*3*5*7)/(2*4*6*8),...
  Impair = 5;      // = 5, 7, 9, 11, ...
  Pair = 6;        // = 6, 8, 10, 12, ...
  x= x*x;          // Eleve x au carr‚
  XPuissance = XPuissance*x;       // Xpuissance = x^5, x^7, x^9, ...
  Total = XPuissance/Impair*Rapport;
  do {
    ValArcSinH = ValArcSinH +Total;
    Rapport = Rapport*Impair/Pair;
    Impair += 2;
    Pair += 2;
    XPuissance = XPuissance*x;

    ValArcSinH = ValArcSinH -XPuissance/Impair*Rapport;
    Rapport = Rapport*Impair/Pair;
    Impair += 2;
    Pair += 2;
    XPuissance = XPuissance*x;

    Total = XPuissance/Impair*Rapport;
  } while (PuissancePrecision<Abs(Total));

  return ValArcSinH;
}

NombreReel ArcTanH (NombreReel x)
{
  NombreReel
    ValArcTanH,XPuissance,Quotient;
  NombreEntier d;

  // ARCTANH(x) = 1/2* ln( (1+x)/(1-x) )
  //            = x +x^3/3 +x^5/5 +x^7/7 +...
  if ((x<=-1) || (1<=x)) {
    printf ("x doit appartenir à ]-1,1[ !!!\n");
    return 0;
  }

  ValArcTanH = x;
  XPuissance = x*x*x;
  x= x*x;          // Eleve x au carr‚
  d = 3;
  Quotient = XPuissance/d;
  do {
    ValArcTanH = ValArcTanH +Quotient;
    XPuissance = XPuissance*x;
    d += 2;
    Quotient = XPuissance/d;
  } while (PuissancePrecision<Abs(Quotient));

  return ValArcTanH;
}

//=============== ARCCOSH, ARCSINH, ARCTANH ====================================
//=============== ARCCOSH, ARCSINH, ARCTANH ====================================
int main(int argc, char* argv[])
{
  NombreReel x,e;
  printf (">>> Precision = %u (=nombre de chiffres significatifs) <<<\n\n",Precision);

  x = 1; AfficheExemple ("Exp",Exp);
  x = 2; AfficheExemple ("Exp",Exp);
  x = 2; AfficheExemple ("Ln",Ln);
  x = 3; AfficheExemple ("Ln",Ln);
  printf ("\n");
  getchar ();

  // --------------------------------------------
  x = 0.123; AfficheExemple ("Cos",Cos);
  x = 0.456; AfficheExemple ("Cos",Cos);
  x = 0.789; AfficheExemple ("Cos",Cos);
  printf ("\n");
  x = 0.123; AfficheExemple ("Sin",Sin);
  x = 0.456; AfficheExemple ("Sin",Sin);
  x = 0.789; AfficheExemple ("Sin",Sin);
  printf ("\n");
  x = 0.123; AfficheExemple ("Tan",Tan);
  x = 0.456; AfficheExemple ("Tan",Tan);
  x = 0.789; AfficheExemple ("Tan",Tan);
  printf ("\n");
  getchar ();

  // --------------------------------------------
  x = 0.123; AfficheExemple ("ArcCos",ArcCos);
  x = 0.456; AfficheExemple ("ArcCos",ArcCos);
  x = 0.789; AfficheExemple ("ArcCos",ArcCos);
  printf ("\n");
  x = 0.123; AfficheExemple ("ArcSin",ArcSin);
  x = 0.456; AfficheExemple ("ArcSin",ArcSin);
  x = 0.789; AfficheExemple ("ArcSin",ArcSin);
  printf ("\n");
  x = 0.123; AfficheExemple ("ArcTan",ArcTan);
  x = 0.456; AfficheExemple ("ArcTan",ArcTan);
  x = 0.789; AfficheExemple ("ArcTan",ArcTan);
  printf ("\n");
  getchar ();

  // --------------------------------------------
  x = 0.123; AfficheExemple ("CosH",CosH);
  x = 0.456; AfficheExemple ("CosH",CosH);
  x = 0.789; AfficheExemple ("CosH",CosH);
  printf ("\n");
  x = 0.123; AfficheExemple ("SinH",SinH);
  x = 0.456; AfficheExemple ("SinH",SinH);
  x = 0.789; AfficheExemple ("SinH",SinH);
  printf ("\n");
  x = 0.123; AfficheExemple ("TanH",TanH);
  x = 0.456; AfficheExemple ("TanH",TanH);
  x = 0.789; AfficheExemple ("TanH",TanH);
  printf ("\n");
  getchar ();

  // --------------------------------------------
  x = 0.123; AfficheExemple ("ArcCosH",ArcCosH);
  x = 0.456; AfficheExemple ("ArcCosH",ArcCosH);
  x = 0.789; AfficheExemple ("ArcCosH",ArcCosH);
  printf ("\n");
  x = 0.123; AfficheExemple ("ArcSinH",ArcSinH);
  x = 0.456; AfficheExemple ("ArcSinH",ArcSinH);
  x = 0.789; AfficheExemple ("ArcSinH",ArcSinH);
  printf ("\n");
  x = 0.123; AfficheExemple ("ArcTanH",ArcTanH);
  x = 0.456; AfficheExemple ("ArcTanH",ArcTanH);
  x = 0.789; AfficheExemple ("ArcTanH",ArcTanH);
  printf ("\n");
  getchar();

  // Racine carré
  x = 2; AfficheExemple ("Racine",Racine);
  x = 3; AfficheExemple ("Racine",Racine);
  x = 4; AfficheExemple ("Racine",Racine);
  printf ("\n");

  // Calcul de PI
  printf ("PI = %.12g (= ArcTan(1/2)+ArcTan(1/3))\n",(4*(ArcTan ((NombreReel)1/2)+ArcTan((NombreReel)1/3))));
  getchar ();
  return 0;
}