Algorithms_in_C++ 1.0.0
Set of algorithms implemented in C++.
Loading...
Searching...
No Matches
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}
#define ACCURACY
Definition durand_kerner_roots.cpp:45

◆ 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}
T snprintf(T... args)
T imag(T... args)
T real(T... args)
Here is the call graph for this function:

◆ 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}
bool check_termination(long double delta)
Definition durand_kerner_roots.cpp:91
const char * complex_str(const std::complex< double > &x)
Definition durand_kerner_roots.cpp:76
std::complex< double > poly_function(const std::valarray< double > &coeffs, std::complex< double > x)
Definition durand_kerner_roots.cpp:53
T exit(T... args)
T is_open(T... args)
T isinf(T... args)
T isnan(T... args)
T max(T... args)
T open(T... args)
T perror(T... args)
Here is the call graph for this function:

◆ main()

int main ( int argc,
char ** argv )
283 {
284 /* initialize random seed: */
285 std::srand(std::time(nullptr));
286
287 if (argc < 2) {
288 test1(); // run tests when no input is provided
289 test2(); // and skip tests when input polynomial is provided
290 std::cout << "Please pass the coefficients of the polynomial as "
291 "commandline "
292 "arguments.\n";
293 return 0;
294 }
295
296 int n, degree = argc - 1; // detected polynomial degree
297 std::valarray<double> coeffs(degree); // create coefficiencts array
298
299 // number of roots = degree - 1
301
302 std::cout << "Computing the roots for:\n\t";
303 for (n = 0; n < degree; n++) {
304 coeffs[n] = strtod(argv[n + 1], nullptr);
305 if (n < degree - 1 && coeffs[n] != 0)
306 std::cout << "(" << coeffs[n] << ") x^" << degree - n - 1 << " + ";
307 else if (coeffs[n] != 0)
308 std::cout << "(" << coeffs[n] << ") x^" << degree - n - 1
309 << " = 0\n";
310
311 /* initialize root approximations with random values */
312 if (n < degree - 1) {
313 s0[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
314 s0[n] -= 50.f;
315 s0[n] /= 50.f;
316 }
317 }
318
319 // numerical errors less when the first coefficient is "1"
320 // hence, we normalize the first coefficient
321 {
322 double tmp = coeffs[0];
323 coeffs /= tmp;
324 }
325
326 clock_t end_time, start_time = clock();
327 auto result = durand_kerner_algo(coeffs, &s0, true);
328 end_time = clock();
329
330 std::cout << "\nIterations: " << result.first << "\n";
331 for (n = 0; n < degree - 1; n++)
332 std::cout << "\t" << complex_str(s0[n]) << "\n";
333 std::cout << "absolute average change: " << result.second << "\n";
334 std::cout << "Time taken: "
335 << static_cast<double>(end_time - start_time) / CLOCKS_PER_SEC
336 << " sec\n";
337
338 return 0;
339}
T clock(T... args)
void test2()
Definition durand_kerner_roots.cpp:242
void test1()
Definition durand_kerner_roots.cpp:207
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
uint64_t result(uint64_t n)
Definition fibonacci_sum.cpp:76
T rand(T... args)
T srand(T... args)
T strtod(T... args)
T time(T... args)

◆ 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}
T pow(T... args)
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: