[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[lmi-commits] [lmi] master d5bd712 06/13: Add an instrumented reference
From: |
Greg Chicares |
Subject: |
[lmi-commits] [lmi] master d5bd712 06/13: Add an instrumented reference implementation |
Date: |
Sun, 4 Jul 2021 19:04:44 -0400 (EDT) |
branch: master
commit d5bd7126c935c39c6481199eb2434e4dba8c332b
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>
Add an instrumented reference implementation
---
zero.hpp | 136 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 134 insertions(+), 2 deletions(-)
diff --git a/zero.hpp b/zero.hpp
index c77c6de..da1ae65 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -467,7 +467,7 @@ root_type decimal_root
);
}
-/// A C++ translation of Brent's algol60 reference implementation.
+/// An instrumented translation of Brent's reference implementation.
template<typename FunctionalType>
double brent_zero
@@ -475,6 +475,137 @@ double brent_zero
,double a
,double b
,double t
+ ,std::ostream& os_trace = null_stream()
+ )
+{
+ // Returns a zero of the function f in the given interval [a,b],
+ // to within a tolerance 6ϵ|ζ| + 2t, where ϵ is the relative
+ // machine precision and t is a positive tolerance. Assumes
+ // that f(a) and f(b) have different signs.
+
+ int n_iter {0};
+ interpolation_technique technique {interpolate_initialization};
+
+ os_trace
+ << "#eval"
+ << " a fa"
+ << " b fb"
+ << " c fc"
+ << '\n'
+ ;
+
+ double c, d, e, fa, fb, fc, tol, m, p, q, r, s;
+
+ auto expatiate = [&]()
+ {
+ os_trace
+ << std::setw(3) << n_iter
+ << ' ' << "IBLQb"[technique]
+ << ' ' << std::setw(12) << a << ' ' << std::setw(12) << fa
+ << ' ' << std::setw(12) << b << ' ' << std::setw(12) << fb
+ << ' ' << std::setw(12) << c << ' ' << std::setw(12) << fc
+ << std::endl
+ ;
+ };
+
+ fa = f(a);
+ ++n_iter;
+ fb = f(b);
+ ++n_iter;
+ c = fc = 0.0; // Zero-initialize before calling expatiate().
+ expatiate();
+ interpolate:
+ c = a; fc = fa; d = e = b - a;
+ extrapolate:
+ if(std::fabs(fc) < std::fabs(fb))
+ {
+ a = b; b = c; c = a;
+ fa = fb; fb = fc; fc = fa;
+ }
+ tol = 2.0 * DBL_EPSILON * std::fabs(b) + t;
+ m = 0.5 * (c - b);
+ if(tol < std::fabs(m) && 0.0 != fb)
+ {
+ // See if a bisection is forced.
+ if(std::fabs(e) < tol || std::fabs(fa) <= std::fabs(fb))
+ {
+ technique = interpolate_bisection0;
+ d = e = m;
+ }
+ else
+ {
+ s = fb / fa;
+ 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;
+ 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.0 < p)
+ {
+ q = -q;
+ }
+ else
+ {
+ p = -p;
+ }
+ s = e;
+ e = d;
+ if
+ ( 2.0 * p < 3.0 * m * q - std::fabs(tol * q)
+ && p < std::fabs(0.5 * s * q)
+ )
+ {
+ d = p / q;
+ }
+ else
+ {
+ technique = interpolate_bisection1;
+ d = e = m;
+ }
+ }
+ a = b; fa = fb;
+ if(tol < std::fabs(d))
+ {
+ b += d;
+ }
+ else if(0.0 < m)
+ {
+ b += tol;
+ }
+ else
+ {
+ b -= tol;
+ }
+ fb = f(b);
+ ++n_iter;
+ expatiate();
+ if((0.0 < fb) == (0.0 < fc))
+ {goto interpolate;}
+ else
+ {goto extrapolate;}
+ }
+ return b;
+}
+
+/// A C++ translation of Brent's algol60 reference implementation.
+
+template<typename FunctionalType>
+double brent_zero_reference
+ (FunctionalType& f
+ ,double a
+ ,double b
+ ,double t
)
{
// Returns a zero of the function f in the given interval [a,b],
@@ -482,7 +613,8 @@ double brent_zero
// machine precision and t is a positive tolerance. Assumes
// that f(a) and f(b) have different signs.
double c, d, e, fa, fb, fc, tol, m, p, q, r, s;
- fa = f(a); fb = f(b);
+ fa = f(a);
+ fb = f(b);
interpolate:
c = a; fc = fa; d = e = b - a;
extrapolate:
- [lmi-commits] [lmi] master updated (a03c499 -> 55873f1), Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master bfd60ed 03/13: Rename primary root finder, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master aaadfe1 05/13: Translate the ALGOL scrolls more literally, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master e68be69 07/13: Test lmi_root() against reference implementation, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master 41d1896 08/13: Move some tests out of a hodgepodge, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master 0e29fb4 10/13: Augment a block of unit tests, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master bf65f24 12/13: Remove duplicative tests, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master 28955d1 01/13: Refactor for flexibility, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master 8a1404f 02/13: Describe the preceding commit, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master e6b40ee 04/13: Write root-finder arguments in a more natural order, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master d5bd712 06/13: Add an instrumented reference implementation,
Greg Chicares <=
- [lmi-commits] [lmi] master 9be9b1d 09/13: Refactor, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master a1a1901 11/13: Refactor, Greg Chicares, 2021/07/04
- [lmi-commits] [lmi] master 55873f1 13/13: Test an uncommon code path, Greg Chicares, 2021/07/04