#define float double

void version(void);
int findroot(float (*f[MAX])(), float x[MAX]);
#define X_INIT_CORR 1
#define NUM_OF_STEPS 100000
#define NUM_OF_TRIES 10
#define HOW_FAST 1
#define HOW_FAST_BASE 2
#define HOW_FAST_FLAG 1
#define PRECISION 1e-10
#define DEBUG 0
int eqsolv(float how_fast, float (*f[MAX])(), float x[MAX]);
int solve(float c[MAX][MAX], float b[MAX], float x[MAX]);
float df(int i, float x[MAX], float (*f)());
#define INCREMENT 1e-10


void version(void)
{
  static int calls=0;
  if(!calls)  printf("FindRoot v1.2\n"), calls=1;
}


int findroot(float (*f[MAX])(), float x[MAX])
{
  int i,j,n;
  int error;
  int num_of_err;
  int tot_steps;
  float how_fast;
  float x_i[MAX];
  float fi, f_i2, f_f2;

  version();

  /*
    This subroutine is called FindRoot. It implements precision control
    and avoids some traps in the calculation. All you need to do is to
    define an array of pointers to functions which are to be set to zero,
    and an array of starting values of unknowns. It has several define
    options most of which are transparent. If you define DEBUG 1, it will
    print out error flag which indicates which special case is encounted
    and how it's being treated. If too many error codes other than zero
    are encounted, it might meen that a change of some defines may be
    approapriate for the case. FindRoot returns a value, which is zero,
    if the solution has been found, one, if it cannot be found for these
    starting values of unknowns (even though it varies this value - see
    define X_INIT_CORR 1), and two if the required accuracy has not been
    achieved, even though there were no outstanding errors. If you get 2,
    debugging might help.
  */

  f_i2=0;
  for(i=0;i<MAX;i++) fi=(*f[i])(x), f_i2+=fi*fi;

  tot_steps=0, num_of_err=0, how_fast=HOW_FAST;
  for(i=0;i<NUM_OF_STEPS;i++){

    if(DEBUG) for(j=0;j<MAX;j++)
      printf("x%i=%g\t\tf%i=%g\n", j, x[j], j, (*f[j])(x));

    if(f_i2<(PRECISION*PRECISION)) return 0;
    if(!HOW_FAST_FLAG) how_fast=HOW_FAST;
    for(j=0;j<NUM_OF_TRIES;j++){
      for(n=0;n<MAX;n++) x_i[n]=x[n];
      error=eqsolv(how_fast, f, x_i), tot_steps++;
      if(error) break;
      f_f2=0;
      for(n=0;n<MAX;n++) fi=(*f[n])(x_i), f_f2+=fi*fi;
      if(f_f2<f_i2) break;
      else error=MAX+1, how_fast/=HOW_FAST_BASE;

      if(DEBUG) printf("Current error code is %i, repeating with \
smaller value of how_fast\n", error);

    }

    if(DEBUG) printf("Step number %i: Current value of error code \
is %i\n", tot_steps, error);

    if(!error) num_of_err=0;
    if(num_of_err==MAX) return 1;
    if(!error){
      for(j=0;j<MAX;j++) x[j]=x_i[j]; f_i2=f_f2;
    }
    else{
      x[num_of_err]+=X_INIT_CORR, num_of_err++;
      f_i2=0;
      for(j=0;j<MAX;j++) fi=(*f[j])(x), f_i2+=fi*fi;
    }
  }
  return 2;
}


int eqsolv(float how_fast, float (*f[MAX])(), float x[MAX])
{
  int i,j;
  float c[MAX][MAX], b[MAX];

  version();

  if(how_fast>1)
    printf("warning: how_fast=%g>1, setting to 1\n", how_fast), how_fast=1;
  if(how_fast<=0)
    printf("warning: how_fast=%g<=0, setting to 1\n", how_fast), how_fast=1;

  /*
    This function solves system of equations f[i](x[j])=0.
    Each call gets better approximation. You have to initialize
    array of unknown variables x[j], though. HOW_FAST=1 implements
    Newton's method, smaller value helps avoid divergencies or
    locate multiple roots (so does init_x)
  */

  for(i=0;i<MAX;i++){
    b[i]=-how_fast*(*f[i])(x);
    for (j=0;j<MAX;j++) b[i]+=df(j, x, (*f[i]))*x[j];
  }

  for(i=0;i<MAX;i++) for(j=0;j<MAX;j++) c[i][j]=df(j, x, f[i]);

  return solve(c, b, x);
}


int solve(float c[MAX][MAX], float b[MAX], float x[MAX])
{
  int i,j,n,k;
  int flag1, flag2, error;
  float div;

  version();

  /*
    This is the actual calculation - by Gauss' method. You need to
    supply the matrix and the colomn of your linear equation. Another
    parameter (x) gives the location at which the roots are stored.
  */

  error=0;
  for(j=0;j<MAX;j++){
    flag1=0, flag2=0;
    for(i=0;i<MAX;i++) if(c[i][j]!=0){ flag1=1;
      for(k=0;k<j;k++) if(c[i][k]!=0){flag2=1; break;}
      if(flag2) flag1=0,flag2=0;
      if(flag1) break;
    }

    if(!flag1){ error++; continue;}
    div=c[i][j];
    for(k=j;k<MAX;k++) c[i][k]=c[i][k]/div;
    b[i]=b[i]/div;

    for(n=0;n<MAX;n++){
      if(n!=i){
        div=c[n][j];
        for(k=j;k<MAX;k++) c[n][k]-=c[i][k]*div;
        b[n]-=b[i]*div;
      }
    }
  }

  /*
    This is where the found roots are written in their order
  */

  if(!error) for(i=0;i<MAX;i++) for(j=0;j<MAX;j++) if(c[i][j]!=0) x[j]=b[i];

  return error;
}


float df(int i, float x[MAX], float (*f)())
{
  int j;
  static float inc=INCREMENT;
  float x1[MAX], x2[MAX];

  /*
    Finding derivative of a function f - no comment
  */

  for(j=0;j<MAX;j++){
    if(j==i) x1[j]=x[j]-inc, x2[j]=x[j]+inc;
    else x1[j]=x[j], x2[j]=x[j];
  }
  return ((*f)(x2)-(*f)(x1))/(2*inc);
}