/* L-25 MCS 572 Monday 13 March 2006 : Newton's method */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "dcmplx.h"

dcmplx eval ( int n, dcmplx f[n], dcmplx x );
/*
 * DESCRIPTION :
 *   Applies Horner's method to evaluate the polynomial f at x.
 *
 * ON ENTRY :
 *   n        #coefficients of the polynomial equals one plus its degree;
 *   f        complex coefficients of the polynomial,
 *            f[i] is the coefficient with the monomial x^i;
 *   x        where to evaluate the polynomial. 
 *
 * ON RETURN :
 *   y        y = eval() the function value of f at x. */

dcmplx* diff ( int n, dcmplx f[n] );
/*
 * DESCRIPTION :
 *   Allocates memory and computes the coefficients of the derivative of f. */

int newton ( int n, dcmplx f[n], dcmplx *df, dcmplx *x, float eps, int m,
             int verbose );
/*
 * DESCRIPTION :
 *   Performs at most m steps of Newton's method to find a root of f.
 *
 * ON ENTRY :
 *   n       number of coefficients in f, equals one plus the degree of f;
 *   f       complex coefficients of a polynomial in one variable;
 *   df      derivative, given as array of n-1 complex coefficients;
 *   x       initial approximation for the root;
 *   eps     accuracy requirement;
 *   m       maximal number of iterations;
 *   verbose should be 1 if output written to screen, 0 for no output.
 *
 * ON RETURN :
 *   x       improved approximation for the root;
 *   v       v = newton() return value of newton is number of iterations,
 *           if v = m+1, then the desired accuracy was not reached. */

int update_roots ( int n, dcmplx r[n], int *cnt, dcmplx x, float tol,
                   int verbose );
/*
 * DESCRIPTION :
 *   Updates r with a new approximate root x.
 *
 * ON ENTRY :
 *   n       degree of the polynomial, total size of r;
 *   r       roots currently been found;
 *   cnt     number of roots currently in r;
 *   x       new approximate root;
 *   tol     to decide whether two numbers are equal;
 *   verbose should be 1 if intermediate output is wanted, 0 otherwise.
 *
 * ON RETURN :
 *   r       updated sequence of roots;
 *   cnt     updated counter;
 *   v       v = update_roots() is index of x in r. */

int main ( int argc, char *argv[] )
{
   int n,i,max,nit,cnt;
   dcmplx *r,*f,*df,x,y;
   float tol;
   char ans;

   srand(time(NULL));
   printf("Give the degree of the polynomial : "); scanf("%d",&n);
   r = (dcmplx*)calloc(n,sizeof(dcmplx));    /* space for roots */
   f = (dcmplx*)calloc(++n,sizeof(dcmplx));  /* coefficient vector */
   for(i=0; i<n; i++) f[i] = random_dcmplx1();
   printf("The random coefficients : \n");
   for(i=0; i<n; i++,printf("\n")) write_dcmplx(f[i]);
   df = diff(n,f);
   printf("Give the maximal number of steps : "); scanf("%d",&max);
   printf("Give the required accuracy : "); scanf("%f",&tol);
   cnt = 0;
   do
   {
      x = random_dcmplx1();
      nit = newton(n,f,df,&x,tol,max,1);
      printf("The approximate root after %d Newton steps :\n",nit);
      write_dcmplx(x); printf("\n");
      y = eval(n,f,x);
      printf("The approximate root evaluated at the polynomial :\n");
      write_dcmplx(y); printf("\n");
      scanf("%c",&ans); /* skip new line symbol */
      printf("Accept this root ? (y/n) "); scanf("%c",&ans);
      if(ans == 'y') i = update_roots(n-1,r,&cnt,x,tol,1);
      scanf("%c",&ans); /* skip new line symbol */
      printf("Found %d roots, continue ? (y/n) ",cnt); scanf("%c",&ans);
   }
   while (ans == 'y');

   return 0;
}

dcmplx eval ( int n, dcmplx f[n], dcmplx x )
{
   dcmplx y = f[n-1];
   int i;

   for(i=n-2; i>=0; i--)
   {
      y = mul_dcmplx(y,x);         /* y = y*x */
      y = add_dcmplx(y,f[i]);      /* y = y + f[i] */
   }
   return y;
}

dcmplx* diff ( int n, dcmplx f[n] )
{
   dcmplx *y = (dcmplx*)calloc(n-1,sizeof(dcmplx));
   int i;

   for(i=1; i<n; i++)
      y[i-1] = mul_double(f[i],((double) i));

   return y;
}

int newton ( int n, dcmplx f[n], dcmplx *df, dcmplx *x, float eps, int m,
             int verbose )
{
   int i;
   dcmplx y,dy,dx;

   y = eval(n,f,*x);            /* y = f(x) */

   for(i=1; i<=m; i++)
   {
      dy = eval(n-1,df,*x);     /* dy = f'(x) */
      dx = div_dcmplx(y,dy);    /* dx = f(x)/f'(x) */
      *x = sub_dcmplx(*x,dx);   /* x = x - f(x)/f'(x) */
      y = eval(n,f,*x);         /* y = f(x), for stop criterium */
      if(verbose > 0)
      {
         printf("************ results at step %d ******************\n",i);
         printf("    Newton update : "); write_dcmplx(dx); printf("\n");
         printf("new approximation : "); write_dcmplx(*x); printf("\n");
         printf("   function value : "); write_dcmplx(y);  printf("\n");
      }
      if((modulus(y)<=eps) || (modulus(dx) <= eps)) return i;
   }

   return i;
}

int update_roots ( int n, dcmplx r[n], int *cnt, dcmplx x, float tol,
                   int verbose )
{
   int i;

   for(i=0; i<*cnt; i++)
      if(equal_dcmplx(r[i],x,tol) == 1)
      {
         if(verbose > 0)
         {
            printf("Found "); write_dcmplx(r[i]); printf(" at %d.\n",i);
         }
         return i;
      }
   r[(*cnt)++] = x;
   if(verbose > 0)
   {
      printf("Added "); write_dcmplx(x); printf(" as root %d.\n",*cnt);
   }
   return (*cnt-1);
}