Improve the accuracy of tgamma (BZ #26983)

With this patch, the maximal known error for tgamma is now reduced to 9 ulps
for dbl-64, for all rounding modes. Since exhaustive testing is not possible
for dbl-64, it might be that there are still cases with an error larger than
9 ulps, but all known cases are fixed (intensive tests were done to find cases
with large errors).

Tested on x86_64 and powerpc.
---
 math/auto-libm-test-in                  |   6 +-
 math/auto-libm-test-out-tgamma          | 276 ++++++++++++++++++++++++
 math/mul_split.h                        |  74 +++++++
 sysdeps/ieee754/dbl-64/e_gamma_r.c      |  39 +++-
 sysdeps/ieee754/dbl-64/gamma_product.c  |   4 +
 sysdeps/ieee754/ldbl-96/gamma_product.c |   5 +
 sysdeps/x86_64/fpu/libm-test-ulps       |   6 +-
 7 files changed, 395 insertions(+), 15 deletions(-)

On 26/03/2021 12:34, Paul Zimmermann wrote:
> With this patch, the maximal known error for tgamma is now reduced to 9 ulps
> for dbl-64, for all rounding modes. Since exhaustive testing is not possible
> for dbl-64, it might be that there are still cases with an error larger than
> 9 ulps, but all known cases are fixed (intensive tests were done to find cases
> with large errors).
> 
> Tested on x86_64 and powerpc.

I have also check on aarch64, arm, s390x, sparc, and i686 and saw no
regressions (besides the libm-test-ulps update). Some comments below.

> ---
>  math/auto-libm-test-in                  |   6 +-
>  math/auto-libm-test-out-tgamma          | 276 ++++++++++++++++++++++++
>  math/mul_split.h                        |  74 +++++++
>  sysdeps/ieee754/dbl-64/e_gamma_r.c      |  39 +++-
>  sysdeps/ieee754/dbl-64/gamma_product.c  |   4 +
>  sysdeps/ieee754/ldbl-96/gamma_product.c |   5 +
>  sysdeps/x86_64/fpu/libm-test-ulps       |   6 +-
>  7 files changed, 395 insertions(+), 15 deletions(-)
> 
> diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
> index 4edaaa8ee1..af4af15cf4 100644
> --- a/math/auto-libm-test-in
> +++ b/math/auto-libm-test-in
> @@ -8245,8 +8245,12 @@ tgamma -0x6.ce9158p+0
>  tgamma -0xd.cbf53d0e7d06p+0
>  # the next value generates larger error bounds on x86_64 (binary32)
>  tgamma -0x3.0aa534p+0
> -# the next value generates larger error bounds on x86_64 (binary64)
> +# the next values generate large error bounds on x86_64 (binary64)
>  tgamma -0x1.62b8c36c7180bp+4
> +tgamma -0x1.62c4d519e8677p+3
> +tgamma -0x1.c033cc426752fp+2
> +tgamma -0x1.62cfd0d34ade2p+3
> +tgamma -0x1.8814da6eb7dbp+5
>  
>  y0 0.125
>  y0 0.75

Ok.

> diff --git a/math/auto-libm-test-out-tgamma b/math/auto-libm-test-out-tgamma
> index 032ad3a7d7..20a9916afa 100644
> --- a/math/auto-libm-test-out-tgamma
> +++ b/math/auto-libm-test-out-tgamma
> @@ -27397,3 +27397,279 @@ tgamma -0x1.62b8c36c7180bp+4
>  = tgamma tonearest ibm128 -0x1.62b8c36c7180bp+4 : -0xf.3fb73d327d6a801abb3039b104p-72 : inexact-ok
>  = tgamma towardzero ibm128 -0x1.62b8c36c7180bp+4 : -0xf.3fb73d327d6a801abb3039b1p-72 : inexact-ok
>  = tgamma upward ibm128 -0x1.62b8c36c7180bp+4 : -0xf.3fb73d327d6a801abb3039b1p-72 : inexact-ok
> +tgamma -0x1.62c4d519e8677p+3
> += tgamma downward binary32 -0xb.1626ap+0 : 0x3.fac968p-24 : inexact-ok
> += tgamma tonearest binary32 -0xb.1626ap+0 : 0x3.fac96cp-24 : inexact-ok
> += tgamma towardzero binary32 -0xb.1626ap+0 : 0x3.fac968p-24 : inexact-ok
> += tgamma upward binary32 -0xb.1626ap+0 : 0x3.fac96cp-24 : inexact-ok
> += tgamma downward binary64 -0xb.1626ap+0 : 0x3.fac96a100e60ep-24 : inexact-ok
> += tgamma tonearest binary64 -0xb.1626ap+0 : 0x3.fac96a100e61p-24 : inexact-ok
> += tgamma towardzero binary64 -0xb.1626ap+0 : 0x3.fac96a100e60ep-24 : inexact-ok
> += tgamma upward binary64 -0xb.1626ap+0 : 0x3.fac96a100e61p-24 : inexact-ok
> += tgamma downward intel96 -0xb.1626ap+0 : 0x3.fac96a100e60f1f8p-24 : inexact-ok
> += tgamma tonearest intel96 -0xb.1626ap+0 : 0x3.fac96a100e60f1fcp-24 : inexact-ok
> += tgamma towardzero intel96 -0xb.1626ap+0 : 0x3.fac96a100e60f1f8p-24 : inexact-ok
> += tgamma upward intel96 -0xb.1626ap+0 : 0x3.fac96a100e60f1fcp-24 : inexact-ok
> += tgamma downward m68k96 -0xb.1626ap+0 : 0x3.fac96a100e60f1f8p-24 : inexact-ok
> += tgamma tonearest m68k96 -0xb.1626ap+0 : 0x3.fac96a100e60f1fcp-24 : inexact-ok
> += tgamma towardzero m68k96 -0xb.1626ap+0 : 0x3.fac96a100e60f1f8p-24 : inexact-ok
> += tgamma upward m68k96 -0xb.1626ap+0 : 0x3.fac96a100e60f1fcp-24 : inexact-ok
> += tgamma downward binary128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfb6p-24 : inexact-ok
> += tgamma tonearest binary128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfb8p-24 : inexact-ok
> += tgamma towardzero binary128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfb6p-24 : inexact-ok
> += tgamma upward binary128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfb8p-24 : inexact-ok
> += tgamma downward ibm128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfp-24 : inexact-ok
> += tgamma tonearest ibm128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcep-24 : inexact-ok
> += tgamma towardzero ibm128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfp-24 : inexact-ok
> += tgamma upward ibm128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcep-24 : inexact-ok
> += tgamma downward binary32 -0xb.1626bp+0 : 0x3.fac6p-24 : inexact-ok
> += tgamma tonearest binary32 -0xb.1626bp+0 : 0x3.fac6p-24 : inexact-ok
> += tgamma towardzero binary32 -0xb.1626bp+0 : 0x3.fac6p-24 : inexact-ok
> += tgamma upward binary32 -0xb.1626bp+0 : 0x3.fac604p-24 : inexact-ok
> += tgamma downward binary64 -0xb.1626bp+0 : 0x3.fac600633d74p-24 : inexact-ok
> += tgamma tonearest binary64 -0xb.1626bp+0 : 0x3.fac600633d742p-24 : inexact-ok
> += tgamma towardzero binary64 -0xb.1626bp+0 : 0x3.fac600633d74p-24 : inexact-ok
> += tgamma upward binary64 -0xb.1626bp+0 : 0x3.fac600633d742p-24 : inexact-ok
> += tgamma downward intel96 -0xb.1626bp+0 : 0x3.fac600633d7414fp-24 : inexact-ok
> += tgamma tonearest intel96 -0xb.1626bp+0 : 0x3.fac600633d7414f4p-24 : inexact-ok
> += tgamma towardzero intel96 -0xb.1626bp+0 : 0x3.fac600633d7414fp-24 : inexact-ok
> += tgamma upward intel96 -0xb.1626bp+0 : 0x3.fac600633d7414f4p-24 : inexact-ok
> += tgamma downward m68k96 -0xb.1626bp+0 : 0x3.fac600633d7414fp-24 : inexact-ok
> += tgamma tonearest m68k96 -0xb.1626bp+0 : 0x3.fac600633d7414f4p-24 : inexact-ok
> += tgamma towardzero m68k96 -0xb.1626bp+0 : 0x3.fac600633d7414fp-24 : inexact-ok
> += tgamma upward m68k96 -0xb.1626bp+0 : 0x3.fac600633d7414f4p-24 : inexact-ok
> += tgamma downward binary128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d97364607cp-24 : inexact-ok
> += tgamma tonearest binary128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d97364607cp-24 : inexact-ok
> += tgamma towardzero binary128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d97364607cp-24 : inexact-ok
> += tgamma upward binary128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d97364607ep-24 : inexact-ok
> += tgamma downward ibm128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d973646p-24 : inexact-ok
> += tgamma tonearest ibm128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d973646p-24 : inexact-ok
> += tgamma towardzero ibm128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d973646p-24 : inexact-ok
> += tgamma upward ibm128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d9736461p-24 : inexact-ok
> += tgamma downward binary64 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7ep-24 : inexact-ok
> += tgamma tonearest binary64 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7ep-24 : inexact-ok
> += tgamma towardzero binary64 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7ep-24 : inexact-ok
> += tgamma upward binary64 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d8p-24 : inexact-ok
> += tgamma downward intel96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e69p-24 : inexact-ok
> += tgamma tonearest intel96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e694p-24 : inexact-ok
> += tgamma towardzero intel96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e69p-24 : inexact-ok
> += tgamma upward intel96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e694p-24 : inexact-ok
> += tgamma downward m68k96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e69p-24 : inexact-ok
> += tgamma tonearest m68k96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e694p-24 : inexact-ok
> += tgamma towardzero m68k96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e69p-24 : inexact-ok
> += tgamma upward m68k96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e694p-24 : inexact-ok
> += tgamma downward binary128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429ecp-24 : inexact-ok
> += tgamma tonearest binary128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429ecp-24 : inexact-ok
> += tgamma towardzero binary128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429ecp-24 : inexact-ok
> += tgamma upward binary128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429eep-24 : inexact-ok
> += tgamma downward ibm128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429p-24 : inexact-ok
> += tgamma tonearest ibm128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b13942ap-24 : inexact-ok
> += tgamma towardzero ibm128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429p-24 : inexact-ok
> += tgamma upward ibm128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b13942ap-24 : inexact-ok
> +tgamma -0x1.c033cc426752fp+2
> += tgamma downward binary32 -0x7.00cf3p+0 : 0xf.f7014p-8 : inexact-ok
> += tgamma tonearest binary32 -0x7.00cf3p+0 : 0xf.f7014p-8 : inexact-ok
> += tgamma towardzero binary32 -0x7.00cf3p+0 : 0xf.f7014p-8 : inexact-ok
> += tgamma upward binary32 -0x7.00cf3p+0 : 0xf.f7015p-8 : inexact-ok
> += tgamma downward binary64 -0x7.00cf3p+0 : 0xf.f7014642d9848p-8 : inexact-ok
> += tgamma tonearest binary64 -0x7.00cf3p+0 : 0xf.f7014642d9848p-8 : inexact-ok
> += tgamma towardzero binary64 -0x7.00cf3p+0 : 0xf.f7014642d9848p-8 : inexact-ok
> += tgamma upward binary64 -0x7.00cf3p+0 : 0xf.f7014642d985p-8 : inexact-ok
> += tgamma downward intel96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok
> += tgamma tonearest intel96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok
> += tgamma towardzero intel96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok
> += tgamma upward intel96 -0x7.00cf3p+0 : 0xf.f7014642d9849fdp-8 : inexact-ok
> += tgamma downward m68k96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok
> += tgamma tonearest m68k96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok
> += tgamma towardzero m68k96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok
> += tgamma upward m68k96 -0x7.00cf3p+0 : 0xf.f7014642d9849fdp-8 : inexact-ok
> += tgamma downward binary128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2ae8p-8 : inexact-ok
> += tgamma tonearest binary128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2afp-8 : inexact-ok
> += tgamma towardzero binary128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2ae8p-8 : inexact-ok
> += tgamma upward binary128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2afp-8 : inexact-ok
> += tgamma downward ibm128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea28p-8 : inexact-ok
> += tgamma tonearest ibm128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2cp-8 : inexact-ok
> += tgamma towardzero ibm128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea28p-8 : inexact-ok
> += tgamma upward ibm128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2cp-8 : inexact-ok
> += tgamma downward binary32 -0x7.00cf38p+0 : 0xf.f6627p-8 : inexact-ok
> += tgamma tonearest binary32 -0x7.00cf38p+0 : 0xf.f6628p-8 : inexact-ok
> += tgamma towardzero binary32 -0x7.00cf38p+0 : 0xf.f6627p-8 : inexact-ok
> += tgamma upward binary32 -0x7.00cf38p+0 : 0xf.f6628p-8 : inexact-ok
> += tgamma downward binary64 -0x7.00cf38p+0 : 0xf.f6627d248f8ep-8 : inexact-ok
> += tgamma tonearest binary64 -0x7.00cf38p+0 : 0xf.f6627d248f8ep-8 : inexact-ok
> += tgamma towardzero binary64 -0x7.00cf38p+0 : 0xf.f6627d248f8ep-8 : inexact-ok
> += tgamma upward binary64 -0x7.00cf38p+0 : 0xf.f6627d248f8e8p-8 : inexact-ok
> += tgamma downward intel96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25dp-8 : inexact-ok
> += tgamma tonearest intel96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ep-8 : inexact-ok
> += tgamma towardzero intel96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25dp-8 : inexact-ok
> += tgamma upward intel96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ep-8 : inexact-ok
> += tgamma downward m68k96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25dp-8 : inexact-ok
> += tgamma tonearest m68k96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ep-8 : inexact-ok
> += tgamma towardzero m68k96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25dp-8 : inexact-ok
> += tgamma upward m68k96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ep-8 : inexact-ok
> += tgamma downward binary128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcbp-8 : inexact-ok
> += tgamma tonearest binary128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcb8p-8 : inexact-ok
> += tgamma towardzero binary128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcbp-8 : inexact-ok
> += tgamma upward binary128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcb8p-8 : inexact-ok
> += tgamma downward ibm128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcp-8 : inexact-ok
> += tgamma tonearest ibm128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcp-8 : inexact-ok
> += tgamma towardzero ibm128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcp-8 : inexact-ok
> += tgamma upward ibm128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81586p-8 : inexact-ok
> += tgamma downward binary64 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad9999808p-8 : inexact-ok
> += tgamma tonearest binary64 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999981p-8 : inexact-ok
> += tgamma towardzero binary64 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad9999808p-8 : inexact-ok
> += tgamma upward binary64 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999981p-8 : inexact-ok
> += tgamma downward intel96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok
> += tgamma tonearest intel96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok
> += tgamma towardzero intel96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok
> += tgamma upward intel96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4bp-8 : inexact-ok
> += tgamma downward m68k96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok
> += tgamma tonearest m68k96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok
> += tgamma towardzero m68k96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok
> += tgamma upward m68k96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4bp-8 : inexact-ok
> += tgamma downward binary128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814ecp-8 : inexact-ok
> += tgamma tonearest binary128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814ecp-8 : inexact-ok
> += tgamma towardzero binary128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814ecp-8 : inexact-ok
> += tgamma upward binary128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814ec8p-8 : inexact-ok
> += tgamma downward ibm128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814cp-8 : inexact-ok
> += tgamma tonearest ibm128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608815p-8 : inexact-ok
> += tgamma towardzero ibm128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814cp-8 : inexact-ok
> += tgamma upward ibm128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608815p-8 : inexact-ok
> +tgamma -0x1.62cfd0d34ade2p+3
> += tgamma downward binary32 -0xb.167e8p+0 : 0x3.e85604p-24 : inexact-ok
> += tgamma tonearest binary32 -0xb.167e8p+0 : 0x3.e85608p-24 : inexact-ok
> += tgamma towardzero binary32 -0xb.167e8p+0 : 0x3.e85604p-24 : inexact-ok
> += tgamma upward binary32 -0xb.167e8p+0 : 0x3.e85608p-24 : inexact-ok
> += tgamma downward binary64 -0xb.167e8p+0 : 0x3.e85606c693e9ap-24 : inexact-ok
> += tgamma tonearest binary64 -0xb.167e8p+0 : 0x3.e85606c693e9cp-24 : inexact-ok
> += tgamma towardzero binary64 -0xb.167e8p+0 : 0x3.e85606c693e9ap-24 : inexact-ok
> += tgamma upward binary64 -0xb.167e8p+0 : 0x3.e85606c693e9cp-24 : inexact-ok
> += tgamma downward intel96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok
> += tgamma tonearest intel96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok
> += tgamma towardzero intel96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok
> += tgamma upward intel96 -0xb.167e8p+0 : 0x3.e85606c693e9bf28p-24 : inexact-ok
> += tgamma downward m68k96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok
> += tgamma tonearest m68k96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok
> += tgamma towardzero m68k96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok
> += tgamma upward m68k96 -0xb.167e8p+0 : 0x3.e85606c693e9bf28p-24 : inexact-ok
> += tgamma downward binary128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8c3ap-24 : inexact-ok
> += tgamma tonearest binary128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8c3cp-24 : inexact-ok
> += tgamma towardzero binary128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8c3ap-24 : inexact-ok
> += tgamma upward binary128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8c3cp-24 : inexact-ok
> += tgamma downward ibm128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8cp-24 : inexact-ok
> += tgamma tonearest ibm128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8cp-24 : inexact-ok
> += tgamma towardzero ibm128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8cp-24 : inexact-ok
> += tgamma upward ibm128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8dp-24 : inexact-ok
> += tgamma downward binary32 -0xb.167e9p+0 : 0x3.e852b8p-24 : inexact-ok
> += tgamma tonearest binary32 -0xb.167e9p+0 : 0x3.e852b8p-24 : inexact-ok
> += tgamma towardzero binary32 -0xb.167e9p+0 : 0x3.e852b8p-24 : inexact-ok
> += tgamma upward binary32 -0xb.167e9p+0 : 0x3.e852bcp-24 : inexact-ok
> += tgamma downward binary64 -0xb.167e9p+0 : 0x3.e852b838ddf22p-24 : inexact-ok
> += tgamma tonearest binary64 -0xb.167e9p+0 : 0x3.e852b838ddf24p-24 : inexact-ok
> += tgamma towardzero binary64 -0xb.167e9p+0 : 0x3.e852b838ddf22p-24 : inexact-ok
> += tgamma upward binary64 -0xb.167e9p+0 : 0x3.e852b838ddf24p-24 : inexact-ok
> += tgamma downward intel96 -0xb.167e9p+0 : 0x3.e852b838ddf2348cp-24 : inexact-ok
> += tgamma tonearest intel96 -0xb.167e9p+0 : 0x3.e852b838ddf2349p-24 : inexact-ok
> += tgamma towardzero intel96 -0xb.167e9p+0 : 0x3.e852b838ddf2348cp-24 : inexact-ok
> += tgamma upward intel96 -0xb.167e9p+0 : 0x3.e852b838ddf2349p-24 : inexact-ok
> += tgamma downward m68k96 -0xb.167e9p+0 : 0x3.e852b838ddf2348cp-24 : inexact-ok
> += tgamma tonearest m68k96 -0xb.167e9p+0 : 0x3.e852b838ddf2349p-24 : inexact-ok
> += tgamma towardzero m68k96 -0xb.167e9p+0 : 0x3.e852b838ddf2348cp-24 : inexact-ok
> += tgamma upward m68k96 -0xb.167e9p+0 : 0x3.e852b838ddf2349p-24 : inexact-ok
> += tgamma downward binary128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd745531ep-24 : inexact-ok
> += tgamma tonearest binary128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd745531ep-24 : inexact-ok
> += tgamma towardzero binary128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd745531ep-24 : inexact-ok
> += tgamma upward binary128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd745532p-24 : inexact-ok
> += tgamma downward ibm128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd74553p-24 : inexact-ok
> += tgamma tonearest ibm128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd74553p-24 : inexact-ok
> += tgamma towardzero ibm128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd74553p-24 : inexact-ok
> += tgamma upward ibm128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd74554p-24 : inexact-ok
> += tgamma downward binary64 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa6p-24 : inexact-ok
> += tgamma tonearest binary64 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa8p-24 : inexact-ok
> += tgamma towardzero binary64 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa6p-24 : inexact-ok
> += tgamma upward binary64 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa8p-24 : inexact-ok
> += tgamma downward intel96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok
> += tgamma tonearest intel96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok
> += tgamma towardzero intel96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok
> += tgamma upward intel96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e88p-24 : inexact-ok
> += tgamma downward m68k96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok
> += tgamma tonearest m68k96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok
> += tgamma towardzero m68k96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok
> += tgamma upward m68k96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e88p-24 : inexact-ok
> += tgamma downward binary128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e7829836ap-24 : inexact-ok
> += tgamma tonearest binary128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e7829836ap-24 : inexact-ok
> += tgamma towardzero binary128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e7829836ap-24 : inexact-ok
> += tgamma upward binary128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e7829836cp-24 : inexact-ok
> += tgamma downward ibm128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e782983p-24 : inexact-ok
> += tgamma tonearest ibm128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e782983p-24 : inexact-ok
> += tgamma towardzero ibm128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e782983p-24 : inexact-ok
> += tgamma upward ibm128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e782984p-24 : inexact-ok
> +tgamma -0x1.8814da6eb7dbp+5
> += tgamma downward binary32 -0x3.1029b4p+4 : 0x0p+0 : inexact-ok underflow errno-erange
> += tgamma tonearest binary32 -0x3.1029b4p+4 : 0x0p+0 : inexact-ok underflow errno-erange
> += tgamma towardzero binary32 -0x3.1029b4p+4 : 0x0p+0 : inexact-ok underflow errno-erange
> += tgamma upward binary32 -0x3.1029b4p+4 : 0x8p-152 : inexact-ok underflow errno-erange-ok
> += tgamma downward binary64 -0x3.1029b4p+4 : 0x3.fd9063a22653ep-204 : inexact-ok
> += tgamma tonearest binary64 -0x3.1029b4p+4 : 0x3.fd9063a22653ep-204 : inexact-ok
> += tgamma towardzero binary64 -0x3.1029b4p+4 : 0x3.fd9063a22653ep-204 : inexact-ok
> += tgamma upward binary64 -0x3.1029b4p+4 : 0x3.fd9063a22654p-204 : inexact-ok
> += tgamma downward intel96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok
> += tgamma tonearest intel96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok
> += tgamma towardzero intel96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok
> += tgamma upward intel96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e8p-204 : inexact-ok
> += tgamma downward m68k96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok
> += tgamma tonearest m68k96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok
> += tgamma towardzero m68k96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok
> += tgamma upward m68k96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e8p-204 : inexact-ok
> += tgamma downward binary128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5bcp-204 : inexact-ok
> += tgamma tonearest binary128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5bcp-204 : inexact-ok
> += tgamma towardzero binary128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5bcp-204 : inexact-ok
> += tgamma upward binary128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5bep-204 : inexact-ok
> += tgamma downward ibm128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5p-204 : inexact-ok
> += tgamma tonearest ibm128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a6p-204 : inexact-ok
> += tgamma towardzero ibm128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5p-204 : inexact-ok
> += tgamma upward ibm128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a6p-204 : inexact-ok
> += tgamma downward binary32 -0x3.1029b8p+4 : 0x0p+0 : inexact-ok underflow errno-erange
> += tgamma tonearest binary32 -0x3.1029b8p+4 : 0x0p+0 : inexact-ok underflow errno-erange
> += tgamma towardzero binary32 -0x3.1029b8p+4 : 0x0p+0 : inexact-ok underflow errno-erange
> += tgamma upward binary32 -0x3.1029b8p+4 : 0x8p-152 : inexact-ok underflow errno-erange-ok
> += tgamma downward binary64 -0x3.1029b8p+4 : 0x3.fd2a955e20f56p-204 : inexact-ok
> += tgamma tonearest binary64 -0x3.1029b8p+4 : 0x3.fd2a955e20f58p-204 : inexact-ok
> += tgamma towardzero binary64 -0x3.1029b8p+4 : 0x3.fd2a955e20f56p-204 : inexact-ok
> += tgamma upward binary64 -0x3.1029b8p+4 : 0x3.fd2a955e20f58p-204 : inexact-ok
> += tgamma downward intel96 -0x3.1029b8p+4 : 0x3.fd2a955e20f57628p-204 : inexact-ok
> += tgamma tonearest intel96 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762cp-204 : inexact-ok
> += tgamma towardzero intel96 -0x3.1029b8p+4 : 0x3.fd2a955e20f57628p-204 : inexact-ok
> += tgamma upward intel96 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762cp-204 : inexact-ok
> += tgamma downward m68k96 -0x3.1029b8p+4 : 0x3.fd2a955e20f57628p-204 : inexact-ok
> += tgamma tonearest m68k96 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762cp-204 : inexact-ok
> += tgamma towardzero m68k96 -0x3.1029b8p+4 : 0x3.fd2a955e20f57628p-204 : inexact-ok
> += tgamma upward m68k96 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762cp-204 : inexact-ok
> += tgamma downward binary128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f01716p-204 : inexact-ok
> += tgamma tonearest binary128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f01718p-204 : inexact-ok
> += tgamma towardzero binary128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f01716p-204 : inexact-ok
> += tgamma upward binary128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f01718p-204 : inexact-ok
> += tgamma downward ibm128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f017p-204 : inexact-ok
> += tgamma tonearest ibm128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f017p-204 : inexact-ok
> += tgamma towardzero ibm128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f017p-204 : inexact-ok
> += tgamma upward ibm128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f018p-204 : inexact-ok
> += tgamma downward binary64 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5cp-204 : inexact-ok
> += tgamma tonearest binary64 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5cp-204 : inexact-ok
> += tgamma towardzero binary64 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5cp-204 : inexact-ok
> += tgamma upward binary64 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5ep-204 : inexact-ok
> += tgamma downward intel96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d4p-204 : inexact-ok
> += tgamma tonearest intel96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d8p-204 : inexact-ok
> += tgamma towardzero intel96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d4p-204 : inexact-ok
> += tgamma upward intel96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d8p-204 : inexact-ok
> += tgamma downward m68k96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d4p-204 : inexact-ok
> += tgamma tonearest m68k96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d8p-204 : inexact-ok
> += tgamma towardzero m68k96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d4p-204 : inexact-ok
> += tgamma upward m68k96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d8p-204 : inexact-ok
> += tgamma downward binary128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d82ep-204 : inexact-ok
> += tgamma tonearest binary128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d82ep-204 : inexact-ok
> += tgamma towardzero binary128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d82ep-204 : inexact-ok
> += tgamma upward binary128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d83p-204 : inexact-ok
> += tgamma downward ibm128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d8p-204 : inexact-ok
> += tgamma tonearest ibm128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d8p-204 : inexact-ok
> += tgamma towardzero ibm128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d8p-204 : inexact-ok
> += tgamma upward ibm128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d9p-204 : inexact-ok

Ok.

> diff --git a/math/mul_split.h b/math/mul_split.h
> index 8cc60ec630..24138c8c39 100644
> --- a/math/mul_split.h
> +++ b/math/mul_split.h
> @@ -47,4 +47,78 @@ mul_split (double *hi, double *lo, double x, double y)
>  #endif
>  }
>  
> +/* Define WANT_MUL_EXPANSION to use mul_expansion().
> +   Define WANT_DIV_EXPANSION to use div_expansion().  */
> +

This is error prone and no scalable (it requires define multiple macro for
each new symbol).  Just define the function as static inline or if the
routines are not hotspots and/or large move to its own object (for
this patch I think making then static inline should be better).

> +#if defined(WANT_MUL_EXPANSION) || defined(WANT_DIV_EXPANSION)
> +/* Add a + b exactly, such that *hi + *lo = a + b.
> +   Assumes |a| >= |b| and rounding to nearest.  */
> +static void
> +fast_two_sum (double *hi, double *lo, double a, double b)
> +{
> +  double e;
> +
> +  *hi = a + b;
> +  e = *hi - a; /* exact */
> +  *lo = b - e; /* exact */

Double space before '*/'.

> +  /* Now *hi + *lo = a + b exactly.  */
> +}
> +#endif
> +
> +#ifdef WANT_MUL_EXPANSION
> +/* Multiplication of two floating-point expansions: *hi + *lo is an
> +   approximation of (h1+l1)*(h2+l2), assuming |l1| <= 1/2*ulp(h1)
> +   and |l2| <= 1/2*ulp(h2) and rounding to nearest.  */
> +static void
> +mul_expansion (double *hi, double *lo, double h1, double l1,
> +	       double h2, double l2)
> +{
> +  double r, e;
> +
> +  mul_split (hi, lo, h1, h2);
> +  r = h1 * l2 + h2 * l1;
> +  /* Now add r to (hi,lo) using fast two-sum, where we know |r| < |hi|.  */
> +  fast_two_sum (hi, &e, *hi, r);
> +  *lo -= e;
> +}
> +#endif
> +
> +#ifdef WANT_DIV_EXPANSION
> +/* Calculate X / Y and store the approximate result in *HI + *LO.  It is
> +   assumed that Y is not zero, that no overflow nor underflow occurs, and
> +   rounding is to nearest.  */
> +static void
> +div_split (double *hi, double *lo, double x, double y)
> +{
> +  double a, b;
> +
> +  *hi = x / y;
> +  mul_split (&a, &b, *hi, y);
> +  /* a + b = hi*y, which should be near x.  */
> +  a = x - a; /* huge cancellation */
> +  a = a - b;
> +  /* Now x ~ hi*y + a thus x/y ~ hi + a/y.  */
> +  *lo = a / y;
> +}
> +
> +/* Division of two floating-point expansions: *hi + *lo is an
> +   approximation of (h1+l1)/(h2+l2), assuming |l1| <= 1/2*ulp(h1)
> +   and |l2| <= 1/2*ulp(h2), h2+l2 is not zero, and rounding to nearest.  */
> +static void
> +div_expansion (double *hi, double *lo, double h1, double l1,
> +	       double h2, double l2)
> +{
> +  double r, e;
> +
> +  div_split (hi, lo, h1, h2);
> +  /* (h1+l1)/(h2+l2) ~ h1/h2 + (l1*h2 - l2*h1)/h2^2 */
> +  r = (l1 * h2 - l2 * h1) / (h2 * h2);
> +  /* Now add r to (hi,lo) using fast two-sum, where we know |r| < |hi|.  */
> +  fast_two_sum (hi, &e, *hi, r);
> +  *lo += e;
> +  /* Renormalize since |lo| might be larger than 0.5 ulp(hi).  */
> +  fast_two_sum (hi, lo, *hi, *lo);
> +}
> +#endif
> +
>  #endif /* _MUL_SPLIT_H */
> diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c
> index fe69920c76..7c8dc800cf 100644
> --- a/sysdeps/ieee754/dbl-64/e_gamma_r.c
> +++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c
> @@ -24,6 +24,9 @@
>  #include <math-underflow.h>
>  #include <float.h>
>  #include <libm-alias-finite.h>
> +#define WANT_MUL_EXPANSION
> +#define WANT_DIV_EXPANSION
> +#include <mul_split.h>
>  
>  /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
>     approximation to gamma function.  */
> @@ -88,7 +91,6 @@ gamma_positive (double x, int *exp2_adj)
>  	 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
>  	 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
>  	 factored out.  */
> -      double exp_adj = -eps;
>        double x_adj_int = round (x_adj);
>        double x_adj_frac = x_adj - x_adj_int;
>        int x_adj_log2;
> @@ -99,18 +101,22 @@ gamma_positive (double x, int *exp2_adj)
>  	  x_adj_mant *= 2.0;
>  	}
>        *exp2_adj = x_adj_log2 * (int) x_adj_int;
> -      double ret = (__ieee754_pow (x_adj_mant, x_adj)
> -		    * __ieee754_exp2 (x_adj_log2 * x_adj_frac)
> -		    * __ieee754_exp (-x_adj)
> -		    * sqrt (2 * M_PI / x_adj)
> -		    / prod);
> -      exp_adj += x_eps * __ieee754_log (x_adj);
> +      double h1, l1, h2, l2;
> +      mul_split (&h1, &l1, __ieee754_pow (x_adj_mant, x_adj),
> +			   __ieee754_exp2 (x_adj_log2 * x_adj_frac));
> +      mul_split (&h2, &l2, __ieee754_exp (-x_adj), sqrt (2 * M_PI / x_adj));
> +      mul_expansion (&h1, &l1, h1, l1, h2, l2);
> +      /* Divide by prod * (1 + eps).  */
> +      div_expansion (&h1, &l1, h1, l1, prod, prod * eps);
> +      double exp_adj = x_eps * __ieee754_log (x_adj);
>        double bsum = gamma_coeff[NCOEFF - 1];
>        double x_adj2 = x_adj * x_adj;
>        for (size_t i = 1; i <= NCOEFF - 1; i++)
>  	bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
>        exp_adj += bsum / x_adj;
> -      return ret + ret * __expm1 (exp_adj);
> +      /* Now return (h1+l1) * exp(exp_adj), where exp_adj is small.  */
> +      l1 += h1 * __expm1 (exp_adj);
> +      return h1 + l1;
>      }
>  }
>  

Ok.

> @@ -188,9 +194,20 @@ __ieee754_gamma_r (double x, int *signgamp)
>  			       ? __sin (M_PI * frac)
>  			       : __cos (M_PI * (0.5 - frac)));
>  	      int exp2_adj;
> -	      double tret = M_PI / (-x * sinpix
> -				    * gamma_positive (-x, &exp2_adj));
> -	      ret = __scalbn (tret, -exp2_adj);
> +	      double h1, l1, h2, l2;
> +	      h2 = gamma_positive (-x, &exp2_adj);
> +	      mul_split (&h1, &l1, sinpix, h2);
> +	      /* sinpix*gamma_positive(.) = h1 + l1 */
> +	      mul_split (&h2, &l2, h1, x);
> +	      /* h1*x = h2 + l2 */
> +	      /* (h1 + l1) * x = h1*x + l1*x = h2 + l2 + l1*x */
> +	      l2 += l1 * x;
> +	      /* x*sinpix*gamma_positive(.) ~ h2 + l2 */
> +	      h1 = 0x3.243f6a8885a3p+0;   /* binary64 approximation of Pi */
> +	      l1 = 0x8.d313198a2e038p-56; /* |h1+l1-Pi| < 3e-33 */
> +	      /* Now we divide h1 + l1 by h2 + l2.  */
> +	      div_expansion (&h1, &l1, h1, l1, h2, l2);
> +	      ret = __scalbn (-h1, -exp2_adj);
>  	      math_check_force_underflow_nonneg (ret);
>  	    }
>  	}

Ok.

> diff --git a/sysdeps/ieee754/dbl-64/gamma_product.c b/sysdeps/ieee754/dbl-64/gamma_product.c
> index 643df11b53..2e3b5148dc 100644
> --- a/sysdeps/ieee754/dbl-64/gamma_product.c
> +++ b/sysdeps/ieee754/dbl-64/gamma_product.c
> @@ -16,6 +16,10 @@
>     License along with the GNU C Library; if not, see
>     <https://www.gnu.org/licenses/>.  */
>  
> +/* Note: this function is only used on an architecture that doesn't use
> +   ldbl-96 (i.e., anything other than x86_64, i386, ia64 or m68k).
> +   On other architectures the one in ldbl-96 is used.  */
> +
>  #include <math.h>
>  #include <math_private.h>
>  #include <fenv_private.h>

I think there is no need to add such comment, it is implict from
sysdep path.

> diff --git a/sysdeps/ieee754/ldbl-96/gamma_product.c b/sysdeps/ieee754/ldbl-96/gamma_product.c
> index 9e00e71452..f881649801 100644
> --- a/sysdeps/ieee754/ldbl-96/gamma_product.c
> +++ b/sysdeps/ieee754/ldbl-96/gamma_product.c
> @@ -16,6 +16,11 @@
>     License along with the GNU C Library; if not, see
>     <https://www.gnu.org/licenses/>.  */
>  
> +/* Note: this function is used on an architecture that uses
> +   ldbl-96 (i.e., x86_64, i386, ia64 or m68k), and in particular
> +   it is called from dbl-64/e_gamma_r.c.
> +   On other architectures the one in dbl-64 is used.  */
> +
>  #include <math.h>
>  #include <math-narrow-eval.h>
>  #include <math_private.h>

Same as before.

> diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
> index bd1fa63702..e0a362a347 100644
> --- a/sysdeps/x86_64/fpu/libm-test-ulps
> +++ b/sysdeps/x86_64/fpu/libm-test-ulps
> @@ -1734,16 +1734,16 @@ float128: 4
>  ldouble: 5
>  
>  Function: "tgamma_downward":
> -double: 8
> +double: 9
>  float: 7
>  float128: 5
> -ldouble: 5
> +ldouble: 6
>  
>  Function: "tgamma_towardzero":
>  double: 9
>  float: 7
>  float128: 5
> -ldouble: 5
> +ldouble: 6
>  
>  Function: "tgamma_upward":
>  double: 9
> 

Ok.

Dear Adhemerval,

thank you for your review. I will send a new version addressing your comments.
For the notes in {dbl-64,ldbl-96}/gamma_product.c, I had added them since for
a newcomer like me, it is not trivial to know which version is used on which
target. Anyway I have removed them in the new version.

Paul