# Computing irrational numbers

I am genuinely curious, how do people compute decimal digits of irrational numbers in general, and $\pi$ or nth roots of integers in particular? How do they reach arbitrary accuracy?

• Algorithms or infinite series are two popular ones. May 23, 2015 at 13:19
• Different problems often require different methods. May 23, 2015 at 13:26
• @jm324354 In the end, all programming languages are essentially equivalent. I believe Mathematica itself is mostly written in C/C++. There's no magic in Mathematica: some programmers actually built up the libraries to handle thousands of thousands of digits of accuracy using C, so of course it is possible. May 23, 2015 at 16:17
• I should mention that most irrational numbers (i.e. all except a countable set) are not computable, one example being Chaitin's constant. May 23, 2015 at 18:07
• @JeffreyL. Any computational process is an algorithm. So saying that people compute the decimal expansions of rational numbers using algorithms is like saying that people solve mathematical problems "using math". May 23, 2015 at 22:14

I applaud your genuine curiosity. There is no general method. To compute the decimal digits for any particular irrational number you start with a definition of the number.

For $\pi$ there are lots of methods known. @GregoryGrant 's answer points to one of the oldest and most famous. Archimedes started with inscribed and circumscribed regular hexagons, then doubled the number of sides several times to increase the accuracy of his estimate. (See http://en.wikipedia.org/wiki/Pi#Polygon_approximation_era).

Searching for "computing digits of pi" will lead you to information about modern methods.

For roots of polynomials (like nth roots of integers) Newton's method is efficient.

For many irrationals, continued fractions are a good bet, as @Arteom notes.

For some particular irrationals, clever tricks (based on deep theory) give you lots of digits quickly. For example, the fractions 3/2, 7/5, 17/12, 41/29, ... are better and better approximations to $\sqrt{2}$. You should be able to guess the pattern.

• what is the pattern? May 25, 2015 at 21:30
• @bnosnehpets I'd like to leave it as a puzzle rather than providing an answer. Here's a hint. I added 41/29 to the list. There's a simple way to find the next numerator and denominator from the current numerator and denominator. (The denominator is a little easier.) May 25, 2015 at 21:55

Here's how Archimedes did it 2,500 years ago. $\pi$ is the circumference of a circle with diameter one. What Archimedes did was to inscribe a polygon with $n$ sides inside the circle. Like this for $n=5$.

He then calculated the circumference of the polygon. This gives an approximation to $\pi$ for each $n$ and as $n$ grows it gets closer and closer to $\pi$.

• Great geometric answer for a lower bound. As I recall, he ended up using a polygon with something like 96 sides. The other perk of this method is that if you circumscribe the circle with another polygon (so that the circle is tangent to the sides of this polygon), you obtain an upper bound on the value of $\pi$. May 24, 2015 at 23:36

You generally need to use arbitrary precision arithmetic to compute large numbers of digits of typical irrational numbers. The exception is oddball things like the Champernowne constant 0.12345678910111213141516… :)

There are various arbitrary precision arithmetic packages available. Internally, they use big arrays or lists to hold groups of the digits of the numbers they manipulate. The computer can do simple arithmetic on numbers represented this way using algorithms similar to those taught in grade school which allow you to do arithmetic on multi-digit numbers given the ability to do arithmetic on single digit numbers.

However, for very large numbers the naive multiplication & division algorithms are rather inefficient, so more subtle techniques are used. Multiplication of huge numbers can be performed efficiently using the Fast Fourier transform. Division can be performed by multiplying by the reciprocal of the divisor; reciprocals can be calculated using Newton's method.

Once you have such algorithms at your disposal it's fairly straightforward to calculate $n$th roots using Newton's method, as Lucian mentioned in the comments.

However, it's not necessary to have a full arbitrary precision arithmetic package to calculate large numbers of digits of some of the well-known irrational numbers. It's trivial in most programming languages to add or subtract large numbers that are stored as an array of digit blocks. It's also quite easy to multiply such numbers by a normal-sized number; division by a normal-sized number is similarly easy to implement. And that's all you need to calculate large numbers of terms of many Taylor series to huge precision.

One somewhat famous example of this technique is this gem written in C by the late Dik T. Winter of CWI, which computes $\pi$ to 800 decimal digits.

int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,
f[b]=d%--g,d/=g--,--b;d*=b);}


Unfortunately, that cryptic little program cannot be extended to yield more (correct) digits, due to overflow problems.

But here's a program I wrote a few years ago that can generate large numbers of digits of $e$; it can calculate 10,000 digits in a couple of seconds. The original version of this program was written in C, this version is in Python 2; as you might guess the C version is significantly faster :). Hopefully, the comments explain the algorithm, even if you aren't familiar with Python. The core idea is to treat the Taylor series of $e$ as a number written in factorial base and to simply convert that number into decimal notation.

#! /usr/bin/env python

"""
ECalc - Big e calculator

Generate digits of e, using factorial base notation
Inspired by Dik T. Winter

Created by PM 2Ring 2003.12.26
Converted from C to Python 2007.12.10
"""

import sys, math

# Global constants
N = 50                            # Default number of digits
W = 5                             # Digits per block
A = 10 ** W
LR2P = .5 * math.log(2.*math.pi)  # Log Root 2 Pi
L10 = math.log(10.)               # Log 10

def invlfact(n):
""" Inverse Stirling's approx for log n factorial, using Newton's method """
x = y = n * L10 - LR2P
for i in xrange(3):
x = (y + x) / math.log(x)
return int(round(x))

def bigE(n):
""" Generate n digits of e using factorial base notation """
kmax = invlfact(n) + 1
print "Calculating e to %d places, using %d cells" % (n, kmax)

# Table for factorial base digits, initialized with series for e-2 = 1/2! + 1/3! + ...
# We ignore first two entries
d = [1] * kmax
print "e ~= 2."

j = 1
while n>0:
# Get the next W digits by firstly multiplying d by A,
# propagating carries back down the array,
# then printing the integer part & removing it.
kmax = invlfact(n)    # Number of cells needed for this loop
c = 0                 # Clear carry
for k in xrange(kmax, 1, -1):
d[k] = d[k] * A + c
c = d[k] // k
d[k] -=  c*k

# Print block & field separator. May need modifying if you make W large
jj = (j%10 and 1  or  j%50 and 2  or  j%200 and 3  or  4) - 1
print "%0*d%s" % (W, c, "\n" * jj),

j += 1; n -= W

def main():
# Get number of digits
n = len(sys.argv) > 1 and int(sys.argv[1]) or N
bigE(n)

if __name__ == '__main__':
main()


According to MathWorld

Around 1966, MIT hacker Eric Jensen wrote a very concise program (requiring less than a page of assembly language) that computed e by converting from factorial base to decimal.

I expect that Eric's algorithm was substantially the same as the one in my code, which I derived from an obfuscated C program written by the late Dik Winter.

• Request a source for getting the Eric Jensen code and other details provided by him for that assembly language program. May 14, 2018 at 23:36
• @jiten Sorry, I have no idea, and a question asking for that information would be off-topic here, but maybe someone on Retrocomputing SE could help. I assume you're just asking this for historical interest, you don't need that info to write your own version. May 16, 2018 at 3:09
• Here's a Python version of the Dik Winter Pi code. It can calculate more digits correctly because Python integers don't overflow. Apr 14, 2021 at 20:30

For some irrational numbers, like $\pi$, there are convenient infinite series that converge to them. So for example $$\sum_{n=1}^\infty\frac{1}{n^2}=\frac{\pi^2}{6}$$ By adding up more and more terms of this series you get closer and closer to $\pi^2/6$. You can take an estimate for some value of $n$, multiply by $6$ and then take the square root. That will give you an approximation of $\pi$. This illustrates the idea, but in practice nobody uses this series because there are series that converge to $\pi$ much faster.

Using inequalities: Let's say $x$ is the irrational number you want determine. Then by adding decimals and evaluating if is less or greater, to obtain the desirable accuracy.

So if $x = \sqrt{3} = 1.732\dots$, by squaring $1.7$ and $1.8$ and comparing with $3$ and repeating the process again:

\begin{align} &(1.7)^2<3<(1.8)^2 \\ \text{and then } &(1.73)^2<3<(1.74)^2 \\ \text{and then } &(1.732)^2<3<(1.733)^2 \end{align}

Continued fractions are commonly used in order to get the decimal value with enough accuracy. For example, the A001203 sequence in the Sloan's online encyclopedia represents $\pi$ in a form of continued fraction (actually, all infinite, rational numbers can be decomposed into a finite continued fraction, whereas irrational numbers cannot).

• This is just trading one problem for another. How do you calculate the continued fraction expansion for $\pi$? (This method does work for numbers like $\sqrt 2$ or $e$, though, because the continued fraction is predictable.) May 24, 2015 at 6:59
• @MarioCarneiro Sure, but at least we can make an approximation.We can start with the reccurence: $a_{0}=\lfloor{x} \rfloor, x_{0}=x-a_{0}$, $a_{n}=\lfloor {\frac{1}{x_{n-1}}} \rfloor, x_{n}=\frac{1}{x_{n-1}}-a_{n}$. A continued fraction $[a_{0}, a_{1}, \ldots, a_{n}, \ldots]$ will give us an approximate value of $x$. Regarding square roots - the expansion of the continued fraction is periodic, so it enables to establish an exact form. May 24, 2015 at 9:11
• Calculating those floors requires increasingly accurate approximations or rational bounds on the given number $x$, and conversely the sequence $a_n$ gives you increasingly accurate approximations to $x$. It's seriously no better than changing a number from one base to another, or trading one convergent series for another. For a "general" irrational number for which the continued fraction shows no patterns, you simply have to store the list of convergents, in which case you may as well just store the sequence of decimal digits to equal effect (and you can't really claim to have "computed" it). May 24, 2015 at 13:29
• @MarioCarneiro - the ordinary continued fraction expansion for pi can be difficult, but there are generalised continued fraction expansions for pi that are quite neat. For instance, see en.wikipedia.org/wiki/… May 25, 2015 at 3:55
• Thanks - that makes sense now. I hope my edit gives the right meaning. May 25, 2015 at 20:28

I don't know why nobody has mentioned the fast algorithms that are actually used for calculating $π$ to extremely high precision. Here are some:

• I believe versions of Mathematica uses a variant (the Chudnovsky algorithm) of the first method you mentioned. May 25, 2015 at 9:55
• @TitoPiezasIII: Interesting. Thanks for mentioning! May 25, 2015 at 10:13

how do people compute decimal digits of irrational numbers in general?

In general, you can't. Almost all real numbers are incomputable, only a countable number of real numbers are computable, $\pi$ and $e$ among them as well as almost every other number you can name.

An example of an incomputable number is Chaitin's Constant.

Because most of the numbers we use are computable, it is easy to fall into the trap that most numbers are computable; this simply isn't so.

There is no general way to compute the decimal expansion for an arbitrary member of the set of irrational numbers, I. :(

--Proof Start--

Here's a wishy washy proof (The statements that most obviously lack support are marked with a "*", so if anyone's interested, you should look into them yourself):

Consider R to be the set of real numbers. I is a proper subset of R, and both are uncountably infinite.* (A set that is uncountably infinite is defined by being "too big" to have a 1-to-1 correspondence with the set of integers, which is know to be countably infinite.)

Now, consider D to be the set of definable numbers, that is, numbers that have definitions. D is at most countable simply because there is a 1-to-1 correspondence between definable numbers and a subset of definitions.* (Informally, a "definition" may be thought of as a finite permutation of a finite alphabet of characters -- think dictionary -- and thus the entire set of definitions would be at most countable.*)

Now, consider C to be the set of computable numbers. Because of the way computers work -- think Turing machine -- in order for a number to be computed, it must be defined, thus, C is a subset of D.*

Since C is a subset of D, and D is at most countable, and I is know to be uncountable, there are infinitely more irrational numbers than computable numbers. Therefore, one could not hope to find a general way to "compute" any property on the set of irrational numbers.

--Proof End--

As for the second question, there is an intersection between the computable numbers and the irrational numbers. Many of these numbers can be defined as infinite sums. One can calculate decimal expansions of arbitrary accuracy for these numbers because each term in the sum is generally decreasing and taking a finite sum of the first few terms can give you an idea of the first few digits of an irrational number (say 3.1415 of pi) while the remaining infinite terms can still be calculated, but are insignificant.

One way to do it which i like the most are continued fractions because they make you feel that you can write an irrational number in rational terms that is $\frac{p}{q}$ form but actually its not a rational form.

A continued fraction is of the form :$$a_0+\cfrac{1}{a_1+\cfrac{1}{a_2+\cfrac{1}{a_3+\cdots}}}$$

So a continued fraction for $\pi$ would be something like this:$$3+\cfrac{1}{7+\cfrac{1}{15+\cfrac{1}{1+\cfrac{1}{292+\cdots}}}}$$ So as you can see if you take the first few terms that is $$3+\frac{1}{7}$$ you get $\frac{22}{7}$ that is the famous value of $\pi$ used in high school. And one more fact which is very nice about these continued fractions are as you consider more and more fractions more accurate your answer starts becoming. There are many such continued fractions for other transcendental and irrational number which are really nice,also these continued fractions give you nice bounds for values.

• Continued fractions are a great way, but not for calculating $\pi$. Because there's no nice way of calculating the coefficients $3, 7, 15, 1, 292, \dots$ May 25, 2015 at 12:20
• @user38034 there are other continued fractions for pi that do have a repeating pattern May 25, 2015 at 13:08

One simple way to compute nth roots of integers is to use Khovanskii's matrix method. A good explanation of this process can be found here.

To compute the $n^{th}$ root of the integer $N$, you can use Newton's iterations for the equation $$x^n-N=0:$$ $$x_{k+1}=x_k-\frac{x_k^n-N}{nx_k^{n-1}}=\left(1-\frac1n\right)x_k+\frac{N}{nx_k^{n-1}}.$$

You will need to implement all four arithmetic operations for long numbers; powers are computed using repeated squarings.

For a good starting approximation, just use the floating-point estimate.

As the iterations are "self correcting", you can increase the accuracy as the computation progresses.