/*******************************************************
*             Equations of degree 2 to 4               *
* ---------------------------------------------------- *
* This program calculates the real or complex roots of *
* algebraic equations of degree 2, 3 and 4.            *
* ---------------------------------------------------- *
* Reference: Mathématiques et statitiques - Programmes *
*            en BASIC. Editions du P.S.I., 1981        *
*            [BIBLI 13].                               *
*                                                      *
*                    C++ version by J-P Moreau, Paris  *
*                          (www.jpmoreau.fr)           *
* ---------------------------------------------------- *
* SAMPLE RUN:                                          *
*                                                      *
* ROOTS OF EQUATIONS OF DEGREE 2, 3 OU 4               *
* Example: X^4+4X^2-3X+3 = 0                           *
*                                                      *
* INSTRUCTIONS:                                        *
*                                                      *
* 1. Input degree of equation                          *
* 2. Input coefficients beginning with highest exponent*
* 3. Solutions are given under the form:               *
*     (real part) + I (imaginary part)                 *
*                                                      *
* Degree of equation (2 to 4): 4                       *
*                                                      *
* Coefficient of X^4 = 1                               *
* Coefficient of X^3 = 0                               *
* Coefficient of X^2 = 4                               *
* Coefficient of X^1 = -3                              *
* Coefficient of X^0 = 3                               *
*                                                      *
* Root 1 = 0.449702657 + I 0.731070322                 *
* Root 2 = 0.449702657 - I 0.731070322                 *
* Root 1 = -0.449702657 + I 1.967231848                *
* Root 2 = -0.449702657 - I 1.967231848                *
*                                                      *
*******************************************************/
#include <stdio.h>
#include <math.h>

double a[5][5], r[5];
int    i,n;
double aa,b,c,d,ii,im,k,l,m,q,rr,s,sw;


double f(double x) {
  return x * x * x + a[3][2] * x * x + a[3][1] * x + a[3][0];
}


/******************************************
* This Subroutine calculates the roots of *
*    X^2 + A(2,1)*X + A(2,0) = 0          *
*    W = Determinant of equation          *
******************************************/
void Root_2()  {
  double q1, q2, w;
  w = a[2][1] * a[2][1] - 4 * a[2][0];
  q1 = -a[2][1] / 2.0;
  if (w < 0) {
    q2 = sqrt(-w) / 2.0;
    im = q2;
    r[1] = q1; r[2] = q1;
  }
  else if (w==0) {
    r[1] = q1; r[2] = q1; im = 0.0;
  }
  else if (w > 0) {
    q2 = sqrt(w) / 2.0; im = 0.0;
    r[1] = q1 + q2; r[2] = q1 - q2;
  }
}


/*******************************************
* This subroutine calculates the roots of  *
* X^3 + A(3,2)*X^2 + A(3,1)*X + A(3,0) = 0 *
*******************************************/
void Root_3()  {
  //Labels e100,e200,e500
  double am, er, te, tt, xa, xb, xc, y1, y2;
  // one root equals zero
  if (a[3][0] == 0)  {
    xc = 0; goto e500;
  }
  // looking for( maximum coefficient in fabsolute magnitude
  am = fabs(a[3][0]);
  for (i = 1; i<4; i++) {
    tt = fabs(a[3][i]);
    if (am < tt) am = tt;
  }
  //Define interval where a real root exists
  // according to sign of A(3,0)  
  if (a[3][0] > 0) {
    aa = -am - 1;
    b = 0;
    goto e100;
  }
  aa = 0; b = am + 1;

e100:
  // Searching for( xc = real root in interval (a,b)
  //   (by Bisection method)
  
  // Define segment (xa,xb) including the root
  xa = aa; xb = b; y1 = f(aa); y2 = f(b);

  // Desired précision
  er = 0.000001;
  if (y1 == 0) {
    xc = aa; goto e500;
  }
  if (y2 == 0) {
    xc = b; goto e500;
  }
  // xc : middle of segment
  e200: xc = (xa + xb) / 2.0; te = f(xc);
  if (te == 0) goto e500;
  if ((xb - xa) < er) goto e500;
  if (f(xa) * te > 0) {
    xa = xc; goto e200;
  }
  xb = xc; goto e200;
  // r[3] is a real root
e500: r[3] = xc;
  if (sw = -1) return;
  // Calculates the roots of remaining quadratic equation
  // Define the equation coefficients
  a[2][1] = a[3][2] + xc;
  a[2][0] = a[3][1] + a[3][2] * xc + xc * xc;
  Root_2();
  return;
}

// Root search main subroutine
void Root_4()  {
  // Labels: e100,e200,e300
  // Normalization of coefficients
  for (i = 0; i<n; i++) {
    a[n][i] /= a[n][n];
  }
  // Branching according to degree n
  if (n==2) {
    Root_2();
    return;
  }
  else if (n==3) {
    Root_3();
    return;
  }
  else if (n==4) {
    aa = a[4][3];
    b = a[4][2];
    c = a[4][1];
    d = a[4][0];
    q = b - (3.0 * aa * aa / 8.0);
    rr = c - (aa * b / 2) + (aa * aa * aa / 8.0);
    s = d - (aa * c / 4) + (aa * aa * b / 16.0) - (3 * aa * aa * aa * aa / 256.0);
    // Défine coefficients of cubic equation
    a[3][2] = q / 2.0;
    a[3][1] = (q * q - 4 * s) / 16.0;
    a[3][0] = -(rr * rr / 64.0);
    // Calculate a real root of this equation
    if ((rr != 0) || (a[3, 1] >= 0)) goto e100;
    // Particular case when this equation is of 2nd order
    a[2][1] = a[3][2]; a[2][0] = a[3][1];
    Root_2();
    // One takes the positive root
    r[3] = r[1];
    goto e200;
    // Calling Root_3 with sw=-1 to calculate one root only
e100: sw = -1;
    Root_3();
    // real root of above equation
e200: k = sqrt(r[3]);
    // Calculate L and M if k=0
    if (k == 0) {
      rr = sqrt(q * q - 4 * s); goto e300;
    }
    q = q + (4 * r[3]);
    rr = rr / (2 * k);
e300: l = (q - rr) / 2.0; m = (q + rr) / 2.0;
    // Solving two equations of degree 2
    a[2][1] = 2 * k; a[2][0] = l;
    // 1st equation
    Root_2();
    // Transferring solutions in r[3], r[4], ii
    r[3] = r[1] - (a[4][3] / 4); r[4] = r[2] - (a[4][3] / 4.0); ii = im;
    a[2][1] = (double) -2 * k; a[2][0] = m;
    // 2nd equation
    Root_2();
    r[2] = r[2] - (a[4][3] / 4.0); r[1] = r[1] - (a[4][3] / 4.0);
  }
}


void main()  {
  printf("\n ROOTS OF EQUATIONS OF DEGREE 2, 3 OR 4\n");
  printf(" EXAMPLE : X^4+4X^2-3X+3 = 0\n\n");
  printf(" INSTRUCTIONS:\n\n");
  printf(" 1. Input degree of equation\n");
  printf(" 2. Input coefficients beginning by highest exponent\n");
  printf(" 3. Solutions are given under the form:\n");
  printf("    (real part) + I (imaginary part)\n\n");

  printf(" DEGREE OF EQUATION [2 to 4]: "); scanf("%d",&n);
  if (n < 2) n = 2;
  if (n > 4) n = 4;
  printf("\n");
  for (i = n; i>-1; i--) {
    printf(" COEFFICIENT OF X^%d = ",i); scanf("%lf",&a[n][i]);
  }
  printf("\n");
  // calling root search main subroutine
  Root_4();
  // writing results
  if (n==2)  {
    //case n=2
    if (im == 0) {
      printf(" Root 1 = %12.9f\n",r[1]);
      printf(" Root 2 = %12.9f\n",r[2]);
    }
    else {
      printf(" Root 1 = %12.9f +I %12.9f\n",r[1],im);
      printf(" Root 2 = %12.9f -I %12.9f\n",r[2],im);
    }
  }
  else if (n==3)  {
    //case n=3
    printf(" Root 1 = %12.9f\n",r[3]);
    if (im == 0) { 
      printf(" Root 2 = %12.9f\n",r[1]);
      printf(" Root 3 = %12.9f\n",r[2]);
    }
    else {
      printf(" Root 2 = %12.9f +I %12.9f\n",r[1],im);
      printf(" Root 3 = %12.9f -I %12.9f\n",r[2],im);
    }
  }
  else if (n==4) {
    //case n=4
    if (im==0) {
      printf(" Root 1 = %12.9f\n",r[1]);
      printf(" Root 2 = %12.9f\n",r[2]);
    }
    else {
      printf(" Root 1 = %12.9f +I %12.9f\n",r[1],im);
      printf(" Root 2 = %12.9f -I %12.9f\n",r[2],im);
    }
    if (ii==0) {
      printf(" Root 3 = %12.9f\n",r[3]);
      printf(" Root 4 = %12.9f\n",r[4]);
    }
    else {
      printf(" Root 3 = %12.9f +I %12.9f\n",r[3],ii);
      printf(" Root 4 = %12.9f -I %12.9f\n",r[4],ii);
    }
  }
  printf("\n");
}

// End of file root4.cpp