EIC code displayed by LXR master/​include/​opencascade/​MathRoot_Brent.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_Brent_HeaderFile
 #define _MathRoot_Brent_HeaderFile
 
 #include <MathUtils_Types.hxx>
 #include <MathUtils_Config.hxx>
 #include <MathUtils_Core.hxx>
 
 #include <cmath>
 #include <utility>
 
 namespace MathRoot
 {
 using namespace MathUtils;
 
 //! Brent's method for root finding.
 //! Combines bisection, secant, and inverse quadratic interpolation.
 //! Guaranteed to converge if a valid bracket is provided.
 //!
 //! Algorithm:
 //! 1. Start with bracket [a, b] where f(a) * f(b) < 0
 //! 2. At each step, try inverse quadratic interpolation
 //! 3. If interpolation step is rejected, use bisection
 //! 4. Acceptance criteria ensure superlinear convergence when possible
 //!
 //! @tparam Function type with Value(double theX, double& theF) method
 //! @param theFunc function to find root of
 //! @param theLower lower bound of bracket (f(theLower) and f(theUpper) must have opposite signs)
 //! @param theUpper upper bound of bracket
 //! @param theConfig solver configuration
 //! @return result containing root location and convergence status
 template <typename Function>
 MathUtils::ScalarResult Brent(Function&                theFunc,
                               double                   theLower,
                               double                   theUpper,
                               const MathUtils::Config& theConfig = MathUtils::Config())
 {
   MathUtils::ScalarResult aResult;
 
   double aA  = theLower;
   double aB  = theUpper;
   double aFa = 0.0;
   double aFb = 0.0;
 
   // Evaluate at endpoints
   if (!theFunc.Value(aA, aFa))
   {
     aResult.Status = MathUtils::Status::NumericalError;
     return aResult;
   }
   if (!theFunc.Value(aB, aFb))
   {
     aResult.Status = MathUtils::Status::NumericalError;
     return aResult;
   }
 
   // Check that bracket is valid (sign change)
   if (aFa * aFb > 0.0)
   {
     aResult.Status = MathUtils::Status::InvalidInput;
     return aResult;
   }
 
   // Ensure |f(a)| >= |f(b)| (b is the better approximation)
   if (std::abs(aFa) < std::abs(aFb))
   {
     std::swap(aA, aB);
     std::swap(aFa, aFb);
   }
 
   double aC  = aA; // Previous iterate
   double aFc = aFa;
   double aD  = aB - aA; // Step size
   double aE  = aD;      // Previous step size
 
   for (int anIter = 0; anIter < theConfig.MaxIterations; ++anIter)
   {
     aResult.NbIterations = anIter + 1;
 
     const double aTol = 2.0 * MathUtils::THE_EPSILON * std::abs(aB) + 0.5 * theConfig.XTolerance;
     const double aM   = 0.5 * (aC - aB);
 
     if (std::abs(aFb) < theConfig.FTolerance || aFb == 0.0 || std::abs(aM) <= aTol)
     {
       aResult.Status = MathUtils::Status::OK;
       aResult.Root   = aB;
       aResult.Value  = aFb;
       return aResult;
     }
 
     double aS = 0.0; // New approximation
 
     // Try inverse quadratic interpolation if we have three distinct points
     if (std::abs(aFa - aFc) > MathUtils::THE_ZERO_TOL
         && std::abs(aFb - aFc) > MathUtils::THE_ZERO_TOL)
     {
       // Inverse quadratic interpolation
       aS = aA * aFb * aFc / ((aFa - aFb) * (aFa - aFc))
            + aB * aFa * aFc / ((aFb - aFa) * (aFb - aFc))
            + aC * aFa * aFb / ((aFc - aFa) * (aFc - aFb));
     }
     else
     {
       // Secant method
       aS = aB - aFb * (aB - aA) / (aFb - aFa);
     }
 
     // Decide whether to accept the interpolation step
     bool aUseInterp = false;
 
     // Check if s is between (3a+b)/4 and b
     const double aBound1 = (3.0 * aA + aB) / 4.0;
     if ((aS > std::min(aBound1, aB) && aS < std::max(aBound1, aB)))
     {
       // Accept interpolation if step is smaller than half the previous step
       // (ensures convergence rate). Minimum step is enforced later.
       if (std::abs(aS - aB) < std::abs(aE) / 2.0)
       {
         aUseInterp = true;
       }
     }
 
     if (!aUseInterp)
     {
       // Bisection step
       aS = aB + aM;
       aE = aM;
       aD = aM;
     }
     else
     {
       aE = aD;
       aD = aS - aB;
     }
 
     // Update previous values
     aA  = aB;
     aFa = aFb;
 
     // Compute new point, ensuring minimum step
     if (std::abs(aD) > aTol)
     {
       aB = aS;
     }
     else
     {
       aB += (aM > 0.0) ? aTol : -aTol;
     }
 
     // Evaluate function at new point
     if (!theFunc.Value(aB, aFb))
     {
       aResult.Status = MathUtils::Status::NumericalError;
       aResult.Root   = aB;
       return aResult;
     }
 
     // Update bracket
     if (aFb * aFc > 0.0)
     {
       aC  = aA;
       aFc = aFa;
       aD  = aB - aA;
       aE  = aD;
     }
     else if (std::abs(aFc) < std::abs(aFb))
     {
       // Swap b and c if c is better (use std::swap to avoid overwriting)
       aA  = aB;
       aFa = aFb;
       std::swap(aB, aC);
       std::swap(aFb, aFc);
     }
   }
 
   // Maximum iterations reached
   aResult.Status = MathUtils::Status::MaxIterations;
   aResult.Root   = aB;
   aResult.Value  = aFb;
   return aResult;
 }
 
 } // namespace MathRoot
 
 #endif // _MathRoot_Brent_HeaderFile

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