Add ModAB root finding algorithm

.gitignore

@@ -50,3 +50,4 @@ documentation/doxygen/*.dot
 
 # Pycache
 __pycache__
+/out

math/mathcore/CMakeLists.txt

@@ -32,6 +32,7 @@ set(HEADERS
   Math/BrentMethods.h
   Math/BrentMinimizer1D.h
   Math/BrentRootFinder.h
+  Math/ModABRootFinder.h
   Math/ChebyshevPol.h
   Math/Delaunay2D.h
   Math/DistFuncMathCore.h
@@ -125,6 +126,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(MathCore
     src/BrentMethods.cxx
     src/BrentMinimizer1D.cxx
     src/BrentRootFinder.cxx
+    src/ModABRootFinder.cxx
     src/ChebyshevPol.cxx
     src/DataRange.cxx
     src/Delaunay2D.cxx

math/mathcore/inc/Math/ModABRootFinder.h

@@ -0,0 +1,98 @@
+// @(#)root/mathcore:$Id$
+// Authors: Nedelcho Ganchovski    03/03/2026
+
+/**********************************************************************
+ *                                                                    *
+ * Copyright (c) 2026  CERN                                           *
+ * All rights reserved.                                               *
+ *                                                                    *
+ * For the licensing terms see $ROOTSYS/LICENSE.                      *
+ * For the list of contributors see $ROOTSYS/README/CREDITS.          *
+ *                                                                    *
+ **********************************************************************/
+
+// Header for the RootFinder
+//
+// Created by: Nedelcho Ganchovski  : Wed March 03 2026
+//
+
+#ifndef ROOT_Math_ModABRootFinder
+#define ROOT_Math_ModABRootFinder
+
+#include "Math/IFunction.h"
+
+#include "Math/IRootFinderMethod.h"
+
+namespace ROOT {
+namespace Math {
+
+//___________________________________________________________________________________________
+/**
+   Class for finding the root of a one dimensional function using the ModAB algorithm.
+   It is based on the Modified Anderson-Björck method (2022 Ganchovski & Traykov) that
+   adaptively switches between bisection and regula-falsi with
+   side-correction, yielding superlinear convergence on well-behaved
+   functions while retaining the robustness of bisection.
+
+   @ingroup RootFinders
+
+ */
+
+class ModABRootFinder : public IRootFinderMethod {
+public:
+   /** Default Constructor. */
+   ModABRootFinder();
+
+   /** Default Destructor. */
+   ~ModABRootFinder() override {}
+
+   /** Set function to solve and the interval in where to look for the root.
+
+       \@param f Function to be minimized.
+       \@param xlow Lower bound of the search interval.
+       \@param xup Upper bound of the search interval.
+   */
+   using IRootFinderMethod::SetFunction;
+   bool SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup) override;
+
+   /** Returns the X value corresponding to the function value fy for (xmin<x<xmax).
+       Method:
+       Modified Anderson-Björck method is applied on the bracketed interval.
+
+       \@param maxIter maximum number of iterations.
+       \@param absTol desired absolute error in the minimum position.
+       \@param absTol desired relative error in the minimum position.
+   */
+   bool Solve(int maxIter = 100, double absTol = 1E-8, double relTol = 1E-10) override;
+
+   /** Returns root value. Need to call first Solve(). */
+   double Root() const override { return fRoot; }
+
+   /** Returns status of last estimate. If = 0 is OK */
+   int Status() const override { return fStatus; }
+
+   /** Return number of iteration used to find minimum */
+   int Iterations() const override { return fNIter; }
+
+   /** Return name of root finder algorithm ("ModABRootFinder"). */
+   const char *Name() const override;
+
+   // static function used to modify the default parameters
+
+   /** set number of default Npx used at construction time (when SetNpx is not called)
+       Default value is 100
+    */
+
+private:
+   const IGenFunction *fFunction; ///< Pointer to the function.
+   int fNIter;                    ///< Number of iterations needed for the last estimation.
+   int fStatus;                   ///< Status of code of the last estimate
+   double fXMin;                  ///< Lower bound of the search interval.
+   double fXMax;                  ///< Upper bound of the search interval
+   double fRoot;                  ///< Current estimation of the function root.
+};
+
+} // namespace Math
+} // namespace ROOT
+
+#endif /* ROOT_Math_ModABRootFinder */

math/mathcore/inc/Math/RootFinder.h

@@ -75,6 +75,7 @@ namespace ROOT {
       public:
 
          enum EType { kBRENT,                                   // Methods from MathCore
+                     kMODAB,                                    // Modified A&B method added in MathCore
                      kGSL_BISECTION, kGSL_FALSE_POS, kGSL_BRENT, // GSL Normal
                      kGSL_NEWTON, kGSL_SECANT, kGSL_STEFFENSON   // GSL Derivatives
          };

math/mathcore/src/BrentRootFinder.cxx

@@ -22,8 +22,8 @@ namespace ROOT {
 namespace Math {
 
 
-static int gDefaultNpx = 100; // default nunmber of points used in the grid to bracked the root
-static int gDefaultNSearch = 10;  // number of time the iteration (bracketing -Brent ) is repeted
+static int gDefaultNpx = 100; // default number of points used in the grid to bracket the root
+static int gDefaultNSearch = 10;  // number of time the iteration (bracketing -Brent ) is repeated
 
    BrentRootFinder::BrentRootFinder() : fFunction(nullptr),
                                         fLogScan(false), fNIter(0),
@@ -47,15 +47,14 @@ bool BrentRootFinder::SetFunction(const ROOT::Math::IGenFunction& f, double xlow
    // invalid previous status
    fStatus = -1;
 
-   if (xlow >= xup)
+   if (xlow > xup)
    {
-      double tmp = xlow;
-      xlow = xup;
-      xup = tmp;
+      fXMin = xup;
+      fXMax = xlow;
+   } else {
+      fXMin = xlow;
+      fXMax = xup;
    }
-   fXMin = xlow;
-   fXMax = xup;
-
    return true;
 }
 

math/mathcore/src/ModABRootFinder.cxx

@@ -0,0 +1,158 @@
+// @(#)root/mathcore:$Id$
+// Authors: Nedelcho Ganchovski    03/03/2026
+
+/**********************************************************************
+ *                                                                    *
+ * Copyright (c) 2026  CERN                                           *
+ * All rights reserved.                                               *
+ *                                                                    *
+ * For the licensing terms see $ROOTSYS/LICENSE.                      *
+ * For the list of contributors see $ROOTSYS/README/CREDITS.          *
+ *                                                                    *
+ **********************************************************************/
+
+#include "Math/ModABRootFinder.h"
+#include "Math/IFunctionfwd.h"
+#include <cmath>
+
+#include "Math/Error.h"
+
+namespace ROOT {
+namespace Math {
+
+ModABRootFinder::ModABRootFinder()
+   : fFunction(nullptr), fNIter(0), fStatus(-1), fXMin(0), fXMax(0), fRoot(0)
+{
+
+}
+
+bool ModABRootFinder::SetFunction(const ROOT::Math::IGenFunction &f, double xmin, double xmax)
+{
+   // Set function to solve and the interval in where to look for the root.
+
+   fFunction = &f;
+   // invalid previous status
+   fStatus = -1;
+   if (xmin > xmax) {
+      fXMin = xmax;
+      fXMax = xmin;
+   }
+   else {
+      fXMin = xmin;
+      fXMax = xmax;
+   }
+   return true;
+}
+
+const char *ModABRootFinder::Name() const
+{
+   return "ModABRootFinder";
+}
+
+int sign(double x) {
+   return (x > 0) - (x < 0);
+}
+
+bool ModABRootFinder::Solve(int maxIter, double absTol, double relTol)
+{
+   // Returns the X value corresponding to the function value fy for (xmin<x<xmax).
+
+   if (!fFunction) {
+      MATH_ERROR_MSG("ModABRootFinder::Solve", "Function has not been set");
+      return false;
+   }
+   double x1 = fXMin, y1 = (*fFunction)(x1);   
+   if (y1 == 0.0) {
+      fRoot = x1;
+      fStatus = 0;
+      return true;
+   }
+   double x2 = fXMax, y2 = (*fFunction)(x2);   
+   if (y2 == 0.0) {
+      fRoot = x2;
+      fStatus = 0;
+      return true;
+   }
+   if (sign(y1) == sign(y2)) {
+      MATH_ERROR_MSG("ModABRootFinder::Solve", "Function values at the interval endpoints have the same sign");
+      return false;
+   }
+   int side = 0;
+   bool bisection = true;
+   double threshold = x2 - x1;
+   const double C = 16;
+   for (int i = 1; i <= maxIter; ++i) {
+      double x3, y3;
+      if (bisection) {
+         x3 = (x1 + x2) / 2.0;
+         y3 = (*fFunction)(x3);
+         double ym = (y1 + y2) / 2.0;
+         double r = 1 - std::fabs(ym / (y2 - y1));
+         double k = r * r;
+         if (std::fabs(ym - y3) < k * (std::fabs(y3) + std::fabs(ym))) {
+            bisection = false;
+            threshold = (x2 - x1) * C;
+         }
+      } else {
+         x3 = (x1 * y2 - y1 * x2) / (y2 - y1);
+         // Clamp x3 when floating point round-off errors shoots it out of the bracketing interval. Assign also y values to avoid redundant re-evaluation
+         if (x3 <= x1)
+         {
+             x3 = x1;
+             y3 = y1;
+         }
+         else if (x3 >= x2)
+         {
+             x3 = x2;
+             y3 = y2;
+         }
+         else
+             y3 = (*fFunction)(x3);
+
+         threshold /= 2.0;
+      }
+      fRoot = x3;
+      fNIter = i;
+      double eps = absTol + relTol * std::abs(x3);
+      if (std::fabs(y3) == 0.0 || x2 - x1 <= eps) {
+         fStatus = 0;
+         return true;
+      }
+      if (y1 * y3 > 0.0) {
+         if (side == 1) {
+            double m = 1 - y3 / y1;
+            if (m <= 0)
+               y2 /= 2;
+            else
+               y2 *= m;
+         } else if (!bisection)
+            side = 1;
+
+         x1 = x3;
+         y1 = y3;
+      } else {
+         if (side == -1) {
+            double m = 1 - y3 / y2;
+            if (m <= 0)
+               y1 /= 2;
+            else
+               y1 *= m;
+         } else if (!bisection)
+            side = -1;
+
+         x2 = x3;
+         y2 = y3;
+      }
+      if (x2 - x1 > threshold)
+      {
+          bisection = true;
+          side = 0;
+      }
+   }
+   fNIter = maxIter;
+   MATH_ERROR_MSG("ModABRootFinder::Solve", "Search didn't converge");
+   fStatus = -2;
+   return false;
+}
+} // namespace Math
+} // namespace ROOT

math/mathcore/src/RootFinder.cxx

@@ -14,6 +14,7 @@
 #include "Math/RootFinder.h"
 #include "Math/IRootFinderMethod.h"
 #include "Math/BrentRootFinder.h"
+#include "Math/ModABRootFinder.h"
 
 #include "RConfigure.h"
 
@@ -57,6 +58,12 @@ bool RootFinder::SetMethod(RootFinder::EType type)
       return true;
    }
 
+   if ( type == RootFinder::kMODAB )
+   {
+      fSolver = new ModABRootFinder();
+      return true;
+   }
+
 #ifdef MATH_NO_PLUGIN_MANAGER    // no PM available
 #ifdef R__HAS_MATHMORE
 

math/mathcore/test/testRootFinder.cxx

@@ -287,6 +287,12 @@ TEST(Parabola,BrentRootFinder)
    runTestBrent(0);      
 }
 
+TEST(Parabola,ModABRootFinder)
+{
+   // test parabola using ModAB Root Finder
+   runTestBrent(0, ROOT::Math::RootFinder::kMODAB);
+}
+
 #ifdef R__HAS_MATHMORE
 TEST(Parabola, GSLBrent)
 {
@@ -318,8 +324,14 @@ TEST(LogParabola,TF1GetX)
 
 TEST(LogParabola,BrentRootFinder)
 {
-   // test parabola using Brent Root Finder 
-   runTestBrent(0);      
+   // test log-parabola using Brent Root Finder 
+   runTestBrent(0);  //Brent fails with runTestBrent(1)!!! Changed back to 0  
+}
+
+TEST(LogParabola,ModABRootFinder)
+{
+   // test log-parabola using ModAB Root Finder
+   runTestBrent(1, ROOT::Math::RootFinder::kMODAB);
 }
 
 #ifdef R__HAS_MATHMORE
@@ -345,6 +357,12 @@ TEST(GammaCDF,BrentRootFinder)
    runTestBrent(2);      
 }
 
+TEST(GammaCDF, ModABRootFinder)
+{
+   // test gamma cdf using ModAB Root Finder
+   runTestBrent(2, ROOT::Math::RootFinder::kMODAB);
+}
+
 #ifdef R__HAS_MATHMORE
 TEST(GammaCDF,GSLBrent)
 {