problems with libm

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

problems with libm

Moritz Buhl-3
Hi,

while testing arm hardware on OpenBSD I noticed that some floating point
operations cause failures of other tests.
In fact the current libm is incorrect according to the ISO C Std Annex
G. I found this out after porting some FreeBSD lib msun tests. Many edge
cases for complex floating point operations are not covered at all.

My question on this is what I should do about this. Port the FreeBSD
msun library? Ignore the problem? Patch the current implementation?

I attached a test that is currently breaking. There are many more. To
fix these I just copied the FreeBSD file that implements the failing
function. But I am not sure if this is a good approach.

mbuhl

Index: regress/lib/libm/msun/Makefile
===================================================================
RCS file: /cvs/src/regress/lib/libm/msun/Makefile,v
retrieving revision 1.1.1.1
diff -u -p -r1.1.1.1 Makefile
--- regress/lib/libm/msun/Makefile 21 Feb 2019 16:14:03 -0000 1.1.1.1
+++ regress/lib/libm/msun/Makefile 31 May 2019 19:50:32 -0000
@@ -1,6 +1,7 @@
 # $OpenBSD: Makefile,v 1.1.1.1 2019/02/21 16:14:03 bluhm Exp $
 
 TESTS =
+TESTS += cexp_test
 TESTS += conj_test
 TESTS += fenv_test
 TESTS += lrint_test
Index: regress/lib/libm/msun/cexp_test.c
===================================================================
RCS file: regress/lib/libm/msun/cexp_test.c
diff -N regress/lib/libm/msun/cexp_test.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ regress/lib/libm/msun/cexp_test.c 1 Jun 2019 05:40:51 -0000
@@ -0,0 +1,326 @@
+/*-
+ * Copyright (c) 2008-2011 David Schultz <[hidden email]>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/*
+ * Tests for corner cases in cexp*().
+ */
+
+#include <sys/cdefs.h>
+
+#include <sys/param.h>
+
+#include <assert.h>
+#include <complex.h>
+#include <fenv.h>
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "test-utils.h"
+
+#pragma STDC FENV_ACCESS ON
+#pragma STDC CX_LIMITED_RANGE OFF
+
+/*
+ * Test that a function returns the correct value and sets the
+ * exception flags correctly. The exceptmask specifies which
+ * exceptions we should check. We need to be lenient for several
+ * reasons, but mainly because on some architectures it's impossible
+ * to raise FE_OVERFLOW without raising FE_INEXACT. In some cases,
+ * whether cexp() raises an invalid exception is unspecified.
+ *
+ * These are macros instead of functions so that assert provides more
+ * meaningful error messages.
+ *
+ * XXX The volatile here is to avoid gcc's bogus constant folding and work
+ *     around the lack of support for the FENV_ACCESS pragma.
+ */
+#define test(func, z, result, exceptmask, excepts, checksign) do { \
+ volatile long double complex _d = z; \
+ assert(feclearexcept(FE_ALL_EXCEPT) == 0); \
+ assert(cfpequal_cs((func)(_d), (result), (checksign))); \
+ assert(((void)(func), fetestexcept(exceptmask) == (excepts))); \
+} while (0)
+
+#define test_c(t, func, z, result, exceptmask, excepts, checksign) do { \
+        volatile t complex r = result;                                  \
+        test(func, z, r, exceptmask, excepts, checksign);               \
+} while (0)
+
+#define test_f(func, z, result, exceptmask, excepts, checksign) do {    \
+        test_c(float, func, z, result, exceptmask, excepts, checksign); \
+} while (0)
+
+/* Test within a given tolerance. */
+#define test_tol(func, z, result, tol) do { \
+ volatile long double complex _d = z; \
+ assert(cfpequal_tol((func)(_d), (result), (tol), \
+    FPE_ABS_ZERO | CS_BOTH)); \
+} while (0)
+
+/* Test all the functions that compute cexp(x). */
+#define testall(x, result, exceptmask, excepts, checksign) do { \
+ test(cexp, x, result, exceptmask, excepts, checksign); \
+ test_f(cexpf, x, result, exceptmask, excepts, checksign); \
+} while (0)
+
+/*
+ * Test all the functions that compute cexp(x), within a given tolerance.
+ * The tolerance is specified in ulps.
+ */
+#define testall_tol(x, result, tol) do { \
+ test_tol(cexp, x, result, tol * DBL_ULP()); \
+ test_tol(cexpf, x, result, tol * FLT_ULP()); \
+} while (0)
+
+/* Various finite non-zero numbers to test. */
+static const float finites[] =
+{ -42.0e20, -1.0, -1.0e-10, -0.0, 0.0, 1.0e-10, 1.0, 42.0e20 };
+
+
+/* Tests for 0 */
+static void
+test_zero(void)
+{
+
+ /* cexp(0) = 1, no exceptions raised */
+ testall(0.0, 1.0, ALL_STD_EXCEPT, 0, 1);
+ testall(-0.0, 1.0, ALL_STD_EXCEPT, 0, 1);
+ testall(CMPLXL(0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1);
+ testall(CMPLXL(-0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1);
+}
+
+/*
+ * Tests for NaN.  The signs of the results are indeterminate unless the
+ * imaginary part is 0.
+ */
+static void
+test_nan(void)
+{
+ unsigned i;
+
+ /* cexp(x + NaNi) = NaN + NaNi and optionally raises invalid */
+ /* cexp(NaN + yi) = NaN + NaNi and optionally raises invalid (|y|>0) */
+ for (i = 0; i < nitems(finites); i++) {
+ printf("# Run %d..\n", i);
+ testall(CMPLXL(finites[i], NAN), CMPLXL(NAN, NAN),
+ ALL_STD_EXCEPT & ~FE_INVALID, 0, 0);
+ if (finites[i] == 0.0)
+ continue;
+ /* XXX FE_INEXACT shouldn't be raised here */
+ testall(CMPLXL(NAN, finites[i]), CMPLXL(NAN, NAN),
+ ALL_STD_EXCEPT & ~(FE_INVALID | FE_INEXACT), 0, 0);
+ }
+
+ /* cexp(NaN +- 0i) = NaN +- 0i */
+ testall(CMPLXL(NAN, 0.0), CMPLXL(NAN, 0.0), ALL_STD_EXCEPT, 0, 1);
+ testall(CMPLXL(NAN, -0.0), CMPLXL(NAN, -0.0), ALL_STD_EXCEPT, 0, 1);
+
+ /* cexp(inf + NaN i) = inf + nan i */
+ testall(CMPLXL(INFINITY, NAN), CMPLXL(INFINITY, NAN),
+ ALL_STD_EXCEPT, 0, 0);
+ /* cexp(-inf + NaN i) = 0 */
+ testall(CMPLXL(-INFINITY, NAN), CMPLXL(0.0, 0.0),
+ ALL_STD_EXCEPT, 0, 0);
+ /* cexp(NaN + NaN i) = NaN + NaN i */
+ testall(CMPLXL(NAN, NAN), CMPLXL(NAN, NAN),
+ ALL_STD_EXCEPT, 0, 0);
+}
+
+static void
+test_inf(void)
+{
+ unsigned i;
+
+ /* cexp(x + inf i) = NaN + NaNi and raises invalid */
+ for (i = 0; i < nitems(finites); i++) {
+ printf("# Run %d..\n", i);
+ testall(CMPLXL(finites[i], INFINITY), CMPLXL(NAN, NAN),
+ ALL_STD_EXCEPT, FE_INVALID, 1);
+ }
+ /* cexp(-inf + yi) = 0 * (cos(y) + sin(y)i) */
+ /* XXX shouldn't raise an inexact exception */
+ testall(CMPLXL(-INFINITY, M_PI_4), CMPLXL(0.0, 0.0),
+ ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ testall(CMPLXL(-INFINITY, 3 * M_PI_4), CMPLXL(-0.0, 0.0),
+ ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ testall(CMPLXL(-INFINITY, 5 * M_PI_4), CMPLXL(-0.0, -0.0),
+ ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ testall(CMPLXL(-INFINITY, 7 * M_PI_4), CMPLXL(0.0, -0.0),
+ ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ testall(CMPLXL(-INFINITY, 0.0), CMPLXL(0.0, 0.0),
+ ALL_STD_EXCEPT, 0, 1);
+ testall(CMPLXL(-INFINITY, -0.0), CMPLXL(0.0, -0.0),
+ ALL_STD_EXCEPT, 0, 1);
+ /* cexp(inf + yi) = inf * (cos(y) + sin(y)i) (except y=0) */
+ /* XXX shouldn't raise an inexact exception */
+ testall(CMPLXL(INFINITY, M_PI_4), CMPLXL(INFINITY, INFINITY),
+ ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ testall(CMPLXL(INFINITY, 3 * M_PI_4), CMPLXL(-INFINITY, INFINITY),
+ ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ testall(CMPLXL(INFINITY, 5 * M_PI_4), CMPLXL(-INFINITY, -INFINITY),
+ ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ testall(CMPLXL(INFINITY, 7 * M_PI_4), CMPLXL(INFINITY, -INFINITY),
+ ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ /* cexp(inf + 0i) = inf + 0i */
+ testall(CMPLXL(INFINITY, 0.0), CMPLXL(INFINITY, 0.0),
+ ALL_STD_EXCEPT, 0, 1);
+ testall(CMPLXL(INFINITY, -0.0), CMPLXL(INFINITY, -0.0),
+ ALL_STD_EXCEPT, 0, 1);
+}
+
+static void
+test_reals(void)
+{
+ unsigned i;
+
+ for (i = 0; i < nitems(finites); i++) {
+ /* XXX could check exceptions more meticulously */
+ printf("# Run %d..\n", i);
+ test(cexp, CMPLXL(finites[i], 0.0),
+     CMPLXL(exp(finites[i]), 0.0),
+     FE_INVALID | FE_DIVBYZERO, 0, 1);
+ test(cexp, CMPLXL(finites[i], -0.0),
+     CMPLXL(exp(finites[i]), -0.0),
+     FE_INVALID | FE_DIVBYZERO, 0, 1);
+ test_f(cexpf, CMPLXL(finites[i], 0.0),
+     CMPLXL(expf(finites[i]), 0.0),
+     FE_INVALID | FE_DIVBYZERO, 0, 1);
+ test_f(cexpf, CMPLXL(finites[i], -0.0),
+     CMPLXL(expf(finites[i]), -0.0),
+     FE_INVALID | FE_DIVBYZERO, 0, 1);
+ }
+}
+
+static void
+test_imaginaries(void)
+{
+ unsigned i;
+
+ for (i = 0; i < nitems(finites); i++) {
+ printf("# Run %d..\n", i);
+ test(cexp, CMPLXL(0.0, finites[i]),
+     CMPLXL(cos(finites[i]), sin(finites[i])),
+     ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ test(cexp, CMPLXL(-0.0, finites[i]),
+     CMPLXL(cos(finites[i]), sin(finites[i])),
+     ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ test_f(cexpf, CMPLXL(0.0, finites[i]),
+     CMPLXL(cosf(finites[i]), sinf(finites[i])),
+     ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ test_f(cexpf, CMPLXL(-0.0, finites[i]),
+     CMPLXL(cosf(finites[i]), sinf(finites[i])),
+     ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
+ }
+}
+
+static void
+test_small(void)
+{
+ static const double tests[] = {
+     /* csqrt(a + bI) = x + yI */
+     /* a b x y */
+ 1.0, M_PI_4, M_SQRT2 * 0.5 * M_E, M_SQRT2 * 0.5 * M_E,
+ -1.0, M_PI_4, M_SQRT2 * 0.5 / M_E, M_SQRT2 * 0.5 / M_E,
+ 2.0, M_PI_2, 0.0, M_E * M_E,
+ M_LN2, M_PI, -2.0, 0.0,
+ };
+ double a, b;
+ double x, y;
+ unsigned i;
+
+ for (i = 0; i < nitems(tests); i += 4) {
+ printf("# Run %d..\n", i);
+ a = tests[i];
+ b = tests[i + 1];
+ x = tests[i + 2];
+ y = tests[i + 3];
+ test_tol(cexp, CMPLXL(a, b), CMPLXL(x, y), 3 * DBL_ULP());
+
+ /* float doesn't have enough precision to pass these tests */
+ if (x == 0 || y == 0)
+ continue;
+ test_tol(cexpf, CMPLXL(a, b), CMPLXL(x, y), 1 * FLT_ULP());
+        }
+}
+
+/* Test inputs with a real part r that would overflow exp(r). */
+static void
+test_large(void)
+{
+
+ test_tol(cexp, CMPLXL(709.79, 0x1p-1074),
+ CMPLXL(INFINITY, 8.94674309915433533273e-16), DBL_ULP());
+ test_tol(cexp, CMPLXL(1000, 0x1p-1074),
+ CMPLXL(INFINITY, 9.73344457300016401328e+110), DBL_ULP());
+ test_tol(cexp, CMPLXL(1400, 0x1p-1074),
+ CMPLXL(INFINITY, 5.08228858149196559681e+284), DBL_ULP());
+ test_tol(cexp, CMPLXL(900, 0x1.23456789abcdep-1020),
+ CMPLXL(INFINITY, 7.42156649354218408074e+83), DBL_ULP());
+ test_tol(cexp, CMPLXL(1300, 0x1.23456789abcdep-1020),
+ CMPLXL(INFINITY, 3.87514844965996756704e+257), DBL_ULP());
+
+ test_tol(cexpf, CMPLXL(88.73, 0x1p-149),
+ CMPLXL(INFINITY, 4.80265603e-07), 2 * FLT_ULP());
+ test_tol(cexpf, CMPLXL(90, 0x1p-149),
+ CMPLXL(INFINITY, 1.7101492622e-06f), 2 * FLT_ULP());
+ test_tol(cexpf, CMPLXL(192, 0x1p-149),
+ CMPLXL(INFINITY, 3.396809344e+38f), 2 * FLT_ULP());
+ test_tol(cexpf, CMPLXL(120, 0x1.234568p-120),
+ CMPLXL(INFINITY, 1.1163382522e+16f), 2 * FLT_ULP());
+ test_tol(cexpf, CMPLXL(170, 0x1.234568p-120),
+ CMPLXL(INFINITY, 5.7878851079e+37f), 2 * FLT_ULP());
+}
+
+int
+main(void)
+{
+
+ printf("1..7\n");
+
+ test_zero();
+ printf("ok 1 - cexp zero\n");
+
+ test_nan();
+ printf("ok 2 - cexp nan\n");
+
+ test_inf();
+ printf("ok 3 - cexp inf\n");
+
+ test_reals();
+ printf("ok 4 - cexp reals\n");
+
+ test_imaginaries();
+ printf("ok 5 - cexp imaginaries\n");
+
+ test_small();
+ printf("ok 6 - cexp small\n");
+
+ test_large();
+ printf("ok 7 - cexp large\n");
+
+ return (0);
+}

Reply | Threaded
Open this post in threaded view
|

Re: problems with libm

Daniel Dickman


> On Jul 1, 2019, at 2:50 AM, Moritz Buhl <[hidden email]> wrote:
>
> Hi,
>
> while testing arm hardware on OpenBSD I noticed that some floating point
> operations cause failures of other tests.
> In fact the current libm is incorrect according to the ISO C Std Annex
> G. I found this out after porting some FreeBSD lib msun tests. Many edge
> cases for complex floating point operations are not covered at all.

You’re saying that OpenBSD complex functions in libm may have some rough edges? Sounds right to me.

In the past I’ve found running the numpy regress tests to be a decent suite of tests for libm. I remember the last time I looked at the remaining failures they were mostly related to complex math from what I could tell.

I can share my wip update of numpy 1.16.4 if you’d be interesting in these tests on arm. In fact I would be interested in the results from that platform...

>
> My question on this is what I should do about this. Port the FreeBSD
> msun library? Ignore the problem? Patch the current implementation?

Patches to correct the current code would be very welcome. (In my case especially if it improves numpy test results).

>
> I attached a test that is currently breaking. There are many more. To
> fix these I just copied the FreeBSD file that implements the failing
> function. But I am not sure if this is a good approach.

I’ll take a look at this on some of the platforms I have access to, but arm is not a platform I own. In my experience, the failures can be platform dependent so we’ll see if I can reproduce elsewhere or not.

In general it would be useful to send minimal code illustrating failures and then show the outputs on the platform you’re testing.

Here’s an example libm problem report you could use as a template for reporting problems if you’d like:

https://marc.info/?l=openbsd-tech&m=140167886413618&w=2

>
> mbuhl
>
> Index: regress/lib/libm/msun/Makefile
> ===================================================================
> RCS file: /cvs/src/regress/lib/libm/msun/Makefile,v
> retrieving revision 1.1.1.1
> diff -u -p -r1.1.1.1 Makefile
> --- regress/lib/libm/msun/Makefile    21 Feb 2019 16:14:03 -0000    1.1.1.1
> +++ regress/lib/libm/msun/Makefile    31 May 2019 19:50:32 -0000
> @@ -1,6 +1,7 @@
> # $OpenBSD: Makefile,v 1.1.1.1 2019/02/21 16:14:03 bluhm Exp $
>
> TESTS =
> +TESTS +=    cexp_test
> TESTS +=    conj_test
> TESTS +=    fenv_test
> TESTS +=    lrint_test
> Index: regress/lib/libm/msun/cexp_test.c
> ===================================================================
> RCS file: regress/lib/libm/msun/cexp_test.c
> diff -N regress/lib/libm/msun/cexp_test.c
> --- /dev/null    1 Jan 1970 00:00:00 -0000
> +++ regress/lib/libm/msun/cexp_test.c    1 Jun 2019 05:40:51 -0000
> @@ -0,0 +1,326 @@
> +/*-
> + * Copyright (c) 2008-2011 David Schultz <[hidden email]>
> + * All rights reserved.
> + *
> + * Redistribution and use in source and binary forms, with or without
> + * modification, are permitted provided that the following conditions
> + * are met:
> + * 1. Redistributions of source code must retain the above copyright
> + *    notice, this list of conditions and the following disclaimer.
> + * 2. Redistributions in binary form must reproduce the above copyright
> + *    notice, this list of conditions and the following disclaimer in the
> + *    documentation and/or other materials provided with the distribution.
> + *
> + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
> + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
> + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
> + * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
> + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
> + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
> + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
> + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
> + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
> + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
> + * SUCH DAMAGE.
> + */
> +
> +/*
> + * Tests for corner cases in cexp*().
> + */
> +
> +#include <sys/cdefs.h>
> +
> +#include <sys/param.h>
> +
> +#include <assert.h>
> +#include <complex.h>
> +#include <fenv.h>
> +#include <float.h>
> +#include <math.h>
> +#include <stdio.h>
> +
> +#include "test-utils.h"
> +
> +#pragma STDC FENV_ACCESS    ON
> +#pragma    STDC CX_LIMITED_RANGE    OFF
> +
> +/*
> + * Test that a function returns the correct value and sets the
> + * exception flags correctly. The exceptmask specifies which
> + * exceptions we should check. We need to be lenient for several
> + * reasons, but mainly because on some architectures it's impossible
> + * to raise FE_OVERFLOW without raising FE_INEXACT. In some cases,
> + * whether cexp() raises an invalid exception is unspecified.
> + *
> + * These are macros instead of functions so that assert provides more
> + * meaningful error messages.
> + *
> + * XXX The volatile here is to avoid gcc's bogus constant folding and work
> + *     around the lack of support for the FENV_ACCESS pragma.
> + */
> +#define    test(func, z, result, exceptmask, excepts, checksign)    do {    \
> +    volatile long double complex _d = z;                \
> +    assert(feclearexcept(FE_ALL_EXCEPT) == 0);            \
> +    assert(cfpequal_cs((func)(_d), (result), (checksign)));        \
> +    assert(((void)(func), fetestexcept(exceptmask) == (excepts)));    \
> +} while (0)
> +
> +#define test_c(t, func, z, result, exceptmask, excepts, checksign) do { \
> +        volatile t complex r = result;                                  \
> +        test(func, z, r, exceptmask, excepts, checksign);               \
> +} while (0)
> +
> +#define test_f(func, z, result, exceptmask, excepts, checksign) do {    \
> +        test_c(float, func, z, result, exceptmask, excepts, checksign); \
> +} while (0)
> +
> +/* Test within a given tolerance. */
> +#define    test_tol(func, z, result, tol)                do {    \
> +    volatile long double complex _d = z;                \
> +    assert(cfpequal_tol((func)(_d), (result), (tol),        \
> +        FPE_ABS_ZERO | CS_BOTH));                    \
> +} while (0)
> +
> +/* Test all the functions that compute cexp(x). */
> +#define    testall(x, result, exceptmask, excepts, checksign)    do {    \
> +    test(cexp, x, result, exceptmask, excepts, checksign);        \
> +    test_f(cexpf, x, result, exceptmask, excepts, checksign);        \
> +} while (0)
> +
> +/*
> + * Test all the functions that compute cexp(x), within a given tolerance.
> + * The tolerance is specified in ulps.
> + */
> +#define    testall_tol(x, result, tol)                do {    \
> +    test_tol(cexp, x, result, tol * DBL_ULP());            \
> +    test_tol(cexpf, x, result, tol * FLT_ULP());            \
> +} while (0)
> +
> +/* Various finite non-zero numbers to test. */
> +static const float finites[] =
> +{ -42.0e20, -1.0, -1.0e-10, -0.0, 0.0, 1.0e-10, 1.0, 42.0e20 };
> +
> +
> +/* Tests for 0 */
> +static void
> +test_zero(void)
> +{
> +
> +    /* cexp(0) = 1, no exceptions raised */
> +    testall(0.0, 1.0, ALL_STD_EXCEPT, 0, 1);
> +    testall(-0.0, 1.0, ALL_STD_EXCEPT, 0, 1);
> +    testall(CMPLXL(0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1);
> +    testall(CMPLXL(-0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1);
> +}
> +
> +/*
> + * Tests for NaN.  The signs of the results are indeterminate unless the
> + * imaginary part is 0.
> + */
> +static void
> +test_nan(void)
> +{
> +    unsigned i;
> +
> +    /* cexp(x + NaNi) = NaN + NaNi and optionally raises invalid */
> +    /* cexp(NaN + yi) = NaN + NaNi and optionally raises invalid (|y|>0) */
> +    for (i = 0; i < nitems(finites); i++) {
> +        printf("# Run %d..\n", i);
> +        testall(CMPLXL(finites[i], NAN), CMPLXL(NAN, NAN),
> +            ALL_STD_EXCEPT & ~FE_INVALID, 0, 0);
> +        if (finites[i] == 0.0)
> +            continue;
> +        /* XXX FE_INEXACT shouldn't be raised here */
> +        testall(CMPLXL(NAN, finites[i]), CMPLXL(NAN, NAN),
> +            ALL_STD_EXCEPT & ~(FE_INVALID | FE_INEXACT), 0, 0);
> +    }
> +
> +    /* cexp(NaN +- 0i) = NaN +- 0i */
> +    testall(CMPLXL(NAN, 0.0), CMPLXL(NAN, 0.0), ALL_STD_EXCEPT, 0, 1);
> +    testall(CMPLXL(NAN, -0.0), CMPLXL(NAN, -0.0), ALL_STD_EXCEPT, 0, 1);
> +
> +    /* cexp(inf + NaN i) = inf + nan i */
> +    testall(CMPLXL(INFINITY, NAN), CMPLXL(INFINITY, NAN),
> +        ALL_STD_EXCEPT, 0, 0);
> +    /* cexp(-inf + NaN i) = 0 */
> +    testall(CMPLXL(-INFINITY, NAN), CMPLXL(0.0, 0.0),
> +        ALL_STD_EXCEPT, 0, 0);
> +    /* cexp(NaN + NaN i) = NaN + NaN i */
> +    testall(CMPLXL(NAN, NAN), CMPLXL(NAN, NAN),
> +        ALL_STD_EXCEPT, 0, 0);
> +}
> +
> +static void
> +test_inf(void)
> +{
> +    unsigned i;
> +
> +    /* cexp(x + inf i) = NaN + NaNi and raises invalid */
> +    for (i = 0; i < nitems(finites); i++) {
> +        printf("# Run %d..\n", i);
> +        testall(CMPLXL(finites[i], INFINITY), CMPLXL(NAN, NAN),
> +            ALL_STD_EXCEPT, FE_INVALID, 1);
> +    }
> +    /* cexp(-inf + yi) = 0 * (cos(y) + sin(y)i) */
> +    /* XXX shouldn't raise an inexact exception */
> +    testall(CMPLXL(-INFINITY, M_PI_4), CMPLXL(0.0, 0.0),
> +        ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +    testall(CMPLXL(-INFINITY, 3 * M_PI_4), CMPLXL(-0.0, 0.0),
> +        ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +    testall(CMPLXL(-INFINITY, 5 * M_PI_4), CMPLXL(-0.0, -0.0),
> +        ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +    testall(CMPLXL(-INFINITY, 7 * M_PI_4), CMPLXL(0.0, -0.0),
> +        ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +    testall(CMPLXL(-INFINITY, 0.0), CMPLXL(0.0, 0.0),
> +        ALL_STD_EXCEPT, 0, 1);
> +    testall(CMPLXL(-INFINITY, -0.0), CMPLXL(0.0, -0.0),
> +        ALL_STD_EXCEPT, 0, 1);
> +    /* cexp(inf + yi) = inf * (cos(y) + sin(y)i) (except y=0) */
> +    /* XXX shouldn't raise an inexact exception */
> +    testall(CMPLXL(INFINITY, M_PI_4), CMPLXL(INFINITY, INFINITY),
> +        ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +    testall(CMPLXL(INFINITY, 3 * M_PI_4), CMPLXL(-INFINITY, INFINITY),
> +        ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +    testall(CMPLXL(INFINITY, 5 * M_PI_4), CMPLXL(-INFINITY, -INFINITY),
> +        ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +    testall(CMPLXL(INFINITY, 7 * M_PI_4), CMPLXL(INFINITY, -INFINITY),
> +        ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +    /* cexp(inf + 0i) = inf + 0i */
> +    testall(CMPLXL(INFINITY, 0.0), CMPLXL(INFINITY, 0.0),
> +        ALL_STD_EXCEPT, 0, 1);
> +    testall(CMPLXL(INFINITY, -0.0), CMPLXL(INFINITY, -0.0),
> +        ALL_STD_EXCEPT, 0, 1);
> +}
> +
> +static void
> +test_reals(void)
> +{
> +    unsigned i;
> +
> +    for (i = 0; i < nitems(finites); i++) {
> +        /* XXX could check exceptions more meticulously */
> +        printf("# Run %d..\n", i);
> +        test(cexp, CMPLXL(finites[i], 0.0),
> +             CMPLXL(exp(finites[i]), 0.0),
> +             FE_INVALID | FE_DIVBYZERO, 0, 1);
> +        test(cexp, CMPLXL(finites[i], -0.0),
> +             CMPLXL(exp(finites[i]), -0.0),
> +             FE_INVALID | FE_DIVBYZERO, 0, 1);
> +        test_f(cexpf, CMPLXL(finites[i], 0.0),
> +             CMPLXL(expf(finites[i]), 0.0),
> +             FE_INVALID | FE_DIVBYZERO, 0, 1);
> +        test_f(cexpf, CMPLXL(finites[i], -0.0),
> +             CMPLXL(expf(finites[i]), -0.0),
> +             FE_INVALID | FE_DIVBYZERO, 0, 1);
> +    }
> +}
> +
> +static void
> +test_imaginaries(void)
> +{
> +    unsigned i;
> +
> +    for (i = 0; i < nitems(finites); i++) {
> +        printf("# Run %d..\n", i);
> +        test(cexp, CMPLXL(0.0, finites[i]),
> +             CMPLXL(cos(finites[i]), sin(finites[i])),
> +             ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +        test(cexp, CMPLXL(-0.0, finites[i]),
> +             CMPLXL(cos(finites[i]), sin(finites[i])),
> +             ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +        test_f(cexpf, CMPLXL(0.0, finites[i]),
> +             CMPLXL(cosf(finites[i]), sinf(finites[i])),
> +             ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +        test_f(cexpf, CMPLXL(-0.0, finites[i]),
> +             CMPLXL(cosf(finites[i]), sinf(finites[i])),
> +             ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);
> +    }
> +}
> +
> +static void
> +test_small(void)
> +{
> +    static const double tests[] = {
> +         /* csqrt(a + bI) = x + yI */
> +         /* a    b    x            y */
> +         1.0,    M_PI_4,    M_SQRT2 * 0.5 * M_E,    M_SQRT2 * 0.5 * M_E,
> +        -1.0,    M_PI_4,    M_SQRT2 * 0.5 / M_E,    M_SQRT2 * 0.5 / M_E,
> +         2.0,    M_PI_2,    0.0,            M_E * M_E,
> +         M_LN2,    M_PI,    -2.0,            0.0,
> +    };
> +    double a, b;
> +    double x, y;
> +    unsigned i;
> +
> +    for (i = 0; i < nitems(tests); i += 4) {
> +        printf("# Run %d..\n", i);
> +        a = tests[i];
> +        b = tests[i + 1];
> +        x = tests[i + 2];
> +        y = tests[i + 3];
> +        test_tol(cexp, CMPLXL(a, b), CMPLXL(x, y), 3 * DBL_ULP());
> +
> +        /* float doesn't have enough precision to pass these tests */
> +        if (x == 0 || y == 0)
> +            continue;
> +        test_tol(cexpf, CMPLXL(a, b), CMPLXL(x, y), 1 * FLT_ULP());
> +        }
> +}
> +
> +/* Test inputs with a real part r that would overflow exp(r). */
> +static void
> +test_large(void)
> +{
> +
> +    test_tol(cexp, CMPLXL(709.79, 0x1p-1074),
> +         CMPLXL(INFINITY, 8.94674309915433533273e-16), DBL_ULP());
> +    test_tol(cexp, CMPLXL(1000, 0x1p-1074),
> +         CMPLXL(INFINITY, 9.73344457300016401328e+110), DBL_ULP());
> +    test_tol(cexp, CMPLXL(1400, 0x1p-1074),
> +         CMPLXL(INFINITY, 5.08228858149196559681e+284), DBL_ULP());
> +    test_tol(cexp, CMPLXL(900, 0x1.23456789abcdep-1020),
> +         CMPLXL(INFINITY, 7.42156649354218408074e+83), DBL_ULP());
> +    test_tol(cexp, CMPLXL(1300, 0x1.23456789abcdep-1020),
> +         CMPLXL(INFINITY, 3.87514844965996756704e+257), DBL_ULP());
> +
> +    test_tol(cexpf, CMPLXL(88.73, 0x1p-149),
> +         CMPLXL(INFINITY, 4.80265603e-07), 2 * FLT_ULP());
> +    test_tol(cexpf, CMPLXL(90, 0x1p-149),
> +         CMPLXL(INFINITY, 1.7101492622e-06f), 2 * FLT_ULP());
> +    test_tol(cexpf, CMPLXL(192, 0x1p-149),
> +         CMPLXL(INFINITY, 3.396809344e+38f), 2 * FLT_ULP());
> +    test_tol(cexpf, CMPLXL(120, 0x1.234568p-120),
> +         CMPLXL(INFINITY, 1.1163382522e+16f), 2 * FLT_ULP());
> +    test_tol(cexpf, CMPLXL(170, 0x1.234568p-120),
> +         CMPLXL(INFINITY, 5.7878851079e+37f), 2 * FLT_ULP());
> +}
> +
> +int
> +main(void)
> +{
> +
> +    printf("1..7\n");
> +
> +    test_zero();
> +    printf("ok 1 - cexp zero\n");
> +
> +    test_nan();
> +    printf("ok 2 - cexp nan\n");
> +
> +    test_inf();
> +    printf("ok 3 - cexp inf\n");
> +
> +    test_reals();
> +    printf("ok 4 - cexp reals\n");
> +
> +    test_imaginaries();
> +    printf("ok 5 - cexp imaginaries\n");
> +
> +    test_small();
> +    printf("ok 6 - cexp small\n");
> +
> +    test_large();
> +    printf("ok 7 - cexp large\n");
> +
> +    return (0);
> +}
>
Reply | Threaded
Open this post in threaded view
|

Re: problems with libm

Theo de Raadt-2
Daniel Dickman <[hidden email]> wrote:


> > msun library? Ignore the problem? Patch the current implementation?
>
> Patches to correct the current code would be very welcome. (In my case especially if it improves numpy test results).

Agree.  Don't replace the code.

Replacements are likely to contain other bugs, in if pulled from a repo
that supports fewer plaforms.

So it's same process as other code:  find bugs, fix bugs, repeat...

Reply | Threaded
Open this post in threaded view
|

Re: problems with libm

Ingo Feinerer-2
In reply to this post by Moritz Buhl-3
Moritz Buhl wrote:
> ... I noticed that some floating point operations cause failures of other tests.
> ...
> Many edge cases for complex floating point operations are not covered at all.

Hi,

https://marc.info/?l=openbsd-tech&m=150737856618497&w=2 is another example of
an edge case for complex floating point operations.

https://github.com/wch/r-source/blob/trunk/src/main/complex.c#L452-L455 gives
a solution by checking if the imaginary part of the input complex number is
too large (as otherwise sinh() is called which grows exponentially (see e.g.
https://www.wolframalpha.com/input/?i=sinh(x) ) resulting in an overflow.)

Note that the ctan() implementation in R is under GPL, so I am unsure if the
check can be taken as is and committed to OpenBSD.

s_ctanf.c probably needs a similar treatment.

Best regards,
Ingo

Index: s_ctan.c
===================================================================
RCS file: /cvs/src/lib/libm/src/s_ctan.c,v
retrieving revision 1.7
diff -u -p -r1.7 s_ctan.c
--- s_ctan.c 12 Sep 2016 19:47:02 -0000 1.7
+++ s_ctan.c 11 Jul 2019 12:31:41 -0000
@@ -135,9 +135,11 @@ double complex
 ctan(double complex z)
 {
  double complex w;
- double d;
+ double d, wy, x, y;
 
- d = cos (2.0 * creal (z)) + cosh (2.0 * cimag (z));
+ x = 2.0 * creal(z);
+ y = 2.0 * cimag(z);
+ d = cos(x) + cosh(y);
 
  if (fabs(d) < 0.25)
  d = _ctans (z);
@@ -148,7 +150,12 @@ ctan(double complex z)
  return (w);
  }
 
- w = sin (2.0 * creal(z)) / d + (sinh (2.0 * cimag(z)) / d) * I;
+ if (isnan(y) || fabs(y) < 50.0)
+ wy = sinh(y) / d;
+ else
+ wy = (y < 0 ? -1.0 : 1.0);
+
+ w = sin(x) / d + wy * I;
  return (w);
 }
 DEF_STD(ctan);

Reply | Threaded
Open this post in threaded view
|

Re: problems with libm

Moritz Buhl-3
In reply to this post by Moritz Buhl-3
Hi,

I made the FreeBSD msun regression tests compile on OpenBSD.
https://github.com/moritzbuhl/msun-regress

3 out of 19 test files pass. 14 files die after the first error case.
Two files (ctrig_test.c and trig_test.c) use atf and after some hacks
they report all error cases. 840 for ctrig_test and 88 for trig_test.

These test files should be reviewed carefully as I know for sure that many
don't work on i386 (adding some volatile keywords usually helps).

I believe all these errors paint a good picture. I will be looking into
fixing what I can.