Bug 3976 - libm rounding modes do not work correctly for many archs
Summary: libm rounding modes do not work correctly for many archs
Status: RESOLVED FIXED
Alias: None
Product: glibc
Classification: Unclassified
Component: math (show other bugs)
Version: unspecified
: P2 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords:
: 6869 15010 (view as bug list)
Depends on:
Blocks:
 
Reported: 2007-02-06 13:14 IST by Pierre Habouzit
Modified: 2016-02-26 15:07 IST (History)
6 users (show)

See Also:
Host:
Target:
Build:
Last reconfirmed:
fweimer: security-


Attachments
Test a math function in the 4 rounding modes. (834 bytes, text/x-csrc)
2012-02-29 02:35 IST, Vincent Lefèvre
Details

Note You need to log in before you can comment on or make changes to this bug.
Description Pierre Habouzit 2007-02-06 13:14:10 IST
It seems that changing the rounding modes make some functions (like exp, or 
sin, or ...) buggy on many archs. The testing code is:

===============================================================================
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fenv.h>

static int rnd[4] = { FE_TONEAREST, FE_TOWARDZERO, FE_DOWNWARD, FE_UPWARD };
static char rc[4] = "NZDU";

int main (int argc, char *argv[])
{
  int i;

  for (i = 1; i < argc; i++)
    {
      int r;
      double x;
      char *end;

      x = strtod (argv[i], &end);
      if (*end != '\0')
        exit (EXIT_FAILURE);

      for (r = 0; r < 4; r++)
        {
          double y;

          fesetround (rnd[r]);
          y = exp (x);
          printf ("%c: exp(%.17g) = %.17g\n", rc[r], x, y);
        }
    }

  return 0;
}
===============================================================================

I've checked on my amd64 that it indeed works in 32bits mode (and it also work 
on an i386 machine) but it does not in 64bits mode:

[madcoder mad] gcc -m32 -lm -o a a.c ; ./a 1 2
N: exp(1) = 2.7182818284590451
Z: exp(1) = 2.7182818284590451
D: exp(1) = 2.7182818284590451
U: exp(1) = 2.7182818284590455
N: exp(2) = 7.3890560989306504
Z: exp(2) = 7.3890560989306495
D: exp(2) = 7.3890560989306495
U: exp(2) = 7.3890560989306504

[madcoder mad] gcc -m64 -lm -o a a.c ; ./a 1 2
N: exp(1) = 2.7182818284590451
Z: exp(1) = 2.7182818284590451
D: exp(1) = 2.7182818284590451
U: exp(1) = 0.04788398250919005
N: exp(2) = 7.3890560989306504
Z: exp(2) = 4.0037745305985499
D: exp(2) = 4.0037745305985499
U: exp(2) = 7.3890560989306504

I tested it on other archs, here is the summary:
 - i386, m68k, ia64 have been tested OK.
 - arm: only FE_ROUNDTONEAREST exists (which is ok as per C standard)
 - amd64, mips, ppc, hppa, s390, sparc produce wrong results, in their own 
unique way.

I've not been able to test alpha yet though.
Comment 1 Pierre Habouzit 2007-02-06 22:13:19 IST
(In reply to comment #0)
> I've not been able to test alpha yet though.

alpha gives curious results as all values are exactly the same, which is quite 
reasonnable, but may hide a bug too.

It gave (for 1 2 3 4 5):

N: exp(1) = 2.7182818284590451
Z: exp(1) = 2.7182818284590451
D: exp(1) = 2.7182818284590451
U: exp(1) = 2.7182818284590451
N: exp(2) = 7.3890560989306504
Z: exp(2) = 7.3890560989306504
D: exp(2) = 7.3890560989306504
U: exp(2) = 7.3890560989306504
N: exp(3) = 20.085536923187668
Z: exp(3) = 20.085536923187668
D: exp(3) = 20.085536923187668
U: exp(3) = 20.085536923187668
N: exp(4) = 54.598150033144236
Z: exp(4) = 54.598150033144236
D: exp(4) = 54.598150033144236
U: exp(4) = 54.598150033144236
N: exp(5) = 148.4131591025766
Z: exp(5) = 148.4131591025766
D: exp(5) = 148.4131591025766
U: exp(5) = 148.4131591025766
Comment 2 Vincent Lefèvre 2008-02-25 13:11:22 IST
You can see other results in the bug I originally reported on the Debian BTS:
  http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=153022
In short, sin can even give out-of-range results, and pow can segfault.
Comment 3 Christian Keil 2008-10-27 08:17:04 IST
At least on my Core2 Duo it's also not working with still other values (32bit
version is ok):

gcc -m64 -lm -o a a.c; ./a 1 2
N: exp(1) = 2.7182818284590451
Z: exp(1) = 2.7182818284590451
D: exp(1) = 2.7182818284590451
U: exp(1) = 7.1387612927397726
N: exp(2) = 7.3890560989306504
Z: exp(2) = 4.0037745305985499
D: exp(2) = 4.0037745305985499
U: exp(2) = 7.3890560989306504

As far as I dug into it, the 32bit and 64bit versions use other code. 64bit
comes from the IBM Accurate Mathematical Library. For exp
sysdeps/ieee754/dbl-64/e_exp.c produces the wrong results.
Comment 4 Jakub Jelinek 2008-10-27 08:25:06 IST
The functions have defined behavior only in the default rounding mode
(round-to-even), anything else is undefined behavior and completely programmer's
fault for calling the functions in those rounding modes.
Comment 5 Vincent Lefèvre 2008-10-27 09:33:27 IST
(In reply to comment #4)
> The functions have defined behavior only in the default rounding mode
> (round-to-even), anything else is undefined behavior and completely programmer's
> fault for calling the functions in those rounding modes.

No, it isn't. There's no undefined behavior there. The C standard says (F.9):
"Whether the functions honor the rounding direction mode is implementation-defined."

So, this just means that an implementation doesn't need to return a result
rounded to the correct direction (this is implementation-defined). Still it
should return an approximate value whatever the rounding direction mode is, and
not behave erratically.
Comment 6 joseph@codesourcery.com 2008-10-27 12:42:01 IST
Subject: Re:  libm rounding modes do not work correctly for
 many archs

On Mon, 27 Oct 2008, vincent+libc at vinc17 dot org wrote:

> So, this just means that an implementation doesn't need to return a result
> rounded to the correct direction (this is implementation-defined). Still it
> should return an approximate value whatever the rounding direction mode is, and
> not behave erratically.

Furthermore, some functions whose implementations only work in 
round-to-nearest mode do save the mode then set round-to-nearest using 
feholdexcept and fesetround (e.g. sysdeps/ieee754/dbl-64/e_exp2.c).  I 
think this is the best thing to do in such cases.  (Regarding the issue of 
rounding modes in signal handlers, I think the best thing is for ABI 
documents to make explicit that library functions may temporarily save and 
restore rounding modes; that's what is being done in the Power 
Architecture ABI document being worked on, to document what the de facto 
ABI is right now.)

Comment 7 Vincent Lefèvre 2010-03-21 22:24:35 IST
*** Bug 6869 has been marked as a duplicate of this bug. ***
Comment 8 Joseph Myers 2012-02-29 01:14:47 IST
Problem with exp confirmed on x86_64 with current sources.

Patch for exp proposed: http://sourceware.org/ml/libc-alpha/2012-02/msg00748.html

Given that patch applied:

* I cannot confirm the problem with sin or cos on x86_64 (though tests should be added to the testsuite).

* pow (1.6, 1.6) does not segfault, but the result in round-upward mode is substantially inaccurate; pow will need a similar fix (and test in the testsuite).

If other functions have problems in current sources, it would be best to have a separate bug for each function with problems (identifying clearly the platform on which the problem occurs).
Comment 9 Vincent Lefèvre 2012-02-29 02:35:24 IST
Created attachment 6256 [details]
Test a math function in the 4 rounding modes.

(In reply to comment #8)
> * I cannot confirm the problem with sin or cos on x86_64 (though tests should
> be added to the testsuite).

I still get the bug on the argument 100 under Debian (glibc 2.13).

> * pow (1.6, 1.6) does not segfault, but the result in round-upward mode is
> substantially inaccurate;

I confirm, but pow(1.01,1.1) crashes:

N: pow(1.01,1.1000000000000001) = 1.0110054835779234
Z: pow(1.01,1.1000000000000001) = 1.0110054835779232
D: pow(1.01,1.1000000000000001) = 1.0110054835779232
zsh: segmentation fault (core dumped)  ./tfct-4rm 1.01 1.1

> pow will need a similar fix (and test in the testsuite).

Yes, like the other functions.

> If other functions have problems in current sources, [...]

I would say that each function probably has the same problem.

I did the tests with

gcc -std=c99 tfct-4rm.c -o tfct-4rm -lm -DFCT=exp
gcc -std=c99 tfct-4rm.c -o tfct-4rm -lm -DFCT=sin
gcc -std=c99 tfct-4rm.c -o tfct-4rm -lm -DFCT=cos
gcc -std=c99 tfct-4rm.c -o tfct-4rm -lm -DFCT=pow -DTWOARGS

using the attached code.
Comment 10 Vincent Lefèvre 2012-02-29 02:54:10 IST
(In reply to comment #9)
> I would say that each function probably has the same problem.

Actually log and atan seem to behave correctly (and even honor the rounding mode).

But tan, sinh and cosh are buggy:

$ ./tfct-4rm 1.6
N: tan(1.6000000000000001) = -34.232532735557314
Z: tan(1.6000000000000001) = -34.232532735557307
D: tan(1.6000000000000001) = -34.232532735557307
U: tan(1.6000000000000001) = -1.9000553624671392

$ ./tfct-4rm 100
N: sinh(100) = 1.3440585709080678e+43
Z: sinh(100) = 1.3440585709080676e+43
D: sinh(100) = 1.3440585709080676e+43
U: sinh(100) = -5.7576991416613074e+42

$ ./tfct-4rm 100
N: cosh(100) = 1.3440585709080678e+43
Z: cosh(100) = 1.3440585709080676e+43
D: cosh(100) = 1.3440585709080676e+43
U: cosh(100) = -5.7576991416613074e+42

tanh seems accurate, but it doesn't honor the rounding mode:

$ ./tfct-4rm 0.8
N: tanh(0.80000000000000004) = 0.66403677026784913
Z: tanh(0.80000000000000004) = 0.6640367702678488
D: tanh(0.80000000000000004) = 0.66403677026784902
U: tanh(0.80000000000000004) = 0.66403677026784891
Comment 11 Joseph Myers 2012-03-02 15:16:39 IST
exp fixed by:

commit 28afd92dbdb4fef4358051aad5cb944a9527a4b5
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Fri Mar 2 15:12:53 2012 +0000

    Fix exp in non-default rounding modes (bug 3976).

This may fix some cases of sinh and cosh (where they use exp internally) though tests still need adding for those functions and other functions still need to be fixed.
Comment 12 Joseph Myers 2012-03-02 16:04:12 IST
sinh and cosh appear to have been fixed by the exp fix.  I confirm what you report for tan and pow.  For sin and cos, I don't see the problem for 100 with current sources, but do for 7:

N: cos(7) = 0.7539022543433046
Z: cos(7) = 0.7539022543433046
D: cos(7) = 0.7539022543433046
U: cos(7) = -3.5593280928702544e+244
N: sin(7) = 0.65698659871878906
Z: sin(7) = 0.65698659871878906
D: sin(7) = 0.65698659871878906
U: sin(7) = 6.5993918533387462e+246
Comment 13 Joseph Myers 2012-03-02 16:08:45 IST
Examples for pow involving arguments that are exactly representable in all floating-point types (always nice for the testsuite):

N: pow(1.0625,1.125) = 1.0705822930287614
Z: pow(1.0625,1.125) = 1.0705822930287612
D: pow(1.0625,1.125) = 1.0705822930287612
U: pow(1.0625,1.125) = 0.032633984947502387
N: pow(1.5,1.03125) = 1.5191270987147432
Z: pow(1.5,1.03125) = 1.0003045181943615
D: pow(1.5,1.03125) = 1.0003045181943615
U: pow(1.5,1.03125) = 1.5191270987147434
Comment 14 Joseph Myers 2012-03-02 20:54:41 IST
sin, cos, tan fixed by:

commit 804360ed837e3347c9cd9738f25345d2587a1242
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Fri Mar 2 20:51:39 2012 +0000

    Fix sin, cos, tan in non-default rounding modes (bug 3976).
Comment 15 Joseph Myers 2012-03-05 12:27:36 IST
cosh, sinh tests (some errors in round-upward mode up to 27ulp) added by:

commit ca811b2256d2e48c7288219e9e11dcbab3000f19
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Mon Mar 5 12:20:24 2012 +0000

    Test cosh, sinh in non-default rounding modes (bug 3976).

pow fixed by:

commit b7cd39e8f8c5cf2844f20eb03f545d19c4c25987
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Mon Mar 5 12:22:46 2012 +0000

    Fix pow in non-default rounding modes (bug 3976).

So all the reported issues are fixed here.  If cases are found of functions that still (in any rounding mode, for any of float, double, long double) produce wild results or crash with glibc after my patches - or of functions that should always produce correctly rounded results such as rint, nextafter, fma, sqrt, that fail to do so (in any rounding mode, for any of float, double, long double), please file them as separate bugs for each function found to have a problem.
Comment 16 Vincent Lefèvre 2012-03-05 13:32:24 IST
(In reply to comment #10)
> (In reply to comment #9)
> > I would say that each function probably has the same problem.
> 
> Actually log and atan seem to behave correctly (and even honor the rounding
> mode).

I forgot that since Ziv's strategy is used to provide correct rounding, my tests where incomplete (I tested only on a few random values). To understand the problem, here's what Ziv's strategy does:

1. Compute an approximation with slightly more precision than the target precision, typically by representing it as an exact sum of two floating-point values a+b, where |b| is much smaller than |a|. An error bound e1 is determined and a rounding test is performed. Most of the time, the correct rounding can be guaranteed by this test, so that the computed result can be returned. But there's a low probability of failure (problem known as the "Table Maker's Dilemma"), in which case a second step is necessary.

2. Ditto with even more precision and a lower error bound e2. In case of failure, a third step is necessary.

3. Computation with even more precision. IIRC, a 768-bit precision is used here in IBM's library (that's the mp*.c files), which assumes that this is sufficient.

So, for functions that are believed to behave correctly in the directed rounding modes, and for which the active rounding mode is not changed temporarily in the function, one should consider the following remarks:

One would need to test cases that fall in (2) and in (3). The probabilities of failure depend on the algorithm and implied error bounds e1 and e2, and possibly on the magnitude of the input (or some other criteria). To get cases in (2) from random tests, one may need to test several dozens of thousands of inputs (analyzing the code can give a better idea of what is needed -- it doesn't seem to be much documented). Ideally, tools like gcov to get code coverage would be useful. From my work on the Table Maker's Dilemma, I also have a list of the hardest-to-round cases ("worst cases") for some functions in some domains, which should probably fall in (3).

AFAIK, in IBM's library, (3) uses integer arithmetic, so that it shouldn't depend on the active rounding mode (though I'm not completely sure because FP may also be used, e.g. at the end).

Note that testing cases for (3) only is not sufficient, because with such cases, the results from (2) are discarded.

It is also possible that in directed rounding modes, (1) and/or (2) always fail, meaning that a slower step would be performed. This wouldn't be incorrect, but rather inefficient. This should be checked too...
Comment 17 Andreas Jaeger 2013-01-10 19:22:58 IST
*** Bug 15010 has been marked as a duplicate of this bug. ***
Comment 18 Jackie Rosen 2014-02-16 19:37:05 IST Comment hidden (spam)