include/opencascade/MathPoly_Quadratic.hxx

0001 // Copyright (c) 2025 OPEN CASCADE SAS
0002 //
0003 // This file is part of Open CASCADE Technology software library.
0004 //
0005 // This library is free software; you can redistribute it and/or modify it under
0006 // the terms of the GNU Lesser General Public License version 2.1 as published
0007 // by the Free Software Foundation, with special exception defined in the file
0008 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0009 // distribution for complete text of the license and disclaimer of any warranty.
0010 //
0011 // Alternatively, this file may be used under the terms of Open CASCADE
0012 // commercial license or contractual agreement.
0013
0014 #ifndef _MathPoly_Quadratic_HeaderFile
0015 #define _MathPoly_Quadratic_HeaderFile
0016
0017 #include <MathUtils_Types.hxx>
0018 #include <MathUtils_Core.hxx>
0019
0020 #include <cmath>
0021
0022 //! Polynomial root finding algorithms.
0023 namespace MathPoly
0024 {
0025 using namespace MathUtils;
0026
0027 //! Solve linear equation: a*x + b = 0
0028 //! Handles degenerate cases (a = 0).
0029 //! @param theA coefficient of x
0030 //! @param theB constant term
0031 //! @return result containing 0 or 1 root, or infinite solutions flag
0032 #ifdef _MSC_VER
0033   #pragma warning(push)
0034   #pragma warning(disable : 4723) // potential divide by 0 - guarded by IsZero() check
0035 #endif
0036 inline MathUtils::PolyResult Linear(double theA, double theB)
0037 {
0038   MathUtils::PolyResult aResult;
0039
0040   if (MathUtils::IsZero(theA))
0041   {
0042     if (MathUtils::IsZero(theB))
0043     {
0044       aResult.Status = MathUtils::Status::InfiniteSolutions;
0045     }
0046     else
0047     {
0048       aResult.Status = MathUtils::Status::NoSolution;
0049     }
0050   }
0051   else
0052   {
0053     aResult.Status   = MathUtils::Status::OK;
0054     aResult.NbRoots  = 1;
0055     aResult.Roots[0] = -theB / theA;
0056   }
0057
0058   return aResult;
0059 }
0060 #ifdef _MSC_VER
0061   #pragma warning(pop)
0062 #endif
0063
0064 //! Solve quadratic equation: a*x^2 + b*x + c = 0
0065 //! Uses numerically stable formulas to avoid catastrophic cancellation.
0066 //!
0067 //! Algorithm:
0068 //! 1. Handle linear case when a = 0
0069 //! 2. Compute discriminant D = b^2 - 4ac
0070 //! 3. For D < 0: no real roots
0071 //! 4. For D = 0: one double root
0072 //! 5. For D > 0: two roots using stable formula
0073 //!
0074 //! @param theA coefficient of x^2
0075 //! @param theB coefficient of x
0076 //! @param theC constant term
0077 //! @return result containing 0, 1, or 2 real roots (sorted in ascending order)
0078 inline MathUtils::PolyResult Quadratic(double theA, double theB, double theC)
0079 {
0080   MathUtils::PolyResult aResult;
0081
0082   // Linear case: b*x + c = 0
0083   if (MathUtils::IsZero(theA))
0084   {
0085     return Linear(theB, theC);
0086   }
0087
0088   // Scale coefficients for numerical stability
0089   const double aScale = std::max({std::abs(theA), std::abs(theB), std::abs(theC)});
0090   if (aScale < MathUtils::THE_ZERO_TOL)
0091   {
0092     aResult.Status = MathUtils::Status::InfiniteSolutions;
0093     return aResult;
0094   }
0095
0096   const double aA = theA / aScale;
0097   const double aB = theB / aScale;
0098   const double aC = theC / aScale;
0099
0100   // Compute discriminant
0101   const double aDisc = aB * aB - 4.0 * aA * aC;
0102
0103   // Tolerance for discriminant based on coefficient magnitudes
0104   const double aDiscTol = MathUtils::THE_ZERO_TOL * (aB * aB + std::abs(4.0 * aA * aC));
0105
0106   if (aDisc < -aDiscTol)
0107   {
0108     // No real roots (complex conjugate pair)
0109     aResult.Status  = MathUtils::Status::OK;
0110     aResult.NbRoots = 0;
0111     return aResult;
0112   }
0113
0114   if (std::abs(aDisc) <= aDiscTol)
0115   {
0116     // Double root
0117     aResult.Status   = MathUtils::Status::OK;
0118     aResult.NbRoots  = 1;
0119     aResult.Roots[0] = -aB / (2.0 * aA);
0120     return aResult;
0121   }
0122
0123   // Two distinct real roots
0124   // Use numerically stable formula to avoid catastrophic cancellation
0125   const double aSqrtDisc = std::sqrt(aDisc);
0126
0127   // q = -0.5 * (b + sign(b) * sqrt(discriminant))
0128   const double aQ = -0.5 * (aB + MathUtils::SignTransfer(aSqrtDisc, aB));
0129
0130   aResult.Status   = MathUtils::Status::OK;
0131   aResult.NbRoots  = 2;
0132   aResult.Roots[0] = aQ / aA;
0133   aResult.Roots[1] = aC / aQ;
0134
0135   // Sort roots in ascending order
0136   if (aResult.Roots[0] > aResult.Roots[1])
0137   {
0138     std::swap(aResult.Roots[0], aResult.Roots[1]);
0139   }
0140
0141   return aResult;
0142 }
0143
0144 } // namespace MathPoly
0145
0146 #endif // _MathPoly_Quadratic_HeaderFile