OCCT

// 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_Secant_HeaderFile
#define _MathRoot_Secant_HeaderFile

#include <MathUtils_Types.hxx>
#include <MathUtils_Config.hxx>
#include <MathUtils_Core.hxx>
#include <MathUtils_Convergence.hxx>
#include <MathRoot_Utils.hxx>

#include <cmath>

namespace MathRoot
{
using namespace MathUtils;

//! Secant method for root finding.
//! Does not require derivative, uses finite difference approximation.
//! Converges superlinearly (order ~1.618, the golden ratio).
//!
//! Algorithm:
//! 1. Start with two initial points x0, x1
//! 2. Approximate derivative: f'(x) ~ (f(x1) - f(x0)) / (x1 - x0)
//! 3. Newton-like update: x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
//! 4. Repeat with x0 = x1, x1 = x2
//!
//! @tparam Function type with Value(double theX, double& theF) method
//! @param theFunc function object providing only value
//! @param theX0 first initial point
//! @param theX1 second initial point (different from theX0)
//! @param theConfig solver configuration
//! @return result containing root location and convergence status
template <typename Function>
MathUtils::ScalarResult Secant(Function&                theFunc,
                               double                   theX0,
                               double                   theX1,
                               const MathUtils::Config& theConfig = MathUtils::Config())
{
  MathUtils::ScalarResult aResult;

  if (!std::isfinite(theX0) || !std::isfinite(theX1) || theX0 == theX1
      || !Utils::IsValidConfig(theConfig))
  {
    aResult.Status = MathUtils::Status::InvalidInput;
    return aResult;
  }

  double aX0 = theX0;
  double aX1 = theX1;
  double aF0 = 0.0;
  double aF1 = 0.0;

  // Evaluate at initial points
  if (!theFunc.Value(aX0, aF0))
  {
    aResult.Status = MathUtils::Status::CallbackError;
    return aResult;
  }
  if (!std::isfinite(aF0))
  {
    aResult.Status = MathUtils::Status::NumericalError;
    return aResult;
  }
  if (!theFunc.Value(aX1, aF1))
  {
    aResult.Status = MathUtils::Status::CallbackError;
    return aResult;
  }
  if (!std::isfinite(aF1))
  {
    aResult.Status = MathUtils::Status::NumericalError;
    return aResult;
  }

  if (std::abs(aF0) <= theConfig.FTolerance)
  {
    aResult.Status = MathUtils::Status::OK;
    aResult.Root   = aX0;
    aResult.Value  = aF0;
    return aResult;
  }

  for (uint32_t anIter = 0; anIter < theConfig.MaxIterations; ++anIter)
  {
    aResult.NbIterations = anIter + 1;

    // Check convergence
    if (std::abs(aF1) <= theConfig.FTolerance)
    {
      aResult.Status = MathUtils::Status::OK;
      aResult.Root   = aX1;
      aResult.Value  = aF1;
      return aResult;
    }

    // Secant step
    const double aDenom = aF1 - aF0;
    if (!std::isfinite(aDenom) || aDenom == 0.0)
    {
      aResult.Status = MathUtils::Status::NumericalError;
      aResult.Root   = aX1;
      aResult.Value  = aF1;
      return aResult;
    }

    const double aXNew = aX1 - (aF1 / aDenom) * (aX1 - aX0);
    if (!std::isfinite(aXNew))
    {
      aResult.Status = MathUtils::Status::NumericalError;
      aResult.Root   = aX1;
      aResult.Value  = aF1;
      return aResult;
    }

    if (MathUtils::IsXConverged(aX1, aXNew, theConfig.XTolerance))
    {
      double aFNew = 0.0;
      if (!theFunc.Value(aXNew, aFNew))
      {
        aResult.Status = MathUtils::Status::CallbackError;
        aResult.Root   = aX1;
        aResult.Value  = aF1;
        return aResult;
      }
      if (!std::isfinite(aFNew))
      {
        aResult.Status = MathUtils::Status::NumericalError;
        aResult.Root   = aX1;
        aResult.Value  = aF1;
        return aResult;
      }
      aResult.Status = (std::abs(aFNew) <= theConfig.FTolerance) ? MathUtils::Status::OK
                                                                 : MathUtils::Status::NotConverged;
      aResult.Root   = aXNew;
      aResult.Value  = aFNew;
      return aResult;
    }

    // Update for next iteration
    aX0 = aX1;
    aF0 = aF1;
    aX1 = aXNew;

    if (!theFunc.Value(aX1, aF1))
    {
      aResult.Status = MathUtils::Status::CallbackError;
      aResult.Root   = aX0;
      aResult.Value  = aF0;
      return aResult;
    }
    if (!std::isfinite(aF1))
    {
      aResult.Status = MathUtils::Status::NumericalError;
      aResult.Root   = aX0;
      aResult.Value  = aF0;
      return aResult;
    }
  }

  aResult.Status = MathUtils::Status::MaxIterations;
  aResult.Root   = aX1;
  aResult.Value  = aF1;
  return aResult;
}

//! Secant method with automatic initial points.
//! Creates second point by small perturbation of the initial guess.
//!
//! @tparam Function type with Value(double theX, double& theF) method
//! @param theFunc function object providing only value
//! @param theX0 initial guess
//! @param theConfig solver configuration
//! @return result containing root location and convergence status
template <typename Function>
MathUtils::ScalarResult SecantAuto(Function&                theFunc,
                                   double                   theX0,
                                   const MathUtils::Config& theConfig = MathUtils::Config())
{
  if (!std::isfinite(theX0))
  {
    MathUtils::ScalarResult aResult;
    aResult.Status = MathUtils::Status::InvalidInput;
    return aResult;
  }
  // Create second point by small perturbation
  const double aDelta = (std::abs(theX0) > 1.0) ? 0.01 * theX0 : 0.01;
  if (!std::isfinite(theX0 + aDelta) || theX0 + aDelta == theX0)
  {
    MathUtils::ScalarResult aResult;
    aResult.Status = MathUtils::Status::InvalidInput;
    return aResult;
  }
  return Secant(theFunc, theX0, theX0 + aDelta, theConfig);
}

} // namespace MathRoot

#endif // _MathRoot_Secant_HeaderFile