Obtaining integral of multiparametrical functions in C++

All replies

• 

  

  0

  Sign in to vote

  Example of using

  testfunc.cpp

  #include <stdio.h>
  #include <stdlib.h>
  #include <windows.h>
  #include <conio.h>
  #include <math.h>
  //#include "horndiagf.h"
  //#include "integral1.h"
  #include "cusparse.h"

  
  int main(void)

  {

  double psi0e,psi0h, resnuah,resnuae,a,aroz,b,broz,psih,psie,RH,RE  ;
  //horn parameters
  psih=1.07;
  psie=23.75;
  a=0.0721;
  aroz=0.080;

  RH=2.140;
  RE=0.25;
  b=0.034;
  broz=0.220;

  
  psi0e=20.62; //psi0dop
  psi0h=69;    //psi0osn

  double lyambda0=0.1;
  double lyambdamin=0.10309;
  double lyambdamax=0.0971;
  //0.856 0.863 0.870 0.86 0.88 0.89
  printf("\n Main unit  entry OK ");

  printf ("\n MAIN: psi0e=%lf deg", psi0e);
  printf ("\n MAIN: broz=%lf m", broz);
  printf ("\n MAIN: Lyambda=%lf m", lyambda0);

  printf("\n Run get_nuah(double psi0E)  ");

  
  double FE1PSI0DOP=FOE(psi0e*M_PI/180,broz,lyambdamin);
  double FE2PSI0DOP=FOE(psi0e*M_PI/180,broz,lyambda0);
  double FE3PSI0DOP=FOE(psi0e*M_PI/180,broz,lyambdamax);
  printf("\n psi0dop= %lf, FE1(psi0dop)=%lf, lyambdamin=%lf m ",psi0e,FE1PSI0DOP,lyambdamin );
  printf("\n psi0dop= %lf, FE2(psi0dop)=%lf, lyambdamin=%lf m ",psi0e,FE2PSI0DOP,lyambda0 );
  printf("\n psi0dop= %lf, FE3(psi0dop)=%lf, lyambdamin=%lf m ",psi0e,FE3PSI0DOP,lyambdamax);

  double FH1PSI0DOP=FOH(psi0h*M_PI/180,aroz,lyambdamin);
  double FH2PSI0DOP=FOH(psi0h*M_PI/180,aroz,lyambda0);
  double FH3PSI0DOP=FOH(psi0h*M_PI/180,aroz,lyambdamax);
  printf("\n psi0osn= %lf, FH1(psi0osn)=%lf, lyambdamin=%lf m ",psi0h,FH1PSI0DOP,lyambdamin);
  printf("\n psi0osn= %lf, FH2(psi0osn)=%lf, lyambda0=%lf m ",psi0h,FH2PSI0DOP,lyambda0 );
  printf("\n psi0osn= %lf, FH3(psi0dop)=%lf, lyambdamax=%lf m ",psi0h,FH3PSI0DOP,lyambdamax );

  
  double resnuah1=get_nuah( psi0h,aroz, lyambdamin);

  double resnuah2=get_nuah( psi0h,aroz, lyambda0);

  double resnuah3=get_nuah( psi0h,aroz, lyambdamax);
  printf("\n result : nua h =%lf , lyambdamin=%lf m ",resnuah1,lyambdamin);
  printf("\n result : nua h =%lf , lyambda0  =%lf m ",resnuah2,lyambda0);
  printf("\n result : nua h =%lf , lyambdamax=%lf m ",resnuah3,lyambdamax);

  
  double resnuae1=get_nuae( psi0e,broz, lyambdamin);
  double resnuae2=get_nuae( psi0e,broz, lyambda0);
  double resnuae3=get_nuae( psi0e,broz, lyambdamax);
  printf("\n result : nua e 1 =%lf , lyambdamin=%lf m ",resnuae1,lyambdamin);
  printf("\n result : nua e 2 =%lf , lyambda0  =%lf m ",resnuae2,lyambda0);
  printf("\n result : nua e 3 =%lf , lyambdamax=%lf m  ",resnuae3,lyambdamax);
  getch();
  return 0;

  
  }

  cusparse.h

  #include <stdio.h>
  #include <stdlib.h>
  #include <windows.h>
  #include <conio.h>
  #include <math.h>
  #include "horndiagf.h"
  #include "integral1.h"

   

  double funcint1_h(double Teta, double ar, double lw)
  {
   // ar for H
  return   FOH(Teta,ar,lw)*tan(0.5*Teta);
  }

  double funcint2_h(double Teta, double ar, double lw)
  {
    // ar for H
   //psi0H,aroz,cos in FOH
  return    pow( FOH(Teta,ar, lw),2)*sin(Teta);
  }

  
  double get_nuah(double Psi0H, double ar, double lw )
  {
   
  double eps=0.0000000001;
  printf ("\n Obtaining nu a H  \n");
  Psi0H=Psi0H*M_PI/180;
   double psimin ;
   psimin=0;
  return  2*pow(tan(0.5*Psi0H),-2)*pow(fabs(IntRomberg(psimin,Psi0H,eps,ar,lw,funcint1_h ) ),2)/ IntRomberg(psimin,Psi0H,eps,ar,lw,funcint2_h ) ;     
   
  }

  double funcint1_e(double Psi, double br, double Lw)
  {
  //for FE(Teta)  : br, Teta, lw , psi is Teta
  //for integral (0;Psi0E; f(Psi); dPsi) 
  //f(Psi)=FoE(Psi)*tg(0.5*Psi) 
  //printf("F1e"); 

  
  return  FOE(Psi ,br,Lw)*tan(0.5*Psi);
  }

  double funcint2_e(double Psi, double br, double lw)
  {
  //for FE(Teta)  : br, Teta, lw , psi is Teta
  //for integral (0;Psi0; f(Psi); dPsi) 
  //f(Psi)=(FoE(Psi)^2)*sin( Psi)
    
  return pow( FOE(Psi,  br, lw),2)*sin(Psi);
  }

  
  double get_nuae(double Psi0E, double br, double LW )

  {
  //LW-lyambda for calculation
   
  double eps=0.0000000001;
  Psi0E=Psi0E*M_PI/180;
  printf ("\n Obtaining nu a E  \n");
   double funcint1_e(double, double,double);
   double funcint2_e(double, double,double);
  // double IntRomberg(double ,double, double ,double ,double , double(*)(double,double,double) );
   double psimin ;
   psimin=1e-15;
  return    2*pow(tan(0.5*Psi0E),-2)*fabs(pow( IntRomberg(psimin,Psi0E,eps,br,LW,funcint1_e ),2))  / IntRomberg(psimin,Psi0E,eps,br,LW,funcint2_e ) ;     
   
  }

  double subfunctionf0psie(double Psi,double  razm,double lw)
  { 
  return pow(FOE(Psi, razm,lw),2)*sin(0.5*Psi);
  }

  double subfunctionf0psih(double Psi, double  razm,double lw)
  { 
  return pow(FOH(Psi, razm,lw),2)*sin(0.5*Psi);
  }

  double getnu0e(double Psi0, double br, double LW )
  {
  Psi0=Psi0*M_PI/180; 
  double eps=0.000000001; 
  double psimin ;
   psimin=1e-15;
  return IntRomberg(psimin,Psi0,eps,br,LW,subfunctionf0psie) /IntRomberg(psimin,M_PI,eps,br,LW,subfunctionf0psie);  
  }

  double getnu0h(double Psi0, double ar, double LW )
  {
  Psi0=Psi0*M_PI/180; 
  double eps=0.000000001; 
  double psimin ;
   psimin=1e-15;
  return IntRomberg(psimin,Psi0,eps,ar,LW,subfunctionf0psih) /IntRomberg(psimin,M_PI,eps,ar,LW,subfunctionf0psih);  
  }

  horndiagf.h

  #include <stdio.h>
  #include <stdlib.h>
  #include <windows.h>
  #include <conio.h>
  #include <math.h>

  
  //DS OPROMINUVACHA
  //DS oprominuvacha v osnovniy ploschini

  
  double  FOH(double psi, double ar, double lw )
  {
   
  //FH(Psi) - normalised directivity pattern of the horn for vector H
  //   H
  // ar  -  horizontal  size of  aperture of horn 0,080  m,
  // lw  -  lyambda obt., m
  //teta is psi0H
    
  return  fabs(0.5*(1+cos(psi))*cos(M_PI*ar*sin(psi)/lw)/(1-pow((2*ar*sin(psi)/lw),2)));
  }

  //v dopomizhniy ploschini
  double  FOE(double teta, double br, double lw )
  {
  //FE(teta) -normalised directivity pattern of the horn for vector E
  //    teta in radians
  // br-vert. size aperture of aperture of horn ,  0,220  m,
  // lw  -lyambda  , m

  return fabs(0.5*(1+cos(teta))*(   sin(M_PI*br*sin(teta)/lw) / (M_PI*br*sin(teta )/lw )  ));
  }

  // example
  //psi0e=20.62*M_PI/180; //psi0dop
  //psi0h=69*M_PI/180;    //psi0osn
  //lyambda0=0.1;
  //ap = 0.08; bp = 0.22
  //  FOE(psi0dop, br,  lyambdamin );
  //  FOE(psi0dop, br,  lyambda0 );           // =0.25865
  //  FOE(psi0dop, br,  lyambdamax );

  //  FOH(psi0osn, ar,  lyambdamin );     
  //  FOH(psi0osn, ar,  lyambda0 );
  //  FOH(psi0osn, ar,  lyambdamax  ); FH1(psi0osn) = 0.37146

  
  //FE1(0.1597) = 0.79426
  //FE2(0.1647) = 0.7939
  //FE3(0.1695) = 0.79399

  //wDN2h :=2*(180/M_PI)*root[(FH2(psi) - 0.707),psi, 0, 1.5]=76.06224
  //wDN2e :=2*(180/M_PI)*root[(FE2(psi) - 0.707),psi, 0.000001, 1.5]=22.9172
  //FH(psi,ar,lyambda)=0.707 ->psi1, psi2 Width of directivity pattern
  //FE(teta,br,lyambda)=0.707-> Width of directivity pattern

  Sunday, September 27, 2015 7:58 PM
• 

  

  0

  Sign in to vote

  Simply rootfinder

  #include <stdio.h>
  #include <conio.h>
  #include <stdlib.h>
  #include <math.h>

  double function1(double X)
  {
  return X*X-4;

  }

  double rootsvect[100];

   void root(double (*f)(double), double step, double argmin , double argmax, double ksi,  int *nroots) 
  {
   double arg;
  // finds roots (x0,x1,xn) of f(x)=0
  int n =0;
  arg=argmin;
  while(arg<argmax)
   {
    if( (fabs(f(arg) ))<=ksi)
    {
   
    rootsvect[n]=arg;
    n++;
   
    }
   arg=arg+step;
   }

  *nroots=n;
  return;
  }

  int main (void)
  {
  printf("\n Finding roots of function1(x)=0 " );
  int nr,i;
  double step=1e-4;
  printf("\n Calculating " );
  root(function1, step, -3 , 3,step/10,  &nr);
   
  printf ("\n Found %d root(s) ",nr);

  for (i=0;i<nr;i++)
  {
  printf ("\n   %d -root(s) is %le ",i,rootsvect[i] );
  }

  
  getch();
  return 0;

  }

  How to rebuilt it to f(x)=x*x-4.261  for quick calculating of zeros with need treshold   ?

  use           -ksi<=f(x)<=ksi   if((condition 1)&&(condition 2))          /диапазон , исключить ложные корни/
  

  • Edited by USERPC01 Sunday, September 27, 2015 9:07 PM

  Sunday, September 27, 2015 9:04 PM
• 

  

  0

  Sign in to vote

   Examples from Internet  for secant method

  #include <stdio.h>  //#include <cstdio.h
  #include <conio.h>  //
  #include <math.h>   //#include <cmath.h>
  float f(float x)
  {
      return(x*x*x-4.025); // f(x)= x^3-4
  }
  int  main(void)
  {
      float a,b,c,d,e;
      int count=1,n;
      printf("\n\nEnter the values of a and b:\n"); //(a,b) must contain the solution.
      scanf("%f%f",&a,&b);
      printf("Enter the values of allowed error and maximun number of iterations:\n");
      scanf("%f %d",&e,&n);
      do
      {
          if(f(a)==f(b))
          {
              printf("\nSolution cannot be found as the values of a and b are same.\n");
          return 0;
          }
          c=(a*f(b)-b*f(a))/(f(b)-f(a));
          a=b;
          b=c;
          printf("Iteration No-%d    x=%f\n",count,c);
          count++;
          if(count==n)
          {
          break;
          }
      } while(fabs(f(c))>e);
      printf("\n The required solution is %f\n",c);
   return 0;
  }

  #include <stdio.h>
  #include <conio.h>
  #include <math.h>
  #define ESP 0.0001
  #define F(x) (x)*(x) - 4*(x) - 10
  int  main(void)
  {
    float x1,x2,x3,f1,f2,t;
    //clrscr();
    printf("\nEnter the value of x1: ");
    scanf("%f",&x1);
    printf("\nEnter the value of x2: ");
    scanf("%f",&x2);
    printf("\n______________________________________________\n");
    printf("\n    x1\t  x2\t  x3\t     f(x1)\t f(x2)");
    printf("\n______________________________________________\n");
    do
    {
    f1=F(x1);
    f2=F(x2);
    x3=x2-((f2*(x2-x1))/(f2-f1));
    printf("\n%f   %f   %f   %f   %f",x1,x2,x3,f1,f2);
    x1=x2;
    x2=x3;
    if(f2<0)
      t=fabs(f2);
    else
      t=f2;
    }while(t>ESP);
  printf("\n______________________________________________\n");
  printf("\n\nApp.root = %f",x3);
  getch();
  return 0;
  }

  /*
   OUT PUT
  ---------

  
  Enter the value of x1: 4

  Enter the value of x2: 2

  ___________________________________________________________

      x1        x2         x3         f(x1)       f(x2)
  ___________________________________________________________

  4.000000   2.000000   9.000000   -10.000000   -14.000000
  2.000000   9.000000   4.000000   -14.000000    35.000000
  9.000000   4.000000   5.111111    35.000000   -10.000000
  4.000000   5.111111   5.956522   -10.000000   -4.320987
  5.111111   5.956522   5.722488   -4.320987     1.654063
  5.956522   5.722488   5.741121    1.654063    -0.143084
  5.722488   5.741121   5.741659   -0.143084    -0.004015
  5.741121   5.741659   5.741657   -0.004015     0.000010
  ___________________________________________________________

  
  App.root = 5.741657

  */

  Newton-Raphson    method

  #include <conio.h>
  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>

  int user_power,i=0,cnt=0,flag=0;
  int coef[10]={0};
  float x1=0,x2=0,t=0;
  float fx1=0,fdx1=0;

  int main(void)
  {

      //clrscr();

      printf("\n\n\t\t\t PROGRAM FOR NEWTON RAPHSON GENERAL");

      printf("\n\n\n\tENTER THE TOTAL NO. OF POWER:::: ");
      scanf("%d",&user_power);

      for(i=0;i<=user_power;i++)
      {
          printf("\n\t x^%d::",i);
          scanf("%d",&coef[i]);
      }

      printf("\n");

      printf("\n\t THE POLYNOMIAL IS ::: ");
      for(i=user_power;i>=0;i--)//printing coeff.
      {
          printf(" %dx^%d",coef[i],i);
      }

      printf("\n\tINTIAL X1---->");
      scanf("%f",&x1);

      printf("\n ******************************************************");
      printf("\n ITERATION    X1    FX1    F'X1  ");
      printf("\n **********************************************************");

      do
      {
              cnt++;
              fx1=fdx1=0;
              for(i=user_power;i>=1;i--)
              {
                  fx1+=coef[i] * (pow(x1,i)) ;
              }
              fx1+=coef[0];
              for(i=user_power;i>=0;i--)
              {
                  fdx1+=coef[i]* (i*pow(x1,(i-1)));
              }
              t=x2;
              x2=(x1-(fx1/fdx1));

              x1=x2;

              printf("\n %d         %.3f  %.3f  %.3f ",cnt,x2,fx1,fdx1);

      }while((fabs(t - x1))>=0.0001);
      printf("\n\t THE ROOT OF EQUATION IS %f",x2);
      getch();
      return 0;
  }

  /*******************************OUTPUT***********************************/

  /*
              PROGRAM FOR NEWTON RAPHSON GENERAL

  
      ENTER THE TOTAL NO. OF POWER:::: 3

       x^0::-3

       x^1::-1

       x^2::0

       x^3::1

  
       THE POLYNOMIAL IS :::  1x^3 0x^2 -1x^1 -3x^0

       INTIAL X1---->3

   **************************************
   ITERATION    X1    FX1    F'X1
   **************************************
   1         2.192  21.000 26.000
   2         1.794  5.344  13.419
   3         1.681  0.980  8.656
   4         1.672  0.068  7.475
   5         1.672  0.000  7.384
   **************************************

       THE ROOT OF EQUATION IS 1.671700
  */

  
  

  • Edited by USERPC01 Saturday, October 3, 2015 6:40 PM

  Saturday, October 3, 2015 6:39 PM
• 

  

  0

  Sign in to vote

  http://www.mcs.csueastbay.edu/~malek/TeX/root.pdf

  Saturday, October 3, 2015 6:44 PM
• 

  

  0

  Sign in to vote

                

            

  <meta content="6.42" itemprop="price" /> A First Course in Numerical Analysis: Second Edition (Dover Books on Mathematics)      Anthony Ralston; Philip Rabinowitz ,1978     http://www.abebooks.co.uk/book-search/title/a-first-course-in-numerical-analysis/author/ralston-a-rabinowitz-philip/ 

  
  

  • Edited by USERPC01 Saturday, October 3, 2015 6:55 PM

  Saturday, October 3, 2015 6:52 PM
• 

  

  0

  Sign in to vote

  #include  <stdio.h>
  #include  <conio.h>
  #include  <stdlib.h>
  #include  <math.h>
  //Ridder's rule ;from site  exponenta.ru Численные методы на Mathcad'е Ю. Ю. Тарасевич ,

  
  double f( double x)
  {
  //x=[1,4]
  //return x*x-log(x)-2*cos(x)-1;
  //x=[-4,4]
  return  (x+0.25)*(x+0.25)- 1  ;
  }

  double sign(double x)
  {
   if (x==0) return 0;
   if (x<1)  return -1;
   return 1;
  }

  void ridder1(double a, double b, double eps_user)
  {
  double   eps,x,c, i_x[200] ;
  int i;
  //fa=f(a);
  //fb=f(b);
  i=0;
  eps=fabs(a-b);
  if (f(a)*f(b)<0)
   {
    while(eps>eps_user)
    {
     i=i+1;
     c=0.5*(a+b);
     //Finding  poinds of  intersection  g(x) and OX 
     x=c+(c-a)*
  (sign(f(a)-f(b))*f(c))/pow(( f(c)*f(c)-(f(a)*f(b))),0.5 );

     if(f(a)*f(x)<0)
      {
      // left root localisation  interval
      a=a;
      b=x;
      }
     if(f(x)*f(b)<0)
      {
      // right root localisation  interval
      a=x;
      b=b;
      }

     i_x[i]=x;
     eps=fabs(i_x[i]-i_x[i-1]);

     //Creating of roots matrix

     //Usermatrix[i][0]=eps;
     //Usermatrix[i][1]=x;
     //Usermatrix[i][2]=f(x);

     // printf("\n Eps[%d]=%le ",i,eps);
     // printf(";  x[%d]=%le ",i,x);
    //  printf(";  f(x[%d])=%le  ",i,f(x));
     }
      printf("\n Eps[%d]=%le ;  x[%d]=%le ;  f(x[%d])=%le",i,eps,i,x,i,f(x));
     
      
    
    // printf("\n  Iteration with |eps|<|eps_user| i=%d   " ,i );
   
    return ;
   }

  
   printf("\n Error in localisation interval");
   return;
  }

  
  void ridder(double a, double b, double eps_user)
  {
  if ((a<0)&&(b>0))
  {
  ridder1( a,   0,   eps_user);
  ridder1( 0,   b,   eps_user);
  }

  if (((a<=0)&&(b<=0))||((a>=0)&&(b>=0)))
  {
  ridder1( a,   b,   eps_user);
  }

  return;
  }

  int  main(void)
  {
  double  a,b,Epsilon;
  a=-4;
  b=4;
  Epsilon=0.0000001;
   
  ridder( a,  b,   Epsilon);

  
  getch();
  return 0;
  }

  
  

  • Edited by USERPC01 Friday, October 9, 2015 12:26 AM

  Friday, October 9, 2015 12:25 AM
• 

  

  0

  Sign in to vote

  #include<stdio.h>

  #include<math.h>

  #define I 2

  float f(float x)

  {

      return cos(x) - x*exp(x);

  }

  main ()

  {

      int i, itr, maxmitr;

      float x[4], li, di, mu, s, l, allerr;

      printf("\nEnter the three initial guesses:\n");

      for (i=I-2; i<3; i++)

          scanf("%f", &x[i]);

      printf("Enter allowed error and maximum iterations:\n");

      scanf("%f %d", &allerr, &maxmitr);

      for (itr=1; itr<=maxmitr; itr++)

      {

          li = (x[I] - x[I-1]) / (x[I-1] - x[I-2]);

          di = (x[I] - x[I-2]) / (x[I-1] - x[I-2]);

          mu = f(x[I-2])*li*li - f(x[I-1])*di*di + f(x[I])*(di+li);

          s = sqrt ((mu*mu - 4*f(x[I])*di*li*(f(x[I-2])*li - f(x[I-1])*di + f(x[I]))));

          if (mu < 0)

              l = (2*f(x[I])*di)/(-mu+s);

          else

              l = (2*f(x[I])*di)/(-mu-s);

          x[I+1] = x[I] + l*(x[I] - x[I-1]);

          printf("At iteration no. %3d, x = %7.5f\n", itr, x[I+1]);

          if (fabs (x[I+1] - x[I]) < allerr)

          {

          printf("After %3d iterations, the required root is %6.4f\n", itr, x[I+1]);

          return 0;

          }

          for (i=I-2; i<3; i++)

              x[i] = x[i+1];

      }

  printf("The required solution does not converge or iterations are insufficient\n");

  return 1;

  }

  Friday, October 9, 2015 12:27 AM
• 

  

  0

  Sign in to vote

  #include<stdio.h>

  #include<conio.h>

  float f(float x)

  {

      return(1/(1+x));

  }

  int main()

  {

      int i,n;

      float x0,xn,h,y[20],so,se,ans,x[20];

      printf("\n Enter values of x0,xn,h: ");

      scanf("%f%f%f",&x0,&xn,&h);

      n=(xn-x0)/h;

      if(n%2==1)

      {

          n=n+1;

      }

      h=(xn-x0)/n;

      printf("\n Refined value of n and h are:%d %f\n",n,h);

      printf("\n Y values: \n");

      for(i=0; i<=n; i++)

      {

          x[i]=x0+i*h;

          y[i]=f(x[i]);

          printf("\n %f\n",y[i]);

      }

      so=0;

      se=0;

      for(i=1; i<n; i++)

      {

          if(i%2==1)

          {

              so=so+y[i];

          }

          else

          {

              se=se+y[i];

          }

   

      }

      ans=h/3*(y[0]+y[n]+4*so+2*se);

      printf("\n Final integration is %f",ans);

   

      getch();

  }

  Friday, October 9, 2015 12:27 AM
• 

  

  0

  Sign in to vote

  Example for Brent 's algorithm from internet:

  # include <cstdlib>
  # include <iostream>
  # include <iomanip>
  # include <cmath>
  # include <ctime>
  # include <string>
  //# include "brent.hpp"
  # include "brent.cpp"
  using namespace std;
  using namespace brent;

  const double check_tolerance(1e-15);

  void test_zero_all ( );
  void test_zero_rc_all ( );
  void test_local_min_all ( );
  void test_local_min_rc_all ( );
  void test_glomin_all ( );

  
  void test_zero_one ( double a, double b, double t, double f ( double x ),
    string title );
  void test_zero_rc_one ( double a, double b, double t, double f ( double x ),
    string title );
  template <class T>
  void test_local_min_one ( double a, double b, double t, T f,
    string title );
  void test_local_min_rc_one ( double a, double b, double t,
    double f ( double x ), string title );
  void test_glomin_one ( double a, double b, double c, double m, double e,
    double t, double f ( double x ), string title );

  double f_01 ( double x );
  double f_02 ( double x );
  double f_03 ( double x );
  double f_04 ( double x );
  double f_05 ( double x );
  double g_01 ( double x );
  double g_02 ( double x );
  double g_03 ( double x );
  double g_04 ( double x );
  double g_05 ( double x );
  double h_01 ( double x );
  double h_02 ( double x );
  double h_03 ( double x );
  double h_04 ( double x );
  double h_05 ( double x );

  //****************************************************************************80

  int main ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    MAIN is the main program for BRENT_PRB.
  //
  //  Discussion:
  //
  //    BRENT_PRB tests the BRENT library.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    John Burkardt
  //
  {
    timestamp ( );
    cout << "\n";
    cout << "BRENT_PRB\n";
    cout << "  C++ version\n";
    cout << "  Test the BRENT library.\n";

    test_zero_all ( );
    test_zero_rc_all ( );
    test_local_min_all ( );
    test_local_min_rc_all ( );
    test_glomin_all ( );
  //
  //  Terminate.
  //
    cout << "\n";
    cout << "BRENT_PRB\n";
    cout << "  Normal end of execution.\n";
    cout << "\n";
    timestamp ( );

    return 0;
  }
  //****************************************************************************80

  void test_zero_all ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_ZERO_ALL tests Brent's zero finding routine on all test functions.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    John Burkardt
  //
  {
    double a;
    double b;
    double t;

    cout << "\n";
    cout << "TEST_ZERO_ALL\n";
    cout << "  Test the Brent ZERO routine, which seeks\n";
    cout << "  a root of a function F(X)\n";
    cout << "  in an interval [A,B].\n";

    t = r8_epsilon ( );

    a = 1.0;
    b = 2.0;

    test_zero_one ( a, b, t, f_01,
      "f_01(x) = sin ( x ) - x / 2" );

    a = 0.0;
    b = 1.0;

    test_zero_one ( a, b, t, f_02,
      "f_02(x) = 2 * x - exp ( - x )" );

    a = -1.0;
    b =  0.5;

    test_zero_one ( a, b, t, f_03,
      "f_03(x) = x * exp ( - x )" );

    a =  0.0001;
    b =  20.0;

    test_zero_one ( a, b, t, f_04,
      "f_04(x) = exp ( x ) - 1 / ( 100 * x * x )" );

    a = -5.0;
    b =  2.0;

    test_zero_one ( a, b, t, f_05,
      "f_05(x) = (x+3) * (x-1) * (x-1)" );

    return;
  }
  //****************************************************************************80

  void test_zero_rc_all ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_ZERO_RC_ALL tests ZERO_RC on all test functions.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 October 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  {
    double a;
    double b;
    double t;

    cout << "\n";
    cout << "TEST_ZERO_RC_ALL\n";
    cout << "  Test the ZERO_RC routine, which seeks\n";
    cout << "  a root of a function F(X)\n";
    cout << "  in an interval [A,B].\n";

    t = r8_epsilon ( );

    a = 1.0;
    b = 2.0;

    test_zero_rc_one ( a, b, t, f_01,
      "f_01(x) = sin ( x ) - x / 2" );

    a = 0.0;
    b = 1.0;

    test_zero_rc_one ( a, b, t, f_02,
      "f_02(x) = 2 * x - exp ( - x )" );

    a = -1.0;
    b =  0.5;

    test_zero_rc_one ( a, b, t, f_03,
      "f_03(x) = x * exp ( - x )" );

    a =  0.0001;
    b =  20.0;

    test_zero_rc_one ( a, b, t, f_04,
      "f_04(x) = exp ( x ) - 1 / ( 100 * x * x )" );

    a = -5.0;
    b =  2.0;

    test_zero_rc_one ( a, b, t, f_05,
      "f_05(x) = (x+3) * (x-1) * (x-1)" );

    return;
  }
  //****************************************************************************80

  void test_local_min_all ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_LOCAL_MIN_ALL tests Brent"s local minimizer on all test functions.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    John Burkardt
  //
  {
    double a;
    double b;
    double t;

    cout << "\n";
    cout << "TEST_LOCAL_MIN_ALL\n";
    cout << "  Test the Brent LOCAL_MIN routine, which seeks\n";
    cout << "  a local minimizer of a function F(X)\n";
    cout << "  in an interval [A,B].\n";

    t = r8_epsilon ( );

    a = 0.0;
    b = 3.141592653589793;

    test_local_min_one ( a, b, t, g_01,
      "g_01(x) = ( x - 2 ) * ( x - 2 ) + 1" );

    a = 0.0;
    b = 1.0;

    test_local_min_one ( a, b, t, g_02,
      "g_02(x) = x * x + exp ( - x )" );

    a = -2.0;
    b =  2.0;

  // coefficients in ascending order, starting with scalar term:
    const int degree(4);
    double coefflist[1+degree] = {3, 1, 2, 0, 1};
    Poly quartic(coefflist, degree);

    test_local_min_one ( a, b, t, quartic,
      "g_03(x) = x^4 + 2x^2 + x + 3" );

    a =  0.0001;
    b =  1.0;

    test_local_min_one ( a, b, t, g_04,
      "g_04(x) = exp ( x ) + 1 / ( 100 x )" );

    a =  0.0002;
    b = 2.0;

    test_local_min_one ( a, b, t, g_05,
      "g_05(x) = exp ( x ) - 2x + 1/(100x) - 1/(1000000x^2)" );

    return;
  }
  //****************************************************************************80

  void test_local_min_rc_all ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_LOCAL_MIN_RC_ALL tests LOCAL_MIN_RC on all test functions.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    16 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  {
    double a;
    double b;
    double t;

    cout << "\n";
    cout << "TEST_LOCAL_MIN_RC_ALL\n";
    cout << "  Test the reverse communication version of\n";
    cout << "  the Brent LOCAL_MIN routine, which seeks\n";
    cout << "  a local minimizer of a function F(X)\n";
    cout << "  in an interval [A,B].\n";

    t = 10.0 * sqrt ( r8_epsilon ( ) );

    a = 0.0;
    b = 3.141592653589793;

    test_local_min_rc_one ( a, b, t, g_01,
      "g_01(x) = ( x - 2 ) * ( x - 2 ) + 1" );

    a = 0.0;
    b = 1.0;

    test_local_min_rc_one ( a, b, t, g_02,
      "g_02(x) = x * x + exp ( - x )" );

    a = -2.0;
    b =  2.0;

    test_local_min_rc_one ( a, b, t, g_03,
      "g_03(x) = x^4 + 2x^2 + x + 3" );

    a =  0.0001;
    b =  1.0;

    test_local_min_rc_one ( a, b, t, g_04,
      "g_04(x) = exp ( x ) + 1 / ( 100 x )" );

    a =  0.0002;
    b = 2.0;

    test_local_min_rc_one ( a, b, t, g_05,
      "g_05(x) = exp ( x ) - 2x + 1/(100x) - 1/(1000000x^2)" );

    return;
  }
  //****************************************************************************80

  void test_glomin_all ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_GLOMIN_ALL tests the Brent global minimizer on all test functions.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  {
    double a;
    double b;
    double c;
    double e;
    double m;
    double t;

    cout << "\n";
    cout << "TEST_GLOMIN_ALL\n";
    cout << "  Test the Brent GLOMIN routine, which seeks\n";
    cout << "  a global minimizer of a function F(X)\n";
    cout << "  in an interval [A,B],\n";
    cout << "  given some upper bound M \n";
    cout << "  for the second derivative of F.\n";

    e = sqrt ( r8_epsilon ( ) );
    t = sqrt ( r8_epsilon ( ) );

    a = 7.0;
    b = 9.0;
    c = ( a + b ) / 2.0;
    m = 0.0;

    test_glomin_one ( a, b, c, m, e, t, h_01,
      "h_01(x) = 2 - x, M = 0" );

    a = 7.0;
    b = 9.0;
    c = ( a + b ) / 2.0;
    m = 100.0;

    test_glomin_one ( a, b, c, m, e, t, h_01,
      "h_01(x) = 2 - x, M = 100" );

    a = -1.0;
    b = +2.0;
    c = ( a + b ) / 2.0;
    m = 2.0;

    test_glomin_one ( a, b, c, m, e, t, h_02,
      "h_02(x) = x * x, M = 2" );

    a = -1.0;
    b = +2.0;
    c = ( a + b ) / 2.0;
    m = 2.1;

    test_glomin_one ( a, b, c, m, e, t, h_02,
      "h_02(x) = x * x, M = 2.1" );

    a = -0.5;
    b =  +2.0;
    c = ( a + b ) / 2.0;
    m = 14.0;

    test_glomin_one ( a, b, c, m, e, t, h_03,
      "h_03(x) = x^3 + x^2, M = 14" );

    a = -0.5;
    b =  +2.0;
    c = ( a + b ) / 2.0;
    m = 28.0;

    test_glomin_one ( a, b, c, m, e, t, h_03,
      "h_03(x) = x^3 + x^2, M = 28" );

    a = -10.0;
    b = +10.0;
    c = ( a + b ) / 2.0;
    m = 72.0;

    test_glomin_one ( a, b, c, m, e, t, h_04,
      "h_04(x) = ( x + sin(x) ) * exp(-x*x), M = 72" );

    a = -10.0;
    b = +10.0;
    c = ( a + b ) / 2.0;
    m = 72.0;

    test_glomin_one ( a, b, c, m, e, t, h_05,
      "h_05(x) = ( x - sin(x) ) * exp(-x*x), M = 72" );

    return;
  }
  //****************************************************************************80

  void test_zero_one ( double a, double b, double t, double f ( double x ),
    string title )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_ZERO_ONE tests Brent's zero finding routine on one test function.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    13 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double A, B, the two endpoints of the change of sign
  //    interval.
  //
  //    Input, double T, a positive error tolerance.
  //
  //    Input, double F ( double x ), the name of a user-supplied
  //    function which evaluates the function whose zero is being sought.
  //
  //    Input, string TITLE, a title for the problem.
  //
  {
    double fa;
    double fb;
    double fz;
    double z;

    z = zero ( a, b, t, f );
    fz = f ( z );
    fa = f ( a );
    fb = f ( b );

    cout << "\n";
    cout << "  " << title << "\n";
    cout << "\n";
    cout << "           A                 Z             B\n";
    cout << "         F(A)              F(Z)          F(B)\n";
    cout << "  " << setw(14) << a
         << "  " << setw(14) << z
         << "  " << setw(14) << b << "\n";
    cout << "  " << setw(14) << fa
         << "  " << setw(14) << fz
         << "  " << setw(14) << fb << "\n";

    if (abs(fz) > check_tolerance) {
      cerr << "*** error ***" << endl;
      cerr << "fz " << fz
          << " exceeds check_tolerance " << check_tolerance << endl;
      exit(1);
    }

    return;
  }
  //****************************************************************************80

  void test_zero_rc_one ( double a, double b, double t, double f ( double x ),
    string title )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_ZERO_RC_ONE tests ZERO_RC on one test function.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 October 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double A, B, the two endpoints of the change of sign
  //    interval.
  //
  //    Input, double MACHEP, an estimate for the relative machine
  //    precision.
  //
  //    Input, double T, a positive error tolerance.
  //
  //    Input, double F ( double x ), the name of a user-supplied
  //    function which evaluates the function whose zero is being sought.
  //
  //    Input, string TITLE, a title for the problem.
  //
  {
    double arg;
    int status;
    double value;

    cout << "\n";
    cout << "  " << title << "\n";
    cout << "\n";
    cout << "    STATUS      X               F(X)\n";
    cout << "\n";

    status = 0;

    for ( ; ; )
    {
      zero_rc ( a, b, t, arg, status, value );

      if ( status < 0 )
      {
        cout << "\n";
        cout << "  ZERO_RC returned an error flag!\n";
        break;
      }

      value = f ( arg );

      cout << "  " << setw(8) << status
           << "  " << setw(14) << arg
           << "  " << setw(14) << value << "\n";

      if ( status == 0 )
      {
        break;
      }
    }
    if (abs(value) > check_tolerance) {
      cerr << "*** error ***" << endl;
      cerr << "final value " << value
          << " exceeds check_tolerance " << check_tolerance << endl;
      exit(1);
    }

    return;
  }

  //****************************************************************************80

  template <class T>
  void test_local_min_one ( double a, double b, double t, T f,
    string title )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_LOCAL_MIN_ONE tests Brent's local minimizer on one test function.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the interval.
  //
  //    Input, double T, a positive absolute error tolerance.
  //
  //    Input, double F ( double x ), the name of a user-supplied
  //    function, whose local minimum is being sought.
  //
  //    Input, string TITLE, a title for the problem.
  //
  {
    double fa;
    double fb;
    double fx;
    double x;

    fx = local_min ( a, b, t, f, x );
    fa = f ( a );
    fb = f ( b );

    cout << "\n";
    cout << "  " << title << "\n";
    cout << "\n";
    cout << "           A                 X             B\n";
    cout << "         F(A)              F(X)          F(B)\n";
    cout << "  " << setw(14) << a
         << "  " << setw(14) << x
         << "  " << setw(14) << b << "\n";
    cout << "  " << setw(14) << fa
         << "  " << setw(14) << fx
         << "  " << setw(14) << fb << "\n";

    return;
  }
  //****************************************************************************80

  void test_local_min_rc_one ( double a, double b, double t,
    double f ( double x ), string title )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_LOCAL_MIN_RC_ONE tests LOCAL_MIN_RC on one test function.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    16 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the interval.
  //
  //    Input, double T, a positive absolute error tolerance.
  //
  //    Input, double F ( double x ), the name of a user-supplied
  //    function, whose local minimum is being sought.
  //
  //    Input, string TITLE, a title for the problem.
  //
  {
    double a2;
    double arg;
    double b2;
    int status;
    int step;
    double value;

    cout << "\n";
    cout << "  " << title << "\n";
    cout << "\n";
    cout << "  Step      X                          F(X)\n";
    cout << "\n";
    step = 0;

    arg = a;
    value = f ( arg );
    cout << "  " << setw(4) << step
         << "  " << setprecision(16) << setw(24) << arg
         << "  " << setprecision(16) << setw(24) << value << "\n";

    arg = b;
    value = f ( arg );
    cout << "  " << setw(4) << step
         << "  " << setprecision(16) << setw(24) << arg
         << "  " << setprecision(16) << setw(24) << value << "\n";

    a2 = a;
    b2 = b;
    status = 0;

    for ( ; ; )
    {
      arg = local_min_rc ( a2, b2, status, value );

      if ( status < 0 )
      {
        cout << "\n";
        cout << "TEST_LOCAL_MIN_RC_ONE - Fatal error!\n";
        cout << "  LOCAL_MIN_RC returned negative status.\n";
        break;
      }

      value = f ( arg );

      step = step + 1;
      cout << "  " << setw(4) << step
           << "  " << setprecision(16) << setw(24) << arg
           << "  " << setprecision(16) << setw(24) << value << "\n";

      if ( 50 < step )
      {
        cout << "\n";
        cout << "TEST_LOCAL_MIN_RC_ONE - Fatal error!\n";
        cout << "  Too many steps!\n";
        break;
       }
      if ( status == 0 )
      {
        break;
      }
    }

    return;
  }
  //****************************************************************************80

  void test_glomin_one ( double a, double b, double c, double m,
    double e, double t, double f ( double x ), string title )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TEST_GLOMIN_ONE tests the Brent global minimizer on one test function.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the interval.
  //
  //    Input, double C, an initial guess for the global
  //    minimizer.  If no good guess is known, C = A or B is acceptable.
  //
  //    Input, double M, the bound on the second derivative.
  //
  //    Input, double E, a positive tolerance, a bound for the
  //    absolute error in the evaluation of F(X) for any X in [A,B].
  //
  //    Input, double T, a positive absolute error tolerance.
  //
  //    Input, double F ( double x ), the name of a user-supplied
  //    function whose global minimum is being sought.
  //
  //    Input, string TITLE, a title for the problem.
  //
  {
    double fa;
    double fb;
    double fx;
    double x;

    fx = glomin ( a, b, c, m, e, t, f, x );
    fa = f ( a );
    fb = f ( b );

    cout << "\n";
    cout << "  " << title << "\n";
    cout << "\n";
    cout << "           A                 X             B\n";
    cout << "         F(A)              F(X)          F(B)\n";
    cout << "  " << setprecision(6) << setw(14) << a
         << "  " << setw(14) << x
         << "  " << setw(14) << b << "\n";
    cout << "  " << setw(14) << fa
         << "  " << setw(14) << fx
         << "  " << setw(14) << fb << "\n";

    return;
  }
  //****************************************************************************80

  double f_01 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    F_01 evaluates sin ( x ) - x / 2.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    13 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double F_01, the value of the function at X.
  //
  {
    double value;

    value = sin ( x ) - 0.5 * x;

    return value;
  }
  //****************************************************************************80

  double f_02 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    F_02 evaluates 2*x-exp(-x).
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    13 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double F_02, the value of the function at X.
  //
  {
    double value;

    value = 2.0 * x - exp ( - x );

    return value;
  }
  //****************************************************************************80

  double f_03 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    F_03 evaluates x*exp(-x).
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    13 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double F_03, the value of the function at X.
  //
  {
    double value;

    value = x * exp ( - x );

    return value;
  }
  //****************************************************************************80

  double f_04 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    F_04 evaluates exp(x) - 1 / (100*x*x).
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    13 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double F_04, the value of the function at X.
  //
  {
    double value;

    value = exp ( x ) - 1.0 / 100.0 / x / x;

    return value;
  }
  //****************************************************************************80

  double f_05 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    F_05 evaluates (x+3)*(x-1)*(x-1).
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    13 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double F_05, the value of the function at X.
  //
  {
    double value;

    value = ( x + 3.0 ) * ( x - 1.0 ) * ( x - 1.0 );

    return value;
  }
  //****************************************************************************80

  double g_01 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    G_01 evaluates (x-2)^2 + 1.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double G_01, the value of the function at X.
  //
  {
    double value;

    value = ( x - 2.0 ) * ( x - 2.0 ) + 1.0;

    return value;
  }
  //****************************************************************************80

  double g_02 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    G_02 evaluates x^2 + exp ( - x ).
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double G_02, the value of the function at X.
  //
  {
    double value;

    value = x * x + exp ( - x );

    return value;
  }
  //****************************************************************************80

  double g_03 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    G_03 evaluates x^4+2x^2+x+3.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double G_03, the value of the function at X.
  //
  {
    double value;

    value = ( ( x * x + 2.0 ) * x + 1.0 ) * x + 3.0;

    return value;
  }
  //****************************************************************************80

  double g_04 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    G_04 evaluates exp(x)+1/(100X)
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double G_04, the value of the function at X.
  //
  {
    double value;

    value = exp ( x ) + 0.01 / x;

    return value;
  }
  //****************************************************************************80

  double g_05 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    G_05 evaluates exp(x) - 2x + 1/(100x) - 1/(1000000x^2)
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double G_05, the value of the function at X.
  //
  {
    double value;

    value = exp ( x ) - 2.0 * x + 0.01 / x - 0.000001 / x / x;

    return value;
  }
  //****************************************************************************80

  double h_01 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    H_01 evaluates 2 - x.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double H_01, the value of the function at X.
  //
  {
    double value;

    value = 2.0 - x;

    return value;
  }
  //****************************************************************************80

  double h_02 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    H_02 evaluates x^2.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double H_02, the value of the function at X.
  //
  {
    double value;

    value = x * x;

    return value;
  }
  //****************************************************************************80

  double h_03 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    H_03 evaluates x^3+x^2.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double H_03, the value of the function at X.
  //
  {
    double value;

    value = x * x * ( x + 1.0 );

    return value;
  }
  //****************************************************************************80

  double h_04 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    H_04 evaluates ( x + sin ( x ) ) * exp ( - x * x ).
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double H_04, the value of the function at X.
  //
  {
    double value;

    value = ( x + sin ( x ) ) * exp ( - x * x );

    return value;
  }
  //****************************************************************************80

  double h_05 ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    H_05 evaluates ( x - sin ( x ) ) * exp ( - x * x ).
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    14 April 2008
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the point at which F is to be evaluated.
  //
  //    Output, double H_05, the value of the function at X.
  //
  {
    double value;

    value = ( x - sin ( x ) ) * exp ( - x * x );

    return value;
  }

   

  Friday, October 9, 2015 12:35 AM
• 

  

  0

  Sign in to vote

  //brent.hpp

  #include <vector>

  namespace brent
  {

  
  class func_base
  {

  public:
   
    virtual double operator() (double) = 0;

  };

  
  class monicPoly :
   
    public func_base
   {

   public:
    std::vector<double> coeff;
    virtual double operator() (double x);

  // constructors:

    monicPoly(const size_t degree)
     : coeff(degree) {}

    monicPoly(const std::vector<double>& v)
     : coeff(v) {}

    monicPoly(const double* c, size_t degree)
     : coeff(std::vector<double>(c, c+degree)) {}
  };

  class Poly : public func_base {
  public:
    std::vector<double> coeff;
      // a vector of size nterms i.e. 1+degree
    virtual double operator() (double x);

      // constructors:

    Poly(const size_t degree)
     : coeff(1+degree) {}
    Poly(const std::vector<double>& v)
     : coeff(v) {}
    Poly(const double* c, size_t degree)
     : coeff(std::vector<double>(c, 1+c+degree)) {}
  };

  double glomin ( double a, double b, double c, double m, double e, double t,
    func_base& f, double &x );
  double local_min ( double a, double b, double t, func_base& f,
    double &x );
  double local_min_rc ( double &a, double &b, int &status, double value );
  double r8_abs ( double x );
  double r8_epsilon ( );
  double r8_max ( double x, double y );
  double r8_sign ( double x );
  void timestamp ( );
  double zero ( double a, double b, double t, func_base& f );
  void zero_rc ( double a, double b, double t, double &arg, int &status,
    double value );

  // === simple wrapper functions
  // === for convenience and/or compatibility
  double glomin ( double a, double b, double c, double m, double e, double t,
    double f ( double x ), double &x );
  double local_min ( double a, double b, double t, double f ( double x ),
    double &x );
  double zero ( double a, double b, double t, double f ( double x ) );

  }

   

  Friday, October 9, 2015 12:41 AM
• 

  

  0

  Sign in to vote

  //brent.cpp

  #include <cstdlib>
  //#include <windows.h>
  #include <iostream>

  #include <cmath>

  #include <ctime>

  
  using namespace std;

  
  # include "brent.hpp"

  
  namespace brent{

  
  //****************************************************************************80

  double glomin ( double a, double b, double c, double m, double e, double t,
    func_base& f, double &x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    GLOMIN seeks a global minimum of a function F(X) in an interval [A,B].
  //
  //  Discussion:
  //
  //    This function assumes that F(X) is twice continuously differentiable
  //    over [A,B] and that F''(X) <= M for all X in [A,B].
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    Original FORTRAN77 version by Richard Brent.
  //    C++ version by John Burkardt.
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the interval.
  //    It must be the case that A < B.
  //
  //    Input, double C, an initial guess for the global
  //    minimizer.  If no good guess is known, C = A or B is acceptable.
  //
  //    Input, double M, the bound on the second derivative.
  //
  //    Input, double E, a positive tolerance, a bound for the
  //    absolute error in the evaluation of F(X) for any X in [A,B].
  //
  //    Input, double T, a positive error tolerance.
  //
  //    Input, func_base& F, a user-supplied c++ functor whose
  //    global minimum is being sought.  The input and output
  //    of F() are of type double.
  //
  //    Output, double &X, the estimated value of the abscissa
  //    for which F attains its global minimum value in [A,B].
  //
  //    Output, double GLOMIN, the value F(X).
  //
  {
    double a0;
    double a2;
    double a3;
    double d0;
    double d1;
    double d2;
    double h;
    int k;
    double m2;
    double macheps;
    double p;
    double q;
    double qs;
    double r;
    double s;
    double sc;
    double y;
    double y0;
    double y1;
    double y2;
    double y3;
    double yb;
    double z0;
    double z1;
    double z2;
   
    a0 = b;
    x = a0;
    a2 = a;
    y0 = f ( b );
    yb = y0;
    y2 = f ( a );
    y = y2;

    if ( y0 < y )
    {
      y = y0;
    }
    else
    {
      x = a;
    }

    if ( m <= 0.0 || b <= a )
    {
      return y;
    }

    macheps = r8_epsilon ( );

    m2 = 0.5 * ( 1.0 + 16.0 * macheps ) * m;

    if ( c <= a || b <= c )
    {
      sc = 0.5 * ( a + b );
    }
    else
    {
      sc = c;
    }

    y1 = f ( sc );
    k = 3;
    d0 = a2 - sc;
    h = 9.0 / 11.0;

    if ( y1 < y )
    {
      x = sc;
      y = y1;
    }
  //
  //  Loop.
  //
    for ( ; ; )
    {
      d1 = a2 - a0;
      d2 = sc - a0;
      z2 = b - a2;
      z0 = y2 - y1;
      z1 = y2 - y0;
      r = d1 * d1 * z0 - d0 * d0 * z1;
      p = r;
      qs = 2.0 * ( d0 * z1 - d1 * z0 );
      q = qs;

      if ( k < 1000000 || y2 <= y )
      {
        for ( ; ; )
        {
          if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
            z2 * m2 * r * ( z2 * q - r ) )
          {
            a3 = a2 + r / q;
            y3 = f ( a3 );

            if ( y3 < y )
            {
              x = a3;
              y = y3;
            }
          }
          k = ( ( 1611 * k ) % 1048576 );
          q = 1.0;
          r = ( b - a ) * 0.00001 * ( double ) ( k );

          if ( z2 <= r )
          {
            break;
          }
        }
      }
      else
      {
        k = ( ( 1611 * k ) % 1048576 );
        q = 1.0;
        r = ( b - a ) * 0.00001 * ( double ) ( k );

        while ( r < z2 )
        {
          if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
            z2 * m2 * r * ( z2 * q - r ) )
          {
            a3 = a2 + r / q;
            y3 = f ( a3 );

            if ( y3 < y )
            {
              x = a3;
              y = y3;
            }
          }
          k = ( ( 1611 * k ) % 1048576 );
          q = 1.0;
          r = ( b - a ) * 0.00001 * ( double ) ( k );
        }
      }

      r = m2 * d0 * d1 * d2;
      s = sqrt ( ( ( y2 - y ) + t ) / m2 );
      h = 0.5 * ( 1.0 + h );
      p = h * ( p + 2.0 * r * s );
      q = q + 0.5 * qs;
      r = - 0.5 * ( d0 + ( z0 + 2.01 * e ) / ( d0 * m2 ) );

      if ( r < s || d0 < 0.0 )
      {
        r = a2 + s;
      }
      else
      {
        r = a2 + r;
      }

      if ( 0.0 < p * q )
      {
        a3 = a2 + p / q;
      }
      else
      {
        a3 = r;
      }

      for ( ; ; )
      {
        a3 = r8_max ( a3, r );

        if ( b <= a3 )
        {
          a3 = b;
          y3 = yb;
        }
        else
        {
          y3 = f ( a3 );
        }

        if ( y3 < y )
        {
          x = a3;
          y = y3;
        }

        d0 = a3 - a2;

        if ( a3 <= r )
        {
          break;
        }

        p = 2.0 * ( y2 - y3 ) / ( m * d0 );

        if ( ( 1.0 + 9.0 * macheps ) * d0 <= r8_abs ( p ) )
        {
          break;
        }

        if ( 0.5 * m2 * ( d0 * d0 + p * p ) <= ( y2 - y ) + ( y3 - y ) + 2.0 * t )
        {
          break;
        }
        a3 = 0.5 * ( a2 + a3 );
        h = 0.9 * h;
      }

      if ( b <= a3 )
      {
        break;
      }

      a0 = sc;
      sc = a2;
      a2 = a3;
      y0 = y1;
      y1 = y2;
      y2 = y3;
    }

    return y;
  }
  //****************************************************************************80

  double local_min ( double a, double b, double t, func_base& f,
    double &x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    LOCAL_MIN seeks a local minimum of a function F(X) in an interval [A,B].
  //
  //  Discussion:
  //
  //    The method used is a combination of golden section search and
  //    successive parabolic interpolation.  Convergence is never much slower
  //    than that for a Fibonacci search.  If F has a continuous second
  //    derivative which is positive at the minimum (which is not at A or
  //    B), then convergence is superlinear, and usually of the order of
  //    about 1.324....
  //
  //    The values EPS and T define a tolerance TOL = EPS * abs ( X ) + T.
  //    F is never evaluated at two points closer than TOL.
  //
  //    If F is a unimodal function and the computed values of F are always
  //    unimodal when separated by at least SQEPS * abs ( X ) + (T/3), then
  //    LOCAL_MIN approximates the abscissa of the global minimum of F on the
  //    interval [A,B] with an error less than 3*SQEPS*abs(LOCAL_MIN)+T.
  //
  //    If F is not unimodal, then LOCAL_MIN may approximate a local, but
  //    perhaps non-global, minimum to the same accuracy.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    Original FORTRAN77 version by Richard Brent.
  //    C++ version by John Burkardt.
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the interval.
  //
  //    Input, double T, a positive absolute error tolerance.
  //
  //    Input, func_base& F, a user-supplied c++ functor whose
  //    local minimum is being sought.  The input and output
  //    of F() are of type double.
  //
  //    Output, double &X, the estimated value of an abscissa
  //    for which F attains a local minimum value in [A,B].
  //
  //    Output, double LOCAL_MIN, the value F(X).
  //
  {
    double c;
    double d;
    double e;
    double eps;
    double fu;
    double fv;
    double fw;
    double fx;
    double m;
    double p;
    double q;
    double r;
    double sa;
    double sb;
    double t2;
    double tol;
    double u;
    double v;
    double w;
  //
  //  C is the square of the inverse of the golden ratio.
  //
    c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );

    eps = sqrt ( r8_epsilon ( ) );

    sa = a;
    sb = b;
    x = sa + c * ( b - a );
    w = x;
    v = w;
    e = 0.0;
    fx = f ( x );
    fw = fx;
    fv = fw;

    for ( ; ; )
    {
      m = 0.5 * ( sa + sb ) ;
      tol = eps * r8_abs ( x ) + t;
      t2 = 2.0 * tol;
  //
  //  Check the stopping criterion.
  //
      if ( r8_abs ( x - m ) <= t2 - 0.5 * ( sb - sa ) )
      {
        break;
      }
  //
  //  Fit a parabola.
  //
      r = 0.0;
      q = r;
      p = q;

      if ( tol < r8_abs ( e ) )
      {
        r = ( x - w ) * ( fx - fv );
        q = ( x - v ) * ( fx - fw );
        p = ( x - v ) * q - ( x - w ) * r;
        q = 2.0 * ( q - r );
        if ( 0.0 < q )
        {
          p = - p;
        }
        q = r8_abs ( q );
        r = e;
        e = d;
      }

      if ( r8_abs ( p ) < r8_abs ( 0.5 * q * r ) &&
           q * ( sa - x ) < p &&
           p < q * ( sb - x ) )
      {
  //
  //  Take the parabolic interpolation step.
  //
        d = p / q;
        u = x + d;
  //
  //  F must not be evaluated too close to A or B.
  //
        if ( ( u - sa ) < t2 || ( sb - u ) < t2 )
        {
          if ( x < m )
          {
            d = tol;
          }
          else
          {
            d = - tol;
          }
        }
      }
  //
  //  A golden-section step.
  //
      else
      {
        if ( x < m )
        {
          e = sb - x;
        }
        else
        {
          e = sa - x;
        }
        d = c * e;
      }
  //
  //  F must not be evaluated too close to X.
  //
      if ( tol <= r8_abs ( d ) )
      {
        u = x + d;
      }
      else if ( 0.0 < d )
      {
        u = x + tol;
      }
      else
      {
        u = x - tol;
      }

      fu = f ( u );
  //
  //  Update A, B, V, W, and X.
  //
      if ( fu <= fx )
      {
        if ( u < x )
        {
          sb = x;
        }
        else
        {
          sa = x;
        }
        v = w;
        fv = fw;
        w = x;
        fw = fx;
        x = u;
        fx = fu;
      }
      else
      {
        if ( u < x )
        {
          sa = u;
        }
        else
        {
          sb = u;
        }

        if ( fu <= fw || w == x )
        {
          v = w;
          fv = fw;
          w = u;
          fw = fu;
        }
        else if ( fu <= fv || v == x || v == w )
        {
          v = u;
          fv = fu;
        }
      }
    }
    return fx;
  }
  //****************************************************************************80

  double local_min_rc ( double &a, double &b, int &status, double value )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    LOCAL_MIN_RC seeks a minimizer of a scalar function of a scalar variable.
  //
  //  Discussion:
  //
  //    This routine seeks an approximation to the point where a function
  //    F attains a minimum on the interval (A,B).
  //
  //    The method used is a combination of golden section search and
  //    successive parabolic interpolation.  Convergence is never much
  //    slower than that for a Fibonacci search.  If F has a continuous
  //    second derivative which is positive at the minimum (which is not
  //    at A or B), then convergence is superlinear, and usually of the
  //    order of about 1.324...
  //
  //    The routine is a revised version of the Brent local minimization
  //    algorithm, using reverse communication.
  //
  //    It is worth stating explicitly that this routine will NOT be
  //    able to detect a minimizer that occurs at either initial endpoint
  //    A or B.  If this is a concern to the user, then the user must
  //    either ensure that the initial interval is larger, or to check
  //    the function value at the returned minimizer against the values
  //    at either endpoint.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //    David Kahaner, Cleve Moler, Steven Nash,
  //    Numerical Methods and Software,
  //    Prentice Hall, 1989,
  //    ISBN: 0-13-627258-4,
  //    LC: TA345.K34.
  //
  //  Parameters
  //
  //    Input/output, double &A, &B.  On input, the left and right
  //    endpoints of the initial interval.  On output, the lower and upper
  //    bounds for an interval containing the minimizer.  It is required
  //    that A < B.
  //
  //    Input/output, int &STATUS, used to communicate between
  //    the user and the routine.  The user only sets STATUS to zero on the first
  //    call, to indicate that this is a startup call.  The routine returns STATUS
  //    positive to request that the function be evaluated at ARG, or returns
  //    STATUS as 0, to indicate that the iteration is complete and that
  //    ARG is the estimated minimizer.
  //
  //    Input, double VALUE, the function value at ARG, as requested
  //    by the routine on the previous call.
  //
  //    Output, double LOCAL_MIN_RC, the currently considered point.
  //    On return with STATUS positive, the user is requested to evaluate the
  //    function at this point, and return the value in VALUE.  On return with
  //    STATUS zero, this is the routine's estimate for the function minimizer.
  //
  //  Local parameters:
  //
  //    C is the squared inverse of the golden ratio.
  //
  //    EPS is the square root of the relative machine precision.
  //
  {
    static double arg;
    static double c;
    static double d;
    static double e;
    static double eps;
    static double fu;
    static double fv;
    static double fw;
    static double fx;
    static double midpoint;
    static double p;
    static double q;
    static double r;
    static double tol;
    static double tol1;
    static double tol2;
    static double u;
    static double v;
    static double w;
    static double x;
  //
  //  STATUS (INPUT) = 0, startup.
  //
    if ( status == 0 )
    {
      if ( b <= a )
      {
        cout << "\n";
        cout << "LOCAL_MIN_RC - Fatal error!\n";
        cout << "  A < B is required, but\n";
        cout << "  A = " << a << "\n";
        cout << "  B = " << b << "\n";
        status = -1;
        exit ( 1 );
      }
      c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );

      eps = sqrt ( r8_epsilon ( ) );
      tol = r8_epsilon ( );

      v = a + c * ( b - a );
      w = v;
      x = v;
      e = 0.0;

      status = 1;
      arg = x;

      return arg;
    }
  //
  //  STATUS (INPUT) = 1, return with initial function value of FX.
  //
    else if ( status == 1 )
    {
      fx = value;
      fv = fx;
      fw = fx;
    }
  //
  //  STATUS (INPUT) = 2 or more, update the data.
  //
    else if ( 2 <= status )
    {
      fu = value;

      if ( fu <= fx )
      {
        if ( x <= u )
        {
          a = x;
        }
        else
        {
          b = x;
        }
        v = w;
        fv = fw;
        w = x;
        fw = fx;
        x = u;
        fx = fu;
      }
      else
      {
        if ( u < x )
        {
          a = u;
        }
        else
        {
          b = u;
        }

        if ( fu <= fw || w == x )
        {
          v = w;
          fv = fw;
          w = u;
          fw = fu;
        }
        else if ( fu <= fv || v == x || v == w )
        {
          v = u;
          fv = fu;
        }
      }
    }
  //
  //  Take the next step.
  //
    midpoint = 0.5 * ( a + b );
    tol1 = eps * r8_abs ( x ) + tol / 3.0;
    tol2 = 2.0 * tol1;
  //
  //  If the stopping criterion is satisfied, we can exit.
  //
    if ( r8_abs ( x - midpoint ) <= ( tol2 - 0.5 * ( b - a ) ) )
    {
      status = 0;
      return arg;
    }
  //
  //  Is golden-section necessary?
  //
    if ( r8_abs ( e ) <= tol1 )
    {
      if ( midpoint <= x )
      {
        e = a - x;
      }
      else
      {
        e = b - x;
      }
      d = c * e;
    }
  //
  //  Consider fitting a parabola.
  //
    else
    {
      r = ( x - w ) * ( fx - fv );
      q = ( x - v ) * ( fx - fw );
      p = ( x - v ) * q - ( x - w ) * r;
      q = 2.0 * ( q - r );
      if ( 0.0 < q )
      {
        p = - p;
      }
      q = r8_abs ( q );
      r = e;
      e = d;
  //
  //  Choose a golden-section step if the parabola is not advised.
  //
      if (
        ( r8_abs ( 0.5 * q * r ) <= r8_abs ( p ) ) ||
        ( p <= q * ( a - x ) ) ||
        ( q * ( b - x ) <= p ) )
      {
        if ( midpoint <= x )
        {
          e = a - x;
        }
        else
        {
          e = b - x;
        }
        d = c * e;
      }
  //
  //  Choose a parabolic interpolation step.
  //
      else
      {
        d = p / q;
        u = x + d;

        if ( ( u - a ) < tol2 )
        {
          d = tol1 * r8_sign ( midpoint - x );
        }

        if ( ( b - u ) < tol2 )
        {
          d = tol1 * r8_sign ( midpoint - x );
        }
      }
    }
  //
  //  F must not be evaluated too close to X.
  //
    if ( tol1 <= r8_abs ( d ) )
    {
      u = x + d;
    }
    if ( r8_abs ( d ) < tol1 )
    {
      u = x + tol1 * r8_sign ( d );
    }
  //
  //  Request value of F(U).
  //
    arg = u;
    status = status + 1;

    return arg;
  }
  //****************************************************************************80

  double r8_abs ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    R8_ABS returns the absolute value of an R8.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    07 May 2006
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the quantity whose absolute value is desired.
  //
  //    Output, double R8_ABS, the absolute value of X.
  //
  {
    double value;

    if ( 0.0 <= x )
    {
      value = x;
    }
    else
    {
      value = - x;
    }
    return value;
  }
  //****************************************************************************80

  double r8_epsilon ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    R8_EPSILON returns the R8 roundoff unit.
  //
  //  Discussion:
  //
  //    The roundoff unit is a number R which is a power of 2 with the
  //    property that, to the precision of the computer's arithmetic,
  //      1 < 1 + R
  //    but
  //      1 = ( 1 + R / 2 )
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    01 September 2012
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Output, double R8_EPSILON, the R8 round-off unit.
  //
  {
    const double value = 2.220446049250313E-016;

    return value;
  }
  //****************************************************************************80

  double r8_max ( double x, double y )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    R8_MAX returns the maximum of two R8's.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    18 August 2004
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, Y, the quantities to compare.
  //
  //    Output, double R8_MAX, the maximum of X and Y.
  //
  {
    double value;

    if ( y < x )
    {
      value = x;
    }
    else
    {
      value = y;
    }
    return value;
  }
  //****************************************************************************80

  double r8_sign ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    R8_SIGN returns the sign of an R8.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    18 October 2004
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the number whose sign is desired.
  //
  //    Output, double R8_SIGN, the sign of X.
  //
  {
    double value;

    if ( x < 0.0 )
    {
      value = -1.0;
    }
    else
    {
      value = 1.0;
    }
    return value;
  }
  //****************************************************************************80

  void timestamp ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TIMESTAMP prints the current YMDHMS date as a time stamp.
  //
  //  Example:
  //
  //    31 May 2001 09:45:54 AM
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    24 September 2003
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    None
  //
  {
    const int TIME_SIZE(40);

    static char time_buffer[TIME_SIZE];
    const struct tm *tm;
    time_t now;

    now = time ( NULL );
    tm = localtime ( &now );

    strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

    cout << time_buffer << "\n";

    return;
  }
  //****************************************************************************80

  double zero ( double a, double b, double t, func_base& f )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    ZERO seeks the root of a function F(X) in an interval [A,B].
  //
  //  Discussion:
  //
  //    The interval [A,B] must be a change of sign interval for F.
  //    That is, F(A) and F(B) must be of opposite signs.  Then
  //    assuming that F is continuous implies the existence of at least
  //    one value C between A and B for which F(C) = 0.
  //
  //    The location of the zero is determined to within an accuracy
  //    of 6 * MACHEPS * r8_abs ( C ) + 2 * T.
  //
  //    Thanks to Thomas Secretin for pointing out a transcription error in the
  //    setting of the value of P, 11 February 2013.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    11 February 2013
  //
  //  Author:
  //
  //    Original FORTRAN77 version by Richard Brent.
  //    C++ version by John Burkardt.
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the change of sign interval.
  //
  //    Input, double T, a positive error tolerance.
  //
  //    Input, func_base& F, the name of a user-supplied c++ functor
  //    whose zero is being sought.  The input and output
  //    of F() are of type double.
  //
  //    Output, double ZERO, the estimated value of a zero of
  //    the function F.
  //
  {
    double c;
    double d;
    double e;
    double fa;
    double fb;
    double fc;
    double m;
    double macheps;
    double p;
    double q;
    double r;
    double s;
    double sa;
    double sb;
    double tol;
  //
  //  Make local copies of A and B.
  //
    sa = a;
    sb = b;
    fa = f ( sa );
    fb = f ( sb );

    c = sa;
    fc = fa;
    e = sb - sa;
    d = e;

    macheps = r8_epsilon ( );

    for ( ; ; )
    {
      if ( r8_abs ( fc ) < r8_abs ( fb ) )
      {
        sa = sb;
        sb = c;
        c = sa;
        fa = fb;
        fb = fc;
        fc = fa;
      }

      tol = 2.0 * macheps * r8_abs ( sb ) + t;
      m = 0.5 * ( c - sb );

      if ( r8_abs ( m ) <= tol || fb == 0.0 )
      {
        break;
      }

      if ( r8_abs ( e ) < tol || r8_abs ( fa ) <= r8_abs ( fb ) )
      {
        e = m;
        d = e;
      }
      else
      {
        s = fb / fa;

        if ( sa == c )
        {
          p = 2.0 * m * s;
          q = 1.0 - s;
        }
        else
        {
          q = fa / fc;
          r = fb / fc;
          p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
          q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
        }

        if ( 0.0 < p )
        {
          q = - q;
        }
        else
        {
          p = - p;
        }

        s = e;
        e = d;

        if ( 2.0 * p < 3.0 * m * q - r8_abs ( tol * q ) &&
          p < r8_abs ( 0.5 * s * q ) )
        {
          d = p / q;
        }
        else
        {
          e = m;
          d = e;
        }
      }
      sa = sb;
      fa = fb;

      if ( tol < r8_abs ( d ) )
      {
        sb = sb + d;
      }
      else if ( 0.0 < m )
      {
        sb = sb + tol;
      }
      else
      {
        sb = sb - tol;
      }

      fb = f ( sb );

      if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
      {
        c = sa;
        fc = fa;
        e = sb - sa;
        d = e;
      }
    }
    return sb;
  }
  //****************************************************************************80

  void zero_rc ( double a, double b, double t, double &arg, int &status,
    double value )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    ZERO_RC seeks the root of a function F(X) using reverse communication.
  //
  //  Discussion:
  //
  //    The interval [A,B] must be a change of sign interval for F.
  //    That is, F(A) and F(B) must be of opposite signs.  Then
  //    assuming that F is continuous implies the existence of at least
  //    one value C between A and B for which F(C) = 0.
  //
  //    The location of the zero is determined to within an accuracy
  //    of 6 * MACHEPS * r8_abs ( C ) + 2 * T.
  //
  //    The routine is a revised version of the Brent zero finder
  //    algorithm, using reverse communication.
  //
  //    Thanks to Thomas Secretin for pointing out a transcription error in the
  //    setting of the value of P, 11 February 2013.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    11 February 2013
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the change of sign interval.
  //
  //    Input, double T, a positive error tolerance.
  //
  //    Output, double &ARG, the currently considered point.  The user
  //    does not need to initialize this value.  On return with STATUS positive,
  //    the user is requested to evaluate the function at ARG, and return
  //    the value in VALUE.  On return with STATUS zero, ARG is the routine's
  //    estimate for the function's zero.
  //
  //    Input/output, int &STATUS, used to communicate between
  //    the user and the routine.  The user only sets STATUS to zero on the first
  //    call, to indicate that this is a startup call.  The routine returns STATUS
  //    positive to request that the function be evaluated at ARG, or returns
  //    STATUS as 0, to indicate that the iteration is complete and that
  //    ARG is the estimated zero
  //
  //    Input, double VALUE, the function value at ARG, as requested
  //    by the routine on the previous call.
  //
  {
    static double c;
    static double d;
    static double e;
    static double fa;
    static double fb;
    static double fc;
    double m;
    static double macheps;
    double p;
    double q;
    double r;
    double s;
    static double sa;
    static double sb;
    double tol;
  //
  //  Input STATUS = 0.
  //  Initialize, request F(A).
  //
    if ( status == 0 )
    {
      macheps = r8_epsilon ( );

      sa = a;
      sb = b;
      e = sb - sa;
      d = e;

      status = 1;
      arg = a;
      return;
    }
  //
  //  Input STATUS = 1.
  //  Receive F(A), request F(B).
  //
    else if ( status == 1 )
    {
      fa = value;
      status = 2;
      arg = sb;
      return;
    }
  //
  //  Input STATUS = 2
  //  Receive F(B).
  //
    else if ( status == 2 )
    {
      fb = value;

      if ( 0.0 < fa * fb )
      {
        status = -1;
        return;
      }
      c = sa;
      fc = fa;
    }
    else
    {
      fb = value;

      if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
      {
        c = sa;
        fc = fa;
        e = sb - sa;
        d = e;
      }
    }
  //
  //  Compute the next point at which a function value is requested.
  //
    if ( r8_abs ( fc ) < r8_abs ( fb ) )
    {
      sa = sb;
      sb = c;
      c = sa;
      fa = fb;
      fb = fc;
      fc = fa;
    }

    tol = 2.0 * macheps * r8_abs ( sb ) + t;
    m = 0.5 * ( c - sb );

    if ( r8_abs ( m ) <= tol || fb == 0.0 )
    {
      status = 0;
      arg = sb;
      return;
    }

    if ( r8_abs ( e ) < tol || r8_abs ( fa ) <= r8_abs ( fb ) )
    {
      e = m;
      d = e;
    }
    else
    {
      s = fb / fa;

      if ( sa == c )
      {
        p = 2.0 * m * s;
        q = 1.0 - s;
      }
      else
      {
        q = fa / fc;
        r = fb / fc;
        p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
        q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
      }

      if ( 0.0 < p )
      {
        q = - q;
      }
      else
      {
        p = - p;
      }
      s = e;
      e = d;

      if ( 2.0 * p < 3.0 * m * q - r8_abs ( tol * q ) &&
           p < r8_abs ( 0.5 * s * q ) )
      {
        d = p / q;
      }
      else
      {
        e = m;
        d = e;
      }
    }

    sa = sb;
    fa = fb;

    if ( tol < r8_abs ( d ) )
    {
      sb = sb + d;
    }
    else if ( 0.0 < m )
    {
      sb = sb + tol;
    }
    else
    {
      sb = sb - tol;
    }

    arg = sb;
    status = status + 1;

    return;
  }

  // ======================================================================
  // === Simple wrapper functions
  // === for convenience and/or compatibility.
  //
  // === The three functions are the same as above,
  // === except that they take a plain function F
  // === instead of a c++ functor.  In all cases, the
  // === input and output of F() are of type double.

  typedef double DoubleOfDouble (double);

  class func_wrapper : public func_base {
    DoubleOfDouble* func;
  public:
    func_wrapper(DoubleOfDouble* f) {
      func = f;
    }
    virtual double operator() (double x){
      return func(x);
    }
  };

  //****************************************************************************80

  double glomin ( double a, double b, double c, double m, double e,
           double t, double f ( double x ), double &x ){
    func_wrapper foo(f);
    return glomin(a, b, c, m, e, t, foo, x);
  }

  //****************************************************************************80

  double local_min ( double a, double b, double t, double f ( double x ),
    double &x ){
    func_wrapper foo(f);
    return local_min(a, b, t, foo, x);
  }

  //****************************************************************************80

  double zero ( double a, double b, double t, double f ( double x ) ){
    func_wrapper foo(f);
    return zero(a, b, t, foo);
  }

  // ======================================================================
  // Generally useful functor to evaluate a monic polynomial.
  // For details, see class definition in brent.hpp

  double monicPoly::operator()(double x){
    double rslt(1);
    for (int ii = coeff.size()-1; ii >= 0; ii--){
      rslt *= x;
      rslt += coeff[ii];
    }
    return rslt;
  }

  // Similarly, evaluate a general polynomial (not necessarily monic):
  double Poly::operator()(double x){
    double rslt(0);
    for (int ii = coeff.size()-1; ii >= 0; ii--){
      rslt *= x;
      rslt += coeff[ii];
    }
    return rslt;
  }

  } // end namespace brent

  Friday, October 9, 2015 12:43 AM
• 

  

  0

  Sign in to vote

  //brent.cpp

  #include <cstdlib>
  //#include <windows.h>
  #include <iostream>

  #include <cmath>

  #include <ctime>

  
  using namespace std;

  
  # include "brent.hpp"

  
  namespace brent{

  
  //****************************************************************************80

  double glomin ( double a, double b, double c, double m, double e, double t,
    func_base& f, double &x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    GLOMIN seeks a global minimum of a function F(X) in an interval [A,B].
  //
  //  Discussion:
  //
  //    This function assumes that F(X) is twice continuously differentiable
  //    over [A,B] and that F''(X) <= M for all X in [A,B].
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    Original FORTRAN77 version by Richard Brent.
  //    C++ version by John Burkardt.
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the interval.
  //    It must be the case that A < B.
  //
  //    Input, double C, an initial guess for the global
  //    minimizer.  If no good guess is known, C = A or B is acceptable.
  //
  //    Input, double M, the bound on the second derivative.
  //
  //    Input, double E, a positive tolerance, a bound for the
  //    absolute error in the evaluation of F(X) for any X in [A,B].
  //
  //    Input, double T, a positive error tolerance.
  //
  //    Input, func_base& F, a user-supplied c++ functor whose
  //    global minimum is being sought.  The input and output
  //    of F() are of type double.
  //
  //    Output, double &X, the estimated value of the abscissa
  //    for which F attains its global minimum value in [A,B].
  //
  //    Output, double GLOMIN, the value F(X).
  //
  {
    double a0;
    double a2;
    double a3;
    double d0;
    double d1;
    double d2;
    double h;
    int k;
    double m2;
    double macheps;
    double p;
    double q;
    double qs;
    double r;
    double s;
    double sc;
    double y;
    double y0;
    double y1;
    double y2;
    double y3;
    double yb;
    double z0;
    double z1;
    double z2;
   
    a0 = b;
    x = a0;
    a2 = a;
    y0 = f ( b );
    yb = y0;
    y2 = f ( a );
    y = y2;

    if ( y0 < y )
    {
      y = y0;
    }
    else
    {
      x = a;
    }

    if ( m <= 0.0 || b <= a )
    {
      return y;
    }

    macheps = r8_epsilon ( );

    m2 = 0.5 * ( 1.0 + 16.0 * macheps ) * m;

    if ( c <= a || b <= c )
    {
      sc = 0.5 * ( a + b );
    }
    else
    {
      sc = c;
    }

    y1 = f ( sc );
    k = 3;
    d0 = a2 - sc;
    h = 9.0 / 11.0;

    if ( y1 < y )
    {
      x = sc;
      y = y1;
    }
  //
  //  Loop.
  //
    for ( ; ; )
    {
      d1 = a2 - a0;
      d2 = sc - a0;
      z2 = b - a2;
      z0 = y2 - y1;
      z1 = y2 - y0;
      r = d1 * d1 * z0 - d0 * d0 * z1;
      p = r;
      qs = 2.0 * ( d0 * z1 - d1 * z0 );
      q = qs;

      if ( k < 1000000 || y2 <= y )
      {
        for ( ; ; )
        {
          if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
            z2 * m2 * r * ( z2 * q - r ) )
          {
            a3 = a2 + r / q;
            y3 = f ( a3 );

            if ( y3 < y )
            {
              x = a3;
              y = y3;
            }
          }
          k = ( ( 1611 * k ) % 1048576 );
          q = 1.0;
          r = ( b - a ) * 0.00001 * ( double ) ( k );

          if ( z2 <= r )
          {
            break;
          }
        }
      }
      else
      {
        k = ( ( 1611 * k ) % 1048576 );
        q = 1.0;
        r = ( b - a ) * 0.00001 * ( double ) ( k );

        while ( r < z2 )
        {
          if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
            z2 * m2 * r * ( z2 * q - r ) )
          {
            a3 = a2 + r / q;
            y3 = f ( a3 );

            if ( y3 < y )
            {
              x = a3;
              y = y3;
            }
          }
          k = ( ( 1611 * k ) % 1048576 );
          q = 1.0;
          r = ( b - a ) * 0.00001 * ( double ) ( k );
        }
      }

      r = m2 * d0 * d1 * d2;
      s = sqrt ( ( ( y2 - y ) + t ) / m2 );
      h = 0.5 * ( 1.0 + h );
      p = h * ( p + 2.0 * r * s );
      q = q + 0.5 * qs;
      r = - 0.5 * ( d0 + ( z0 + 2.01 * e ) / ( d0 * m2 ) );

      if ( r < s || d0 < 0.0 )
      {
        r = a2 + s;
      }
      else
      {
        r = a2 + r;
      }

      if ( 0.0 < p * q )
      {
        a3 = a2 + p / q;
      }
      else
      {
        a3 = r;
      }

      for ( ; ; )
      {
        a3 = r8_max ( a3, r );

        if ( b <= a3 )
        {
          a3 = b;
          y3 = yb;
        }
        else
        {
          y3 = f ( a3 );
        }

        if ( y3 < y )
        {
          x = a3;
          y = y3;
        }

        d0 = a3 - a2;

        if ( a3 <= r )
        {
          break;
        }

        p = 2.0 * ( y2 - y3 ) / ( m * d0 );

        if ( ( 1.0 + 9.0 * macheps ) * d0 <= r8_abs ( p ) )
        {
          break;
        }

        if ( 0.5 * m2 * ( d0 * d0 + p * p ) <= ( y2 - y ) + ( y3 - y ) + 2.0 * t )
        {
          break;
        }
        a3 = 0.5 * ( a2 + a3 );
        h = 0.9 * h;
      }

      if ( b <= a3 )
      {
        break;
      }

      a0 = sc;
      sc = a2;
      a2 = a3;
      y0 = y1;
      y1 = y2;
      y2 = y3;
    }

    return y;
  }
  //****************************************************************************80

  double local_min ( double a, double b, double t, func_base& f,
    double &x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    LOCAL_MIN seeks a local minimum of a function F(X) in an interval [A,B].
  //
  //  Discussion:
  //
  //    The method used is a combination of golden section search and
  //    successive parabolic interpolation.  Convergence is never much slower
  //    than that for a Fibonacci search.  If F has a continuous second
  //    derivative which is positive at the minimum (which is not at A or
  //    B), then convergence is superlinear, and usually of the order of
  //    about 1.324....
  //
  //    The values EPS and T define a tolerance TOL = EPS * abs ( X ) + T.
  //    F is never evaluated at two points closer than TOL.
  //
  //    If F is a unimodal function and the computed values of F are always
  //    unimodal when separated by at least SQEPS * abs ( X ) + (T/3), then
  //    LOCAL_MIN approximates the abscissa of the global minimum of F on the
  //    interval [A,B] with an error less than 3*SQEPS*abs(LOCAL_MIN)+T.
  //
  //    If F is not unimodal, then LOCAL_MIN may approximate a local, but
  //    perhaps non-global, minimum to the same accuracy.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    Original FORTRAN77 version by Richard Brent.
  //    C++ version by John Burkardt.
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the interval.
  //
  //    Input, double T, a positive absolute error tolerance.
  //
  //    Input, func_base& F, a user-supplied c++ functor whose
  //    local minimum is being sought.  The input and output
  //    of F() are of type double.
  //
  //    Output, double &X, the estimated value of an abscissa
  //    for which F attains a local minimum value in [A,B].
  //
  //    Output, double LOCAL_MIN, the value F(X).
  //
  {
    double c;
    double d;
    double e;
    double eps;
    double fu;
    double fv;
    double fw;
    double fx;
    double m;
    double p;
    double q;
    double r;
    double sa;
    double sb;
    double t2;
    double tol;
    double u;
    double v;
    double w;
  //
  //  C is the square of the inverse of the golden ratio.
  //
    c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );

    eps = sqrt ( r8_epsilon ( ) );

    sa = a;
    sb = b;
    x = sa + c * ( b - a );
    w = x;
    v = w;
    e = 0.0;
    fx = f ( x );
    fw = fx;
    fv = fw;

    for ( ; ; )
    {
      m = 0.5 * ( sa + sb ) ;
      tol = eps * r8_abs ( x ) + t;
      t2 = 2.0 * tol;
  //
  //  Check the stopping criterion.
  //
      if ( r8_abs ( x - m ) <= t2 - 0.5 * ( sb - sa ) )
      {
        break;
      }
  //
  //  Fit a parabola.
  //
      r = 0.0;
      q = r;
      p = q;

      if ( tol < r8_abs ( e ) )
      {
        r = ( x - w ) * ( fx - fv );
        q = ( x - v ) * ( fx - fw );
        p = ( x - v ) * q - ( x - w ) * r;
        q = 2.0 * ( q - r );
        if ( 0.0 < q )
        {
          p = - p;
        }
        q = r8_abs ( q );
        r = e;
        e = d;
      }

      if ( r8_abs ( p ) < r8_abs ( 0.5 * q * r ) &&
           q * ( sa - x ) < p &&
           p < q * ( sb - x ) )
      {
  //
  //  Take the parabolic interpolation step.
  //
        d = p / q;
        u = x + d;
  //
  //  F must not be evaluated too close to A or B.
  //
        if ( ( u - sa ) < t2 || ( sb - u ) < t2 )
        {
          if ( x < m )
          {
            d = tol;
          }
          else
          {
            d = - tol;
          }
        }
      }
  //
  //  A golden-section step.
  //
      else
      {
        if ( x < m )
        {
          e = sb - x;
        }
        else
        {
          e = sa - x;
        }
        d = c * e;
      }
  //
  //  F must not be evaluated too close to X.
  //
      if ( tol <= r8_abs ( d ) )
      {
        u = x + d;
      }
      else if ( 0.0 < d )
      {
        u = x + tol;
      }
      else
      {
        u = x - tol;
      }

      fu = f ( u );
  //
  //  Update A, B, V, W, and X.
  //
      if ( fu <= fx )
      {
        if ( u < x )
        {
          sb = x;
        }
        else
        {
          sa = x;
        }
        v = w;
        fv = fw;
        w = x;
        fw = fx;
        x = u;
        fx = fu;
      }
      else
      {
        if ( u < x )
        {
          sa = u;
        }
        else
        {
          sb = u;
        }

        if ( fu <= fw || w == x )
        {
          v = w;
          fv = fw;
          w = u;
          fw = fu;
        }
        else if ( fu <= fv || v == x || v == w )
        {
          v = u;
          fv = fu;
        }
      }
    }
    return fx;
  }
  //****************************************************************************80

  double local_min_rc ( double &a, double &b, int &status, double value )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    LOCAL_MIN_RC seeks a minimizer of a scalar function of a scalar variable.
  //
  //  Discussion:
  //
  //    This routine seeks an approximation to the point where a function
  //    F attains a minimum on the interval (A,B).
  //
  //    The method used is a combination of golden section search and
  //    successive parabolic interpolation.  Convergence is never much
  //    slower than that for a Fibonacci search.  If F has a continuous
  //    second derivative which is positive at the minimum (which is not
  //    at A or B), then convergence is superlinear, and usually of the
  //    order of about 1.324...
  //
  //    The routine is a revised version of the Brent local minimization
  //    algorithm, using reverse communication.
  //
  //    It is worth stating explicitly that this routine will NOT be
  //    able to detect a minimizer that occurs at either initial endpoint
  //    A or B.  If this is a concern to the user, then the user must
  //    either ensure that the initial interval is larger, or to check
  //    the function value at the returned minimizer against the values
  //    at either endpoint.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    17 July 2011
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //    David Kahaner, Cleve Moler, Steven Nash,
  //    Numerical Methods and Software,
  //    Prentice Hall, 1989,
  //    ISBN: 0-13-627258-4,
  //    LC: TA345.K34.
  //
  //  Parameters
  //
  //    Input/output, double &A, &B.  On input, the left and right
  //    endpoints of the initial interval.  On output, the lower and upper
  //    bounds for an interval containing the minimizer.  It is required
  //    that A < B.
  //
  //    Input/output, int &STATUS, used to communicate between
  //    the user and the routine.  The user only sets STATUS to zero on the first
  //    call, to indicate that this is a startup call.  The routine returns STATUS
  //    positive to request that the function be evaluated at ARG, or returns
  //    STATUS as 0, to indicate that the iteration is complete and that
  //    ARG is the estimated minimizer.
  //
  //    Input, double VALUE, the function value at ARG, as requested
  //    by the routine on the previous call.
  //
  //    Output, double LOCAL_MIN_RC, the currently considered point.
  //    On return with STATUS positive, the user is requested to evaluate the
  //    function at this point, and return the value in VALUE.  On return with
  //    STATUS zero, this is the routine's estimate for the function minimizer.
  //
  //  Local parameters:
  //
  //    C is the squared inverse of the golden ratio.
  //
  //    EPS is the square root of the relative machine precision.
  //
  {
    static double arg;
    static double c;
    static double d;
    static double e;
    static double eps;
    static double fu;
    static double fv;
    static double fw;
    static double fx;
    static double midpoint;
    static double p;
    static double q;
    static double r;
    static double tol;
    static double tol1;
    static double tol2;
    static double u;
    static double v;
    static double w;
    static double x;
  //
  //  STATUS (INPUT) = 0, startup.
  //
    if ( status == 0 )
    {
      if ( b <= a )
      {
        cout << "\n";
        cout << "LOCAL_MIN_RC - Fatal error!\n";
        cout << "  A < B is required, but\n";
        cout << "  A = " << a << "\n";
        cout << "  B = " << b << "\n";
        status = -1;
        exit ( 1 );
      }
      c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );

      eps = sqrt ( r8_epsilon ( ) );
      tol = r8_epsilon ( );

      v = a + c * ( b - a );
      w = v;
      x = v;
      e = 0.0;

      status = 1;
      arg = x;

      return arg;
    }
  //
  //  STATUS (INPUT) = 1, return with initial function value of FX.
  //
    else if ( status == 1 )
    {
      fx = value;
      fv = fx;
      fw = fx;
    }
  //
  //  STATUS (INPUT) = 2 or more, update the data.
  //
    else if ( 2 <= status )
    {
      fu = value;

      if ( fu <= fx )
      {
        if ( x <= u )
        {
          a = x;
        }
        else
        {
          b = x;
        }
        v = w;
        fv = fw;
        w = x;
        fw = fx;
        x = u;
        fx = fu;
      }
      else
      {
        if ( u < x )
        {
          a = u;
        }
        else
        {
          b = u;
        }

        if ( fu <= fw || w == x )
        {
          v = w;
          fv = fw;
          w = u;
          fw = fu;
        }
        else if ( fu <= fv || v == x || v == w )
        {
          v = u;
          fv = fu;
        }
      }
    }
  //
  //  Take the next step.
  //
    midpoint = 0.5 * ( a + b );
    tol1 = eps * r8_abs ( x ) + tol / 3.0;
    tol2 = 2.0 * tol1;
  //
  //  If the stopping criterion is satisfied, we can exit.
  //
    if ( r8_abs ( x - midpoint ) <= ( tol2 - 0.5 * ( b - a ) ) )
    {
      status = 0;
      return arg;
    }
  //
  //  Is golden-section necessary?
  //
    if ( r8_abs ( e ) <= tol1 )
    {
      if ( midpoint <= x )
      {
        e = a - x;
      }
      else
      {
        e = b - x;
      }
      d = c * e;
    }
  //
  //  Consider fitting a parabola.
  //
    else
    {
      r = ( x - w ) * ( fx - fv );
      q = ( x - v ) * ( fx - fw );
      p = ( x - v ) * q - ( x - w ) * r;
      q = 2.0 * ( q - r );
      if ( 0.0 < q )
      {
        p = - p;
      }
      q = r8_abs ( q );
      r = e;
      e = d;
  //
  //  Choose a golden-section step if the parabola is not advised.
  //
      if (
        ( r8_abs ( 0.5 * q * r ) <= r8_abs ( p ) ) ||
        ( p <= q * ( a - x ) ) ||
        ( q * ( b - x ) <= p ) )
      {
        if ( midpoint <= x )
        {
          e = a - x;
        }
        else
        {
          e = b - x;
        }
        d = c * e;
      }
  //
  //  Choose a parabolic interpolation step.
  //
      else
      {
        d = p / q;
        u = x + d;

        if ( ( u - a ) < tol2 )
        {
          d = tol1 * r8_sign ( midpoint - x );
        }

        if ( ( b - u ) < tol2 )
        {
          d = tol1 * r8_sign ( midpoint - x );
        }
      }
    }
  //
  //  F must not be evaluated too close to X.
  //
    if ( tol1 <= r8_abs ( d ) )
    {
      u = x + d;
    }
    if ( r8_abs ( d ) < tol1 )
    {
      u = x + tol1 * r8_sign ( d );
    }
  //
  //  Request value of F(U).
  //
    arg = u;
    status = status + 1;

    return arg;
  }
  //****************************************************************************80

  double r8_abs ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    R8_ABS returns the absolute value of an R8.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    07 May 2006
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the quantity whose absolute value is desired.
  //
  //    Output, double R8_ABS, the absolute value of X.
  //
  {
    double value;

    if ( 0.0 <= x )
    {
      value = x;
    }
    else
    {
      value = - x;
    }
    return value;
  }
  //****************************************************************************80

  double r8_epsilon ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    R8_EPSILON returns the R8 roundoff unit.
  //
  //  Discussion:
  //
  //    The roundoff unit is a number R which is a power of 2 with the
  //    property that, to the precision of the computer's arithmetic,
  //      1 < 1 + R
  //    but
  //      1 = ( 1 + R / 2 )
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    01 September 2012
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Output, double R8_EPSILON, the R8 round-off unit.
  //
  {
    const double value = 2.220446049250313E-016;

    return value;
  }
  //****************************************************************************80

  double r8_max ( double x, double y )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    R8_MAX returns the maximum of two R8's.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    18 August 2004
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, Y, the quantities to compare.
  //
  //    Output, double R8_MAX, the maximum of X and Y.
  //
  {
    double value;

    if ( y < x )
    {
      value = x;
    }
    else
    {
      value = y;
    }
    return value;
  }
  //****************************************************************************80

  double r8_sign ( double x )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    R8_SIGN returns the sign of an R8.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    18 October 2004
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    Input, double X, the number whose sign is desired.
  //
  //    Output, double R8_SIGN, the sign of X.
  //
  {
    double value;

    if ( x < 0.0 )
    {
      value = -1.0;
    }
    else
    {
      value = 1.0;
    }
    return value;
  }
  //****************************************************************************80

  void timestamp ( )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    TIMESTAMP prints the current YMDHMS date as a time stamp.
  //
  //  Example:
  //
  //    31 May 2001 09:45:54 AM
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    24 September 2003
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Parameters:
  //
  //    None
  //
  {
    const int TIME_SIZE(40);

    static char time_buffer[TIME_SIZE];
    const struct tm *tm;
    time_t now;

    now = time ( NULL );
    tm = localtime ( &now );

    strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

    cout << time_buffer << "\n";

    return;
  }
  //****************************************************************************80

  double zero ( double a, double b, double t, func_base& f )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    ZERO seeks the root of a function F(X) in an interval [A,B].
  //
  //  Discussion:
  //
  //    The interval [A,B] must be a change of sign interval for F.
  //    That is, F(A) and F(B) must be of opposite signs.  Then
  //    assuming that F is continuous implies the existence of at least
  //    one value C between A and B for which F(C) = 0.
  //
  //    The location of the zero is determined to within an accuracy
  //    of 6 * MACHEPS * r8_abs ( C ) + 2 * T.
  //
  //    Thanks to Thomas Secretin for pointing out a transcription error in the
  //    setting of the value of P, 11 February 2013.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    11 February 2013
  //
  //  Author:
  //
  //    Original FORTRAN77 version by Richard Brent.
  //    C++ version by John Burkardt.
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the change of sign interval.
  //
  //    Input, double T, a positive error tolerance.
  //
  //    Input, func_base& F, the name of a user-supplied c++ functor
  //    whose zero is being sought.  The input and output
  //    of F() are of type double.
  //
  //    Output, double ZERO, the estimated value of a zero of
  //    the function F.
  //
  {
    double c;
    double d;
    double e;
    double fa;
    double fb;
    double fc;
    double m;
    double macheps;
    double p;
    double q;
    double r;
    double s;
    double sa;
    double sb;
    double tol;
  //
  //  Make local copies of A and B.
  //
    sa = a;
    sb = b;
    fa = f ( sa );
    fb = f ( sb );

    c = sa;
    fc = fa;
    e = sb - sa;
    d = e;

    macheps = r8_epsilon ( );

    for ( ; ; )
    {
      if ( r8_abs ( fc ) < r8_abs ( fb ) )
      {
        sa = sb;
        sb = c;
        c = sa;
        fa = fb;
        fb = fc;
        fc = fa;
      }

      tol = 2.0 * macheps * r8_abs ( sb ) + t;
      m = 0.5 * ( c - sb );

      if ( r8_abs ( m ) <= tol || fb == 0.0 )
      {
        break;
      }

      if ( r8_abs ( e ) < tol || r8_abs ( fa ) <= r8_abs ( fb ) )
      {
        e = m;
        d = e;
      }
      else
      {
        s = fb / fa;

        if ( sa == c )
        {
          p = 2.0 * m * s;
          q = 1.0 - s;
        }
        else
        {
          q = fa / fc;
          r = fb / fc;
          p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
          q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
        }

        if ( 0.0 < p )
        {
          q = - q;
        }
        else
        {
          p = - p;
        }

        s = e;
        e = d;

        if ( 2.0 * p < 3.0 * m * q - r8_abs ( tol * q ) &&
          p < r8_abs ( 0.5 * s * q ) )
        {
          d = p / q;
        }
        else
        {
          e = m;
          d = e;
        }
      }
      sa = sb;
      fa = fb;

      if ( tol < r8_abs ( d ) )
      {
        sb = sb + d;
      }
      else if ( 0.0 < m )
      {
        sb = sb + tol;
      }
      else
      {
        sb = sb - tol;
      }

      fb = f ( sb );

      if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
      {
        c = sa;
        fc = fa;
        e = sb - sa;
        d = e;
      }
    }
    return sb;
  }
  //****************************************************************************80

  void zero_rc ( double a, double b, double t, double &arg, int &status,
    double value )

  //****************************************************************************80
  //
  //  Purpose:
  //
  //    ZERO_RC seeks the root of a function F(X) using reverse communication.
  //
  //  Discussion:
  //
  //    The interval [A,B] must be a change of sign interval for F.
  //    That is, F(A) and F(B) must be of opposite signs.  Then
  //    assuming that F is continuous implies the existence of at least
  //    one value C between A and B for which F(C) = 0.
  //
  //    The location of the zero is determined to within an accuracy
  //    of 6 * MACHEPS * r8_abs ( C ) + 2 * T.
  //
  //    The routine is a revised version of the Brent zero finder
  //    algorithm, using reverse communication.
  //
  //    Thanks to Thomas Secretin for pointing out a transcription error in the
  //    setting of the value of P, 11 February 2013.
  //
  //  Licensing:
  //
  //    This code is distributed under the GNU LGPL license.
  //
  //  Modified:
  //
  //    11 February 2013
  //
  //  Author:
  //
  //    John Burkardt
  //
  //  Reference:
  //
  //    Richard Brent,
  //    Algorithms for Minimization Without Derivatives,
  //    Dover, 2002,
  //    ISBN: 0-486-41998-3,
  //    LC: QA402.5.B74.
  //
  //  Parameters:
  //
  //    Input, double A, B, the endpoints of the change of sign interval.
  //
  //    Input, double T, a positive error tolerance.
  //
  //    Output, double &ARG, the currently considered point.  The user
  //    does not need to initialize this value.  On return with STATUS positive,
  //    the user is requested to evaluate the function at ARG, and return
  //    the value in VALUE.  On return with STATUS zero, ARG is the routine's
  //    estimate for the function's zero.
  //
  //    Input/output, int &STATUS, used to communicate between
  //    the user and the routine.  The user only sets STATUS to zero on the first
  //    call, to indicate that this is a startup call.  The routine returns STATUS
  //    positive to request that the function be evaluated at ARG, or returns
  //    STATUS as 0, to indicate that the iteration is complete and that
  //    ARG is the estimated zero
  //
  //    Input, double VALUE, the function value at ARG, as requested
  //    by the routine on the previous call.
  //
  {
    static double c;
    static double d;
    static double e;
    static double fa;
    static double fb;
    static double fc;
    double m;
    static double macheps;
    double p;
    double q;
    double r;
    double s;
    static double sa;
    static double sb;
    double tol;
  //
  //  Input STATUS = 0.
  //  Initialize, request F(A).
  //
    if ( status == 0 )
    {
      macheps = r8_epsilon ( );

      sa = a;
      sb = b;
      e = sb - sa;
      d = e;

      status = 1;
      arg = a;
      return;
    }
  //
  //  Input STATUS = 1.
  //  Receive F(A), request F(B).
  //
    else if ( status == 1 )
    {
      fa = value;
      status = 2;
      arg = sb;
      return;
    }
  //
  //  Input STATUS = 2
  //  Receive F(B).
  //
    else if ( status == 2 )
    {
      fb = value;

      if ( 0.0 < fa * fb )
      {
        status = -1;
        return;
      }
      c = sa;
      fc = fa;
    }
    else
    {
      fb = value;

      if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
      {
        c = sa;
        fc = fa;
        e = sb - sa;
        d = e;
      }
    }
  //
  //  Compute the next point at which a function value is requested.
  //
    if ( r8_abs ( fc ) < r8_abs ( fb ) )
    {
      sa = sb;
      sb = c;
      c = sa;
      fa = fb;
      fb = fc;
      fc = fa;
    }

    tol = 2.0 * macheps * r8_abs ( sb ) + t;
    m = 0.5 * ( c - sb );

    if ( r8_abs ( m ) <= tol || fb == 0.0 )
    {
      status = 0;
      arg = sb;
      return;
    }

    if ( r8_abs ( e ) < tol || r8_abs ( fa ) <= r8_abs ( fb ) )
    {
      e = m;
      d = e;
    }
    else
    {
      s = fb / fa;

      if ( sa == c )
      {
        p = 2.0 * m * s;
        q = 1.0 - s;
      }
      else
      {
        q = fa / fc;
        r = fb / fc;
        p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
        q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
      }

      if ( 0.0 < p )
      {
        q = - q;
      }
      else
      {
        p = - p;
      }
      s = e;
      e = d;

      if ( 2.0 * p < 3.0 * m * q - r8_abs ( tol * q ) &&
           p < r8_abs ( 0.5 * s * q ) )
      {
        d = p / q;
      }
      else
      {
        e = m;
        d = e;
      }
    }

    sa = sb;
    fa = fb;

    if ( tol < r8_abs ( d ) )
    {
      sb = sb + d;
    }
    else if ( 0.0 < m )
    {
      sb = sb + tol;
    }
    else
    {
      sb = sb - tol;
    }

    arg = sb;
    status = status + 1;

    return;
  }

  // ======================================================================
  // === Simple wrapper functions
  // === for convenience and/or compatibility.
  //
  // === The three functions are the same as above,
  // === except that they take a plain function F
  // === instead of a c++ functor.  In all cases, the
  // === input and output of F() are of type double.

  typedef double DoubleOfDouble (double);

  class func_wrapper : public func_base {
    DoubleOfDouble* func;
  public:
    func_wrapper(DoubleOfDouble* f) {
      func = f;
    }
    virtual double operator() (double x){
      return func(x);
    }
  };

  //****************************************************************************80

  double glomin ( double a, double b, double c, double m, double e,
           double t, double f ( double x ), double &x ){
    func_wrapper foo(f);
    return glomin(a, b, c, m, e, t, foo, x);
  }

  //****************************************************************************80

  double local_min ( double a, double b, double t, double f ( double x ),
    double &x ){
    func_wrapper foo(f);
    return local_min(a, b, t, foo, x);
  }

  //****************************************************************************80

  double zero ( double a, double b, double t, double f ( double x ) ){
    func_wrapper foo(f);
    return zero(a, b, t, foo);
  }

  // ======================================================================
  // Generally useful functor to evaluate a monic polynomial.
  // For details, see class definition in brent.hpp

  double monicPoly::operator()(double x){
    double rslt(1);
    for (int ii = coeff.size()-1; ii >= 0; ii--){
      rslt *= x;
      rslt += coeff[ii];
    }
    return rslt;
  }

  // Similarly, evaluate a general polynomial (not necessarily monic):
  double Poly::operator()(double x){
    double rslt(0);
    for (int ii = coeff.size()-1; ii >= 0; ii--){
      rslt *= x;
      rslt += coeff[ii];
    }
    return rslt;
  }

  } // end namespace brent

  Friday, October 9, 2015 12:45 AM
• 

  

  0

  Sign in to vote

  //montecarlo.cpp

  #include <stdlib.h> // or <cstdlib>
  #include <stdio.h> // or <cstdio>
  #include <time.h>  // or <ctime>

  #include <math.h> // or <cmath>

  /*

  #include <cstdlib>
  #include <cstdio>
  #include <ctime>

  */
   
  double func(double x)
  {
         // return 1-x*x;
         return x*x;   // y=x*x
  }
   
  int main()
  {
          int point_a = 1; //-1
          int point_b = 5; //1
          int number_of_random=0;
          double s=0;
          printf ("how much random points /100000000/? ");
          scanf ("%ld",&number_of_random);
          srand ((unsigned) time(NULL));
          int i;
          for (  i=0; i<number_of_random; i++)
          {
               s+=func (point_a+((double)rand()/RAND_MAX*(point_b-point_a)));
          }
          s=s/(double)number_of_random*(point_b-point_a);
          printf ("\nintegral %f \n",s);
          return 0;
  }

  Friday, October 9, 2015 2:28 AM