/* * Newton's method is an iterative method to approximate * the roots of a real function taking the following form: * x_(n+1) = x_n - f(x_n) / f'(x_n) * * For an arbitrary number y and it's cube x := y^3, * we know that x^(1/3) = y. Therefore we define our function * f(x) = x^3 - y * y being the number we want the cube root for * and x being it's cube root, and it's derivative * f'(x) = 3x^2 * * We'll therefore find x using the formula * x_(n+1) = x_n - (x_n^3 - y) / (3x_n^2) * We could shorten it, but we want to have as little divisions * as possible, as they are slow and prone to rounding errors. * * We need a starting value x0 that's at least as close to 0 as possible. * To find this value, we will use argmin f(x) for an integer x. */ #include double abs(const double x) { return x >= 0 ? x : -x; } double f(const double x, const double y) { return x*x*x - y; } double f_prime(const double x) { return 3*x*x; } int find_x0(const double y) { int argmin = 0; int min = abs(f(argmin, y)); while(true) { if (abs(f(argmin - 1, y)) < min) --argmin; else if (abs(f(argmin + 1, y)) < min) ++argmin; else return argmin; min = abs(f(argmin, y)); } } int main() { const double y = 100; const double precision = 0.0001; const int max_iterations = 9; int iterations = 0; double x_n = find_x0(y); double x_n_1 = x_n; do { x_n = x_n_1; std::cout << "x_" << iterations << " = " << x_n << std::endl; x_n_1 = x_n - f(x_n, y) / f_prime(x_n); if (++iterations > max_iterations) { std::cout << "Newton's method diverged or the wanted precision was not reached.\n"; break; } } while (abs(x_n_1 - x_n) > precision); std::cout << "\nThe cube root of " << y << " is " << x_n; return 0; }

Cube Root (Newton's method)