#define VERSION 0.01

#define PRECISION 1e-4
#define PRECISION_I 1e-4
#define MAX_STEPS_F 1000
#define MAX_STEPS_I 1000

#define float double

float df(int i, float x[MAX], float (*f)());
#define INCREMENT 1e-5
void version(void);


float fy(float (*f)(), float x[MAX], float n[MAX], float y)
{
  int i;
  float x_aux[MAX];

  for(i=0;i<MAX;i++) x_aux[i]=x[i]+y*n[i];
  return (*f)(x_aux);
}

int findmin(float (*f)(), float x[MAX], float x_min[MAX], float x_max[MAX])
{
  int i, j, k;
  int flag;
  float n[MAX];
  float x_inc[MAX];
  float error;
  float y_mid, y_min, y_max, y_aux, y_min_aux, y_max_aux;
  float f_mid, f_min, f_max, f_aux;

  for(i=0;i<MAX;i++) if((x[i]<x_min[i])||(x[i]>x_max[i])) break;
  if(i!=MAX) return -1;

  for(i=0;i<MAX;i++) if(x_min[i]==x_max[i]) break;
  if(i!=MAX) return -2;

  for(i=0;i<MAX_STEPS_F;i++){
    flag=0;
    for(j=0;j<MAX;j++){
      n[j]=-df(j, x, f);
      printf("n[%i]=%g\n", j, n[j]);
      if(n[j]!=0) flag=1;
    }
    if(!flag) return 1;

    y_max=y_min=y_mid=y_aux=0;

    for(j=0;j<MAX;j++){
      if(n[j]!=0){
	y_min_aux=(x_min[j]-x[j])/n[j];
	y_max_aux=(x_max[j]-x[j])/n[j];
	if(y_min_aux>y_max_aux){
	  y_aux=y_max_aux;
	  y_max_aux=y_min_aux;
	  y_min_aux=y_aux;
	}
	if((y_min==0)||(y_min_aux>y_min)) y_min=y_min_aux;
	if((y_max==0)||(y_max_aux<y_max)) y_max=y_max_aux;
      }
    }

    printf("HELLO %g %g %g %g\n", y_min, y_mid, y_max, y_aux);

    f_mid=fy(f, x, n, y_mid);
    f_min=fy(f, x, n, y_min);
    f_max=fy(f, x, n, y_max);
    if((f_min<f_mid)||(f_max<f_mid)) return -3;
    for(j=0;j<MAX_STEPS_I;j++){
      if(j%2){
	y_aux=(y_mid+y_min)/2;
	f_aux=fy(f, x, n, y_aux);
	if(f_aux<f_mid){
	  y_max=y_mid, f_max=f_mid;
	  y_mid=y_aux, f_mid=f_aux;
	}
	else y_min=y_aux, f_min=f_aux;
      }
      else{
	y_aux=(y_mid+y_max)/2;
	f_aux=fy(f, x, n, y_aux);
	if(f_aux<f_mid){
	  y_min=y_mid, f_min=f_mid;
	  y_mid=y_aux, f_mid=f_aux;
	}
	else y_max=y_aux, f_max=f_aux;
      }

      if(fabs(y_max-y_min)<PRECISION) break;


      printf("Step # %i\ty_min=%g, y_mid=%g, y_max=%g, y_aux=%g\n", j, y_min, y_mid, y_max, y_aux);
      printf("Step # %i\tf_min=%g, f_mid=%g, f_max=%g, f_aux=%g\n", j, f_min, f_mid, f_max, f_aux);

      for(k=0;k<MAX;k++) printf("x%i = %g\t\t", k, x[k]+y_mid*n[k]);
      printf("\nf(x)=%g\n\n", f_mid);
    }


    printf("j=%i\n", j);
    if(j==MAX_STEPS_I) printf("Warning: number of cycles in linear minimization exceeded %i\n", MAX_STEPS_I);

    error=0;
    for(j=0;j<MAX;j++){
      x[j]+=x_inc[j]=y_mid*n[j];
      error+=x_inc[j]*x_inc[j];
      printf("%g\t\t", x_inc[j]);
    }
    error=sqrt(error);
    printf("Error=%g\n", error);
    if(error<PRECISION) break;
  }
  printf("i=%i\n", i);

  if(i==MAX_STEPS_F) return 1;
  return 0;
}


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

  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);
}


void version(void)
{
  static int calls=0;
  if(!calls)  printf("FindRoot v%g\n",VERSION), calls=1;
}