From 114f0eff2724ccc05ca7061cea76a24559537dfc Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Sun, 11 Jul 2021 13:14:59 +0200
Subject: [PATCH] Root Finding: Newton-Raphson to check for zero after each
 eval to avoid NaN. Closes #774.

---
 src/Numerics/RootFinding/NewtonRaphson.cs       |  5 +++++
 src/Numerics/RootFinding/RobustNewtonRaphson.cs | 15 +++++++++++++++
 2 files changed, 20 insertions(+)

diff --git a/src/Numerics/RootFinding/NewtonRaphson.cs b/src/Numerics/RootFinding/NewtonRaphson.cs
index 9d5e7e67..ba8650c9 100644
--- a/src/Numerics/RootFinding/NewtonRaphson.cs
+++ b/src/Numerics/RootFinding/NewtonRaphson.cs
@@ -99,6 +99,11 @@ namespace MathNet.Numerics.RootFinding
             {
                 // Evaluation
                 double fx = f(root);
+                if (fx == 0.0)
+                {
+                    return true;
+                }
+
                 double dfx = df(root);
 
                 // Netwon-Raphson step
diff --git a/src/Numerics/RootFinding/RobustNewtonRaphson.cs b/src/Numerics/RootFinding/RobustNewtonRaphson.cs
index 02e1674c..3135b0b3 100644
--- a/src/Numerics/RootFinding/RobustNewtonRaphson.cs
+++ b/src/Numerics/RootFinding/RobustNewtonRaphson.cs
@@ -91,6 +91,11 @@ namespace MathNet.Numerics.RootFinding
 
             root = 0.5*(lowerBound + upperBound);
             double fx = f(root);
+            if (fx == 0.0)
+            {
+                return true;
+            }
+
             double lastStep = Math.Abs(upperBound - lowerBound);
             for (int i = 0; i < maxIterations; i++)
             {
@@ -119,6 +124,11 @@ namespace MathNet.Numerics.RootFinding
                     // Bisection
                     root = 0.5*(upperBound + lowerBound);
                     fx = f(root);
+                    if (fx == 0.0)
+                    {
+                        return true;
+                    }
+
                     lastStep = 0.5*Math.Abs(upperBound - lowerBound);
                     if (Math.Sign(fx) == Math.Sign(fmin))
                     {
@@ -146,6 +156,11 @@ namespace MathNet.Numerics.RootFinding
 
                 // Evaluation
                 fx = f(root);
+                if (fx == 0.0)
+                {
+                    return true;
+                }
+
                 lastStep = step;
 
                 // Update bounds