calculating the Gamma function for negative real values (C++, Boost)

Question

The following infinite serious require the calculation of factorial for non-integer, negative, real numbers:

(it is a way to calculate the circumference of an ellipse, a and b are the semi-major and semi minor axis and h is defined as:
h = (a-b)^2/(a+b)^2)

The factorial function can be extended to negative values via the Gamma function which is defined for all the real numbers that are not negative integers.

While coding the serious I tried boost::math::factorial and boost::math::tgamma which gives results only down to -1 (not included) -1.5 for example give an error.

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

    int main()
    {
        double x;
        double f;
        double tg;

        x = -0.5;
        f = boost::math::factorial<double>(x);
        tg = boost::math::tgamma<double>(x);
        cout << "factorial of " << x << " = " << f << endl;
        cout << "tgamma of " << x << " = " << tg << endl << endl;

        x = -1.5;
        f = boost::math::factorial<double>(x);
        tg = boost::math::tgamma<double>(x);
        cout << "factorial of " << x << " = " << f << endl;
        cout << "tgamma of " << x << " = " << tg << endl << endl;
    
        return 0;
    }

the output:

factorial of -0.5 = 1
tgamma of -0.5 = -3.54491
terminate called after throwing an instance of 'boost::exception_detail::clone_implboost::exception_detail::error_info_injector<std::domain_error >' what(): Error in function boost::math::tgamma(long double): Evaluation of tgamma at a negative integer 0. Aborted (core dumped)

boost factorial: boost factorial
boost tgamma: boost tgamma

My questions:

Is there an alternative to boost that can calculate the gamma function for its negative domain?
I could not find in the boost documentation linked above mention of the implemented domain for factorials and tgamma functions. In fact I might simply be using them wrong. What is the way to ascertain what is indeed the domain/correct usage?

Thanks.

Just for curiosity, why did not you use `std::tgamma` ? https://en.cppreference.com/w/cpp/numeric/math/tgamma — Damien, Jul 23 '20 at 16:08
In the last boost tgamma documentation, a figure is provided, which shows that you should not have problem in calculating it for x = -1.5. https://www.boost.org/doc/libs/1_73_0/libs/math/doc/html/math_toolkit/sf_gamma/tgamma.html — Damien, Jul 23 '20 at 16:39
@Damien thank you, I did not know that std::tgamma exists. and yes, the image from the boost documentation apply that boost::tgamma function is indeed defined for the negative domain, but the domain is not defined explicitly. — sivan shani, Jul 24 '20 at 14:03

score 2 · Answer 1 · answered Jul 24 '20 at 14:07

Gamma functions are part of the standard library since C++11 [1].

The usage is as follows:

#include <cmath>

std::tgamma(-0.5) # -3.5449077
std::lgamma(-0.5) # 1.2655121

You might as well use tgammal and tgammaf for long double and float types correspondingly.

https://en.cppreference.com/w/cpp/numeric/math

score 2 · Accepted Answer · answered Jul 24 '20 at 14:12

I understand what's going wrong. The boost::math::factorial function takes an unsigned integer by definition:

template <class T>
inline T factorial(unsigned i)
{
   return factorial<T>(i, policies::policy<>());
}

This means that if you call it with a double, it will get implicitly converted to unsigned. That's not what you want. Also, factorial ends up using tgamma internally, so you get this:

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

void foo(long double x) {
    using namespace boost::math;
    try {
        auto f = factorial<long double>(x);
        std::cout << "factorial of " << static_cast<unsigned>(x) << " = " << f << "\n";
    } catch(std::exception const& e) {
        std::cout << "error at " << static_cast<unsigned>(x) << ": " << std::quoted(e.what()) << "\n";
    }
}

int main() {
    std::cout << std::unitbuf;
    foo(-2);
}

Will end up doing this:

#0  boost::math::tgamma<long double, boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy> > (a=4294967295, z=...)
    at /home/sehe/custom/boost_1_73_0/boost/math/special_functions/gamma.hpp:1994
No locals.
#1  0x0000555555558eb3 in boost::math::factorial<long double, boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy> > (i=4294967294, pol=...)
    at /home/sehe/custom/boost_1_73_0/boost/math/special_functions/factorials.hpp:44
        result = -0.667762310955655363645
#2  0x0000555555558674 in boost::math::factorial<long double> (i=4294967294)
    at /home/sehe/custom/boost_1_73_0/boost/math/special_functions/factorials.hpp:53
No locals.
#3  0x0000555555557792 in foo (x=-2) at /home/sehe/Projects/stackoverflow/test.cpp:7
        f = <invalid float value>
#4  0x000055555555791f in main () at /home/sehe/Projects/stackoverflow/test.cpp:16
No locals.

So it's attempting to give you boost::math::factorial<long double> (i=4294967294)

Fix

Don't use factorials for other than non-negative integers.

Live On Compiler Explorer

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

void foo(long double x) {
    using namespace boost::math;
    try {
        auto tg = tgamma<long double>(x);
        std::cout << "tgamma    of " << x << " = " << tg << "\n" << std::endl;
    } catch(std::exception const& e) {
        std::cout << "error at " << x << ": " << std::quoted(e.what()) << std::endl;
    }
}

int main() {
    for (auto x : { 1., 2., 3., 4., 5., -.2, -2., -.5, -1.5 })
        foo(x);
}

Prints:

tgamma    of 1 = 1

tgamma    of 2 = 1

tgamma    of 3 = 2

tgamma    of 4 = 6

tgamma    of 5 = 24

tgamma    of -0.2 = -5.82115

error at -2: "Error in function boost::math::tgamma<long double>(long double): Evaluation of tgamma at a negative integer -2."
tgamma    of -0.5 = -3.54491

tgamma    of -1.5 = 2.36327

It understandably overflows at -2, but that's correct.

thank you very much. I tried both the std::tgamma and also your code and both works. much appreciated. — sivan shani, Jul 24 '20 at 14:32

user14717 · Answer 3 · 2020-07-24T13:49:49.120

1

I just tried to reproduce your problem using g++-10 on Mac, using boost.math commit hash 87929356790ad0 (current develop), but I get:

factorial of -0.5 = 1
tgamma of -0.5 = -3.54491

factorial of -1.5 = 1
tgamma of -1.5 = 2.36327

So more information is necessary, Boost version, OS, compiler, so on.

edited Jul 24 '20 at 13:49

answered Jul 23 '20 at 19:46

user14717

4,757
2
44
68

Maybe you can supply _your_ platform details as well? – sehe Jul 24 '20 at 13:46
@user14717 thank you very much, at least now I know that in general it supposed to work even if not in the version I have, my platform details: Ubuntu 18.04 gcc version 9.3.0 (Ubuntu 9.3.0-11ubuntu0~18.04.1) libboost-dev version 1.65.1.0ubuntu1 – sivan shani Jul 24 '20 at 14:07

calculating the Gamma function for negative real values (C++, Boost)

3 Answers3

Fix