Write a function or complete program that takes a positive number n and performs n steps of an iterative algorithm for calculating π that has quadratic convergence (i.e. it approximately doubles the number of accurate digits at every iteration) then returns or prints out 2n correct digits (including the beginning 3). One such algorithm is the Gauss–Legendre algorithm, but you are free to use a different algorithm if you prefer.


input 1 → output 3.1
input 2 → output 3.141
input 5 → output 3.1415926535897932384626433832795


  • Each iteration of the algorithm must perform a constant number of basic operations such as addition, subtraction, multiplication, division, power and root (with integer exponent/degree) - each such operation on "big" integer/decimal numbers is counted as one even if it involves one or more loops internally. To be clear, trigonometric functions and powers involving complex numbers are not basic operations.
  • The algorithm is expected to have an initialization step which must also have a constant number of operations.
  • If the algorithm needs 1 or 2 more iterations to get to 2n correct digits, you can perform up to n+2 iterations instead of just n.
  • If it wasn't clear enough, after the correct 2n digits, your program must not print anything else (such as more correct digits, wrong digits or the complete works of Shakespeare).
  • Your program must support values of n from 1 to at least 20.
  • Your program should not take more than an hour for n=20 on a modern computer (not a hard rule, but try to keep it reasonable).
  • The program must not obtain more than 20 accurate digits after the initialization and first iteration of the algorithm.
  • The program must be runnable in Linux using freely available software.
  • The source code must use only ASCII characters.


Straightforward code golf, shortest code wins.


The winner is Digital Trauma, I finally finished running his code on n=20 (just kidding). Special prize goes to primo for his very fast python solution and different algorithm :)

  • 1
    \$\begingroup\$ Quadratic convergence is error ~ N^(1/2). What you describe is exponential convergence error ~ 2^(-N). \$\endgroup\$
    – yo'
    Commented Mar 18, 2015 at 11:26
  • \$\begingroup\$ @yo' are you saying that wikipedia is wrong? \$\endgroup\$ Commented Mar 18, 2015 at 19:55
  • 1
    \$\begingroup\$ Misleading, at least to say: "quadratic convergence" is ~q^(n^2) according to the 1st section there and ~q^2 according to the 2nd section there. \$\endgroup\$
    – yo'
    Commented Mar 18, 2015 at 20:03
  • 1
    \$\begingroup\$ I don't understand codegolf: surely anyone could just write their own programming language specifically for a single task like this, then write a program of, say, 0 bytes? \$\endgroup\$
    – minseong
    Commented Mar 18, 2015 at 20:52
  • 2
    \$\begingroup\$ @theonlygusti that would be a standard loophole and would get disqualified \$\endgroup\$ Commented Mar 18, 2015 at 20:54

5 Answers 5


dc, 99 bytes



With whitespace and comments for "readability":

2?dsi               # Push 2. push input n, duplicate and store in i
1+^k                # Set calculation precision to 2^(n+1)
1dddsa              # Push four 1s. Store 1st in a
2v/sb               # Store 1/sqrt(2) in b
4/st                # Store 1/4 in t
sp                  # Store 1 in p
[                   # Start iteration loop macro
lalb*v              # Save sqrt(a*b) on stack
lalb+2/d            # Save a[i+1] = (a[i]+b[i])/2 on stack and duplicate
la-d*lp*ltr-        # Save t-p(a[i]-a[i+1])^2 on the stack
st                  # Store t result from stack
sa                  # Store a result from stack
sb                  # Store b result from stack
lp2*sp              # Store 2p in p
li1-dsi0<m]         # Decrement iteration counter i; recurse into macro if < 0
dsmx                # Duplicate, store and run macro
K2/1-k              # Set display precision to 2^n-1
lalb+d*4lt*/        # Save (a+b)^2/4t on stack
p                   # Print result

dc needs to be told how many digits of precision should be used. The calculation precision needs to be higher than the final display precision, so the calculation precision is set to be 2^(n+1) digits. I have verified the accuracy of the output with n=10 against http://www.angio.net/pi/digits/pi1000000.txt.

This slows down dramatically for larger n; n=12 takes 1.5 mins on my VM. Running a few different samples shows the time complexity is O(e^n) (not surprising). Extrapolating this to n=20 would have a runtime of 233 days. Oh well. Better than heat-death-of-the-universe at least.

This is moderately golfed - the stack is used to eliminate temporary variables during the calculations of each iteration, but there is possibly more use of the stack to shorten this more.

$ dc glpi.dc <<< 1
$ dc glpi.dc <<< 2
$ dc glpi.dc <<< 5
$ time dc glpi.dc <<< 7

real    0m0.048s
user    0m0.039s
sys 0m0.000s

If you don't like dc wrapping output at 70 chars, you can set the environment variable DC_LINE_LENGTH to 0:

$ DC_LINE_LENGTH=0 dc glpi.dc <<< 8
  • 2
    \$\begingroup\$ Haha, "readability." Doesn't really apply to dc. ;) \$\endgroup\$
    – Alex A.
    Commented Mar 17, 2015 at 18:23
  • \$\begingroup\$ It seems to print a lot more than 32 digits for input 5 \$\endgroup\$ Commented Mar 18, 2015 at 0:12
  • \$\begingroup\$ I added a rule for that, plus another one about running time (but not really strict). I also don't like how your output is split into multiple lines with backslashes, is that a limitation of dc? \$\endgroup\$ Commented Mar 18, 2015 at 1:09
  • \$\begingroup\$ I'm afraid the output is wrong for n=6 \$\endgroup\$ Commented Mar 18, 2015 at 3:35
  • 1
    \$\begingroup\$ Great, and you got it under 100 too :) Could you also post the actual golfed 99-char program with no whitespace/comments? \$\endgroup\$ Commented Mar 18, 2015 at 20:17

R, 156 bytes

Let's get this party started... with the absolute naivest implementation of the Gauss-Legendre algorithm ever.

for(i in 1:scan()){if(i<2){a=p=Rmpfr::mpfr(1,2e6);t=a/4;b=t^t}else{x=(a+b)/2;b=(a*b)^.5;t=t-p*(a-x)^2;a=x;p=2*p};o=(a+b)^2/(4*t)};cat(Rmpfr::format(o,2^i))

Ungolfed + explanation:

# Generate n approximations of pi, where n is read from stdin
for (i in 1:scan()) {

    # Initial values on the first iteration
    if (i < 2) {
        a <- p <- Rmpfr::mpfr(1, 1e7)
        t <- a/4
        b <- t^t
    } else {
        # Compute new values
        x <- (a + b) / 2
        b <- (a*b)^0.5
        t <- t - p*(a - x)^2

        # Store values for next iteration
        a <- x
        p <- 2*p

    # Approximate pi 
    o <- (a + b)^2 / (4*t)

# Print the result with 2^n digits
cat(Rmpfr::format(o, 2^i))

The mpfr() function is part of the Rmpfr package. It creates an mpfr object using the first argument as the value and the second argument as the number of bits of precision. We assign a and p to 1, and by defining t based on a (and b based on t), the mpfr type propogates to all four variables, thereby maintaining precision throughout.

As mentioned, this requires the R package Rmpfr, which is an acronym for "R Multiple Precision Floating point Reliable." The package uses GMP in the background. Unfortunately base R does not have the ability to do high-precision arithmetic, hence the package dependency.

Don't have Rmpfr? No sweat. install.packages("Rmpfr") and all of your dreams will come true.

Notice that 2e6 was specified as the precision. That means we have 2,000,000 bits of precision, which is enough to maintain precision for at least n = 20. (Note: n = 20 takes a long time but less than an hour on my computer.)

The approach here is literally just a regurgitation of the formulas on the Wikipedia page, but hey, we have to start somewhere.

Any input is welcome as always!

I had to rewrite a lot of this but I still have to acknowledge that Peter Taylor helped me knock 70 bytes off of my first score. In the words of DigitalTrauma, "boom."


Python 2, 214 bytes

This challenge presented a good excuse for me to learn the Decimal module. The Decimal numbers have definable precision and have square root support. I have set the precision to a conservative estimate of the accuracy depending on the loop count.


I have updated the program to improve accuracy and speed, at the expense of golfing. By using the decimal sqrt() method and replacing the x**2 usage with x*x, it is now 200 times faster. This means it can now compute 20 loops giving a million digit result in 6.5 hours. The Decimal numbers often have an error in the last digit (caused by operations on the limit of precision), so the program now uses and discards 5 extra digits so only accurate digits are printed.

from decimal import*
for i in[0]*e:f=a;a,b=(a+b)/j,(a*b).sqrt();c=f-a;t-=c*c*p;p+=p
print str(l*l/g/t)[:-5]

Sample run:

$ echo 1 | python min.py 
$ echo 2 | python min.py 
$ echo 3 | python min.py 
$ echo 5 | python min.py 
$ echo 12 | python min.py

The ungolfed code:

from decimal import *
d = Decimal

loops = input()
# this is a conservative estimate for precision increase with each loop:
getcontext().prec = 5 + (1<<loops)

# constants:
one = d(1)
two = d(2)
four = two*two
half = one/two

# initialise:
a = one
b = one / two.sqrt()
t = one / four
p = one

for i in [0]*loops :
    temp = a;
    a, b = (a+b)/two, (a*b).sqrt();
    pterm = temp-a;
    t -= pterm*pterm * p;
    p += p

ab = a+b
print str(ab*ab / four / t)[:-5]
  • 4
    \$\begingroup\$ Heh half = one/two \$\endgroup\$ Commented Mar 17, 2015 at 19:18
  • \$\begingroup\$ It seems that you're not printing the correct number of digits. And I wonder if the slowness is due to the unnecessary use of **. \$\endgroup\$ Commented Mar 18, 2015 at 0:19
  • 1
    \$\begingroup\$ @aditsu, I have reduced the accuracy to expected digit count (but throwing away perfectly good accuracy from an iteration is making my teeth itch). Good suggestion on the ** effect. I found lots of speed by getting rid of them. I can't meet the 20 loops in 1 hour though. Perhaps with pypy or Cython? Hmmm. I will consider that. \$\endgroup\$ Commented Mar 18, 2015 at 9:45
  • \$\begingroup\$ Much better :) For this problem, throwing away good accuracy is less evil than continuing into bad accuracy. The 1 hour limit is based on my cjam/java test code run with java 8. Maybe python doesn't have efficient multiplication/division/sqrt for large numbers (Karatsuba & co)? \$\endgroup\$ Commented Mar 18, 2015 at 20:14
  • \$\begingroup\$ @aditsu: I believe integers have karatsuba (and just that) -- but with 32-bit limb size rather than 64-bit limb size. Who knows about Decimal. \$\endgroup\$
    – user16488
    Commented Mar 21, 2015 at 22:24

Python (2.7) - 131 bytes

from gmpy import*

Update: Now using gmpy rather than gmpy2. For whatever reason, in gmpy2 setting the precision on a single value doesn't propagate to other values. The result of any calculation reverts back to the precision of the current context. Precision does propagate in gmpy, which seems more intuitive to me. It's also considerably less verbose.

Using one of the many algorithms devised by Borwein and Borwein, slightly refactored. n=20 takes about 11 seconds on my box. Not the most efficient method, but still not bad.


The original algorithm was the following:

Refactoring was done incrementally, but the end result is that

The major simplification happens in pn+1. It's also slightly faster due to having eliminated a division.

The current implementation pushes an back one iteration in the calculation of pn+1, which allows for a different initialization of p0 (2√2), but is otherwise identical.

Sample Usage

$ echo 1 | python pi-borwein.py

$ echo 2 | python pi-borwein.py

$ echo 5 | python pi-borwein.py

$ echo 10 | python pi-borwein.py
  • \$\begingroup\$ Great, but you're missing a digit for n=7 \$\endgroup\$ Commented Mar 18, 2015 at 22:42
  • \$\begingroup\$ Also, is it this algorithm? \$\endgroup\$ Commented Mar 18, 2015 at 22:51
  • \$\begingroup\$ @aditsu fixed, and yes it is. \$\endgroup\$
    – primo
    Commented Mar 18, 2015 at 22:52
  • \$\begingroup\$ Now the last digit is wrong for n=5 \$\endgroup\$ Commented Mar 18, 2015 at 22:54
  • 1
    \$\begingroup\$ @aditsu pip install gmpy worked for me; gmpy and gmpy2 are separate packages. However, it does rely on the deprecated distutils. \$\endgroup\$
    – primo
    Commented Mar 25, 2015 at 7:26

bc and Newton's method, 43 bytes

Newton's method for finding zeros of any function converges quadratically and the algorithm is way simpler than for Gauss-Legendre. It basically boils down to:

xnew = xold - f(xold)/f'(xold)

So here goes the according snippet:


A bit more readable:

/* desired number of iterations */
n = 20;

/* starting estimate for pi */
x = 3;

/* set precision to 2^n */
scale = 2^n;

/* perform n iteration steps */
  // f:=sin, f'=cos
  x -= s(x)/c(x)

To test this, run bc -l in a shell and paste the above snippet. Be prepared to wait for a while; n=20 has been running for about 5min now and no end in sight, yet. n=10 takes about 40s.

  • 4
    \$\begingroup\$ Not sure if sine and cosine qualify as "basic operations such as addition, subtraction, multiplication, division and power (including roots)". However, if you rolled your own sine/cosine, that would probably be acceptable. \$\endgroup\$
    – primo
    Commented Mar 20, 2015 at 8:13
  • 1
    \$\begingroup\$ Neat formula, though -- it says that π is a fixed point of f(x) = x - tan(x) \$\endgroup\$
    – Casey Chu
    Commented Mar 22, 2015 at 7:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.