// P-1 MCS 572 due Fri 27 Sep 2024 : muller.c

/* Illustrates the method of Muller to compute a root of a polynomial.
 *
 * Compile with -O3 and -lm options, i.e.:
 * gcc -o muller -O3 muller.c -lm 
 * To redirect the output to a file, type
 * ./muller > output.txt
 * then output.txt contains the file with the dimension
 * and all numbers separated by commas. */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <complex.h>

void evalpoly
 ( int deg, double complex *p, double complex x, double complex *y );
/*
 * Evaluates the polynomial p at x, y = p(x).
 *
 * ON ENTRY :
 *   deg      degree of the polynomial p
 *   p        coefficients of a polynomial in one variable 
 *   x        argument for p 
 *
 * ON RETURN :
 *   y        values of p at x. */

int muller
 ( int deg, double complex *p,
   double complex x0, double complex x1, double complex x2,
   double complex *root, double *absdx, double *abspx, int *nbrit,
   double dxtol, double pxtol, int maxit, int verbose );
/*
 * Applies the method of Muller to compute a root of p.
 *
 * ON ENTRY : 
 *   deg      degree of the polynomial p
 *   p        coefficients of a polynomial in one variable 
 *   x0       start point
 *   x1       point distinct from the start point
 *   x2       point distinct from x0 and x1
 *   root     an approximation for the root 
 *   dxtol    the tolerance on the forward error 
 *   pxtol    the tolerance on the backward error 
 *   maxit    the maximum number of iterations
 *   verbose  the verbose flag, if true, writes one line each step
 *
 * ON RETURN :
 *   root     an approximation for the root 
 *   absdx    the estimated forward error 
 *   abspx    the estimated backward error 
 *   nbrit    the number of iterations
 *   fail     is the integer returned by the function,
 *            is 1 if tolerances not reached, 0 otherwise. */

void simple_test ( int d, int verbose);
/*
 * Runs a simple test on a polyomial of degree d. */

void make_matrix ( void );
/*
 * Makes a matrix for p = x^8 - 1 to represent the attraction basins
 * of the method of Laguerre. */

int rank ( int n, double complex *r, double complex z, double tol );
/*
 * Returns the rank of z in the array r of size n.
 *
 * ON ENTRY :
 *   n        number of elements in r
 *   r        n complex numbers
 *   z        some complex number
 *   tol      to decide if two complex numbers are the same.
 *
 * ON RETURN :
 *   the index of z in r, or -1 if z does not appear in r. */

int main ( int argc, char *argv[] )
{
   // simple_test(6, 1);
   make_matrix();

   return 0;
}

void evalpoly
 ( int deg, double complex *p, double complex x, double complex *y )
{
   int i;

   *y = p[deg];

   for(i=deg-1; i>=0; i--) *y = (*y)*x + p[i];
}

int muller
 ( int deg, double complex *p,
   double complex x0, double complex x1, double complex x2,
   double complex *root, double *absdx, double *abspx, int *nbrit,
   double dxtol, double pxtol, int maxit, int verbose )
{
   int i;
   double complex fx0,fx1,fx2,proot,dx,h0,h1,d0,d1;
   double complex a,b,c,rad,den;

   evalpoly(deg,p,x0,&fx0);
   evalpoly(deg,p,x1,&fx1);
   evalpoly(deg,p,x2,&fx2);
   *root = x2;
   proot = fx2;
   dx = 1.0;
   *absdx = cabs(dx);
   *abspx = cabs(proot);

   if(verbose > 0)
   {
      printf("step :        real(root)               imag(root)");
      printf("          |dx|     |p(x)|\n");
      printf(" %3d : %23.16e  %23.16e", 0, creal(*root), cimag(*root));
      printf("  %.2e  %.2e\n", *absdx, *abspx);
   }
   for(i=1; i<=maxit; i++)
   {
      // the quadric q passes through (x0, fx0), (x1, fx1), (x2, fx2)
      h0 = x1 - x0;
      h1 = x2 - x1;
      d0 = (fx1 - fx0)/h0;
      d1 = (fx2 - fx1)/h1;
      // q = a*x^2 + b*x + c
      a = (d1 - d0)/(h1 + h0);
      b = a*h1 + d1;
      c = fx2;
      // apply the quadratic formula
      rad = csqrt(b*b - 4*a*c);
      // the closest root has the largest denominator
      if(abs(b + rad) > abs(b - rad))
         den = b + rad;
      else
         den = b - rad;
      dx = -2.0*c/den;
      *root = x2 + dx;
      evalpoly(deg,p,*root,&proot);
      *absdx = cabs(dx);
      *abspx = cabs(proot);
      if(verbose > 0)
      {
         printf(" %3d : %23.16e  %23.16e", i, creal(*root), cimag(*root));
         printf("  %.2e  %.2e\n", *absdx, *abspx);
      }
      if((*absdx < dxtol) | (*abspx < pxtol))
      {
         if(verbose > 0) printf("succeeded after %d step(s)\n", i);
         *nbrit = i;
         return 0;
      }
      else
      {
         x0 = x1;    fx0 = fx1;
         x1 = x2;    fx1 = fx2;
         x2 = *root; fx2 = proot;
      }
   }
   if(verbose > 0) printf("failed after %d step(s)\n", maxit);

   return 1;
}

void simple_test ( int d, int verbose )
{
   double complex *cff;
   double complex x0,x1,x2,root;
   double absdx,abspx,rnd;
   int i,nit,fail;

   srand(time(NULL));

   cff = (double complex*)calloc(d,sizeof(double complex));
   for(i=0; i<d; i++)
      cff[i] = ((double) rand())/RAND_MAX + (((double) rand())/RAND_MAX)*I;

   printf("the coefficients of a random polynomial :\n");
   for(i=0; i<d; i++) 
      printf("  %d : %.16e + %.16e\n", i, creal(cff[i]), cimag(cff[i]));

   x0 = ((double) rand())/RAND_MAX + (((double) rand())/RAND_MAX)*I;
   x1 = ((double) rand())/RAND_MAX + (((double) rand())/RAND_MAX)*I;
   x1 = x0 + 1.0e-04*x1;
   x2 = ((double) rand())/RAND_MAX + (((double) rand())/RAND_MAX)*I;
   x2 = x0 + 1.0e-04*x2;

   printf("the start points : \n");
   printf("  x0 : %.16e + %.16e\n", creal(x0), cimag(x0));
   printf("  x1 : %.16e + %.16e\n", creal(x1), cimag(x1));
   printf("  x2 : %.16e + %.16e\n", creal(x2), cimag(x2));

   fail = muller(d,cff,x0,x1,x2,&root,&absdx,&abspx,&nit,1.0e-8,1.0e-8,10,
                 verbose);

   printf(" approximation  : %.16e  %.16e\n", creal(root), cimag(root));
   printf("  forward error : %.2e\n", absdx);
   printf(" backward error : %.2e\n", abspx);
   printf("number of steps : %d\n", nit);
   printf("         failed : %d\n", fail);
}

int rank ( int n, double complex *r, double complex z, double tol )
{
   int i;

   for(i=0; i<n; i++)
      if(cabs(r[i] - z) < tol) return i;

   return -1;
}

void make_matrix ( void )
{
   double complex *cff,*roots,z0,z1,z2,rt;
   double absdx,abspx;
   const double sq2 = sqrt(2.0)/2.0;
   const int dim = 10001; // dimension of the matrix
   const double left = -1.1;
   const double right = 1.1;
   const double dz = (right - left)/(dim-1);
   int **mat;
   int i,j,nit;

   cff = (double complex*)calloc(9,sizeof(double complex));
   roots = (double complex*)calloc(8,sizeof(double complex));
   mat = (int**)calloc(dim,sizeof(int*));
   for(i=0; i<dim; i++)
   {
      mat[i] = (int*) calloc(dim,sizeof(int));
      for(j=0; j<dim; j++) mat[i][j] = 0;
   }

   cff[0] =  1.0 + 0.0*I; cff[1] = 0.0 + 0.0*I; cff[2] = 0.0 + 0.0*I;
   cff[3] =  0.0 + 0.0*I; cff[4] = 0.0 + 0.0*I; cff[5] = 0.0 + 0.0*I;
   cff[6] =  0.0 + 0.0*I; cff[7] = 0.0 + 0.0*I; cff[8] = -1.0 + 0.0*I;
 
   roots[0] =  1.0 + 0.0*I; roots[1] =  sq2 + sq2*I; roots[2] =  0.0 + 1.0*I;
   roots[3] = -sq2 + sq2*I; roots[4] = -1.0 + 0.0*I; roots[5] = -sq2 - sq2*I;
   roots[6] =  0.0 - 1.0*I; roots[7] =  sq2 - sq2*I;
/*
   printf("the roots :\n");
   for(i=0; i<8; i++)
      printf("  %23.16e  %23.16e\n",creal(roots[i]), cimag(roots[i]));
 */
   for(i=0; i<dim; i++)
      for(j=0; j<dim; j++)
      {
         z0 = (left + i*dz) + (left + j*dz)*I;
         if(rank(8,roots,z0,0.01) != -1)
            mat[i][j] = 8; 
         else
         {
            z1 = ((double) rand())/RAND_MAX + (((double) rand())/RAND_MAX)*I;
            z1 = z0 + 1.0e-04*z1;
            z2 = ((double) rand())/RAND_MAX + (((double) rand())/RAND_MAX)*I;
            z2 = z0 + 1.0e-04*z2;
            muller(8,cff,z0,z1,z2,&rt,&absdx,&abspx,&nit,1.0e-8,1.0e-8,10,0);
            mat[i][j] = rank(8,roots,rt,1.0e-4);
         }
      }

   printf("%d\n",dim);
   for(i=0; i<dim; i++) // write comma separated integers
   {
      printf("%d", mat[i][0]);

      for(j=1; j<dim; j++) printf(",%d", mat[i][j]);
      
      printf("\n");
   }
}