// méthodes issues de quartic.c
// on a laissé les informations des  programmes originaux
// l'ensemble est organisée sous forme d'une classe pour encapsuler

//********************************************************************************
//       l'écriture initiale a été modifiée -> C++
//
//       calcul :
//               - racine d'un polynome du second degré
//               - racine d'un polynome du 3ieme degré
//               - racine d'un polynome du 4ième degré avec différentes méthodes
//                      .  Descartes' algorithm
//                      .  Ferrari's algorithm
//                      .  Neumark's algorithm
//                      .  Yacoub & Fraidenraich's algorithm
//                      .  Christianson's algorithm
//********************************************************************************



#ifndef QUARTIC_H
#define QUARTIC_H

#include <stdio.h>

class Quartic

{ public:


/*   quartic.c
     version 54

     a test of relative accuracy of various quartic routines.

     use:
       quartic [-a] [-c n] [-q n] [-d n]
       (no parameters)  : a loop through 10,000 prechosen quartic coefficients
       -a               : improve cubic roots (default: no iteration)
       -c  n            : keep reading n cubic roots from standard input
                          (if n=0 read coefficients), terminate with 'q'.
       -d  n            : set debug value to 'n' (default: 5)
       -q  n            : keep reading n quartic roots from standard input
                          (if n=0 read coefficients), terminate with 'q'.

     debug <= 1 : print cases used

     method-
      http://linus.socs.uts.edu.au/~don/pubs/solving.html

      4 Apr 2004 bringing yacfraid notation into line with solving.html
      4 Apr 2004 fixed bug in final root calculation in yacfraid
      8 Mar 2004 corrected modified yacfraid algorithm
      8 Mar 2004 modified yacfraid algorithm
     31 Jan 2004 printing results when number of roots vary 
     30 Jan 2004 printing error of chosen cubic if debug<1
     27 Jan 2004 seeking worst coefficients for each algorithm
     27 Jan 2004 choosing best route for k in chris
     26 Jan 2004 conforming variables to solving.html
     26 Jan 2004 fixed bug in yacfraid for multiplicity 3
     25 Jan 2004 speeding up chris
     24 Jan 2004 fixed chris error (A <-> B)
     23 Jan 2004 use debug<1 for diagnostics
     21 Jan 2004 solve cubic in neumark, yacfraid and chris if d=0
     20 Jan 2004 3 roots to cubic in chris if e1=0
     19 Jan 2004 debugging Christianson routine
     14 Jan 2004 adding Christianson, cut determinant in quadratic call
      5 Jan 2004 improving cubic roots by Newton-Raphson iteration
     23 Dec 2003 seeking snag in yacfraid
     21 Dec 2003 putting diagnostic printout in yacfraid
     18 Dec 2003 allowing input of equation coefficients
     18 Dec 2003 recording cubic root giving least worst error
     17 Dec 2003 recording cubic root giving the most quartic roots
     16 Dec 2003 using the cubic root giving the most quartic roots
     16 Dec 2003 trying consistency of 3 cubic roots where available
     15 Dec 2003 trying all 3 cubic roots where available
     13 Dec 2003 trying to fix Neumark
     13 Dec 2003 cleaning up diagnostic format
     12 Dec 2003 initialising n,m,po3 in cubic
     12 Dec 2003 allow cubic to return 3 zero roots if p=q=r=0
     10 Dec 2003 added setargs
      2 Dec 2003 finding worst cases
      2 Dec 2003 negating j if p>0 in cubic
      1 Dec 2003 changing v in cubic from (sinsqk > doub0) to (sinsqk >= doub0)
      1 Dec 2003 test changing v in cubic from po3sq+po3sq to doub2*po3sq
     30 Nov 2003 counting cases in all solving routines
     29 Nov 2003 testing wsq >= doub0
     29 Nov 2003 mult by doub2 for v in cubic
     29 Nov 2003 cut po3cu from cubic
     29 Nov 2003 better quadratic
     29 Nov 2003 count agreements
     17 Nov 2003 count special cases
     16 Nov 2003 option of loop or read coefficients from input
     15 Nov 2003 fixed cubic() bug
     11 Nov 2003 added Brown and Yacoub & Fraidenraich's algorithms
     21 Jan 1989 quartic selecting Descartes, Ferrari, or Neumark algorithms 
     16 Jul 1981 Don Herbison-Evans

"Solving Quartics and Cubics for Graphics", D. Herbison-Evans,
Graphics Gems V (ed.: A. Paeth) Academic Press (Chesnut Hill), pp. 3-15 (1995). 

"Solving Quartics and Cubics for Graphics", D. Herbison-Evans,
Research Report CS-86-56, Department of Computer Science, University of Waterloo (1986)
 
"Caterpillars and the Inaccurate Solution of Cubic and Quartic Equations", 
D. Herbison-Evans, Australian Computer Science Communications, 
Vol. 5, No. 1, pp. 80-90 (1983)

   subroutines:
      errors      - find errors in a set of roots
      acos3       - find arcos(x/3)
      curoot      - find cube root
      quadratic   - solve a quadratic
      cubic       - solve a cubic
      cubnewton   - improve cubic roots by iteration
      quartic     - solve a quartic
      descartes   - use Descartes' algorithm to solve a quartic  
      ferrari     - use Ferrari's algorithm to solve a quartic 
      neumark     - use Neumark's algorithm to solve a quartic 
      yacfraid    - use Yacoub & Fraidenraich's algorithm to solve a quartic 
      chris       - use Christianson's algorithm to solve a quartic
*/


/***********************************/

// constructeur permettant d'initialiser les variables
Quartic();

// modification du paramètre de débug et d'itération pour cubic
void Change(int debuge, bool iteratee) {debug=debuge;iterate=iteratee;};

//****************************************************
//
//   find the errors
//
//   called by quartictest, docoeff, compare,
//            chris, descartes, ferrari, neumark, yacfraid.
//****************************************************

double errors(double a,double b,double c,double d,double rts[4],double rterr[4],int nrts);
//double a,b,c,d,rts[4],rterr[4];
//int nrts;

/**********************************************/
//     find cos(acos(x)/3) 
//    
//     16 Jul 1981   Don Herbison-Evans 
//
//     called by cubic . 
/**********************************************/

double acos3(double x);
//   double x ;

/***************************************/
//     find cube root of x.
//
//     30 Jan 1989   Don Herbison-Evans 
//
//     called by cubic . 
/***************************************/

double curoot(double x);
//   double x ;

/****************************************************/
//     solve the quadratic equation - 
//
//         x**2 + b*x + c = 0 
//
//     14 Jan 2004   cut determinant in quadratic call
//     29 Nov 2003   improved
//     16 Jul 1981   Don Herbison-Evans
//
//     called by  cubic,quartic,chris,descartes,ferrari,neumark.
/****************************************************/

int quadratic(double b,double c,double rts[4]);
//   double b,c,rts[4];


/**************************************************/
//   find the real roots of the cubic - 
//       x**3 + p*x**2 + q*x + r = 0 
//
//     12 Dec 2003 initialising n,m,po3
//     12 Dec 2003 allow return of 3 zero roots if p=q=r=0
//      2 Dec 2003 negating j if p>0
//      1 Dec 2003 changing v from (sinsqk > doub0) to (sinsqk >= doub0)
//      1 Dec 2003 test changing v from po3sq+po3sq to doub2*po3sq
//     16 Jul 1981 Don Herbison-Evans
//
//   input parameters - 
//     p,q,r - coeffs of cubic equation. 
//     iterate : booléen indiquant si l'on veut améliorer les racines via du newton (4 fois)
//     debug : si <1 impression d'un paquet de paramètre ??
//
//   output- 
//     the number of real roots
//     v3 - the roots. 
//
//   global constants -
//     rt3 - sqrt(3) 
//     inv3 - 1/3 
//     doubmax - square root of largest number held by machine 
//
//     method - 
//     see D.E. Littlewood, "A University Algebra" pp.173 - 6 
//
//     15 Nov 2003 output 3 real roots: Don Herbison-Evans
//        Apr 1981 initial version: Charles Prineas
//
//     called by  cubictest,quartic,chris,yacfraid,neumark,descartes,ferrari.
//     calls      quadratic,acos3,curoot,cubnewton. 
/**************************************************/

int cubic(double p,double q,double r, double v3[4]) ;
 


/****************************************************/
//    improve roots of cubic by Newton-Raphson iteration
//
//    5 Jan 2004  Don Herbison-Evans
//
//    called by cubic.
/****************************************************/

void cubnewton(double p,double q,double r,int n3,double v3[4]);
//int n3;
//double p,q,r,v3[4];
//

//****************************************************
//   Solve quartic equation using either
//   quadratic, Ferrari's or Neumark's algorithm.
//
//   called by 
//   calls  descartes, ferrari, neumark, yacfraid.
//
//     15 Dec 2003  added yacfraid
//     10 Dec 2003  added descartes with neg coeffs
//     21 Jan 1989  Don Herbison-Evans
//****************************************************
int quartic(double a,double b,double c,double d,double rts[4]);
//double a,b,c,d,rts[4];

//*****************************************/
//   Solve quartic equation using
//   Descartes-Euler-Cardano algorithm
//
//   called by quartic, compare, quartictest.
//
//   Strong, T. "Elemementary and Higher Algebra"
//      Pratt and Oakley, p. 469 (1859)
//
//   16 Jul 1981  Don Herbison-Evans
//*****************************************/
int descartes(double a,double b,double c,double d,double rts[4]);
//double a,b,c,d,rts[4];


//****************************************************/
//     solve the quartic equation - 
//
//     x**4 + a*x**3 + b*x**2 + c*x + d = 0 
//
//     calls     cubic, quadratic.
//
//     input - 
//     a,b,c,e - coeffs of equation. 
//
//     output - 
//     n4 - number of real roots. 
//     rts - array of root values. 
//
//     method :  Ferrari - Lagrange
//     Theory of Equations, H.W. Turnbull p. 140 (1947)
//
//     16 Jul 1981 Don Herbison-Evans
//
//     calls  cubic, quadratic 
//****************************************************/
int ferrari(double a,double b,double c,double d,double rts[4]);
//   double a,b,c,d,rts[4];
 


//*****************************************/
//     solve the quartic equation - 
//
//     x**4 + a*x**3 + b*x**2 + c*x + d = 0 
//
//    calls     cubic, quadratic.
//
//    input parameters - 
//     a,b,c,e - coeffs of equation. 
//
//     output parameters - 
//     n4 - number of real roots. 
//     rts - array of root values. 
//
//     method -  S. Neumark 
//        "Solution of Cubic and Quartic Equations" - Pergamon 1965 
//
//      1 Dec 1985   translated to C with help of Shawn Neely
//     16 Jul 1981   Don Herbison-Evans
//*****************************************/
int neumark(double a,double b,double c,double d,double rts[4]);
//   double a,b,c,d,rts[4];
 


//****************************************************/ 
//     solve the quartic equation - 
//     x**4 + a*x**3 + b*x**2 + c*x + d = 0 
//
//     calls     cubic, quadratic.
//
//     input parameters - 
//     a,b,c,e - coeffs of equation. 
//
//     output parameters - 
//     n4 - number of real roots. 
//     rts - array of root values. 
//
//     method - 
//        K.S. Brown 
//        Reducing Quartics to Cubics,
//        http://www.seanet.com/~ksbrown/kmath296.htm (1967)
// 
//        Michael Daoud Yacoub & Gustavo Fraidenraich
//        "A new simple solution of the general quartic equation"
//        Revised 16 Feb 2004
//
//     14 Nov 2003 Don Herbison-Evans
//*****************************************/
int yacfraid(double a,double b,double c,double d,double rts[4]);
//   double a,b,c,d,rts[4];


//*****************************************/
//   Solve quartic equation using
//   Christianson's algorithm
//   calls errors, quadratic, cubic.
//
//   Bruce Christianson, Solving Quartics Using Palindromes,
//      Mathematical Gazette, Vol. 75, pp. 327-328 (1991)
//
//   14 Jan 2004  Don Herbison-Evans

int chris(double a,double b,double c,double d,double rts[4]);
//double a,b,c,d,rts[4];
 
 // données propres
 protected:
 #define NCASES 40

     double d3o8,d3o256;       /* double precision constants */
     double doub0;
     double doub1,doub2;
     double doub3,doub4;
     double doub5,doub6;
     double doub8,doub9,doub12;
     double doub16,doub24;
     double doub27,doub64;
     double doubmax;           /* approx square root of max double number */
     double doubmin;           /* smallest double number */
     double doubtol;           /* tolerance of double numbers */
     double inv2,inv3,inv4;
     double inv8,inv16,inv32,inv64,inv128;
     double qrts[4][3];        /* quartic roots for each cubic root */
     double rt3;               /* square root of 3 */
     double rterc[4],rterd[4],rterf[4],rtern[4],rterq[4],rtery[4]; /* errors of roots */
     double worst3[3];         /* worst error for a given cubic root */
     int debug;                /* <1 for lots of diagnostics, >5 for none */   
     bool iterate; // indique si l'on veut itérer on pas 
     int j3;
     int n1,n2,n3,n4[3];
     int nqud[NCASES];
     int ncub[NCASES];
     int nchr[NCASES];
     int ndes[NCASES];
     int nfer[NCASES];
     int nneu[NCASES];
     int nyac[NCASES];
};

 
#endif  