[<prev] [next>] [thread-next>] [day] [month] [year] [list]
Message-ID: <1726645.2jRaziT4Rs@omega>
Date: Mon, 25 Jan 2021 21:55:40 +0100
From: Bruno Haible <bruno@...sp.org>
To: musl@...ts.openwall.com
Subject: expl function imprecise on arm64 and s390x
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/
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
