[MathMore] allow using several methods of Polynomial.h in const context

    - Remove unused `fRoots` data member (it was cleared before each root
      computation and provided no persistent state).
    - Return polynomial roots by value instead of via internal storage.
    - Mark the affected query/utility methods as `const` to make the class
      usable in const contexts.
    - Apply clang-format to the modified sections.

Author: meiyasan

math/mathmore/inc/Math/Polynomial.h

@@ -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();
+   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();
+   std::vector<std::complex<double>> FindNumRoots() const;
 
    /**
       Order of Polynomial
@@ -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

math/mathmore/src/Polynomial.cxx

@@ -147,43 +147,40 @@ IGenFunction * Polynomial::Clone() const {
     return f;
 }
 
+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--;
-    }
-
-    fRoots.clear();
-    fRoots.reserve(n);
+   }
 
+   std::vector<std::complex<double>> roots;
 
-    if (n == 0) {
-      return fRoots;
-    }
-    else if (n == 1 ) {
-      if ( Parameters()[1] == 0) return fRoots;
+   if (n == 0) {
+      return roots;
+   } else if (n == 1) {
+      if (Parameters()[1] == 0)
+         return roots;
       double r = - Parameters()[0]/ Parameters()[1];
-      fRoots.push_back( std::complex<double> ( r, 0.0) );
-    }
-    // quadratic equations
-    else if (n == 2 ) {
+      roots.push_back(std::complex<double>(r, 0.0));
+   }
+   // quadratic equations
+   else if (n == 2) {
       gsl_complex z1, z2;
       int nr = gsl_poly_complex_solve_quadratic(Parameters()[2], Parameters()[1], Parameters()[0], &z1, &z2);
       if ( nr != 2) {
 #ifdef DEBUG
          std::cout << "Polynomial quadratic ::-  FAILED to find roots" << std::endl;
 #endif
-         return fRoots;
+         return roots;
       }
-      fRoots.push_back(std::complex<double>(z1.dat[0],z1.dat[1]) );
-      fRoots.push_back(std::complex<double>(z2.dat[0],z2.dat[1]) );
-    }
-    // cubic equations
-    else if (n == 3 ) {
+      roots.push_back(std::complex<double>(z1.dat[0], z1.dat[1]));
+      roots.push_back(std::complex<double>(z2.dat[0], z2.dat[1]));
+   }
+   // cubic equations
+   else if (n == 3) {
       gsl_complex  z1, z2, z3;
       // renormmalize params in this case
       double w = Parameters()[3];
@@ -195,15 +192,14 @@ const std::vector< std::complex <double> > &  Polynomial::FindRoots(){
 #ifdef DEBUG
          std::cout << "Polynomial  cubic::-  FAILED to find roots" << std::endl;
 #endif
-         return fRoots;
-
+         return roots;
       }
-      fRoots.push_back(std::complex<double> (z1.dat[0],z1.dat[1]) );
-      fRoots.push_back(std::complex<double> (z2.dat[0],z2.dat[1]) );
-      fRoots.push_back(std::complex<double> (z3.dat[0],z3.dat[1]) );
-    }
-    // quartic equations
-    else if (n == 4 ) {
+      roots.push_back(std::complex<double>(z1.dat[0], z1.dat[1]));
+      roots.push_back(std::complex<double>(z2.dat[0], z2.dat[1]));
+      roots.push_back(std::complex<double>(z3.dat[0], z3.dat[1]));
+   }
+   // quartic equations
+   else if (n == 4) {
       // quartic eq.
       gsl_complex  z1, z2, z3, z4;
       // renormalize params in this case
@@ -217,60 +213,55 @@ const std::vector< std::complex <double> > &  Polynomial::FindRoots(){
 #ifdef DEBUG
          std::cout << "Polynomial quartic ::-  FAILED to find roots" << std::endl;
 #endif
-         return fRoots;
+         return roots;
       }
-      fRoots.push_back(std::complex<double> (z1.dat[0],z1.dat[1]) );
-      fRoots.push_back(std::complex<double> (z2.dat[0],z2.dat[1]) );
-      fRoots.push_back(std::complex<double> (z3.dat[0],z3.dat[1]) );
-      fRoots.push_back(std::complex<double> (z4.dat[0],z4.dat[1]) );
-    }
-    else {
-    // for higher order polynomial use numerical fRoots
-      FindNumRoots();
-    }
-
-    return fRoots;
-
-  }
-
-
-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;
+      roots.push_back(std::complex<double>(z1.dat[0], z1.dat[1]));
+      roots.push_back(std::complex<double>(z2.dat[0], z2.dat[1]));
+      roots.push_back(std::complex<double>(z3.dat[0], z3.dat[1]));
+      roots.push_back(std::complex<double>(z4.dat[0], z4.dat[1]));
+   } else {
+      // for higher order polynomial use numerical roots
+      return FindNumRoots();
+   }
+
+   return roots;
 }
-const std::vector< std::complex <double> > &  Polynomial::FindNumRoots(){
-
-
-    // check if order is correct
-    unsigned int n = fOrder;
-    while ( Parameters()[n] == 0 ) {
-      n--;
-    }
-    fRoots.clear();
-    fRoots.reserve(n);
-
 
-    if (n == 0) {
-      return fRoots;
-    }
-
-    gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc( n + 1);
-    std::vector<double> z(2*n);
-    int status = gsl_poly_complex_solve ( Parameters(), n+1, w, &z.front() );
-    gsl_poly_complex_workspace_free(w);
-    if (status != GSL_SUCCESS) return fRoots;
-    for (unsigned int i = 0; i < n; ++i)
-      fRoots.push_back(std::complex<double> (z[2*i],z[2*i+1] ) );
+std::vector<double> Polynomial::FindRealRoots() const
+{
+   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 fRoots;
+   return realRoots;
 }
+std::vector<std::complex<double>> Polynomial::FindNumRoots() const
+{
 
+   // check if order is correct
+   unsigned int n = fOrder;
+   while (Parameters()[n] == 0) {
+      n--;
+   }
+
+   std::vector<std::complex<double>> roots;
+   if (n == 0) {
+      return roots;
+   }
+
+   gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(n + 1);
+   std::vector<double> z(2 * n);
+   int status = gsl_poly_complex_solve(Parameters(), n + 1, w, &z.front());
+   gsl_poly_complex_workspace_free(w);
+   if (status != GSL_SUCCESS)
+      return roots;
+   for (unsigned int i = 0; i < n; ++i)
+      roots.push_back(std::complex<double>(z[2 * i], z[2 * i + 1]));
+
+   return roots;
+}
 
 } // namespace Math
 } // namespace ROOT