[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[lmi-commits] [lmi] odd/brent 8ccfa28: Test Chandrupatla's root-finding
|
From: |
Greg Chicares |
|
Subject: |
[lmi-commits] [lmi] odd/brent 8ccfa28: Test Chandrupatla's root-finding method |
|
Date: |
Thu, 17 Jun 2021 20:56:03 -0400 (EDT) |
branch: odd/brent
commit 8ccfa282bdebb3fa0f704dc56695c540dc0ae033
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>
Test Chandrupatla's root-finding method
It is claimed that
"Chandrupatlas's method is more efficient than Dekker's and Brent's,
especially for higher order roots"
on page 97 of
Philipp O. J. Scherer
Computational Physics (2nd ed.)
ISBN 978-3-319-00400-6
If that is true, then this alternative ought to improve the speed of lmi
solves, which involve finding a zero of a polynomial of high order:
months_per_year * number_of_years
which is typically 600 and not uncommonly 1200. But does it really?
* chandrupatla.txt: rough pseudocode for Chandrupatla's method as
presented in Scherer's book
* zero3.hpp: the same, implemented in C++
* zero2.hpp: a version of Brent's method, adapted from Scherer's book,
with many errors corrected
* zero.hpp: adapted for this investigation by adding Chandrupatla's
criterion as presented by Scherer
* zero_test.cpp: adapted to test functions in 'zero*.hpp'
Future commits will discuss this similar investigation that refers to
Scherer's book:
https://www.embeddedrelated.com/showarticle/855.php
Interesting graphs may be found here:
https://github.com/SimpleArt/solver/wiki/Methods
Brent has made his book (oddly rotated) available here:
https://maths-people.anu.edu.au/~brent/pd/rpb011i.pdf
The portion relevant here is the fourth chapter. He has also published
this four-page abridgement, which is less daunting than the full book:
https://maths-people.anu.edu.au/~brent/pd/rpb005.pdf
This commit finds a zero of f(x)=x^19 by various methods, labelled thus
in the unit-test output:
- genuine brent: lmi's translation from Brent's ALGOL 60 to C++, which
is believed to be faithful (the number of iterations differs from
Brent's Table 4.1, presumably due to hardware differences)
- brent: Scherer's paraphrase of Brent, with numerous corrections (this
takes fewer iterations than the "genuine" version, presumably due
to gratuitous differences not corrected)
- chand: Chandrupatla's method as given by Scherer (this requires far
fewer iterations, but every one of them is a bisection)
- decimal_root: lmi's production code, which is the "genuine brent"
algorithm with documented changes, especially for rounding (the
number of iterations is slightly lower than "genuine brent", but
that is of no relevance)
Questions to be investigated:
- Does the apparent superiority of Chandrupatla's method to Brent's
in a handful of examples stem from its apparent preference for
bisection, which seems optimal for some of those examples?
- Is the acceptance criterion for inverse quadratic interpolation (IQI)
given by Chandrupatla superior to the one Brent used, such that
grafting it into Brent's algorithm would be an improvement?
---
chandrupatla.txt | 48 +++++++++++++
zero.hpp | 114 ++++++++++++++++++++++++++++---
zero2.hpp | 202 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
zero3.hpp | 144 +++++++++++++++++++++++++++++++++++++++
zero_test.cpp | 39 ++++++++---
5 files changed, 528 insertions(+), 19 deletions(-)
diff --git a/chandrupatla.txt b/chandrupatla.txt
new file mode 100644
index 0000000..d5be528
--- /dev/null
+++ b/chandrupatla.txt
@@ -0,0 +1,48 @@
+void chandrupatla()
+{
+start with interval [x0,x1] such that f(x0) * f(x1) <= 0.0
+b = x0;
+a = c = x1;
+fb = f(b);
+fa = fc = f(a);
+t = 0.5;
+
+for(;;)
+ {
+ xt = a + t * (b - a);
+ ft = f(xt);
+ if(sgn(ft) == sgn(fa))
+ {
+ c = a; fc = fa;
+ a = xt; fa = fx;
+ }
+ else
+ {
+ c = b; fc = fb;
+ b = a; fb = fa;
+ a = xt; fa = ft;
+ }
+ xm = a;
+ fm = fa;
+ if(fabs(fb) < fabs(fa))
+ {
+ xm = b; fm = fb;
+ }
+ tol = 2.0 * epsilon * fabs(xm) + epsilon_a;
+ tl = tol / fabs(b - c);
+ if(tl > 0.5 || 0 == fm) {exit();}
+ xi = (a - b) / (c - b);
+ phi = (fa - fb) / (fc - fb);
+ if(1.0 - sqrt(1 - xi) < phi && phi < sqrt(xi))
+ {
+ t = (fa * fc) / ((fb - fa) * (fb - fc))
+ + (fa * fb) / ((fc - fa) * (fc - fb))
+ * (c - a ) / (b - a)
+ ;
+ }
+ else {t = 0.5;}
+ if(t < tl) {t = tl;}
+ if(t > 1 - tl) {t = 1 - tl;}
+ }
+}
+#endif // zero_hpp
diff --git a/zero.hpp b/zero.hpp
index 205329d..3a58f82 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -29,6 +29,7 @@
#include <cfloat> // DECIMAL_DIG
#include <cmath> // fabs(), pow()
+#include <iostream>
#include <limits>
#include <ostream>
@@ -42,6 +43,7 @@ enum interpolation_technique
{interpolate_initialization
,interpolate_bisection0
,interpolate_linear
+ ,interpolate_linear2
,interpolate_inverse_quadratic
,interpolate_bisection1 // bisection when quadratic rejected
};
@@ -67,13 +69,17 @@ inline void expatiate
,interpolation_technique technique
,double x
,double fx
+ ,bool cond_0
+ ,bool cond_b
+ ,bool cond_c
)
{
if(os_trace.good())
{
os_trace
<< "iteration " << n_iter
- << " " << "IBLQb"[technique]
+ << " " << "IBLcQb"[technique]
+ << " " << cond_0 << cond_b << cond_c
<< " iterand " << x
<< " value " << fx
<< std::endl
@@ -262,7 +268,7 @@ root_type decimal_root
,std::ostream& os_trace = null_stream()
)
{
- constexpr double epsilon {std::numeric_limits<double>::epsilon()};
+ static constexpr double epsilon {std::numeric_limits<double>::epsilon()};
round_to<double> const round_dec {decimals, r_to_nearest};
@@ -282,14 +288,14 @@ root_type decimal_root
}
double fa = static_cast<double>(f(a));
- detail::expatiate(os_trace, n_iter++, technique, a, fa);
+ detail::expatiate(os_trace, n_iter++, technique, a, fa, false, false,
false);
if(0.0 == fa) // Note 0.
{
return {a, root_is_valid, n_iter};
}
double fb = static_cast<double>(f(b));
- detail::expatiate(os_trace, n_iter++, technique, b, fb);
+ detail::expatiate(os_trace, n_iter++, technique, b, fb, false, false,
false);
if(0.0 == fb) // Note 0 [bis].
{
return {b, root_is_valid, n_iter};
@@ -308,18 +314,23 @@ root_type decimal_root
for(;;)
{
+ bool cond_0 = false;
+ bool cond_c = false;
+ bool cond_b = false;
+//os_trace << " old: " << a << ' ' << b << ' ' << c << "\n " << fa << ' '
<< fb << ' ' << fc << std::endl;
if((0.0 < fb) == (0.0 < fc))
{
c = a;
fc = fa;
d = e = b - a;
+//os_trace << " int: " << a << ' ' << b << ' ' << c << "\n " << fa << ' '
<< fb << ' ' << fc << std::endl;
}
- // If 'c' is a closer approximant than 'b', then swap them,
- // discarding the old value of 'a'.
if(std::fabs(fc) < std::fabs(fb))
{
+//os_trace << " was: " << a << ' ' << b << ' ' << c << "\n " << fa << ' '
<< fb << ' ' << fc << std::endl;
a = b; b = c; c = a;
fa = fb; fb = fc; fc = fa;
+//os_trace << " ext: " << a << ' ' << b << ' ' << c << "\n " << fa << ' '
<< fb << ' ' << fc << std::endl;
}
double tol = 2.0 * epsilon * std::fabs(b) + t;
double m = 0.5 * (c - b);
@@ -356,6 +367,9 @@ root_type decimal_root
technique = interpolate_linear;
p = 2.0 * m * s;
q = 1.0 - s;
+ // NEVER USE SECANT ?
+// technique = interpolate_bisection1;
+// d = e = m;
}
else
{
@@ -364,6 +378,21 @@ root_type decimal_root
double r = fb / fc;
p = s * (2.0 * m * q * (q - r) - (b - a) * (r - 1.0));
q = (q - 1.0) * (r - 1.0) * (s - 1.0);
+#if 0
+ double Xi = (b - a) / (c - a);
+ double Phi = (fb - fa) / (fc - fa);
+ cond_c =
+ (Phi * Phi) < Xi
+ && ((1.0 - Phi) * (1.0 - Phi)) < (1.0 - Xi)
+ ;
+ if(!cond_c) // interpolate linearly instead?
+ {
+os_trace << " chandrupatla..." << std::endl;
+ technique = interpolate_linear2;
+ p = 2.0 * m * s;
+ q = 1.0 - s;
+ }
+#endif // 0
}
if(0.0 < p)
{
@@ -375,12 +404,70 @@ root_type decimal_root
}
s = e;
e = d;
+
+// double xi = (a - b) / (c - b);
+// double phi = (fa - fb) / (fc - fb);
+// no...b is the latest, a is the previous, and c the oldest
+// double xi = (b - a) / (c - a);
+// double phi = (fb - fa) / (fc - fa);
+// no...try again...WAIT, ISN'T THIS THE CORRECT ONE?
+// double xi = ( b - c) / ( a - c);
+// double phi = (fb - fc) / (fa - fc);
+// no...try again
+ double xi = (fb - fc) / (fa - fc);
+ double phi = ( b - c) / ( a - c);
+#if 0
+ cond_c =
+ (phi * phi) < xi
+ && ((1.0 - phi) * (1.0 - phi)) < (1.0 - xi)
+// && interpolate_inverse_quadratic == technique
+ ;
+#endif // 0
+ // Chandrupatla acceptance criterion
+ cond_c =
+ (phi * phi) < xi
+ && ((1.0 - phi) * (1.0 - phi)) < (1.0 - xi)
+ ;
+//os_trace << phi << ' ' << xi << std::endl;
+
+ cond_b = p < 1.5 * m * q - std::fabs(tol * q);
+ cond_0 = p < std::fabs(0.5 * s * q);
if
- ( p < 1.5 * m * q - std::fabs(tol * q)
- && p < std::fabs(0.5 * s * q)
+ ( p < std::fabs(0.5 * s * q)
+#if 1
+ && p < 1.5 * m * q - std::fabs(tol * q)
+#else
+ && (interpolate_inverse_quadratic != technique || cond_c)
+#endif // 1
)
{
+//std::cout << cond_c << cond_b << cond_0 << std::endl;
d = p / q;
+ if(interpolate_inverse_quadratic == technique && !cond_c)
+ {
+#if 0
+ os_trace
+ << " chandrupatla would reject this:\n"
+ << " " << b + p / q << " in favor of\n"
+ << " " << b + m
+ << " based on these values:\n"
+ << " " << p << " p\n"
+ << " " << 1.5 * m * q - std::fabs(tol * q) << " 1.5 *
m * q - std::fabs(tol * q)\n"
+ << " " << a << " a\n"
+ << " " << b << " b\n"
+ << " " << c << " c\n"
+ << " " << fa << " fa\n"
+ << " " << fb << " fb\n"
+ << " " << fc << " fc\n"
+ << " " << xi << " xi\n"
+ << " " << phi << " phi\n"
+ << " " << ((phi * phi) < xi) << " (phi * phi) < xi\n"
+ << " " << (((1.0 - phi) * (1.0 - phi)) < (1.0 - xi))
<< " ((1.0 - phi) * (1.0 - phi)) < (1.0 - xi)\n"
+ << std::endl;
+#endif // 0
+ technique = interpolate_bisection1;
+ d = e = m;
+ }
}
else
{
@@ -415,7 +502,7 @@ root_type decimal_root
else
{
fb = static_cast<double>(f(b));
- detail::expatiate(os_trace, n_iter++, technique, b, fb);
+ detail::expatiate(os_trace, n_iter++, technique, b, fb, cond_0,
cond_b, cond_c);
}
}
}
@@ -430,7 +517,9 @@ double brent_zero
,FunctionalType& f
)
{
- constexpr double epsilon {std::numeric_limits<double>::epsilon()};
+ static constexpr double epsilon {std::numeric_limits<double>::epsilon()};
+ interpolation_technique technique {interpolate_initialization};
+std::cout.precision(21);
// Returns a zero of the function f in [a,b] to within a tolerance
// 6 * epsilon * |z| + 2 * t
@@ -442,6 +531,7 @@ double brent_zero
double d = b - a;
double e = d;
+int j = 0;
for(;;)
{
if((0.0 < fb) == (0.0 < fc))
@@ -464,6 +554,7 @@ double brent_zero
if(std::fabs(e) < tol || std::fabs(fa) <= std::fabs(fb))
{
// Bisection.
+ technique = interpolate_bisection0;
d = e = m;
}
else
@@ -473,12 +564,14 @@ double brent_zero
if(a == c)
{
// Linear interpolation.
+ technique = interpolate_linear;
p = 2.0 * m * s;
q = 1.0 - s;
}
else
{
// Inverse quadratic interpolation.
+ technique = interpolate_inverse_quadratic;
q = fa / fc;
double r = fb / fc;
p = s * (2.0 * m * q * (q - r) - (b - a) * (r - 1.0));
@@ -520,6 +613,7 @@ double brent_zero
{
b -= tol;
}
+std::cout << ++j << " ze " << "IBLcQb"[technique] << ' ' << b << std::endl;
fb = f(b);
}
}
diff --git a/zero2.hpp b/zero2.hpp
new file mode 100644
index 0000000..50cdaf9
--- /dev/null
+++ b/zero2.hpp
@@ -0,0 +1,202 @@
+// Root finding by Brent's method.
+//
+// Copyright (C) 2021 Gregory W. Chicares.
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License version 2 as
+// published by the Free Software Foundation.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
+//
+// https://savannah.nongnu.org/projects/lmi
+// email: <gchicares@sbcglobal.net>
+// snail: Chicares, 186 Belle Woods Drive, Glastonbury CT 06033, USA
+
+// Adapted from pseudocode on pages 93-94 of
+// Philipp O. J. Scherer
+// Computational Physics (2nd ed.)
+// ISBN 978-3-319-00400-6
+// and corrected against Brent's original.
+
+#ifndef zero2_hpp
+#define zero2_hpp
+
+#include "config.hpp"
+
+#include "zero.hpp"
+
+#include <cfloat> // DECIMAL_DIG
+#include <cmath> // fabs(), pow()
+#include <iostream>
+#include <limits>
+
+template<typename T>
+inline T signum(T t)
+{
+ return (0 == t) ? 0 : std::signbit(t) ? -1 : 1;
+}
+
+/// A C++ equivalent of Brent's algol60 original, for reference only.
+#if 0
+/// A C++ equivalent of Brent's algol60 original, for reference only.
+
+template<typename FunctionalType>
+double brent_zero
+ (double a
+ ,double b
+ ,double t
+ ,FunctionalType& f
+ )
+#endif // 0
+template<typename FunctionalType>
+double brent
+ (double a
+ ,double b
+ ,double t
+ ,FunctionalType& f
+ )
+{
+ static constexpr double epsilon {std::numeric_limits<double>::epsilon()};
+ interpolation_technique technique {interpolate_initialization};
+std::cout.precision(21);
+
+ double fa = f(a);
+ double fb = f(b);
+#if 0 // RECTIFIED
+ // Brent: at top of loop (not bottom, as here):
+ // if fb and fc have the same sign, then swap a and c
+ // this code seems to do something equivalent
+ double fc = fa; // Brent: = fb
+ double c = a; // likewise
+#endif // 0
+#if 1 // RECTIFIED
+ double fc = fb;
+ double c = b;
+#endif // 1
+
+ double d = b - a;
+ double e = d;
+
+//std::cout << " br " << ' ' << a << ' ' << b << ' ' << c << ' ' << fa << ' '
<< fb << ' ' << fc << std::endl;
+
+int j = 0;
+ for(;;)
+ {
+#if 1 // RECTIFIED
+ // Brent does this at the top
+ if((0.0 < fb) == (0.0 < fc))
+// if(signum(fb) == signum(fc)) // should be equivalent
+ {
+ c = a;
+ fc = fa;
+ d = e = b - a;
+ }
+#endif // 1
+ // if c is a better approximant than b, then exchange
+ // Brent does the same, but it's confusing:
+ // it looks like a three-way shift
+ // but it actually wipes out the old 'a'
+ if(std::fabs(fc) < std::fabs(fb))
+ {
+ a = b; b = c; c = a;
+ fa = fb; fb = fc; fc = fa;
+ }
+ // Brent: set 'tol' carefully--not just epsilon
+// double tol = epsilon; // RECTIFIED
+ double tol = 2.0 * epsilon * std::fabs(b) + t;
+ double m = 0.5 * (c - b);
+ // Brent: '<=' rather than '<'
+// if(0.0 == fb || std::fabs(m) < tol) // RECTIFIED
+ if(0.0 == fb || std::fabs(m) <= tol)
+ {
+ return b;
+ }
+ // bisect if previous iteration took too small a step or
+ // produced no improvement
+ // Brent: '<=' rather than '<' in second condition
+// if(std::fabs(e) < tol || std::fabs(fa) < std::fabs(fb)) // RECTIFIED
+ if(std::fabs(e) < tol || std::fabs(fa) <= std::fabs(fb))
+ {
+ technique = interpolate_bisection0;
+ d = e = m;
+ }
+ // else try interpolating
+ else
+ {
+ double p, q;
+ // Brent: store fb/fa for later use--avoid numerical problems?
+ if(a == c)
+ // linear
+ {
+ technique = interpolate_linear;
+ p = 2.0 * m * fb / fa;
+ // Brent: opposite sign
+// q = (fb - fa) / fa; // PARTLY RECTIFIED
+ q = (fa - fb) / fa;
+ }
+ else
+ // inverse quadratic
+ {
+ technique = interpolate_inverse_quadratic;
+ // Brent: arrange differently--avoid numerical problems?
+ p = 2.0 * m * fb * (fa - fb) / (fc * fc) - (b - a) * fb * (fb
- fc) / (fa * fc);
+ q = (fa / fc - 1.0) * (fb / fc - 1.0) * (fb / fa - 1.0);
+ }
+ // force 'p' to be positive
+ // Brent: '<' rather than '<='
+ if(0.0 <= p)
+ {q = -q;}
+ else
+ {p = -p;}
+ // store previous step
+ double s = e;
+ e = d;
+ if
+ // Brent: absolute value |tol * q|
+ ( 2.0 * p < 3.0 * m * q - tol * q
+ && p < std::fabs(0.5 * s * q)
+ )
+ {
+ d = p / q;
+ }
+ else
+ {
+ d = e = m;
+ }
+ } // closing brace not evident in textbook
+
+ a = b;
+ fa = fb;
+ if(tol < std::fabs(d))
+ {
+ b += d;
+ }
+ else
+ {
+ // Brent: adds 'tol' if 0==m
+ b += tol * signum(m);
+ }
+ fb = f(b);
+//std::cout << ++j << " br " << ' ' << a << ' ' << b << ' ' << c << ' ' <<
"IBLcQb"[technique] << ' ' << b << std::endl;
+std::cout << ++j << " br " << "IBLcQb"[technique] << ' ' << b << std::endl;
+#if 0 // RECTIFIED
+ // Brent does this at the top
+ if((0.0 < fb) == (0.0 < fc))
+// if(signum(fb) == signum(fc)) // should be equivalent
+ {
+ c = a;
+ fc = fa;
+ e = d = b - a;
+ }
+#endif // 0
+ }
+}
+
+#endif // zero2_hpp
diff --git a/zero3.hpp b/zero3.hpp
new file mode 100644
index 0000000..b4c28cc
--- /dev/null
+++ b/zero3.hpp
@@ -0,0 +1,144 @@
+// Root finding by Chandrupatla's method.
+//
+// Copyright (C) 2021 Gregory W. Chicares.
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License version 2 as
+// published by the Free Software Foundation.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
+//
+// https://savannah.nongnu.org/projects/lmi
+// email: <gchicares@sbcglobal.net>
+// snail: Chicares, 186 Belle Woods Drive, Glastonbury CT 06033, USA
+
+// Adapted from pseudocode on pages 95-97 of
+// Philipp O. J. Scherer
+// Computational Physics (2nd ed.)
+// ISBN 978-3-319-00400-6
+// Scherer recommends this method over Brent's, though the evidence he
+// gives is not strong. His graphs do show that it outperforms Brent's
+// method when pure bisection also does.
+
+#ifndef zero3_hpp
+#define zero3_hpp
+
+#include "config.hpp"
+
+#include "zero.hpp"
+
+#include <cfloat> // DECIMAL_DIG
+#include <cmath> // fabs(), pow()
+#include <iostream>
+#include <limits>
+
+#if 0
+template<typename T>
+inline T signum(T t)
+{
+ return (0 == t) ? 0 : std::signbit(t) ? -1 : 1;
+}
+#endif // 0
+
+template<typename FunctionalType>
+double chand
+ (double a
+ ,double b
+ ,double // t
+ ,FunctionalType& f
+ )
+{
+ static constexpr double epsilon {std::numeric_limits<double>::epsilon()};
+ interpolation_technique technique {interpolate_initialization};
+std::cout.precision(21);
+
+ // textbook swaps 'a' and 'b' here compared to its "brent" routine
+ // but that probably doesn't matter
+ double fa = f(a);
+ double fb = f(b);
+ double c = a;
+ double fc = fa;
+ double t = 0.5;
+
+int j = 0;
+ for(;;)
+ {
+ // 't' is coefficient for interpolating (linearly) by any technique
+ double x = a + t * (b - a);
+ double fx = f(x);
+//std::cout << ++j << " ch " << "IBLcQb"[technique] << ' ' << x << std::endl;
+std::cout << ++j << " ch " << "IBLcQb"[technique] << ' ' << t << ' ' << x << '
' << fx << ' ' << a << ' ' << b << ' ' << c << ' ' << b << std::endl;
+ if(signum(fx) == signum(fa))
+ {
+ c = a;
+ fc = fa;
+ a = x;
+ fa = fx; // textbook has weird capitalization on RHS
+ }
+ else
+ {
+ c = b;
+ b = a;
+ a = x;
+ fc = fb;
+ fb = fa;
+ fa = fx;
+ }
+ double xm = a;
+ double fm = fa;
+ if(std::fabs(fb) < std::fabs(fa))
+ {
+ xm = b;
+ fm = fb;
+ }
+ // textbook: 'epsM' and 'epsA' not defined
+ double epsM = epsilon;
+ double epsA = epsilon;
+ double tol = 2.0 * epsM * std::fabs(xm) + epsA;
+// Brent uses this:
+// double tol = 2.0 * epsilon * std::fabs(b) + t;
+ double tl = tol / std::fabs(b - c);
+ if(0.5 < tl || 0.0 == fm)
+ {
+ return xm;
+ }
+// specimen values:
+// 0.0659374999999999961142 a
+// 0.0406250000000000013878 b this is one side of the bracket
+// 0.0912499999999999977796 c
+// 0.0447250871687336973292 true value
+ double xi = (a - b) / (c - b);
+ double phi = (fa - fb) / (fc - fb);
+ if
+ ( 1.0 - std::sqrt(1.0 - xi) < phi
+ && phi < std::sqrt(xi)
+ )
+ {
+ technique = interpolate_inverse_quadratic;
+ // textbook: really use 't'?
+ t =
+ (fa / (fb - fa)) * (fc / (fb - fc))
+ + ((c - a) / (b - a)) * (fa / (fc - fa)) * (fb / (fc - fb))
+ ;
+ }
+ else
+ {
+ technique = interpolate_bisection0;
+ t = 0.5;
+ }
+ // constrain t to (tl, 1-tl)
+ if(t < tl)
+ {t = tl;}
+ if(1.0 - tl < t)
+ {t = 1.0 - tl;}
+ }
+}
+
+#endif // zero3_hpp
diff --git a/zero_test.cpp b/zero_test.cpp
index 176806d..d18d1f0 100644
--- a/zero_test.cpp
+++ b/zero_test.cpp
@@ -22,6 +22,8 @@
#include "pchfile.hpp"
#include "zero.hpp"
+#include "zero2.hpp"
+#include "zero3.hpp"
#include "materially_equal.hpp"
#include "test_tools.hpp"
@@ -111,7 +113,11 @@ struct e_functor
struct e_nineteenth
{
- double operator()(double z) {return std::pow(z, 19);}
+ double operator()(double z) {return std::pow(z, 19);} // -1.0 , 4.0
+// double operator()(double z) {return std::cos(z) - 0.999;} // -0.01, 0.8
+// double operator()(double z) {return std::pow(z - 1.7, 17.0);} // 0.0 , 2.0
+// double operator()(double z) {return std::pow((z - 1.0), 3);} // 0.0 , 1.8
+// double operator()(double z) {return std::pow(z, 2.0) - 2.0;} // 0.0 , 2.0
};
/// A function that's unfriendly to the secant method.
@@ -193,7 +199,7 @@ int test_main(int, char*[])
std::ostringstream oss;
r = decimal_root(0.5, 5.0, bias_none, 9, e_function, oss);
- std::cout << oss.str() << std::endl;
+// std::cout << oss.str() << std::endl;
// Test use with function object.
@@ -269,10 +275,23 @@ int test_main(int, char*[])
// 195 Brent's table 4.1 (IBM 360)
// 171 x86_64 brent_zero (IEEE 754)
// 169 x86_64 decimal_root (differs slightly due to rounding)
- double d = brent_zero(-1.0, 4.0, 1.0e-20, e_19);
+ double lo = -1.00;
+ double hi = 4.0 ;
+std::cout << "test genuine brent" << std::endl;
+ double d = brent_zero(lo, hi, 1.0e-16, e_19);
LMI_TEST(std::fabs(d) <= epsilon);
- r = decimal_root(-1.0, 4.0, bias_none, 20, e_19);
+std::cout << "test brent" << std::endl;
+ d = brent(lo, hi, 1.0e-16, e_19);
+ std::cout << d << std::endl;
+
+std::cout << "test chand" << std::endl;
+ d = chand(lo, hi, 1.0e-16, e_19);
+ std::cout << d << std::endl;
+
+std::cout << "test decimal_root" << std::endl;
+std::ostringstream os2;
+ r = decimal_root(lo, hi, bias_none, 16, e_19, os2);
LMI_TEST(root_is_valid == r.validity);
LMI_TEST(std::fabs(r.root) <= epsilon);
// Assertions labelled 'weak' test the number of iterations
@@ -280,12 +299,14 @@ int test_main(int, char*[])
// rather than a theoretical maximum. Perhaps they'll always
// succeed, because floating-point behavior is determinate;
// but small variations betoken no catastrophe.
- LMI_TEST(169 <= r.n_iter && r.n_iter <= 173); // weak
+std::cout << r.n_iter << std::endl;
+// LMI_TEST(169 <= r.n_iter && r.n_iter <= 173); // weak
+std::cout << os2.str() << std::endl;
- d = brent_zero(-100.0, 100.0, 0.5, eq_2_1);
+// d = brent_zero(-100.0, 100.0, 0.5, eq_2_1);
double zeta = -100.0;
double eq_2_1_upper = zeta + max_err(zeta, 0.5);
- LMI_TEST(-100.0 <= d && d <= eq_2_1_upper);
+// LMI_TEST(-100.0 <= d && d <= eq_2_1_upper);
r = decimal_root(-100.0, 100.0, bias_none, 0, eq_2_1);
LMI_TEST(root_is_valid == r.validity);
@@ -293,7 +314,7 @@ int test_main(int, char*[])
LMI_TEST(10 == max_n_iter_bolzano(-100.0, 100.0, 0.5, -100.0));
LMI_TEST(98 == max_n_iter_brent (-100.0, 100.0, 0.5, -100.0));
LMI_TEST(r.n_iter <= 98);
- LMI_TEST_EQUAL(20, r.n_iter); // weak
+// LMI_TEST_EQUAL(20, r.n_iter); // weak
// Number of iterations required:
// 23 for brent_zero() [above]
// 20 for decimal_root()
@@ -306,7 +327,7 @@ int test_main(int, char*[])
LMI_TEST( 53 == max_n_iter_bolzano(-100.0, 100.0, 0.0, -100.0));
LMI_TEST(2807 == max_n_iter_brent (-100.0, 100.0, 0.0, -100.0));
LMI_TEST(r.n_iter <= 2807);
- LMI_TEST_EQUAL(66, r.n_iter); // weak
+// LMI_TEST_EQUAL(66, r.n_iter); // weak
r = decimal_root(-1.0, 1.0, bias_none, 13, signum_offset);
LMI_TEST(root_is_valid == r.validity);