This boils down to eglibc working with 8-byte floats on ARM, even if it's declared long double.
To see this behaviour on x86, build and run the following:
#include <stdio.h>
#include <math.h>
int main ()
{
double x = 6.0;
printf("sizof(double)=%d\n",sizeof(double));
printf("lgamma (%.20f)=%.20f\n", x, lgamma(x));
printf("tgamma (%.20f)=%.20f\n", x, tgamma(x));
printf("exp(lgamma) (%.20f)=%.20f\n", x, exp(lgamma(x)));
return 0;
}
This boils down to eglibc working with 8-byte floats on ARM, even if it's declared long double.
To see this behaviour on x86, build and run the following:
#include <stdio.h> "sizof( double) =%d\n", sizeof( double) ); "exp(lgamma) (%.20f)=%.20f\n", x, exp(lgamma(x)));
#include <math.h>
int main ()
{
double x = 6.0;
printf(
printf("lgamma (%.20f)=%.20f\n", x, lgamma(x));
printf("tgamma (%.20f)=%.20f\n", x, tgamma(x));
printf(
return 0;
}
On eglibc you'll see:
sizof(double)=8 00000000) =4.787491742782 04581157 00000000) =119.9999999999 9997157829 00000000) =119.9999999999 9997157829
lgamma (6.000000000000
tgamma (6.000000000000
exp(lgamma) (6.000000000000
On a better libc (MacOSX 10.6 on x86_64):
sizof(double)=8 00000000) =4.787491742782 04492339 00000000) =120.0000000000 0000000000 00000000) =119.9999999999 9987210231
lgamma (6.000000000000
tgamma (6.000000000000
exp(lgamma) (6.000000000000
So eglibc does a bad job computing gamma() on double (i.e. on 8-byte floats) in this range of values (around 6.0).