/*****************************************************
*   ROOTS OF A SECOND DEGREE EQUATION WITH COMPLEX   *
*   COEFFICIENTS  AZ2 + BZ + C = 0                   *
* -------------------------------------------------- *
* SAMPLE RUN:                                        *
* (Find complex roots of equation:                   *
* (-4+i)z2 + (2-3i)z +(5-i) = 0).                    *
*                                                    *
* SOLVING A COMPLEX 2ND DEGREE EQUATION              *
*     az2 + bz + c = 0                               *
*                                                    *
*  a real part      = -4                             *
*  a imaginary part = 1                              *
*  b real part      = 2                              *
*  b imaginary part = -3                             *
*  c real part      = 5                              *
*  c imaginary part = -1                             *
*                                                    *
*  Root 1 =  1.444644 - 0.352759 i                   *
*  Root 2 = -0.797586 - 0.235476 i                   *
*                                                    *
* -------------------------------------------------- *
* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp  *
* and A. Reverchon (vol 2), Eyrolles, Paris, 1988"   *
* [BIBLI 03].                                        *
*                                                    *
*                 C++ Release By J-P Moreau, Paris.  *
*                       (www.jpmoreau.fr)            *
*****************************************************/
#include <stdio.h>
#include <math.h>

typedef struct {
  double x,y;
} complex;

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

void main() {
  complex a,b,c,s1,s2;
  printf("\n SOLVING A COMPLEX 2ND DEGREE EQUATION\n");
  printf("      az2 + bz + c = 0\n\n");
  printf("    a      real part = "); scanf("%lf",&a.x);
  printf("    a imaginary part = "); scanf("%lf",&a.y);
  printf("    b      real part = "); scanf("%lf",&b.x);
  printf("    b imaginary part = "); scanf("%lf",&b.y);
  printf("    c      real part = "); scanf("%lf",&c.x);
  printf("    c imaginary part = "); scanf("%lf",&c.y);

  equa2c(&a,&b,&c,&s1,&s2);

  printf("\n");
  if (s1.x > 0) printf(" Root 1 =  %f ", s1.x);
  else printf(" Root 1 = %f ", s1.x);
  if (s1.y > 0) printf("+ %f i\n", fabs(s1.y));
  else printf("- %f i\n", fabs(s1.y));
  if (s2.x > 0) printf(" Root 2 =  %f ", s2.x);
  else printf(" Root 2 = %f ", s2.x);
  if (s2.y > 0) printf("+ %f i\n", fabs(s2.y));
  else printf("- %f i\n", fabs(s2.y));
  printf("\n");
}

//end of file tequa2.cpp