0

Advertisement

Solutions are welcome in any language. :-) 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`

d 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.

Question author Chris-jester-young | Source

# Answer

1

Advertisement

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,

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

Answer author Nlucaroni

Advertisement

Tickanswer.com is providing the only single recommended solution of the question What is the fastest way to get the value of π? under the categories i.e performance , algorithm , language-agnostic , unix , pi , . Our team of experts filter the best solution for you.