Exact Integer Square and Cube Roots

Here I present simple algorithms BigIntSqrt and BigIntCbrt for computing exact integer values (truncated towards zero) of square and cube roots of arbitrary integers using only integer arithmetic. These can be used to implement BigInt equivalents of Math.sqrt and Math.cbrt in ECMAScript. Both algorithms are efficient, taking only log(log(n)) operations to compute square or cube roots of n. I'm also including proofs that these algorithms compute the exact results for all inputs:

Runnable code is in roots.js.

The goal is to get BigIntSqrt and BigIntCbrt included in the ECMAScript language under some suitable name in the effort to extend Math library functions such as Math.abs to also cover BigInt (the exact form and whether they'll be overloads or separate functions are to be decided). This is preferable to having users blindly copy snippets they find “on the net”. For example, at the time of this writing the top Google search result for "javascript bigint sqrt" produces a link to a buggy and needlessly slow ECMAScript algorithm that sometimes produces incorrect answers, such when computing the square root of 4.

— Waldemar Horwat — May 2022

Notation

All math operations have their usual mathematical meanings on real numbers.

As is the convention in mathematics, a positive real number or integer means that it's greater than zero; similarly a negative real number or integer means that it's less than zero. The integer or real number zero is neither positive nor negative.
⌊x⌋ is the floor of the real number x (i.e. truncated towards -∞ to an integer):
⌊7.1⌋ = 7, ⌊–3.2⌋ = –4, ⌊5⌋ = 5, ⌊–2⌋ = –2.
[x] is the real number x truncated towards 0 to an integer:
[7.1] = 7, [–3.2] = –3, [5] = 5, [–2] = –2.

We can combine [x] with division to denote integer division truncating towards 0:

is the quotient of x divided by y truncated towards 0 to an integer:
[17/5] = [3.4] = 3, [–7/2] = [–3.5] = –3, [10/2] = 5.
When x and y are integers, this is the same as ECMAScript's BigInt division of x and y.

When x ≥ 0 and y > 0, the result of x/y is nonnegative, so truncating it towards 0 is the same as truncating it towards -∞.

In such nonnegative cases we'll sometimes use instead of . In those cases they're equivalent and both denote ECMAScript's BigInt division of x and y.

ECMAScript Algorithms

Bit-Size

We'll need a helper function BigIntLog2 that, given a positive integer n, returns the position of its most significant set bit when expressed in binary. For example, BigIntLog2(1) = 0, BigIntLog2(2) = 1, BigIntLog2(255) = 7, BigIntLog2(256) = 8.

Mathematically BigIntLog2 is defined as

To put it another way, BigIntLog2(n) is the integer w that satisfies 2^w ≤ n < 2^w+1.

We'll assume we have an efficient ECMAScript implementation of BigIntLog2 that takes a positive BigInt and returns a BigInt greater than or equal to zero.

BigIntSqrt

We implement BigIntSqrt in ECMAScript as follows. For simplicity we assume that the argument n has already been checked to be a BigInt.

function BigIntSqrt(n) {
  if (n < 0n)
    throw RangeError("Square root of negative BigInt");
  if (n === 0n)
    return 0n;
  const w = BigIntLog2(n);  // BigIntLog2 returns a BigInt
  let x = 1n << (w >> 1n);  // x is the initial guess x0 here
  let next = (x + n/x) >> 1n;
  do {
    x = next;
  } while ((next = (x + n/x) >> 1n) < x);
  return x;
}

All numbers are nonnegative, so the right-shifts by 1 are equivalent to dividing by 2.

Examples

If we step through BigIntSqrt(123456n), we get the following values of x, w, and the final next:

w = 16n
x₀ = 256n
x₁ = 369n
x₂ = 351n
next = 351n

and the result is x₂ = 351n.

Occasionally the final next can be greater than the last value of x, as in BigIntSqrt(80n):

w = 6n
x₀ = 8n
x₁ = 9n
x₂ = 8n
next = 9n

where the result is x₂ = 8n.

The algorithm converges rapidly, using log(log(n)) iterations. Let's take the square root of a googol:

BigIntSqrt(10n**100n)

w = 332n
x₀ = 93536104789177786765035829293842113257979682750464n
x₁ = 100223346596432806305328643989950694205293778977332n
x₂ = 100000248862684355295361037376723795016205636438280n
x₃ = 100000000000309662407688436639554331697921065069655n
x₄ = 100000000000000000000000479454033675513027179143720n
x₅ = 100000000000000000000000000000000000000000000001149n
x₆ = 100000000000000000000000000000000000000000000000000n
next = 100000000000000000000000000000000000000000000000000n

with the correct answer x₆ = 10⁵⁰.

Let's figure out the exact values of the first 101 decimal digits of the square root of 2:

BigIntSqrt(2n * 10n**200n)

w = 665n
x₀ = 8749002899132047697490008908470485461412677723572849745703082425639811996797503692894052708092215296n
x₁ = 15804375362388773670902839937288652325225438810014411485002799645096232357466133803201506572647827944n
x₂ = 14229549421664044022085597366571792103664648001192143121732144276497668172208742442918179650436359389n
x₃ = 14142404120358730362124525665217767401060790557815075015694913950869126908913392090386697015371402453n
x₄ = 14142135626279684017183484419478198323023092160355018411871574810862373478578236444214290814824160615n
x₅ = 14142135623730950488246557030817812967256722016462383289181725954672117314556926820527863606311500271n
x₆ = 14142135623730950488016887242096980785698583684684248602001325910890328013825461394333000393858369222n
x₇ = 14142135623730950488016887242096980785696718753769480731766797379907324784621070511468589068012098747n
x₈ = 14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727n
next = 14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727n

The result x₈ is ⌊10¹⁰⁰ √2⌋, so √2 = 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727…

Small value optimization

We can add an optimization that computes BigIntSqrt directly using floating-point arithmetic for n small enough that Math.sqrt on n converted to a Number is accurate to within a few least significant bits. The threshold of what counts as "small enough" varies depending on the implementation accuracy of Math.sqrt, but 2⁴⁴ should be safe for any reasonable implementation.

function BigIntSqrt(n) {
  if (n < 0n)
    throw RangeError("Square root of negative BigInt");
  if (n < 0x100000000000n) {  // 2^44
    return BigInt(Math.floor(Math.sqrt(Number(n) + 0.5)));
  }
  ... rest as before

Adding 0.5 to n avoids potential problems if n is a perfect square of an integer s² but Math.sqrt is slightly inaccurate and produces an approximation that's slightly below s.

BigIntCbrt

We implement BigIntCbrt as follows, which produces the cube root of n truncated towards zero. Once again we assume that the argument n has already been checked to be a BigInt.

function BigIntCbrt(n) {
  if (n < 0n)
    return -BigIntCbrt(-n);
  if (n === 0n)
    return 0n;
  const w = BigIntLog2(n);  // BigIntLog2 returns a BigInt
  let x = 1n << (w / 3n);  // x is the initial guess x0 here
  let next = (2n*x + n/(x*x)) / 3n;
  do {
    x = next;
  } while ((next = (2n*x + n/(x*x)) / 3n) < x);
  return x;
}

Examples

If we step through BigIntCbrt(125n), we get the following values of x, w, and the final next:

w = 6n
x₀ = 4n
x₁ = 5n
next = 5n

and the result is x₁ = 5n.

Let's take the cube root of a googol:

BigIntCbrt(10n**100n)

w = 332n
x₀ = 1298074214633706907132624082305024n
x₁ = 2843626090122429343812008194362090n
x₂ = 2307975193491828189478923582864817n
x₃ = 2164422626317660162834416053331641n
x₄ = 2154480709428036015136058273653333n
x₅ = 2154434691014844364863943029148908n
x₆ = 2154434690031883722207769275611571n
x₇ = 2154434690031883721759293566519350n
next = 2154434690031883721759293566519350n

with the answer x₇ = 2154434690031883721759293566519350, which is the first 34 digits of the decimal representation of the cube root of 10.

Small value optimization

As with BigIntSqrt, we can do an optimization for small values of n:

function BigIntCbrt(n) {
  if (n < 0n)
    return -BigIntCbrt(-n);
  if (n < 0x100000000000n) {  // 2^44
    return BigInt(Math.floor(Math.cbrt(Number(n) + 0.5)));
  }
  ... rest as before

Explanation

Now let's delve into how and why the BigIntSqrt and BigIntCbrt algorithms work.

Newton's Method on Real Numbers

The basic approach of computing the square or cube root of n is based on solving the equation x² – n = 0 or x³ – n = 0 for a real number x. We can do this by starting with an initial guess x₀ and then using a variant of Newton's method to refine it to produce successive approximations x₁, x₂, and so on until we find the desired answer. Let's first take a look at how Newton's method works on real numbers:

Given an approximation x_i to a root of the equation f(x) = 0, Newton's method produces the next approximation

For computing square roots we're looking for roots of f(x) = x² – n so Newton's method becomes

For cube roots we're looking for roots of f(x) = x³ – n, in which case Newton's method becomes

These will produce an infinite series of successively more accurate real number approximations of the square or cube root of n.

Newton's Method on Integers

The standard Newton's method uses real numbers and produces an infinite series of approximations. Let's modify it to use only integer arithmetic to find integer square or cube roots truncated towards 0. Later we'll show that we'll arrive at the exact answer in finitely many operations. A similar algorithm for square roots is described on Wikipedia's entry on integer square roots but without the detailed proof of correctness.

Square Root Algorithm

The algorithm that implements BigIntSqrt(n) where n is a nonnegative integer can be stated mathematically as follows:

If n = 0, then return 0. Otherwise, let

For i = 0, 1, 2, 3, … compute the series

until we find the lowest k > 0 such that x_k+1 ≥ x_k. Return x_k.

Later we will prove that our search for such a k terminates and that x_k satisfies

Cube Root Algorithm

The algorithm that implements BigIntCbrt(n) where n is any integer can be stated mathematically as follows:

If n = 0, then return 0.

If n < 0, then return –BigIntCbrt(–n).

Otherwise n is positive. Let

For i = 0, 1, 2, 3, … compute the series

until we find the lowest k > 0 such that x_k+1 ≥ x_k. Return x_k.

Later we will prove that our search for such a k terminates and that x_k satisfies

Choice of Initial Guess

The initial guess x₀ must be positive and should be close to the result. As proven below, the code presented here and the above formulas pick x₀ to always be within a factor of 2 of the final result, which results in quick convergence. Were we to pick the initial guess x₀ = 1 (as some web sites recommend), the algorithm would still work but converge slowly, for example taking 172 iterations to compute the square root of 10¹⁰⁰ instead of the 6 iterations using the algorithm presented here.

Proofs

Lemmas

Let's start with a few lemmas about the floor function.

Given a real number x, ⌊x⌋ is the unique integer i such that x = i + f and 0 ≤ f < 1. From that we can derive the following lemmas.

Lemma 1

This is obvious from the definition of ⌊x⌋.

Lemma 2

To prove this, let x = i + f where i is an integer and 0 ≤ f < 1. Then let i = jn + k where j and k are integers and 0 ≤ k ≤ n – 1. Then we have 0 ≤ k + f < n and thus

Lemma 3

If all of a, b, c, and n are integers with b ≠ 0 and c > 0, we can eliminate the inner floor in:

Proof: By lemmas 1 and 2 we have

BigIntSqrt Proof

Recall that the series for computing BigIntSqrt(n) when integer n > 0 consists of

Also let's define

Given n ≥ 1, we have s ≥ 1 and x₀ is an integer greater than 0. All subsequent terms of the series are also integers due to the definition of x_i+1. Next we'll show by induction that all terms of the series after the zeroth one (i.e. x_i with i > 0) are greater than or equal to s.

Lower bound on BigIntSqrt series terms

Suppose x_i ≥ 1. We'll show that x_i+1 ≥ s.

The square of any real number is nonnegative, so we have

We can divide both sides by the positive quantity 2x_i and simplify to get

Taking the floor of both sides and then using lemma 3 we get

Thus x_i+1 ≥ s ≥ 1, which completes the induction.

Upper bound on BigIntSqrt series terms

We already know that x_i ≥ s for all i > 0. Now suppose x_i > s for some i. We'll show that x_i+1 < x_i so the series is strictly decreasing as long as terms are greater than s.

x_i and s are integers, so x_i > s implies

By the definition of s, we get

Combining the above two inequalities yields

Squaring both sides (which is valid because the function f(a) = a² is monotonically increasing for positive a) and then adding x_i² to both sides results in

Dividing both sides by the positive value 2x_i we get

Applying lemma 3 yields the upper bound on x_i+1.

Combining the upper bound with the lower bound, we get

There are only finitely many integers between s and x_i so the series must decrease on each step (other than the zeroth because we don't necessarily have x₀ > s) and eventually reach x_k = s for some k. At that point the series cannot decrease further, so we can detect x_k = s by looking for k > 0 and x_k+1 ≥ x_k.

Running time

Now let's show that the worst-case running time is O(log(log(n))) integer operations for positive n.

For each term in the x_i series above let's introduce the concept of a corresponding relative error e_i which indicates how much x_i differs from √n. Define

or equivalently

Zeroth iteration

First let's calculate the possible range of e₀.

Applying lemma 2 to the definition of x₀ have

By the properties of the floor function, we have

Exponentiation with a positive base is monotonically increasing, so we can exponentiate all three formulas while preserving the inequalities

which we can simplify using the properties of exponentiation and the definition of x₀ to

Thus x₀ is always within a factor of 2 of √n.

To compute bounds on e₀ we divide all sides by √n and subtract 1 to get

and once again simplify and substitute the definition of e₀ to get the bounds

Subsequent iterations

Let's compute e_i+1 in terms of e_i. Combining the definition of e_i+1, the recursion formula for x_i+1, and x_i = √n(1+e_i) we have

Thus

Let's name the right side f(e_i).

Its derivative is

so f(e_i) is monotonically increasing for e_i ≥ 0 and monotonically decreasing for –1 < e_i ≤ 0 (f has a pole at e_i = –1).

From the bounds on e₀ we can compute the bounds on f(e₀) and thus e₁.

f(e₀) is monotonically decreasing on that range so

Since e_i+1 ≤ f(e_i), we have

Termination

The algorithm terminates when we get e_i ≤ 0 with i > 0. That's because e_i ≤ 0 implies x_i ≤ √n, but we earlier showed the lower bound x_i ≥ s for i > 0, so x_i is an integer such that

which requires x_i = s.

If we get 0 < e_i ≤ 1/√n with i > 0 then the algorithm must terminate on the next iteration. That's because 0 < e_i implies s < x_i while e_i ≤ 1/√n implies

so x_i is an integer such that

s is an integer so the only possibility is x_i = s + 1, but we know from the lower and upper bounds that in this case

so x_i + 1 = s

Number of iterations

Recall the recursion relation

When e_i is much greater than 1 this series decreases exponentially. For example, if e_i = 2¹⁰⁰, then e_i+1 ≈ 2⁹⁹, e_i+2 ≈ 2⁹⁸, and so on, taking roughly 100 iterations to reach 1. This is what would happen had we chosen a bad initial guess such as x₀ = 1. Fortunately the algorithm presented here does not do that.

When e_i is positive and less than 1, the series decreases doubly-exponentially (the exponential of an exponential). For example, if e_i = 2^–100, then e_i+1 ≈ 2^–201, e_i+2 ≈ 2^–403, and so on. Our algorithm operates in this regime. Let's formalize it.

We showed earlier that

We can prove by induction that, as long as the series hasn't reached the termination condition by going to zero or negative earlier, for i > 0 we have

The base case e₁ is above. For the induction case,

If n = 1, x₀ = x₁ = 1 and we're done immediately. If n > 1, we're guaranteed that e_i < 1/√n for some positive i no higher than log₂(log₂(√n)), and the algorithm must terminate either on that or the subsequent iteration. Thus, including the zeroth and the last iteration, the maximum running time is log₂(log₂(√n)) + 2 iterations, which is O(log(log(n)).

BigIntCbrt Proof

The series for computing BigIntCbrt(n) when integer n > 0 consists of

Also let's define

Given n ≥ 1, we have s ≥ 1 and x₀ is an integer greater than 0. All subsequent terms of the series are also integers due to the definition of x_i+1. Next we'll show by induction that all terms of the series after the zeroth one (i.e. x_i with i > 0) are greater than or equal to s.

Lower bound on BigIntCbrt series terms

Suppose x_i ≥ 1. We'll show that x_i+1 ≥ s.

The left factor is always positive and the right factor is the square of a real number which is always nonnegative, so we have

Expanding the products and collecting terms yields

We can divide both sides by the positive quantity 3x_i² and simplify to get

Taking the floor of both sides and then using lemma 3 we get

Thus x_i+1 ≥ s ≥ 1, which completes the induction.

Upper bound on BigIntCbrt series terms

We already know that x_i ≥ s for all i > 0. Now suppose x_i > s for some i. We'll show that x_i+1 < x_i so the series is strictly decreasing as long as terms are greater than s.

x_i and s are integers, so x_i > s implies

By the definition of s, we get

Combining the above two inequalities yields

Cubing both sides (which is valid because the function f(a) = a³ is monotonically increasing) and then adding 2x_i³ to both sides results in

Dividing both sides by the positive value 3x_i² we get

Applying lemma 3 yields the upper bound on x_i+1.

Combining the upper bound with the lower bound, we get

There are only finitely many integers between s and x_i so the series must decrease on each step (other than the zeroth because we don't necessarily have x₀ > s) and eventually reach x_k = s for some k. At that point the series cannot decrease further, so we can detect x_k = s by looking for k > 0 and x_k+1 ≥ x_k.

QED

Running time

The cube root running time analysis is analogous to the square root case, albeit with more complicated recursion formulas. The algorithm is still O(log(log(n))) with slightly different constants.

Name		Name	Last commit message	Last commit date
Latest commit History 41 Commits
formulas		formulas
LICENSE		LICENSE
README.md		README.md
roots.js		roots.js
scratchpad.js		scratchpad.js

License

waldemarhorwat/integer-roots

Folders and files

Latest commit

History

Repository files navigation

Exact Integer Square and Cube Roots

Notation

ECMAScript Algorithms

Bit-Size

BigIntSqrt

Examples

Small value optimization

BigIntCbrt

Examples

Small value optimization

Explanation

Newton's Method on Real Numbers

Newton's Method on Integers

Square Root Algorithm

Cube Root Algorithm

Choice of Initial Guess

Proofs

Lemmas

Lemma 1

Lemma 2

Lemma 3

BigIntSqrt Proof

Lower bound on BigIntSqrt series terms

Upper bound on BigIntSqrt series terms

Running time

Zeroth iteration

Subsequent iterations

Termination

Number of iterations

BigIntCbrt Proof

Lower bound on BigIntCbrt series terms

Upper bound on BigIntCbrt series terms

Running time

About

Resources

License

Stars

Watchers

Forks

Languages