[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[lmi-commits] [lmi] master 27cdaa7 01/11: Refactor
From: |
Greg Chicares |
Subject: |
[lmi-commits] [lmi] master 27cdaa7 01/11: Refactor |
Date: |
Thu, 15 Jul 2021 14:57:10 -0400 (EDT) |
branch: master
commit 27cdaa7abd3366aba8516a54183800bd84107879
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>
Refactor
Brent's method uses bisection in four different circumstances, which
will usefully be distinguished soon in the optional trace.
---
zero.hpp | 31 ++++++++++++++++++++++++++-----
1 file changed, 26 insertions(+), 5 deletions(-)
diff --git a/zero.hpp b/zero.hpp
index a860644..9c8b5a5 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -365,7 +365,12 @@ root_type lmi_root
; // Do nothing.
}
}
- if(std::fabs(e) < tol || std::fabs(fa) <= std::fabs(fb))
+ if(std::fabs(e) < tol)
+ {
+ technique = interpolate_bisection0;
+ d = e = m;
+ }
+ else if(std::fabs(fa) <= std::fabs(fb))
{
technique = interpolate_bisection0;
d = e = m;
@@ -398,10 +403,26 @@ root_type lmi_root
}
s = e;
e = d;
- if
- ( 2.0 * p < 3.0 * m * q - std::fabs(tol * q)
- && p < std::fabs(0.5 * s * q)
- )
+ // Use the criteria in Brent's ALGOL, which differ
+ // slightly from their descriptions in his text.
+ //
+ // AfMWD says on page 51:
+ // "we reject i [i.e., b + p/q] if 2|p| ≥ 3|mq|"
+ // Difference: the ALGOL subtracts tol×|q| [i.e., δ|q|]
+ bool const k0 = 2.0 * p < 3.0 * m * q - std::fabs(tol * q);
+ // AfMWD says on page 50:
+ // "Let e be the value of p/q at the step before the
+ // last one."
+ // (That value is 's', both above and in the ALGOL.)
+ // "If |e| < δ or |p/q| ≥ ½|e| then we do a bisection"
+ // Difference: the ALGOL tests |e| < δ elsewhere
+ bool const k1 = p < std::fabs(0.5 * s * q);
+ // Do not attempt to invert these conditions, e.g.
+ // - if(a < b) x() else y();
+ // + if(b <= a) y() else x();
+ // because NaNs break such reasoning; instead, make sure
+ // the 'else' clause performs bisection.
+ if(k0 && k1)
{
d = p / q;
}
- [lmi-commits] [lmi] master updated (debc275 -> a277ed6), Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master 0090739 03/11: Rename an enumeration, Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master fb3d854 02/11: Refactor instrumented reference implementation, Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master 0d9b36b 05/11: Make enumerators print a useful value directly, Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master 27cdaa7 01/11: Refactor,
Greg Chicares <=
- [lmi-commits] [lmi] master ff2f404 04/11: Use a more informative unit-test macro, Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master 0b44a84 07/11: Reenumerate root-finding activities, Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master dce3edc 08/11: Show variable shifts in optional trace, Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master 5a3bd1d 10/11: Regularize whitespace, Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master 70f1126 09/11: Find a root that coincides with an input bound, Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master b24bcce 06/11: Improve indentation, Greg Chicares, 2021/07/15
- [lmi-commits] [lmi] master a277ed6 11/11: Uniformly test two functions in parallel, Greg Chicares, 2021/07/15