6972

42
+60  Q:

## Fastest way to get value of pi

Solutions welcome in any language. :-) I'm looking for the fastest way to obtain the value of pi, as a personal challenge. More specifically I'm using ways that don't involve using #defined constants like M_PI, or hard-coding the number in.

The program below tests the various ways I know of. The inline assembly version is, in theory, the fastest option, though clearly not portable; I've included it as a baseline to compare the other versions against. In my tests, with built-ins, the 4 * atan(1) version is fastest on GCC 4.2, because it auto-folds the atan(1) into a constant. With -fno-builtin specified, the atan2(0, -1) version is fastest.

Here's the main testing program (pitimes.c):

#include <math.h>
#include <stdio.h>
#include <time.h>

#define ITERS 10000000
#define TESTWITH(x) {                                                       \
diff = 0.0;                                                             \
time1 = clock();                                                        \
for (i = 0; i < ITERS; ++i)                                             \
diff += (x) - M_PI;                                                 \
time2 = clock();                                                        \
printf("%s\t=> %e, time => %f\n", #x, diff, diffclock(time2, time1));   \
}

static inline double
diffclock(clock_t time1, clock_t time0)
{
return (double) (time1 - time0) / CLOCKS_PER_SEC;
}

int
main()
{
int i;
clock_t time1, time2;
double diff;

/* Warmup. The atan2 case catches GCC's atan folding (which would
* optimise the 4 * atan(1) - M_PI'' to a no-op), if -fno-builtin
* is not used. */
TESTWITH(4 * atan(1))
TESTWITH(4 * atan2(1, 1))

#if defined(__GNUC__) && (defined(__i386__) || defined(__amd64__))
extern double fldpi();
TESTWITH(fldpi())
#endif

/* Actual tests start here. */
TESTWITH(atan2(0, -1))
TESTWITH(acos(-1))
TESTWITH(2 * asin(1))
TESTWITH(4 * atan2(1, 1))
TESTWITH(4 * atan(1))

return 0;
}


And the inline assembly stuff (fldpi.c), noting that it will only work for x86 and x64 systems:

double
fldpi()
{
double pi;
asm("fldpi" : "=t" (pi));
return pi;
}


And a build script that builds all the configurations I'm testing (build.sh):

#!/bin/sh
gcc -O3 -Wall -c           -m32 -o fldpi-32.o fldpi.c
gcc -O3 -Wall -c           -m64 -o fldpi-64.o fldpi.c

gcc -O3 -Wall -ffast-math  -m32 -o pitimes1-32 pitimes.c fldpi-32.o
gcc -O3 -Wall              -m32 -o pitimes2-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -fno-builtin -m32 -o pitimes3-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -ffast-math  -m64 -o pitimes1-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall              -m64 -o pitimes2-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall -fno-builtin -m64 -o pitimes3-64 pitimes.c fldpi-64.o -lm


Apart from testing between various compiler flags (I've compared 32-bit against 64-bit too, because the optimisations are different), I've also tried switching the order of the tests around. The atan2(0, -1) version still comes out top every time, though.

I'm keen to hear what results you have, as well as improvements to the testing process. :-)

+20  A:

Here's a general description of a technique for calculating pi that i learnt in high-school.

I only share this because I think it is simple enough that anyone can remember it, indefinitely, plus it teaches you the concept of "Monte-Carlo" methods -- which are statistical methods of arriving at answers that don't immediately appear to be deducible through random processes.

Draw a square, and inscribe a quadrant (one quarter of a semi-circle) inside that square (a quadrant with radius equal to the side of the square, so it fills as much of the square as possible)

Now throw a dart at the square, and record where it lands -- that is, choose a random point anywhere inside the square. Of course it landed inside the square, but is it inside the semi-circle? record this fact.

Repeat this process many times -- and you will find there is a ratio of the number of points inside the semi-circle versus the total number thrown, call this ratio x.

Since the area of the square is r times r, you can deduce that the area of the semi circle is x times r times r (that is, x times r squared). Hence x times 4 will give you pi.

This is not a quick method to use. But it's a nice example of a monte carlo method. And if you look around, you may find that many problems otherwise outside your computational skills can be solved by such methods.

This is the method we used to calculate Pi in a java project in school. Just used a randomizer to come up with the x,y coordinates and the more 'darts' we threw the closer to Pi we came.
+4  A:

I apologize, this is my fault. The original question got pulled because our inappropriate soft-delete threshold was set absurdly low (TWO!)

Now it's set to 10.

Carry on! And be nice this time! Sheesh!

+46  A:

The monte-carlo method, as mentioned, applies some great concepts but it is, clearly, not the fastest --not by a long shot, not by any reasonable usefulness. Also, it all depends on what kind of accuracy you are looking for. The fastest pi I know of is the digits hard coded. Looking at Pi and Pi[PDF], there are a lot of formulas.

Here is a method that converges quickly (~14digits per iteration). The current fastest application, PiFast, uses this formula with the FFT. I'll just write the formula, since the code is straight forward. This formula was almost found by Ramanujan and discovered by Chudnovsky. It is actually how he calculated several billion digits of the number --so it isn't a method to disregard. The formula will overflow quickly since we are dividing factorials, it would be advantageous then to delay such calculating to remove terms. , where,

k1=545140134; k2=13591409; k3=640320; k4=100100025; k5=327843840; k6=53360;


Below is the Brent–Salamin algorithm. Wikipedia mentions that when a and b are 'close enough' then (a+b)^2/4t will be an approximation of pi. I'm not sure what 'close enough' means, but from my tests, one iteration got 2digits, two got 7, and three had 15, of course this is with doubles, so it might have error based on it's representation and the 'true' calculation could be more accurate.

let pi_2 iters =
let rec loop_ a b t p i =
if i = 0 then a,b,t,p
else
let a_n = (a +. b) /. 2.0
and b_n = sqrt (a*.b)
and p_n = 2.0 *. p in
let t_n = t -. (p *. (a -. a_n) *. (a -. a_n)) in
loop_ a_n b_n t_n p_n (i - 1)
in
let a,b,t,p = loop_ (1.0) (1.0 /. (sqrt 2.0)) (1.0/.4.0) (1.0) iters in
(a +. b) *. (a +. b) /. (4.0 *. t)


Lastly, how about some pi golf (800digits)? 160characters!

int a=10000,b,c=2800,d,e,f,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);}

Woot! Very nice. :-)
I know! all I wanted was a nice complete answer. I just thought of things late, and I can't spell right!
Assuming you are trying to implement the first one yourself, wouldn't sqr(k3) be a problem? I'm pretty sure it would end up an irrational number that you will have to estimate (IIRC, all roots that are not whole numbers are irrational). Everything else looks pretty straight-forward if you are using infinite precision arithmetic but that square root is a deal breaker. The second one includes a sqrt as well.
+11  A:

Got a physics engine handy?

http://www.wikihow.com/Calculate-Pi-by-Throwing-Frozen-Hot-Dogs

Granted, this is slightly non-optimal, but it's way more interesting.

Hm. My Trend Micro OfficeScan has that url listed as malicious.
Not sure what to tell you... It seems fine to me.
+1  A:

@Ryan: Just slightly non-optimal and way more interesting, like Bogosort. :-P Very neat!

+1  A:

@Chris: Don't forget Bogosort's cousin, Bozosort.

A:

@Ryan: Can we consider Bozosort and/or Bogosort "sorting by Monte Carlo"? You can have a metric of sorted-ness, based on the proportion of pair-wise comparisons are in the desired order, and tally up these metrics as the simulation goes....

+1  A:

If by fastest you mean fastest to type in the code, here's the golfscript solution:

;''6666,-2%{2+.2/@*\/10.3??2*+}*1000<~\;
A:

Check out Earl F Glynn's computer lab site... specifically the Buffon's Needles page, which describes the maths and theory behind the Monte Carlo simulation used, and where you can download some source code in Delphi.

+7  A:

3.14

3.14159265358979323846That's accurate enough to calculate the area of the universe to an accuracy of the width of an atom. It should work fine for most uses.:)
Not even close. The largest atom, Cs, has a width of about 0.3nm. The universe is 15 billion light years across (a very conservative estimate, I'm told). A light year is 9.46*10^15 meters. You'd need about 12 more digits of Pi to make that calculation.
Bill: I wish I could upvote your comment
@Bill: 3.1415926535897932384626433832795 So this could be enough, couldn't it? :P I'd like to upvote your comment too :P
@Bill @DrJokepu @Andrea Why not upvote my response ;-)
@Andrea: That's only 11 more digits (and I'm not checking accuracy), so you might want to add a few more just to account for smaller atoms and bigger universes. :) @jake: I had you upvoted long ago. :)
3.14159265358 is enough to calculate the circumference of the earth with a precision of a millimeter
@lubos: assuming the Earth is a sphere, which it isn't.
The Earth is an oblate spheroid (or an oblate ellipsoid)
@Andrea and Bill: There are a lot of misconceptions about the size of the universe (perhaps it's related to how mindboggingly big it is? Or how everything is just speculation?). 46.5 billion ly is the current best estimate of the size of the *observable* universe. According to some theories, the size of the *entire* universe can be at least 10^26 times bigger than the observable universe. So, you might still need more than 26 more digits :P. I'll refer you to http://en.wikipedia.org/wiki/Size_of_the_universe.
See, even I got that wrong. 46.5 billion ly is an estimate for the radius of the observable universe. That gives a diameter close to 100 billion ly, so add one more digit, please.
I thought the observable universe was related to the Big Bang, so if that happened 13.7 billion years ago (estimate from Wikipedia) the radius of the observable universe would be 13.7 billion light-years, since anything farther away wouldn't have reached us yet.
@David: no because space itself is expanding.
This should be enough to calculate the size of the universe to the width of an atom: http://www.exploratorium.edu/pi/Pi10-6.html
+19  A:
assert(inIndiana && bill246passed)
TESTWITH(16/5)


Much quicker ;)

LOL ROFL LOL LOL LOL
+7  A:

There's actually a whole book dedicated (amongst other things) to fast methods for the computation of \pi: 'Pi and the AGM', by Jonathan and Peter Borwein (available on Amazon).

I studied the AGM and related algorithms quite a bit: it's quite interesting (though sometimes non-trivial).

Note that to implement most modern algorithms to compute \pi, you will need a multiprecision arithmetic library (GMP is quite a good choice, though it's been a while since I last used it).

The time-complexity of the best algorithms is in O(M(n)log(n)), where M(n) is the time-complexity for the multiplication of two n-bit integers (M(n)=O(n log(n) log(log(n))) using FFT-based algorithms, which are usually needed when computing digits of \pi, and such an algorithm is implemented in GMP).

Note that even though the mathematics behind the algorithms might not be trivial, the algorithms themselves are usually a few lines of pseudo-code, and their implementation is usually very straightforward (if you chose not to write your own multiprecision arithmetic :-) ). Here you can find a sample implementation in Java (and easily recover the algorithm form it!).

+6  A:

The BBP formula allows you to compute the nth digit - in base 2 (or 16) - without having to even bother with the previous n-1 digits first :)

+24  A:

I really like this program, which approximates pi by looking at its own area :-)

IOCCC 1998 : westley.c

Ow. That hurts my brain.
If you replace _ with -F<00||--F-OO-- it should be easier to follow :-)
or, if you replace _ with "if (previous character is '-') { OO--; } F--;"
I knew the guy who wrote it. Very nice guy, quiet, not the sort of person you'd expect to come up with lots of weird C like that.
Gah! That's utterly insane.
it prints 0.25 here -.-
+1 litb, it prints 0.25 for me too
This program was great in 1998, but was broken because modern preprocessors are more liberal with inserting spaces around macro expansions to prevent things like this from working. It is a relic, unfortunately.
Pass --traditional-cpp to *cpp* to get the intended behavior.
+4  A:

Pick a better algorithm.
This one is more work, but converges fast.

+1  A:

## Calculate PI at compile-time with D.

( Copied from DSource.org )

/** Calculate pi at compile time
*
* Compile with dmd -c pi.d
*/
module calcpi;

import meta.math;
import meta.conv;

/** real evaluateSeries!(real x, real metafunction!(real y, int n) term)
*
* Evaluate a power series at compile time.
*
* Given a metafunction of the form
*  real term!(real y, int n),
* which gives the nth term of a convergent series at the point y
* (where the first term is n==1), and a real number x,
* this metafunction calculates the infinite sum at the point x
* by adding terms until the sum doesn't change any more.
*/
template evaluateSeries(real x, alias term, int n=1, real sumsofar=0.0)
{
static if (n>1 && sumsofar == sumsofar + term!(x, n+1)) {
const real evaluateSeries = sumsofar;
} else {
const real evaluateSeries = evaluateSeries!(x, term, n+1, sumsofar + term!(x, n));
}
}

/*** Calculate atan(x) at compile time.
*
* Uses the Maclaurin formula
*  atan(z) = z - z^3/3 + Z^5/5 - Z^7/7 + ...
*/
template atan(real z)
{
const real atan = evaluateSeries!(z, atanTerm);
}

template atanTerm(real x, int n)
{
const real atanTerm =  (n & 1 ? 1 : -1) * pow!(x, 2*n-1)/(2*n-1);
}

/// Machin's formula for pi
/// pi/4 = 4 atan(1/5) - atan(1/239).
pragma(msg, "PI = " ~ fcvt!(4.0 * (4*atan!(1/5.0) - atan!(1/239.0))) );

Unfortunately, tangents are arctangents are based on pi, somewhat invalidating this calculation.
I just copied it from DSource.org
+4  A:

This is a "classic" method, very easy to implement. This implementation, in python (not so fast language) does it:

from math import pi
from time import time

precision = 10**6 # higher value -> higher precision
# lower  value -> higher speed

t = time()

calc = 0
for k in xrange(0, precision):
calc += ((-1)**k) / (2*k+1.)
calc *= 4. # this is just a little optimization

t = time()-t

print "Calculated: %.40f" % calc
print "Costant pi: %.40f" % pi
print "Difference: %.40f" % abs(calc-pi)
print "Time elapsed: %s" % repr(t)


Anyway the fastest way to get a precise as-much-as-you-want value of pi in python is:

from gmpy import pi
print pi(3000) # the rule is the same as
# the precision on the previous code


here is the piece of source for the gmpy pi method, I don't think the code is as much useful as the comment in this case:

static char doc_pi[]="\
pi(n): returns pi with n bits of precision in an mpf object\n\
";

/* This function was originally from netlib, package bmp, by
* Richard P. Brent. Paulo Cesar Pereira de Andrade converted
* it to C and used it in his LISP interpreter.
*
*
*   sets mp pi = 3.14159... to the available precision.
*   uses the gauss-legendre algorithm.
*   this method requires time o(ln(t)m(t)), so it is slower
*   than mppi if m(t) = o(t**2), but would be faster for
*   large t if a faster multiplication algorithm were used
*   for a description of the method, see - multiple-precision
*   zero-finding and the complexity of elementary function
*   evaluation (by r. p. brent), in analytic computational
*   complexity (edited by j. f. traub), academic press, 1976, 151-176.
*   rounding options not implemented, no guard digits used.
*/
static PyObject *
Pygmpy_pi(PyObject *self, PyObject *args)
{
PympfObject *pi;
int precision;
mpf_t r_i2, r_i3, r_i4;
mpf_t ix;

ONE_ARG("pi", "i", &precision);
if(!(pi = Pympf_new(precision))) {
return NULL;
}

mpf_set_si(pi->f, 1);

mpf_init(ix);
mpf_set_ui(ix, 1);

mpf_init2(r_i2, precision);

mpf_init2(r_i3, precision);
mpf_set_d(r_i3, 0.25);

mpf_init2(r_i4, precision);
mpf_set_d(r_i4, 0.5);
mpf_sqrt(r_i4, r_i4);

for (;;) {
mpf_set(r_i2, pi->f);
mpf_div_ui(pi->f, pi->f, 2);
mpf_mul(r_i4, r_i2, r_i4);
mpf_sub(r_i2, pi->f, r_i2);
mpf_mul(r_i2, r_i2, r_i2);
mpf_mul(r_i2, r_i2, ix);
mpf_sub(r_i3, r_i3, r_i2);
mpf_sqrt(r_i4, r_i4);
mpf_mul_ui(ix, ix, 2);
/* Check for convergence */
if (!(mpf_cmp_si(r_i2, 0) &&
mpf_get_prec(r_i2) >= (unsigned)precision)) {
mpf_mul(pi->f, pi->f, r_i4);
mpf_div(pi->f, pi->f, r_i3);
break;
}
}

mpf_clear(ix);
mpf_clear(r_i2);
mpf_clear(r_i3);
mpf_clear(r_i4);

return (PyObject*)pi;
}


EDIT: I had some problem with cut and paste and identation, anyway you can find the source here.

+31  A:
1. Go to Amazon.com;
2. Search for pi book (containing only digits);
3. Filter out books about cooking;
4. Order the book;
5. Wait for delivery;
6. Start typing the digits from the book into your program;
7. List item;
8. Get bored;
9. Order a scanner;
10. Wait for scanner delivery;
11. Realize that you also need an OCR program;
12. Order OCR program, yes you will get a floppy in the mail;
13. Add a PDF feature to your program so that it can read the PDF file containing PI with thousands of decimals;
14. Celebrate with coffee.
Note: You can speed up this algorithm by clearing the "Super Saver Shipping" flag in step 5. Setting this flag does *not* optimize for speed.
Or you could order it on Kindle but then if somebody in Indiana objects to the value amazon will just delete it.
"Filter out books about cooking" - ha ha :-D
... and books about tigers in rowboats (http://www.amazon.com/Life-Pi-Yann-Martel/dp/B0017HV05E/)
+4  A:
<joke>22/7</joke>

Nice, that's a very E4X-compatible answer! :-P
Good grief, 22/7 is merely the upper bound found by Archimedes 2000 years ago! Please, [use 355/113.](http://betterexplained.com/articles/prehistoric-calculus-discovering-pi/) ;)
+1  A:

I once thought that I could write something to do one of these Monte Carlo simulations. I didn't get very far before I discovered that in order to decide if something is inside or outside the circle, you had to describe the circle itself--something for which you need to know the value of pi! So that ended that. Was I wrong about this?

YES! You actually measure the distance of a point to the origin and if that is < 1 then you are within the circle. Distance is sqrt(x*x+y*y), no use of pi involved. This assumes you are generating 2 numbers (x and y) from -1 to 1.
@nlucaroni: And remember that sqrt(1) = 1, so you can skip that. Test if x*x + y*y is less than or equal to 1; if so, it's in the circle.
Thanks, guys! Now I can go worry if the random number generator used pi to make the numbers :)
+14  A:

Easy, the value of Pi is exactly 3, says so in the Bible

:-^

Did some googling on this premise and there is some fascinating debate on the subject, although it seems silly. To assume the Bible says pi is exactly 3 you'd have to assume "circular in shape" = a perfect circle or a cubit (defined as the length one's forearm) is a precise unit of measurement.
This baffles me, it's a light-hearted joke based on the well know quote that the bible says it's 3, and some people mark it down. Like get a sense of humour guys!
You got voted down (not by me) because 1) your "well know[n] quote" is frequently used to disparage the religious, which makes it less than lighthearted; and 2) it isn't even accurate. Hint: read the the whole passage, run the numbers again, and don't leave out the bowl's thickness this time.
Disparage the religious and my upvotes will be with you.
No, Pi is exactly 3.2. See House Bill #246 in the Indiana House of Representatives.
+2  A:

This version (in Delphi) is nothing special, but it is at least faster than the version Nick Hodge posted on his blog :). On my machine, it takes about 16 seconds to do a billion iterations, giving a value of 3.1415926525879 (the accurate part is in bold).

program calcpi;

{\$APPTYPE CONSOLE}

uses
SysUtils;

var
start, finish: TDateTime;

function CalculatePi(iterations: integer): double;
var
numerator, denominator, i: integer;
sum: double;
begin
{
PI may be approximated with this formula:
4 * (1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 .......)
//}
numerator := 1;
denominator := 1;
sum := 0;
for i := 1 to iterations do begin
sum := sum + (numerator/denominator);
denominator := denominator + 2;
numerator := -numerator;
end;
Result := 4 * sum;
end;

begin
try
start := Now;
WriteLn(FloatToStr(CalculatePi(StrToInt(ParamStr(1)))));
finish := Now;
WriteLn('Seconds:' + FormatDateTime('hh:mm:ss.zz',finish-start));
except
on E:Exception do
Writeln(E.Classname, ': ', E.Message);
end;
end.

A:

Take 314000 and divide by 100000.

+2  A:

Just came across this one that should be here for completeness:

calculate PI in Piet

It has the rather nice property that the precision can be improved making the program bigger.

Here's some insight into the language itself

+1  A:

Let's not forget the infamous macro injection attack:

#define fastPi i=ITERS);\
printf("fastPi\t=> 0.000000, time => 0.000000\n");\
return 0;\
(0

TESTWITH(fastPi())


Gives:

fastPi        => 0.000000, time => 0.000000
+2  A:

How I need a drink, alcoholic of course, after the heavy sessions involving quantum mechanics...

Someone who didn't get it edited "for grammar" and broke it. Not only that, but it was wrong all along, using "session" which is 7 letters where pi has 8. I've edited to fix. Also see http://en.wikipedia.org/wiki/Piphilology for more variations on the same theme.
+2  A:

Since the question is for the fastest way to get the value, rather than the fastest way to calculate the value:

Go to http://www.google.com/search?q=value+of+pi and select one or more of the result pages which appear likely to contain a value accurate enough to meet your needs.

If that's not quick enough, you can bypass the google step and go directly to http://www.dbooth.net/internerd/pifinders.cfm to get the value of pi to far more digits than I care to count at the moment.

For greater precision, http://ja0hxv.calico.jp/pai/epivalue.html has it to 100 billion digits, but they're split up into 1000 separate zip files, which would require the additional steps of retrieving, unzipping, and concatenating the files before you had your value.

+1  A:

Back in the old days, with small word sizes and slow or non-existent floating-point operations, we used to do stuff like this:

/* Return approximation of n * PI; n is integer */
#define pi_times(n) (((n) * 22) / 7)


For applications that don't require a lot of precision (video games, for example), this is very fast and is accurate enough.

For more accuracy use 355 / 113. Very accurate for the size of numbers involved.
+9  A:

I'm looking for the fastest way to obtain the value of pi

A reasonably fast algorithm:

curl http://www.google.com/search?q=pi


Note that it only results in 9 significant digits, and substantial work has to be done on the back end to supply more.

However, there's already a web interface, and it may even be accessible through a SOAP or similar API.

Try it out here.

Another option is

curl http://www.wolframalpha.com/input/?i=pi


Which provides more digits, and an optional interface to get even more digits beyond those initially provided. Further, it provides it in alternative forms, such as a continued fraction.

You can try this one out immediately as well.

A related one that's a little more fun than google --> http://www.wolframalpha.com/input/?i=pi
Thanks! Edited.
would've given you an upvote if I had one
+2  A:

Pi is exactly 3! [Prof. Frink (Simpsons)]

Joke, but here's one in C# (.NET-Framework required).

using System;
using System.Text;

class Program {
static void Main(string[] args) {
int Digits = 100;

BigNumber x = new BigNumber(Digits);
BigNumber y = new BigNumber(Digits);
x.ArcTan(16, 5);
y.ArcTan(4, 239);
x.Subtract(y);
string pi = x.ToString();
Console.WriteLine(pi);
}
}

public class BigNumber {
private UInt32[] number;
private int size;
private int maxDigits;

public BigNumber(int maxDigits) {
this.maxDigits = maxDigits;
this.size = (int)Math.Ceiling((float)maxDigits * 0.104) + 2;
number = new UInt32[size];
}
public BigNumber(int maxDigits, UInt32 intPart)
: this(maxDigits) {
number = intPart;
for (int i = 1; i < size; i++) {
number[i] = 0;
}
}
private void VerifySameSize(BigNumber value) {
if (Object.ReferenceEquals(this, value))
throw new Exception("BigNumbers cannot operate on themselves");
if (value.size != this.size)
throw new Exception("BigNumbers must have the same size");
}

VerifySameSize(value);

int index = size - 1;
while (index >= 0 && value.number[index] == 0)
index--;

UInt32 carry = 0;
while (index >= 0) {
UInt64 result = (UInt64)number[index] +
value.number[index] + carry;
number[index] = (UInt32)result;
if (result >= 0x100000000U)
carry = 1;
else
carry = 0;
index--;
}
}
public void Subtract(BigNumber value) {
VerifySameSize(value);

int index = size - 1;
while (index >= 0 && value.number[index] == 0)
index--;

UInt32 borrow = 0;
while (index >= 0) {
UInt64 result = 0x100000000U + (UInt64)number[index] -
value.number[index] - borrow;
number[index] = (UInt32)result;
if (result >= 0x100000000U)
borrow = 0;
else
borrow = 1;
index--;
}
}
public void Multiply(UInt32 value) {
int index = size - 1;
while (index >= 0 && number[index] == 0)
index--;

UInt32 carry = 0;
while (index >= 0) {
UInt64 result = (UInt64)number[index] * value + carry;
number[index] = (UInt32)result;
carry = (UInt32)(result >> 32);
index--;
}
}
public void Divide(UInt32 value) {
int index = 0;
while (index < size && number[index] == 0)
index++;

UInt32 carry = 0;
while (index < size) {
UInt64 result = number[index] + ((UInt64)carry << 32);
number[index] = (UInt32)(result / (UInt64)value);
carry = (UInt32)(result % (UInt64)value);
index++;
}
}
public void Assign(BigNumber value) {
VerifySameSize(value);
for (int i = 0; i < size; i++) {
number[i] = value.number[i];
}
}

public override string ToString() {
BigNumber temp = new BigNumber(maxDigits);
temp.Assign(this);

StringBuilder sb = new StringBuilder();
sb.Append(temp.number);
sb.Append(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.CurrencyDecimalSeparator);

int digitCount = 0;
while (digitCount < maxDigits) {
temp.number = 0;
temp.Multiply(100000);
sb.AppendFormat("{0:D5}", temp.number);
digitCount += 5;
}

return sb.ToString();
}
public bool IsZero() {
foreach (UInt32 item in number) {
if (item != 0)
return false;
}
return true;
}

public void ArcTan(UInt32 multiplicand, UInt32 reciprocal) {
BigNumber X = new BigNumber(maxDigits, multiplicand);
X.Divide(reciprocal);
reciprocal *= reciprocal;

this.Assign(X);

BigNumber term = new BigNumber(maxDigits);
UInt32 divisor = 1;
bool subtractTerm = true;
while (true) {
X.Divide(reciprocal);
term.Assign(X);
divisor += 2;
term.Divide(divisor);
if (term.IsZero())
break;

if (subtractTerm)
this.Subtract(term);
else
subtractTerm = !subtractTerm;
}
}
}

+4  A:

instead of defining pi as a constant, I always use acos(-1).

cos(-1), or acos(-1)? :-P That (the latter) is one of the test cases in my original code. It's among my preferred (along with atan2(0, -1), which really is the same as acos(-1), except that acos is usually implemented in terms of atan2), but some compilers optimise for 4 * atan(1)!
+1  A:

Brent's method posted above by Chris is very good; Brent generally is a giant in the field of arbitrary-precision arithmetic.

If all you want is the Nth digit, the famous BBP formula is useful in hex

The Brent method wasn't posted by me; it was posted by Andrea, and I just happened to be the last person who edited the post. :-) But I agree, that post deserves an upvote.
+1  A:

If you are willing to use an approximation, 355 / 113 is good for 6 decimal digits, and has the added advantage of being usable with integer expressions. That's not as important these days, as "floating point math co-processor" ceased to have any meaning, but it was quite important once.

+1  A:

Pi is irrational. In any language, the value of pi is precisely the mathematical constant π. You can't get its value any more accurately than that.

Except, it's not "precisely" π, if you're talking "any language" from the point of view of programming languages. Rather, it's a floating-approximation of π, which is as close as floating-point can get, but the difference still gets propagated in future calculations. (See http://stackoverflow.com/questions/1053 for more on this topic.)
You're talking about a rational approximation of pi, not the value of pi. The value of pi is just pi, and if you keep it that way, then the precision gets propagated to future calculations. Multiply by 2 and you get the expression 2π, which is precise. Divide it by 3 and you get π/3, which is also precise. Prematurely approximating to a rational number with floating point representation is a mistake.
+1  A:

Uhh ...

#define PI (3.141592653589793238464)


If you need more digits there are complicated algorithms for producing them as have been posted here. But in general no applications really need that.

As stated in the question: "More specifically I'm using ways that don't involve using #defined constants like M_PI, **or hard-coding the number in**."
+2  A:

Since PI is transcendental, surely all algorithms take forever to generate the value of PI?

:-P

Yes, under an eager evaluation system. :-P Although, most people only need about a double-precision float's worth of a pi, so, when you evaluate under a lazy evaluation system, you just stop evaluating after the right number of digits have been obtained. :-P
Thankfully, someone saying that you cannot calculate PI. You can only approximate it.
+7  A:

In the interests of completeness, a C++ template version, which in an optimised build will compute PI at compile time and will inline to a single value.

#include <iostream>

template<int I>
struct sign
{
enum {value = (I % 2) == 0 ? 1 : -1};
};

template<int I, int J>
struct pi_calc
{
inline static double value ()
{
return (pi_calc<I-1, J>::value () + pi_calc<I-1, J+1>::value ()) / 2.0;
}
};

template<int J>
struct pi_calc<0, J>
{
inline static double value ()
{
return (sign<J>::value * 4.0) / (2.0 * J + 1.0) + pi_calc<0, J-1>::value ();
}
};

template<>
struct pi_calc<0, 0>
{
inline static double value ()
{
return 4.0;
}
};

template<int I>
struct pi
{
inline static double value ()
{
return pi_calc<I, I>::value ();
}
};

int main ()
{
std::cout.precision (12);

const double pi_value = pi<10>::value ();

std::cout << "pi ~ " << pi_value << std::endl;

return 0;
}


Note for I > 10, optimised builds can be slow, as can non-optimised runs. I think for 12 iterations there are around 80k calls to value() (without memoization).

Thank you very much. The question needed one of these.
I run this and get "pi ~ 3.14159265383"
Well, that's accurate to 9dp's. Are you objecting to something or just making an observation?
A:

If you want to compute an approximation of the value of π (for some reason), you should try a binary extraction algorithm. Bellard's improvement of BBP gives does PI in O(N^2).

If you want to obtain an approximation of the value of π to do calculations, then:

PI = 3.141592654


Granted, that's only an approximation, and not entirely accurate. It's off by a little more than 0.00000000004102. (four ten-trillionths, about 4/10,000,000,000).

If you want to do math with π, then get yourself a pencil and paper or a computer algebra package, and use π's exact value, π.

If you really want a formula, this one is fun:

## π = -i ln(-1)

Your formula depends on how you define ln in the complex plane. It has to be non-contiguous along one line in the complex plane, and it's quite common for that line to be the negative real axis.
I was just being "cute" with Euler's formula :P
+2  A:

If this article is true, then the algorithm that Bellard has created could be one of the speediest available. He has created pi to 2.7 TRILLION digits using a DESKTOP PC!

...and he has published his work here

Good work Bellard, You are a pioneer!

http://www.theregister.co.uk/2010/01/06/very_long_pi/

Bellard pioneered in many, many ways...first there was LZEXE, quite possibly the first executable compressor (think what UPX does, then flip back in time to the '80s), and of course now, both QEMU and FFMPEG are widely used. Oh, and his IOCCC entry.... :-P
+5  A:

Back in the late '90's, I wrote an applet for my homepage that would calculate the nth digit of pi -- with an algorithm that ran in O(1) time!

Sadly, to achieve this blazing speed, the algorithm had some rather stark limitations, as I explained on another page that I strongly encouraged visitors to take a look at. It was non-deterministic, and only had a 10% chance of returning the correct answer.

The hate mail that thing generated still makes me smile....

+2  A:

Is there a practical application to calculating pi at runtime instead of hard-coding it?

Like, say, to write code that is portable to alternate universes?

Yes, it means someone auditing your code doesn't actually have to check all the digits, besides making it portable to other precision levels (e.g., if the program is ported to use 128-bit floating point, or decimal64, or whatever).
I see. That's a good point. Of course the ideal thing would be to have the value calculated at compile time, but most compilers can't do that.
@chris: what if you made a mistake in your algo - particularly if you use floating point it is very tricky? It is more likely you have to change your algo. than "lookup" value, when you change the machine. IMHO, any practical engineer should just look up. That said your question is good one, but practically of no use to most us, other than guys with supercomputers.
@Fakrudeen: I didn't say you have to use a very complicated algorithm---acos(-1), atan2(0, -1), and 4 * atan(1) are all very well-known ways to get pi in floating point. My question really is: which one of those, or perhaps some other one, is preferable? For example, in GCC, atan(1) is constant-folded to become pi/4, so that's probably considered preferable if you're programming with GCC.
+1  A:

with doubles:

4.0 * (4.0 * Math.Atan(0.2) - Math.Atan(1.0 / 239.0))
`

this will be accurate up to 14 decimal places, enough to fill a double
(the inaccuracy is probably because the rest of the decimals in the arc tangents are truncated)

also seth, it's

3.141592653589793238463, not 64