Durand-Kerner Method Polynomial Roots

Hello world! I've been learning C++ for a few weeks now. I already know MATLAB and python, so I'm not a complete beginner. I've written the following code using the Durand-Kerner method to find all complex roots of a polynomial. Thoughts? Suggestions? I'm new to C++, but not to programming. Are there any things I could be doing more efficiently? Thanks!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include <iostream>
#include <cmath>
#include <complex>
#include <cstdlib>
using namespace std;
typedef complex<double> dcomp;

// "poly" evaluates at x a polynomial of the form:
// f(x) = pow(x, n) + a1*pow(x, n - 1) + a2*pow(x, n - 2) + . . . + a(n - 2)*x*x + a(n - 1)*x + a(n)
// where the vector A = {a1, a2, a3, . . . , a(n - 2), a(n - 1), a(n)}
dcomp poly(double A[], int n, dcomp x)
{
    dcomp y = pow(x, n);
    for (int i = 0; i < n; i++)
        y += A[i]*pow(x, (n - i - 1) );
    return y;
}

// polyroot uses the Durand-Kerner method to find all roots (real and complex) of a polynomial of the form:
// f(x) = pow(x, n) + a1*pow(x, n - 1) + a2*pow(x, n - 2) + . . . + a(n - 2)*x*x + a(n - 1)*x + a(n)
// where the vector A = {a1, a2, a3, . . . , a(n - 2), a(n - 1), a(n)}
dcomp * polyroot(double A[], int n)
{
    int iterations = 1000;
    dcomp z = dcomp (0.4, 0.9);
    int size = sizeof(z);
    dcomp * R;
    R = (dcomp *) malloc(size*n);
    for (int i = 0; i < n; i++)
        R[i] = pow(z, i);
    for (int i = 0; i < iterations; i++)
    {
        for (int j = 0; j < n; j ++)
        {
            dcomp B = poly(A, n, R[j]);
            for (int k = 0; k < n; k++)
            {
                if (k != j)
                B /= R[j] - R[k];
            }
            R[j] -= B;
        }
    }
    return R;
}

// TEST
int main()
{
    double A[6] = { 0, -29, -58, 54, 216, 216 };
    int n = 6;
    dcomp * R = polyroot(A, n);
    for (int i = 0; i < 6; i++)
        cout << "R[" << i << "] = " << R[i] << endl;
    return 0;
}
Last edited on
std::complex<double> isn't a type safe to use with malloc(). Use new, instead.
Also, you're not freeing the pointer returned by polyroot().
Last edited on
Thanks for the advice Hyperion! I didn't realize new[] could be used with variable-sized arrays.
BrentSpinor wrote:
I didn't realize new[] could be used with variable-sized arrays.


new is just C++'s version of C's malloc(). It is most certainly for dynamically sized arrays ;)
Topic archived. No new replies allowed.