Mathematical Constants in C++

None of C++ standards has defined π number. So every programmer who needs to do some math including, for example, trigonomerty or statistics, has to define π herself. Luckily header files have defined π (See M_PI in math.h). Though it is not ultimately safe and portable to use these constants. It might casue numeric errors in critical applications. In this blog post, I am going to show a better way to reach a safe π.

Fortunately it is possible to calculate π number at compile time, having zero overhead in runtime, while complying to IEEE 754 completely. Well, the latter part requires your compile flags to be set properly, but that’s not the point. Using this method your π will always be the same on all architectures that implement IEEE 754.

This implementation, uses C++’s meta-programming features (that’s a fancy word for templates. Not really, but almost). There are many ways to converge a value towards π. One of them is this nice serie called Bailey–Borwein–Plouffe formula, discovered in 1995 by, David H. Bailey, Peter Borwein, and Simon Plouffe. It has been proven that:

$$\pi=\sum_{n=0}^{\infty}{\left[\frac{1}{16^{n}}\left(\frac{4}{8n+1}-\frac{2}{8n +4}-\frac{1}{8n+5}-\frac{1}{8n+6}\right)\right]}$$

Thus one can calculate π number, to whatever precision she wants.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

#include <cstdint>
#include <limits>
#include <iostream>

template <uint64_t b, uint64_t e>
struct pow {
    static const uint64_t result = b * pow < b, e - 1 >::result;
};

template <uint64_t b>
struct pow<b, 0> {
    static const uint64_t result = 1;
};

template <uint64_t n>
struct bbp {
    constexpr static double pi
            = (1.0 / pow<16, n>::result)
            * (   4.0 / (8 * n + 1.0) - 2.0 / (8 * n + 4.0)
                - 1.0 / (8 * n + 5.0) - 1.0 / (8 * n + 6.0))
            + bbp < n - 1 >::pi;
};

template <>
struct bbp < -1 > {
    constexpr static double pi = 0;
};

int main(int, char* []) {
    std::cout.precision(std::numeric_limits<double>::digits10);
    std::cout << "pi = " << bbp<1>::pi << std::endl;
    std::cout << "pi = " << bbp<5>::pi << std::endl;
    std::cout << "pi = " << bbp<10>::pi << std::endl;
    return 0;
}

Results:

pi: 3.14142246642247
pi: 3.14159265322809
pi: 3.14159265358979
pi: 3.14159265358979

A Note on Other Series

There are many series that converge towards π. For example Newton’s method, approximates π using it’s geometrical properties with this formula:

$$\pi=\frac{3}{4}\sqrt{3}+24\int_{0}^{\frac{1}{4}}{\sqrt{x-x^{2}}},dx$$

This one is computationally intensive and it can’t be simply calculated at compile time.

There is also this one which moves towards π slower compared to BBP method:

$$\pi=8\sum _{n=0}^{\infty }{\frac {1}{(4n+1)(4n+3)}}$$

And there is also this formula which is both computationally intensive and slowly moving towards π:

$$\pi=\sum_{n=0}^{\infty}{\frac{\left({n!}^{2}2^{n-1}\right)}{\left(2n+1\right)! }}$$

Inspecting lots of formulas to reach π, I found BBP method the easiest to implement with less errors. Leave a comment if you know about a better method (:

Other constants

π is not the only irrational number used in math. There are others like e (Euler’s number) which is irrational and one can move towards it using this series:

$$e=\sum \limits _{n=0}^{\infty }{\frac {1}{n!}}={\frac {1}{1}}+{\frac {1}{1}}+{\frac {1}{1\cdot 2}}+{\frac {1}{1\cdot 2\cdot 3}}+\cdots$$

But there’s no need to implement this one. e is already defined in C and C++ standards, via exp function.