boost:math:factorial2 throws an error while computing double factorial of -1?

Question

The official documentation of boost library in C++ confirms that -1!! is defined. However, when I try to compute the double factorial of -1, it throws the following error "Error in function boost::math::tgamma result of gamma is too large to represent". I can implement a code based on iteration to compute the same (if it comes to that), but would like to use optimized libraries wherever possible. Any suggestions on rectifying this issue?

#include <iostream>
#include  <boost/math/special_functions/factorials.hpp>

int main(int argc, char const *argv[]) {
  double t  = boost::math::double_factorial<double>(-1);
  std::cout << t <<"\n";
  return 0;
}

This is the minimal working example.

Answer 1

This is how boost::math::double_factorial is declared:

namespace boost{ namespace math{

template <class T>
T double_factorial(unsigned i);

template <class T, class Policy>
T double_factorial(unsigned i, const Policy&);

}} // namespaces

According to the documentation for boost::math::double_factorial

The argument to double_factorial is type unsigned even though technically -1!! is defined.

This means that, although mathematically -1!! is well defined, the function does not support negative inputs. In fact, -1 will be silently converted to an unsigned, which will result in a large number-- hence the overflow.

As a workaround, you could declare an additional double_factorial function inside the namespace boost::math

Example, where for simplicity I have allowed for just -1!!.

#include <iostream>
#include <exception>
#include  <boost/math/special_functions/factorials.hpp>

namespace boost { namespace math {
    template <class T>
    T double_factorial(int i) { if (i == -1) return T{1}; else throw(std::runtime_error("Unknown double factorial")); }
} }

int main(int argc, char const *argv[]) {
  double t  = boost::math::double_factorial<double>(-1);
  std::cout << t <<"\n";
  return 0;
}

See it live on Coliru

Warning

This solution will work, because the compiler will choose your new overload, inside the same namespace, when you call double_factorial with a negative argument. It is however potentially dangerous: if a future version of boost allows for negative values, the code may not compile anymore, or ignore the new boost version.

A better solution is to wrap around double_gamma to a function you define, like as follows

template <class T>
T my_double_factorial(int i) {
    if (i > 0)
       return boost::math::double_factorial<T>(i);
    else if (i == -1)
       return T{1};
    else
       throw(std::runtime_error("Unknown double factorial"));
}