root finding - function won't exit?

[sexyzebra19]

I'm trying to implement the algorithm described here under 'definition', http://en.wikipedia.org/wiki/Laguerre%27s_method

To summarise it takes an initial guess 0, and then iterates to find the root.

In lines 68 to 71 I was hoping this would exit the function, i.e. when the root is found...but it's still carrying on.

Can anyone please explain why?

Thanks in advance!

#include <iostream>
#include <cmath>
#include <vector>
#include <iomanip>

using namespace std;

//Declare functions
void poly(vector<double> A, double x, double& p, double& px, double& pxx, double a, double G, double H, double denom1, double denom2, int degree, int max);

//---------------------------------------------------------------------------------------

int main()
{
    vector<double> A;                  // vector of coefficients a^degree...a^0
    int degree;                        // highest degree
    int max = 0;                       // maximum number of iterations
    double x;                          // value of x, initial guess
    double p;                          // p(x)
    double px;                         // p'(x)
    double pxx;                        // p''(x)
    double a;
    double G;
    double H;
    double denom1;
    double denom2;

    poly(A, x, p, px, pxx, a, G, H, denom1, denom2, degree, max);
    return 0;
}

//---------------------------------------------------------------------------------------

void poly(vector<double> A, double x, double& p, double& px, double& pxx, double a, double G, double H, double denom1, double denom2, int degree, int max)
{
    cout << "For polynomial p(x) = a0 + a1x + a2x^2 + ... + anx^n" << endl;
    cout << "Enter the polynomial degree, n" << endl;
    cout << "Example:  2 for a quadratic, 3 for a cubic..." << endl;
    cout << "Degree: ";
    cin >> degree;

    A.resize (degree + 1);

    for (int i = degree; i >= 0; i--)
    {
        cout << "Enter coefficient a" << i << ": ";
        cin >> A[i];
    }

    cout << "Enter the initial guess x: ";
    cin >> x;
    cout << "Enter the maximum number of iterations allowed: ";
    cin >> max;


    for (int i = 1; i <= max; i++)
    {

    //Calculations for p(x)

    p = A[degree];
    for (int i = (degree - 1); i >= 0; i--)
    {
        p = p*x;
        p = p+A[i];
    }

    if (p==0)
    {
        return;
    }

    //First differentiation on the coefficients
    for (int i = 1; i <= degree; i++)
    {
        A[i - 1] = i * A[i];
    }

    //Horner's Method - calculations for p'(x)
    px = A[degree - 1];
    for (int i = (degree - 2); i >= 0; i--)
    {
        px = px*x;
        px = px+A[i];
    }

    //Second differentiation on the coefficients
    for (int i = 1; i <= degree-1 ; i++)
    {
        A[i - 1] = i * A[i];
    }

    //Horner's Method - calculations for p''(x)
    pxx = A[degree - 2];
    for (int i = (degree - 3); i >= 0; i--)
    {
        pxx = pxx*x;
        pxx = pxx+A[i];
    }

    //Laguerre's Method
    G = px/p;
    H = G*G - pxx/p;
    denom1 = (G + sqrt((degree-1)*(degree*H - G*G)));
    denom2 = (G - sqrt((degree-1)*(degree*H - G*G)));

    if ( abs(denom1) >= abs(denom2) )
    {
        a = degree/denom1;
    }
    else
    {
        a = degree/denom2;
    }

    cout << "The root approximation for x = " <<fixed<<setprecision(7)<< x << " is " <<fixed<<setprecision(7)<< a << endl;

    x = x - a;
}
return;
}

[jonsca]

This time you don't need all that precision. Take the fixed and setprecision statements off of the cout(I only told you 7 digits last time so you could see error at 10^-6). P is missing zero by a minuscule amount and continuing. There may be other errors because the approx is off by about 1. Also, you should check your denominator in your Laguerre's calculation against zero to make sure that won't crash the program.
I'm not sure why changing the precision would affect the calculation?? aside from the output. One thing though that this brings up is that it is difficult if not impossible to use == with floats or doubles due to these conditions where you can be so close to zero and not quite zero. So you can do things like myvariable <= 0.00000000001 or something (or take the difference between your variable and an ideal myvariable - ideal < 0.000000001 etc.).

[Clinton Portis]

perhaps you could floor() the variable prior to testing for zero.