What is the fastest way to get the value of π?

Question

I'm looking for the fastest way to obtain the value of π, as a personal challenge. More specifically, I'm using ways that don't involve using #define 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 against the other versions. 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) that 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 optimizations are different), I've also tried switching the order of the tests around. But still, the atan2(0, -1) version still comes out on top every time.

Solution 1

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 measure. Also, it all depends on what kind of accuracy you are looking for. The fastest π I know of is the one with the digits hard coded. Looking at Pi and Pi[PDF], there are a lot of formulae.

Here is a method that converges quickly — about 14 digits per iteration. PiFast, the current fastest application, uses this formula with the FFT. I'll just write the formula, since the code is straightforward. 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 and, since we are dividing factorials, it would be advantageous then to delay such calculations to remove terms.

where,

Below is the Brent–Salamin algorithm. Wikipedia mentions that when a and b are "close enough" then (a + b)² / 4t will be an approximation of π. I'm not sure what "close enough" means, but from my tests, one iteration got 2 digits, two got 7, and three had 15, of course this is with doubles, so it might have an error based on its 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 (800 digits)? 160 characters!

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);}

Solution 2

I really like this program, because it approximates π by looking at its own area.

IOCCC 1988 : westley.c

#define _ -F<00||--F-OO--;
int F=00,OO=00;main(){F_OO();printf("%1.3f\n",4.*-F/OO/OO);}F_OO()
{
            _-_-_-_
       _-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
        _-_-_-_-_-_-_-_
            _-_-_-_
}

Solution 3

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.

Solution 4

In the interests of completeness, a C++ template version, which, for an optimised build, will compute an approximation of 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, likewise for non-optimised runs. For 12 iterations I believe there are around 80k calls to value() (in the absence of memoisation).

Solution 5

The following answers precisely how to do this in the fastest possible way -- with the least computing effort. Even if you don't like the answer, you have to admit that it is indeed the fastest way to get the value of PI.

The FASTEST way to get the value of Pi is:

1. chose your favourite programming language
2. load its Math library
3. and find that Pi is already defined there -- ready for use!

In case you don't have a Math library at hand..

The SECOND FASTEST way (more universal solution) is:

look up Pi on the Internet, e.g. here:

http://www.eveandersson.com/pi/digits/1000000 (1 million digits .. what's your floating point precision? )

or here:

http://3.141592653589793238462643383279502884197169399375105820974944592.com/

or here:

http://en.wikipedia.org/wiki/Pi

It's really fast to find the digits you need for whatever precision arithmetic you would like to use, and by defining a constant, you can make sure that you don't waste precious CPU time.

Not only is this a partly humorous answer, but in reality, if anybody would go ahead and compute the value of Pi in a real application .. that would be a pretty big waste of CPU time, wouldn't it? At least I don't see a real application for trying to re-compute this.

Also consider that NASA only uses 15 digits of Pi for calculating interplanetary travel:

• TL;DR: https://twitter.com/Rainmaker1973/status/1463477499434835968
• JPL Explanation: https://www.jpl.nasa.gov/edu/news/2016/3/16/how-many-decimals-of-pi-do-we-really-need/

Dear Moderator: please note that the OP asked: "Fastest Way to get the value of PI"

View more solutions

Comments

[krusty.ar]

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

[erikkallen]

Why do you consider using atan(1) different from using M_PI? I'd understand why you want to do this if you only used arithmetic operations, but with atan I don't see the point.

[C. K. Young]

@erik: Not all languages have a built-in constant like M_PI. I was trying to find an "authoritative" way to get a (floating-point) value of pi that (in theory) works across a variety of languages (and/or their built-in libraries). My current preferred method is using atan2(0, -1), but perhaps there are better ways.

[Mark Cooper]

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! theregister.co.uk/2010/01/06/very_long_pi

[C. K. Young]

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

[Tilo]

the question is: why would you not want to use a constant? e.g. either defined by a library or by yourself? Computing Pi is a waste of CPU cycles, as this problem has been solved over and over and over again to a number of significant digits much greater than needed for daily computations

[ern0]

There is only one solution which is faster than pre-calculate constant PI: pre-calculate all the values appear in formulas, e.g. when circumference needed, you may pre-calculate 2*PI instead of multiplying every time the PI by 2 in runtime.

[Zeus]

@ChrisJester-Young Nope. That's not the intention.. I recently started reading more into the rules.. and thought I would do my part on the threads that I visit to remind users that might have forgotten coz of the long time frame. I'm in no way trying to be a police here. Apologies if i came across rude.

[C. K. Young]

@Zeus In this specific case, my question was actually intended to be a "fun" micro-optimisation question (which, these days, would probably be a better fit for Programming Puzzles & Code Golf), but the general premise of "fastest way to calculate pi" seemed to be useful enough to keep this question here. So, at some stage, I will probably reevaluate whether I should just accept the best algorithmic answer (probably nlucaroni's one), without regard to whether it's related to micro-optimisation.

[C. K. Young]

@HopelessN00b In the dialect of English I speak, "optimisation" is spelt with an "s", not a "z" (which is pronounced as "zed", BTW, not "zee" ;-)). (This is not the first time I've had to revert this sort of edit, too, if you look at the review history.)

[C. K. Young]

nlucaroni's answer has reached 100 upvotes (congrats), so it's probably a good point to green-tick it. Enjoy! (Though, since it's community wiki and all, it is generating no rep, so not even sure if nlucaroni will even notice this.)

[C. K. Young]

@Pacerier See en.wiktionary.org/wiki/boggle and also en.wiktionary.org/wiki/mindboggling.

[Admin]

This should ideally get the attention of Mysticial, since he's the world record holder of computing pi to the largest number of digits. stackoverflow.com/questions/14283270/…

[signonsridhar]

9801/(1103√8)..gives six decimal places..this is the fastest way to calculate PI? = 3.14159273

[C. K. Young]

@signonsridhar No, we're only talking about computation methods that give the exact same result as M_PI when truncated to double-precision.

[signonsridhar]

@Chris Jester-Young Well I just saw a video on Ramanujan who gave this way to calculate PI. So I just shared it :>

[Admin]

You can start off by learning how the formulas of getting π, then try to write it in a C file. Using the Brent-Salamin Formula, you can try to figure out its algorithm and convert it somehow as a program. I don't know how to do it yet because I don't understand a few things...

[Pat]

If you replace _ with -F<00||--F-OO-- it should be easier to follow :-)

[Grant Johnson]

Unfortunately, tangents are arctangents are based on pi, somewhat invalidating this calculation.

[FryGuy]

or, if you replace _ with "if (previous character is '-') { OO--; } F--;"

[Jeff Keslinke]

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.

[C. K. Young]

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)!

[David Thornley]

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

[visual_learner]

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.

[erikkallen]

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.

[Rahul Das]

Pass --traditional-cpp to cpp to get the intended behavior.

[Beweelam]

I run this and get "pi ~ 3.14159265383"

[jon-hanson]

Well, that's accurate to 9dp's. Are you objecting to something or just making an observation?

[Bill K]

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.

[Stephen]

in my experience, 'close enough' usually means there's a taylor series approximation involved.

[Sebastião Miranda]

what is the name of the algorithm used here to compute PI ?

[jon-hanson]

@sebastião-miranda Leibniz's formula, with an averaging acceleration improve convergence. pi_calc<0, J> calculates each successive term from the formula and the non-specialised pi_calc<I, J> calculates the average.

[Petter Friberg]

@Pat if you wounder why I edited it was because I saw this answer in the LQP queue stackoverflow.com/review/low-quality-posts/16750528, hence to avoid deletion I added the code in the link to the answer.

[Max]

Dear Tilo: please note that the OP said: "I'm looking for the fastest way to obtain the value of π, as a personal challenge. More specifically, I'm using ways that don't involve using #define constants like M_PI, or hard-coding the number in.

[Tilo]

Dear @Max: please note that the OP edited their original question after I answered it - that's hardly my fault ;) My solution is still the fastest way, and solves the problem with any desired floating point precision and no CPU cycles elegantly :)

[Max]

Oh sorry, I didn't realise. Just a thought, wouldn't the hard coded constants have less precision than calculating pi? I guess it depends on what language it is and how willing the creator is to put all of the digits in :-)

[Yin Cognyto]

I realize that you answered this in the most honest and funny way possible, but I also realize that there are many people taking it seriously and using the idea as a way of life - the number of upvotes on this proves it: don't do anything to use your brain, because someone else did, does or will do it for you. After all, folks already send already made birthday wishes to friends from their phone cause they can't come up with a couple of original words to express that...

[Tilo]

vm.tiktok.com/ZTdsuAWQK/?k=1