[musl] expl function imprecise on arm64 and s390x

mailing list of musl libc
 help / color / mirror / code / Atom feed

* [musl] expl function imprecise on arm64 and s390x
@ 2021-01-25 20:55 Bruno Haible
  2021-01-25 22:31 ` Rich Felker
  0 siblings, 1 reply; 3+ messages in thread
From: Bruno Haible @ 2021-01-25 20:55 UTC (permalink / raw)
  To: musl

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/


^ permalink raw reply	[flat|nested] 3+ messages in thread

* Re: [musl] expl function imprecise on arm64 and s390x
  2021-01-25 20:55 [musl] expl function imprecise on arm64 and s390x Bruno Haible
@ 2021-01-25 22:31 ` Rich Felker
  2021-01-31 12:36   ` [musl] expm1l, logl, log1pl, log2l, log10l, remainderl functions " Bruno Haible
  0 siblings, 1 reply; 3+ messages in thread
From: Rich Felker @ 2021-01-25 22:31 UTC (permalink / raw)
  To: Bruno Haible; +Cc: musl

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

^ permalink raw reply	[flat|nested] 3+ messages in thread

* [musl] expm1l, logl, log1pl, log2l, log10l, remainderl functions imprecise on arm64 and s390x
  2021-01-25 22:31 ` Rich Felker
@ 2021-01-31 12:36   ` Bruno Haible
  0 siblings, 0 replies; 3+ messages in thread
From: Bruno Haible @ 2021-01-31 12:36 UTC (permalink / raw)
  To: Rich Felker; +Cc: musl

Rich Felker wrote:
> 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. ...
> Hopefully we'll get the rest in better shape soon.

Indeed. Here are more test cases. Hope you can use them. I computed the
"correct values" using the arbitrary-precision arithmetic of GNU clisp.


Test case for expm1l:
===============================================================================
#ifndef __NO_MATH_INLINES
# define __NO_MATH_INLINES 1 /* for glibc */
#endif
#include <float.h>
#include <math.h>
#include <stdio.h>
#undef expm1l
extern
#ifdef __cplusplus
"C"
#endif
long double expm1l (long double);
static long double dummy (long double x) { return 0; }
int main (int argc, char *argv[])
{
  long double (* volatile my_expm1l) (long double) = argc ? expm1l : 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_expm1l (x);
    long double z = my_expm1l (- x);
    volatile long double t = (1.0L + y) * z;
    long double err = (y + t) * 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 |= 1;
  }
  return result;
}
===============================================================================
Output on musl/arm64 and musl/s390x:

LDBL_MANT_DIG = 113
y = 85646.90127933183975983411073684692
z = -0.9999883242906707492281270788225811
err = -1.11952e+22

Correct values:
y = 85646.9012793317906249282906550266078499527412871
z = -0.999988324290670724047210016651430966314251470277


Test case for logl:
===============================================================================
#ifndef __NO_MATH_INLINES
# define __NO_MATH_INLINES 1 /* for glibc */
#endif
#include <float.h>
#include <math.h>
#include <stdio.h>
#undef logl
extern
#ifdef __cplusplus
"C"
#endif
long double logl (long double);
static long double dummy (long double x) { return 0; }
int main (int argc, char *argv[])
{
  long double (* volatile my_logl) (long double) = argc ? logl : 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 = 16.981137113807045L;
    long double y = my_logl (x);
    long double z = my_logl (1.0L / x);
    long double err = (y + z) * 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 |= 1;
  }

  return result;
}
===============================================================================
Output on musl/arm64 and musl/s390x:
LDBL_MANT_DIG = 113
y = 2.832103146474115096253854062524624
z = -2.83210314647411554034306391258724
err = -2.30584e+18

Correct values:
y = 2.83210314647411532843886058554627196598416729434
z = -2.83210314647411532843886058554627196598416729434


Test case for log1pl:
===============================================================================
#ifndef __NO_MATH_INLINES
# define __NO_MATH_INLINES 1 /* for glibc */
#endif
#include <float.h>
#include <math.h>
#include <stdio.h>
#undef log1pl
extern
#ifdef __cplusplus
"C"
#endif
long double log1pl (long double);
static long double dummy (long double x) { return 0; }
int main (int argc, char *argv[])
{
  long double (* volatile my_log1pl) (long double) = argc ? log1pl : 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_log1pl (x);
    long double z = my_log1pl (- x / (1.0L + x));
    long double err = (y + z) * 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 >= -900.0L && err <= 900.0L))
      result |= 1;
  }

  return result;
}
===============================================================================
Output on musl/arm64 and musl/s390x:
LDBL_MANT_DIG = 113
y = 2.514303626638787925173801340861246
z = -2.51430362663878748108459149079863
err = 2.30584e+18

Correct values:
y = 2.514303626638787808991104205181310474969595858254
z = -2.514303626638787808991104205181310474969595858254


Test case for log2l:
===============================================================================
#ifndef __NO_MATH_INLINES
# define __NO_MATH_INLINES 1 /* for glibc */
#endif
#include <float.h>
#include <math.h>
#include <stdio.h>
#undef log2l
extern
#ifdef __cplusplus
"C"
#endif
long double log2l (long double);
static long double dummy (long double x) { return 0; }
int main (int argc, char *argv[])
{
  long double (* volatile my_log2l) (long double) = argc ? log2l : 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_log2l (x);
    long double z = my_log2l (1.0L / x);
    long double err = (y + z) * 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 >= -10000.0L && err <= 10000.0L))
      result |= 1;
  }
  if (DBL_MAX_EXP < LDBL_MAX_EXP)
    {
      long double x = ldexpl (1.0L, DBL_MAX_EXP);
      long double y = my_log2l (x);
      if (y > 0 && y + y == y) /* infinite? */
        result |= 2;
    }

  return result;
}
===============================================================================
Output and exit code on musl/arm64 and musl/s390x:
LDBL_MANT_DIG = 113
y = 3.50563691176277458794174890499562
z = -3.505636911762774143852539054933004
err = 2.30584e+18
3

Correct values:
y = 3.505636911762774324741549890160677088025417867858
z = -3.505636911762774324741549890160677088025417867858


Test case for log10l:
===============================================================================
#ifndef __NO_MATH_INLINES
# define __NO_MATH_INLINES 1 /* for glibc */
#endif
#include <float.h>
#include <math.h>
#include <stdio.h>
#undef log10l
extern
#ifdef __cplusplus
"C"
#endif
long double log10l (long double);
static long double dummy (long double x) { return 0; }
int main (int argc, char *argv[])
{
  long double (* volatile my_log10l) (long double) = argc ? log10l : 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 = 7.90097792256024576L;
    long double y = my_log10l (x);
    long double z = my_log10l (1.0L / x);
    long double err = (y + z) * 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 |= 1;
  }

  return result;
}
===============================================================================
Output on musl/arm64 and musl/s390x:
LDBL_MANT_DIG = 113
y = 0.8976808482634930363985859003150836
z = -0.8976808482634929253762834377994295
err = 5.76461e+17

Correct values:
y = 0.897680848263492989548770123475961808070106674916
z = -0.897680848263492989548770123475961808070106674916


Test case for remainderl:
===============================================================================
#ifndef __NO_MATH_INLINES
# define __NO_MATH_INLINES 1 /* for glibc */
#endif
#include <float.h>
#include <math.h>
#include <stdio.h>
#undef remainderl
extern
#ifdef __cplusplus
"C"
#endif
long double remainderl (long double, long double);
static long double dummy (long double x, long double y) { return 0; }
int main (int argc, char *argv[])
{
  long double (* volatile my_remainderl) (long double, long double) = argc ? remainderl : 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.357996098166760874793827872740874983168L;
    long double y = 0.486497838502717923110029188864352615388L;
    long double z = my_remainderl (x, y) - (x - 23 * y);
    long double err = z * TWO_LDBL_MANT_DIG;
    printf ("LDBL_MANT_DIG = %d\nz = %.*Lg\nerr = %g\n",
            LDBL_MANT_DIG,
            (int)(LDBL_MANT_DIG * 0.30103), z,
            (double)err);
    if (!(err >= -32.0L && err <= 32.0L))
      result |= 1;
  }

  return result;
}
===============================================================================
Output on musl/arm64 and musl/s390x:
LDBL_MANT_DIG = 113
z = -6.067178037652900601980525052256813e-15
err = -3.15026e+19

Correct value: z = 0



^ permalink raw reply	[flat|nested] 3+ messages in thread

end of thread, other threads:[~2021-01-31 12:37 UTC | newest]

Thread overview: 3+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2021-01-25 20:55 [musl] expl function imprecise on arm64 and s390x Bruno Haible
2021-01-25 22:31 ` Rich Felker
2021-01-31 12:36   ` [musl] expm1l, logl, log1pl, log2l, log10l, remainderl functions " Bruno Haible

Code repositories for project(s) associated with this public inbox

	https://git.vuxu.org/mirror/musl/

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).