Calculating $\sqrt[3]{x}$ with Babylonian method

Question

Consider my attempt to implement the Babylonian method in C:

int sqrt3(int x) {
    double abs_err = 1.0;
    double xold = x;
    double xnew = 0;

    while(abs_err > 1e-8) {
        xnew =  (2 * xold + x/(xold* xold))/3;
        abs_err= xnew-xold;
        if (abs_err < 0) abs_err = -abs_err;
        xold=xnew;
    }
    return xnew;
} 
int main() {
    int a;

    scanf("%d", &a);
    printf(" Result is: %f",sqrt3(a));
    return 0;   
}

Result is for x=27: 0.0000? Where is my mistake?

`1/3` performs an integer (truncating) division, so the result is `0`. Either do something like `1.0/3`, or move `/3` to the end of the expression (so the lhs of `/` is a `double` rather than `int`). — HolyBlackCat, May 26 '20 at 16:01
Please note that SO site doesn't provide MathJax, so that it doesn't interpret LaTeX. BTW, it seems to be a *cube* root, rather than a sqrt3. — Bob__, May 26 '20 at 16:07
Yes sorry. It should be the a cube root. But why does my programm return 0 0.0000? — Steven Hunt, May 26 '20 at 16:14
you are quite liberal with number types. Maybe it is a good idea to first implement the function in a more type safe language. I also have my doubts about the line where you update ``xnew`` AND I think the initial value of xold is usually x and xnew starts with 1.0, not 0.0. This Julia function seems to work: `` function cubeRoot( a ) x = a y = 1.0 precision = 0.000001 while (abs(x-y) > precision) x = (x+y) / 2 y = a / x / x end x end`` — BitTickler, May 26 '20 at 16:24
@BitTickler. Thank you for your advice. But where exactly do you see my problem? For exmaple: /3 is numerical problematic if you conisder floating point numbers, isn't it? — Steven Hunt, May 26 '20 at 16:30
It starts with the naming. In my Julia code above, x is the upper bound, y the lower bound and they move towards each other. ``xold, xnew`` do not clarify that point. Then, if you intend to do floating point arithmetic, your constants should be 2.0 instead of 2 and 3.0 instead of 3 to eliminate any possibility for erroneous implicit conversions. More readable and less to worry about. — BitTickler, May 26 '20 at 16:39

score 0 · Accepted Answer · answered May 26 '20 at 16:19

0

While the function returns an int, that value is printed with the wrong format specifier, %f instead of %d.

Change the signature (and the name, if I may) into something like this

double cube_root(double x) { ... }

Or change the format specifier, if you really want an int.

answered May 26 '20 at 16:19

Bob__

12,361
3
28
42

Thank you :) I did not see that my return type is int :) – Steven Hunt May 26 '20 at 16:26

score 0 · Answer 2 · answered May 26 '20 at 17:23

Following the explanation from tutorialspoint, which states, that the basic idea is to implement the Newton Raphson method for solving nonlinear equations, IMHO, the code below displays this fact more clearly. Since there is already an accepted answer, I add this answer just for future reference.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

double rootCube( double a)
{
    double x = a;
    double y = 1.0;
    const double precision = 0.0000001;
    while(fabs(x-y) > precision)
    {
        x = (x + y) / 2.0;
        y = a / x / x;
    }
    return x;
}

int main(int argc, const char* argv[])
{
    if(argc > 1)
    {
        double a = 
            strtod(argv[1],NULL); 
        printf("cubeRoot(%f) = %f\n", a, rootCube(a));
    }

    return 0;
}

Here, in contrast to the original code of the question, it is more obvious, that x and y are the bounds, which are being improved until a sufficiently accurate solution is found.

With modification of the line in the while block, where y is being updated, this code can also be used to solve similar equations. For finding the square root, for example, this line would look like this: y = a / x.

Please note that while there are several way to improve OP's code (the estimation of the error for starters), the algorithm showed in this tutorialspoint link doens't really look as a Newton's method. At least it doesn't seem to have a quadratic convergence rate (see e.g. [here](https://godbolt.org/z/i6arkE)), it's more like a bisection method. — Bob__, May 26 '20 at 19:23
@Bob__ I really should avoid that tutorialspoint site. This is not the first time I found wrong information there. Thanks for pointing it out. — BitTickler, May 27 '20 at 18:03

Calculating $\sqrt[3]{x}$ with Babylonian method

2 Answers2