lmi-commits
[Top][All Lists]

## [lmi-commits] [lmi] master a041329 08/23: Consider bounds not to bracket

 From: Greg Chicares Subject: [lmi-commits] [lmi] master a041329 08/23: Consider bounds not to bracket a root if value at either is NaN Date: Tue, 27 Jul 2021 21:59:51 -0400 (EDT)

```branch: master
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>

Consider bounds not to bracket a root if value at either is NaN
---
zero.hpp      |  6 +++---
zero_test.cpp | 35 ++++++++++++++++++++++++-----------
2 files changed, 27 insertions(+), 14 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index f602b1b..1f26878 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -29,7 +29,7 @@
#include "round_to.hpp"

#include <cfloat>                       // DBL_EPSILON, DECIMAL_DIG
-#include <cmath>                        // fabs(), isfinite(), pow()
+#include <cmath>                        // fabs(), isfinite(), isnan(), pow()
#include <cstdint>                      // uint64_t
#include <cstring>                      // memcpy()
#include <functional>                   // function(), identity()
@@ -465,8 +465,8 @@ root_type lmi_root
return {b, root_is_valid, n_iter};
}

-    // f(a) and f(b) must have different signs.
-    if((0.0 < fa) == (0.0 < fb))
+    // f(a) and f(b) must have different signs; neither may be a NaN.
+    if(std::isnan(fa) || std::isnan(fb) || (0.0 < fa) == (0.0 < fb))
{
recapitulate();
return {0.0, root_not_bracketed, n_iter};
diff --git a/zero_test.cpp b/zero_test.cpp
index 1350fed..b6fec8b 100644
--- a/zero_test.cpp
+++ b/zero_test.cpp
@@ -517,6 +517,15 @@ void test_NaNs()
int const max_n_iter = max_n_iter_brent(-1.0e100, 1.0e100, 5.0e-1, pi);
LMI_TEST_RELATION(r.n_iter,<=,max_n_iter);
LMI_TEST(materially_equal(1.0e100, std::fabs(r.root)));
+
+    // If the function's value is a NaN at either input bound, then
+    // no root is bracketed.
+    r = lmi_root(NaN_signed,  -1.0e100, 2.0 * pi, 5.0e-1);
+    LMI_TEST_EQUAL(root_not_bracketed, r.validity);
+    r = lmi_root(NaN_signed, -2.0 * pi,  1.0e100, 5.0e-1);
+    LMI_TEST_EQUAL(root_not_bracketed, r.validity);
+    r = lmi_root(NaN_signed, -2.0 * pi, 2.0 * pi, 5.0e-1);
+    LMI_TEST_EQUAL(root_not_bracketed, r.validity);
}

/// Find a root that coincides with one or both bounds.
@@ -658,17 +667,21 @@ void test_biases()

// Various tests--see function-template definition.

-    test_bias(-1.0e100, 4.0e100, -100, e, std::exp(1.0));
-    test_bias(-1.0    , 4.0    ,    0, e, std::exp(1.0));
-    test_bias( 0.5    , 5.0    ,    1, e, std::exp(1.0));
-    test_bias( 0.5    , 5.0    ,    2, e, std::exp(1.0));
-    test_bias( 0.5    , 5.0    ,    3, e, std::exp(1.0));
-    test_bias( 0.5    , 5.0    ,    4, e, std::exp(1.0));
-    test_bias( 0.5    , 5.0    ,    5, e, std::exp(1.0));
-    test_bias( 0.5    , 5.0    ,    6, e, std::exp(1.0));
-    test_bias( 0.5    , 5.0    ,    7, e, std::exp(1.0));
-    test_bias( 0.5    , 5.0    ,    8, e, std::exp(1.0));
-    test_bias(-1.0    , 4.0    ,  100, e, std::exp(1.0));
+    // Rounding to -100 decimals makes the maximum error 1e+100,
+    // which probably isn't useful in practice.
+    test_bias(0.0, 4.0e100, -100, e, std::exp(1.0));
+    test_bias(0.0, 4.0    ,    0, e, std::exp(1.0));
+    test_bias(0.5, 5.0    ,    1, e, std::exp(1.0));
+    test_bias(0.5, 5.0    ,    2, e, std::exp(1.0));
+    test_bias(0.5, 5.0    ,    3, e, std::exp(1.0));
+    test_bias(0.5, 5.0    ,    4, e, std::exp(1.0));
+    test_bias(0.5, 5.0    ,    5, e, std::exp(1.0));
+    test_bias(0.5, 5.0    ,    6, e, std::exp(1.0));
+    test_bias(0.5, 5.0    ,    7, e, std::exp(1.0));
+    test_bias(0.5, 5.0    ,    8, e, std::exp(1.0));
+    // Rounding to 100 decimals shouldn't round at all; the
+    // effective maximum error is 6ϵ × e = 3.62148e-15 .
+    test_bias(0.0, 4.0    ,  100, e, std::exp(1.0));
}

/// Test the worked-out example given here:

```