0027252: Implicit-implicit intersection (Cylinder-Plane) loses intersection curve
[occt.git] / src / math / math_BrentMinimum.cxx
1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15
16 #include <math_BrentMinimum.hxx>
17 #include <math_Function.hxx>
18 #include <StdFail_NotDone.hxx>
19
20 static const Standard_Real CGOLD = 0.3819660; //0.5*(3 - sqrt(5));
21
22
23 //=======================================================================
24 //function : SHFT
25 //purpose  : Shifts arguments
26 //=======================================================================
27 inline void SHFT(Standard_Real &theA, Standard_Real &theB,
28                  Standard_Real &theC, Standard_Real &theD)
29 {
30   theA = theB;
31   theB = theC;
32   theC = theD;
33 }
34
35 //=======================================================================
36 //function : math_BrentMinimum
37 //purpose  : Constructor
38 //=======================================================================
39 math_BrentMinimum::math_BrentMinimum(const Standard_Real    theTolX,
40                                      const Standard_Integer theNbIterations,
41                                      const Standard_Real    theZEPS)
42 : a   (0.0),
43   b   (0.0),
44   x   (0.0),
45   fx  (0.0),
46   fv  (0.0),
47   fw  (0.0),
48   XTol(theTolX),
49   EPSZ(theZEPS),
50   Done   (Standard_False),
51   iter   (0),
52   Itermax(theNbIterations),
53   myF    (Standard_False)
54 {
55 }
56
57 //=======================================================================
58 //function : math_BrentMinimum
59 //purpose  : Constructor
60 //=======================================================================
61 math_BrentMinimum::math_BrentMinimum(const Standard_Real    theTolX,
62                                      const Standard_Real    theFbx,
63                                      const Standard_Integer theNbIterations,
64                                      const Standard_Real    theZEPS)
65 : a   (0.0),
66   b   (0.0),
67   x   (0.0),
68   fx  (theFbx),
69   fv  (0.0),
70   fw  (0.0),
71   XTol(theTolX),
72   EPSZ(theZEPS),
73   Done   (Standard_False),
74   iter   (0),
75   Itermax(theNbIterations),
76   myF    (Standard_True)
77 {
78 }
79
80 //=======================================================================
81 //function : ~math_BrentMinimum
82 //purpose  : Destructor
83 //=======================================================================
84 math_BrentMinimum::~math_BrentMinimum()
85 {
86 }
87
88 //=======================================================================
89 //function : Perform
90 //purpose  : 
91 //=======================================================================
92 void math_BrentMinimum::Perform(math_Function& F,
93                                 const Standard_Real    ax,
94                                 const Standard_Real    bx,
95                                 const Standard_Real    cx)
96 {
97   Standard_Boolean OK;  
98   Standard_Real etemp, fu, p, q, r;  
99   Standard_Real tol1, tol2, u, v, w, xm;
100   Standard_Real e = 0.0;
101   Standard_Real d = RealLast();
102
103   a = ((ax < cx) ? ax : cx);
104   b = ((ax > cx) ? ax : cx);
105   x = w = v = bx;
106   if (!myF) {
107     OK = F.Value(x, fx);
108     if (!OK) return;
109   }
110   fw = fv = fx;
111   for (iter = 1; iter <= Itermax; iter++) {
112     xm = 0.5 * (a + b);
113     tol1 = XTol * fabs(x) + EPSZ;
114     tol2 = 2.0 * tol1;
115     if (IsSolutionReached(F)) {
116       Done = Standard_True;
117       return;
118     }
119     if (fabs(e) > tol1) {
120       r = (x - w) * (fx - fv);
121       q = (x - v) * (fx - fw);
122       p = (x - v) * q - (x - w) * r;
123       q = 2.0 * (q - r);
124       if (q > 0.0) p = -p;
125       q = fabs(q);
126       etemp = e;
127       e = d;
128       if (fabs(p) >= fabs(0.5 * q * etemp)
129         || p <= q * (a - x) || p >= q * (b - x)) {
130         e = (x >= xm ? a - x : b - x);
131         d = CGOLD * e;
132       }
133       else {
134         d = p / q;
135         u = x + d;
136         if (u - a < tol2 || b - u < tol2) d = Sign(tol1, xm - x);
137       }
138     }
139     else {
140       e = (x >= xm ? a - x : b - x);
141       d = CGOLD * e;
142     }
143     u = (fabs(d) >= tol1 ? x + d : x + Sign(tol1, d));
144     OK = F.Value(u, fu);
145     if (!OK) return;
146     if (fu <= fx) {
147       if (u >= x) a = x; else b = x;
148       SHFT(v, w, x, u);
149       SHFT(fv, fw, fx, fu);
150     }
151     else {
152       if (u < x) a = u; else b = u;
153       if (fu <= fw || w == x) {
154         v = w;
155         w = u;
156         fv = fw;
157         fw = fu;
158       }
159       else if (fu <= fv || v == x || v == w) {
160         v = u;
161         fv = fu;
162       }
163     }
164   }
165   Done = Standard_False;
166   return;
167 }
168
169 //=======================================================================
170 //function : Dump
171 //purpose  : 
172 //=======================================================================
173 void math_BrentMinimum::Dump(Standard_OStream& o) const
174 {
175   o << "math_BrentMinimum ";
176   if(Done) {
177     o << " Status = Done \n";
178     o << " Location value = " << x <<"\n";
179     o << " Minimum value = " << fx << "\n";
180     o << " Number of iterations = " << iter <<"\n";
181   }
182   else {
183     o << " Status = not Done \n";
184   }
185 }