by vogt@linux.vnet.ibm.com:

This is a copy of the bug report from the gcc bugtracker.  See
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63269 for the mentioned patch.

--


Dominik Vogt 2014-09-15 14:10:37 UTC
------------------------------------
There are several flaws in the math library test code.  On s390x, the test TestLog2
fails for three reasons:

1) This part of the test is too strict for i == 0:

for i := -1074; i <= 1023; i++ { 
  f := Ldexp(1, i) 
  l := Log2(f) 
  if l != float64(i) { 
    t.Errorf("Log2(2**%d) = %g, want %d", i, l, i) 
  } 
} 

It should really use the veryclose function() as there is no guarantee that the result
of this calculation from log10.go is exact:

  return Log(frac)*(1/Ln2) + float64(exp)

With frac == 0.5 and exp == 1, the result may be just close to 0, not exactly zero:

* If the compiler generated two instructions, one to mulitply and one to add, the result
is (with a small error delta):

  Product: Log(0.5)*(1/Ln2) = -1 + delta -> rounded to -1
  Sum: -1 + 1 = 0 -> rounded to 0
  Result: 0

* If the compiler generates a single multiply and add instruction, only the final result
is rouced, i.e.

  Product: Log(0.5)*(1/Ln2) = -1 + delta (not rounded)
  Sum: -1 + delta + 1 = delta -> rounded to delta'
  Result: delta' != 0

So, the fix is to use the "veryclose" error tolerance.

2) There's a bug in the tolerance() function:

        if a != 0 { 
                e = e * a 
                if e < 0 { 
                        e = -e 
                } 
        } 

This snippet uses the defective value a but should use the expected value b instead:

        if b != 0 { 
                e = e * b 
                if e < 0 { 
                        e = -e 
                } 
        } 

Otherwise, bad things happen for an expected value 0 when the defective value is very
close, i.e. they never match.

3) The alike() function does allow an error tolerance.  This kicks in on s390x in the
case that the expected value is 0 and the defective value is not identical.

        switch { 
        case IsNaN(a) && IsNaN(b): 
                return true 
        case a == b: 
                return Signbit(a) == Signbit(b) 
        } 


A possible fix would be to detect this one situation and return true, and rely on more
specific tests to chack that the error is small ebough.

        switch { 
        case IsNaN(a) && IsNaN(b): 
                return true 
        case a == 0 && !IsNaN(b) && !IsInf(b, 0): 
                // allow deviations when the expected value is zero
                return true 
        case a == b: 
                return Signbit(a) == Signbit(b) 
        }


Comment 3 Ian Lance Taylor 2014-11-05 03:54:14 UTC
--------------------------------------------------
First, let me say that this code is in the Go master library and must be fixed there. 
It might be more effective to discuss it on the Go issue tracker at
http://golang.org/issue.

I don't agree with your argument for item 1.  You say that the wrong result happens when
using a fused multiply-add instruction.  The math library is compiled with
-ffp-contract=off (see MATH_FLAG in configure.ac and Makefile.am), so the compiler
should not be generating a fused multiply-add instruction.

I'm not entirely persuaded by your argument for item 2.  Zero is a special value.  When
we expect a zero, we should get a zero, not something close to zero.  I don't think this
change is correct in general.  It may be correct for some specific cases, but then we
need to investigate those.

Item 3 is the same sort of thing: when we expect zero, we should, in general, get
exactly zero.


Comment 4 Dominik Vogt 2014-11-05 14:34:38 UTC
----------------------------------------------
regarding 2)

> I'm not entirely persuaded by your argument for item 2. ...

Hm, good that you doubted it, because the actual mistake is somehwere else:  The
unpatched code has

  if l != float64(i)

but if you want to use a tolerance here this must become

  if !veryclose(float64(i), l) {

With the argument reversed.  This could/should be cleaned up by renaming the arguments
of the tolerance() function, e.g. a -> expected, b -> result, e -> maxerr.

> Zero is a special
> value.  When we expect a zero, we should get a zero, not something close to
> zero.  I don't think this change is correct in general.  It may be correct for
> some specific cases, but then we need to investigate those.

Actually, this has nothing to do with 0 being special here, abut with scaling of the
allowed error: Multiplying it by 0 yields zero error tolerance, so the tolerance()
function does not do that.

--> This chunk is not necessary, but a (separate) cleanup patch might help to avoid
future confusion.


Comment 5 Dominik Vogt 2014-11-05 16:47:03 UTC
----------------------------------------------

regarding 1)

My earlier explanation of the problem was wrong.  Multiply and add is not generated; it
probably only was in the artificial test case that I made and certainly did not compile
with -ffp-contract=off.

In this calculation in log2(),

  Log(frac)*(1/Ln2) + float64(exp)

Gcc does constant folding for (1/Ln2) and generates a multiply instruction and then adds
the second term.  Same result if you write "*Log2E" instead of
"*(1/Ln2)").  But with

  Log(frac)/Ln2 + float64(exp)

it generates a divide instruction.  The multiplication and the division yield results
that differ in the least significant bit, and I don't see how this could be prevented in
general; it's just an artifact of the floating point format.  I've verified that the
constants Ln2, 1/Ln2 and Log2E are bit correct.

The "easy" way to fix this is increasing the allowed tolerance as my patch
does (note that the arguments of the veryclose() call need to be swapped, see previous
comment).

The "right" way to fix this is to calculate platform specific ULPs for all the
algorithms from the math library and use these.  That's what glibc does.