Calculate π with quadratic convergence 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.

Examples:

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

Requirements:


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.


Scoring:

Straightforward code golf, shortest code wins.

Winner:

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 :)
dc, 99 bytes
Golfed:
2?dsi1+^k1dddsa2v/sb4/stsp[lalb*vlalb+2/dla-d*lp*ltr-stsasblp2*spli1-dsi0&lt;m]dsmxK2/1-klalb+d*4lt*/p

With whitespace and comments for &quot;readability&quot;:
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&lt;m]         # Decrement iteration counter i; recurse into macro if &lt; 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 &lt;&lt;&lt; 1
3.1
$ dc glpi.dc &lt;&lt;&lt; 2
3.141
$ dc glpi.dc &lt;&lt;&lt; 5
3.1415926535897932384626433832795
$ time dc glpi.dc &lt;&lt;&lt; 7
3.1415926535897932384626433832795028841971693993751058209749445923078\
164062862089986280348253421170679821480865132823066470938446

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 &lt;&lt;&lt; 8
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648
$