EIC code displayed by LXR master/​include/​opencascade/​MathRoot_MultipleUtils.hxx	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-09-28 09:20:52

 // Copyright (c) 2025 OPEN CASCADE SAS
 //
 // This file is part of Open CASCADE Technology software library.
 //
 // This library is free software; you can redistribute it and/or modify it under
 // the terms of the GNU Lesser General Public License version 2.1 as published
 // by the Free Software Foundation, with special exception defined in the file
 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
 // distribution for complete text of the license and disclaimer of any warranty.
 //
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
 #ifndef _MathRoot_MultipleUtils_HeaderFile
 #define _MathRoot_MultipleUtils_HeaderFile
 
 #include <MathUtils_Types.hxx>
 #include <MathUtils_Config.hxx>
 #include <MathRoot_Brent.hxx>
 #include <math_Vector.hxx>
 
 #include <NCollection_DynamicArray.hxx>
 
 #include <algorithm>
 #include <cmath>
 
 //! @file MathRoot_MultipleUtils.hxx
 //! @brief Internal utilities for FindAllRoots / FindAllRootsWithDerivative.
 //!
 //! Contains result/config types, helper functions, functor adapters and the
 //! shared core implementation used by MathRoot_Multiple.hxx.
 
 namespace MathRoot
 {
 using namespace MathUtils;
 
 // ============================================================================
 //  Result and configuration types
 // ============================================================================
 
 //! Result for multiple root finding.
 //! Contains all found roots sorted in ascending order.
 struct MultipleResult
 {
   MathUtils::Status                Status = MathUtils::Status::NotConverged; //!< Computation status
   size_t                           NbIterations = 0; //!< Total iterations across all roots
   NCollection_DynamicArray<double> Roots;            //!< Found roots (sorted)
   NCollection_DynamicArray<double> Values;           //!< Function values at roots
   bool IsAllNull = false; //!< True if function is essentially zero in range
 
   //! Returns true if computation succeeded.
   bool IsDone() const { return Status == MathUtils::Status::OK; }
 
   //! Conversion to bool for convenient checking.
   explicit operator bool() const { return IsDone(); }
 
   //! Returns the number of roots found.
   int NbRoots() const { return Roots.Length(); }
 
   //! Access root by index (0-based).
   double operator[](int theIndex) const { return Roots.Value(theIndex); }
 };
 
 //! Configuration for multiple root finding.
 struct MultipleConfig
 {
   int    NbSamples     = 100;   //!< Number of sample points for initial search
   double XTolerance    = 1e-10; //!< Tolerance on X for convergence
   double FTolerance    = 1e-10; //!< Tolerance on F(X) for convergence
   double NullTolerance = 1e-12; //!< Tolerance to consider function as null
   int    MaxIterations = 100;   //!< Max iterations per root refinement
   double Offset        = 0.0;   //!< Find roots of f(x) - Offset = 0
 };
 
 // ============================================================================
 //  Helper functions
 // ============================================================================
 
 //! In-place insertion sort of roots and corresponding values by ascending root value.
 inline void SortRoots(MultipleResult& theResult)
 {
   for (int i = 1; i < theResult.Roots.Length(); ++i)
   {
     const double aKeyRoot = theResult.Roots[i];
     const double aKeyVal  = theResult.Values[i];
     int          j        = i - 1;
     while (j >= 0 && theResult.Roots[j] > aKeyRoot)
     {
       theResult.Roots[j + 1]  = theResult.Roots[j];
       theResult.Values[j + 1] = theResult.Values[j];
       --j;
     }
     theResult.Roots[j + 1]  = aKeyRoot;
     theResult.Values[j + 1] = aKeyVal;
   }
 }
 
 //! Helper to add a root if it is not a duplicate of an already found root.
 inline void AddRoot(MultipleResult& theResult, double theEpsX, double theRoot, double theValue)
 {
   for (int k = 0; k < theResult.Roots.Length(); ++k)
   {
     if (std::abs(theRoot - theResult.Roots.Value(k)) < theEpsX)
     {
       return;
     }
   }
   theResult.Roots.Append(theRoot);
   theResult.Values.Append(theValue);
 }
 
 //! Compute the minimal X tolerance compatible with the established multi-root behavior.
 inline double EffectiveXTolerance(double theLower, double theUpper, double theXTolerance)
 {
   const double aMinEpsX = 1.0e-10 * (std::abs(theLower) + std::abs(theUpper));
   return std::max(theXTolerance, aMinEpsX);
 }
 
 // ============================================================================
 //  Functor adapters for Value-only interface
 // ============================================================================
 
 //! Samples a Value-only function and stores f(x)-offset into a math_Vector.
 //! @tparam Function type with Value(double theX, double& theF) method
 template <typename Function>
 struct MultipleSampleValueFn
 {
   Function&    myFunc;
   math_Vector& mySamples;
   const double myOffset;
 
   bool operator()(int theIndex, double theX) const
   {
     double aF = 0.0;
     if (!myFunc.Value(theX, aF))
       return false;
     mySamples(theIndex) = aF - myOffset;
     return true;
   }
 };
 
 //! Returns the sampled value at a given index from a math_Vector.
 struct MultipleGetValueFn
 {
   const math_Vector& mySamples;
 
   double operator()(int theIndex) const { return mySamples(theIndex); }
 };
 
 //! Brent wrapper that adapts a Value-only function for offset root finding.
 //! @tparam Function type with Value(double theX, double& theF) method
 template <typename Function>
 struct MultipleBrentValueWrapper
 {
   Function& myFunc;
   double    myOffset;
 
   bool Value(double theX, double& theY) const
   {
     if (!myFunc.Value(theX, theY))
       return false;
     theY -= myOffset;
     return true;
   }
 };
 
 //! Evaluates original (non-offset) function value at a root point via Value interface.
 //! @tparam Function type with Value(double theX, double& theF) method
 template <typename Function>
 struct MultipleGetRootValueFn
 {
   Function& myFunc;
 
   double operator()(double theX) const
   {
     double aF = 0.0;
     myFunc.Value(theX, aF);
     return aF;
   }
 };
 
 // ============================================================================
 //  Functor adapters for Values (with derivative) interface
 // ============================================================================
 
 //! Wrapper exposing a function derivative through the Value() contract required by Brent.
 //! @tparam Function type with Values(double theX, double& theF, double& theDF) method
 template <typename Function>
 struct MultipleDerivativeValueWrapper
 {
   Function& myFunc;
 
   bool Value(double theX, double& theY) const
   {
     double aF = 0.0;
     return myFunc.Values(theX, aF, theY);
   }
 };
 
 //! Evaluate the original function value.
 //! @tparam Function type with Value(double theX, double& theF) method
 template <typename Function>
 bool EvaluateValue(Function& theFunc, double theX, double& theValue)
 {
   return theFunc.Value(theX, theValue);
 }
 
 //! Evaluate the function value shifted by the requested offset.
 //! @tparam Function type with Value(double theX, double& theF) method
 template <typename Function>
 bool EvaluateShiftedValue(Function& theFunc, double theX, double theOffset, double& theValue)
 {
   if (!theFunc.Value(theX, theValue))
   {
     return false;
   }
 
   theValue -= theOffset;
   return true;
 }
 
 //! Evaluate the function value and derivative, then shift the value by the requested offset.
 //! @tparam Function type with Values(double theX, double& theF, double& theDF) method
 template <typename Function>
 bool EvaluateShiftedValues(Function& theFunc,
                            double    theX,
                            double    theOffset,
                            double&   theValue,
                            double&   theDerivative)
 {
   if (!theFunc.Values(theX, theValue, theDerivative))
   {
     return false;
   }
 
   theValue -= theOffset;
   return true;
 }
 
 //! Refine a bracketed sign change using the same Brent/Newton sequence as math_FunctionRoots.
 //! @tparam Function type with Value() and Values() methods
 template <typename Function>
 bool RefineBracketedRoot(Function&       theFunc,
                          double          theOffset,
                          double          theX1,
                          double          theY1,
                          double          theX2,
                          double          theY2,
                          double          theTolerance,
                          double          theEpsX,
                          MultipleResult& theResult)
 {
   constexpr int    THE_MAX_ITERATIONS = 100;
   constexpr double THE_EPS2           = 2.0e-14;
   constexpr double THE_DERIV_EPS      = 1.0e-10;
 
   int    anIter = 0;
   double aTol2  = 0.5 * theTolerance;
   double aA     = theX1;
   double aB     = theX2;
   double aC     = theX2;
   double aD     = 0.0;
   double anE    = 0.0;
   double aFa    = theY1;
   double aFb    = theY2;
   double aFc    = theY2;
 
   for (anIter = 1; anIter <= THE_MAX_ITERATIONS; ++anIter)
   {
     if ((aFb > 0.0 && aFc > 0.0) || (aFb < 0.0 && aFc < 0.0))
     {
       aC  = aA;
       aFc = aFa;
       anE = aD = aB - aA;
     }
 
     if (std::abs(aFc) < std::abs(aFb))
     {
       const double aPrevA  = aA;
       const double aPrevFa = aFa;
       aA                   = aB;
       aB                   = aC;
       aC                   = aPrevA;
       aFa                  = aFb;
       aFb                  = aFc;
       aFc                  = aPrevFa;
     }
 
     const double aTol1 = THE_EPS2 * std::abs(aB) + aTol2;
     const double aXm   = 0.5 * (aC - aB);
     if (std::abs(aXm) < aTol1 || aFb == 0.0)
     {
       double aNewtonX = aB;
       for (int aNewtonIter = 0; aNewtonIter < 5; ++aNewtonIter)
       {
         double aY    = 0.0;
         double aDfdx = 0.0;
         if (!EvaluateShiftedValues(theFunc, aNewtonX, theOffset, aY, aDfdx))
         {
           return false;
         }
 
         if (std::abs(aDfdx) <= THE_DERIV_EPS)
         {
           break;
         }
 
         aNewtonX -= aY / aDfdx;
         if (aNewtonX < theX1 || aNewtonX > theX2)
         {
           break;
         }
 
         if (!EvaluateShiftedValue(theFunc, aNewtonX, theOffset, aY))
         {
           return false;
         }
 
         if (std::abs(aY) < std::abs(aFb))
         {
           aB  = aNewtonX;
           aFb = aY;
         }
       }
 
       double aRootValue = 0.0;
       if (!EvaluateValue(theFunc, aB, aRootValue))
       {
         return false;
       }
 
       theResult.NbIterations += static_cast<size_t>(anIter);
       AddRoot(theResult, theEpsX, aB, aRootValue);
       return true;
     }
 
     if (std::abs(anE) >= aTol1 && std::abs(aFa) > std::abs(aFb))
     {
       double       aP = 0.0;
       double       aQ = 0.0;
       const double aS = aFb / aFa;
       if (aA == aC)
       {
         aP = 2.0 * aXm * aS;
         aQ = 1.0 - aS;
       }
       else
       {
         aQ              = aFa / aFc;
         const double aR = aFb / aFc;
         aP              = aS * (2.0 * aXm * aQ * (aQ - aR) - (aB - aA) * (aR - 1.0));
         aQ              = (aQ - 1.0) * (aR - 1.0) * (aS - 1.0);
       }
 
       if (aP > 0.0)
       {
         aQ = -aQ;
       }
 
       aP                 = std::abs(aP);
       const double aMin1 = 3.0 * aXm * aQ - std::abs(aTol1 * aQ);
       const double aMin2 = std::abs(anE * aQ);
       if (2.0 * aP < std::min(aMin1, aMin2))
       {
         anE = aD;
         aD  = aP / aQ;
       }
       else
       {
         aD  = aXm;
         anE = aD;
       }
     }
     else
     {
       aD  = aXm;
       anE = aD;
     }
 
     aA  = aB;
     aFa = aFb;
     if (std::abs(aD) > aTol1)
     {
       aB += aD;
     }
     else
     {
       aB += (aXm >= 0.0) ? std::abs(aTol1) : -std::abs(aTol1);
     }
 
     if (!EvaluateShiftedValue(theFunc, aB, theOffset, aFb))
     {
       return false;
     }
   }
 
   theResult.NbIterations += THE_MAX_ITERATIONS;
   return true;
 }
 
 //! Derivative-aware implementation mirroring the proven math_FunctionRoots heuristics
 //! while operating on the modern callable-based MathRoot API.
 //! @tparam Function type with Value() and Values() methods
 template <typename Function>
 MultipleResult FindAllRootsWithDerivativeImpl(Function&             theFunc,
                                               double                theLower,
                                               double                theUpper,
                                               const MultipleConfig& theConfig)
 {
   MultipleResult aResult;
   aResult.Status = MathUtils::Status::OK;
 
   const double aLower         = std::min(theLower, theUpper);
   const double aUpper         = std::max(theLower, theUpper);
   const int    aNbSamples     = std::max(2 * theConfig.NbSamples, 20);
   const double aDx            = (aUpper - aLower) / aNbSamples;
   const double aEpsX          = EffectiveXTolerance(aLower, aUpper, theConfig.XTolerance);
   const double aRawXTolerance = theConfig.XTolerance;
   const double aMajorDx       = 5.0 * aDx;
 
   math_Vector aValues(0, aNbSamples);
   double      aX = aLower;
   for (int i = 0; i <= aNbSamples; ++i, aX += aDx)
   {
     if (aX > aUpper)
     {
       aX = aUpper;
     }
 
     if (!EvaluateShiftedValue(theFunc, aX, theConfig.Offset, aValues(i)))
     {
       aResult.Status = MathUtils::Status::NumericalError;
       return aResult;
     }
   }
 
   aResult.IsAllNull = true;
   for (int i = 0; i <= aNbSamples; ++i)
   {
     if (aValues(i) > theConfig.NullTolerance || aValues(i) < -theConfig.NullTolerance)
     {
       aResult.IsAllNull = false;
       break;
     }
   }
   if (aResult.IsAllNull)
   {
     return aResult;
   }
 
   const double aTolerance = aEpsX;
   double       aX1        = aLower;
   for (int i = 0, anIp1 = 1; i < aNbSamples; ++i, ++anIp1, aX1 += aDx)
   {
     double aX2 = aX1 + aDx;
     if (aX2 > aUpper)
     {
       aX2 = aUpper;
     }
 
     if ((aValues(i) < 0.0 && aValues(anIp1) > 0.0) || (aValues(i) > 0.0 && aValues(anIp1) < 0.0))
     {
       if (!RefineBracketedRoot(theFunc,
                                theConfig.Offset,
                                aX1,
                                aValues(i),
                                aX2,
                                aValues(anIp1),
                                aTolerance,
                                aEpsX,
                                aResult))
       {
         aResult.Status = MathUtils::Status::NumericalError;
         return aResult;
       }
     }
   }
 
   for (int i = 0; i <= aNbSamples; ++i)
   {
     if (aValues(i) != 0.0)
     {
       continue;
     }
 
     const double aZeroX  = std::min(aLower + i * aDx, aUpper);
     double       aLeftX  = aZeroX - 0.5 * aDx;
     double       aRightX = aZeroX + 0.5 * aDx;
     if (aLeftX < aLower)
     {
       aLeftX = aLower;
     }
     if (aLeftX > aUpper)
     {
       aLeftX = aUpper;
     }
     if (aRightX < aLower)
     {
       aRightX = aLower;
     }
     if (aRightX > aUpper)
     {
       aRightX = aUpper;
     }
 
     double aLeftY  = 0.0;
     double aRightY = 0.0;
     if (!EvaluateShiftedValue(theFunc, aLeftX, theConfig.Offset, aLeftY)
         || !EvaluateShiftedValue(theFunc, aRightX, theConfig.Offset, aRightY))
     {
       aResult.Status = MathUtils::Status::NumericalError;
       return aResult;
     }
 
     if (aLeftY * aRightY < 0.0)
     {
       if (!RefineBracketedRoot(theFunc,
                                theConfig.Offset,
                                aLeftX,
                                aLeftY,
                                aRightX,
                                aRightY,
                                aTolerance,
                                aEpsX,
                                aResult))
       {
         aResult.Status = MathUtils::Status::NumericalError;
         return aResult;
       }
     }
     else if (aLeftY != 0.0 || aRightY != 0.0)
     {
       double aRootValue = 0.0;
       if (!EvaluateValue(theFunc, aZeroX, aRootValue))
       {
         aResult.Status = MathUtils::Status::NumericalError;
         return aResult;
       }
       AddRoot(aResult, aEpsX, aZeroX, aRootValue);
     }
   }
 
   if (aValues(0) <= theConfig.FTolerance && aValues(0) >= -theConfig.FTolerance)
   {
     double aRootValue = 0.0;
     if (!EvaluateValue(theFunc, aLower, aRootValue))
     {
       aResult.Status = MathUtils::Status::NumericalError;
       return aResult;
     }
     AddRoot(aResult, aEpsX, aLower, aRootValue);
   }
 
   if (aValues(aNbSamples) <= theConfig.FTolerance && aValues(aNbSamples) >= -theConfig.FTolerance)
   {
     double aRootValue = 0.0;
     if (!EvaluateValue(theFunc, aUpper, aRootValue))
     {
       aResult.Status = MathUtils::Status::NumericalError;
       return aResult;
     }
     AddRoot(aResult, aEpsX, aUpper, aRootValue);
   }
 
   int    anIm1 = 0;
   int    anIp1 = 2;
   double aMidX = aLower + aDx;
   for (int i = 1; i < aNbSamples; ++i, ++anIm1, ++anIp1, aMidX += aDx)
   {
     if (aMidX > aUpper)
     {
       aMidX = aUpper;
     }
 
     bool isRediscretize = false;
     if (aValues(i) > 0.0)
     {
       if (aValues(anIm1) > aValues(i) && aValues(anIp1) > aValues(i))
       {
         double aProbeX  = std::max(aLower, aMidX - aDx);
         double aProbeY  = 0.0;
         double aProbeDy = 0.0;
         if (!EvaluateShiftedValues(theFunc, aProbeX, theConfig.Offset, aProbeY, aProbeDy))
         {
           aResult.Status = MathUtils::Status::NumericalError;
           return aResult;
         }
 
         if (std::abs(aProbeDy) > 1.0e-10)
         {
           const double aStep = aProbeY / aProbeDy;
           if (aStep < aMajorDx && aStep > -aMajorDx)
           {
             isRediscretize = true;
           }
         }
 
         if (!isRediscretize)
         {
           aProbeX = std::min(aUpper, aMidX + aDx);
           if (!EvaluateShiftedValues(theFunc, aProbeX, theConfig.Offset, aProbeY, aProbeDy))
           {
             aResult.Status = MathUtils::Status::NumericalError;
             return aResult;
           }
 
           if (std::abs(aProbeDy) > 1.0e-10)
           {
             const double aStep = aProbeY / aProbeDy;
             if (aStep < aMajorDx && aStep > -aMajorDx)
             {
               isRediscretize = true;
             }
           }
         }
       }
     }
     else if (aValues(i) < 0.0)
     {
       if (aValues(anIm1) < aValues(i) && aValues(anIp1) < aValues(i))
       {
         double aProbeX  = std::max(aLower, aMidX - aDx);
         double aProbeY  = 0.0;
         double aProbeDy = 0.0;
         if (!EvaluateShiftedValues(theFunc, aProbeX, theConfig.Offset, aProbeY, aProbeDy))
         {
           aResult.Status = MathUtils::Status::NumericalError;
           return aResult;
         }
 
         if (std::abs(aProbeDy) > 1.0e-10)
         {
           const double aStep = aProbeY / aProbeDy;
           if (aStep < aMajorDx && aStep > -aMajorDx)
           {
             isRediscretize = true;
           }
         }
 
         if (!isRediscretize)
         {
           aProbeX = std::min(aUpper, aMidX + aDx);
           if (!EvaluateShiftedValues(theFunc, aProbeX, theConfig.Offset, aProbeY, aProbeDy))
           {
             aResult.Status = MathUtils::Status::NumericalError;
             return aResult;
           }
 
           if (std::abs(aProbeDy) > 1.0e-10)
           {
             const double aStep = aProbeY / aProbeDy;
             if (aStep < aMajorDx && aStep > -aMajorDx)
             {
               isRediscretize = true;
             }
           }
         }
       }
     }
 
     if (!isRediscretize)
     {
       continue;
     }
 
     double aX0      = std::max(aLower, aMidX - aDx);
     double aX3      = std::min(aUpper, aMidX + aDx);
     double aRoot1   = 0.0;
     double aRoot2   = 0.0;
     double aVal1    = 0.0;
     double aVal2    = 0.0;
     double aDer1    = 0.0;
     double aDer2    = 0.0;
     bool   hasRoot1 = false;
     bool   hasRoot2 = false;
 
     MultipleDerivativeValueWrapper<Function> aDerivativeWrapper{theFunc};
     MathUtils::Config                        aDerivativeConfig;
     aDerivativeConfig.XTolerance    = aRawXTolerance;
     aDerivativeConfig.FTolerance    = 0.0;
     aDerivativeConfig.MaxIterations = theConfig.MaxIterations;
 
     MathUtils::ScalarResult aDerivativeRoot =
       Brent(aDerivativeWrapper, aX0, aX3, aDerivativeConfig);
     aResult.NbIterations += aDerivativeRoot.NbIterations;
     if (aDerivativeRoot.IsDone() && aDerivativeRoot.Root.has_value())
     {
       aRoot1                = *aDerivativeRoot.Root;
       double aOriginalValue = 0.0;
       if (!EvaluateValue(theFunc, aRoot1, aOriginalValue))
       {
         aResult.Status = MathUtils::Status::NumericalError;
         return aResult;
       }
       aVal1 = std::abs(aOriginalValue - theConfig.Offset);
       if (aVal1 < theConfig.FTolerance)
       {
         hasRoot1 = true;
         if (!EvaluateShiftedValues(theFunc, aRoot1, theConfig.Offset, aOriginalValue, aDer1))
         {
           aResult.Status = MathUtils::Status::NumericalError;
           return aResult;
         }
       }
     }
 
     double           aXProbe1             = 0.0;
     double           aXProbe2             = 0.0;
     constexpr double THE_GOLDEN_INV_RATIO = 1.0 / MathUtils::THE_GOLDEN_RATIO;
     constexpr double THE_GOLDEN_INV_COMP  = 1.0 - THE_GOLDEN_INV_RATIO;
     const double     aTolCR               = aEpsX * 10.0;
     const double     aLocalTolX           = 0.001 * aEpsX;
     double           aF0                  = aValues(anIm1);
     double           aF3                  = aValues(anIp1);
     const bool       isSearchMinimum      = (aF0 > 0.0);
 
     if (std::abs(aX3 - aMidX) > std::abs(aX0 - aMidX))
     {
       aXProbe1 = aMidX;
       aXProbe2 = aMidX + THE_GOLDEN_INV_COMP * (aX3 - aMidX);
     }
     else
     {
       aXProbe2 = aMidX;
       aXProbe1 = aMidX - THE_GOLDEN_INV_COMP * (aMidX - aX0);
     }
 
     double aF1 = 0.0;
     double aF2 = 0.0;
     if (!EvaluateShiftedValue(theFunc, aXProbe1, theConfig.Offset, aF1)
         || !EvaluateShiftedValue(theFunc, aXProbe2, theConfig.Offset, aF2))
     {
       aResult.Status = MathUtils::Status::NumericalError;
       return aResult;
     }
 
     while (std::abs(aX3 - aX0) > aTolCR * (std::abs(aXProbe1) + std::abs(aXProbe2))
            && std::abs(aXProbe1 - aXProbe2) > aLocalTolX)
     {
       if (isSearchMinimum)
       {
         if (aF2 < aF1)
         {
           aX0      = aXProbe1;
           aXProbe1 = aXProbe2;
           aXProbe2 = THE_GOLDEN_INV_RATIO * aXProbe1 + THE_GOLDEN_INV_COMP * aX3;
           aF0      = aF1;
           aF1      = aF2;
           if (!EvaluateShiftedValue(theFunc, aXProbe2, theConfig.Offset, aF2))
           {
             aResult.Status = MathUtils::Status::NumericalError;
             return aResult;
           }
         }
         else
         {
           aX3      = aXProbe2;
           aXProbe2 = aXProbe1;
           aXProbe1 = THE_GOLDEN_INV_RATIO * aXProbe2 + THE_GOLDEN_INV_COMP * aX0;
           aF3      = aF2;
           aF2      = aF1;
           if (!EvaluateShiftedValue(theFunc, aXProbe1, theConfig.Offset, aF1))
           {
             aResult.Status = MathUtils::Status::NumericalError;
             return aResult;
           }
         }
       }
       else
       {
         if (aF2 > aF1)
         {
           aX0      = aXProbe1;
           aXProbe1 = aXProbe2;
           aXProbe2 = THE_GOLDEN_INV_RATIO * aXProbe1 + THE_GOLDEN_INV_COMP * aX3;
           aF0      = aF1;
           aF1      = aF2;
           if (!EvaluateShiftedValue(theFunc, aXProbe2, theConfig.Offset, aF2))
           {
             aResult.Status = MathUtils::Status::NumericalError;
             return aResult;
           }
         }
         else
         {
           aX3      = aXProbe2;
           aXProbe2 = aXProbe1;
           aXProbe1 = THE_GOLDEN_INV_RATIO * aXProbe2 + THE_GOLDEN_INV_COMP * aX0;
           aF3      = aF2;
           aF2      = aF1;
           if (!EvaluateShiftedValue(theFunc, aXProbe1, theConfig.Offset, aF1))
           {
             aResult.Status = MathUtils::Status::NumericalError;
             return aResult;
           }
         }
       }
 
       if (aF1 * aF0 < 0.0)
       {
         if (!RefineBracketedRoot(theFunc,
                                  theConfig.Offset,
                                  aX0,
                                  aF0,
                                  aXProbe1,
                                  aF1,
                                  aTolerance,
                                  aEpsX,
                                  aResult))
         {
           aResult.Status = MathUtils::Status::NumericalError;
           return aResult;
         }
       }
 
       if (aF2 * aF3 < 0.0)
       {
         if (!RefineBracketedRoot(theFunc,
                                  theConfig.Offset,
                                  aXProbe2,
                                  aF2,
                                  aX3,
                                  aF3,
                                  aTolerance,
                                  aEpsX,
                                  aResult))
         {
           aResult.Status = MathUtils::Status::NumericalError;
           return aResult;
         }
       }
     }
 
     if ((isSearchMinimum && aF1 < aF2) || (!isSearchMinimum && aF1 > aF2))
     {
       if (std::abs(aF1) < theConfig.FTolerance)
       {
         hasRoot2 = true;
         aRoot2   = aXProbe1;
         aVal2    = std::abs(aF1);
       }
     }
     else if (std::abs(aF2) < theConfig.FTolerance)
     {
       hasRoot2 = true;
       aRoot2   = aXProbe2;
       aVal2    = std::abs(aF2);
     }
 
     if (hasRoot1 && hasRoot2)
     {
       if (aVal2 - aVal1 > theConfig.FTolerance)
       {
         double aRootValue = 0.0;
         if (!EvaluateValue(theFunc, aRoot1, aRootValue))
         {
           aResult.Status = MathUtils::Status::NumericalError;
           return aResult;
         }
         AddRoot(aResult, aEpsX, aRoot1, aRootValue);
       }
       else if (aVal1 - aVal2 > theConfig.FTolerance)
       {
         double aRootValue = 0.0;
         if (!EvaluateValue(theFunc, aRoot2, aRootValue))
         {
           aResult.Status = MathUtils::Status::NumericalError;
           return aResult;
         }
         AddRoot(aResult, aEpsX, aRoot2, aRootValue);
       }
       else
       {
         double aShiftedValue = 0.0;
         if (!EvaluateShiftedValues(theFunc, aRoot2, theConfig.Offset, aShiftedValue, aDer2))
         {
           aResult.Status = MathUtils::Status::NumericalError;
           return aResult;
         }
 
         const double aChosenRoot = (std::abs(aDer1) < std::abs(aDer2)) ? aRoot1 : aRoot2;
         double       aRootValue  = 0.0;
         if (!EvaluateValue(theFunc, aChosenRoot, aRootValue))
         {
           aResult.Status = MathUtils::Status::NumericalError;
           return aResult;
         }
         AddRoot(aResult, aEpsX, aChosenRoot, aRootValue);
       }
     }
     else if (hasRoot1)
     {
       double aRootValue = 0.0;
       if (!EvaluateValue(theFunc, aRoot1, aRootValue))
       {
         aResult.Status = MathUtils::Status::NumericalError;
         return aResult;
       }
       AddRoot(aResult, aEpsX, aRoot1, aRootValue);
     }
     else if (hasRoot2)
     {
       double aRootValue = 0.0;
       if (!EvaluateValue(theFunc, aRoot2, aRootValue))
       {
         aResult.Status = MathUtils::Status::NumericalError;
         return aResult;
       }
       AddRoot(aResult, aEpsX, aRoot2, aRootValue);
     }
   }
 
   SortRoots(aResult);
   return aResult;
 }
 
 // ============================================================================
 //  Interval handlers
 // ============================================================================
 
 //! No-op interval handler for functions without derivative.
 struct MultipleNoExtraHandler
 {
   void operator()(int, double, double, double, double, MultipleResult&, double) const {}
 };
 
 // ============================================================================
 //  Core implementation
 // ============================================================================
 
 //! Core implementation for finding all roots in an interval.
 //! Shared logic for both Value-only and Values (with derivative) interfaces.
 //!
 //! @tparam SampleFn callable (int theIndex, double theX) -> bool
 //! @tparam GetValueFn callable (int theIndex) -> double, returns f(x)-offset at sample
 //! @tparam BrentWrapperT type with Value(double, double&) for Brent root finding
 //! @tparam GetRootValueFn callable (double theX) -> double, returns original f(x)
 //! @tparam IntervalExtraFn callable (int, x0, x1, f0, f1, result, epsX) -> void
 template <typename SampleFn,
           typename GetValueFn,
           typename BrentWrapperT,
           typename GetRootValueFn,
           typename IntervalExtraFn>
 MultipleResult FindAllRootsImpl(double                theLower,
                                 double                theUpper,
                                 const MultipleConfig& theConfig,
                                 SampleFn              theSampleFn,
                                 GetValueFn            theGetValue,
                                 BrentWrapperT&        theBrentWrapper,
                                 GetRootValueFn        theGetRootValue,
                                 IntervalExtraFn       theIntervalExtra)
 {
   MultipleResult aResult;
   aResult.Status = MathUtils::Status::OK;
 
   // Ensure proper ordering
   const double aLower = std::min(theLower, theUpper);
   const double aUpper = std::max(theLower, theUpper);
 
   // Minimum samples
   const int    aNbSamples = std::max(2 * theConfig.NbSamples, 20);
   const double aDx        = (aUpper - aLower) / aNbSamples;
 
   // Ensure EpsX is not too small relative to interval
   const double aEpsX = EffectiveXTolerance(aLower, aUpper, theConfig.XTolerance);
 
   // Sample function values
   math_Vector aXValues(0, aNbSamples);
   for (int i = 0; i <= aNbSamples; ++i)
   {
     double aX = aLower + i * aDx;
     if (aX > aUpper)
       aX = aUpper;
     aXValues(i) = aX;
 
     if (!theSampleFn(i, aX))
     {
       aResult.Status = MathUtils::Status::NumericalError;
       return aResult;
     }
   }
 
   // Check if function is essentially null everywhere
   aResult.IsAllNull = true;
   for (int i = 0; i <= aNbSamples; ++i)
   {
     if (std::abs(theGetValue(i)) > theConfig.NullTolerance)
     {
       aResult.IsAllNull = false;
       break;
     }
   }
 
   if (aResult.IsAllNull)
   {
     return aResult;
   }
 
   // Find sign changes
   for (int i = 0; i < aNbSamples; ++i)
   {
     const double aF0 = theGetValue(i);
     const double aF1 = theGetValue(i + 1);
     const double aX0 = aXValues(i);
     const double aX1 = aXValues(i + 1);
 
     // Exact zero at sample point
     if (std::abs(aF0) < theConfig.FTolerance)
     {
       AddRoot(aResult, aEpsX, aX0, aF0 + theConfig.Offset);
       continue;
     }
 
     // Sign change detected
     if (aF0 * aF1 < 0.0)
     {
       MathUtils::Config aBrentConfig;
       aBrentConfig.XTolerance    = aEpsX;
       aBrentConfig.FTolerance    = theConfig.FTolerance;
       aBrentConfig.MaxIterations = theConfig.MaxIterations;
 
       MathUtils::ScalarResult aBrentResult = Brent(theBrentWrapper, aX0, aX1, aBrentConfig);
       aResult.NbIterations += aBrentResult.NbIterations;
 
       if (aBrentResult.IsDone() && aBrentResult.Root.has_value())
       {
         AddRoot(aResult, aEpsX, *aBrentResult.Root, theGetRootValue(*aBrentResult.Root));
       }
     }
     else
     {
       // Additional per-interval processing (e.g., tangential root detection)
       theIntervalExtra(i, aX0, aX1, aF0, aF1, aResult, aEpsX);
     }
   }
 
   // Check last sample point
   if (std::abs(theGetValue(aNbSamples)) < theConfig.FTolerance)
   {
     AddRoot(aResult, aEpsX, aXValues(aNbSamples), theGetValue(aNbSamples) + theConfig.Offset);
   }
 
   SortRoots(aResult);
   return aResult;
 }
 
 } // namespace MathRoot
 
 #endif // _MathRoot_MultipleUtils_HeaderFile

This page was automatically generated by the 2.3.7 LXR engine. The LXR team