/*
   Корни уравнений. Функция zero()
*/
#include <stdio.h>
#include <math.h>
#include <float.h>
/*
   Функция zero() - Комбинация методов хорд и квадратичной интерполяции.

   Приведенные ниже схемы, иллюстрируют существо метода:

        начальное положение:  ---|----------------------------|----->x
                                 a                            b

          после метода хорд:  ---|--------|-------------------|----->x
                                 a        c                   b

  после квадр. интерполяции:  ---|--------|---------|---------|----->x
                                 a        c         d         b

     итоги одного шага:

 а) f(c),f(d)-разные знаки:       --------|---------|------>x
                                          a         b

 б) f(c),f(d)-одинаковые,
    f(a),f(c)-разные знаки:     -|------------------|------>x
                                 a                  b

 в) f(a),f(c),f(d)-одинаковые
                     знаки:   ---|--------|---------|------>x
                                 a        c         b

 В зависимости от того, сколько точек мы имеем после очередного шага (2 или 3),
 следующий шаг начинается либо с метода хорд (при 2 точках), либо прямо с мето-
 да квадратичной интерполяции.

*/
int zero(float a, float b, float eps, float f(float), float *x, float *fx)
{
  float c,d,fa,fb,fc,fd,a0,a1,a2,ak,bk,ck,sa,d1,d2,buf;

// правильно задан интервал (a,b)?
  fa=f(a); fb=f(b);
  if(fa*fb>0) return 1;  // корня нет
  eps+=FLT_EPSILON; // защищаемся от eps=0
  c=FLT_MAX;
  do
    {
      if(c==FLT_MAX)  // запускаем метод хорд
       { c=a-fa*(b-a)/(fb-fa); fc=f(c); }
      if((c-a<eps)||(b-c<eps)||(fabs(fc)<1e-6)) { *x=c; *fx=fc; return 0; }
// квадратичная интерполяция
      a2=a*(fc-fb)+b*(fa-fc)+c*(fb-fa);
      ak=a*a; bk=b*b; ck=c*c;
      a1=ak*(fb-fc)+bk*(fc-fa)+ck*(fa-fb);
      a0=ak*(b*fc-c*fb)+bk*(c*fa-a*fc)+ck*(a*fb-b*fa);
      sa=sqrt(a1*a1-4*a2*a0);
      d1=-(a1+sa)/(2*a2); d2=-(a1-sa)/(2*a2);
      if((d1>a)&&(d1<b)) d=d1;    // ищем подходящий корень
       else if((d2>a)&&(d2<b)) d=d2;
             else if(fa*fc<0) { a=c; fa=fc; }  // нет подходящего корня!
                   else { b=c; fb=fc; }
      fd=f(d);
      if((d-a<eps)||(b-d<eps)) { *x=d; *fx=fd; return 0; }
      if(c>d) { buf=c; c=d; d=buf; buf=fc; fc=fd; fd=buf; }
      if(fc*fd<0)  // разные знаки
           { a=c; fa=fc; b=d; fb=fd; c=FLT_MAX; }
       else if(fa*fc<0) { b=d; fb=fd; }
             else { a=c; fa=fc; c=d; fc=fd; }
    }
  while((b-a)>eps);
  if(fabs(fa)>fabs(fb)) { *x=b; *fx=fb; }
   else { *x=a; *fx=fa; }
  return 0;  // нашли корень!
}

int nu,kf;

// Функция f() - вычисление значений функций
float f(float x)
{
  float x2,x3,x4,x5,fx;

  kf++;
  x2=x*x; x3=x2*x; x4=x2*x2; x5=x3*x2;
  switch(nu+1)
     {
       case  1: fx=x3+1; break;
       case  2: fx=x3-6*x+2; break;
       case  3: fx=x3+x-10000; break;
       case  4: fx=2*x3+3*x2+4*x+5; break;
       case  5: fx=x4-4*x-1; break;
       case  6: fx=2*x4-3*x2+75*x-10000; break;
       case  7: fx=x5-x-0.2; break;
       case  8: fx=x5+x4+x-1; break;
       case  9: fx=x5-10*x+128; break;
       case 10: fx=x+exp(x); break;
       case 11: fx=x*log(x)-1; break;
       case 12: fx=tan(x)-x; break;
       case 13: fx=x-sin(x)-0.25; break;
       case 14: fx=2*x-exp(-0.1*x); break;
       case 15: fx=x-exp(-0.1*x); break;
       case 16: fx=x+log(x); break;
       case 17: fx=x-5*log(x); break;
       case 18: fx=x+log(x)-0.125*exp(x); break;
       case 19: fx=sin(x)+cos(2*x); break;
       case 20: fx=(x5-2*x4)*x2+7*x-8; break;
       case 21: fx=(x5-2*x4)*x4+9*x-10; break;
       case 22: fx=x4-5*x2+8*x-8;
     }
  return fx;
}

void main()
{
  float a[22]={-2,2,21,-2,-1,8,1,0,-3,-1,1.5,3,0,0.2,0.2,0.5,0.5,2.8,2,0,1,-2};
  float b[22]={0,3,22,-1,1,9,2,1,-2,0,2,4.6,1.5,1.5,2,1.5,1.5,4,4,2,2,2};
  float eps=0.0;
  float x,fx;
  int ir,kff=0;
  FILE *fp;

  fp=fopen("koi_mf.txt","w");
  fprintf(fp,"\n\nур-ие     a       b     kf           x              fx  ");
  for(nu=0; nu<22; nu++)
    {
      kf=0;
      ir=zero(a[nu],b[nu],eps,f,&x,&fx);
      if(ir) fprintf(fp,"\n  %2d\t %3.1g\t %3.1g\t  - нет корня !!! ",
                                                     nu+1,a[nu],b[nu]);
       else fprintf(fp,"\n %2d\t% 3.1f\t% 3.1f\t%2d\t% 8.7f\t% 8.3e ",
                                             nu+1,a[nu],b[nu],kf,x,fx);
      kff+=kf;
    }
  fprintf(fp,"\n Расчет закончен! (kff=%d) ",kff);
  fclose(fp);
}
