Algorithms_in_C++  1.0.0
Set of algorithms implemented in C++.
durand_kerner_roots.cpp File Reference

Compute all possible approximate roots of any given polynomial using Durand Kerner algorithm More...

#include <algorithm>
#include <cassert>
#include <cmath>
#include <complex>
#include <cstdlib>
#include <ctime>
#include <fstream>
#include <iostream>
#include <valarray>
Include dependency graph for durand_kerner_roots.cpp:

Macros

#define ACCURACY   1e-10
 
#define MAX_BUFF_SIZE   50
 

Functions

std::complex< double > poly_function (const std::valarray< double > &coeffs, std::complex< double > x)
 
const char * complex_str (const std::complex< double > &x)
 
bool check_termination (long double delta)
 
std::pair< uint32_t, double > durand_kerner_algo (const std::valarray< double > &coeffs, std::valarray< std::complex< double >> *roots, bool write_log=false)
 
void test1 ()
 
void test2 ()
 
int main (int argc, char **argv)
 

Detailed Description

Compute all possible approximate roots of any given polynomial using Durand Kerner algorithm

Author
Krishna Vedala

Test the algorithm online: https://gist.github.com/kvedala/27f1b0b6502af935f6917673ec43bcd7

Try the highly unstable Wilkinson's polynomial:

./numerical_methods/durand_kerner_roots 1 -210 20615 -1256850 53327946
-1672280820 40171771630 -756111184500 11310276995381 -135585182899530
1307535010540395 -10142299865511450 63030812099294896 -311333643161390640
1206647803780373360 -3599979517947607200 8037811822645051776
-12870931245150988800 13803759753640704000 -8752948036761600000
2432902008176640000

Sample implementation results to compute approximate roots of the equation \(x^4-1=0\):
Error evolution during root approximations computed every
iteration. Roots evolution - shows the initial approximation of the
roots and their convergence to a final approximation along with the iterative
approximations

Macro Definition Documentation

◆ ACCURACY

#define ACCURACY   1e-10

maximum accuracy limit

Function Documentation

◆ check_termination()

bool check_termination ( long double  delta)

check for termination condition

Parameters
[in]deltapoint at which to evaluate the polynomial
Returns
false if termination not reached
true if termination reached
91  {
92  static long double past_delta = INFINITY;
93  if (std::abs(past_delta - delta) <= ACCURACY || delta < ACCURACY)
94  return true;
95  past_delta = delta;
96  return false;
97 }

◆ complex_str()

const char* complex_str ( const std::complex< double > &  x)

create a textual form of complex number

Parameters
[in]xpoint at which to evaluate the polynomial
Returns
pointer to converted string
76  {
77 #define MAX_BUFF_SIZE 50
78  static char msg[MAX_BUFF_SIZE];
79 
80  std::snprintf(msg, MAX_BUFF_SIZE, "% 7.04g%+7.04gj", x.real(), x.imag());
81 
82  return msg;
83 }

◆ durand_kerner_algo()

std::pair<uint32_t, double> durand_kerner_algo ( const std::valarray< double > &  coeffs,
std::valarray< std::complex< double >> *  roots,
bool  write_log = false 
)

Implements Durand Kerner iterative algorithm to compute all roots of a polynomial.

Parameters
[in]coeffscoefficients of the polynomial
[out]rootsthe computed roots of the polynomial
[in]write_logflag whether to save the log file (default = false)
Returns
pair of values - number of iterations taken and final accuracy achieved
111  {
112  long double tol_condition = 1;
113  uint32_t iter = 0;
114  int n;
115  std::ofstream log_file;
116 
117  if (write_log) {
118  /*
119  * store intermediate values to a CSV file
120  */
121  log_file.open("durand_kerner.log.csv");
122  if (!log_file.is_open()) {
123  perror("Unable to create a storage log file!");
124  std::exit(EXIT_FAILURE);
125  }
126  log_file << "iter#,";
127 
128  for (n = 0; n < roots->size(); n++) log_file << "root_" << n << ",";
129 
130  log_file << "avg. correction";
131  log_file << "\n0,";
132  for (n = 0; n < roots->size(); n++)
133  log_file << complex_str((*roots)[n]) << ",";
134  }
135 
136  bool break_loop = false;
137  while (!check_termination(tol_condition) && iter < INT16_MAX &&
138  !break_loop) {
139  tol_condition = 0;
140  iter++;
141  break_loop = false;
142 
143  if (log_file.is_open())
144  log_file << "\n" << iter << ",";
145 
146 #ifdef _OPENMP
147 #pragma omp parallel for shared(break_loop, tol_condition)
148 #endif
149  for (n = 0; n < roots->size(); n++) {
150  if (break_loop)
151  continue;
152 
153  std::complex<double> numerator, denominator;
154  numerator = poly_function(coeffs, (*roots)[n]);
155  denominator = 1.0;
156  for (int i = 0; i < roots->size(); i++)
157  if (i != n)
158  denominator *= (*roots)[n] - (*roots)[i];
159 
160  std::complex<long double> delta = numerator / denominator;
161 
162  if (std::isnan(std::abs(delta)) || std::isinf(std::abs(delta))) {
163  std::cerr << "\n\nOverflow/underrun error - got value = "
164  << std::abs(delta) << "\n";
165  // return std::pair<uint32_t, double>(iter, tol_condition);
166  break_loop = true;
167  }
168 
169  (*roots)[n] -= delta;
170 
171 #ifdef _OPENMP
172 #pragma omp critical
173 #endif
174  tol_condition = std::max(tol_condition, std::abs(std::abs(delta)));
175  }
176  // tol_condition /= (degree - 1);
177 
178  if (break_loop)
179  break;
180 
181  if (log_file.is_open()) {
182  for (n = 0; n < roots->size(); n++)
183  log_file << complex_str((*roots)[n]) << ",";
184  }
185 
186 #if defined(DEBUG) || !defined(NDEBUG)
187  if (iter % 500 == 0) {
188  std::cout << "Iter: " << iter << "\t";
189  for (n = 0; n < roots->size(); n++)
190  std::cout << "\t" << complex_str((*roots)[n]);
191  std::cout << "\t\tabsolute average change: " << tol_condition
192  << "\n";
193  }
194 #endif
195 
196  if (log_file.is_open())
197  log_file << tol_condition;
198  }
199 
200  return std::pair<uint32_t, long double>(iter, tol_condition);
201 }
Here is the call graph for this function:

◆ poly_function()

std::complex<double> poly_function ( const std::valarray< double > &  coeffs,
std::complex< double >  x 
)

Evaluate the value of a polynomial with given coefficients

Parameters
[in]coeffscoefficients of the polynomial
[in]xpoint at which to evaluate the polynomial
Returns
\(f(x)\)
54  {
55  double real = 0.f, imag = 0.f;
56  int n;
57 
58  // #ifdef _OPENMP
59  // #pragma omp target teams distribute reduction(+ : real, imag)
60  // #endif
61  for (n = 0; n < coeffs.size(); n++) {
63  coeffs[n] * std::pow(x, coeffs.size() - n - 1);
64  real += tmp.real();
65  imag += tmp.imag();
66  }
67 
68  return std::complex<double>(real, imag);
69 }
Here is the call graph for this function:

◆ test1()

void test1 ( )

Self test the algorithm by checking the roots for \(x^2+4=0\) to which the roots are \(0 \pm 2i\)

207  {
208  const std::valarray<double> coeffs = {1, 0, 4}; // x^2 - 2 = 0
211  std::complex<double>(0., 2.),
212  std::complex<double>(0., -2.) // known expected roots
213  };
214 
215  /* initialize root approximations with random values */
216  for (int n = 0; n < roots.size(); n++) {
217  roots[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
218  roots[n] -= 50.f;
219  roots[n] /= 25.f;
220  }
221 
222  auto result = durand_kerner_algo(coeffs, &roots, false);
223 
224  for (int i = 0; i < roots.size(); i++) {
225  // check if approximations are have < 0.1% error with one of the
226  // expected roots
227  bool err1 = false;
228  for (int j = 0; j < roots.size(); j++)
229  err1 |= std::abs(std::abs(roots[i] - expected[j])) < 1e-3;
230  assert(err1);
231  }
232 
233  std::cout << "Test 1 passed! - " << result.first << " iterations, "
234  << result.second << " accuracy"
235  << "\n";
236 }
Here is the call graph for this function:

◆ test2()

void test2 ( )

Self test the algorithm by checking the roots for \(0.015625x^3-1=0\) to which the roots are \((4+0i),\,(-2\pm3.464i)\)

242  {
243  const std::valarray<double> coeffs = {// 0.015625 x^3 - 1 = 0
244  1. / 64., 0., 0., -1.};
246  const std::valarray<std::complex<double>> expected = {
247  std::complex<double>(4., 0.), std::complex<double>(-2., 3.46410162),
248  std::complex<double>(-2., -3.46410162) // known expected roots
249  };
250 
251  /* initialize root approximations with random values */
252  for (int n = 0; n < roots.size(); n++) {
253  roots[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
254  roots[n] -= 50.f;
255  roots[n] /= 25.f;
256  }
257 
258  auto result = durand_kerner_algo(coeffs, &roots, false);
259 
260  for (int i = 0; i < roots.size(); i++) {
261  // check if approximations are have < 0.1% error with one of the
262  // expected roots
263  bool err1 = false;
264  for (int j = 0; j < roots.size(); j++)
265  err1 |= std::abs(std::abs(roots[i] - expected[j])) < 1e-3;
266  assert(err1);
267  }
268 
269  std::cout << "Test 2 passed! - " << result.first << " iterations, "
270  << result.second << " accuracy"
271  << "\n";
272 }
Here is the call graph for this function:
check_termination
bool check_termination(long double delta)
Definition: durand_kerner_roots.cpp:91
std::pair
complex_str
const char * complex_str(const std::complex< double > &x)
Definition: durand_kerner_roots.cpp:76
durand_kerner_algo
std::pair< uint32_t, double > durand_kerner_algo(const std::valarray< double > &coeffs, std::valarray< std::complex< double >> *roots, bool write_log=false)
Definition: durand_kerner_roots.cpp:109
std::isnan
T isnan(T... args)
std::snprintf
T snprintf(T... args)
std::cerr
std::ofstream
STL class.
std::isinf
T isinf(T... args)
std::perror
T perror(T... args)
poly_function
std::complex< double > poly_function(const std::valarray< double > &coeffs, std::complex< double > x)
Definition: durand_kerner_roots.cpp:53
std::valarray< double >
std::complex::real
T real(T... args)
std::ofstream::open
T open(T... args)
ACCURACY
#define ACCURACY
Definition: durand_kerner_roots.cpp:45
std::rand
T rand(T... args)
std::complex
STL class.
std::complex::imag
T imag(T... args)
std::max
T max(T... args)
std::ofstream::is_open
T is_open(T... args)
std::exit
T exit(T... args)
std::pow
T pow(T... args)