[<prev] [next>] [<thread-prev] [day] [month] [year] [list]
Date: Sun, 31 Jan 2021 13:36:59 +0100
From: Bruno Haible <bruno@...sp.org>
To: Rich Felker <dalias@...c.org>
Cc: musl@...ts.openwall.com
Subject: expm1l, logl, log1pl, log2l, log10l, remainderl functions imprecise on arm64 and s390x
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
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.
