More fun with constexpr

pi.cpp

// Copyright Dave Abrahams 2011. Distributed under the Boost
// Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <iostream>
#include <iomanip>
#include <boost/math/constants/constants.hpp>

#ifndef TRACE
# define CONSTEXPR constexpr
#else
# define CONSTEXPR
#endif

//
// Calculate pi via the Van Wijngaarden transformation without relying
// on memoization.
//
// TODO: finish this comment
// π₀,₁ = 4;  π₀,ⱼ = π₀,ⱼ-1 ...

//
// node: a simple list link
//
template <class T>
struct node
{
    constexpr node(T value, node const* next = nullptr)
        : value(value), next(next)
    {}

    constexpr node const&
    nth(unsigned n)
    {
        return n == 0 ? *this : next->nth(n-1);
    }

    constexpr node push(T x)
    {
        return node(x, this);
    }
    
    T value;
    node const* next;
};

// position - describes the values that (i,j) pairs as we fill in the
// chart of values:
//
//  diagonal # d == i+j
//
//  length of diagonal == d/2+1
//
//  visitation order
//
//     j
//     0  1  2  3  4  5  6
//    +-------------------
// i 0|0  1  2  4  6  9  12
//   1|-  3  5  7  10 13
//   2|-  -  8  11 14
//   3|-  -  -  15
//
// index in history of element above
//
//     j
//     0  1  2  3  4  5  6
//    +-------------------
// i 0|-  -  -  -  -  -  -
//   1|-  1  2  2  3  3
//   2|-  -  2  3  3
//   3|-  -  -  3
//
// index in history of element left
//
//     j
//     0  1  2  3  4  5  6
//    +-------------------
// i 0|-  0  0  1  1  2  2
//
struct position
{
    constexpr position(unsigned i, unsigned j)
        : i(i), j(j) {}
    
    constexpr position next()
    {
        return i+1 <= j-1 ? position{i+1,j-1} : position{0,i+j+1};
    }

    unsigned i, j;
};

inline std::ostream& operator<<(std::ostream& s, position p)
{
    s << "( " << p.i << ", " << p.j << ")";
}

template <class T>
CONSTEXPR
T calculate_pi(position pos, node<T> const& history);

template <class T>
CONSTEXPR
T calculate_pi_(position pos, node<T> const& history, T next_value)
{
#ifdef TRACE
    std::cout << "π" << pos << " = " << std::setprecision(100) << next_value << std::endl;
#endif
    // the problem with this method is that it's hard to know where to
    // stop because it doesn't converge monotonically.  (20,22) seems
    // to give maximally precise long doubles on my machine
    return
        pos.i == 20 && pos.j == 22
        ? next_value
        : calculate_pi(pos.next(), history.push(next_value));
}

template <class T>
constexpr T pi0_next(unsigned j, T pi0_prev)
{
    return pi0_prev + (j%2?T(4.0):T(-4.0))/(2*j-1);
}

template <class T>
CONSTEXPR
T calculate_pi(position pos, node<T> const& history)
{
    
    return calculate_pi_(
        pos, history,
        pos.i == 0 ?
        pi0_next(pos.j, history.nth((pos.j-1)/2).value) 
        : (history.value + history.nth((pos.i+pos.j+1)/2).value)/T(2.0)
        );
}

template <class T>
CONSTEXPR
T pi_()
{
    return calculate_pi(position(0,1), node<T>(0.0));
};

#ifndef TRACE
constexpr long double pi = pi_<long double>();
#endif

int main()
{
    using number = std::complex<long double>;
    
    std::cout << std::setprecision(100);

#ifdef TRACE
    std::cout << pi_<long double>() << std::endl;
#else
    std::cout << pi << std::endl;
#endif

    std::cout << number(4.0) * (number(4.0)*atan(number(1.0)/number(5))        -   atan( number(1.0)/number(239) ) ) << std::endl;
    std::cout << boost::math::constants::pi<long double>() << std::endl;
}