Horner-Newton Method

[tswvyim]

Okay, so I'm trying to build a program to find the complex roots of any given polynomials and have run into a dead end. No errors are given when compiling, however, after asking for an initial guess it won't return any other sort of value. I'm assuming its some sort of syntax error because I believe that the actual algorithms are sound. Here's the code

#include <iostream>
#include <complex>
#include <cmath>

using namespace std;

  typedef complex<double> dcmplx; 

int main()
{
  dcmplx a[20],b[20],c[20];
    dcmplx x;
    int n,i;

    cout<<"Enter the degree of equations : ";cin>>n;
    cout<<"Enter all coefficients \n";
    for(i=0;i<=n;i++)
        cin>>a[i];
    cout<<"enter initial guess:";cin>>x;

    while(n>0)
    {
        c[0]=b[0]=a[0];

        for(i=1;i<=n;i++)
            b[i]=a[i]+b[i-1]*x;

        for(i=1;i<n;i++)
            c[i]=b[i]+c[i-1]*x;

        if(abs(c[n-1])==0.0)   
        {
            cout<<"Divide by Zero : Undefined ";
            break;
        }

        x=b[n]/c[n-1];
        if(abs(b[n]/c[n-1])<0.00001)
        {
            cout<<"\nRoot is: "<<x<<"f(x) "<<b[n];
            for(i=1;i<n;i++)
                a[i]=b[i];
            n--;
        }
    }
return 0;
}

Thanks!

[codewalker]

Look at

x=b[n]/c[n-1];
        if(abs(b[n]/c[n-1])<0.00001)

Value of n is not getting modified, thus abs(b[n]/c[n-1]) is the same value all the time and the execution never enters the if block (and hence the value of n never gets modified)

[BrentSpinor]

I thought Horner's method could only be used to find the real roots of a polynomial? Anyways, I've written the following implementation of the Durand-Kerner method:

#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 = new dcomp [size*n];
    //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()
{
// Test the Durand-Kerner method for finding complex roots of a polynomial
    double A[3] = {-3, 3, -5};
    int m = 3;
    dcomp * R = polyroot(A, m);
    for (int i = 0; i < 3; i++)
        cout << "R[" << i << "] = " << R[i] << endl;
    delete R;
    R = 0;

    cout << "\n\n";

    double B[6] = { 0, -29, -58, 54, 216, 216 };
    int n = 6;
    dcomp * S = polyroot(B, n);
    for (int i = 0; i < 6; i++)
        cout << "S[" << i << "] = " << S[i] << endl;
    delete S;
    S = 0;

    return 0;
}

[tswvyim]

Okay, so in order to modify n do I throw in another for loop?