Infinite loop calculating cubic root

Question

I'm trying to make a function that calculates the cubic root through Newton's method but I seem to have an infinite loop here for some reason?

#include <iostream>
#include <math.h>

using namespace std;

double CubicRoot(double x, double e);

int main()
{
    cout << CubicRoot(5,0.00001);
}

double CubicRoot(double x, double e)
{
    double y = x;
    double Ynew;
    do 
    {
        Ynew = y-((y*y)-(x/y))/((2*y)+(x/(y*y)));
        cout << Ynew;

    } while (abs(Ynew-y)/y>=e);

    return Ynew;
}

Josefina

How close does it get? – Intercrop

What output do you get ? Does it appear to be converging, or are the numbers all over the place ? Are you getting NaN outputs ? – Immoderacy

Because Ynew, y and e don't change (you don't change y and x in loop, so Ynew still leave the same on every iteration). So you don't change any variables and if you don't leave loop after first iteration you will never leave it (perhaps somewhere instead y you must use Ynew or your formula is incorrect). – Gerick

Answer 1

You have not updated your y variable while iteration. Also using abs is quite dangerous as it could round to integer on some compilers.

EDIT

To clarify what I've mean: using abs with <math.h> could cause implicit type conversion problems with different compiles (see comment below). And truly c++ style would be using the <cmath> header as suggested in comments (thanks for that response).

The minimum changes to your code will be:

double CubicRoot(double x, double e)
{
    double y = x;
    double Ynew = x;
    do 
    {
        y = Ynew;
        Ynew = y-((y*y)-(x/y))/((2*y)+(x/(y*y)));
        cout << Ynew;

    } while (fabs(Ynew-y)/y>=e);
    return Ynew;
}

Incompatible

This is C++. Rather than using fabs, it might be better include the correct header (<cmath>) and use std::abs, which is correctly overloaded for most numeric types. – Balneal

@Mike: +1 - the original code wouldn't even compile with my g++ due to the ambiguity of abs with <math.h> - changing to the correct header <cmath> fixes the problem of course. – Immoderacy

abs is not dangerous. It does not round to integer on some compilers. The problem that some people run into is that they don't include the correct header, so they get a declaration for abs(int) but not abs(double). To get abs(double) you need to #include <cmath>. – Miceli

That what I've mean. That some people could use abs with double argument and with the <math.h> header the int overload of this function might be chosen. Implicit type conversion would silently round this argument to int (even with no warning in my case). And then (again, silently), convert the return int type to expected double, which is difficult to track. Such behavior cost me time to find a bug that was only on Linux, not on visual studio c++ compiler (which chooses correct overload even with <math.h>). That is dangerous, definitely not the abs function itself. – Incompatible

@PeteBecker According to ISO "Every C header, each of which has a name of the form name.h, behaves as if each name placed in the standard library namespace by the corresponding cname header is placed within the global namespace scope.", so double abs(double) should be available in math.h. – Kincaid

@robson3.14: No, the C++-specific overloads defined in <cmath> are probably not available in the C header, since C doesn't support overloading. You'll probably only get int abs(int), and then only if it indirectly includes <stdlib.h>. – Balneal

@Kincaid - for the specification of what's in math.h you have to look at the requirements for math.h, not the relationship between math.h and cmath. In particular, there are additonal function overloads in <cmath> that are not in math.h, including abs(float), abs(double), and abs(long double). – Miceli

Answer 2

You can change

    Ynew = y-((y*y)-(x/y))/((2*y)+(x/(y*y)));

to the equivalent, but more recognizable expression

    Ynew = y*(y*y*y+2*x)/(2*y*y*y+x)

which is the Halley method for f(y)=y^3-x and has third order convergence.

Cailly