Estimate 𝜋 using Machin's formula

Question

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

Answer 1

Answers to your questions

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.

Answer 2

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;
}