/* newton.c  newton's method to find a root, needs complex.c */
#include <stdio.h>
#include <math.h>
#include "complex.h"

complex f(complex x)
{
  return cxadd(cxmul(cxmul(x,x),cxmul(x,x)),cmplx(1.0,0.0));
}

complex fp(complex x)
{
  return cxmul(cxmul(x,x),cxmul(x,cmplx(4.0,0.0)));
}

int main(int argc, char * argv[])
{
  complex x, y;
  int i;
  
  printf("newton.c running with complex \n");
  x = cmplx(1.0,1.0);

  printf("find f(x)=0 x_next = x - f(x)/f'(x) \n", x);
  printf("expect convergence, sometimes \n");
  printf("f(x)=x^4+1  f'(x)=4*x^3 beware x=0 \n");

  printf("x=%20.15f+%20.15fi \n", x);
  printf("f(x)=%20.15f+%20.15fi \n", f(x));
  printf("fp(x)=%20.15f+%20.15fi \n", fp(x));

  for(i=0; i<12; i++)
  {
    x = cxsub(x,cxdiv(f(x),fp(x))); /* next approximation */
    printf("i=%d, x=%20.15f+%20.15fi \n", i, x);
  }
  printf("sqrt(2)/2=%20.15f \n", sqrt(2.0)/2.0);

  return 0;
}