Newton's Method

Introduction

This article introduces Newton's method for finding approximate solutions of equations. The method was proposed by Newton in the 17th century.

More specifically, given a continuous and monotonic function $f(x)$ on $[a,b]$, we seek an approximate solution to $f(x)=0$.

Explanation

Initially, we start with the given $f(x)$ and an approximate solution $x_0$ (the choice of the initial value is related to Newton fractals; see 3Blue1Brown's Newton fractal).

Suppose that our current approximate solution is $x_i$. Draw the line $l$ tangent to $f(x)$ at $(x_i,f(x_i))$, and let the abscissa of the intersection between $l$ and the $x$-axis be $x_{i+1}$. This is a better approximate solution. Repeat this iteration. From the geometric meaning of the derivative, we obtain:

$$f'(x_i) = \frac{f(x_i)}{x_{i} - x_{i+1}}$$

Rearranging gives the following recurrence:

$$x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}$$

Intuitively, if $f(x)$ is sufficiently smooth, then as the number of iterations increases, $x_i$ gets closer and closer to the solution of the equation.

Newton's method has quadratic convergence, which means that the number of correct digits in the approximate solution doubles after each iteration. For a proof of convergence, see citizendium - Newton method Convergence analysis.

Of course, Newton's method also has drawbacks; see the analysis in Xiaolin Wu - Roots of Equations, pages 18–20.

Computing Square Roots

We try to use Newton's method to compute a square root. Let $f(x)=x^2-n$; an approximate solution of this equation is an approximate value of $\sqrt{n}$. Therefore,

$$x_{i+1}=x_i-\frac{x_i^2-n}{2x_i}=\frac{x_i+\frac{n}{x_i}}{2}$$

In an implementation, choose an appropriate precision. The code is as follows.

Implementation

C++

double sqrt_newton(double n) {
  constexpr static double eps = 1E-15;
  double x = 1;
  while (true) {
    double nx = (x + n / x) / 2;
    if (abs(x - nx) < eps) break;
    x = nx;
  }
  return x;
}

Python

def sqrt_newton(n):
    eps = 1e-15
    x = 1
    while True:
        nx = (x + n / x) / 2
        if abs(x - nx) < eps:
            break
        x = nx
    return x

Computing the Integer Square Root

Although we can call the sqrt() function to obtain a square root, we will also discuss a variant of Newton's method for finding the largest integer solution to the inequality $x^2\le n$. We still consider a process similar to Newton iteration, but make a small change to the boundary condition. If $x$'s approximate solution decreased in the previous iteration but increases in the current iteration, we do not perform the current iteration and exit the loop.

Implementation

C++

int isqrt_newton(int n) {
  int x = 1;
  bool decreased = false;
  for (;;) {
    int nx = (x + n / x) >> 1;
    if (x == nx || (nx > x && decreased)) break;
    decreased = nx < x;
    x = nx;
  }
  return x;
}

Python

def isqrt_newton(n):
    x = 1
    decreased = False
    while True:
        nx = (x + n // x) // 2
        if x == nx or (nx > x and decreased):
            break
        decreased = nx < x
        x = nx
    return x

High-Precision Square Roots

Finally, consider Newton iteration with high-precision numbers. The iteration itself is unchanged, but this time we need to pay attention to the initial approximate solution, namely the value of $x_0$. Since the numbers for which high precision is needed are usually very large, different initial values can have a significant effect on the algorithm's efficiency. A natural idea is to choose $x_0=2^{\left\lfloor\frac{1}{2}\log_2n\right\rfloor}$, which makes $x_0$ quick to compute and gives a reasonably close approximation to the square root.

Implementation

Here is an implementation in Java:

public static BigInteger isqrtNewton(BigInteger n) {
  BigInteger a = BigInteger.ONE.shiftLeft(n.bitLength() / 2);
  boolean p_dec = false;
  for (;;) {
    BigInteger b = n.divide(a).add(a).shiftRight(1);
    if (a.compareTo(b) == 0 || a.compareTo(b) < 0 && p_dec)
      break;
    p_dec = a.compareTo(b) > 0;
    a = b;
  }
  return a;
}

In practice, when $n=10^{1000}$, this algorithm takes 60 ms to run. If we do not optimize the value of $x_0$ and start the iteration directly from $x_0=1$, the running time increases to 120 ms.

Exercises

• UVa 10428 - The Roots

• LeetCode 69. Square Root of x

  This page is mainly translated from the blog post Newton's method (tangent) for finding roots and its English translation Newton's method for finding roots. The Russian version is licensed under Public Domain + Leave a Link; the English version is licensed under CC-BY-SA 4.0.