00001 # include <cstdlib>
00002 # include <iostream>
00003 # include <cmath>
00004 # include <ctime>
00005
00006 using namespace std;
00007
00008 # include "brent.h"
00009
00010 namespace shogun
00011 {
00012
00013
00014
00015 double glomin ( double a, double b, double c, double m, double e, double t,
00016 func_base& f, double &x )
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075 {
00076 double a0;
00077 double a2;
00078 double a3;
00079 double d0;
00080 double d1;
00081 double d2;
00082 double h;
00083 int k;
00084 double m2;
00085 double macheps;
00086 double p;
00087 double q;
00088 double qs;
00089 double r;
00090 double s;
00091 double sc;
00092 double y;
00093 double y0;
00094 double y1;
00095 double y2;
00096 double y3;
00097 double yb;
00098 double z0;
00099 double z1;
00100 double z2;
00101
00102 a0 = b;
00103 x = a0;
00104 a2 = a;
00105 y0 = f ( b );
00106 yb = y0;
00107 y2 = f ( a );
00108 y = y2;
00109
00110 if ( y0 < y )
00111 {
00112 y = y0;
00113 }
00114 else
00115 {
00116 x = a;
00117 }
00118
00119 if ( m <= 0.0 || b <= a )
00120 {
00121 return y;
00122 }
00123
00124 macheps = r8_epsilon ( );
00125
00126 m2 = 0.5 * ( 1.0 + 16.0 * macheps ) * m;
00127
00128 if ( c <= a || b <= c )
00129 {
00130 sc = 0.5 * ( a + b );
00131 }
00132 else
00133 {
00134 sc = c;
00135 }
00136
00137 y1 = f ( sc );
00138 k = 3;
00139 d0 = a2 - sc;
00140 h = 9.0 / 11.0;
00141
00142 if ( y1 < y )
00143 {
00144 x = sc;
00145 y = y1;
00146 }
00147
00148
00149
00150 for ( ; ; )
00151 {
00152 d1 = a2 - a0;
00153 d2 = sc - a0;
00154 z2 = b - a2;
00155 z0 = y2 - y1;
00156 z1 = y2 - y0;
00157 r = d1 * d1 * z0 - d0 * d0 * z1;
00158 p = r;
00159 qs = 2.0 * ( d0 * z1 - d1 * z0 );
00160 q = qs;
00161
00162 if ( k < 1000000 || y2 <= y )
00163 {
00164 for ( ; ; )
00165 {
00166 if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
00167 z2 * m2 * r * ( z2 * q - r ) )
00168 {
00169 a3 = a2 + r / q;
00170 y3 = f ( a3 );
00171
00172 if ( y3 < y )
00173 {
00174 x = a3;
00175 y = y3;
00176 }
00177 }
00178 k = ( ( 1611 * k ) % 1048576 );
00179 q = 1.0;
00180 r = ( b - a ) * 0.00001 * ( double ) ( k );
00181
00182 if ( z2 <= r )
00183 {
00184 break;
00185 }
00186 }
00187 }
00188 else
00189 {
00190 k = ( ( 1611 * k ) % 1048576 );
00191 q = 1.0;
00192 r = ( b - a ) * 0.00001 * ( double ) ( k );
00193
00194 while ( r < z2 )
00195 {
00196 if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
00197 z2 * m2 * r * ( z2 * q - r ) )
00198 {
00199 a3 = a2 + r / q;
00200 y3 = f ( a3 );
00201
00202 if ( y3 < y )
00203 {
00204 x = a3;
00205 y = y3;
00206 }
00207 }
00208 k = ( ( 1611 * k ) % 1048576 );
00209 q = 1.0;
00210 r = ( b - a ) * 0.00001 * ( double ) ( k );
00211 }
00212 }
00213
00214 r = m2 * d0 * d1 * d2;
00215 s = sqrt ( ( ( y2 - y ) + t ) / m2 );
00216 h = 0.5 * ( 1.0 + h );
00217 p = h * ( p + 2.0 * r * s );
00218 q = q + 0.5 * qs;
00219 r = - 0.5 * ( d0 + ( z0 + 2.01 * e ) / ( d0 * m2 ) );
00220
00221 if ( r < s || d0 < 0.0 )
00222 {
00223 r = a2 + s;
00224 }
00225 else
00226 {
00227 r = a2 + r;
00228 }
00229
00230 if ( 0.0 < p * q )
00231 {
00232 a3 = a2 + p / q;
00233 }
00234 else
00235 {
00236 a3 = r;
00237 }
00238
00239 for ( ; ; )
00240 {
00241 a3 = r8_max ( a3, r );
00242
00243 if ( b <= a3 )
00244 {
00245 a3 = b;
00246 y3 = yb;
00247 }
00248 else
00249 {
00250 y3 = f ( a3 );
00251 }
00252
00253 if ( y3 < y )
00254 {
00255 x = a3;
00256 y = y3;
00257 }
00258
00259 d0 = a3 - a2;
00260
00261 if ( a3 <= r )
00262 {
00263 break;
00264 }
00265
00266 p = 2.0 * ( y2 - y3 ) / ( m * d0 );
00267
00268 if ( ( 1.0 + 9.0 * macheps ) * d0 <= r8_abs ( p ) )
00269 {
00270 break;
00271 }
00272
00273 if ( 0.5 * m2 * ( d0 * d0 + p * p ) <= ( y2 - y ) + ( y3 - y ) + 2.0 * t )
00274 {
00275 break;
00276 }
00277 a3 = 0.5 * ( a2 + a3 );
00278 h = 0.9 * h;
00279 }
00280
00281 if ( b <= a3 )
00282 {
00283 break;
00284 }
00285
00286 a0 = sc;
00287 sc = a2;
00288 a2 = a3;
00289 y0 = y1;
00290 y1 = y2;
00291 y2 = y3;
00292 }
00293
00294 return y;
00295 }
00296
00297
00298 double local_min ( double a, double b, double t, func_base& f,
00299 double &x )
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 {
00365 double c;
00366 double d = 0.0;
00367 double e;
00368 double eps;
00369 double fu;
00370 double fv;
00371 double fw;
00372 double fx;
00373 double m;
00374 double p;
00375 double q;
00376 double r;
00377 double sa;
00378 double sb;
00379 double t2;
00380 double tol;
00381 double u;
00382 double v;
00383 double w;
00384
00385
00386
00387 c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );
00388
00389 eps = sqrt ( r8_epsilon ( ) );
00390
00391 sa = a;
00392 sb = b;
00393 x = sa + c * ( b - a );
00394 w = x;
00395 v = w;
00396 e = 0.0;
00397 fx = f ( x );
00398 fw = fx;
00399 fv = fw;
00400
00401 for ( ; ; )
00402 {
00403 m = 0.5 * ( sa + sb ) ;
00404 tol = eps * r8_abs ( x ) + t;
00405 t2 = 2.0 * tol;
00406
00407
00408
00409 if ( r8_abs ( x - m ) <= t2 - 0.5 * ( sb - sa ) )
00410 {
00411 break;
00412 }
00413
00414
00415
00416 r = 0.0;
00417 q = r;
00418 p = q;
00419
00420 if ( tol < r8_abs ( e ) )
00421 {
00422 r = ( x - w ) * ( fx - fv );
00423 q = ( x - v ) * ( fx - fw );
00424 p = ( x - v ) * q - ( x - w ) * r;
00425 q = 2.0 * ( q - r );
00426 if ( 0.0 < q )
00427 {
00428 p = - p;
00429 }
00430 q = r8_abs ( q );
00431 r = e;
00432 e = d;
00433 }
00434
00435 if ( r8_abs ( p ) < r8_abs ( 0.5 * q * r ) &&
00436 q * ( sa - x ) < p &&
00437 p < q * ( sb - x ) )
00438 {
00439
00440
00441
00442 d = p / q;
00443 u = x + d;
00444
00445
00446
00447 if ( ( u - sa ) < t2 || ( sb - u ) < t2 )
00448 {
00449 if ( x < m )
00450 {
00451 d = tol;
00452 }
00453 else
00454 {
00455 d = - tol;
00456 }
00457 }
00458 }
00459
00460
00461
00462 else
00463 {
00464 if ( x < m )
00465 {
00466 e = sb - x;
00467 }
00468 else
00469 {
00470 e = sa - x;
00471 }
00472 d = c * e;
00473 }
00474
00475
00476
00477 if ( tol <= r8_abs ( d ) )
00478 {
00479 u = x + d;
00480 }
00481 else if ( 0.0 < d )
00482 {
00483 u = x + tol;
00484 }
00485 else
00486 {
00487 u = x - tol;
00488 }
00489
00490 fu = f ( u );
00491
00492
00493
00494 if ( fu <= fx )
00495 {
00496 if ( u < x )
00497 {
00498 sb = x;
00499 }
00500 else
00501 {
00502 sa = x;
00503 }
00504 v = w;
00505 fv = fw;
00506 w = x;
00507 fw = fx;
00508 x = u;
00509 fx = fu;
00510 }
00511 else
00512 {
00513 if ( u < x )
00514 {
00515 sa = u;
00516 }
00517 else
00518 {
00519 sb = u;
00520 }
00521
00522 if ( fu <= fw || w == x )
00523 {
00524 v = w;
00525 fv = fw;
00526 w = u;
00527 fw = fu;
00528 }
00529 else if ( fu <= fv || v == x || v== w )
00530 {
00531 v = u;
00532 fv = fu;
00533 }
00534 }
00535 }
00536 return fx;
00537 }
00538
00539
00540 double local_min_rc ( double &a, double &b, int &status, double value )
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624 {
00625 static double arg;
00626 static double c;
00627 static double d;
00628 static double e;
00629 static double eps;
00630 static double fu;
00631 static double fv;
00632 static double fw;
00633 static double fx;
00634 static double midpoint;
00635 static double p;
00636 static double q;
00637 static double r;
00638 static double tol;
00639 static double tol1;
00640 static double tol2;
00641 static double u;
00642 static double v;
00643 static double w;
00644 static double x;
00645
00646
00647
00648 if ( status == 0 )
00649 {
00650 if ( b <= a )
00651 {
00652 cout << "\n";
00653 cout << "LOCAL_MIN_RC - Fatal error!\n";
00654 cout << " A < B is required, but\n";
00655 cout << " A = " << a << "\n";
00656 cout << " B = " << b << "\n";
00657 status = -1;
00658 exit ( 1 );
00659 }
00660 c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );
00661
00662 eps = sqrt ( r8_epsilon ( ) );
00663 tol = r8_epsilon ( );
00664
00665 v = a + c * ( b - a );
00666 w = v;
00667 x = v;
00668 e = 0.0;
00669
00670 status = 1;
00671 arg = x;
00672
00673 return arg;
00674 }
00675
00676
00677
00678 else if ( status == 1 )
00679 {
00680 fx = value;
00681 fv = fx;
00682 fw = fx;
00683 }
00684
00685
00686
00687 else if ( 2 <= status )
00688 {
00689 fu = value;
00690
00691 if ( fu <= fx )
00692 {
00693 if ( x <= u )
00694 {
00695 a = x;
00696 }
00697 else
00698 {
00699 b = x;
00700 }
00701 v = w;
00702 fv = fw;
00703 w = x;
00704 fw = fx;
00705 x = u;
00706 fx = fu;
00707 }
00708 else
00709 {
00710 if ( u < x )
00711 {
00712 a = u;
00713 }
00714 else
00715 {
00716 b = u;
00717 }
00718
00719 if ( fu <= fw || w == x )
00720 {
00721 v = w;
00722 fv = fw;
00723 w = u;
00724 fw = fu;
00725 }
00726 else if ( fu <= fv || v == x || v == w )
00727 {
00728 v = u;
00729 fv = fu;
00730 }
00731 }
00732 }
00733
00734
00735
00736 midpoint = 0.5 * ( a + b );
00737 tol1 = eps * r8_abs ( x ) + tol / 3.0;
00738 tol2 = 2.0 * tol1;
00739
00740
00741
00742 if ( r8_abs ( x - midpoint ) <= ( tol2 - 0.5 * ( b - a ) ) )
00743 {
00744 status = 0;
00745 return arg;
00746 }
00747
00748
00749
00750 if ( r8_abs ( e ) <= tol1 )
00751 {
00752 if ( midpoint <= x )
00753 {
00754 e = a - x;
00755 }
00756 else
00757 {
00758 e = b - x;
00759 }
00760 d = c * e;
00761 }
00762
00763
00764
00765 else
00766 {
00767 r = ( x - w ) * ( fx - fv );
00768 q = ( x - v ) * ( fx - fw );
00769 p = ( x - v ) * q - ( x - w ) * r;
00770 q = 2.0 * ( q - r );
00771 if ( 0.0 < q )
00772 {
00773 p = - p;
00774 }
00775 q = r8_abs ( q );
00776 r = e;
00777 e = d;
00778
00779
00780
00781 if (
00782 ( r8_abs ( 0.5 * q * r ) <= r8_abs ( p ) ) ||
00783 ( p <= q * ( a - x ) ) ||
00784 ( q * ( b - x ) <= p ) )
00785 {
00786 if ( midpoint <= x )
00787 {
00788 e = a - x;
00789 }
00790 else
00791 {
00792 e = b - x;
00793 }
00794 d = c * e;
00795 }
00796
00797
00798
00799 else
00800 {
00801 d = p / q;
00802 u = x + d;
00803
00804 if ( ( u - a ) < tol2 )
00805 {
00806 d = tol1 * r8_sign ( midpoint - x );
00807 }
00808
00809 if ( ( b - u ) < tol2 )
00810 {
00811 d = tol1 * r8_sign ( midpoint - x );
00812 }
00813 }
00814 }
00815
00816
00817
00818 if ( tol1 <= r8_abs ( d ) )
00819 {
00820 u = x + d;
00821 }
00822 if ( r8_abs ( d ) < tol1 )
00823 {
00824 u = x + tol1 * r8_sign ( d );
00825 }
00826
00827
00828
00829 arg = u;
00830 status = status + 1;
00831
00832 return arg;
00833 }
00834
00835
00836 double r8_abs ( double x )
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862 {
00863 double value;
00864
00865 if ( 0.0 <= x )
00866 {
00867 value = x;
00868 }
00869 else
00870 {
00871 value = - x;
00872 }
00873 return value;
00874 }
00875
00876
00877 double r8_epsilon ( )
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909 {
00910 double r;
00911
00912 r = 1.0;
00913
00914 while ( 1.0 < ( double ) ( 1.0 + r ) )
00915 {
00916 r = r / 2.0;
00917 }
00918
00919 return ( 2.0 * r );
00920 }
00921
00922
00923 double r8_max ( double x, double y )
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949 {
00950 double value;
00951
00952 if ( y < x )
00953 {
00954 value = x;
00955 }
00956 else
00957 {
00958 value = y;
00959 }
00960 return value;
00961 }
00962
00963
00964 double r8_sign ( double x )
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990 {
00991 double value;
00992
00993 if ( x < 0.0 )
00994 {
00995 value = -1.0;
00996 }
00997 else
00998 {
00999 value = 1.0;
01000 }
01001 return value;
01002 }
01003
01004
01005 void timestamp ( )
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033 {
01034 # define TIME_SIZE 40
01035
01036 static char time_buffer[TIME_SIZE];
01037 const struct tm *tm;
01038 time_t now;
01039
01040 now = time ( NULL );
01041 tm = localtime ( &now );
01042
01043 strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
01044
01045 cout << time_buffer << "\n";
01046
01047 return;
01048 # undef TIME_SIZE
01049 }
01050
01051
01052 double zero ( double a, double b, double t, func_base& f )
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105 {
01106 double c;
01107 double d;
01108 double e;
01109 double fa;
01110 double fb;
01111 double fc;
01112 double m;
01113 double macheps;
01114 double p;
01115 double q;
01116 double r;
01117 double s;
01118 double sa;
01119 double sb;
01120 double tol;
01121
01122
01123
01124 sa = a;
01125 sb = b;
01126 fa = f ( sa );
01127 fb = f ( sb );
01128
01129 c = sa;
01130 fc = fa;
01131 e = sb - sa;
01132 d = e;
01133
01134 macheps = r8_epsilon ( );
01135
01136 for ( ; ; )
01137 {
01138 if ( r8_abs ( fc ) < r8_abs ( fb ) )
01139 {
01140 sa = sb;
01141 sb = c;
01142 c = sa;
01143 fa = fb;
01144 fb = fc;
01145 fc = fa;
01146 }
01147
01148 tol = 2.0 * macheps * r8_abs ( sb ) + t;
01149 m = 0.5 * ( c - sb );
01150
01151 if ( r8_abs ( m ) <= tol || fb == 0.0 )
01152 {
01153 break;
01154 }
01155
01156 if ( r8_abs ( e ) < tol || r8_abs ( fa ) <= r8_abs ( fb ) )
01157 {
01158 e = m;
01159 d = e;
01160 }
01161 else
01162 {
01163 s = fb / fa;
01164
01165 if ( sa == c )
01166 {
01167 p = 2.0 * m * s;
01168 q = 1.0 - s;
01169 }
01170 else
01171 {
01172 q = fa / fc;
01173 r = fb / fc;
01174 p = s * ( 2.0 * m * a * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
01175 q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
01176 }
01177
01178 if ( 0.0 < p )
01179 {
01180 q = - q;
01181 }
01182 else
01183 {
01184 p = - p;
01185 }
01186
01187 s = e;
01188 e = d;
01189
01190 if ( 2.0 * p < 3.0 * m * q - r8_abs ( tol * q ) &&
01191 p < r8_abs ( 0.5 * s * q ) )
01192 {
01193 d = p / q;
01194 }
01195 else
01196 {
01197 e = m;
01198 d = e;
01199 }
01200 }
01201 sa = sb;
01202 fa = fb;
01203
01204 if ( tol < r8_abs ( d ) )
01205 {
01206 sb = sb + d;
01207 }
01208 else if ( 0.0 < m )
01209 {
01210 sb = sb + tol;
01211 }
01212 else
01213 {
01214 sb = sb - tol;
01215 }
01216
01217 fb = f ( sb );
01218
01219 if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
01220 {
01221 c = sa;
01222 fc = fa;
01223 e = sb - sa;
01224 d = e;
01225 }
01226 }
01227 return sb;
01228 }
01229
01230
01231 void zero_rc ( double a, double b, double t, double &arg, int &status,
01232 double value )
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295 {
01296 static double c;
01297 static double d;
01298 static double e;
01299 static double fa;
01300 static double fb;
01301 static double fc;
01302 double m;
01303 static double macheps;
01304 double p;
01305 double q;
01306 double r;
01307 double s;
01308 static double sa;
01309 static double sb;
01310 double tol;
01311
01312
01313
01314
01315 if ( status == 0 )
01316 {
01317 macheps = r8_epsilon ( );
01318
01319 sa = a;
01320 sb = b;
01321 e = sb - sa;
01322 d = e;
01323
01324 status = 1;
01325 arg = a;
01326 return;
01327 }
01328
01329
01330
01331
01332 else if ( status == 1 )
01333 {
01334 fa = value;
01335 status = 2;
01336 arg = sb;
01337 return;
01338 }
01339
01340
01341
01342
01343 else if ( status == 2 )
01344 {
01345 fb = value;
01346
01347 if ( 0.0 < fa * fb )
01348 {
01349 status = -1;
01350 return;
01351 }
01352 c = sa;
01353 fc = fa;
01354 }
01355 else
01356 {
01357 fb = value;
01358
01359 if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
01360 {
01361 c = sa;
01362 fc = fa;
01363 e = sb - sa;
01364 d = e;
01365 }
01366 }
01367
01368
01369
01370 if ( r8_abs ( fc ) < r8_abs ( fb ) )
01371 {
01372 sa = sb;
01373 sb = c;
01374 c = sa;
01375 fa = fb;
01376 fb = fc;
01377 fc = fa;
01378 }
01379
01380 tol = 2.0 * macheps * r8_abs ( sb ) + t;
01381 m = 0.5 * ( c - sb );
01382
01383 if ( r8_abs ( m ) <= tol || fb == 0.0 )
01384 {
01385 status = 0;
01386 arg = sb;
01387 return;
01388 }
01389
01390 if ( r8_abs ( e ) < tol || r8_abs ( fa ) <= r8_abs ( fb ) )
01391 {
01392 e = m;
01393 d = e;
01394 }
01395 else
01396 {
01397 s = fb / fa;
01398
01399 if ( sa == c )
01400 {
01401 p = 2.0 * m * s;
01402 q = 1.0 - s;
01403 }
01404 else
01405 {
01406 q = fa / fc;
01407 r = fb / fc;
01408 p = s * ( 2.0 * m * a * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
01409 q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
01410 }
01411
01412 if ( 0.0 < p )
01413 {
01414 q = - q;
01415 }
01416 else
01417 {
01418 p = - p;
01419 }
01420 s = e;
01421 e = d;
01422
01423 if ( 2.0 * p < 3.0 * m * q - r8_abs ( tol * q ) &&
01424 p < r8_abs ( 0.5 * s * q ) )
01425 {
01426 d = p / q;
01427 }
01428 else
01429 {
01430 e = m;
01431 d = e;
01432 }
01433 }
01434
01435 sa = sb;
01436 fa = fb;
01437
01438 if ( tol < r8_abs ( d ) )
01439 {
01440 sb = sb + d;
01441 }
01442 else if ( 0.0 < m )
01443 {
01444 sb = sb + tol;
01445 }
01446 else
01447 {
01448 sb = sb - tol;
01449 }
01450
01451 arg = sb;
01452 status = status + 1;
01453
01454 return;
01455 }
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466 typedef double DoubleOfDouble (double);
01467
01468 class func_wrapper : public func_base {
01469 DoubleOfDouble* func;
01470 public:
01471 func_wrapper(DoubleOfDouble* f) {
01472 func = f;
01473 }
01474 virtual double operator() (double x){
01475 return func(x);
01476 }
01477 };
01478
01479
01480
01481 double glomin ( double a, double b, double c, double m, double e,
01482 double t, double f ( double x ), double &x ){
01483 func_wrapper foo(f);
01484 return glomin(a, b, c, m, e, t, foo, x);
01485 }
01486
01487
01488
01489 double local_min ( double a, double b, double t, double f ( double x ),
01490 double &x ){
01491 func_wrapper foo(f);
01492 return local_min(a, b, t, foo, x);
01493 }
01494
01495
01496
01497 double zero ( double a, double b, double t, double f ( double x ) ){
01498 func_wrapper foo(f);
01499 return zero(a, b, t, foo);
01500 }
01501
01502
01503
01504
01505
01506 double monicPoly::operator()(double x){
01507 double rslt(1);
01508 for (int ii = coeff.size()-1; ii >= 0; ii--){
01509 rslt *= x;
01510 rslt += coeff[ii];
01511 }
01512 return rslt;
01513 }
01514
01515
01516 double Poly::operator()(double x){
01517 double rslt(0);
01518 for (int ii = coeff.size()-1; ii >= 0; ii--){
01519 rslt *= x;
01520 rslt += coeff[ii];
01521 }
01522 return rslt;
01523 }
01524
01525 }