391 lines
13 KiB
C
391 lines
13 KiB
C
|
// 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
|