24 #ifndef FunctionTools_h
25 #define FunctionTools_h
65 double m=coeffLinearTerm, b=coeffConstTerm;
88 double coeffConstTerm) {
89 double A=coeffQuadTerm, B=coeffLinearTerm, C=coeffConstTerm;
98 if (disc < 0)
return roots;
99 q = -0.5*(B + (B < 0 ? -1.0 : 1.0)*sqrt(disc));
101 if (q == 0) roots << q/A;
106 if (!roots.contains(temp)) roots << temp;
123 double coeffLinearTerm,
double coeffConstTerm) {
124 double a=coeffQuadTerm, b=coeffLinearTerm, c=coeffConstTerm;
128 if (coeffCubicTerm == 0.0)
133 if ( coeffCubicTerm != 1.0) {
141 Q = (a*a - 3.0*b) / 9.0;
142 R = (2*a*a*a - 9.0*a*b + 27.0*c) / 54.0;
145 if ( a == 0.0 && b == 0.0 ) {
146 roots << (c < 0 ? 1.0 : -1.0)*pow(fabs(c),1.0/3.0);
148 else if (R*R <= Q*Q*Q) {
149 double theta = acos(R/sqrt(Q*Q*Q)),temp;
152 roots << Q*cos(theta/3.0) - a /3.0;
154 temp = Q*cos((theta + TWOPI)/3.0) - a / 3.0;
155 if (!roots.contains(temp)) roots << temp;
156 temp = Q*cos((theta - TWOPI)/3.0) - a / 3.0;
157 if (!roots.contains(temp)) roots << temp;
161 A = (R < 0 ? 1.0 : -1.0 ) * pow(fabs(R) + sqrt(R*R - Q*Q*Q),1.0/3.0);
162 B = A == 0.0 ? 0.0 : Q / A;
164 roots << (A + B) - a / 3.0;
191 template <
typename Functor>
193 const QList<double> pt2,
double tol,
int maxIter,
double &root) {
194 double a=pt1[0], b=pt2[0], c, d=0, fa = pt1[1], fb = pt2[1], fc, p, q, r, s, t,tol1, bnew, fbnew, temp, deltaI,deltaF;
197 double offset = (a + b) / 2.0;
202 if ( (fa > 0) - (fa < 0) == (fb > 0) - (fb < 0) ) {
203 QString msg =
"The function evaluations of two bounding points passed to Brents Method "
204 "have the same sign. Therefore, they don't necessary bound a root. No "
205 "root finding will be attempted.\n";
212 if (fabs(fa) < fabs(fb)) {
226 for (
int iter=0;iter<maxIter; iter++) {
227 tol1 = DBL_EPSILON*2.0*fabs(b);
228 if (a != c && b != c) {
233 p = s*(t*(r-t)*(c-b)-(1-r)*(b-a));
234 q = (t-1)*(r-1)*(s-1);
239 bnew = b - fb * (b - a) / (fb - fa);
247 p = temp < b ? temp : b;
248 q = b > temp ? b : temp;
249 deltaI = fabs(b - bnew);
251 (bnew < p || bnew > q ) ||
255 (mflag && deltaI >= fabs(b-c)/2.0) ||
260 (!mflag && deltaI >= fabs(c-d)/2.0) ||
265 (mflag && fabs(b-c) < tol1) ||
270 (!mflag && fabs(c-d) < tol1) ) {
280 fbnew = func(bnew + offset);
282 QString msg =
"Function evaluation failed at:" +
toString(bnew) +
283 ". Function must be continuous and defined for the entire interval "
284 "inorder to gaurentee brentsRootFinder will work.";
291 if ( (fa > 0) - (fa < 0) == (fbnew > 0) - (fbnew < 0) ) {
293 deltaF = fabs(a - bnew);
299 deltaF = fabs(b - bnew);
304 if (fabs(fa) < fabs(fb)) {
313 if (fabs(fb) < tol) {
318 else if ( deltaF < tol1 && fabs(b) < 100.0*tol) {
328 else if (deltaF < tol1) {