mailing list of musl libc
 help / color / mirror / code / Atom feed
From: Bruno Haible <bruno@clisp.org>
To: musl@lists.openwall.com
Subject: [musl] expl function imprecise on arm64 and s390x
Date: Mon, 25 Jan 2021 21:55:40 +0100	[thread overview]
Message-ID: <1726645.2jRaziT4Rs@omega> (raw)

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/


             reply	other threads:[~2021-01-25 20:55 UTC|newest]

Thread overview: 3+ messages / expand[flat|nested]  mbox.gz  Atom feed  top
2021-01-25 20:55 Bruno Haible [this message]
2021-01-25 22:31 ` Rich Felker
2021-01-31 12:36   ` [musl] expm1l, logl, log1pl, log2l, log10l, remainderl functions " Bruno Haible

Reply instructions:

You may reply publicly to this message via plain-text email
using any one of the following methods:

* Save the following mbox file, import it into your mail client,
  and reply-to-all from there: mbox

  Avoid top-posting and favor interleaved quoting:
  https://en.wikipedia.org/wiki/Posting_style#Interleaved_style

* Reply using the --to, --cc, and --in-reply-to
  switches of git-send-email(1):

  git send-email \
    --in-reply-to=1726645.2jRaziT4Rs@omega \
    --to=bruno@clisp.org \
    --cc=musl@lists.openwall.com \
    /path/to/YOUR_REPLY

  https://kernel.org/pub/software/scm/git/docs/git-send-email.html

* If your mail client supports setting the In-Reply-To header
  via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line before the message body.
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).