/**************************************************************
*            Solve complex equation of degree 4:              *
*            a z^4 + b z^3 + c z^2 + d z + e = 0              *
* ----------------------------------------------------------- *
* SAMPLE RUN:                                                 *
*    a = (1,0) b = (2,-10) c = (-4,1) d = (5,2) e = (3,-7.5)  *
*                                                             *
* Equation az^4 + bz^3 + cz^2 + d z + e = 0                   *
* Solutions z1..z4:                                           *
*                                                             *
* z1 = (-2.139849,9.625142)                                   *
* F(z1) = (0.000000,0.000000)                                 *
* z2 = (0.690192,-0.706380)                                   *
* F(z2) = (0.000000,-0.000000)                                *
* z3 = (-0.875877,0.223149)                                   *
* F(z3) = (0.000000,-0.000000)                                *
* z4 = (0.325534,0.858088)                                    *
* F(z4) = (0.000000,-0.000000)                                *
*                                                             *
* ----------------------------------------------------------- *
* Reference: Mathématiques et statitiques - Programmes en     *
*            BASIC. Editions du P.S.I., 1981.                 *
*                                                             *
*            Adapted to complex Domain By J-P Moreau, Paris.  *
*                            (www.jpmoreau.fr)                *
***************************************************************
 Explanations:
 ------------
 Let us solve complex equation z^4 + a*z^3 + b*z^2 + c*z + d = 0  (1)

 Let us set z = z1 - a/4 and equation (1) becomes:

     (z1²+2ckz1+l)(z1²-2ckz1+m) = 0       (2)

 where:   l + m -4ck² = q  |            q = 3a²/8
          2 ck (m-1)  = r  | (3) with:  r = c-(ab/2)+(a^3/8)
                 lm   = s  |            s = d-(ac/4)+(a²b/16)-(3a^4/256)

 ck² = zz is a solution of complex cubic equation:

          zz^3 + aa zz^2 + bb z + cc = 0   (4)

 with:    aa = q/2, bb = (q²-4s)/16, cc = -r²/64

 So we solve equation (4) using Croot3 (3 roots zz1, zz2, zz3), then
 we calculate complex l and m by solving system (3), finally, we solve
 both 2nd degree equations of product (2). We obtain three sets of four
 complex roots, only one of which is valid, i.e. annulates equation (1). 

----------------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

#define  PI  4*atan(1.0)   // used by complex cubic root

typedef struct {
          double x, y;
		} COMPLEX;

// return ABS(z)
double CABS(COMPLEX z) {
  return (sqrt(z.x*z.x + z.y*z.y));
}

// z = z1 * z2
void CMUL(COMPLEX z1, COMPLEX z2, COMPLEX *z) {
  z->x = z1.x*z2.x - z1.y*z2.y;
  z->y = z1.x*z2.y + z1.y*z2.x;
}

// z = z1 / z2
void CDIV(COMPLEX z1, COMPLEX z2, COMPLEX *z) {
  double d; COMPLEX c;
  d = z2.x*z2.x + z2.y*z2.y;
  if (d<1E-15)
    printf(" Complex Divide by zero!\n");
  else {
    c.x=z2.x; c.y=-z2.y;
    CMUL(z1,c,z);
    z->x=z->x/d; z->y=z->y/d;
  }
}

// z1 = SQRT(z)
void CSQRT(COMPLEX z, COMPLEX *z1) {
  double r;
  r=sqrt(z.x*z.x+z.y*z.y);
  z1->x = sqrt((r+z.x)/2);
  z1->y = sqrt((r-z.x)/2);
  if (z.y<0)  z1->y=-z1->y;
}

double ATAN(double numerator, double denominator) {
// Return a phase between -PI and +PI
  double phase;
  if (fabs(denominator) < 1e-15) {
    if (numerator < 0) return (-PI/2);
    else return (PI/2);
  }
  else {
    phase=atan(numerator/denominator);
    if (denominator < 0)  phase += PI;
	return phase;
  }
} //ATAN()

// z1, z2, z3 = z^(1/3)
void CUBRT(COMPLEX z, COMPLEX *z1, COMPLEX *z2, COMPLEX *z3) {
  double r,t,tt;
  r=sqrt(z.x*z.x+z.y*z.y);
  r=pow(r,1.0/3.0);
  t = ATAN(z.y,z.x);
  t=t/3.0;
  z1->x=r*cos(t);
  z1->y=r*sin(t);
  tt=t-(2*PI/3);
  z2->x=r*cos(tt);
  z2->y=r*sin(tt);
  tt=t+(2*PI/3);
  z3->x=r*cos(tt);
  z3->y=r*sin(tt);
}

// global variables
COMPLEX  a,b,c,d,e,z1,z2,z3,z4,zz;

// return z1 = z^3 + p*z + q (used by Croot3)
void F3(COMPLEX p,COMPLEX q,COMPLEX z, COMPLEX *z1) {
  COMPLEX u,v,w;
  CMUL(z,z,&u); CMUL(z,u,&v); 
  CMUL(p,z,&w);
  z1->x=v.x + w.x + q.x;
  z1->y=v.y + w.y + q.y;
}

// return a*z^4+b*z^3+c*z^2+d^z+e (used by main - optional)
  void F4(COMPLEX a,COMPLEX b,COMPLEX c,COMPLEX d, COMPLEX e,
	      COMPLEX z, COMPLEX *z1)  {
  COMPLEX u,v,w;
    CMUL(z,z,&u); CMUL(u,u,&v); CMUL(a,v,z1);
    CMUL(u,z,&v); CMUL(b,v,&w);
    z1->x=z1->x+w.x; z1->y=z1->y+w.y;
    CMUL(c,u,&v);
    z1->x=z1->x+v.x; z1->y=z1->y+v.y;
    CMUL(d,z,&u);
    z1->x=z1->x+u.x+e.x; z1->y=z1->y+u.y+e.y;
  }

// return complex roots of equation a*z^2+b*z+c = 0 (used by Croot4)
void Equa2c(COMPLEX a,COMPLEX b,COMPLEX c, COMPLEX *s1, COMPLEX *s2) {
  COMPLEX u,v,w; double r;
  u.x = b.x * b.x - b.y * b.y - 4 * a.x * c.x + 4 * a.y * c.y;
  u.y = 2 * b.x * b.y - 4 * a.x * c.y - 4 * a.y * c.x;
  r = sqrt(u.x*u.x + u.y*u.y);
  v.x = sqrt((r+u.x)/2);
  v.y = sqrt((r-u.x)/2);
  if (u.y<0)  v.y=-v.y;
  w.x = (-b.x - v.x)/2;
  w.y = (-b.y - v.y)/2;
  u.x = (-b.x + v.x)/2;
  u.y = (-b.y + v.y)/2;
  r = a.x * a.x + a.y * a.y;
  s1->x = (a.x * w.x + a.y * w.y)/r;
  s1->y = (a.x * w.y - a.y * w.x)/r;
  s2->x = (a.x * u.x + a.y * u.y)/r;
  s2->y = (a.x * u.y - a.y * u.x)/r;
}

// solve complex equation a*z^3+b*z^2+c*z+d=0 and return
// 3 solutions z1,z2,z3. (used by Croot4)
void Croot3(COMPLEX a,COMPLEX b,COMPLEX c,COMPLEX d, COMPLEX *z1, 
			COMPLEX *z2, COMPLEX *z3) {
/*-----------------------------------------------------------------------------------
   z = Z - (b/3a
   p = (c/a) - (b²/3a²)
   q = (2b^3/27a^3) + (d/a) - (bc)/3a²)

  Now the equation is reduced to Z^3 + pZ + q = 0

  with the roots:

  Z[1..3] = cubroot(-q/2 + sqrt(q²/4 + p^3/27)) + cubroot(-q/2 - sqrt(q²/4 + p^3/27))
-------------------------------------------------------------------------------------*/
  COMPLEX p,q,det,sdet,u,u1,u2,v,v1,v2,w,zz;
  COMPLEX z[10]; int i,index;

// calculate p
  CDIV(c,a,&u); CMUL(b,b,&v); CMUL(a,a,&w);
  w.x=3.0*w.x; w.y=3.0*w.y; CDIV(v,w,&zz);
  p.x = u.x - zz.x; p.y = u.y - zz.y;

// calculate q
  CMUL(b,v,&w); w.x=2.0*w.x; w.y=2.0*w.y;
  CMUL(a,a,&u); CMUL(u,a,&v);
  v.x=27.0*v.x; v.y=27.0*v.y;
  CDIV(w,v,&q);
  CDIV(d,a,&u);
  q.x=q.x+u.x; q.y=q.y+u.y;
  CMUL(b,c,&u); CMUL(a,a,&v);
  v.x=3.0*v.x; v.y=3.0*v.y;
  CDIV(u,v,&w);
  q.x=q.x-w.x; q.y=q.y-w.y;

  CMUL(q,q,&u); u.x=u.x/4.0; u.y=u.y/4.0;
  CMUL(p,p,&v); CMUL(p,v,&w); w.x=w.x/27.0; w.y=w.y/27.0;
  det.x=u.x + w.x; det.y=u.y + w.y;

  CSQRT(det,&sdet);  // now sdet contains sqrt(q²/4 + p^3/27))

  v.x = -q.x/2.0 + sdet.x; v.y = -q.y/2.0 + sdet.y;
  CUBRT(v,&u,&u1,&u2);   // 3 cubic roots

  w.x = -q.x/2.0 - sdet.x; w.y = -q.y/2.0 - sdet.y;
  CUBRT(w,&v,&v1,&v2);   // 3 cubic roots

  z[1].x = u.x + v.x; z[1].y = u.y + v.y;
  z[2].x = u.x + v1.x; z[2].y = u.y + v1.y;
  z[3].x = u.x + v2.x; z[3].y = u.y + v2.y;
  z[4].x = u1.x + v.x; z[4].y = u1.y + v.y;
  z[5].x = u1.x + v1.x; z[5].y = u1.y + v1.y;
  z[6].x = u1.x + v2.x; z[6].y = u1.y + v2.y;
  z[7].x = u2.x + v.x;  z[7].y = u2.y + v.y;
  z[8].x = u2.x + v1.x; z[8].y = u2.y + v1.y;
  z[9].x = u2.x + v2.x; z[9].y = u2.y + v2.y;

// Note: only 3 z[i] are correct solutions.

  index=1;
  for (i=1; i<10; i++) {
    F3(p, q, z[i], &zz);
    if (CABS(zz) < 0.001) { //detect if solution is correct
	  if (index==1) {
        z1->x=z[i].x; z1->y=z[i].y; index++;
      }
      else if (index==2) {
        z2->x=z[i].x; z2->y=z[i].y; index++;      }
      else {
        z3->x=z[i].x; z3->y=z[i].y;
      }
    }
  }

  u.x=3.0*a.x; u.y=3.0*a.y;
  CDIV(b,u,&v);

  z1->x = z1->x - v.x; z1->y = z1->y - v.y;
  z2->x = z2->x - v.x; z2->y = z2->y - v.y;
  z3->x = z3->x - v.x; z3->y = z3->y - v.y;

} // croot3


  // return z^4+a*z^3+b*z^2+c^z+d (used by Croot4)
  void F(COMPLEX a,COMPLEX b,COMPLEX c,COMPLEX d, COMPLEX z,
	     COMPLEX *z1)  {
    COMPLEX u,v,w;
	CMUL(z,z,&u); CMUL(u,u,&v); z1->x=v.x; z1->y=v.y;
    CMUL(u,z,&v); CMUL(a,v,&w);
    z1->x=z1->x+w.x; z1->y=z1->y+w.y;
    CMUL(b,u,&v);
    z1->x=z1->x+v.x; z1->y=z1->y+v.y;
    CMUL(c,z,&u);
    z1->x=z1->x+u.x+d.x; z1->y=z1->y+u.y+d.y;
  }

/**************************************************************
* solve complex equation a*z^4+b*z^3+c*z^2+d*z+e=0 and return *
* 4 solutions z1,z2,z3,z4).                                   *
* ----------------------------------------------------------- *
* INPUTS:                                                     *
*         a,b,c,d,e:: complex coefficients.                   *
* OUTPUTS:                                                    *
*         z1,z2,z3,z4: complex roots of equation.             * 
**************************************************************/
void Croot4(COMPLEX a,COMPLEX b,COMPLEX c,COMPLEX d,COMPLEX e, 
			COMPLEX *z1, COMPLEX *z2, COMPLEX *z3, COMPLEX *z4) {
//Labels: e100,e200,e300,e400;
  COMPLEX a1,aa,b1,bb,c1,cc, q,r,s,t,u,v,w;
  COMPLEX ck,l,m;
  COMPLEX z[4];
  COMPLEX zz[13];
  int i;

  t.x=a.x; t.y=a.y;
  CDIV(b,t,&a); CDIV(c,t,&b); CDIV(d,t,&c); CDIV(e,t,&d);
// Now equation is: z^4+a*z^3+b*z^2+c*z+d = 0

// q = b - (3.0 * a * a / 8.0)
  CMUL(a,a,&u); u.x=u.x*3.0/8.0; u.y=u.y*3.0/8.0;
  q.x=b.x - u.x; q.y=b.y - u.y;

// r = c - (a * b / 2) + (a * a * a / 8.0)
  CMUL(a,b,&u);  r.x=-u.x/2.0; r.y=-u.y/2.0;
  CMUL(a,a,&u); CMUL(a,u,&v); v.x=v.x/8.0; v.y=v.y/8.0;
  r.x=r.x+c.x + v.x; r.y=r.y+c.y + v.y;

// s = d - (a * c / 4) + (a * a * b / 16.0) - (3 * a * a * a * a / 256.0)
  CMUL(a,c,&u); u.x=u.x/4.0; u.y=u.y/4.0;
  CMUL(a,a,&v); CMUL(b,v,&w); v.x=w.x/16.0; v.y=w.y/16.0;
  CMUL(a,a,&t); CMUL(t,t,&w); w.x=w.x*3.0/256.0; w.y=w.y*3.0/256.0;
  s.x=d.x - u.x + v.x -w.x; s.y=d.y - u.y + v.y -w.y;

// Défine coefficients of cubic equation z^3+aa*z^2+bb*z+cc=0  (1)
// aa = q / 2.0
  aa.x=q.x/2.0; aa.y=q.y/2.0;
// bb = (q * q - 4 * s) / 16.0
  CMUL(q,q,&u);
  bb.x=(u.x-4.0*s.x)/16.0; bb.y=(u.y-4.0*s.y)/16.0;
// cc = -(r * r / 64.0)
  CMUL(r,r,&u);
  cc.x=-u.x/64.0; cc.y=-u.y/64.0;

// Calculate complex roots equation (1)
  a1.x=1.0; a1.y=0.0;
  if (CABS(r) > 1e-15 || CABS(bb) > 1e-15) goto e100;

// Particular case when equation (1) is of 2nd order
  b1=aa; c1=bb;
  Equa2c(a1,b1,c1,&z[1],&z[2]);
  z[3].x=0.0; z[3].y=0.0;
  goto e200;

e100: Croot3(a1,aa,bb,cc,&z[1],&z[2],&z[3]);

e200: for (i=1; i<4; i++) {
    CSQRT(z[i], &ck);

// Calculate l and m if k=0
    if (CABS(ck) == 0) {
//    r = SQRT(q * q - 4 * s)
      CMUL(q,q,&u); v.x=4.0*s.x;  v.y=4.0*s.y;
      r.x = u.x - v.x; r.y = u.y - v.y;
      goto e300;
    }
//  q = q + (4 * z[i])
    u.x=4.0*z[i].x; u.y=4.0*z[i].y;
    q.x=q.x + u.x; q.y = q.y + u.y;
//  r = r / (2 * ck)
    u.x=2.0*ck.x; u.y=2.0*ck.y;
    CDIV(r,u,&r);
//   l = (q - r) / 2.0, m = (q + r) / 2.0
e300:l.x=(q.x-r.x)/2.0; l.y=(q.y-r.y)/2.0;
    m.x=(q.x+r.x)/2.0; m.y=(q.y+r.y)/2.0;

//  Solving two equations of degree 2
    b1.x = 2.0 * ck.x; b1.y = 2.0 * ck.y;
//  1st equation
    Equa2c(a1,b1,l,&u,&v);

    b1.x=-b1.x; b1.y=-b1.y;
//  2nd equation
    Equa2c(a1,b1,m,&w,&t);

//  Transferring solutions in zz(i)
    switch(i) {
      case 1: zz[1].x=u.x; zz[1].y=u.y;
		      zz[2].x=v.x; zz[2].y=v.y;
			  zz[3].x=w.x; zz[3].y=w.y;
			  zz[4].x=t.x; zz[4].y=t.y;
			  break;
      case 2: zz[5].x=u.x; zz[5].y=u.y;
		      zz[6].x=v.x; zz[6].y=v.y;
			  zz[7].x=w.x; zz[7].y=w.y;
			  zz[8].x=t.x; zz[8].y=t.y;
			  break;
      case 3: zz[9].x=u.x;  zz[9].y=u.y;
		      zz[10].x=v.x; zz[10].y=v.y;
			  zz[11].x=w.x; zz[11].y=w.y;
			  zz[12].x=t.x; zz[12].y=t.y;
	}  //Note: only 4 solutions are correct
  } //i loop

  u.x=a.x/4.0; u.y=a.y/4.0;
  for (i=1; i<=12; i++) { // shift zz by -a/4
    zz[i].x=zz[i].x-u.x; zz[i].y=zz[i].y-u.y;
  }
// select 4 correct solutions.
  i=1;
e400: F(a,b,c,d,zz[i], &t);
  if (CABS(t) < 0.0001)
	switch(i) {
      case 1: z1->x=zz[1].x; z1->y=zz[1].y; z2->x=zz[2].x; z2->y=zz[2].y; 
              z3->x=zz[3].x; z3->y=zz[3].y; z4->x=zz[4].x; z4->y=zz[4].y;
			  break;
      case 5: z1->x=zz[5].x; z1->y=zz[5].y; z2->x=zz[6].x; z2->y=zz[6].y; 
              z3->x=zz[7].x; z3->y=zz[7].y; z4->x=zz[8].x; z4->y=zz[8].y;
			  break;
      case 8: z1->x=zz[9].x; z1->y=zz[9].y; z2->x=zz[10].x; z2->y=zz[10].y; 
              z3->x=zz[11].x; z3->y=zz[11].y; z4->x=zz[12].x; z4->y=zz[12].y;
			  break;
  }
  i += 4;
  if (i<10) goto e400;

// now complex roots are in z1,z2,z3,z4   

}

void main() {

  printf("\n Equation az^4 + bz^3 + cz^2 + d z + e = 0\n");
  printf(" Solutions z1..z4:\n\n");

// Define complex coefficients a, b, c, d, e
  a.x =  1.0; a.y =   0.0;
  b.x =  2.0; b.y = -10.0;
  c.x = -4.0; c.y =   1.0;
  d.x =  5.0; d.y =   2.0;
  e.x =  3.0; e.y =  -7.5;

  Croot4(a,b,c,d,e, &z1,&z2,&z3,&z4);

  printf(" z1 = (%f,%f)\n", z1.x, z1.y);
  F4(a,b,c,d,e,z1,&zz);
  printf(" F(z1) = (%f,%f)\n", zz.x, zz.y);
  printf(" z2 = (%f,%f)\n", z2.x, z2.y);
  F4(a,b,c,d,e,z2,&zz);
  printf(" F(z2) = (%f,%f)\n", zz.x, zz.y);
  printf(" z3 = (%f,%f)\n", z3.x, z3.y);
  F4(a,b,c,d,e,z3,&zz);
  printf(" F(z3) = (%f,%f)\n", zz.x, zz.y);
  printf(" z4 = (%f,%f)\n", z4.x, z4.y);
  F4(a,b,c,d,e,z4,&zz);
  printf(" F(z4) = (%f,%f)\n\n", zz.x, zz.y);

}

// end of file tcroot4.cpp