[armel] Function tgammal has precision issues in version 2.12.1-0ubuntu10.2 on ARM

I forgot to say that it was on an ARM (little-endian) box!

The bug is in fact in lgammal, as tgammal(x) simply takes expl(lgammal(x));
namely, lgammal(6.0) produces 4.78749174278204581157 on armel,
while it should be closer to 4.7874917427820459942477, i.e. armel gives relative error of about 3.8e-17

Status changed to 'Confirmed' because the bug affects multiple users.

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

On eglibc you'll see:

sizof(double)=8
lgamma (6.00000000000000000000)=4.78749174278204581157
tgamma (6.00000000000000000000)=119.99999999999997157829
exp(lgamma) (6.00000000000000000000)=119.99999999999997157829

On a better libc (MacOSX 10.6 on x86_64):

sizof(double)=8
lgamma (6.00000000000000000000)=4.78749174278204492339
tgamma (6.00000000000000000000)=120.00000000000000000000
exp(lgamma) (6.00000000000000000000)=119.99999999999987210231

So eglibc does a bad job computing gamma() on double (i.e. on 8-byte floats) in this range of values (around 6.0).

This is still a problem with 2.13-24ubuntu4.

This seems to be fixed as far as I can tell building a static soft float binary with the Linaro toolchain (which incidentally no longer supports soft float from 4.9 onwards).

It would be great to know in more details how to get the fix mentioned on an armhfl Ubuntu 12.10 (no, we cannot upgrade the system, unfortunately).