// L-9 MCS 507 Mon 11 Sep 2023 : laguerre.c

/* Illustrates the method of Laguerre to compute a root of a polynomial.
 *
 * Compile with -O3 and -lm options, i.e.:
 * gcc -o laguerre -O3 laguerre.c -lm 
 * To redirect the output to a file, type
 * ./laguerre > 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. */

void diffpoly ( int deg, double complex *p, double complex *dp );
/*
 * Computes the derivative of the polynomial p.
 *
 * ON ENTRY :
 *   deg      degree of the polynomial p
 *   p        coefficients of a polynomial in one variable 
 *   dp       space for the coefficients of the derivative of p
 *
 * ON RETURN :
 *   dp       coefficients of the derivative of p. */

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 laguerre
 ( int deg, double complex *p, double complex *d1p, double complex *d2p,
   double complex *root, double *absdx, double *abspx, int *nbrit,
   double dxtol, double pxtol, int maxit, int verbose );
/*
 * Applies the method of Laguerre to compute a root of p.
 *
 * ON ENTRY : 
 *   deg      degree of the polynomial p
 *   p        coefficients of a polynomial in one variable 
 *   d1p      coefficients of the first derivative of p 
 *   d2p      coefficients of the second derivative of p 
 *   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 ( void );
/*
 * Runs a simple test on a quadric. */

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

int main ( int argc, char *argv[] )
{
   // simple_test();
   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];
}

void diffpoly ( int deg, double complex *p, double complex *dp )
{
   int i;

   for(i=1; i<=deg; i++) dp[i-1] = p[i]*i;
}

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;
}
 
int laguerre
 ( int deg, double complex *p, double complex *d1p, double complex *d2p,
   double complex *root, double *absdx, double *abspx, int *nbrit,
   double dxtol, double pxtol, int maxit, int verbose )
{
   const int degm1 = deg-1;
   int i;
   double complex pval,d1val,d2val,Lroot,Mroot,valsqrt;
   double complex yplus,yminus,dx;

   if(verbose > 0)
   {
      printf("step :        real(root)               imag(root)");
      printf("          |dx|     |p(x)|\n");
      printf(" %3d : %23.16e  %23.16e\n", 0, creal(*root), cimag(*root));
   }
   evalpoly(deg,p,*root,&pval);
   *abspx = cabs(pval);
   for(i=1; i<=maxit; i++)
   {
      if(*abspx < pxtol)
      {
         if(verbose > 0) printf("succeeded after %d step(s)\n", i-1);
         return 0;
      }
      evalpoly(degm1,d1p,*root,&d1val);
      evalpoly(deg-2,d2p,*root,&d2val);
      Lroot = d1val/pval;
      Mroot = Lroot*Lroot - d2val/pval;
      valsqrt = csqrt(degm1*(deg*Mroot - Lroot*Lroot));
      yplus = Lroot + valsqrt;
      yminus = Lroot - valsqrt;
      if(cabs(yplus) > cabs(yminus))
         dx = deg/yplus;
      else
         dx = deg/yminus;
      *root = *root - dx;
      evalpoly(deg,p,*root,&pval);
      *absdx = cabs(dx); *abspx = cabs(pval);
      if(verbose > 0)
      {
         printf(" %3d : %23.16e  %23.16e", i, creal(*root), cimag(*root));
         printf("  %.2e  %.2e\n", *absdx, *abspx);
      }
      if(*absdx < dxtol)
      {
         if(verbose > 0) printf("succeeded after %d step(s)\n", i);
         return 0;
      }
   }
   return 1;
}

void simple_test ( void )
{
   double complex *cff,*dp1,*dp2,*roots;
   double complex x,y,z;
   double absdx,abspx,rnd;
   int i,nit;

   cff = (double complex*)calloc(3,sizeof(double complex));
   dp1 = (double complex*)calloc(2,sizeof(double complex));
   dp2 = (double complex*)calloc(1,sizeof(double complex));
   roots = (double complex*)calloc(2,sizeof(double complex));

   cff[0] = 2.0 + 0.0*I;
   cff[1] = -3.0 + 0.0*I;
   cff[2] = 1.0 + 0.0*I;

   x = 1.0 + 0.0*I; roots[0] = x;
   evalpoly(2,cff,x,&y);
   printf("x = %.2e + %.2e\n", creal(x), cimag(x));
   printf("y = %.2e + %.2e\n", creal(y), cimag(y));

   x = 2.0 + 0.0*I; roots[1] = x;
   evalpoly(2,cff,x,&y);
   printf("x = %.2e + %.2e\n", creal(x), cimag(x));
   printf("y = %.2e + %.2e\n", creal(y), cimag(y));

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

   diffpoly(2,cff,dp1);

   printf("coefficients of the first derivative :\n");
   for(i=0; i<2; i++) 
      printf("%d : %.2e + %.2e\n", i, creal(dp1[i]), cimag(dp1[i]));

   diffpoly(1,dp1,dp2);

   printf("coefficients of the second derivative :\n");
   for(i=0; i<1; i++) 
      printf("%d : %.2e + %.2e\n", i, creal(dp2[i]), cimag(dp2[i]));

   srand(time(NULL));

   z = ((double) rand())/RAND_MAX + (((double) rand())/RAND_MAX)*I;

   rnd = ((double) rand())/RAND_MAX;
   printf("rnd : %.2f\n", rnd);
   if(rnd > 0.5) z = 3.0*z;

   printf("a random complex number :\n");
   printf("%23.16e  %23.16e\n", creal(z), cimag(z));

   laguerre(3,cff,dp1,dp2,&z,&absdx,&abspx,&nit,1.0e-8,1.0e-8,20,1);

   printf("rank : %d\n",rank(2,roots,z,1.0e-4));
}

void make_matrix ( void )
{
   double complex *cff,*dp1,*dp2,*roots,z0;
   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));
   dp1 = (double complex*)calloc(8,sizeof(double complex));
   dp2 = (double complex*)calloc(7,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;
 
   diffpoly(8,cff,dp1);
   diffpoly(7,dp1,dp2);

   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
         {

            laguerre(8,cff,dp1,dp2,&z0,&absdx,&abspx,&nit,1.0e-8,1.0e-8,20,0);
            mat[i][j] = rank(8,roots,z0,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");
   }
}