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

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

 // 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 _MathPoly_Cubic_HeaderFile
 #define _MathPoly_Cubic_HeaderFile
 
 #include <MathUtils_Types.hxx>
 #include <MathUtils_Core.hxx>
 #include <MathUtils_Poly.hxx>
 #include <MathPoly_Quadratic.hxx>
 
 #include <cmath>
 
 //! Polynomial root finding algorithms.
 namespace MathPoly
 {
 using namespace MathUtils;
 
 //! Solve cubic equation: a*x^3 + b*x^2 + c*x + d = 0
 //! Uses Cardano's method with Vieta substitution and trigonometric solution.
 //!
 //! Algorithm:
 //! 1. Handle lower degree cases (a = 0 -> quadratic)
 //! 2. Transform to depressed cubic: t^3 + pt + q = 0 via x = t - b/(3a)
 //! 3. Compute discriminant Delta = q^2/4 + p^3/27
 //! 4. For Delta > 0: one real root (Cardano's formula)
 //! 5. For Delta < 0: three real roots (trigonometric method)
 //! 6. For Delta = 0: multiple roots
 //! 7. Apply Newton-Raphson refinement
 //!
 //! @param theA coefficient of x^3
 //! @param theB coefficient of x^2
 //! @param theC coefficient of x
 //! @param theD constant term
 //! @return result containing 1, 2, or 3 real roots (sorted in ascending order)
 inline MathUtils::PolyResult Cubic(double theA, double theB, double theC, double theD)
 {
   MathUtils::PolyResult aResult;
 
   // Reduce to quadratic if leading coefficient is zero
   if (MathUtils::IsZero(theA))
   {
     return Quadratic(theB, theC, theD);
   }
 
   // Scale coefficients for numerical stability
   const double aScale = std::max({std::abs(theA), std::abs(theB), std::abs(theC), std::abs(theD)});
   if (aScale < MathUtils::THE_ZERO_TOL)
   {
     aResult.Status = MathUtils::Status::InfiniteSolutions;
     return aResult;
   }
 
   // Normalize to monic cubic: x^3 + px^2 + qx + r = 0
   const double aP = theB / theA;
   const double aQ = theC / theA;
   const double aR = theD / theA;
 
   // Substitute x = t - p/3 to get depressed cubic: t^3 + at + b = 0
   const double aP3    = aP / 3.0;
   const double aP3_sq = aP3 * aP3;
   const double a      = aQ - 3.0 * aP3_sq;
   const double b      = aR - aP3 * aQ + 2.0 * aP3_sq * aP3;
 
   // Discriminant: Delta = (b/2)^2 + (a/3)^3
   const double aHalfB  = b / 2.0;
   const double aThirdA = a / 3.0;
   const double aDisc   = aHalfB * aHalfB + aThirdA * aThirdA * aThirdA;
 
   // Tolerance for discriminant
   const double aDiscTol =
     MathUtils::THE_ZERO_TOL * std::max(aHalfB * aHalfB, std::abs(aThirdA * aThirdA * aThirdA));
 
   // Store original coefficients for refinement
   const double aCoeffs[4] = {theD, theC, theB, theA};
 
   if (aDisc > aDiscTol)
   {
     // One real root, two complex conjugate roots
     // Cardano's formula: t = cbrt(-b/2 + sqrt(disc)) + cbrt(-b/2 - sqrt(disc))
     const double aSqrtDisc = std::sqrt(aDisc);
     const double aU        = MathUtils::CubeRoot(-aHalfB + aSqrtDisc);
     const double aV        = MathUtils::CubeRoot(-aHalfB - aSqrtDisc);
 
     aResult.Status   = MathUtils::Status::OK;
     aResult.NbRoots  = 1;
     aResult.Roots[0] = aU + aV - aP3;
 
     // Refine root
     aResult.Roots[0] = MathUtils::RefinePolyRoot(aCoeffs, 3, aResult.Roots[0]);
   }
   else if (aDisc < -aDiscTol)
   {
     // Three distinct real roots - use trigonometric method (casus irreducibilis)
     // t = 2*sqrt(-a/3) * cos(theta/3 + 2*k*pi/3), k = 0, 1, 2
     // where cos(theta) = -b / (2 * sqrt((-a/3)^3))
 
     const double aR_val        = std::sqrt(-aThirdA * aThirdA * aThirdA);
     const double aCosArg       = MathUtils::Clamp(-aHalfB / aR_val, -1.0, 1.0);
     const double aTheta        = std::acos(aCosArg);
     const double aTwoSqrtNegA3 = 2.0 * std::sqrt(-aThirdA);
 
     aResult.Status   = MathUtils::Status::OK;
     aResult.NbRoots  = 3;
     aResult.Roots[0] = aTwoSqrtNegA3 * std::cos(aTheta / 3.0) - aP3;
     aResult.Roots[1] = aTwoSqrtNegA3 * std::cos((aTheta + 2.0 * MathUtils::THE_PI) / 3.0) - aP3;
     aResult.Roots[2] = aTwoSqrtNegA3 * std::cos((aTheta + 4.0 * MathUtils::THE_PI) / 3.0) - aP3;
 
     // Refine all roots
     for (int i = 0; i < 3; ++i)
     {
       aResult.Roots[i] = MathUtils::RefinePolyRoot(aCoeffs, 3, aResult.Roots[i]);
     }
 
     // Sort roots
     MathUtils::SortRoots(aResult.Roots.data(), 3);
   }
   else
   {
     // Discriminant is zero - multiple roots
     const double aU = MathUtils::CubeRoot(-aHalfB);
 
     aResult.Status = MathUtils::Status::OK;
 
     if (MathUtils::IsZero(aU))
     {
       // Triple root at x = -p/3
       aResult.NbRoots  = 1;
       aResult.Roots[0] = -aP3;
     }
     else
     {
       // One single root and one double root
       aResult.NbRoots  = 2;
       aResult.Roots[0] = 2.0 * aU - aP3; // Single root
       aResult.Roots[1] = -aU - aP3;      // Double root
 
       // Refine roots
       aResult.Roots[0] = MathUtils::RefinePolyRoot(aCoeffs, 3, aResult.Roots[0]);
       aResult.Roots[1] = MathUtils::RefinePolyRoot(aCoeffs, 3, aResult.Roots[1]);
 
       // Sort roots
       if (aResult.Roots[0] > aResult.Roots[1])
       {
         std::swap(aResult.Roots[0], aResult.Roots[1]);
       }
     }
   }
 
   return aResult;
 }
 
 } // namespace MathPoly
 
 #endif // _MathPoly_Cubic_HeaderFile

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