What is the fastest way to get the value of π?
Asked Answered
Y

23

352

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.

Yacano answered 1/8, 2008 at 5:21 Comment(17)
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 itselfRandolph
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.Willamina
@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.Yacano
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_piAndizhan
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.... :-PYacano
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 computationsFiesole
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.Skuld
@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.Edelmiraedelson
@Edelmiraedelson 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.Yacano
@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.)Yacano
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.)Yacano
@Pacerier See en.wiktionary.org/wiki/boggle and also en.wiktionary.org/wiki/mindboggling.Yacano
This should ideally get the attention of Mysticial, since he's the world record holder of computing pi to the largest number of digits. #14283770Aircool
9801/(1103√8)..gives six decimal places..this is the fastest way to calculate PI? = 3.14159273Witchery
@Witchery No, we're only talking about computation methods that give the exact same result as M_PI when truncated to double-precision.Yacano
@Chris Jester-Young Well I just saw a video on Ramanujan who gave this way to calculate PI. So I just shared it :>Witchery
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...Shipmaster
F
223

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.

enter image description here

enter image description here

where,

enter image description here

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);}
Faddish answered 2/8, 2008 at 18:22 Comment(2)
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.Middelburg
in my experience, 'close enough' usually means there's a taylor series approximation involved.Amido
P
127

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()
{
            _-_-_-_
       _-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
        _-_-_-_-_-_-_-_
            _-_-_-_
}
Peppery answered 2/9, 2008 at 13:28 Comment(5)
If you replace _ with -F<00||--F-OO-- it should be easier to follow :-)Peppery
or, if you replace _ with "if (previous character is '-') { OO--; } F--;"Fremd
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.Doi
Pass --traditional-cpp to cpp to get the intended behavior.Efta
@Peppery 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.Gorgeous
M
81

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.

Machinist answered 1/8, 2008 at 13:37 Comment(1)
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.Nailbrush
R
61

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

Rickety answered 22/12, 2009 at 15:40 Comment(4)
I run this and get "pi ~ 3.14159265383"Meilhac
Well, that's accurate to 9dp's. Are you objecting to something or just making an observation?Rickety
what is the name of the algorithm used here to compute PI ?Solidago
@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.Rickety
F
45

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 :-) ).

Frost answered 24/8, 2008 at 17:14 Comment(0)
F
44

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:

Fiesole answered 28/10, 2011 at 1:2 Comment(5)
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.Vazquez
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 :)Fiesole
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 :-)Vazquez
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...Trochal
vm.tiktok.com/ZTdsuAWQK/?k=1Fiesole
O
28

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

Owens answered 29/8, 2008 at 9:22 Comment(0)
D
24

This is a "classic" method, very easy to implement. This implementation in python (not the fastest 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 "Constant 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 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 problems with cut and paste and indentation, you can find the source here.

Drift answered 2/10, 2008 at 21:27 Comment(0)
A
24

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

Anear answered 8/3, 2009 at 3:2 Comment(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)!Yacano
C
21

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

;''6666,-2%{2+.2/@*\/10.3??2*+}*`1000<~\;
Carabineer answered 6/8, 2008 at 22:54 Comment(0)
P
20

Use the Machin-like formula

176 * arctan (1/57) + 28 * arctan (1/239) - 48 * arctan (1/682) + 96 * arctan(1/12943) 

[; \left( 176 \arctan \frac{1}{57} + 28 \arctan \frac{1}{239} - 48 \arctan \frac{1}{682} + 96 \arctan \frac{1}{12943}\right) ;], for you TeX the World people.

Implemented in Scheme, for instance:

(+ (- (+ (* 176 (atan (/ 1 57))) (* 28 (atan (/ 1 239)))) (* 48 (atan (/ 1 682)))) (* 96 (atan (/ 1 12943))))

Procreate answered 5/2, 2011 at 5:26 Comment(0)
C
19

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;
        }
    }
}
Consubstantiate answered 26/2, 2009 at 19:22 Comment(0)
R
19

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.

Roslyn answered 17/9, 2009 at 16:30 Comment(0)
P
18

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.

Pharsalus answered 28/2, 2010 at 3:52 Comment(0)
M
17

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))) );
Mollie answered 17/9, 2008 at 17:49 Comment(1)
Unfortunately, tangents are arctangents are based on pi, somewhat invalidating this calculation.Yoho
C
15

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.
Cantone answered 12/1, 2009 at 18:24 Comment(0)
K
13

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.

Kalinda answered 20/2, 2009 at 21:21 Comment(1)
For more accuracy use 355 / 113. Very accurate for the size of numbers involved.Kalindi
V
13

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)

Varletry answered 22/12, 2009 at 21:13 Comment(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.Willamina
S
2

Basically the C version of paperclip optimizer's answer, and much more simpilified:

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

double calc_PI(int K) {
    static const int A = 545140134;
    static const int B = 13591409;
    static const int D = 640320;
    const double ID3 = 1.0 / ((double) D * (double) D * (double) D);
    double sum = 0.0;
    double b = sqrt(ID3);
    long long int p = 1;
    long long int a = B;
    sum += (double) p * (double) a * b;
    for (int k = 1; k < K; ++k) {
        a += A;
        b *= ID3;
        p *= (6 * k) * (6 * k - 1) * (6 * k - 2) * (6 * k - 3) * (6 * k - 4) * (6 * k - 5);
        p /= (3 * k) * (3 * k - 1) * (3 * k - 2) * k * k * k;
        p = -p;
        sum += (double) p * (double) a * b;
    }
    return 1.0 / (12 * sum);
}

int main() {
    for (int k = 1; k <= 5; ++k) {
        printf("k = %i, PI = %.16f\n", k, calc_PI(k));
    }
}

But for more simplification, this algorithm takes Chudnovsky's formula, which I can fully simplify if you don't really understand the code.

Summary: We will get a number from 1 to 5 and add it in to a function we will use to get PI. Then 3 numbers are given to you: 545140134 (A), 13591409 (B), 640320 (D). Then we will use D as a double multiplying itself 3 times into another double (ID3). We will then take the square root of ID3 into another double (b) and assign 2 numbers: 1 (p), the value of B (a). Take note that C is case-insensitive. Then a double (sum) will be created by multiplying the value's of p, a and b, all in doubles. Then a loop up until the number given for the function will start and add up A's value to a, b's value gets multiplied by ID3, p's value will be multiplied by multiple values that I hope you can understand and also gets divided by multiple values as well. The sum will add up by p, a and b once again and the loop will repeat until the value of the loop's number is greater or equal to 5. Later, the sum is multiplied by 12 and returned by the function giving us the result of PI.

Okay, that was long, but I guess you will understand it...

Shipmaster answered 18/10, 2020 at 7:29 Comment(0)
M
1

Calculating π from circle area :-)

<input id="range" type="range" min="10" max="960" value="10" step="50" oninput="calcPi()">
<br>
<div id="cont"></div>

<script>
function generateCircle(width) {
    var c = width/2;
    var delta = 1.0;
    var str = "";
    var xCount = 0;
    for (var x=0; x <= width; x++) {
        for (var y = 0; y <= width; y++) {
            var d = Math.sqrt((x-c)*(x-c) + (y-c)*(y-c));
            if (d > (width-1)/2) {
                str += '.';
            }
            else {
                xCount++;
                str += 'o';
            }
            str += "&nbsp;" 
        }
        str += "\n";
    }
    var pi = (xCount * 4) / (width * width);
    return [str, pi];
}

function calcPi() {
    var e = document.getElementById("cont");
    var width = document.getElementById("range").value;
    e.innerHTML = "<h4>Generating circle...</h4>";
    setTimeout(function() {
        var circ = generateCircle(width);
        e.innerHTML  = "<pre>" + "π = " + circ[1].toFixed(2) + "\n" + circ[0] +"</pre>";
    }, 200);
}
calcPi();
</script>
Marmalade answered 3/6, 2017 at 17:13 Comment(0)
D
1

The Chudnovsky algorithm is pretty fast if you don't mind performing a square root and a couple inverses. It converges to double precision in just 2 iterations.

/*
    Chudnovsky algorithm for computing PI
*/

#include <iostream>
#include <cmath>
using namespace std;

double calc_PI(int K=2) {

    static const int A = 545140134;
    static const int B = 13591409;
    static const int D = 640320;

    const double ID3 = 1./ (double(D)*double(D)*double(D));

    double sum = 0.;
    double b   = sqrt(ID3);
    long long int p = 1;
    long long int a = B;

    sum += double(p) * double(a)* b;

    // 2 iterations enough for double convergence
    for (int k=1; k<K; ++k) {
        // A*k + B
        a += A;
        // update denominator
        b *= ID3;
        // p = (-1)^k 6k! / 3k! k!^3
        p *= (6*k)*(6*k-1)*(6*k-2)*(6*k-3)*(6*k-4)*(6*k-5);
        p /= (3*k)*(3*k-1)*(3*k-2) * k*k*k;
        p = -p;

        sum += double(p) * double(a)* b;
    }

    return 1./(12*sum);
}

int main() {

    cout.precision(16);
    cout.setf(ios::fixed);

    for (int k=1; k<=5; ++k) cout << "k = " << k << "   PI = " << calc_PI(k) << endl;

    return 0;
}

Results:

k = 1   PI = 3.1415926535897341
k = 2   PI = 3.1415926535897931
k = 3   PI = 3.1415926535897931
k = 4   PI = 3.1415926535897931
k = 5   PI = 3.1415926535897931
Donaugh answered 8/5, 2020 at 0:52 Comment(0)
P
1

I think the value of pi is the ratio between the circumference and radius of the circle.

It can be simply achieved by a regular math calculation

Playsuit answered 18/10, 2021 at 4:21 Comment(0)
C
0

Faster than GMPY2 and MPmath Built-ins: A Billion in 45 minutes:


I tried several ways; Manchin, AGM and Chudnovsky Bros. Chudnovsky with Binary Split was the fastest:
My github : https://github.com/Overboard-code/Pi-Pourri

My Binary Split Chudnovsky is about twice the speed of the builtin gmpy2.const_pi(). MPmath.mp.pi() took 50 minutes for a billion, so it was almost as fast as the Chudnovsky.

I would greatly appreciate performance tips as well. I am not sure my code is perfect. It is 100% accurate (all formula agree to 100 million) but maybe could be faster?

I tried gmpy2.const_pi() to 100 million digits and it took 300 seconds vs. 150 seconds for the Chudnovsky on the same machine. pi.txt and pi2.txt were the same.

I got to a billion digits on my old i7 16GB laptop in less than an hour.

Here is a snippet of the fastest of the 12 methods I tried:

class PiChudnovsky:
    """Version of Chudnovsky Bros using Binary Splitting 
        So far this is the winner for fastest time to a million digits on my older intel i7
    """
    A = mpz(13591409)
    B = mpz(545140134)
    C = mpz(640320)
    D = mpz(426880)
    E = mpz(10005)
    C3_24  = pow(C, mpz(3)) // mpz(24)
    #DIGITS_PER_TERM = math.log(53360 ** 3) / math.log(10)  #=> 14.181647462725476
    DIGITS_PER_TERM = 14.181647462725476
    MMILL = mpz(1000000)

    def __init__(self,ndigits):
        """ Initialization
        :param int ndigits: digits of PI computation
        """
        self.ndigits = ndigits
        self.n      = mpz(self.ndigits // self.DIGITS_PER_TERM + 1)
        self.prec   = mpz((self.ndigits + 1) * LOG2_10)
        self.one_sq = pow(mpz(10),mpz(2 * ndigits))
        self.sqrt_c = isqrt(self.E * self.one_sq)
        self.iters  = mpz(0)
        self.start_time = 0

    def compute(self):
        """ Computation """
        try:
            self.start_time = time.time()
            logging.debug("Starting {} formula to {:,} decimal places"
                .format(name,ndigits) )
            __, q, t = self.__bs(mpz(0), self.n)  # p is just for recursion
            pi = (q * self.D * self.sqrt_c) // t
            logging.debug('{} calulation Done! {:,} iterations and {:.2f} seconds.'
                .format( name, int(self.iters),time.time() - self.start_time))
            get_context().precision= int((self.ndigits+10) * LOG2_10)
            pi_s = pi.digits() # digits() gmpy2 creates a string 
            pi_o = pi_s[:1] + "." + pi_s[1:]
            return pi_o,int(self.iters),time.time() - self.start_time
        except Exception as e:
            print (e.message, e.args)
            raise

    def __bs(self, a, b):
        """ PQT computation by BSA(= Binary Splitting Algorithm)
        :param int a: positive integer
        :param int b: positive integer
        :return list [int p_ab, int q_ab, int t_ab]
        """
        try:
            self.iters += mpz(1)
            if self.iters % self.MMILL  == mpz(0):
                logging.debug('Chudnovsky ... {:,} iterations and {:.2f} seconds.'
                    .format( int(self.iters),time.time() - self.start_time))
            if a + mpz(1) == b:
                if a == mpz(0):
                    p_ab = q_ab = mpz(1)
                else:
                    p_ab = mpz((mpz(6) * a - mpz(5)) * (mpz(2) * a - mpz(1)) * (mpz(6) * a - mpz(1)))
                    q_ab = pow(a,mpz(3)) * self.C3_24
                t_ab = p_ab * (self.A + self.B * a)
                if a & 1:
                    t_ab *= mpz(-1)
            else:
                m = (a + b) // mpz(2)
                p_am, q_am, t_am = self.__bs(a, m)
                p_mb, q_mb, t_mb = self.__bs(m, b)
                p_ab = p_am * p_mb
                q_ab = q_am * q_mb
                t_ab = q_mb * t_am + p_am * t_mb
            return [p_ab, q_ab, t_ab]
        except Exception as e:
            print (e.message, e.args)
            raise

Here is the output of 1,000,000,000 digits less than 45 minutes:

python pi-pourri.py -v -d 1,000,000,000 -a 10 

[INFO] 2022-10-03 09:22:51,860 <module>: MainProcess Computing π to 1,000,000,000 digits.
[DEBUG] 2022-10-03 09:25:00,543 compute: MainProcess Starting   Chudnovsky brothers  1988 
    π = (Q(0, N) / 12T(0, N) + 12AQ(0, N))**(C**(3/2))
 formula to 1,000,000,000 decimal places
[DEBUG] 2022-10-03 09:25:04,995 __bs: MainProcess Chudnovsky ... 1,000,000 iterations and 4.45 seconds.
[DEBUG] 2022-10-03 09:25:10,836 __bs: MainProcess Chudnovsky ... 2,000,000 iterations and 10.29 seconds.
[DEBUG] 2022-10-03 09:25:18,227 __bs: MainProcess Chudnovsky ... 3,000,000 iterations and 17.68 seconds.
[DEBUG] 2022-10-03 09:25:24,512 __bs: MainProcess Chudnovsky ... 4,000,000 iterations and 23.97 seconds.
[DEBUG] 2022-10-03 09:25:35,670 __bs: MainProcess Chudnovsky ... 5,000,000 iterations and 35.13 seconds.
[DEBUG] 2022-10-03 09:25:41,376 __bs: MainProcess Chudnovsky ... 6,000,000 iterations and 40.83 seconds.
[DEBUG] 2022-10-03 09:25:49,238 __bs: MainProcess Chudnovsky ... 7,000,000 iterations and 48.69 seconds.
[DEBUG] 2022-10-03 09:25:55,646 __bs: MainProcess Chudnovsky ... 8,000,000 iterations and 55.10 seconds.
[DEBUG] 2022-10-03 09:26:15,043 __bs: MainProcess Chudnovsky ... 9,000,000 iterations and 74.50 seconds.
[DEBUG] 2022-10-03 09:26:21,437 __bs: MainProcess Chudnovsky ... 10,000,000 iterations and 80.89 seconds.
[DEBUG] 2022-10-03 09:26:26,587 __bs: MainProcess Chudnovsky ... 11,000,000 iterations and 86.04 seconds.
[DEBUG] 2022-10-03 09:26:34,777 __bs: MainProcess Chudnovsky ... 12,000,000 iterations and 94.23 seconds.
[DEBUG] 2022-10-03 09:26:41,231 __bs: MainProcess Chudnovsky ... 13,000,000 iterations and 100.69 seconds.
[DEBUG] 2022-10-03 09:26:52,972 __bs: MainProcess Chudnovsky ... 14,000,000 iterations and 112.43 seconds.
[DEBUG] 2022-10-03 09:26:59,517 __bs: MainProcess Chudnovsky ... 15,000,000 iterations and 118.97 seconds.
[DEBUG] 2022-10-03 09:27:07,932 __bs: MainProcess Chudnovsky ... 16,000,000 iterations and 127.39 seconds.
[DEBUG] 2022-10-03 09:27:14,036 __bs: MainProcess Chudnovsky ... 17,000,000 iterations and 133.49 seconds.
[DEBUG] 2022-10-03 09:27:51,629 __bs: MainProcess Chudnovsky ... 18,000,000 iterations and 171.09 seconds.
[DEBUG] 2022-10-03 09:27:58,176 __bs: MainProcess Chudnovsky ... 19,000,000 iterations and 177.63 seconds.
[DEBUG] 2022-10-03 09:28:06,704 __bs: MainProcess Chudnovsky ... 20,000,000 iterations and 186.16 seconds.
[DEBUG] 2022-10-03 09:28:13,376 __bs: MainProcess Chudnovsky ... 21,000,000 iterations and 192.83 seconds.
[DEBUG] 2022-10-03 09:28:18,737 __bs: MainProcess Chudnovsky ... 22,000,000 iterations and 198.19 seconds.
[DEBUG] 2022-10-03 09:28:31,095 __bs: MainProcess Chudnovsky ... 23,000,000 iterations and 210.55 seconds.
[DEBUG] 2022-10-03 09:28:37,789 __bs: MainProcess Chudnovsky ... 24,000,000 iterations and 217.25 seconds.
[DEBUG] 2022-10-03 09:28:46,171 __bs: MainProcess Chudnovsky ... 25,000,000 iterations and 225.63 seconds.
[DEBUG] 2022-10-03 09:28:52,933 __bs: MainProcess Chudnovsky ... 26,000,000 iterations and 232.39 seconds.
[DEBUG] 2022-10-03 09:29:13,524 __bs: MainProcess Chudnovsky ... 27,000,000 iterations and 252.98 seconds.
[DEBUG] 2022-10-03 09:29:19,676 __bs: MainProcess Chudnovsky ... 28,000,000 iterations and 259.13 seconds.
[DEBUG] 2022-10-03 09:29:28,196 __bs: MainProcess Chudnovsky ... 29,000,000 iterations and 267.65 seconds.
[DEBUG] 2022-10-03 09:29:34,720 __bs: MainProcess Chudnovsky ... 30,000,000 iterations and 274.18 seconds.
[DEBUG] 2022-10-03 09:29:47,075 __bs: MainProcess Chudnovsky ... 31,000,000 iterations and 286.53 seconds.
[DEBUG] 2022-10-03 09:29:53,746 __bs: MainProcess Chudnovsky ... 32,000,000 iterations and 293.20 seconds.
[DEBUG] 2022-10-03 09:29:59,099 __bs: MainProcess Chudnovsky ... 33,000,000 iterations and 298.56 seconds.
[DEBUG] 2022-10-03 09:30:07,511 __bs: MainProcess Chudnovsky ... 34,000,000 iterations and 306.97 seconds.
[DEBUG] 2022-10-03 09:30:14,279 __bs: MainProcess Chudnovsky ... 35,000,000 iterations and 313.74 seconds.
[DEBUG] 2022-10-03 09:31:31,710 __bs: MainProcess Chudnovsky ... 36,000,000 iterations and 391.17 seconds.
[DEBUG] 2022-10-03 09:31:38,454 __bs: MainProcess Chudnovsky ... 37,000,000 iterations and 397.91 seconds.
[DEBUG] 2022-10-03 09:31:46,437 __bs: MainProcess Chudnovsky ... 38,000,000 iterations and 405.89 seconds.
[DEBUG] 2022-10-03 09:31:53,285 __bs: MainProcess Chudnovsky ... 39,000,000 iterations and 412.74 seconds.
[DEBUG] 2022-10-03 09:32:05,602 __bs: MainProcess Chudnovsky ... 40,000,000 iterations and 425.06 seconds.
[DEBUG] 2022-10-03 09:32:12,220 __bs: MainProcess Chudnovsky ... 41,000,000 iterations and 431.68 seconds.
[DEBUG] 2022-10-03 09:32:20,708 __bs: MainProcess Chudnovsky ... 42,000,000 iterations and 440.17 seconds.
[DEBUG] 2022-10-03 09:32:27,552 __bs: MainProcess Chudnovsky ... 43,000,000 iterations and 447.01 seconds.
[DEBUG] 2022-10-03 09:32:32,986 __bs: MainProcess Chudnovsky ... 44,000,000 iterations and 452.44 seconds.
[DEBUG] 2022-10-03 09:32:53,904 __bs: MainProcess Chudnovsky ... 45,000,000 iterations and 473.36 seconds.
[DEBUG] 2022-10-03 09:33:00,832 __bs: MainProcess Chudnovsky ... 46,000,000 iterations and 480.29 seconds.
[DEBUG] 2022-10-03 09:33:09,198 __bs: MainProcess Chudnovsky ... 47,000,000 iterations and 488.66 seconds.
[DEBUG] 2022-10-03 09:33:16,000 __bs: MainProcess Chudnovsky ... 48,000,000 iterations and 495.46 seconds.
[DEBUG] 2022-10-03 09:33:27,921 __bs: MainProcess Chudnovsky ... 49,000,000 iterations and 507.38 seconds.
[DEBUG] 2022-10-03 09:33:34,778 __bs: MainProcess Chudnovsky ... 50,000,000 iterations and 514.24 seconds.
[DEBUG] 2022-10-03 09:33:43,298 __bs: MainProcess Chudnovsky ... 51,000,000 iterations and 522.76 seconds.
[DEBUG] 2022-10-03 09:33:49,959 __bs: MainProcess Chudnovsky ... 52,000,000 iterations and 529.42 seconds.
[DEBUG] 2022-10-03 09:34:29,294 __bs: MainProcess Chudnovsky ... 53,000,000 iterations and 568.75 seconds.
[DEBUG] 2022-10-03 09:34:36,176 __bs: MainProcess Chudnovsky ... 54,000,000 iterations and 575.63 seconds.
[DEBUG] 2022-10-03 09:34:41,576 __bs: MainProcess Chudnovsky ... 55,000,000 iterations and 581.03 seconds.
[DEBUG] 2022-10-03 09:34:50,161 __bs: MainProcess Chudnovsky ... 56,000,000 iterations and 589.62 seconds.
[DEBUG] 2022-10-03 09:34:56,811 __bs: MainProcess Chudnovsky ... 57,000,000 iterations and 596.27 seconds.
[DEBUG] 2022-10-03 09:35:09,382 __bs: MainProcess Chudnovsky ... 58,000,000 iterations and 608.84 seconds.
[DEBUG] 2022-10-03 09:35:16,206 __bs: MainProcess Chudnovsky ... 59,000,000 iterations and 615.66 seconds.
[DEBUG] 2022-10-03 09:35:24,295 __bs: MainProcess Chudnovsky ... 60,000,000 iterations and 623.75 seconds.
[DEBUG] 2022-10-03 09:35:31,095 __bs: MainProcess Chudnovsky ... 61,000,000 iterations and 630.55 seconds.
[DEBUG] 2022-10-03 09:35:52,139 __bs: MainProcess Chudnovsky ... 62,000,000 iterations and 651.60 seconds.
[DEBUG] 2022-10-03 09:35:58,781 __bs: MainProcess Chudnovsky ... 63,000,000 iterations and 658.24 seconds.
[DEBUG] 2022-10-03 09:36:07,399 __bs: MainProcess Chudnovsky ... 64,000,000 iterations and 666.86 seconds.
[DEBUG] 2022-10-03 09:36:12,847 __bs: MainProcess Chudnovsky ... 65,000,000 iterations and 672.30 seconds.
[DEBUG] 2022-10-03 09:36:19,763 __bs: MainProcess Chudnovsky ... 66,000,000 iterations and 679.22 seconds.
[DEBUG] 2022-10-03 09:36:32,351 __bs: MainProcess Chudnovsky ... 67,000,000 iterations and 691.81 seconds.
[DEBUG] 2022-10-03 09:36:39,078 __bs: MainProcess Chudnovsky ... 68,000,000 iterations and 698.53 seconds.
[DEBUG] 2022-10-03 09:36:47,830 __bs: MainProcess Chudnovsky ... 69,000,000 iterations and 707.29 seconds.
[DEBUG] 2022-10-03 09:36:54,701 __bs: MainProcess Chudnovsky ... 70,000,000 iterations and 714.16 seconds.
[DEBUG] 2022-10-03 09:39:39,357 __bs: MainProcess Chudnovsky ... 71,000,000 iterations and 878.81 seconds.
[DEBUG] 2022-10-03 09:39:46,199 __bs: MainProcess Chudnovsky ... 72,000,000 iterations and 885.66 seconds.
[DEBUG] 2022-10-03 09:39:54,956 __bs: MainProcess Chudnovsky ... 73,000,000 iterations and 894.41 seconds.
[DEBUG] 2022-10-03 09:40:01,639 __bs: MainProcess Chudnovsky ... 74,000,000 iterations and 901.10 seconds.
[DEBUG] 2022-10-03 09:40:14,219 __bs: MainProcess Chudnovsky ... 75,000,000 iterations and 913.68 seconds.
[DEBUG] 2022-10-03 09:40:19,680 __bs: MainProcess Chudnovsky ... 76,000,000 iterations and 919.14 seconds.
[DEBUG] 2022-10-03 09:40:26,625 __bs: MainProcess Chudnovsky ... 77,000,000 iterations and 926.08 seconds.
[DEBUG] 2022-10-03 09:40:35,212 __bs: MainProcess Chudnovsky ... 78,000,000 iterations and 934.67 seconds.
[DEBUG] 2022-10-03 09:40:41,914 __bs: MainProcess Chudnovsky ... 79,000,000 iterations and 941.37 seconds.
[DEBUG] 2022-10-03 09:41:03,218 __bs: MainProcess Chudnovsky ... 80,000,000 iterations and 962.68 seconds.
[DEBUG] 2022-10-03 09:41:10,213 __bs: MainProcess Chudnovsky ... 81,000,000 iterations and 969.67 seconds.
[DEBUG] 2022-10-03 09:41:18,344 __bs: MainProcess Chudnovsky ... 82,000,000 iterations and 977.80 seconds.
[DEBUG] 2022-10-03 09:41:25,261 __bs: MainProcess Chudnovsky ... 83,000,000 iterations and 984.72 seconds.
[DEBUG] 2022-10-03 09:41:37,663 __bs: MainProcess Chudnovsky ... 84,000,000 iterations and 997.12 seconds.
[DEBUG] 2022-10-03 09:41:44,680 __bs: MainProcess Chudnovsky ... 85,000,000 iterations and 1004.14 seconds.
[DEBUG] 2022-10-03 09:41:53,411 __bs: MainProcess Chudnovsky ... 86,000,000 iterations and 1012.87 seconds.
[DEBUG] 2022-10-03 09:41:58,926 __bs: MainProcess Chudnovsky ... 87,000,000 iterations and 1018.38 seconds.
[DEBUG] 2022-10-03 09:42:05,858 __bs: MainProcess Chudnovsky ... 88,000,000 iterations and 1025.32 seconds.
[DEBUG] 2022-10-03 09:42:46,163 __bs: MainProcess Chudnovsky ... 89,000,000 iterations and 1065.62 seconds.
[DEBUG] 2022-10-03 09:42:53,054 __bs: MainProcess Chudnovsky ... 90,000,000 iterations and 1072.51 seconds.
[DEBUG] 2022-10-03 09:43:02,030 __bs: MainProcess Chudnovsky ... 91,000,000 iterations and 1081.49 seconds.
[DEBUG] 2022-10-03 09:43:09,192 __bs: MainProcess Chudnovsky ... 92,000,000 iterations and 1088.65 seconds.
[DEBUG] 2022-10-03 09:43:21,533 __bs: MainProcess Chudnovsky ... 93,000,000 iterations and 1100.99 seconds.
[DEBUG] 2022-10-03 09:43:28,643 __bs: MainProcess Chudnovsky ... 94,000,000 iterations and 1108.10 seconds.
[DEBUG] 2022-10-03 09:43:37,372 __bs: MainProcess Chudnovsky ... 95,000,000 iterations and 1116.83 seconds.
[DEBUG] 2022-10-03 09:43:44,558 __bs: MainProcess Chudnovsky ... 96,000,000 iterations and 1124.02 seconds.
[DEBUG] 2022-10-03 09:44:06,555 __bs: MainProcess Chudnovsky ... 97,000,000 iterations and 1146.01 seconds.
[DEBUG] 2022-10-03 09:44:12,220 __bs: MainProcess Chudnovsky ... 98,000,000 iterations and 1151.68 seconds.
[DEBUG] 2022-10-03 09:44:19,278 __bs: MainProcess Chudnovsky ... 99,000,000 iterations and 1158.74 seconds.
[DEBUG] 2022-10-03 09:44:28,323 __bs: MainProcess Chudnovsky ... 100,000,000 iterations and 1167.78 seconds.
[DEBUG] 2022-10-03 09:44:35,211 __bs: MainProcess Chudnovsky ... 101,000,000 iterations and 1174.67 seconds.
[DEBUG] 2022-10-03 09:44:48,331 __bs: MainProcess Chudnovsky ... 102,000,000 iterations and 1187.79 seconds.
[DEBUG] 2022-10-03 09:44:54,835 __bs: MainProcess Chudnovsky ... 103,000,000 iterations and 1194.29 seconds.
[DEBUG] 2022-10-03 09:45:03,869 __bs: MainProcess Chudnovsky ... 104,000,000 iterations and 1203.33 seconds.
[DEBUG] 2022-10-03 09:45:10,967 __bs: MainProcess Chudnovsky ... 105,000,000 iterations and 1210.42 seconds.
[DEBUG] 2022-10-03 09:46:32,760 __bs: MainProcess Chudnovsky ... 106,000,000 iterations and 1292.22 seconds.
[DEBUG] 2022-10-03 09:46:39,872 __bs: MainProcess Chudnovsky ... 107,000,000 iterations and 1299.33 seconds.
[DEBUG] 2022-10-03 09:46:48,948 __bs: MainProcess Chudnovsky ... 108,000,000 iterations and 1308.41 seconds.
[DEBUG] 2022-10-03 09:46:54,611 __bs: MainProcess Chudnovsky ... 109,000,000 iterations and 1314.07 seconds.
[DEBUG] 2022-10-03 09:47:01,727 __bs: MainProcess Chudnovsky ... 110,000,000 iterations and 1321.18 seconds.
[DEBUG] 2022-10-03 09:47:14,525 __bs: MainProcess Chudnovsky ... 111,000,000 iterations and 1333.98 seconds.
[DEBUG] 2022-10-03 09:47:21,682 __bs: MainProcess Chudnovsky ... 112,000,000 iterations and 1341.14 seconds.
[DEBUG] 2022-10-03 09:47:30,610 __bs: MainProcess Chudnovsky ... 113,000,000 iterations and 1350.07 seconds.
[DEBUG] 2022-10-03 09:47:37,176 __bs: MainProcess Chudnovsky ... 114,000,000 iterations and 1356.63 seconds.
[DEBUG] 2022-10-03 09:47:59,642 __bs: MainProcess Chudnovsky ... 115,000,000 iterations and 1379.10 seconds.
[DEBUG] 2022-10-03 09:48:06,702 __bs: MainProcess Chudnovsky ... 116,000,000 iterations and 1386.16 seconds.
[DEBUG] 2022-10-03 09:48:15,483 __bs: MainProcess Chudnovsky ... 117,000,000 iterations and 1394.94 seconds.
[DEBUG] 2022-10-03 09:48:22,537 __bs: MainProcess Chudnovsky ... 118,000,000 iterations and 1401.99 seconds.
[DEBUG] 2022-10-03 09:48:35,714 __bs: MainProcess Chudnovsky ... 119,000,000 iterations and 1415.17 seconds.
[DEBUG] 2022-10-03 09:48:41,321 __bs: MainProcess Chudnovsky ... 120,000,000 iterations and 1420.78 seconds.
[DEBUG] 2022-10-03 09:48:48,408 __bs: MainProcess Chudnovsky ... 121,000,000 iterations and 1427.87 seconds.
[DEBUG] 2022-10-03 09:48:57,138 __bs: MainProcess Chudnovsky ... 122,000,000 iterations and 1436.60 seconds.
[DEBUG] 2022-10-03 09:49:04,328 __bs: MainProcess Chudnovsky ... 123,000,000 iterations and 1443.79 seconds.
[DEBUG] 2022-10-03 09:49:46,274 __bs: MainProcess Chudnovsky ... 124,000,000 iterations and 1485.73 seconds.
[DEBUG] 2022-10-03 09:49:52,833 __bs: MainProcess Chudnovsky ... 125,000,000 iterations and 1492.29 seconds.
[DEBUG] 2022-10-03 09:50:01,786 __bs: MainProcess Chudnovsky ... 126,000,000 iterations and 1501.24 seconds.
[DEBUG] 2022-10-03 09:50:08,975 __bs: MainProcess Chudnovsky ... 127,000,000 iterations and 1508.43 seconds.
[DEBUG] 2022-10-03 09:50:21,850 __bs: MainProcess Chudnovsky ... 128,000,000 iterations and 1521.31 seconds.
[DEBUG] 2022-10-03 09:50:28,962 __bs: MainProcess Chudnovsky ... 129,000,000 iterations and 1528.42 seconds.
[DEBUG] 2022-10-03 09:50:34,594 __bs: MainProcess Chudnovsky ... 130,000,000 iterations and 1534.05 seconds.
[DEBUG] 2022-10-03 09:50:43,647 __bs: MainProcess Chudnovsky ... 131,000,000 iterations and 1543.10 seconds.
[DEBUG] 2022-10-03 09:50:50,724 __bs: MainProcess Chudnovsky ... 132,000,000 iterations and 1550.18 seconds.
[DEBUG] 2022-10-03 09:51:12,742 __bs: MainProcess Chudnovsky ... 133,000,000 iterations and 1572.20 seconds.
[DEBUG] 2022-10-03 09:51:19,799 __bs: MainProcess Chudnovsky ... 134,000,000 iterations and 1579.26 seconds.
[DEBUG] 2022-10-03 09:51:28,824 __bs: MainProcess Chudnovsky ... 135,000,000 iterations and 1588.28 seconds.
[DEBUG] 2022-10-03 09:51:35,324 __bs: MainProcess Chudnovsky ... 136,000,000 iterations and 1594.78 seconds.
[DEBUG] 2022-10-03 09:51:48,419 __bs: MainProcess Chudnovsky ... 137,000,000 iterations and 1607.88 seconds.
[DEBUG] 2022-10-03 09:51:55,634 __bs: MainProcess Chudnovsky ... 138,000,000 iterations and 1615.09 seconds.
[DEBUG] 2022-10-03 09:52:04,435 __bs: MainProcess Chudnovsky ... 139,000,000 iterations and 1623.89 seconds.
[DEBUG] 2022-10-03 09:52:11,583 __bs: MainProcess Chudnovsky ... 140,000,000 iterations and 1631.04 seconds.
[DEBUG] 2022-10-03 09:52:17,222 __bs: MainProcess Chudnovsky ... 141,000,000 iterations and 1636.68 seconds.
[DEBUG] 2022-10-03 10:02:43,939 compute: MainProcess    Chudnovsky brothers  1988 
    π = (Q(0, N) / 12T(0, N) + 12AQ(0, N))**(C**(3/2))
 calulation Done! 141,027,339 iterations and 2263.39 seconds.
[INFO] 2022-10-03 10:09:07,119 <module>: MainProcess Last 5 digits of π were 45519 as expected at offset 999,999,995
[INFO] 2022-10-03 10:09:07,119 <module>: MainProcess Calculated π to 1,000,000,000 digits using a formula of:
 10     Chudnovsky brothers  1988 
    π = (Q(0, N) / 12T(0, N) + 12AQ(0, N))**(C**(3/2))
 
[INFO] 2022-10-03 10:09:07,120 <module>: MainProcess Calculation took 141,027,339 iterations and 0:44:06.398345.

math_pi.pi(b = 1000000) is faster to a million. About 40 times faster. But it cannot go to a Billion. 1 Million is the most digits.

The GMPY Builtin looks like:

python pi-pourri.py -v -d 1,000,000,000 -a 11
[INFO] 2022-10-03 14:33:34,729 <module>: MainProcess Computing π to 1,000,000,000 digits.
[DEBUG] 2022-10-03 14:33:34,729 compute: MainProcess Starting   const_pi() function from the gmpy2 library formula to 1,000,000,000 decimal places
[DEBUG] 2022-10-03 15:46:46,575 compute: MainProcess    const_pi() function from the gmpy2 library calulation Done! 1 iterations and 4391.85 seconds.
[INFO] 2022-10-03 15:46:46,575 <module>: MainProcess Last 5 digits of π were 45519 as expected at offset 999,999,995
[INFO] 2022-10-03 15:46:46,575 <module>: MainProcess Calculated π to 1,000,000,000 digits using a formula of:
 11     const_pi() function from the gmpy2 library 
[INFO] 2022-10-03 15:46:46,575 <module>: MainProcess Calculation took 1 iterations and 1:13:11.845652.

The MPmath Builtin is almost as fast. About 12% (6 minutes) slower:

python pi-pourri.py -v -a 12 -d 1,000,000,000  
[INFO] 2022-10-04 09:10:37,085 <module>: MainProcess Computing π to 1,000,000,000 digits.
[DEBUG] 2022-10-04 09:10:37,085 compute: MainProcess Starting   mp.pi() function from the mpmath library formula to 1,000,000,000 decimal places
[DEBUG] 2022-10-04 10:01:25,321 compute: MainProcess    mp.pi() function from the mpmath library calulation Done! 1 iterations and 3048.22 seconds.
[INFO] 2022-10-04 10:01:25,338 <module>: MainProcess Last 5 digits of π were 45519 as expected at offset 999,999,995
[INFO] 2022-10-04 10:01:25,340 <module>: MainProcess Calculated π to 1,000,000,000 digits using a formula of:
 12     mp.pi() function from the mpmath library 
[INFO] 2022-10-04 10:01:25,343 <module>: MainProcess Calculation took 1 iterations and 0:50:48.250337.
Cymric answered 19/9, 2022 at 18:6 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.