The implementation of std::tr1::tgamma looks wrong to me.
Keep in mind that the gamma function is equal to the factorial function for
integer values.
The code blow prints 4 columns:
1) a number, 2) its factorial, computed by hand 3) output using a corrected
version of "Numerical Recipes in C"'s log-gamma function 4) output from
std::tr1::tgamma() function from tr1/cmath
The corrected implementation of NR's gammln() computes the log-gamma of x+1
rather than that of x.
The program's output is then:
n n! exp(NR::gammln(n)) std::tr1::tgamma(n)
1 1 1 1
2 2 2 1
3 6 6 2
4 24 24 6
5 120 120 24
6 720 720 120
7 5040 5040 720
8 40320 40320 5040
9 362880 362880 40320
>From the table, one can see that std::tr1::tgamma(n) lacks the "x+1"
correction.
----------- cut here -------------
#include <cassert>
#include <iostream>
#include <tr1/cmath>
template <typename T>
T gammln(const T &x)
{
static const T cof [] = {76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.1208650973866179e-2, -0.5395239384953e-5};
const T y = x + 1.0; // Correction to NR
const T y55 = y+5.5;
const T tmp = y55 - (y+0.5) * std::log( y55 );
T ser = 1.90015e-10;
for (unsigned i = 0; i < 6; ++i)
ser += cof[i] / (y + 1.0 + static_cast<T>(i));
return -tmp + std::log(2.506628274631005 * (1.0 + ser) / y);
}
template <typename T>
T factorial(const T &x)
{
T r = 1;
for (unsigned i = 2; i <= x; ++i)
r *= i;
return r;
}
int main()
{
std::cout << " n\tn!\texp(NR::gammln(n))\tstd::tr1::tgamma(n)"
<< std::endl;
for (unsigned i = 1u; i < 10u; ++i)
{
std::cout << i << '\t' << factorial(double(i)) << "\t\t" <<
std::exp(gammln(double(i))) << "\t\t" <<
std::tr1::tgamma(double(i)) << std::endl;
}
return 0;
}
----------- cut here -------------
My gcc version:
Using built-in specs.
Target: x86_64-linux-gnu
Configured with: ../src/configure -v
--enable-languages=c,c++,fortran,objc,obj-c++,treelang --prefix=/usr
--enable-shared --with-system-zlib --libexecdir=/usr/lib
--without-included-gettext --enable-threads=posix --enable-nls
--with-gxx-include-dir=/usr/include/c++/4.1.3 --program-suffix=-4.1
--enable-__cxa_atexit --enable-clocale=gnu --enable-libstdcxx-debug
--enable-mpfr --enable-checking=release x86_64-linux-gnu
Thread model: posix
gcc version 4.1.3 20070812 (prerelease) (Ubuntu 4.1.2-15ubuntu2)
[EMAIL PROTECTED]:~/work/src/cvmlcpp/trunk/testing$ gcc-4.2 -v
Using built-in specs.
Target: x86_64-linux-gnu
Configured with: ../src/configure -v
--enable-languages=c,c++,fortran,objc,obj-c++,treelang --prefix=/usr
--enable-shared --with-system-zlib --libexecdir=/usr/lib
--without-included-gettext --enable-threads=posix --enable-nls
--with-gxx-include-dir=/usr/include/c++/4.2 --program-suffix=-4.2
--enable-clocale=gnu --enable-libstdcxx-debug --enable-mpfr --disable-werror
--enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu
--target=x86_64-linux-gnu
Thread model: posix
gcc version 4.2.1 (Ubuntu 4.2.1-2ubuntu1)
Best,
Fokko Beekhof
--
Summary: std::tr1::tgamma produces wrong results [for (x-1) in
stead of for x]
Product: gcc
Version: 4.2.1
Status: UNCONFIRMED
Severity: major
Priority: P3
Component: libstdc++
AssignedTo: unassigned at gcc dot gnu dot org
ReportedBy: fpbeekhof at gmail dot com
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33060
