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/
next 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).