* [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).