15 double glomin (
double a,
double b,
double c,
double m,
double e,
double t,
16 func_base& f,
double &x )
119 if ( m <= 0.0 || b <= a )
126 m2 = 0.5 * ( 1.0 + 16.0 * macheps ) * m;
128 if ( c <= a || b <= c )
130 sc = 0.5 * ( a + b );
157 r = d1 * d1 * z0 - d0 * d0 * z1;
159 qs = 2.0 * ( d0 * z1 - d1 * z0 );
162 if ( k < 1000000 || y2 <= y )
166 if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
167 z2 * m2 * r * ( z2 * q - r ) )
178 k = ( ( 1611 * k ) % 1048576 );
180 r = ( b - a ) * 0.00001 * (
double ) ( k );
190 k = ( ( 1611 * k ) % 1048576 );
192 r = ( b - a ) * 0.00001 * (
double ) ( k );
196 if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
197 z2 * m2 * r * ( z2 * q - r ) )
208 k = ( ( 1611 * k ) % 1048576 );
210 r = ( b - a ) * 0.00001 * (
double ) ( k );
214 r = m2 * d0 * d1 * d2;
215 s = sqrt ( ( ( y2 - y ) + t ) / m2 );
216 h = 0.5 * ( 1.0 + h );
217 p = h * ( p + 2.0 * r * s );
219 r = - 0.5 * ( d0 + ( z0 + 2.01 * e ) / ( d0 * m2 ) );
221 if ( r < s || d0 < 0.0 )
266 p = 2.0 * ( y2 - y3 ) / ( m * d0 );
268 if ( ( 1.0 + 9.0 * macheps ) * d0 <=
r8_abs ( p ) )
273 if ( 0.5 * m2 * ( d0 * d0 + p * p ) <= ( y2 - y ) + ( y3 - y ) + 2.0 * t )
277 a3 = 0.5 * ( a2 + a3 );
298 double local_min (
double a,
double b,
double t, func_base& f,
387 c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );
393 x = sa + c * ( b - a );
403 m = 0.5 * ( sa + sb ) ;
404 tol = eps *
r8_abs ( x ) + t;
409 if (
r8_abs ( x - m ) <= t2 - 0.5 * ( sb - sa ) )
422 r = ( x - w ) * ( fx - fv );
423 q = ( x - v ) * ( fx - fw );
424 p = ( x - v ) * q - ( x - w ) * r;
436 q * ( sa - x ) < p &&
447 if ( ( u - sa ) < t2 || ( sb - u ) < t2 )
477 if ( tol <=
r8_abs ( d ) )
522 if ( fu <= fw || w == x )
529 else if ( fu <= fv || v == x || v== w )
540 double local_min_rc (
double &a,
double &b,
int &status,
double value )
634 static double midpoint;
653 cout <<
"LOCAL_MIN_RC - Fatal error!\n";
654 cout <<
" A < B is required, but\n";
655 cout <<
" A = " << a <<
"\n";
656 cout <<
" B = " << b <<
"\n";
660 c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );
665 v = a + c * ( b - a );
678 else if ( status == 1 )
687 else if ( 2 <= status )
719 if ( fu <= fw || w == x )
726 else if ( fu <= fv || v == x || v == w )
736 midpoint = 0.5 * ( a + b );
737 tol1 = eps *
r8_abs ( x ) + tol / 3.0;
742 if (
r8_abs ( x - midpoint ) <= ( tol2 - 0.5 * ( b - a ) ) )
750 if (
r8_abs ( e ) <= tol1 )
767 r = ( x - w ) * ( fx - fv );
768 q = ( x - v ) * ( fx - fw );
769 p = ( x - v ) * q - ( x - w ) * r;
783 ( p <= q * ( a - x ) ) ||
784 ( q * ( b - x ) <= p ) )
804 if ( ( u - a ) < tol2 )
806 d = tol1 *
r8_sign ( midpoint - x );
809 if ( ( b - u ) < tol2 )
811 d = tol1 *
r8_sign ( midpoint - x );
818 if ( tol1 <=
r8_abs ( d ) )
822 if (
r8_abs ( d ) < tol1 )
914 while ( 1.0 < (
double ) ( 1.0 + r ) )
1034 # define TIME_SIZE 40
1037 const struct tm *tm;
1040 now = time ( NULL );
1041 tm = localtime ( &now );
1043 strftime ( time_buffer,
TIME_SIZE,
"%d %B %Y %I:%M:%S %p", tm );
1045 cout << time_buffer <<
"\n";
1052 double zero (
double a,
double b,
double t, func_base& f )
1148 tol = 2.0 * macheps *
r8_abs ( sb ) + t;
1149 m = 0.5 * ( c - sb );
1151 if (
r8_abs ( m ) <= tol || fb == 0.0 )
1174 p = s * ( 2.0 * m * a * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
1175 q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
1190 if ( 2.0 * p < 3.0 * m * q -
r8_abs ( tol * q ) &&
1191 p <
r8_abs ( 0.5 * s * q ) )
1204 if ( tol <
r8_abs ( d ) )
1219 if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
1231 void zero_rc (
double a,
double b,
double t,
double &arg,
int &status,
1303 static double macheps;
1332 else if ( status == 1 )
1343 else if ( status == 2 )
1347 if ( 0.0 < fa * fb )
1359 if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
1380 tol = 2.0 * macheps *
r8_abs ( sb ) + t;
1381 m = 0.5 * ( c - sb );
1383 if (
r8_abs ( m ) <= tol || fb == 0.0 )
1408 p = s * ( 2.0 * m * a * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
1409 q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
1423 if ( 2.0 * p < 3.0 * m * q -
r8_abs ( tol * q ) &&
1424 p <
r8_abs ( 0.5 * s * q ) )
1438 if ( tol <
r8_abs ( d ) )
1452 status = status + 1;
1474 virtual double operator() (
double x){
1481 double glomin (
double a,
double b,
double c,
double m,
double e,
1482 double t,
double f (
double x ),
double &x ){
1484 return glomin(a, b, c, m, e, t, foo, x);
1489 double local_min (
double a,
double b,
double t,
double f (
double x ),
1497 double zero (
double a,
double b,
double t,
double f (
double x ) ){
1499 return zero(a, b, t, foo);
1506 double monicPoly::operator()(
double x){
1508 for (
int ii = coeff.size()-1; ii >= 0; ii--){
1516 double Poly::operator()(
double x){
1518 for (
int ii = coeff.size()-1; ii >= 0; ii--){