Skip to content

Const polynomial roots #20539

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 3 commits into from
+28 −32
Closed

Const polynomial roots #20539

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ea56fb6
[MathMore] allow using several methods of Polynomial.h in const context
meiyasan Nov 26, 2025
5dfb494
Merge branch 'root-project:master' into polynomial-with-const
meiyasan Nov 26, 2025
3eec0d2
Merge branch 'root-project:master' into polynomial-with-const
meiyasan Nov 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 3 additions & 8 deletions math/mathmore/inc/Math/Polynomial.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,20 +116,20 @@ class Polynomial : public ParamFunction<IParamGradFunction>,
equation is very small. In that case it might be more robust to use the numerical method, by calling directly FindNumRoots()

*/
const std::vector<std::complex <double> > & FindRoots();
const std::vector<std::complex<double>> FindRoots() const;

/**
Find the only the real polynomial roots.
For n <= 4, the roots are found analytically while for larger order an iterative numerical method is used
The numerical method used is from GSL (see <A HREF="https://www.gnu.org/software/gsl/doc/html/poly.html">documentation</A> )
*/
std::vector<double > FindRealRoots();
std::vector<double> FindRealRoots() const;

/**
Find the polynomial roots using always an iterative numerical methods
The numerical method used is from GSL (see <A HREF="https://www.gnu.org/software/gsl/doc/html/poly.html">documentation</A> )
*/
const std::vector<std::complex <double> > & FindNumRoots();
const std::vector<std::complex<double>> FindNumRoots() const;

/**
Order of Polynomial
Expand Down Expand Up @@ -164,11 +164,6 @@ class Polynomial : public ParamFunction<IParamGradFunction>,

// cache Parameters for Gradient
mutable std::vector<double> fDerived_params;

// roots

std::vector< std::complex < double > > fRoots;

};

} // namespace Math
Expand Down
49 changes: 25 additions & 24 deletions math/mathmore/src/Polynomial.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,16 +147,16 @@ IGenFunction * Polynomial::Clone() const {
return f;
}

const std::vector<std::complex<double>> Polynomial::FindRoots() const
{

const std::vector< std::complex <double> > & Polynomial::FindRoots(){


// check if order is correct
unsigned int n = fOrder;
while ( Parameters()[n] == 0 ) {
// check if order is correct
unsigned int n = fOrder;
while (Parameters()[n] == 0) {
n--;
}

std::vector<std::complex<double>> fRoots;
fRoots.clear();
fRoots.reserve(n);

Expand Down Expand Up @@ -226,32 +226,34 @@ const std::vector< std::complex <double> > & Polynomial::FindRoots(){
}
else {
// for higher order polynomial use numerical fRoots
FindNumRoots();
return FindNumRoots();
}

return fRoots;
}

}

std::vector<double> Polynomial::FindRealRoots() const
{

std::vector< double > Polynomial::FindRealRoots(){
FindRoots();
std::vector<double> roots;
roots.reserve(fOrder);
for (unsigned int i = 0; i < fOrder; ++i) {
if (fRoots[i].imag() == 0)
roots.push_back( fRoots[i].real() );
}
return roots;
std::vector<std::complex<double>> fRoots = FindRoots();
std::vector<double> roots;
roots.reserve(fOrder);
for (unsigned int i = 0; i < fOrder; ++i) {
if (fRoots[i].imag() == 0)
roots.push_back(fRoots[i].real());
}
Comment on lines +238 to +244
Copy link
Member

@vepadulano vepadulano Nov 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nitpick, consider the following perhaps a bit more idiomatic code

auto roots = FindRoots();
std::vector<double> realRoots;
realRoots.reserve(roots.size());
for (const auto &r: roots)
   if (r.imag() == 0)
      realRoots.push_back(r.real());

Copy link
Member

@vepadulano vepadulano Nov 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Now that I look at this twice, I see the for loop is stopping after fOrder iterations. Although it's not clear from the code whether the vector returned by FindRoots will have a size greater (or even smaller which would be UB!) than fOrder

Copy link
Author

@meiyasan meiyasan Nov 27, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

well.. I agree the historical implementation sounds risky; what do you think about comparing double and strict 0 integers, shouldn't we implement an epsilon limit ?
additionally, considering we are removing std::vector::reverse() in other cases, we could also drop realRoots.reserve(roots.size()); here.

here would be my preference for simplicity purpose.

auto roots = FindRoots();
std::vector<double> realRoots;
for (const auto &r: roots)
   if (std::abs(r.imag()) < std::numeric_limits<double>::epsilon())
      realRoots.push_back(r.real());

return roots;
}
const std::vector< std::complex <double> > & Polynomial::FindNumRoots(){

const std::vector<std::complex<double>> Polynomial::FindNumRoots() const
{

// check if order is correct
unsigned int n = fOrder;
while ( Parameters()[n] == 0 ) {
// check if order is correct
unsigned int n = fOrder;
while (Parameters()[n] == 0) {
n--;
}

std::vector<std::complex<double>> fRoots;
fRoots.clear();
fRoots.reserve(n);
Comment on lines 257 to 258
Copy link
Member

@vepadulano vepadulano Nov 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
fRoots.clear();
fRoots.reserve(n);

These are not needed now


Expand All @@ -271,6 +273,5 @@ const std::vector< std::complex <double> > & Polynomial::FindNumRoots(){
return fRoots;
}


} // namespace Math
} // namespace ROOT