Re: problems with libm

[Moritz Buhl]

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.

Re: problems with libm

[Ingo Feinerer]

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=150737856618497=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.c12 Sep 2016 19:47:02 -  1.7
+++ s_ctan.c11 Jul 2019 12:31:41 -
@@ -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);

Re: problems with libm

[Theo de Raadt]

Daniel Dickman  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...

Re: problems with libm

[Daniel Dickman]

> On Jul 1, 2019, at 2:50 AM, Moritz Buhl  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=140167886413618=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/Makefile21 Feb 2019 16:14:03 -1.1.1.1
> +++ regress/lib/libm/msun/Makefile31 May 2019 19:50:32 -
> @@ -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/null1 Jan 1970 00:00:00 -
> +++ regress/lib/libm/msun/cexp_test.c1 Jun 2019 05:40:51 -
> @@ -0,0 +1,326 @@
> +/*-
> + * Copyright (c) 2008-2011 David Schultz 
> + * 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 
> +
> +#include 
> +
> +#include 
> +#include 
> +#include 
> +#include 
> +#include 
> +#include 
> +
> +#include "test-utils.h"
> +
> +#pragma STDC FENV_ACCESSON
> +#pragmaSTDC CX_LIMITED_RANGEOFF
> +
> +/*
> + * 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
> + *