diff mbox

libgo patch committed: Compile math library with -ffp-contract=off

Message ID 87zjkp4sm5.fsf@canonical.com
State New
Headers show

Commit Message

Michael-Doyle Hudson March 16, 2014, 11:43 p.m. UTC 
Michael Hudson-Doyle <michael.hudson@linaro.org> writes:

> Ian Lance Taylor <iant@google.com> writes:
>
>> On Thu, Mar 13, 2014 at 6:27 PM, Michael Hudson-Doyle
>> <michael.hudson@linaro.org> wrote:
>>> Ian Lance Taylor <iant@google.com> writes:
>>>
>>>> The bug report http://golang.org/issue/7074 shows that math.Log2(1)
>>>> produces the wrong result on Aarch64, because the Go math package is
>>>> compiled to use a fused multiply-add instruction.  This patch to the
>>>> libgo configure script will use -ffp-contract=off when compiling the
>>>> math package on processors other than x86.  Bootstrapped and ran Go
>>>> testsuite on x86_64-unknown-linux-gnu, not that that tests much.
>>>> Committed to mainline.
>>>
>>> Thanks for this!  If you are willing to go into battle enough to argue
>>> that libgcc should also be compiled with -ffp-contract=off (I did not
>>> have the stomach for that fight) then we'll be down to 1 check-go
>>> failure on aarch64 (which is peano -- due to the absence of
>>> split/copyable stacks and should probably xfail).
>>
>> Hmmm, what is it that fails with libgcc?  Is there a bug report for
>> it?
>
>     https://code.google.com/p/go/issues/detail?id=7066
>
> and then
>
>     http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714
>
> I wanted to propose a version using "Kahan's algorithm for the
> determinant" as described in
>
> http://hal-ens-lyon.archives-ouvertes.fr/docs/00/78/57/86/PDF/Jeannerod_Louvet_Muller_final.pdf
>
> but I haven't gotten around to it...

I got bored / distracted and hacked up this:


it would be better to do this in libgcc of course but I think that's
awkward because libgcc can't link to libm and so on...  It's probably a
little slower than the libgcc version (although this is straight line
code) but I don't really care about that :-)

Cheers,
mwh
diff mbox

Patch

diff --git a/libgo/runtime/go-cdiv.c b/libgo/runtime/go-cdiv.c
index 0a81e45..b96576a 100644
--- a/libgo/runtime/go-cdiv.c
+++ b/libgo/runtime/go-cdiv.c
@@ -13,6 +13,8 @@ 
    the the whole number is Inf, but an operation involving NaN ought
    to result in NaN, not Inf.  */
 
+#include <math.h>
+
 __complex float
 __go_complex64_div (__complex float a, __complex float b)
 {
@@ -29,6 +31,48 @@  __go_complex64_div (__complex float a, __complex float b)
   return a / b;
 }
 
+#ifdef FP_FAST_FMA
+
+// This implements what is sometimes called "Kahan's compensated algorithm for
+// 2 by 2 determinants".  It returns a*b + c*d to a high degree of precision
+// but depends on a fused-multiply add operation that rounds once.
+//
+// The obvious algorithm has problems when a*b and c*d nearly cancel, but the
+// trick is the calculation of 'e': "a*b = w + e" is exact when the operands
+// are considered as real numbers.  So if c*d nearly cancels out w, e restores
+// the result to accuracy.
+double
+Kahan(double a, double b, double c, double d)
+{
+  double w, e, f;
+  w = b * a;
+  e = fma(b, a, -w);
+  f = fma(d, c, w);
+  return f + e;
+}
+
+__complex double
+__go_complex128_div (__complex double a, __complex double b)
+{
+  double r, i, denom;
+  if (__builtin_expect (b == 0+0i, 0))
+    {
+      if (!__builtin_isinf (__real__ a)
+	  && !__builtin_isinf (__imag__ a)
+	  && (__builtin_isnan (__real__ a) || __builtin_isnan (__imag__ a)))
+	{
+	  /* Pass "1" to nan to match math/bits.go.  */
+	  return __builtin_nan("1") + __builtin_nan("1")*1i;
+	}
+    }
+  r = Kahan(__real__ a, __real__ b, __imag__ a, __imag__ b);
+  i = Kahan(__imag__ a, __real__ b, - __real__ a, __imag__ b);
+  denom = (__real__ b)*(__real__ b) + (__imag__ b)*(__imag__ b);
+  return r/denom + i*1.0i/denom;
+}
+
+#else
+
 __complex double
 __go_complex128_div (__complex double a, __complex double b)
 {
@@ -44,3 +88,5 @@  __go_complex128_div (__complex double a, __complex double b)
     }
   return a / b;
 }
+
+#endif