/**********************************************************
*          Solve complex equation of degree 3:            *
*             a z^3 + b z^2 + c z + d = 0                 *
* ------------------------------------------------------- *
* SAMPLE RUN:                                             *
*             a = (1,0) b = (2,-10) c = (-4,1) d = (5,2)  *
*                                                         *
* Equation az^3 + bz^2 + cz + d = 0                       *
* Solutions z:                                            *
* z = ( -0.4703,  0.6685)                                 *
* F(z) = ( -0.0000,  0.0000)                              *
* z = (  0.6018, -0.2911)                                 *
* F(z) = (  0.0000, -0.0000)                              *
* z = ( -2.1315,  9.6226)                                 *
* F(z) = ( -0.0000,  0.0000)                              *
*                                                         *
*                                     J-P Moreau, Paris.  *
*                                     (www.jpmoreau.fr)   *
***********************************************************
  NOTE:
  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))
-------------------------------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

#define  PI  4*atan(1.0)

typedef struct {
          double x, y;
	} COMPLEX;

COMPLEX a,b,c,d,u,v,w;

// 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);
}

// z1 = a*z^3 + b*z2 + c*z + d
void F(COMPLEX z, COMPLEX *z1) {
  CMUL(a,z,&u); CMUL(z,u,&v); CMUL(z,v,&w);
  z1->x=w.x; z1->y=w.y;
  CMUL(b,z,&u); CMUL(u,z,&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;
}


void main() {

  COMPLEX det, p, q, sdet, u1, u2, v1, v2, z1;
  COMPLEX z[10];
  int i;


// seek three complex roots of complex equation:
//  a z^3 + b z^2 + c z + d = 0 

  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;

// calculate p = c/a - b²/3a²
  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,&z1);
  p.x = u.x - z1.x; p.y = u.y - z1.y;

// calculate q = 2b^3/27a^3 + d/a - bc/3a² 
  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;

// calculate det = q²/4 + p^3/27
  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

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

  for (i=1; i<10; i++) {  //substract b/3a
    z[i].x = z[i].x - v.x;
    z[i].y = z[i].y - v.y;
  }

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

  for (i=1; i<10; i++) {
    F(z[i],&z1);
    if (CABS(z1) < 0.001) {  // detect if solution is correct
      printf(" z = (%8.4f,%8.4f)\n", z[i].x, z[i].y);
      printf(" F(z) = (%8.4f,%8.4f)\n", z1.x, z1.y);
	}
  }
  printf("\n");
}

// end of file croot3.cpp