/*
   raiz.Complex.c
   ==============

   Este programa encontra uma aproximacao para uma raiz quadrada de um numero
   complexo dado.  O programa usa o metodo de Newton.  Muito parecido com o
   raiz.c 

 */

#include <stdio.h>
#include <math.h>

#define EPSILON .0000000001

/*
  Prototipos e definicoes necessarias para a manipulacao
  de numeros complexos, do livro de Sedgewick (veja abaixo)
 */

typedef struct { double Re; double Im; } Complex;
Complex COMPLEXinit(double, double);
 double Re(Complex);
 double Im(Complex);
Complex COMPLEXmult(Complex, Complex);

/* Fim do codigo de Sedgewick */

/* Codigo adicionado */

Complex COMPLEXdiv(Complex, Complex);
Complex COMPLEXdiff(Complex, Complex);
Complex COMPLEXadd(Complex, Complex);
Complex COMPLEXconj(Complex);

   void COMPLEXprint(Complex);
   void COMPLEXscan(Complex *);

    int COMPLEXclose(Complex, Complex);

/************************************************************/

Complex meu_sqrt(Complex z, Complex z0, int n);

/************************************************************/

int main()
{
  Complex z, z0, w;
  int n;

  printf("z: ");
  COMPLEXscan(&z);
  printf("z0: ");
  COMPLEXscan(&z0);
  printf("n: ");
  scanf("%d", &n);
  
  w = meu_sqrt(z, z0, n);
  printf("Raiz que encontrei: ");
  COMPLEXprint(w);
  putchar('\n');

  printf("Quadrado: ");
  COMPLEXprint(COMPLEXmult(w, w));
  putchar('\n');

  return 0;
}

/*
   meu_sqrt(): devolve uma aproximacao da raiz quadrada
   de x.  Usa o metodo de Newton; n 'e o numero maximo
   de vezes que a recorrencia do metodo de Newton 'e
   aplicada.  Caso a diferenca entre duas aproximacoes
   sucessivas for pequena, paramos a iteracao e
   devolvemos o valor atual.
 */

Complex meu_sqrt(Complex z, Complex z0, int n)
{
  Complex y = z0, yy = z0, t, meio = COMPLEXinit(.5,0);

  while (n > 0) {
    t = COMPLEXadd(y, COMPLEXdiv(z, y));
    yy = COMPLEXmult(meio, t);
    n--;
    if (COMPLEXclose(yy, y))
      break;
    y = yy;
  }

  return yy;
}

/*
  Implementacoes das funcoes para manipulacao de numeros
  complexos.  Este codigo 'e de um excelente livro, adequado
  para um segundo (e terceiro) curso de programacao:

    This code is from "Algorithms in C, Third Edition," 
    by Robert Sedgewick, Addison Wesley Longman, 1998." 
    (com algumas edicoes e adicoes)
 */

Complex COMPLEXinit(double Re, double Im)
  { Complex t; t.Re = Re; t.Im = Im; return t; }
double Re(Complex z)
  { return z.Re; }
double Im(Complex z)
  { return z.Im; }
Complex COMPLEXmult(Complex a, Complex b)
  { Complex t; 
    t.Re = a.Re*b.Re - a.Im*b.Im; 
    t.Im = a.Re*b.Im + a.Im*b.Re; 
    return t; 
  }

/* Fim do codigo de Sedgewick */

/* Codigo adicionado */

Complex COMPLEXdiv(Complex a, Complex b)
{
  Complex t;
  double r = b.Re*b.Re + b.Im*b.Im;
  t.Re = (a.Re*b.Re + a.Im*b.Im)/r;
  t.Im = (-a.Re*b.Im + a.Im*b.Re)/r;
  return t; 
}

Complex COMPLEXdiff(Complex z, Complex w)
{
  Complex t;
  t.Re=z.Re-w.Re;
  t.Im=z.Im-w.Im;
  return t;
}

Complex COMPLEXadd(Complex z, Complex w)
{
  Complex t;
  t.Re=z.Re+w.Re;
  t.Im=z.Im+w.Im;
  return t;
}

Complex COMPLEXconj(Complex z)
{
  Complex t;
  t.Re=z.Re;
  t.Im=-z.Im;
  return t;
}

void COMPLEXprint(Complex z)
{
  printf("%6.3f %6.3f", z.Re, z.Im);
}

void COMPLEXscan(Complex *z)
{
  scanf("%lf %lf", &((*z).Re), &((*z).Im));
}

int COMPLEXclose(Complex z, Complex zz)
{
  return
    (Re(COMPLEXdiff(z,zz)) < EPSILON)&&(-Re(COMPLEXdiff(z,zz)) < EPSILON)
    && 
    (Im(COMPLEXdiff(z,zz)) < EPSILON)&&(-Im(COMPLEXdiff(z,zz)) < EPSILON);
}