Created
October 21, 2021 17:49
-
-
Save marcusmueller/37a6247319b8313b9504019f6a7b8992 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// Copyright Marcus Müller | |
// SPDX-License-Identifier: LGPL-3.0-or-later | |
#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); } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment