# Obtaining integral of multiparametrical functions in C++ • ### Question

• testintegral.cpp

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

#include "integral1.h"

double function1(double x, double y, double z)
{

return x*x;
}

int main(void)
{
double a=1;
double b=5;
double eps=0.00000001;
double res;
res=IntRomberg(a,b,eps,0,1,function1);
printf("\n a=%lf   \nb= %lf \n eps= %lf  \n res=%lf  ", a,b,eps, res);

getch();
return 0;

}

integral1.h

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

double IntRomberg(double a,double b,double eps,double razm,double lyambda, double(*f) (double , double ,double  ) )
{

double h,s,s0,s1,sn    ;
int i,n;
s=1;
sn=101;
n=4;
s0=(f(a,razm,lyambda)+f(b,razm,lyambda))/2;
s1=f((a+b)/2,razm,lyambda);
while (fabs(s-sn)>eps)
{
sn=s;
h=(b-a)/n;
for (i=0; i<n/2; i++)
s1+=f(  (a+(2*i+1)*h),razm,lyambda );
s=h*(s0+s1);
n*=2;
}
return s;
}

Sunday, September 27, 2015 7:56 PM

### All replies

• 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
// 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
• Simply rootfinder

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

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

}

double rootsvect;

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 Sunday, September 27, 2015 9:07 PM
Sunday, September 27, 2015 9:04 PM
• 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={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;
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 Saturday, October 3, 2015 6:40 PM
Saturday, October 3, 2015 6:39 PM
• http://www.mcs.csueastbay.edu/~malek/TeX/root.pdf

Saturday, October 3, 2015 6:44 PM
• ## <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 Saturday, October 3, 2015 6:55 PM
Saturday, October 3, 2015 6:52 PM
• #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 ;
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]=eps;
//Usermatrix[i]=x;
//Usermatrix[i]=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 Friday, October 9, 2015 12:26 AM
Friday, October 9, 2015 12:25 AM
• #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, 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
• #include<stdio.h>

#include<conio.h>

float f(float x)

{

return(1/(1+x));

}

int main()

{

int i,n;

float x0,xn,h,y,so,se,ans,x;

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+y[n]+4*so+2*se);

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

getch();

}

Friday, October 9, 2015 12:27 AM
• 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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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
• //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
• //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:
//
//
//  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
//
//    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:
//
//
//  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
//
//    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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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.
//
else if ( status == 1 )
{
fa = value;
status = 2;
arg = sb;
return;
}
//
//  Input STATUS = 2
//
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
• //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:
//
//
//  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
//
//    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:
//
//
//  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
//
//    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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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:
//
//
//  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.
//
else if ( status == 1 )
{
fa = value;
status = 2;
arg = sb;
return;
}
//
//  Input STATUS = 2
//
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
• //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