Mapa de calor de la raíz polinomial Littlewood -- ++ camp codereview Relacionados El problema

Littlewood polynomial root heatmap


6
vote

problema

Español

Un poliNomial de littlewood es un polinomio donde cada coeficiente es éter -1 o 1 y cuando las raíces complejas de ellos producen una buena imagen. Como resultado, decidí hacer un programa en C ++ que parcelas un mapa de calor de las raíces complejas de los polinomios de littlewood elegidos aleatoriamente. Además, al principio pensé que sería tan usando GSL (Biblioteca Científica GNU), ya que ya tenía un solucionador polinomial y libpng. Sin embargo, después de implementar una versión básica de él utilizando GSL, me di cuenta de que GSL era lento para usar. En consecuencia, eso significaba que tuve que encontrar otra biblioteca polinomial y después de un poco de búsqueda, encontré esto https://www.codeproject.com/articles/674149/a-real-polinomial-class-with-root-finder . Luego, después de limpiar esa biblioteca y aprender a usar Libpng, fue mayormente suave navegando. Sin embargo, mis principales preocupaciones son que podría seguir siendo formas de mejorar el desempeño del código que no conozco y que aún podría mejorar la calidad del código de manera que no conozco.

png.hh

  #ifndef PNG_HH #define PNG_HH #include <png.h>  namespace png {     void write_image(char const *filename, std::uint8_t const *image_data, std::uint32_t image_width, std::uint32_t image_height)     {         /* create a zeroed out png_image struct */         png_image output_png;         std::memset(&output_png, 0, sizeof(output_png));         output_png.version = PNG_IMAGE_VERSION;         output_png.format = PNG_FORMAT_GRAY;         output_png.width = image_width;         output_png.height = image_height;          /* write the png file */         png_image_write_to_file(&output_png, filename, 0, image_data, image_height, nullptr);          /* cleanup */         png_image_free(&output_png);     } } #endif   

PolynomialRootFinder.hh

  //======================================================================= // Copyright (C) 2003-2013 William Hallahan // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without restriction, // including without limitation the rights to use, copy, modify, merge, // publish, distribute, sublicense, and/or sell copies of the Software, // and to permit persons to whom the Software is furnished to do so, // subject to the following conditions: // // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. // // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT // HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, // WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR // OTHER DEALINGS IN THE SOFTWARE. //=======================================================================  //********************************************************************** //  File: PolynomialRootFinder.h //  Author: Bill Hallahan //  Date: January 30, 2003 // //  Abstract: // //    This file contains the definition for class PolynomialRootFinder. // //**********************************************************************  #ifndef POLYNOMIALROOTFINDER_H #define POLYNOMIALROOTFINDER_H  #include <array>  //====================================================================== //  Class definition. //======================================================================   template<std::int32_t degree> struct PolynomialRootFinder {     std::array<double, degree + 1> m_p_vector;     std::array<double, degree + 1> m_qp_vector;     std::array<double, degree + 1> m_k_vector;     std::array<double, degree + 1> m_qk_vector;     std::array<double, degree + 1> m_svk_vector;      std::int32_t m_n;     std::int32_t m_n_plus_one;     double m_real_s;     double m_imag_s;     double m_u;     double m_v;     double m_a;     double m_b;     double m_c;     double m_d;     double m_a1;     double m_a2;     double m_a3;     double m_a6;     double m_a7;     double m_e;     double m_f;     double m_g;     double m_h;     double m_real_sz;     double m_imag_sz;     double m_real_lz;     double m_imag_lz;     double m_are;     double m_mre;      enum class RootStatus_T     {         SUCCESS,         LEADING_COEFFICIENT_IS_ZERO,         SCALAR_VALUE_HAS_NO_ROOTS,         FAILED_TO_CONVERGE     };      PolynomialRootFinder::RootStatus_T FindRoots(double *coefficient_ptr,         double *real_zero_vector_ptr,         double *imaginary_zero_vector_ptr,         std::int32_t *number_of_roots_found_ptr = 0);      std::int32_t Fxshfr(std::int32_t l2var);      std::int32_t QuadraticIteration(double uu, double vv);      std::int32_t RealIteration(double &sss, std::int32_t &flag);      std::int32_t CalcSc();      void NextK(std::int32_t itype);      void Newest(std::int32_t itype, double &uu, double &vv);      void QuadraticSyntheticDivision(std::int32_t n_plus_one,         double u,         double v,         double *p_ptr,         double *q_ptr,         double &a,         double &b);      void SolveQuadraticEquation(double a,         double b,         double c,         double &sr,         double &si,         double &lr,         double &li); };  #include <cmath> #include <float.h>  namespace {     //------------------------------------------------------------------     //  The following machine constants are used in this method.     //     //    f_BASE  The base of the floating postd::int32_t number system used.     //     //    f_ETA  The maximum relative representation error which     //           can be described as the smallest positive floating     //           postd::int32_t number such that 1.0 + f_ETA is greater than 1.0.     //     //    f_MAXIMUM_FLOAT  The largest floating postd::int32_t number.     //     //    f_MINIMUM_FLOAT  The smallest positive floating postd::int32_t number.     //     //------------------------------------------------------------------      constexpr float f_BASE = 2.0;     constexpr float f_ETA = FLT_EPSILON;     constexpr float f_ETA_N = (10.0f) * f_ETA;     constexpr float f_ETA_N_SQUARED = (100.0f) * f_ETA;     constexpr float f_MAXIMUM_FLOAT = FLT_MAX;     constexpr float f_MINIMUM_FLOAT = FLT_MIN;     constexpr float f_XX_INITIAL_VALUE = (0.70710678f);     constexpr float f_COSR_INITIAL_VALUE = (-0.069756474f);     constexpr float f_SINR_INITIAL_VALUE = (0.99756405f); };  //====================================================================== //  Member Function: PolynomialRootFinder::FindRoots // //  Abstract: // //    This method determines the roots of a polynomial which //    has real coefficients. This code is based on FORTRAN //    code published in reference [1]. The method is based on //    an algorithm the three-stage algorithm described in //    Jenkins and Traub [2]. // // 1. "Collected Algorithms from ACM, Volume III", Algorithms 493-545 //    1983. (The root finding algorithms is number 493) // // 2. Jenkins, M. A. and Traub, J. F., "A three-stage algorithm for //    real polynomials using quadratic iteration", SIAM Journal of //    Numerical Analysis, 7 (1970), 545-566 // // 3. Jenkins, M. A. and Traub, J. F., "Principles for testing //    polynomial zerofinding programs", ACM TOMS 1, //    1 (March 1975), 26-34 // // //  Input: // //    All vectors below must be at least a length equal to degree + 1. // //    coefficicent_ptr           A double precision vector that contains //                               the polynomial coefficients in order //                               of increasing power. // //    degree                     The degree of the polynomial. // //    real_zero_vector_ptr       A double precision vector that will //                               contain the real parts of the roots //                               of the polynomial when this method //                               returns. // //    imaginary_zero_vector_ptr  A double precision vector that will //                               contain the real parts of the roots //                               of the polynomial when this method //                               returns. // //    number_of_roots_found_ptr  A postd::int32_ter to an std::int32_teger that will //                               equal the number of roots found when //                               this method returns. If the method //                               returns SUCCESS then this value will //                               always equal the degree of the //                               polynomial. // //  Return Value: // //    The function returns an enum value of type //    'PolynomialRootFinder::RootStatus_T'. // //======================================================================  template<std::int32_t degree> typename PolynomialRootFinder<degree>::RootStatus_T PolynomialRootFinder<degree>::FindRoots(     double *coefficient_vector_ptr,     double *real_zero_vector_ptr,     double *imaginary_zero_vector_ptr,     std::int32_t *number_of_roots_found_ptr) {     //--------------------------------------------------------------     //  The algorithm fails if the polynomial is not at least     //  degree on or the leading coefficient is zero.     //--------------------------------------------------------------      PolynomialRootFinder::RootStatus_T status;     //--------------------------------------------------------------     //  Allocate temporary vectors used to find the roots.     //--------------------------------------------------------------       std::array<double, degree + 1> temp_vector;     std::array<double, degree + 1> pt_vector;      //--------------------------------------------------------------     //  m_are and m_mre refer to the unit error in + and *     //  respectively. they are assumed to be the same as     //  f_ETA.     //--------------------------------------------------------------      m_are = f_ETA;     m_mre = f_ETA;     double lo = f_MINIMUM_FLOAT / f_ETA;      //--------------------------------------------------------------     // Initialization of constants for shift rotation.     //--------------------------------------------------------------      double xx = f_XX_INITIAL_VALUE;     double yy = -xx;     double cosr = f_COSR_INITIAL_VALUE;     double sinr = f_SINR_INITIAL_VALUE;     m_n = degree;     m_n_plus_one = m_n + 1;      //--------------------------------------------------------------     //  Make a copy of the coefficients in reverse order.     //--------------------------------------------------------------      std::int32_t ii = 0;      for (ii = 0; ii < m_n_plus_one; ++ii) {         m_p_vector[m_n - ii] = coefficient_vector_ptr[ii];     }      //--------------------------------------------------------------     //  Assume failure. The status is set to SUCCESS if all     //  the roots are found.     //--------------------------------------------------------------      status = PolynomialRootFinder::RootStatus_T::FAILED_TO_CONVERGE;      //--------------------------------------------------------------     //  If there are any zeros at the origin, remove them.     //--------------------------------------------------------------      std::int32_t jvar = 0;      while (m_p_vector[m_n] == 0.0) {         jvar = degree - m_n;         real_zero_vector_ptr[jvar] = 0.0;         imaginary_zero_vector_ptr[jvar] = 0.0;         m_n_plus_one = m_n_plus_one - 1;         m_n = m_n - 1;     }      //--------------------------------------------------------------     //  Loop and find polynomial zeros. In the original algorithm     //  this loop was an infinite loop. Testing revealed that the     //  number of main loop iterations to solve a polynomial of a     //  particular degree is usually about half the degree.     //  We loop twice that to make sure the solution is found.     //  (This should be revisited as it might preclude solving     //  some large polynomials.)     //--------------------------------------------------------------       for (std::int32_t count = 0; count < degree; ++count) {         //----------------------------------------------------------         //  Check for less than 2 zeros to finish the solution.         //----------------------------------------------------------          if (m_n <= 2) {             if (m_n > 0) {                 //--------------------------------------------------                 //  Calculate the final zero or pair of zeros.                 //--------------------------------------------------                  if (m_n == 1) {                     real_zero_vector_ptr[degree - 1] =                         -m_p_vector[1] / m_p_vector[0];                      imaginary_zero_vector_ptr[degree - 1] = 0.0;                 } else {                     SolveQuadraticEquation(                         m_p_vector[0],                         m_p_vector[1],                         m_p_vector[2],                         real_zero_vector_ptr[degree - 2],                         imaginary_zero_vector_ptr[degree - 2],                         real_zero_vector_ptr[degree - 1],                         imaginary_zero_vector_ptr[degree - 1]);                 }             }              m_n = 0;             status = PolynomialRootFinder::RootStatus_T::SUCCESS;             break;         } else {             //------------------------------------------------------             //  Find largest and smallest moduli of coefficients.             //------------------------------------------------------              double max = 0.0;             double min = f_MAXIMUM_FLOAT;             double xvar;              for (ii = 0; ii < m_n_plus_one; ++ii) {                 xvar = (double)(::fabs((double)(m_p_vector[ii])));                  if (xvar > max) {                     max = xvar;                 }                  if ((xvar != 0.0) && (xvar < min)) {                     min = xvar;                 }             }              //------------------------------------------------------             //  Scale if there are large or very small coefficients.             //  Computes a scale factor to multiply the coefficients             //  of the polynomial. The scaling is done to avoid             //  overflow and to avoid undetected underflow from             //  std::int32_terfering with the convergence criterion.             //  The factor is a power of the base.             //------------------------------------------------------              bool do_scaling_flag = false;             double sc = lo / min;              if (sc <= 1.0) {                 do_scaling_flag = f_MAXIMUM_FLOAT / sc < max;             } else {                 do_scaling_flag = max < 10.0;                  if (!do_scaling_flag) {                     if (sc == 0.0) {                         sc = f_MINIMUM_FLOAT;                     }                 }             }              //------------------------------------------------------             //  Conditionally scale the data.             //------------------------------------------------------              if (do_scaling_flag) {                 std::int32_t lvar = (std::int32_t)(::log(sc) / ::log(f_BASE) + 0.5);                 double factor = ::pow((double)(f_BASE * 1.0), double(lvar));                  if (factor != 1.0) {                     for (ii = 0; ii < m_n_plus_one; ++ii) {                         m_p_vector[ii] = factor * m_p_vector[ii];                     }                 }             }              //------------------------------------------------------             //  Compute lower bound on moduli of zeros.             //------------------------------------------------------              for (ii = 0; ii < m_n_plus_one; ++ii) {                 pt_vector[ii] = (double)(::fabs((double)(m_p_vector[ii])));             }              pt_vector[m_n] = -pt_vector[m_n];              //------------------------------------------------------             //  Compute upper estimate of bound.             //------------------------------------------------------              xvar = (double)                 (::exp((::log(-pt_vector[m_n]) - ::log(pt_vector[0]))                     / (double)(m_n)));              //------------------------------------------------------             //  If newton step at the origin is better, use it.             //------------------------------------------------------              double xm;              if (pt_vector[m_n - 1] != 0.0) {                 xm = -pt_vector[m_n] / pt_vector[m_n - 1];                  if (xm < xvar) {                     xvar = xm;                 }             }              //------------------------------------------------------             //  Chop the std::int32_terval (0, xvar) until ff <= 0             //------------------------------------------------------              double ff;              for (;;) {                 xm = (double)(xvar * 0.1);                 ff = pt_vector[0];                  for (ii = 1; ii < m_n_plus_one; ++ii) {                     ff = ff * xm + pt_vector[ii];                 }                 if (ff <= 0.0) {                     ff = 0.0;                     break;                 }                  xvar = xm;             }              double dx = xvar;              //------------------------------------------------------             //  Do newton iteration until xvar converges to two             //  decimal places.             //------------------------------------------------------              for (;;) {                 if ((double)(::fabs(dx / xvar)) <= 0.005) {                     break;                 }                  ff = pt_vector[0];                  double df = ff;                  for (ii = 1; ii < m_n; ++ii) {                     ff = ff * xvar + pt_vector[ii];                      df = df * xvar + ff;                  }                  ff = ff * xvar + pt_vector[m_n];                  dx = ff / df;                   xvar = xvar - dx;             }              double bnd = xvar;              //------------------------------------------------------             //  Compute the derivative as the std::int32_tial m_k_vector             //  polynomial and do 5 steps with no shift.             //------------------------------------------------------              std::int32_t n_minus_one = m_n - 1;              for (ii = 1; ii < m_n; ++ii) {                 m_k_vector[ii] =                     (double)(m_n - ii) * m_p_vector[ii] / (double)(m_n);             }              m_k_vector[0] = m_p_vector[0];             double aa = m_p_vector[m_n];             double bb = m_p_vector[m_n - 1];             bool zerok_flag = m_k_vector[m_n - 1] == 0.0;              std::int32_t jj = 0;              for (jj = 1; jj <= 5; ++jj) {                 double cc = m_k_vector[m_n - 1];                  if (zerok_flag) {                     //----------------------------------------------                     //  Use unscaled form of recurrence.                     //----------------------------------------------                      for (jvar = n_minus_one; jvar > 0; --jvar) {                         m_k_vector[jvar] = m_k_vector[jvar - 1];                     }                      m_k_vector[0] = 0.0;                     zerok_flag = m_k_vector[m_n - 1] == 0.0;                 } else {                     //----------------------------------------------                     //  Use scaled form of recurrence if value                     //  of m_k_vector at 0 is nonzero.                     //----------------------------------------------                      double tvar = -aa / cc;                      for (jvar = n_minus_one; jvar > 0; --jvar) {                         m_k_vector[jvar] =                             tvar * m_k_vector[jvar - 1] + m_p_vector[jvar];                     }                      m_k_vector[0] = m_p_vector[0];                     zerok_flag =                         ::fabs(m_k_vector[m_n - 1]) <= ::fabs(bb) * f_ETA_N;                 }             }              //------------------------------------------------------             //  Save m_k_vector for restarts with new shifts.             //------------------------------------------------------              for (ii = 0; ii < m_n; ++ii) {                 temp_vector[ii] = m_k_vector[ii];             }              //------------------------------------------------------             //  Loop to select the quadratic corresponding to             //   each new shift.             //------------------------------------------------------              std::int32_t cnt = 0;              for (cnt = 1; cnt <= 20; ++cnt) {                 //--------------------------------------------------                 //  Quadratic corresponds to a double shift to a                 //  non-real postd::int32_t and its complex conjugate. The                 //  postd::int32_t has modulus 'bnd' and amplitude rotated                 //  by 94 degrees from the previous shift.                 //--------------------------------------------------                  double xxx = cosr * xx - sinr * yy;                 yy = sinr * xx + cosr * yy;                 xx = xxx;                 m_real_s = bnd * xx;                 m_imag_s = bnd * yy;                 m_u = -2.0 * m_real_s;                 m_v = bnd;                  //--------------------------------------------------                 //  Second stage calculation, fixed quadratic.                 //  Variable nz will contain the number of                 //   zeros found when function Fxshfr() returns.                 //--------------------------------------------------                  std::int32_t nz = Fxshfr(20 * cnt);                  if (nz != 0) {                     //----------------------------------------------                     //  The second stage jumps directly to one of                     //  the third stage iterations and returns here                     //  if successful. Deflate the polynomial,                     //  store the zero or zeros and return to the                     //  main algorithm.                     //----------------------------------------------                      jvar = degree - m_n;                     real_zero_vector_ptr[jvar] = m_real_sz;                     imaginary_zero_vector_ptr[jvar] = m_imag_sz;                     m_n_plus_one = m_n_plus_one - nz;                     m_n = m_n_plus_one - 1;                      for (ii = 0; ii < m_n_plus_one; ++ii) {                         m_p_vector[ii] = m_qp_vector[ii];                     }                      if (nz != 1) {                         real_zero_vector_ptr[jvar + 1] = m_real_lz;                         imaginary_zero_vector_ptr[jvar + 1] = m_imag_lz;                     }                      break;                      //----------------------------------------------                     //  If the iteration is unsuccessful another                     //  quadratic is chosen after restoring                     //  m_k_vector.                     //----------------------------------------------                 }                  for (ii = 0; ii < m_n; ++ii) {                     m_k_vector[ii] = temp_vector[ii];                 }             }         }     }      //--------------------------------------------------------------     //  If no convergence with 20 shifts then adjust the degree     //  for the number of roots found.     //--------------------------------------------------------------      if (number_of_roots_found_ptr != 0) {         *number_of_roots_found_ptr = degree - m_n;     }      return status; }  //====================================================================== //  Computes up to l2var fixed shift m_k_vector polynomials, //  testing for convergence in the linear or quadratic //  case. initiates one of the variable shift //  iterations and returns with the number of zeros //  found. // //    l2var  An std::int32_teger that is the limit of fixed shift steps. // //  Return Value: //    nz   An std::int32_teger that is the number of zeros found. //======================================================================  template<std::int32_t degree> std::int32_t PolynomialRootFinder<degree>::Fxshfr(std::int32_t l2var) {     //------------------------------------------------------------------     //  Evaluate polynomial by synthetic division.     //------------------------------------------------------------------      QuadraticSyntheticDivision(m_n_plus_one,         m_u,         m_v,         m_p_vector.data(),         m_qp_vector.data(),         m_a,         m_b);      std::int32_t itype = CalcSc();      std::int32_t nz = 0;     float betav = 0.25;     float betas = 0.25;     float oss = (float)(m_real_s);     float ovv = (float)(m_v);     float ots;     float otv;     double ui = 0.0;     double vi = 0.0;     double svar;       for (std::int32_t jvar = 1; jvar <= l2var; ++jvar) {         //--------------------------------------------------------------         //  Calculate next m_k_vector polynomial and estimate m_v.         //--------------------------------------------------------------          NextK(itype);         itype = CalcSc();         Newest(itype, ui, vi);         float vv = (float)(vi);          //--------------------------------------------------------------         //  Estimate svar         //--------------------------------------------------------------          float ss = 0.0;          if (m_k_vector[m_n - 1] != 0.0) {             ss = (float)(-m_p_vector[m_n] / m_k_vector[m_n - 1]);         }          float tv = 1.0;         float ts = 1.0;          if ((jvar != 1) && (itype != 3)) {             //----------------------------------------------------------             //  Compute relative measures of convergence of             //  svar and m_v sequences.             //----------------------------------------------------------              if (vv != 0.0) {                 tv = (float)(::fabs((vv - ovv) / vv));             }              if (ss != 0.0) {                 ts = (float)(::fabs((ss - oss) / ss));             }              //----------------------------------------------------------             //  If decreasing, multiply two most recent convergence             //  measures.             //----------------------------------------------------------              float tvv = 1.0;              if (tv < otv) {                 tvv = tv * otv;             }              float tss = 1.0;              if (ts < ots) {                 tss = ts * ots;             }              //----------------------------------------------------------             //  Compare with convergence criteria.             //----------------------------------------------------------              bool vpass_flag = tvv < betav;             bool spass_flag = tss < betas;              if (spass_flag || vpass_flag) {                 //------------------------------------------------------                 //  At least one sequence has passed the convergence                 //  test. Store variables before iterating.                 //------------------------------------------------------                  double svu = m_u;                 double svv = m_v;                 std::int32_t ii = 0;                  for (ii = 0; ii < m_n; ++ii) {                     m_svk_vector[ii] = m_k_vector[ii];                 }                  svar = ss;                  //------------------------------------------------------                 //  Choose iteration according to the fastest                 //  converging sequence.                 //------------------------------------------------------                  bool vtry_flag = false;                 bool stry_flag = false;                 bool exit_outer_loop_flag = false;                  bool start_with_real_iteration_flag =                     (spass_flag && ((!vpass_flag) || (tss < tvv)));                  do {                     if (!start_with_real_iteration_flag) {                         nz = QuadraticIteration(ui, vi);                          if (nz > 0) {                             exit_outer_loop_flag = true;                             break;                         }                          //----------------------------------------------                         //  Quadratic iteration has failed. flag                         //  that it has been tried and decrease                         //  the convergence criterion.                         //----------------------------------------------                          vtry_flag = true;                         betav = (float)(betav * 0.25);                     }                      //--------------------------------------------------                     //  Try linear iteration if it has not been                     //  tried and the svar sequence is converging.                     //--------------------------------------------------                      if (((!stry_flag) && spass_flag)                         || start_with_real_iteration_flag) {                         if (!start_with_real_iteration_flag) {                             for (ii = 0; ii < m_n; ++ii) {                                 m_k_vector[ii] = m_svk_vector[ii];                             }                         } else {                             start_with_real_iteration_flag = false;                         }                          std::int32_t iflag = 0;                          nz = RealIteration(svar, iflag);                          if (nz > 0) {                             exit_outer_loop_flag = true;                             break;                         }                          //----------------------------------------------                         //  Linear iteration has failed. Flag that                         //  it has been tried and decrease the                         //  convergence criterion.                         //----------------------------------------------                          stry_flag = true;                         betas = (float)(betas * 0.25);                          if (iflag != 0) {                             //------------------------------------------                             //  If linear iteration signals an almost                             //  double real zero attempt quadratic                             //  iteration.                             //------------------------------------------                              ui = -(svar + svar);                             vi = svar * svar;                              continue;                         }                     }                      //--------------------------------------------------                     //  Restore variables                     //--------------------------------------------------                      m_u = svu;                     m_v = svv;                      for (ii = 0; ii < m_n; ++ii) {                         m_k_vector[ii] = m_svk_vector[ii];                     }                      //----------------------------------------------                     //  Try quadratic iteration if it has not been                     //  tried and the m_v sequence is converging.                     //----------------------------------------------                 } while (vpass_flag && (!vtry_flag));                  if (exit_outer_loop_flag) {                     break;                 }                  //------------------------------------------------------                 //  Recompute m_qp_vector and scalar values to                 //  continue the second stage.                 //------------------------------------------------------                  QuadraticSyntheticDivision(m_n_plus_one,                     m_u,                     m_v,                     m_p_vector.data(),                     m_qp_vector.data(),                     m_a,                     m_b);                  itype = CalcSc();             }         }          ovv = vv;         oss = ss;         otv = tv;         ots = ts;     }      return nz; }  //====================================================================== //  Variable-shift m_k_vector-polynomial iteration for //  a quadratic factor converges only if the zeros are //  equimodular or nearly so. // //    uu  Coefficients of starting quadratic //    vv  Coefficients of starting quadratic // //  Return value: //    nz  The number of zeros found. //======================================================================  template<std::int32_t degree> std::int32_t PolynomialRootFinder<degree>::QuadraticIteration(double uu, double vv) {     //------------------------------------------------------------------     //  Main loop     //------------------------------------------------------------------      double ui = 0.0;     double vi = 0.0;     float omp = 0.0F;     float relstp = 0.0F;     std::int32_t itype = 0;     bool tried_flag = false;     std::int32_t jvar = 0;     std::int32_t nz = 0;     m_u = uu;     m_v = vv;      for(;;) {         SolveQuadraticEquation(1.0,             m_u,             m_v,             m_real_sz,             m_imag_sz,             m_real_lz,             m_imag_lz);          //--------------------------------------------------------------         //  Return if roots of the quadratic are real and not close         //  to multiple or nearly equal and  of opposite sign.         //--------------------------------------------------------------          if (::fabs(::fabs(m_real_sz) - ::fabs(m_real_lz)) > 0.01 * ::fabs(m_real_lz)) {             break;         }          //--------------------------------------------------------------         //  Evaluate polynomial by quadratic synthetic division.         //------------------------------------------------------------------          QuadraticSyntheticDivision(m_n_plus_one,             m_u,             m_v,             m_p_vector.data(),             m_qp_vector.data(),             m_a,             m_b);          float mp = (float)(::fabs(m_a - m_real_sz * m_b) + ::fabs(m_imag_sz * m_b));          //--------------------------------------------------------------         //  Compute a rigorous  bound on the rounding error in         //  evaluting m_p_vector.         //--------------------------------------------------------------          float zm = (float)(::sqrt((float)(::fabs((float)(m_v)))));         float ee = (float)(2.0 * (float)(::fabs((float)(m_qp_vector[0]))));         float tvar = (float)(-m_real_sz * m_b);         std::int32_t ii = 0;          for (ii = 1; ii < m_n; ++ii) {             ee = ee * zm + (float)(::fabs((float)(m_qp_vector[ii])));         }          ee = ee * zm + (float)(::fabs((float)(m_a)+tvar));         ee = (float)((5.0 * m_mre + 4.0 * m_are) * ee             - (5.0 * m_mre + 2.0 * m_are) * ((float)(::fabs((float)(m_a)+tvar)) + (float)(::fabs((float)(m_b))) * zm)             + 2.0 * m_are * (float)(::fabs(tvar)));          //--------------------------------------------------------------         //  Iteration has converged sufficiently if the polynomial         //  value is less than 20 times this bound.         //--------------------------------------------------------------          if (mp <= 20.0 * ee) {             nz = 2;             break;         }          jvar = jvar + 1;          //--------------------------------------------------------------         //  Stop iteration after 20 steps.         //--------------------------------------------------------------          if (jvar > 20) {             break;         }          if ((jvar >= 2) && ((relstp <= 0.01)             && (mp >= omp) && (!tried_flag))) {             //----------------------------------------------------------             //  A cluster appears to be stalling the convergence.             //  Five fixed shift steps are taken with a m_u, m_v             //  close to the cluster.             //----------------------------------------------------------              if (relstp < f_ETA) {                 relstp = f_ETA;             }              relstp = (float)(::sqrt(relstp));             m_u = m_u - m_u * relstp;             m_v = m_v + m_v * relstp;              QuadraticSyntheticDivision(m_n_plus_one,                 m_u,                 m_v,                 m_p_vector.data(),                 m_qp_vector.data(),                 m_a,                 m_b);              for (ii = 0; ii < 5; ++ii) {                 itype = CalcSc();                 NextK(itype);             }              tried_flag = true;             jvar = 0;         }          omp = mp;          //--------------------------------------------------------------         //  Calculate next m_k_vector polynomial and         //  new m_u and m_v.         //--------------------------------------------------------------          itype = CalcSc();         NextK(itype);         itype = CalcSc();         Newest(itype, ui, vi);          //--------------------------------------------------------------         //  If vi is zero the iteration is not converging.         //--------------------------------------------------------------          if (vi == 0.0) {             break;         }          relstp = (float)(::fabs((vi - m_v) / vi));         m_u = ui;         m_v = vi;     }      return nz; }  //====================================================================== //  Variable-shift h polynomial iteration for a real zero. // //    sss      Starting iterate //    flag     Flag to indicate a pair of zeros near real axis. // //  Return Value: //     Number of zero found. //======================================================================  template<std::int32_t degree> std::int32_t PolynomialRootFinder<degree>::RealIteration(double &sss, std::int32_t &flag) {     //------------------------------------------------------------------     //  Main loop     //------------------------------------------------------------------      double tvar = 0.0;     float omp = 0.0F;     std::int32_t nz = 0;     flag = 0;     std::int32_t jvar = 0;     double svar = sss;      for(;;) {         double pv = m_p_vector[0];          //--------------------------------------------------------------         //  Evaluate m_p_vector at svar         //--------------------------------------------------------------          m_qp_vector[0] = pv;         std::int32_t ii = 0;          for (ii = 1; ii < m_n_plus_one; ++ii) {             pv = pv * svar + m_p_vector[ii];             m_qp_vector[ii] = pv;         }          float mp = (float)(::fabs(pv));          //--------------------------------------------------------------         //  Compute a rigorous bound on the error in evaluating p         //--------------------------------------------------------------          double ms = (double)(::fabs(svar));         double ee = (m_mre / (m_are + m_mre)) * (double)(::fabs((double)(m_qp_vector[0])));          for (ii = 1; ii < m_n_plus_one; ++ii) {             ee = ee * ms + (float)(::fabs((double)(m_qp_vector[ii])));         }          //--------------------------------------------------------------         //  Iteration has converged sufficiently if the         //  polynomial value is less than 20 times this bound.         //--------------------------------------------------------------          if (mp <= 20.0 * ((m_are + m_mre) * ee - m_mre * mp)) {             nz = 1;             m_real_sz = svar;             m_imag_sz = 0.0;             break;         }          jvar = jvar + 1;          //--------------------------------------------------------------         //  Stop iteration after 10 steps.         //--------------------------------------------------------------          if (jvar > 10) {             break;         }          if ((jvar >= 2)             && ((::fabs(tvar) <= 0.001 * ::fabs(svar - tvar))                 && (mp > omp))) {             //----------------------------------------------------------             //  A cluster of zeros near the real axis has been             //  encountered. Return with flag set to initiate             //  a quadratic iteration.             //----------------------------------------------------------              flag = 1;             sss = svar;             break;         }          //--------------------------------------------------------------         //  Return if the polynomial value has increased significantly.         //--------------------------------------------------------------          omp = mp;          //--------------------------------------------------------------         //  Compute t, the next polynomial, and the new iterate.         //--------------------------------------------------------------          double kv = m_k_vector[0];         m_qk_vector[0] = kv;          for (ii = 1; ii < m_n; ++ii) {             kv = kv * svar + m_k_vector[ii];             m_qk_vector[ii] = kv;         }          if (::fabs(kv) <= ::fabs(m_k_vector[m_n - 1]) * f_ETA_N) {             m_k_vector[0] = 0.0;              for (ii = 1; ii < m_n; ++ii) {                 m_k_vector[ii] = m_qk_vector[ii - 1];             }         } else {             //----------------------------------------------------------             //  Use the scaled form of the recurrence if the             //  value of m_k_vector at svar is non-zero.             //----------------------------------------------------------              tvar = -pv / kv;             m_k_vector[0] = m_qp_vector[0];              for (ii = 1; ii < m_n; ++ii) {                 m_k_vector[ii] = tvar * m_qk_vector[ii - 1] + m_qp_vector[ii];             }         }          //--------------------------------------------------------------         //  Use unscaled form.         //--------------------------------------------------------------          kv = m_k_vector[0];          for (ii = 1; ii < m_n; ++ii) {             kv = kv * svar + m_k_vector[ii];         }          tvar = 0.0;          if (::fabs(kv) > ::fabs(m_k_vector[m_n - 1]) * f_ETA_N) {             tvar = -pv / kv;         }          svar = svar + tvar;     }      return nz; }  //====================================================================== //  This routine calculates scalar quantities used to compute //  the next m_k_vector polynomial and new estimates of the //  quadratic coefficients. // //  Return Value: //    type  std::int32_teger variable set here indicating how the //    calculations are normalized to avoid overflow. //======================================================================  template<std::int32_t degree> std::int32_t PolynomialRootFinder<degree>::CalcSc() {     //------------------------------------------------------------------     //  Synthetic division of m_k_vector by the quadratic 1, m_u, m_v.     //------------------------------------------------------------------      QuadraticSyntheticDivision(m_n,         m_u,         m_v,         m_k_vector.data(),         m_qk_vector.data(),         m_c,         m_d);      std::int32_t itype = 0;      if ((::fabs(m_c) <= ::fabs(m_k_vector[m_n - 1]) * f_ETA_N_SQUARED)         && (::fabs(m_d) <= ::fabs(m_k_vector[m_n - 2]) * f_ETA_N_SQUARED)) {         //--------------------------------------------------------------         //  itype == 3 Indicates the quadratic is almost a         //  factor of m_k_vector.         //--------------------------------------------------------------          itype = 3;     } else if (::fabs(m_d) >= ::fabs(m_c)) {         //--------------------------------------------------------------         //  itype == 2 Indicates that all formulas are divided by m_d.         //--------------------------------------------------------------          itype = 2;         m_e = m_a / m_d;         m_f = m_c / m_d;         m_g = m_u * m_b;         m_h = m_v * m_b;         m_a3 = (m_a + m_g) * m_e + m_h * (m_b / m_d);         m_a1 = m_b * m_f - m_a;         m_a7 = (m_f + m_u) * m_a + m_h;     } else {         //--------------------------------------------------------------         //  itype == 1 Indicates that all formulas are divided by m_c.         //--------------------------------------------------------------          itype = 1;         m_e = m_a / m_c;         m_f = m_d / m_c;         m_g = m_u * m_e;         m_h = m_v * m_b;         m_a3 = m_a * m_e + (m_h / m_c + m_g) * m_b;         m_a1 = m_b - m_a * (m_d / m_c);         m_a7 = m_a + m_g * m_d + m_h * m_f;     }      return itype;  }  //====================================================================== //  Computes the next k polynomials using scalars computed in CalcSc. //====================================================================== template<std::int32_t degree> void PolynomialRootFinder<degree>::NextK(std::int32_t itype) {     std::int32_t ii = 0;      if (itype == 3) {         //--------------------------------------------------------------         //  Use unscaled form of the recurrence if type is 3.         //--------------------------------------------------------------          m_k_vector[0] = 0.0;         m_k_vector[1] = 0.0;          for (ii = 2; ii < m_n; ++ii) {             m_k_vector[ii] = m_qk_vector[ii - 2];         }     } else {         double temp = m_a;          if (itype == 1) {             temp = m_b;         }          if (::fabs(m_a1) <= ::fabs(temp) * f_ETA_N) {             //----------------------------------------------------------             //  If m_a1 is nearly zero then use a special form of             //  the recurrence.             //----------------------------------------------------------              m_k_vector[0] = 0.0;             m_k_vector[1] = -m_a7 * m_qp_vector[0];              for (ii = 2; ii < m_n; ++ii) {                 m_k_vector[ii] = m_a3 * m_qk_vector[ii - 2] - m_a7 * m_qp_vector[ii - 1];             }         } else {             //----------------------------------------------------------             //  Use scaled form of the recurrence.             //----------------------------------------------------------              m_a7 = m_a7 / m_a1;             m_a3 = m_a3 / m_a1;             m_k_vector[0] = m_qp_vector[0];             m_k_vector[1] = m_qp_vector[1] - m_a7 * m_qp_vector[0];              for (ii = 2; ii < m_n; ++ii) {                 m_k_vector[ii] =                     m_a3 * m_qk_vector[ii - 2] - m_a7 * m_qp_vector[ii - 1] + m_qp_vector[ii];             }         }     }      return; }  //====================================================================== //  Compute new estimates of the quadratic coefficients using the //  scalars computed in CalcSc. //======================================================================  template<std::int32_t degree> void PolynomialRootFinder<degree>::Newest(std::int32_t itype, double &uu, double &vv) {     //------------------------------------------------------------------     //  Use formulas appropriate to setting of itype.     //------------------------------------------------------------------      if (itype == 3) {         //--------------------------------------------------------------         //  If itype == 3 the quadratic is zeroed.         //--------------------------------------------------------------          uu = 0.0;         vv = 0.0;     } else {         double a4;         double a5;          if (itype == 2) {             a4 = (m_a + m_g) * m_f + m_h;             a5 = (m_f + m_u) * m_c + m_v * m_d;         } else {             a4 = m_a + m_u * m_b + m_h * m_f;             a5 = m_c + (m_u + m_v * m_f) * m_d;         }          //--------------------------------------------------------------         //  Evaluate new quadratic coefficients.         //--------------------------------------------------------------          double b1 = -m_k_vector[m_n - 1] / m_p_vector[m_n];         double b2 = -(m_k_vector[m_n - 2] + b1 * m_p_vector[m_n - 1]) / m_p_vector[m_n];         double c1 = m_v * b2 * m_a1;         double c2 = b1 * m_a7;         double c3 = b1 * b1 * m_a3;         double c4 = c1 - c2 - c3;         double temp = a5 + b1 * a4 - c4;          if (temp != 0.0) {             uu = m_u - (m_u * (c3 + c2) + m_v * (b1 * m_a1 + b2 * m_a7)) / temp;             vv = m_v * (1.0 + c4 / temp);         }     }      return; }  //====================================================================== //  Divides p by the quadratic  1, u, v placing the quotient in q //  and the remainder in a,b //======================================================================  template<std::int32_t degree> void PolynomialRootFinder<degree>::QuadraticSyntheticDivision(std::int32_t n_plus_one,     double u,     double v,     double *p_ptr,     double *q_ptr,     double &a,     double &b) {     b = p_ptr[0];     q_ptr[0] = b;     a = p_ptr[1] - u * b;     q_ptr[1] = a;      for (std::int32_t ii = 2; ii < n_plus_one; ++ii) {         double c = p_ptr[ii] - u * a - v * b;         q_ptr[ii] = c;         b = a;         a = c;     }      return; }  //====================================================================== //                                          2 //  Calculate the zeros of the quadratic a x + b x + c. //  the quadratic formula, modified to avoid overflow, is used to find //  the larger zero if the zeros are real and both zeros are complex. //  the smaller real zero is found directly from the product of the //  zeros c / a. //======================================================================  template<std::int32_t degree> void PolynomialRootFinder<degree>::SolveQuadraticEquation(double a,     double b,     double c,     double &sr,     double &si,     double &lr,     double &li) {     if (a == 0.0) {         if (b != 0.0) {             sr = -c / b;         } else {             sr = 0.0;         }          lr = 0.0;         si = 0.0;         li = 0.0;     } else if (c == 0.0) {         sr = 0.0;         lr = -b / a;         si = 0.0;         li = 0.0;     } else {         //--------------------------------------------------------------         //  Compute discriminant avoiding overflow.         //--------------------------------------------------------------          double d;         double e;         double bvar = b / 2.0;          if (::fabs(bvar) < ::fabs(c)) {             if (c < 0.0) {                 e = -a;             } else {                 e = a;             }              e = bvar * (bvar / ::fabs(c)) - e;              d = ::sqrt(::fabs(e)) * ::sqrt(::fabs(c));         } else {             e = 1.0 - (a / bvar) * (c / bvar);             d = ::sqrt(::fabs(e)) * ::fabs(bvar);         }          if (e >= 0.0) {             //----------------------------------------------------------             //  Real zeros             //----------------------------------------------------------              if (bvar >= 0.0) {                 d = -d;             }              lr = (-bvar + d) / a;             sr = 0.0;              if (lr != 0.0) {                 sr = (c / lr) / a;             }              si = 0.0;             li = 0.0;         } else {             //----------------------------------------------------------             //  Complex conjugate zeros             //----------------------------------------------------------              sr = -bvar / a;             lr = sr;             si = ::fabs(d / a);             li = -si;         }     }      return; } #endif   

main.cc

  /* standard headers */ #include <cstdint> #include <array> #include <vector> #include <cmath> #include <random> #include <chrono> #include <cstring> #include <cstdio>  /* omp headers */ #include <omp.h>  /* png headers */ #include "png.hh"  /* note: I did not create the polynomial related files I got them from here  * and yes I know this has terrible code quailty but I could not find anything better   https://www.codeproject.com/articles/674149/a-real-polynomial-class-with-root-finder   */ #include "PolynomialRootFinder.hh"  /* constants */ namespace {     constexpr std::uint32_t width = 500;     constexpr std::uint32_t height = 500;     constexpr std::int32_t degree = 24;     constexpr std::int32_t coefficients = degree + 1;     constexpr std::uint64_t total_samples = 1000000;     std::uint64_t const individual_samples = total_samples / omp_get_max_threads();     using roots_t = std::array<double, coefficients * 2>;     using heatmap_t = std::uint32_t; }   std::int32_t generate_roots(roots_t &output) {     static thread_local std::mt19937_64 mt(std::random_device{}());     std::uniform_int_distribution<std::int32_t> dist(0, 1);      std::array<double, coefficients> cofs;     std::transform(cofs.begin(), cofs.end(), cofs.begin(), [&](auto) { return dist(mt) ? 1 : -1; });     PolynomialRootFinder<degree> poly = {};      std::int32_t roots_found;     if (poly.FindRoots(&cofs[0], &output[0], &output[coefficients], &roots_found) == PolynomialRootFinder<degree>::RootStatus_T::SUCCESS) {         return roots_found;     } else {         return 0;     } }  void generate_heatmap(std::vector<heatmap_t> &heatmap, heatmap_t &max_value) {     roots_t roots = {};     auto map_range = [](auto s, decltype(s) a1, decltype(s) a2, decltype(s) b1, decltype(s) b2) {         return b1 + (s - a1) * (b2 - b1) / (a2 - a1);     };      for (std::uint64_t i = 0; i < individual_samples; ++i) {         /* see if we found any roots */         if (std::int32_t roots_found = generate_roots(roots)) {             /* plot all the roots found to the heatmap */             while (--roots_found >= 0) {                 double const real = roots[roots_found];                 double const imag = roots[static_cast<std::size_t>(roots_found) + coefficients];                  std::int32_t const col = static_cast<std::int32_t>(map_range(real, -1.6, 1.6, 0, width));                 std::int32_t const row = static_cast<std::int32_t>(map_range(imag, -1.6, 1.6, 0, height));                 /* only plot roots that are in bounds */                 if (col < 0 || col >= width || row < 0 || row >= height) continue;                 max_value = std::max(++heatmap[static_cast<std::size_t>(row) * width + col], max_value);             }         }     } }  int main() {     /* create a heatmap*/     std::vector<heatmap_t> heatmap(width * height);      /* start a timer */     std::chrono::time_point<std::chrono::high_resolution_clock> const t1 = std::chrono::high_resolution_clock::now();      /* generate heatmap */     heatmap_t max_value = 0; #pragma omp parallel     generate_heatmap(heatmap, max_value);      /* write image */     std::vector<std::uint8_t> image;     image.resize(width * height);     for (std::int32_t i = 0; i < width * height; ++i) {             std::uint8_t color = static_cast<std::uint8_t>((std::log(heatmap[i]) / std::log(max_value)) * 255.0 + 0.55555);             image[i] = color;     }     png::write_image("output.png", image.data(), width, height);      /* print the time it took */     std::chrono::time_point<std::chrono::high_resolution_clock> const t2 = std::chrono::high_resolution_clock::now();     std::chrono::duration<double> const duration =         std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);     double const time_took = duration.count();     std::printf("It took %f %s", time_took, std::array{ "seconds", "second" } [1.0 == time_took]);      /* wait for user input to close */     (void)std::getchar(); }   

Aquí es lo que produce el programa:

ingrese la descripción de la imagen aquí

Original en ingles

a Littlewood polynomial is a polynomial where each coefficient is ether -1 or 1 and when the complex roots of them it produces a nice image. As a result I decided to make a program in c++ that plots a heatmap of the complex roots of randomly chosen Littlewood polynomials. Moreover at first I thought it would be as using GSL(gnu scientific library) as it already had a polynomial solver and libpng. However I after implementing a basic version of it using GSL I realized that GSL was to slow to use. Consequently that meant I had to find another polynomial library and after a little searching I found this https://www.codeproject.com/articles/674149/a-real-polynomial-class-with-root-finder. Then after cleaning up that library and learning how to use libpng it was mostly smooth sailing. However my main concerns are that I could be still be ways of improving the performance of the code that I don't know of and that I could still improve the code quality in ways that I don't know of.

png.hh

#ifndef PNG_HH #define PNG_HH #include <png.h>  namespace png {     void write_image(char const *filename, std::uint8_t const *image_data, std::uint32_t image_width, std::uint32_t image_height)     {         /* create a zeroed out png_image struct */         png_image output_png;         std::memset(&output_png, 0, sizeof(output_png));         output_png.version = PNG_IMAGE_VERSION;         output_png.format = PNG_FORMAT_GRAY;         output_png.width = image_width;         output_png.height = image_height;          /* write the png file */         png_image_write_to_file(&output_png, filename, 0, image_data, image_height, nullptr);          /* cleanup */         png_image_free(&output_png);     } } #endif 

PolynomialRootFinder.hh

//======================================================================= // Copyright (C) 2003-2013 William Hallahan // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without restriction, // including without limitation the rights to use, copy, modify, merge, // publish, distribute, sublicense, and/or sell copies of the Software, // and to permit persons to whom the Software is furnished to do so, // subject to the following conditions: // // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. // // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT // HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, // WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR // OTHER DEALINGS IN THE SOFTWARE. //=======================================================================  //********************************************************************** //  File: PolynomialRootFinder.h //  Author: Bill Hallahan //  Date: January 30, 2003 // //  Abstract: // //    This file contains the definition for class PolynomialRootFinder. // //**********************************************************************  #ifndef POLYNOMIALROOTFINDER_H #define POLYNOMIALROOTFINDER_H  #include <array>  //====================================================================== //  Class definition. //======================================================================   template<std::int32_t degree> struct PolynomialRootFinder {     std::array<double, degree + 1> m_p_vector;     std::array<double, degree + 1> m_qp_vector;     std::array<double, degree + 1> m_k_vector;     std::array<double, degree + 1> m_qk_vector;     std::array<double, degree + 1> m_svk_vector;      std::int32_t m_n;     std::int32_t m_n_plus_one;     double m_real_s;     double m_imag_s;     double m_u;     double m_v;     double m_a;     double m_b;     double m_c;     double m_d;     double m_a1;     double m_a2;     double m_a3;     double m_a6;     double m_a7;     double m_e;     double m_f;     double m_g;     double m_h;     double m_real_sz;     double m_imag_sz;     double m_real_lz;     double m_imag_lz;     double m_are;     double m_mre;      enum class RootStatus_T     {         SUCCESS,         LEADING_COEFFICIENT_IS_ZERO,         SCALAR_VALUE_HAS_NO_ROOTS,         FAILED_TO_CONVERGE     };      PolynomialRootFinder::RootStatus_T FindRoots(double *coefficient_ptr,         double *real_zero_vector_ptr,         double *imaginary_zero_vector_ptr,         std::int32_t *number_of_roots_found_ptr = 0);      std::int32_t Fxshfr(std::int32_t l2var);      std::int32_t QuadraticIteration(double uu, double vv);      std::int32_t RealIteration(double &sss, std::int32_t &flag);      std::int32_t CalcSc();      void NextK(std::int32_t itype);      void Newest(std::int32_t itype, double &uu, double &vv);      void QuadraticSyntheticDivision(std::int32_t n_plus_one,         double u,         double v,         double *p_ptr,         double *q_ptr,         double &a,         double &b);      void SolveQuadraticEquation(double a,         double b,         double c,         double &sr,         double &si,         double &lr,         double &li); };  #include <cmath> #include <float.h>  namespace {     //------------------------------------------------------------------     //  The following machine constants are used in this method.     //     //    f_BASE  The base of the floating postd::int32_t number system used.     //     //    f_ETA  The maximum relative representation error which     //           can be described as the smallest positive floating     //           postd::int32_t number such that 1.0 + f_ETA is greater than 1.0.     //     //    f_MAXIMUM_FLOAT  The largest floating postd::int32_t number.     //     //    f_MINIMUM_FLOAT  The smallest positive floating postd::int32_t number.     //     //------------------------------------------------------------------      constexpr float f_BASE = 2.0;     constexpr float f_ETA = FLT_EPSILON;     constexpr float f_ETA_N = (10.0f) * f_ETA;     constexpr float f_ETA_N_SQUARED = (100.0f) * f_ETA;     constexpr float f_MAXIMUM_FLOAT = FLT_MAX;     constexpr float f_MINIMUM_FLOAT = FLT_MIN;     constexpr float f_XX_INITIAL_VALUE = (0.70710678f);     constexpr float f_COSR_INITIAL_VALUE = (-0.069756474f);     constexpr float f_SINR_INITIAL_VALUE = (0.99756405f); };  //====================================================================== //  Member Function: PolynomialRootFinder::FindRoots // //  Abstract: // //    This method determines the roots of a polynomial which //    has real coefficients. This code is based on FORTRAN //    code published in reference [1]. The method is based on //    an algorithm the three-stage algorithm described in //    Jenkins and Traub [2]. // // 1. "Collected Algorithms from ACM, Volume III", Algorithms 493-545 //    1983. (The root finding algorithms is number 493) // // 2. Jenkins, M. A. and Traub, J. F., "A three-stage algorithm for //    real polynomials using quadratic iteration", SIAM Journal of //    Numerical Analysis, 7 (1970), 545-566 // // 3. Jenkins, M. A. and Traub, J. F., "Principles for testing //    polynomial zerofinding programs", ACM TOMS 1, //    1 (March 1975), 26-34 // // //  Input: // //    All vectors below must be at least a length equal to degree + 1. // //    coefficicent_ptr           A double precision vector that contains //                               the polynomial coefficients in order //                               of increasing power. // //    degree                     The degree of the polynomial. // //    real_zero_vector_ptr       A double precision vector that will //                               contain the real parts of the roots //                               of the polynomial when this method //                               returns. // //    imaginary_zero_vector_ptr  A double precision vector that will //                               contain the real parts of the roots //                               of the polynomial when this method //                               returns. // //    number_of_roots_found_ptr  A postd::int32_ter to an std::int32_teger that will //                               equal the number of roots found when //                               this method returns. If the method //                               returns SUCCESS then this value will //                               always equal the degree of the //                               polynomial. // //  Return Value: // //    The function returns an enum value of type //    'PolynomialRootFinder::RootStatus_T'. // //======================================================================  template<std::int32_t degree> typename PolynomialRootFinder<degree>::RootStatus_T PolynomialRootFinder<degree>::FindRoots(     double *coefficient_vector_ptr,     double *real_zero_vector_ptr,     double *imaginary_zero_vector_ptr,     std::int32_t *number_of_roots_found_ptr) {     //--------------------------------------------------------------     //  The algorithm fails if the polynomial is not at least     //  degree on or the leading coefficient is zero.     //--------------------------------------------------------------      PolynomialRootFinder::RootStatus_T status;     //--------------------------------------------------------------     //  Allocate temporary vectors used to find the roots.     //--------------------------------------------------------------       std::array<double, degree + 1> temp_vector;     std::array<double, degree + 1> pt_vector;      //--------------------------------------------------------------     //  m_are and m_mre refer to the unit error in + and *     //  respectively. they are assumed to be the same as     //  f_ETA.     //--------------------------------------------------------------      m_are = f_ETA;     m_mre = f_ETA;     double lo = f_MINIMUM_FLOAT / f_ETA;      //--------------------------------------------------------------     // Initialization of constants for shift rotation.     //--------------------------------------------------------------      double xx = f_XX_INITIAL_VALUE;     double yy = -xx;     double cosr = f_COSR_INITIAL_VALUE;     double sinr = f_SINR_INITIAL_VALUE;     m_n = degree;     m_n_plus_one = m_n + 1;      //--------------------------------------------------------------     //  Make a copy of the coefficients in reverse order.     //--------------------------------------------------------------      std::int32_t ii = 0;      for (ii = 0; ii < m_n_plus_one; ++ii) {         m_p_vector[m_n - ii] = coefficient_vector_ptr[ii];     }      //--------------------------------------------------------------     //  Assume failure. The status is set to SUCCESS if all     //  the roots are found.     //--------------------------------------------------------------      status = PolynomialRootFinder::RootStatus_T::FAILED_TO_CONVERGE;      //--------------------------------------------------------------     //  If there are any zeros at the origin, remove them.     //--------------------------------------------------------------      std::int32_t jvar = 0;      while (m_p_vector[m_n] == 0.0) {         jvar = degree - m_n;         real_zero_vector_ptr[jvar] = 0.0;         imaginary_zero_vector_ptr[jvar] = 0.0;         m_n_plus_one = m_n_plus_one - 1;         m_n = m_n - 1;     }      //--------------------------------------------------------------     //  Loop and find polynomial zeros. In the original algorithm     //  this loop was an infinite loop. Testing revealed that the     //  number of main loop iterations to solve a polynomial of a     //  particular degree is usually about half the degree.     //  We loop twice that to make sure the solution is found.     //  (This should be revisited as it might preclude solving     //  some large polynomials.)     //--------------------------------------------------------------       for (std::int32_t count = 0; count < degree; ++count) {         //----------------------------------------------------------         //  Check for less than 2 zeros to finish the solution.         //----------------------------------------------------------          if (m_n <= 2) {             if (m_n > 0) {                 //--------------------------------------------------                 //  Calculate the final zero or pair of zeros.                 //--------------------------------------------------                  if (m_n == 1) {                     real_zero_vector_ptr[degree - 1] =                         -m_p_vector[1] / m_p_vector[0];                      imaginary_zero_vector_ptr[degree - 1] = 0.0;                 } else {                     SolveQuadraticEquation(                         m_p_vector[0],                         m_p_vector[1],                         m_p_vector[2],                         real_zero_vector_ptr[degree - 2],                         imaginary_zero_vector_ptr[degree - 2],                         real_zero_vector_ptr[degree - 1],                         imaginary_zero_vector_ptr[degree - 1]);                 }             }              m_n = 0;             status = PolynomialRootFinder::RootStatus_T::SUCCESS;             break;         } else {             //------------------------------------------------------             //  Find largest and smallest moduli of coefficients.             //------------------------------------------------------              double max = 0.0;             double min = f_MAXIMUM_FLOAT;             double xvar;              for (ii = 0; ii < m_n_plus_one; ++ii) {                 xvar = (double)(::fabs((double)(m_p_vector[ii])));                  if (xvar > max) {                     max = xvar;                 }                  if ((xvar != 0.0) && (xvar < min)) {                     min = xvar;                 }             }              //------------------------------------------------------             //  Scale if there are large or very small coefficients.             //  Computes a scale factor to multiply the coefficients             //  of the polynomial. The scaling is done to avoid             //  overflow and to avoid undetected underflow from             //  std::int32_terfering with the convergence criterion.             //  The factor is a power of the base.             //------------------------------------------------------              bool do_scaling_flag = false;             double sc = lo / min;              if (sc <= 1.0) {                 do_scaling_flag = f_MAXIMUM_FLOAT / sc < max;             } else {                 do_scaling_flag = max < 10.0;                  if (!do_scaling_flag) {                     if (sc == 0.0) {                         sc = f_MINIMUM_FLOAT;                     }                 }             }              //------------------------------------------------------             //  Conditionally scale the data.             //------------------------------------------------------              if (do_scaling_flag) {                 std::int32_t lvar = (std::int32_t)(::log(sc) / ::log(f_BASE) + 0.5);                 double factor = ::pow((double)(f_BASE * 1.0), double(lvar));                  if (factor != 1.0) {                     for (ii = 0; ii < m_n_plus_one; ++ii) {                         m_p_vector[ii] = factor * m_p_vector[ii];                     }                 }             }              //------------------------------------------------------             //  Compute lower bound on moduli of zeros.             //------------------------------------------------------              for (ii = 0; ii < m_n_plus_one; ++ii) {                 pt_vector[ii] = (double)(::fabs((double)(m_p_vector[ii])));             }              pt_vector[m_n] = -pt_vector[m_n];              //------------------------------------------------------             //  Compute upper estimate of bound.             //------------------------------------------------------              xvar = (double)                 (::exp((::log(-pt_vector[m_n]) - ::log(pt_vector[0]))                     / (double)(m_n)));              //------------------------------------------------------             //  If newton step at the origin is better, use it.             //------------------------------------------------------              double xm;              if (pt_vector[m_n - 1] != 0.0) {                 xm = -pt_vector[m_n] / pt_vector[m_n - 1];                  if (xm < xvar) {                     xvar = xm;                 }             }              //------------------------------------------------------             //  Chop the std::int32_terval (0, xvar) until ff <= 0             //------------------------------------------------------              double ff;              for (;;) {                 xm = (double)(xvar * 0.1);                 ff = pt_vector[0];                  for (ii = 1; ii < m_n_plus_one; ++ii) {                     ff = ff * xm + pt_vector[ii];                 }                 if (ff <= 0.0) {                     ff = 0.0;                     break;                 }                  xvar = xm;             }              double dx = xvar;              //------------------------------------------------------             //  Do newton iteration until xvar converges to two             //  decimal places.             //------------------------------------------------------              for (;;) {                 if ((double)(::fabs(dx / xvar)) <= 0.005) {                     break;                 }                  ff = pt_vector[0];                  double df = ff;                  for (ii = 1; ii < m_n; ++ii) {                     ff = ff * xvar + pt_vector[ii];                      df = df * xvar + ff;                  }                  ff = ff * xvar + pt_vector[m_n];                  dx = ff / df;                   xvar = xvar - dx;             }              double bnd = xvar;              //------------------------------------------------------             //  Compute the derivative as the std::int32_tial m_k_vector             //  polynomial and do 5 steps with no shift.             //------------------------------------------------------              std::int32_t n_minus_one = m_n - 1;              for (ii = 1; ii < m_n; ++ii) {                 m_k_vector[ii] =                     (double)(m_n - ii) * m_p_vector[ii] / (double)(m_n);             }              m_k_vector[0] = m_p_vector[0];             double aa = m_p_vector[m_n];             double bb = m_p_vector[m_n - 1];             bool zerok_flag = m_k_vector[m_n - 1] == 0.0;              std::int32_t jj = 0;              for (jj = 1; jj <= 5; ++jj) {                 double cc = m_k_vector[m_n - 1];                  if (zerok_flag) {                     //----------------------------------------------                     //  Use unscaled form of recurrence.                     //----------------------------------------------                      for (jvar = n_minus_one; jvar > 0; --jvar) {                         m_k_vector[jvar] = m_k_vector[jvar - 1];                     }                      m_k_vector[0] = 0.0;                     zerok_flag = m_k_vector[m_n - 1] == 0.0;                 } else {                     //----------------------------------------------                     //  Use scaled form of recurrence if value                     //  of m_k_vector at 0 is nonzero.                     //----------------------------------------------                      double tvar = -aa / cc;                      for (jvar = n_minus_one; jvar > 0; --jvar) {                         m_k_vector[jvar] =                             tvar * m_k_vector[jvar - 1] + m_p_vector[jvar];                     }                      m_k_vector[0] = m_p_vector[0];                     zerok_flag =                         ::fabs(m_k_vector[m_n - 1]) <= ::fabs(bb) * f_ETA_N;                 }             }              //------------------------------------------------------             //  Save m_k_vector for restarts with new shifts.             //------------------------------------------------------              for (ii = 0; ii < m_n; ++ii) {                 temp_vector[ii] = m_k_vector[ii];             }              //------------------------------------------------------             //  Loop to select the quadratic corresponding to             //   each new shift.             //------------------------------------------------------              std::int32_t cnt = 0;              for (cnt = 1; cnt <= 20; ++cnt) {                 //--------------------------------------------------                 //  Quadratic corresponds to a double shift to a                 //  non-real postd::int32_t and its complex conjugate. The                 //  postd::int32_t has modulus 'bnd' and amplitude rotated                 //  by 94 degrees from the previous shift.                 //--------------------------------------------------                  double xxx = cosr * xx - sinr * yy;                 yy = sinr * xx + cosr * yy;                 xx = xxx;                 m_real_s = bnd * xx;                 m_imag_s = bnd * yy;                 m_u = -2.0 * m_real_s;                 m_v = bnd;                  //--------------------------------------------------                 //  Second stage calculation, fixed quadratic.                 //  Variable nz will contain the number of                 //   zeros found when function Fxshfr() returns.                 //--------------------------------------------------                  std::int32_t nz = Fxshfr(20 * cnt);                  if (nz != 0) {                     //----------------------------------------------                     //  The second stage jumps directly to one of                     //  the third stage iterations and returns here                     //  if successful. Deflate the polynomial,                     //  store the zero or zeros and return to the                     //  main algorithm.                     //----------------------------------------------                      jvar = degree - m_n;                     real_zero_vector_ptr[jvar] = m_real_sz;                     imaginary_zero_vector_ptr[jvar] = m_imag_sz;                     m_n_plus_one = m_n_plus_one - nz;                     m_n = m_n_plus_one - 1;                      for (ii = 0; ii < m_n_plus_one; ++ii) {                         m_p_vector[ii] = m_qp_vector[ii];                     }                      if (nz != 1) {                         real_zero_vector_ptr[jvar + 1] = m_real_lz;                         imaginary_zero_vector_ptr[jvar + 1] = m_imag_lz;                     }                      break;                      //----------------------------------------------                     //  If the iteration is unsuccessful another                     //  quadratic is chosen after restoring                     //  m_k_vector.                     //----------------------------------------------                 }                  for (ii = 0; ii < m_n; ++ii) {                     m_k_vector[ii] = temp_vector[ii];                 }             }         }     }      //--------------------------------------------------------------     //  If no convergence with 20 shifts then adjust the degree     //  for the number of roots found.     //--------------------------------------------------------------      if (number_of_roots_found_ptr != 0) {         *number_of_roots_found_ptr = degree - m_n;     }      return status; }  //====================================================================== //  Computes up to l2var fixed shift m_k_vector polynomials, //  testing for convergence in the linear or quadratic //  case. initiates one of the variable shift //  iterations and returns with the number of zeros //  found. // //    l2var  An std::int32_teger that is the limit of fixed shift steps. // //  Return Value: //    nz   An std::int32_teger that is the number of zeros found. //======================================================================  template<std::int32_t degree> std::int32_t PolynomialRootFinder<degree>::Fxshfr(std::int32_t l2var) {     //------------------------------------------------------------------     //  Evaluate polynomial by synthetic division.     //------------------------------------------------------------------      QuadraticSyntheticDivision(m_n_plus_one,         m_u,         m_v,         m_p_vector.data(),         m_qp_vector.data(),         m_a,         m_b);      std::int32_t itype = CalcSc();      std::int32_t nz = 0;     float betav = 0.25;     float betas = 0.25;     float oss = (float)(m_real_s);     float ovv = (float)(m_v);     float ots;     float otv;     double ui = 0.0;     double vi = 0.0;     double svar;       for (std::int32_t jvar = 1; jvar <= l2var; ++jvar) {         //--------------------------------------------------------------         //  Calculate next m_k_vector polynomial and estimate m_v.         //--------------------------------------------------------------          NextK(itype);         itype = CalcSc();         Newest(itype, ui, vi);         float vv = (float)(vi);          //--------------------------------------------------------------         //  Estimate svar         //--------------------------------------------------------------          float ss = 0.0;          if (m_k_vector[m_n - 1] != 0.0) {             ss = (float)(-m_p_vector[m_n] / m_k_vector[m_n - 1]);         }          float tv = 1.0;         float ts = 1.0;          if ((jvar != 1) && (itype != 3)) {             //----------------------------------------------------------             //  Compute relative measures of convergence of             //  svar and m_v sequences.             //----------------------------------------------------------              if (vv != 0.0) {                 tv = (float)(::fabs((vv - ovv) / vv));             }              if (ss != 0.0) {                 ts = (float)(::fabs((ss - oss) / ss));             }              //----------------------------------------------------------             //  If decreasing, multiply two most recent convergence             //  measures.             //----------------------------------------------------------              float tvv = 1.0;              if (tv < otv) {                 tvv = tv * otv;             }              float tss = 1.0;              if (ts < ots) {                 tss = ts * ots;             }              //----------------------------------------------------------             //  Compare with convergence criteria.             //----------------------------------------------------------              bool vpass_flag = tvv < betav;             bool spass_flag = tss < betas;              if (spass_flag || vpass_flag) {                 //------------------------------------------------------                 //  At least one sequence has passed the convergence                 //  test. Store variables before iterating.                 //------------------------------------------------------                  double svu = m_u;                 double svv = m_v;                 std::int32_t ii = 0;                  for (ii = 0; ii < m_n; ++ii) {                     m_svk_vector[ii] = m_k_vector[ii];                 }                  svar = ss;                  //------------------------------------------------------                 //  Choose iteration according to the fastest                 //  converging sequence.                 //------------------------------------------------------                  bool vtry_flag = false;                 bool stry_flag = false;                 bool exit_outer_loop_flag = false;                  bool start_with_real_iteration_flag =                     (spass_flag && ((!vpass_flag) || (tss < tvv)));                  do {                     if (!start_with_real_iteration_flag) {                         nz = QuadraticIteration(ui, vi);                          if (nz > 0) {                             exit_outer_loop_flag = true;                             break;                         }                          //----------------------------------------------                         //  Quadratic iteration has failed. flag                         //  that it has been tried and decrease                         //  the convergence criterion.                         //----------------------------------------------                          vtry_flag = true;                         betav = (float)(betav * 0.25);                     }                      //--------------------------------------------------                     //  Try linear iteration if it has not been                     //  tried and the svar sequence is converging.                     //--------------------------------------------------                      if (((!stry_flag) && spass_flag)                         || start_with_real_iteration_flag) {                         if (!start_with_real_iteration_flag) {                             for (ii = 0; ii < m_n; ++ii) {                                 m_k_vector[ii] = m_svk_vector[ii];                             }                         } else {                             start_with_real_iteration_flag = false;                         }                          std::int32_t iflag = 0;                          nz = RealIteration(svar, iflag);                          if (nz > 0) {                             exit_outer_loop_flag = true;                             break;                         }                          //----------------------------------------------                         //  Linear iteration has failed. Flag that                         //  it has been tried and decrease the                         //  convergence criterion.                         //----------------------------------------------                          stry_flag = true;                         betas = (float)(betas * 0.25);                          if (iflag != 0) {                             //------------------------------------------                             //  If linear iteration signals an almost                             //  double real zero attempt quadratic                             //  iteration.                             //------------------------------------------                              ui = -(svar + svar);                             vi = svar * svar;                              continue;                         }                     }                      //--------------------------------------------------                     //  Restore variables                     //--------------------------------------------------                      m_u = svu;                     m_v = svv;                      for (ii = 0; ii < m_n; ++ii) {                         m_k_vector[ii] = m_svk_vector[ii];                     }                      //----------------------------------------------                     //  Try quadratic iteration if it has not been                     //  tried and the m_v sequence is converging.                     //----------------------------------------------                 } while (vpass_flag && (!vtry_flag));                  if (exit_outer_loop_flag) {                     break;                 }                  //------------------------------------------------------                 //  Recompute m_qp_vector and scalar values to                 //  continue the second stage.                 //------------------------------------------------------                  QuadraticSyntheticDivision(m_n_plus_one,                     m_u,                     m_v,                     m_p_vector.data(),                     m_qp_vector.data(),                     m_a,                     m_b);                  itype = CalcSc();             }         }          ovv = vv;         oss = ss;         otv = tv;         ots = ts;     }      return nz; }  //====================================================================== //  Variable-shift m_k_vector-polynomial iteration for //  a quadratic factor converges only if the zeros are //  equimodular or nearly so. // //    uu  Coefficients of starting quadratic //    vv  Coefficients of starting quadratic // //  Return value: //    nz  The number of zeros found. //======================================================================  template<std::int32_t degree> std::int32_t PolynomialRootFinder<degree>::QuadraticIteration(double uu, double vv) {     //------------------------------------------------------------------     //  Main loop     //------------------------------------------------------------------      double ui = 0.0;     double vi = 0.0;     float omp = 0.0F;     float relstp = 0.0F;     std::int32_t itype = 0;     bool tried_flag = false;     std::int32_t jvar = 0;     std::int32_t nz = 0;     m_u = uu;     m_v = vv;      for(;;) {         SolveQuadraticEquation(1.0,             m_u,             m_v,             m_real_sz,             m_imag_sz,             m_real_lz,             m_imag_lz);          //--------------------------------------------------------------         //  Return if roots of the quadratic are real and not close         //  to multiple or nearly equal and  of opposite sign.         //--------------------------------------------------------------          if (::fabs(::fabs(m_real_sz) - ::fabs(m_real_lz)) > 0.01 * ::fabs(m_real_lz)) {             break;         }          //--------------------------------------------------------------         //  Evaluate polynomial by quadratic synthetic division.         //------------------------------------------------------------------          QuadraticSyntheticDivision(m_n_plus_one,             m_u,             m_v,             m_p_vector.data(),             m_qp_vector.data(),             m_a,             m_b);          float mp = (float)(::fabs(m_a - m_real_sz * m_b) + ::fabs(m_imag_sz * m_b));          //--------------------------------------------------------------         //  Compute a rigorous  bound on the rounding error in         //  evaluting m_p_vector.         //--------------------------------------------------------------          float zm = (float)(::sqrt((float)(::fabs((float)(m_v)))));         float ee = (float)(2.0 * (float)(::fabs((float)(m_qp_vector[0]))));         float tvar = (float)(-m_real_sz * m_b);         std::int32_t ii = 0;          for (ii = 1; ii < m_n; ++ii) {             ee = ee * zm + (float)(::fabs((float)(m_qp_vector[ii])));         }          ee = ee * zm + (float)(::fabs((float)(m_a)+tvar));         ee = (float)((5.0 * m_mre + 4.0 * m_are) * ee             - (5.0 * m_mre + 2.0 * m_are) * ((float)(::fabs((float)(m_a)+tvar)) + (float)(::fabs((float)(m_b))) * zm)             + 2.0 * m_are * (float)(::fabs(tvar)));          //--------------------------------------------------------------         //  Iteration has converged sufficiently if the polynomial         //  value is less than 20 times this bound.         //--------------------------------------------------------------          if (mp <= 20.0 * ee) {             nz = 2;             break;         }          jvar = jvar + 1;          //--------------------------------------------------------------         //  Stop iteration after 20 steps.         //--------------------------------------------------------------          if (jvar > 20) {             break;         }          if ((jvar >= 2) && ((relstp <= 0.01)             && (mp >= omp) && (!tried_flag))) {             //----------------------------------------------------------             //  A cluster appears to be stalling the convergence.             //  Five fixed shift steps are taken with a m_u, m_v             //  close to the cluster.             //----------------------------------------------------------              if (relstp < f_ETA) {                 relstp = f_ETA;             }              relstp = (float)(::sqrt(relstp));             m_u = m_u - m_u * relstp;             m_v = m_v + m_v * relstp;              QuadraticSyntheticDivision(m_n_plus_one,                 m_u,                 m_v,                 m_p_vector.data(),                 m_qp_vector.data(),                 m_a,                 m_b);              for (ii = 0; ii < 5; ++ii) {                 itype = CalcSc();                 NextK(itype);             }              tried_flag = true;             jvar = 0;         }          omp = mp;          //--------------------------------------------------------------         //  Calculate next m_k_vector polynomial and         //  new m_u and m_v.         //--------------------------------------------------------------          itype = CalcSc();         NextK(itype);         itype = CalcSc();         Newest(itype, ui, vi);          //--------------------------------------------------------------         //  If vi is zero the iteration is not converging.         //--------------------------------------------------------------          if (vi == 0.0) {             break;         }          relstp = (float)(::fabs((vi - m_v) / vi));         m_u = ui;         m_v = vi;     }      return nz; }  //====================================================================== //  Variable-shift h polynomial iteration for a real zero. // //    sss      Starting iterate //    flag     Flag to indicate a pair of zeros near real axis. // //  Return Value: //     Number of zero found. //======================================================================  template<std::int32_t degree> std::int32_t PolynomialRootFinder<degree>::RealIteration(double &sss, std::int32_t &flag) {     //------------------------------------------------------------------     //  Main loop     //------------------------------------------------------------------      double tvar = 0.0;     float omp = 0.0F;     std::int32_t nz = 0;     flag = 0;     std::int32_t jvar = 0;     double svar = sss;      for(;;) {         double pv = m_p_vector[0];          //--------------------------------------------------------------         //  Evaluate m_p_vector at svar         //--------------------------------------------------------------          m_qp_vector[0] = pv;         std::int32_t ii = 0;          for (ii = 1; ii < m_n_plus_one; ++ii) {             pv = pv * svar + m_p_vector[ii];             m_qp_vector[ii] = pv;         }          float mp = (float)(::fabs(pv));          //--------------------------------------------------------------         //  Compute a rigorous bound on the error in evaluating p         //--------------------------------------------------------------          double ms = (double)(::fabs(svar));         double ee = (m_mre / (m_are + m_mre)) * (double)(::fabs((double)(m_qp_vector[0])));          for (ii = 1; ii < m_n_plus_one; ++ii) {             ee = ee * ms + (float)(::fabs((double)(m_qp_vector[ii])));         }          //--------------------------------------------------------------         //  Iteration has converged sufficiently if the         //  polynomial value is less than 20 times this bound.         //--------------------------------------------------------------          if (mp <= 20.0 * ((m_are + m_mre) * ee - m_mre * mp)) {             nz = 1;             m_real_sz = svar;             m_imag_sz = 0.0;             break;         }          jvar = jvar + 1;          //--------------------------------------------------------------         //  Stop iteration after 10 steps.         //--------------------------------------------------------------          if (jvar > 10) {             break;         }          if ((jvar >= 2)             && ((::fabs(tvar) <= 0.001 * ::fabs(svar - tvar))                 && (mp > omp))) {             //----------------------------------------------------------             //  A cluster of zeros near the real axis has been             //  encountered. Return with flag set to initiate             //  a quadratic iteration.             //----------------------------------------------------------              flag = 1;             sss = svar;             break;         }          //--------------------------------------------------------------         //  Return if the polynomial value has increased significantly.         //--------------------------------------------------------------          omp = mp;          //--------------------------------------------------------------         //  Compute t, the next polynomial, and the new iterate.         //--------------------------------------------------------------          double kv = m_k_vector[0];         m_qk_vector[0] = kv;          for (ii = 1; ii < m_n; ++ii) {             kv = kv * svar + m_k_vector[ii];             m_qk_vector[ii] = kv;         }          if (::fabs(kv) <= ::fabs(m_k_vector[m_n - 1]) * f_ETA_N) {             m_k_vector[0] = 0.0;              for (ii = 1; ii < m_n; ++ii) {                 m_k_vector[ii] = m_qk_vector[ii - 1];             }         } else {             //----------------------------------------------------------             //  Use the scaled form of the recurrence if the             //  value of m_k_vector at svar is non-zero.             //----------------------------------------------------------              tvar = -pv / kv;             m_k_vector[0] = m_qp_vector[0];              for (ii = 1; ii < m_n; ++ii) {                 m_k_vector[ii] = tvar * m_qk_vector[ii - 1] + m_qp_vector[ii];             }         }          //--------------------------------------------------------------         //  Use unscaled form.         //--------------------------------------------------------------          kv = m_k_vector[0];          for (ii = 1; ii < m_n; ++ii) {             kv = kv * svar + m_k_vector[ii];         }          tvar = 0.0;          if (::fabs(kv) > ::fabs(m_k_vector[m_n - 1]) * f_ETA_N) {             tvar = -pv / kv;         }          svar = svar + tvar;     }      return nz; }  //====================================================================== //  This routine calculates scalar quantities used to compute //  the next m_k_vector polynomial and new estimates of the //  quadratic coefficients. // //  Return Value: //    type  std::int32_teger variable set here indicating how the //    calculations are normalized to avoid overflow. //======================================================================  template<std::int32_t degree> std::int32_t PolynomialRootFinder<degree>::CalcSc() {     //------------------------------------------------------------------     //  Synthetic division of m_k_vector by the quadratic 1, m_u, m_v.     //------------------------------------------------------------------      QuadraticSyntheticDivision(m_n,         m_u,         m_v,         m_k_vector.data(),         m_qk_vector.data(),         m_c,         m_d);      std::int32_t itype = 0;      if ((::fabs(m_c) <= ::fabs(m_k_vector[m_n - 1]) * f_ETA_N_SQUARED)         && (::fabs(m_d) <= ::fabs(m_k_vector[m_n - 2]) * f_ETA_N_SQUARED)) {         //--------------------------------------------------------------         //  itype == 3 Indicates the quadratic is almost a         //  factor of m_k_vector.         //--------------------------------------------------------------          itype = 3;     } else if (::fabs(m_d) >= ::fabs(m_c)) {         //--------------------------------------------------------------         //  itype == 2 Indicates that all formulas are divided by m_d.         //--------------------------------------------------------------          itype = 2;         m_e = m_a / m_d;         m_f = m_c / m_d;         m_g = m_u * m_b;         m_h = m_v * m_b;         m_a3 = (m_a + m_g) * m_e + m_h * (m_b / m_d);         m_a1 = m_b * m_f - m_a;         m_a7 = (m_f + m_u) * m_a + m_h;     } else {         //--------------------------------------------------------------         //  itype == 1 Indicates that all formulas are divided by m_c.         //--------------------------------------------------------------          itype = 1;         m_e = m_a / m_c;         m_f = m_d / m_c;         m_g = m_u * m_e;         m_h = m_v * m_b;         m_a3 = m_a * m_e + (m_h / m_c + m_g) * m_b;         m_a1 = m_b - m_a * (m_d / m_c);         m_a7 = m_a + m_g * m_d + m_h * m_f;     }      return itype;  }  //====================================================================== //  Computes the next k polynomials using scalars computed in CalcSc. //====================================================================== template<std::int32_t degree> void PolynomialRootFinder<degree>::NextK(std::int32_t itype) {     std::int32_t ii = 0;      if (itype == 3) {         //--------------------------------------------------------------         //  Use unscaled form of the recurrence if type is 3.         //--------------------------------------------------------------          m_k_vector[0] = 0.0;         m_k_vector[1] = 0.0;          for (ii = 2; ii < m_n; ++ii) {             m_k_vector[ii] = m_qk_vector[ii - 2];         }     } else {         double temp = m_a;          if (itype == 1) {             temp = m_b;         }          if (::fabs(m_a1) <= ::fabs(temp) * f_ETA_N) {             //----------------------------------------------------------             //  If m_a1 is nearly zero then use a special form of             //  the recurrence.             //----------------------------------------------------------              m_k_vector[0] = 0.0;             m_k_vector[1] = -m_a7 * m_qp_vector[0];              for (ii = 2; ii < m_n; ++ii) {                 m_k_vector[ii] = m_a3 * m_qk_vector[ii - 2] - m_a7 * m_qp_vector[ii - 1];             }         } else {             //----------------------------------------------------------             //  Use scaled form of the recurrence.             //----------------------------------------------------------              m_a7 = m_a7 / m_a1;             m_a3 = m_a3 / m_a1;             m_k_vector[0] = m_qp_vector[0];             m_k_vector[1] = m_qp_vector[1] - m_a7 * m_qp_vector[0];              for (ii = 2; ii < m_n; ++ii) {                 m_k_vector[ii] =                     m_a3 * m_qk_vector[ii - 2] - m_a7 * m_qp_vector[ii - 1] + m_qp_vector[ii];             }         }     }      return; }  //====================================================================== //  Compute new estimates of the quadratic coefficients using the //  scalars computed in CalcSc. //======================================================================  template<std::int32_t degree> void PolynomialRootFinder<degree>::Newest(std::int32_t itype, double &uu, double &vv) {     //------------------------------------------------------------------     //  Use formulas appropriate to setting of itype.     //------------------------------------------------------------------      if (itype == 3) {         //--------------------------------------------------------------         //  If itype == 3 the quadratic is zeroed.         //--------------------------------------------------------------          uu = 0.0;         vv = 0.0;     } else {         double a4;         double a5;          if (itype == 2) {             a4 = (m_a + m_g) * m_f + m_h;             a5 = (m_f + m_u) * m_c + m_v * m_d;         } else {             a4 = m_a + m_u * m_b + m_h * m_f;             a5 = m_c + (m_u + m_v * m_f) * m_d;         }          //--------------------------------------------------------------         //  Evaluate new quadratic coefficients.         //--------------------------------------------------------------          double b1 = -m_k_vector[m_n - 1] / m_p_vector[m_n];         double b2 = -(m_k_vector[m_n - 2] + b1 * m_p_vector[m_n - 1]) / m_p_vector[m_n];         double c1 = m_v * b2 * m_a1;         double c2 = b1 * m_a7;         double c3 = b1 * b1 * m_a3;         double c4 = c1 - c2 - c3;         double temp = a5 + b1 * a4 - c4;          if (temp != 0.0) {             uu = m_u - (m_u * (c3 + c2) + m_v * (b1 * m_a1 + b2 * m_a7)) / temp;             vv = m_v * (1.0 + c4 / temp);         }     }      return; }  //====================================================================== //  Divides p by the quadratic  1, u, v placing the quotient in q //  and the remainder in a,b //======================================================================  template<std::int32_t degree> void PolynomialRootFinder<degree>::QuadraticSyntheticDivision(std::int32_t n_plus_one,     double u,     double v,     double *p_ptr,     double *q_ptr,     double &a,     double &b) {     b = p_ptr[0];     q_ptr[0] = b;     a = p_ptr[1] - u * b;     q_ptr[1] = a;      for (std::int32_t ii = 2; ii < n_plus_one; ++ii) {         double c = p_ptr[ii] - u * a - v * b;         q_ptr[ii] = c;         b = a;         a = c;     }      return; }  //====================================================================== //                                          2 //  Calculate the zeros of the quadratic a x + b x + c. //  the quadratic formula, modified to avoid overflow, is used to find //  the larger zero if the zeros are real and both zeros are complex. //  the smaller real zero is found directly from the product of the //  zeros c / a. //======================================================================  template<std::int32_t degree> void PolynomialRootFinder<degree>::SolveQuadraticEquation(double a,     double b,     double c,     double &sr,     double &si,     double &lr,     double &li) {     if (a == 0.0) {         if (b != 0.0) {             sr = -c / b;         } else {             sr = 0.0;         }          lr = 0.0;         si = 0.0;         li = 0.0;     } else if (c == 0.0) {         sr = 0.0;         lr = -b / a;         si = 0.0;         li = 0.0;     } else {         //--------------------------------------------------------------         //  Compute discriminant avoiding overflow.         //--------------------------------------------------------------          double d;         double e;         double bvar = b / 2.0;          if (::fabs(bvar) < ::fabs(c)) {             if (c < 0.0) {                 e = -a;             } else {                 e = a;             }              e = bvar * (bvar / ::fabs(c)) - e;              d = ::sqrt(::fabs(e)) * ::sqrt(::fabs(c));         } else {             e = 1.0 - (a / bvar) * (c / bvar);             d = ::sqrt(::fabs(e)) * ::fabs(bvar);         }          if (e >= 0.0) {             //----------------------------------------------------------             //  Real zeros             //----------------------------------------------------------              if (bvar >= 0.0) {                 d = -d;             }              lr = (-bvar + d) / a;             sr = 0.0;              if (lr != 0.0) {                 sr = (c / lr) / a;             }              si = 0.0;             li = 0.0;         } else {             //----------------------------------------------------------             //  Complex conjugate zeros             //----------------------------------------------------------              sr = -bvar / a;             lr = sr;             si = ::fabs(d / a);             li = -si;         }     }      return; } #endif 

main.cc

/* standard headers */ #include <cstdint> #include <array> #include <vector> #include <cmath> #include <random> #include <chrono> #include <cstring> #include <cstdio>  /* omp headers */ #include <omp.h>  /* png headers */ #include "png.hh"  /* note: I did not create the polynomial related files I got them from here  * and yes I know this has terrible code quailty but I could not find anything better   https://www.codeproject.com/articles/674149/a-real-polynomial-class-with-root-finder   */ #include "PolynomialRootFinder.hh"  /* constants */ namespace {     constexpr std::uint32_t width = 500;     constexpr std::uint32_t height = 500;     constexpr std::int32_t degree = 24;     constexpr std::int32_t coefficients = degree + 1;     constexpr std::uint64_t total_samples = 1000000;     std::uint64_t const individual_samples = total_samples / omp_get_max_threads();     using roots_t = std::array<double, coefficients * 2>;     using heatmap_t = std::uint32_t; }   std::int32_t generate_roots(roots_t &output) {     static thread_local std::mt19937_64 mt(std::random_device{}());     std::uniform_int_distribution<std::int32_t> dist(0, 1);      std::array<double, coefficients> cofs;     std::transform(cofs.begin(), cofs.end(), cofs.begin(), [&](auto) { return dist(mt) ? 1 : -1; });     PolynomialRootFinder<degree> poly = {};      std::int32_t roots_found;     if (poly.FindRoots(&cofs[0], &output[0], &output[coefficients], &roots_found) == PolynomialRootFinder<degree>::RootStatus_T::SUCCESS) {         return roots_found;     } else {         return 0;     } }  void generate_heatmap(std::vector<heatmap_t> &heatmap, heatmap_t &max_value) {     roots_t roots = {};     auto map_range = [](auto s, decltype(s) a1, decltype(s) a2, decltype(s) b1, decltype(s) b2) {         return b1 + (s - a1) * (b2 - b1) / (a2 - a1);     };      for (std::uint64_t i = 0; i < individual_samples; ++i) {         /* see if we found any roots */         if (std::int32_t roots_found = generate_roots(roots)) {             /* plot all the roots found to the heatmap */             while (--roots_found >= 0) {                 double const real = roots[roots_found];                 double const imag = roots[static_cast<std::size_t>(roots_found) + coefficients];                  std::int32_t const col = static_cast<std::int32_t>(map_range(real, -1.6, 1.6, 0, width));                 std::int32_t const row = static_cast<std::int32_t>(map_range(imag, -1.6, 1.6, 0, height));                 /* only plot roots that are in bounds */                 if (col < 0 || col >= width || row < 0 || row >= height) continue;                 max_value = std::max(++heatmap[static_cast<std::size_t>(row) * width + col], max_value);             }         }     } }  int main() {     /* create a heatmap*/     std::vector<heatmap_t> heatmap(width * height);      /* start a timer */     std::chrono::time_point<std::chrono::high_resolution_clock> const t1 = std::chrono::high_resolution_clock::now();      /* generate heatmap */     heatmap_t max_value = 0; #pragma omp parallel     generate_heatmap(heatmap, max_value);      /* write image */     std::vector<std::uint8_t> image;     image.resize(width * height);     for (std::int32_t i = 0; i < width * height; ++i) {             std::uint8_t color = static_cast<std::uint8_t>((std::log(heatmap[i]) / std::log(max_value)) * 255.0 + 0.55555);             image[i] = color;     }     png::write_image("output.png", image.data(), width, height);      /* print the time it took */     std::chrono::time_point<std::chrono::high_resolution_clock> const t2 = std::chrono::high_resolution_clock::now();     std::chrono::duration<double> const duration =         std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);     double const time_took = duration.count();     std::printf("It took %f %s", time_took, std::array{ "seconds", "second" } [1.0 == time_took]);      /* wait for user input to close */     (void)std::getchar(); } 

here is what the program produces:

enter image description here

  
 
 

Lista de respuestas

7
 
vote
vote
La mejor respuesta
 

Nombres de variables

Hay muchas variables con nombres muy cortos. Sé que es común en las fórmulas matemáticas usar nombres de una sola letra para las variables, pero al menos encontrará algún texto acompañante que explica qué significa todas las letras. Al menos, agregaría algunos comentarios al código en el lugar donde declara variabels como double m_a . Esto puede ser una corta epiglainación, o quizás una referencia a un artículo o un libro, incluido el número de la fórmula donde se introduce por primera vez. Alternativamente, dale a las variables un nombre más largo, pero más descriptivo.

Por otro lado, algunas variables son poco largas y pueden ser cortadas. Por ejemplo, las variables que se refieren a matrices o vectores se escriben comúnmente utilizando la forma plural, y no tenemos que repetir el tipo en el nombre. Así, por ejemplo, en lugar de coefficient_vector_ptr , escriba coefficients .

Uso std::complex Para variables complejas

En lugar de declarar dos variables, una para lo real y la otra para la parte imaginaria, considere declarar una sola std::complex Variable. Todavía puede acceder a los dos componentes individualmente si es necesario, pero reduce la cantidad de variables, y también hay muchas funciones matemáticas que pueden funcionar directamente en variables complejas.

Uso const Los punteros donde corresponde

Veo un poco de uso de constexpr , pero casi ninguna ocurrencia de import pandas as pd d = { 'p2': [(1355150022, 3.82), (1355150088, 4.24), (1355150154, 3.94), (1355150216, 3.87), (1355150287, 6.66)], 'p1': [(1355150040, 7.9), (1355150110, 8.03), (1355150172, 8.41), (1355150234, 7.77)], 'p3': [(1355150021, 1.82), (1355150082, 2.24), (1355150153, 3.33), (1355150217, 7.77), (1355150286, 6.66)], } ret = pd.DataFrame({k:pd.Series({a-a%60:b for a,b in v}) for k,v in d.iteritems()}) 0 . Siempre que esté pasando un puntero a algo a una función, y no está modificando los contenidos, lo convierten en un puntero 998877766554443311 . Esto capturará errores si accidentalmente escribe a una variable 99887766655443312 , y podría darle al compilador algunas oportunidades más para optimizar el código. Por ejemplo, import pandas as pd d = { 'p2': [(1355150022, 3.82), (1355150088, 4.24), (1355150154, 3.94), (1355150216, 3.87), (1355150287, 6.66)], 'p1': [(1355150040, 7.9), (1355150110, 8.03), (1355150172, 8.41), (1355150234, 7.77)], 'p3': [(1355150021, 1.82), (1355150082, 2.24), (1355150153, 3.33), (1355150217, 7.77), (1355150286, 6.66)], } ret = pd.DataFrame({k:pd.Series({a-a%60:b for a,b in v}) for k,v in d.iteritems()}) 3 import pandas as pd d = { 'p2': [(1355150022, 3.82), (1355150088, 4.24), (1355150154, 3.94), (1355150216, 3.87), (1355150287, 6.66)], 'p1': [(1355150040, 7.9), (1355150110, 8.03), (1355150172, 8.41), (1355150234, 7.77)], 'p3': [(1355150021, 1.82), (1355150082, 2.24), (1355150153, 3.33), (1355150217, 7.77), (1355150286, 6.66)], } ret = pd.DataFrame({k:pd.Series({a-a%60:b for a,b in v}) for k,v in d.iteritems()}) 4 en , y import pandas as pd d = { 'p2': [(1355150022, 3.82), (1355150088, 4.24), (1355150154, 3.94), (1355150216, 3.87), (1355150287, 6.66)], 'p1': [(1355150040, 7.9), (1355150110, 8.03), (1355150172, 8.41), (1355150234, 7.77)], 'p3': [(1355150021, 1.82), (1355150082, 2.24), (1355150153, 3.33), (1355150217, 7.77), (1355150286, 6.66)], } ret = pd.DataFrame({k:pd.Series({a-a%60:b for a,b in v}) for k,v in d.iteritems()}) 6 en import pandas as pd d = { 'p2': [(1355150022, 3.82), (1355150088, 4.24), (1355150154, 3.94), (1355150216, 3.87), (1355150287, 6.66)], 'p1': [(1355150040, 7.9), (1355150110, 8.03), (1355150172, 8.41), (1355150234, 7.77)], 'p3': [(1355150021, 1.82), (1355150082, 2.24), (1355150153, 3.33), (1355150217, 7.77), (1355150286, 6.66)], } ret = pd.DataFrame({k:pd.Series({a-a%60:b for a,b in v}) for k,v in d.iteritems()}) 7 puede ser Hecho import pandas as pd d = { 'p2': [(1355150022, 3.82), (1355150088, 4.24), (1355150154, 3.94), (1355150216, 3.87), (1355150287, 6.66)], 'p1': [(1355150040, 7.9), (1355150110, 8.03), (1355150172, 8.41), (1355150234, 7.77)], 'p3': [(1355150021, 1.82), (1355150082, 2.24), (1355150153, 3.33), (1355150217, 7.77), (1355150286, 6.66)], } ret = pd.DataFrame({k:pd.Series({a-a%60:b for a,b in v}) for k,v in d.iteritems()}) 8 punteros.

Evite la fundición innecesariamente

Veo muchos castillos que parecen innecesarios. Por ejemplo:

  import pandas as pd  d = {             'p2': [(1355150022, 3.82), (1355150088, 4.24), (1355150154, 3.94), (1355150216, 3.87), (1355150287, 6.66)],             'p1': [(1355150040, 7.9), (1355150110, 8.03), (1355150172, 8.41), (1355150234, 7.77)],             'p3': [(1355150021, 1.82), (1355150082, 2.24), (1355150153, 3.33), (1355150217, 7.77), (1355150286, 6.66)], }   ret = pd.DataFrame({k:pd.Series({a-a%60:b for a,b in v}) for k,v in d.iteritems()}) 9  

¿Por qué los moldes cuando p1 p2 p3 1355149980 NaN 3.82 1.82 1355150040 7.90 4.24 2.24 1355150100 8.03 3.94 3.33 1355150160 8.41 3.87 7.77 1355150220 7.77 NaN NaN 1355150280 NaN 6.66 6.66 0 ya es una matriz de p1 p2 p3 1355149980 NaN 3.82 1.82 1355150040 7.90 4.24 2.24 1355150100 8.03 3.94 3.33 1355150160 8.41 3.87 7.77 1355150220 7.77 NaN NaN 1355150280 NaN 6.66 6.66 1 s, y 99887766555443322 es también un 99887776655443323 ? También evitaría el uso de la versión de la biblioteca C de p1 p2 p3 1355149980 NaN 3.82 1.82 1355150040 7.90 4.24 2.24 1355150100 8.03 3.94 3.33 1355150160 8.41 3.87 7.77 1355150220 7.77 NaN NaN 1355150280 NaN 6.66 6.66 4 , y use p1 p2 p3 1355149980 NaN 3.82 1.82 1355150040 7.90 4.24 2.24 1355150100 8.03 3.94 3.33 1355150160 8.41 3.87 7.77 1355150220 7.77 NaN NaN 1355150280 NaN 6.66 6.66 5 en su lugar:

                p1    p2    p3 1355149980   NaN  3.82  1.82 1355150040  7.90  4.24  2.24 1355150100  8.03  3.94  3.33 1355150160  8.41  3.87  7.77 1355150220  7.77   NaN   NaN 1355150280   NaN  6.66  6.66 6  

También tenga en cuenta que C ++, para bien o para mal, realizará moldes implícitos y proporcione promociones para usted en algunos casos. Generalmente reducen la cantidad de fundición necesaria. TOMAR Por ejemplo:

                p1    p2    p3 1355149980   NaN  3.82  1.82 1355150040  7.90  4.24  2.24 1355150100  8.03  3.94  3.33 1355150160  8.41  3.87  7.77 1355150220  7.77   NaN   NaN 1355150280   NaN  6.66  6.66 7  

Esto se puede reescribir a:

                p1    p2    p3 1355149980   NaN  3.82  1.82 1355150040  7.90  4.24  2.24 1355150100  8.03  3.94  3.33 1355150160  8.41  3.87  7.77 1355150220  7.77   NaN   NaN 1355150280   NaN  6.66  6.66 8  

Tenga en cuenta que no solo es más corto, es aún más eficiente en este caso: 99887776655443329 tiene una sobrecarga para los exponentes enteros, y puede usar un algoritmo mucho más rápido para calcular el resultado en ese caso.

Evite la elevación de los casos degenerados fuera de los bucles

Veo que este patrón repitió muchas veces:

  double m_a0  

Aquí usted trata double m_a1 como un caso especial y lo ha movido fuera del circuito. Pero se puede reescribir a:

  double m_a2  

Probablemente no haya diferencia en la velocidad, pero este último es simplemente un código más simple, y le dice que en realidad no hay nada especial sobre el primer elemento.

Uso double m_a3 Para evitar repetir (largos) tipos

Si bien no usaría double m_a4 para la mayoría de las matemáticas, se puede usar efectivamente dentro de double m_a5 para evitar repetirse. Por ejemplo:

  double m_a6  
 

Variable names

There are a lot of variables with very short names. I know it is common in mathematical formulas to use single-letter names for variables, but you will at least find some accompanying text that explains what all the letters means. I would at least add some comments to the code at the place where you declare variabels like double m_a. This can be a short epxlaination, or perhaps a reference to a paper or book, including the number of the formula where it is first introduced. Alternatively, give the variables a longer, but more descriptive name.

On the other hand, some variables are bit long, and can be shorted. For example, variables referring to arrays or vectors are commonly written using the plural form, and we don't have to repeat the type in the name. So for example, instead of coefficient_vector_ptr, write coefficients.

Use std::complex for complex variables

Instead of declaring two variables, one for the real and the other for the imaginary part, consider declaring a single std::complex variable. You can still access the two components individually if needed, but it reduces the amount of variables, and there are also a lot of mathematical functions that can work directly on complex variables.

Use const pointers where appropriate

I see some use of constexpr, but almost no occurence of const. Whenever you are passing a pointer to something to a function, and you are not modifying the contents, make it a const pointer. This will catch errors if you accidentily do write to a const variable, and it might give the compiler some more opportunities to optimize the code. For example, filename and image_data in write_image(), and coefficient_vector_ptr in FindRoots() can be made const pointers.

Avoid casting unnecessarily

I see a lot of casts that seem unnecessary. For example:

xvar = (double)(::fabs((double)(m_p_vector[ii]))); 

Why the casts when m_p_vector is already an array of doubles, and xvar is also a double? I would also avoid using the C library version of fabs(), and use std::fabs() instead:

xvar = std::fabs(m_p_vector[ii]); 

Also note that C++, for better or for worse, will perform implicit casts and type promotions for you in some cases. They generally reduce the amount of casting necessary. Take for example:

std::int32_t lvar = (std::int32_t)(::log(sc) / ::log(f_BASE) + 0.5); double factor = ::pow((double)(f_BASE * 1.0), double(lvar)); 

This can be rewritten to:

std::int32_t lvar = std::log(sc) / std::log(f_BASE) + 0.5; double factor = std::pow(f_BASE * 1.0, lvar); 

Note that it not only is shorter, it is even more efficient in this case: std::pow() has an overload for integer exponents, and can use a much faster algorithm to calculate the result in that case.

Avoid hoisting degenerate cases out of loops

I see this pattern repeated a lot of times:

double kv = m_k_vector[0]; m_qk_vector[0] = kv;  for (ii = 1; ii < m_n; ++ii) {     kv = kv * svar + m_k_vector[ii];     m_qk_vector[ii] = kv; } 

Here you treat ii = 0 as a special case and have moved it out of the loop. But it can be rewritten to:

double kv = 0;  for (ii = 0; i < m_n; ++ii) {     kv = kv * svar + m_k_vector[ii];     m_qk_vector[ii] = kv; } 

There is probably no difference in speed, but the latter is just simpler code, and it tells you that there is actually nothing special about the first element.

Use auto to avoid repeating (long) types

While I wouldn't use auto for most of the maths, it can be used effectively inside main() to avoid repeating yourself. For example:

auto const t1 = std::chrono::high_resolution_clock::now(); ...     auto color = static_cast<std::uint8_t>((std::log(heatmap[i]) / std::log(max_value)) * 255.0 + 0.55555); ... auto const t2 = std::chrono::high_resolution_clock::now(); auto duration = t2 - t1; 
 
 

Relacionados problema

10  Quicksort  ( Templated quicksort ) 
original quicksort.h #include <algorithm> namespace quicksort { template <typename iterator, typename value_type> struct traits { static iterato...

1  Integración de oscilador de fase perturbada  ( Perturbed phase oscillator integration ) 
Estoy integrando un sistema de osciladores de fase perturbados. Defino el sistema de ecuación y también la matriz jacobiana. Tengo que remodelar el vector dim...

2  Simplifique el código de verificación de banderas de bits  ( Simplify bit flags checking code ) 
Por favor, ayuda a hacer que este código sea más limpio. Es parte del procedimiento de la ventana que notifica a mi renderizado antes de Área de cliente de ...

9  Entrada de usuario y lectura de contenidos de archivo  ( User input and reading contents of file ) 
Para la divulgación completa: esta es una tarea para mi clase de programación y solo quiero consejos o consejos sobre algunos del código. Detalles de asigna...

14  Implementando una lista relacionada adecuada para un entorno profesional  ( Implementing a proper linked list for a professional environment ) 
Tengo algunas preocupaciones: ¿Es normal que la clase tenga al menos un nodo? En otras palabras, esta implementación no puede tener una lista vinculada va...

1  Piedra Papel tijeras  ( Rock paper scissors ) 
Gracias por su tiempo, soy nuevo en la programación y pasé algunos días haciendo este rock, papel y amp; Juego de tijera. ¿Qué otras mejoras posibles podrían ...

0  Sistema de gestión de la biblioteca en C ++  ( Library management system in c ) 
Estoy haciendo un proyecto de principiante C ++ y es simplemente un sistema de administración de la biblioteca donde contratará un estudiante, devolver un lib...

3  Implementación más portátil de Tolower ()  ( More portable tolower implementation ) 
Me estoy desafiando a intentar intentar escribir una función que sea tan eficiente, portátil y a prueba de fallas posible. La función es muy simple y solo con...

4  Simulación simple de red neural en C ++ (Ronda 2)  ( Simple neural network simulation in c round 2 ) 
Intro Ayer He publicado esta pregunta . Desde entonces, he actualizado mi código para incorporar estas sugerencias . También he eliminado la dependencia d...

21  Algoritmo de suma de comprobación personalizada  ( Custom checksum algorithm ) 
A MIENTROS RUENTANDO, INVERSO: diseñé un algoritmo de suma de comprobación de una MMO que se usa para comprobar el Validez de un artículo que está vinculado...




© 2022 respuesta.top Reservados todos los derechos. Centro de preguntas y respuestas reservados todos los derechos