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;
}

What output do you get ? Does it appear to be converging, or are the numbers all over the place ? Are you getting NaN outputs ? — Paul R, Sep 03 '13 at 11:43
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). — user1837009, Sep 03 '13 at 11:48

Pavel K · Accepted Answer · 2013-09-03T13:27:36.303

8

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;
}

edited Sep 03 '13 at 13:27

answered Sep 03 '13 at 11:46

Pavel K

3,541
2
29
44

3

This is C++. Rather than using `fabs`, it might be better include the correct header (``) and use `std::abs`, which is correctly overloaded for most numeric types. – Mike Seymour Sep 03 '13 at 11:51
@Mike: +1 - the original code wouldn't even compile with my g++ due to the ambiguity of abs with `` - changing to the correct header `` fixes the problem of course. – Paul R Sep 03 '13 at 11:55
1

`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 `. – Pete Becker Sep 03 '13 at 12:54
That what I've mean. That some people could use `abs` with `double` argument and with the `` 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 ). **That** is dangerous, definitely not the `abs` function itself. – Pavel K Sep 03 '13 at 13:39
@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`. – robson3.14 Sep 03 '13 at 14:19
2

@robson3.14: No, the C++-specific overloads defined in `` 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 ``. – Mike Seymour Sep 03 '13 at 16:18
2

@robson3.14 - 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 `` that are not in `math.h`, including `abs(float)`, `abs(double)`, and `abs(long double)`. – Pete Becker Sep 03 '13 at 16:39

score 1 · Answer 2 · answered Dec 15 '13 at 15:53

1

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.

answered Dec 15 '13 at 15:53

Lutz Lehmann

25,219
2
22
51

Infinite loop calculating cubic root

2 Answers2