# Overview

Below is a simple program to estimate pi to the full precision of the provided type. It uses Machin's formula, but this question is about coding style / choices. It compares the result to M_PI from <cmath> - best we have in C++17.

# Questions

1. Is my use of enable_if_t appropriate and idiomatic?
2. Is the use of static_cast<T> for one_over_5 and one_over_239 the "best" solution here? Neccesary, and appropriate?
3. How about the use of decltype, ::epsilon and ::digits10?
4. Accuracy: How about sign =* -1 and den += 2.0 and ... = sign * ...? These are integer operations but use the floating point type.
5. My choice to force the user of pi<T>() to specify a T?

Note: I know the templates should probably be in header file, but all in one here for simplicity on CR.

#include <cmath>
#include <iomanip>
#include <iostream>
#include <type_traits>

template <typename T, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
T arctan(T tan) {
// fabs(tan) > 1 will not converge
T num  = tan;
T den  = 1.0;
T sign = 1.0;
T term = num;
T res  = num;
while (std::fabs(term) > std::numeric_limits<T>::epsilon()) {
num *= tan * tan;
den += 2.0;
sign *= -1.0;
term = sign * num / den;
res += term;
}
return res;
}

template <typename T, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
T pi() {
T one_over_5   = static_cast<T>(1.0) / static_cast<T>(5.0);
T one_over_239 = static_cast<T>(1.0) / static_cast<T>(239.0);
return (static_cast<T>(4.0) * arctan<T>(one_over_5) - arctan<T>(one_over_239)) *
static_cast<T>(4.0);
}

int main() {
auto pi_float = pi<float>();
auto pi_double = pi<double>();
auto pi_long_double = pi<long double>();
std::cout << std::setprecision(std::numeric_limits<decltype(pi_float)>::digits10)
<< pi_float << '\n'
<< std::setprecision(std::numeric_limits<decltype(pi_double)>::digits10)
<< pi_double << '\n'
<< std::setprecision(std::numeric_limits<decltype(pi_long_double)>::digits10)
<< pi_long_double << '\n'
<< M_PI << '\n';
return EXIT_SUCCESS;
}


Output - demonstrating that M_PI from <cmath> is actually incorrect at long double precision in the last 2 digits (the correct value is 3.1415926535897932384626433...)

Note (unrelated to question, but for future reference): It turns out that the reason M_PI is "wrong" is because the M_PI macro does not specify a literal L for long double, and hence the value is rounded to double precision before being displayed. If you put a literal 3.14159265358979323846L then it is accurate. Usually that is provided as the macro M_PIl.

3.14159
3.14159265358979
3.14159265358979324
3.14159265358979312


Is my use of enable_if_t appropriate and idiomatic?

Yes, it is appropriate, and I believe idiomatic for C++17 as well. In C++20 the idiomatic way would be to use concepts of course.

Is the use of static_cast for one_over_5 and one_over_239 the "best" solution here? Neccesary, and appropriate?

It's necessary to cast at least a few constants to ensure the expression has the right type. However, in this case I would use the constructor syntax instead, and only do it for the minimum number of constants necessary, to keep the expression concise and legible:

return (4 * arctan(T{1} / 5) - arctan(T{1} / 239)) * 4;


Also note that you don't need to specify <T> here, it can be automatically deduced, and since this is so short there is no need to declare temporary variables one_over_5 and one_over_239.

How about the use of decltype, ::epsilon and ::digits10?

The use of T::digits10 is fine here. However, using T::epsilon() is always very dangerous. It works here only because you are trying to calculate something close to 1, or if you wanted to calculate the arctangent from a very small number you have the advantage that $$\\lim_{x\to 0} \arctan x = x\$$. In general, if you want to calculate a value using some converging expansion and want to stop if further steps are no longer giving you any increase in precision, you can do that by checking if there is any change in the floating point value:

while (res + std::fabs(term) != res) {...}


Or perhaps using std::nextafter(), to keep in the spirit of T::epsilon():

while (res + std::fabs(term) > std::nextafter(res, res + 1)) {...}


Although for std::nextafter(), it might even be hard to come up with the correct second argument.

Again, your method works fine in this case, it produces the same output as both methods above.

Accuracy: How about sign =* -1 and den += 2.0 and ... = sign * ...? These are integer operations but use the floating point type.

Currently those are floating point operations because the left hand side is a floating point type, and if the other side is an integer it would be promoted to the floating point type before the operation is performed. Note that you can replace all the constants in arctan() with integers and the result will still be correct.

You could change the type of sign and den to be an integer type, but it is not necessary. Multiplying floating point sign by -1 will not make it drift away from 1/-1 over time, and since integers up to $$\2^{24}\$$ can be represented exactly by a float, you will not have any problems with den either, unless the loop would run for more than 8388607 iterations, which would mean you have bigger problems.

My choice to force the user of pi() to specify a T?

It forces the user to give the correct type, which might result in less errors and/or more optimal code. You could consider making it double by default:

template <typename T = double, ...>
T pi() {...}


Note that C++20 introduces std::pi_v<> which also doesn't not have a default type.

# Make the functions constexpr

You can and should make the functions constexpr. Interestingly, the compiler would be free to calculate the results of pi_float, pi_double and pi_long_double at compile time regardless, however GCC and Clang don't without constexpr. With constexpr, GCC does calculate the result at compile time, but Clang still emits a lot of code.

• Excellent as always, thank you. What I meant by "These are integer operations" was,:I perceived there is a tradeoff here: Keep them int and the *-1 and +=2 operations are then integer ops, which are both more accurate and possibly faster, but at the cost of a integer => FP conversion on each iteration? In hindsight, integer for accuracy over speed is probably better? Dec 11, 2020 at 21:53
• Ah, instead of T den, sign have int den, sign? It shouldn't matter; sign will be always be precise, and den will only become imprecise after 4194304 iterations if T = float. Dec 11, 2020 at 21:59
• Yes, sign (ie both 1, -1) are non-approximate in Base2 right? I was also more concerned about den How did you figure this? den will only become imprecise after 4194304 iterations if T = float ?? Dec 11, 2020 at 22:05
• I updated the answer. It's after 8388607 iterations actually, integer values up to 16777216 can be represented exactlty by float, see this question Dec 11, 2020 at 22:08
• for a plain is_floating_point condition, clang's error message is actually quite good on an enable_if: candidate template ignored: requirement 'std::is_floating_point_v<int>' was not satisfied [with T = int] constexpr T arctan(T tan) { Dec 12, 2020 at 6:16

Updated code, incorporating @G. Sliepen's feedback.

Also adding a second arctan() which uses Euler's series - old one renamed to arctan_taylor()

The new series converges faster and does not require special casing for abs(tan) > 1

#include <cmath>
#include <iomanip>
#include <iostream>
#include <type_traits>

template <typename T, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
constexpr T arctan_taylor(T tan) {
if (std::fabs(tan) > 1) return (tan < 0.0 ? -M_PI_2f64 : M_PI_2f64) - arctan_taylor(1 / tan);
T num  = tan;
T den  = 1;
T sign = 1;
T term = num;
T res  = num;
while (res + std::fabs(term) != res) {
num *= tan * tan;
den += 2; // will remain precise for T = float up to 2^(mantissa bits + 1) + 1 = 16,777,217
sign = -sign;
term = sign * num / den;
res += term;
}
return res;
}

template <typename T, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
constexpr T arctan(T tan) {
T sum  = 0;
T prod = 1;
for (int n = 0; sum + std::fabs(prod) != sum; ++n) {
prod = 1;
for (int k = 1; k <= n; ++k)
prod *= 2 * k * tan * tan / ((2 * k + 1) * (1 + tan * tan)); //  Eulers Series
sum += prod;
}
return tan / (1 + tan * tan) * sum;
}

template <typename T, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
constexpr T pi() {
return (4 * arctan(T(1) / 5) - arctan(T(1) / 239)) * 4;
}

int main() {
auto pi_float       = pi<float>();
auto pi_double      = pi<double>();
auto pi_long_double = pi<long double>();
std::cout << std::setprecision(std::numeric_limits<decltype(pi_float)>::digits10) << pi_float
<< '\n'
<< std::setprecision(std::numeric_limits<decltype(pi_double)>::digits10) << pi_double
<< '\n'
<< std::setprecision(std::numeric_limits<decltype(pi_long_double)>::digits10)
<< pi_long_double << '\n'
<< M_PI << '\n';

return EXIT_SUCCESS;
}

• Ah, a rare case where std::copysign() can be used: if (std::fabs(tan) > 1) return std::copysign(M_PI_2f64, tan) - arctan_taylor(1 / tan); Dec 14, 2020 at 16:23
• Looks like the inner for-loop in the version using the Euler series is not necessary? Just remove prod = 1 and replace the loop with prod *= 2 * (n + 1) * tan * tan / ((2 * n + 3) * (1 + tan * tan)); Dec 14, 2020 at 16:27
• I was directly implementing the formula from linked WP page. Not sure what you mean. The inner loop may run for more than one iteration for each n? The prod is reset to 1 for each inner loop (ie the "pi in the formula)? Maybe I missed something, I haven't looked at it very hard. But I think it's meant to be a "sum of multi-term-products". The inner loop doesn't run at all for n = 0 . once for n=1, twice for n=2 etc It gives the right answer. I haven't debugged, traced or optimised it other than that. It's already very fast. Dec 14, 2020 at 21:43
• Love the std::copysign() ! That's a new one to me! The new Euler algo actually also benefits from the same if fabs(tan) > 1. It will converge for bigger values, but it's much faster to use the 1/tan property. Dec 14, 2020 at 21:50