Re: [musl] expl function imprecise on arm64 and s390x

[Rich Felker <dalias@libc.org>]

On Mon, Jan 25, 2021 at 09:55:40PM +0100, Bruno Haible wrote:
> The expl() function in musl libc 1.2.2, on arm64 and s390x, returns results
> that are far from as precise as the 'long double' type.
> 
> Test program:
> ===============================================================================
> #ifndef __NO_MATH_INLINES
> # define __NO_MATH_INLINES 1 /* for glibc */
> #endif
> #include <float.h>
> #include <math.h>
> #include <stdio.h>
> #undef expl
> extern
> #ifdef __cplusplus
> "C"
> #endif
> long double expl (long double);
> static long double dummy (long double x) { return 0; }
> int main (int argc, char *argv[])
> {
>   long double (* volatile my_expl) (long double) = argc ? expl : dummy;
>   int result = 0;
>   {
>     const long double TWO_LDBL_MANT_DIG = /* 2^LDBL_MANT_DIG */
>       (long double) (1U << ((LDBL_MANT_DIG - 1) / 5))
>       * (long double) (1U << ((LDBL_MANT_DIG - 1 + 1) / 5))
>       * (long double) (1U << ((LDBL_MANT_DIG - 1 + 2) / 5))
>       * (long double) (1U << ((LDBL_MANT_DIG - 1 + 3) / 5))
>       * (long double) (1U << ((LDBL_MANT_DIG - 1 + 4) / 5));
>     long double x = 11.358L;
>     long double y = my_expl (x);
>     long double z = my_expl (- x);
>     long double err = (y * z - 1.0L) * TWO_LDBL_MANT_DIG;
>     printf ("LDBL_MANT_DIG = %d\ny = %.*Lg\nz = %.*Lg\nerr = %g\n",
>             LDBL_MANT_DIG,
>             (int)(LDBL_MANT_DIG * 0.30103), y,
>             (int)(LDBL_MANT_DIG * 0.30103), z,
>             (double)err);
>     if (!(err >= -100.0L && err <= 100.0L))
>       result |= 4;
>   }
>   return result;
> }
> ===============================================================================
> $ gcc -Wall foo.c -lm
> 
> Precise values:
> y = 85647.901279331790624928290655026607849952741287...
> z = 1.16757093292759527899833485690336857485297224766...e-5
> 
> Output on glibc/x86_64:
> LDBL_MANT_DIG = 64
> y = 85647.90127933179059
> z = 1.167570932927595279e-05
> err = -0.5
> 
> Output on musl libc/arm64:
> LDBL_MANT_DIG = 113
> y = 85647.90127933183975983411073684692
> z = 1.167570932927594569211357522497963e-05
> err = -1.77747e+17
> 
> Output on musl libc/s390x:
> LDBL_MANT_DIG = 113
> y = 85647.90127933183975983411073684692
> z = 1.167570932927594569211357522497963e-05
> err = -1.77747e+17
> 
> You can see that only the first ca. 15 decimal digits are correct.
> 
> Numerical approximation of transcendental functions is a long-studied
> domain of numerical analysis. A recent introduction to the field, for
> the exponential function, is [1].
> 
> Bruno
> 
> [1] https://justinwillmert.com/articles/2020/numerically-computing-the-exponential-function-with-polynomial-approximations/

Thanks. I suspect you may find others; it's a known issue that not all
of the long double functions are appropriately accurate for archs
where long double is quad. sqrtl (a serious error because it's
required to be correctly rounded) was fixed in the last release cycle.
Hopefully we'll get the rest in better shape soon.

Rich