newton.cc

#include <fmt/format.h>

constexpr unsigned N_SUMMANDS = 9;
constexpr unsigned N_ITER = N_SUMMANDS - 2;

constexpr unsigned long long faculty(unsigned long long n)
{
    return n > 1 ? n * faculty(n - 1) : 1;
}

template <class base_t>
constexpr base_t power(base_t base, unsigned exponent)
{
    if (!exponent)
        return 1;
    if (exponent & 0x1)
        return base * power(base, exponent / 2) * power(base, exponent / 2);
    return power(base, exponent / 2) * power(base, exponent / 2);
}
static_assert(power<unsigned int>(2, 10) == 1024);

template <class fp_t>
constexpr fp_t sin_taylor(fp_t x, unsigned n)
{
    return static_cast<fp_t>(power(x, 2 * n + 1)) / faculty(2 * n + 1) *
           ((n & 0x1) ? (-1) : 1);
}

template <class fp_t>
constexpr fp_t cos_taylor(fp_t x, unsigned n)
{
    return static_cast<fp_t>(power(x, 2 * n)) / faculty(2 * n) * ((n & 0x1) ? (-1) : 1);
}

template <class fp_t = double>
constexpr fp_t sin(fp_t x, unsigned int elements = N_SUMMANDS)
{
    fp_t sum = 0;
    for (unsigned idx = 0; idx <= elements; ++idx) {
        sum += sin_taylor(x, idx);
    }
    return sum;
}
template <class fp_t = double>
constexpr fp_t cos(fp_t x, unsigned int elements = N_SUMMANDS)
{
    fp_t sum = 0;
    for (unsigned idx = 0; idx <= elements; ++idx) {
        sum += cos_taylor(x, idx);
    }
    return sum;
}

static_assert(sin(0, 1) == 0);

template <class fp_t = double>
constexpr fp_t f(fp_t x)
{
    return sin(x);
}
template <class fp_t = double>
constexpr fp_t f_prime(fp_t x)
{
    return cos(x);
}
template <class fp_t = double>
constexpr fp_t newtonrecursion(const unsigned int n, const fp_t x)
{
    return n ? newtonrecursion(n - 1, x - f(x) / (f_prime(x))) : x;
}

constexpr double pi = newtonrecursion(N_ITER, 3.0f);

int main() { fmt::print("{:.20f}", pi); }