views:

6972

answers:

42

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.

Leon Bambrick
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.
Jeff Keslinke
+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!

Jeff Atwood
+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.

alt text, 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[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);}
nlucaroni
Woot! Very nice. :-)
Chris Jester-Young
Oh no! All those edits made your post community-moded. :-(
Chris Jester-Young
I know! all I wanted was a nice complete answer. I just thought of things late, and I can't spell right!
nlucaroni
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.
Bill K
+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.

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

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

Chris Jester-Young
+1  A: 

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

Ryan Fox
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....

Chris Jester-Young
+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<~\;
Michiel de Mare
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.

Drew Gibson
+7  A: 

Read this:

3.14
jakemcgraw
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.:)
Matthew Schinckel
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 the Lizard
Bill: I wish I could upvote your comment
DrJokepu
@Bill: 3.1415926535897932384626433832795 So this could be enough, couldn't it? :P I'd like to upvote your comment too :P
Andrea Ambu
@Bill @DrJokepu @Andrea Why not upvote my response ;-)
jakemcgraw
@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. :)
Bill the Lizard
3.14159265358 is enough to calculate the circumference of the earth with a precision of a millimeter
lubos hasko
@lubos: assuming the Earth is a sphere, which it isn't.
Graeme Perrow
The Earth is an oblate spheroid (or an oblate ellipsoid)
Mitch Wheat
@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.
Martinho Fernandes
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.
Martinho Fernandes
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 Thornley
@David: no because space itself is expanding.
jon hanson
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
Chinmay Kanchi
+19  A: 
assert(inIndiana && bill246passed)
TESTWITH(16/5)

Much quicker ;)

LOL ROFL LOL LOL LOL
PEZ
omg! awesome! LOL!... I never knew about this before!
Lazer
+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!).

OysterD
+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 :)

Tyler
+24  A: 

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

IOCCC 1998 : westley.c

Pat
Ow. That hurts my brain.
Herms
If you replace _ with -F<00||--F-OO-- it should be easier to follow :-)
Pat
or, if you replace _ with "if (previous character is '-') { OO--; } F--;"
FryGuy
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.
David Thornley
Gah! That's utterly insane.
Marcus Downing
it prints 0.25 here -.-
Johannes Schaub - litb
+1 litb, it prints 0.25 for me too
flybywire
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.
Chris Lutz
Pass `--traditional-cpp` to *cpp* to get the intended behavior.
Cirno de Bergerac
+4  A: 

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

EvilTeach
+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))) );
Brad Gilbert
Unfortunately, tangents are arctangents are based on pi, somewhat invalidating this calculation.
Grant Johnson
I just copied it from DSource.org
Brad Gilbert
+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)

You can find more information here.

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.
 *
 * Original comments:
 * 
 *   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
 *   (see comments in mpmul).
 *   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_add(pi->f, pi->f, r_i4);
        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.

Andrea Ambu
+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.
Hugo
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.
Bill the Lizard
Or you could order it on Kindle but then if somebody in Indiana objects to the value amazon will just delete it.
Martin Beckett
"Filter out books about cooking" - ha ha :-D
Tor Haugen
... and books about tigers in rowboats (http://www.amazon.com/Life-Pi-Yann-Martel/dp/B0017HV05E/)
Dour High Arch
+4  A: 
<joke>22/7</joke>
John Topley
Nice, that's a very E4X-compatible answer! :-P
Chris Jester-Young
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/) ;)
Daenyth
+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?

gbarry
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
@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.
David Thornley
Thanks, guys! Now I can go worry if the random number generator used pi to make the numbers :)
gbarry
+14  A: 

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

:-^

Cruachan
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.
JohnFx
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!
Cruachan
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.
Iceman
Disparage the religious and my upvotes will be with you.
innaM
No, Pi is exactly 3.2. See House Bill #246 in the Indiana House of Representatives.
Brian
+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.
JosephStyons
A: 

Take 314000 and divide by 100000.

Andrew G. Johnson
+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

krusty.ar
+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
Ates Goral
+2  A: 

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

Vardhan Varma
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.
Daniel Pryden
+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.

Dave Sherohman
+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.

Kristopher Johnson
For more accuracy use `355 / 113`. Very accurate for the size of numbers involved.
David Thornley
+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.

Adam Davis
A related one that's a little more fun than google --> http://www.wolframalpha.com/input/?i=pi
Brian Knoblauch
Thanks! Edited.
Adam Davis
would've given you an upvote if I had one
mikek3332002
+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[0] = 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");
    }

    public void Add(BigNumber value) {
        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[0]);
        sb.Append(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.CurrencyDecimalSeparator);

        int digitCount = 0;
        while (digitCount < maxDigits) {
            temp.number[0] = 0;
            temp.Multiply(100000);
            sb.AppendFormat("{0:D5}", temp.number[0]);
            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
                this.Add(term);
            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)!
Chris Jester-Young
+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

James Youngman
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.
Chris Jester-Young
+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.

Daniel
+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.

Apocalisp
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.)
Chris Jester-Young
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.
Apocalisp
+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.

Paul Hsieh
As stated in the question: "More specifically I'm using ways that don't involve using `#define`d constants like `M_PI`, **or hard-coding the number in**."
Chris Jester-Young
+2  A: 

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

:-P

jon hanson
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
Chris Jester-Young
Thankfully, someone saying that you cannot calculate PI. You can only approximate it.
dmeister
+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).

jon hanson
Thank you very much. The question needed one of these.
David Thornley
I run this and get "pi ~ 3.14159265383"
maxwellb
Well, that's accurate to 9dp's. Are you objecting to something or just making an observation?
jon hanson
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)

Seth
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.
erikkallen
I was just being "cute" with Euler's formula :P
Seth
+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/

Mark Cooper
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
Chris Jester-Young
+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....

BlairHippo
+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?

Deadcode
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).
Chris Jester-Young
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.
Deadcode
@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
@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.
Chris Jester-Young
+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

qwerty01