Re: Replace cos and avoid FPU trigonometry (was: tanf returns NaN for large inputs)
> From: Greg Steuck
> Date: Mon, 10 Jan 2022 20:59:17 -0800
>
> Greg Steuck writes:
>
> > This failure can be reduced to a trivial program which does change
> > its behavior for the worse if s_cos.S is taken out:
> >
> > #include
> > #include
> >
> > int main(int a, char**b) {
> > double y = -0.34061437849088045332;
> > printf("cos(%lf)=%le delta=%e\n", y, cos(y), 0.94254960031831729956 -
> > cos(y));
> > }
> >
> > In HEAD:
> >
> > cos(-0.340614)=9.425496e-01 delta=-1.110223e-16
> >
> > while with the patch below:
> >
> > cos(-0.340614)=9.425496e-01 delta=0.00e+00
>
> As Daniel noted, I swapped the cases. The HEAD is at 0.0 delta whereas
> the patch used to make it worse.
>
> I went looking for why things are better on FreeBSD and they have a
> different (simpler) implementation of cos. I copied it over. Given the
> common provenance, I expect the copyright situation to be unambiguous.
I think you will also need the changes done in FreeBSD commit
4339c67c485f.
> With the two patches things look almost universally better in
> regress/libm. I attached both logs from amd64.
>
> Anybody has ideas for other tests that make sense to do? Maybe people
> can help me run regress on less common platforms?
>
> Thanks
> Greg
>
> >From a0b065bd3f5d48786f77f654dfb53cbf2617b0b3 Mon Sep 17 00:00:00 2001
> From: Greg Steuck
> Date: Mon, 10 Jan 2022 20:22:07 -0800
> Subject: [PATCH 1/2] Copy cos(3) software implementation from FreeBSD-13
>
> The result passes more tests from msun suite. In particular,
> testacc(cos, -0.34061437849088045332L, 0.94254960031831729956L,
> ALL_STD_EXCEPT, FE_INEXACT);
> matches instead of being 1e-16 off.
> ---
> lib/libm/src/k_cos.c | 45 ++--
> lib/libm/src/s_cos.c | 6 +-
> 2 files changed, 23 insertions(+), 28 deletions(-)
>
> diff --git a/lib/libm/src/k_cos.c b/lib/libm/src/k_cos.c
> index 8f3882b6a00..0839243e90c 100644
> --- a/lib/libm/src/k_cos.c
> +++ b/lib/libm/src/k_cos.c
> @@ -36,13 +36,17 @@
> * ~ cos(x) - x*y,
> * a correction term is necessary in cos(x) and hence
> * cos(x+y) = 1 - (x*x/2 - (r - x*y))
> - * For better accuracy when x > 0.3, let qx = |x|/4 with
> - * the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
> - * Then
> - * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
> - * Note that 1-qx and (x*x/2-qx) is EXACT here, and the
> - * magnitude of the latter is at least a quarter of x*x/2,
> - * thus, reducing the rounding error in the subtraction.
> + * For better accuracy, rearrange to
> + * cos(x+y) ~ w + (tmp + (r-x*y))
> + * where w = 1 - x*x/2 and tmp is a tiny correction term
> + * (1 - x*x/2 == w + tmp exactly in infinite precision).
> + * The exactness of w + tmp in infinite precision depends on w
> + * and tmp having the same precision as x. If they have extra
> + * precision due to compiler bugs, then the extra precision is
> + * only good provided it is retained in all terms of the final
> + * expression for cos(). Retention happens in all cases tested
> + * under FreeBSD, so don't pessimize things by forcibly clipping
> + * any extra precision in w.
> */
>
> #include "math.h"
> @@ -60,25 +64,12 @@ C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9,
> 0xBE8838D4 */
> double
> __kernel_cos(double x, double y)
> {
> - double a,hz,z,r,qx;
> - int32_t ix;
> - GET_HIGH_WORD(ix,x);
> - ix &= 0x7fff; /* ix = |x|'s high word*/
> - if(ix<0x3e40) { /* if x < 2**27 */
> - if(((int)x)==0) return one; /* generate inexact */
> - }
> + double hz,z,r,w;
> +
> z = x*x;
> - r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6);
> - if(ix < 0x3FD3) /* if |x| < 0.3 */
> - return one - (0.5*z - (z*r - x*y));
> - else {
> - if(ix > 0x3fe9) { /* x > 0.78125 */
> - qx = 0.28125;
> - } else {
> - INSERT_WORDS(qx,ix-0x0020,0); /* x/4 */
> - }
> - hz = 0.5*z-qx;
> - a = one-qx;
> - return a - (hz - (z*r-x*y));
> - }
> + w = z*z;
> + r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
> + hz = 0.5*z;
> + w = one-hz;
> + return w + (((one-w)-hz) + (z*r-x*y));
> }
> diff --git a/lib/libm/src/s_cos.c b/lib/libm/src/s_cos.c
> index 8b923d5fe61..1406504e9ab 100644
> --- a/lib/libm/src/s_cos.c
> +++ b/lib/libm/src/s_cos.c
> @@ -57,7 +57,11 @@ cos(double x)
>
> /* |x| ~< pi/4 */
> ix &= 0x7fff;
> - if(ix <= 0x3fe921fb) return __kernel_cos(x,z);
> + if(ix <= 0x3fe921fb) {
> + if(ix<0x3e46a09e) /* if x < 2**-27 * sqrt(2) */
> + if(((int)x)==0) return 1.0; /* generate inexact */
> + return __kernel_cos(x,z);
> + }
Replace cos and avoid FPU trigonometry (was: tanf returns NaN for large inputs)
Greg Steuck writes:
> This failure can be reduced to a trivial program which does change
> its behavior for the worse if s_cos.S is taken out:
>
> #include
> #include
>
> int main(int a, char**b) {
> double y = -0.34061437849088045332;
> printf("cos(%lf)=%le delta=%e\n", y, cos(y), 0.94254960031831729956 -
> cos(y));
> }
>
> In HEAD:
>
> cos(-0.340614)=9.425496e-01 delta=-1.110223e-16
>
> while with the patch below:
>
> cos(-0.340614)=9.425496e-01 delta=0.00e+00
As Daniel noted, I swapped the cases. The HEAD is at 0.0 delta whereas
the patch used to make it worse.
I went looking for why things are better on FreeBSD and they have a
different (simpler) implementation of cos. I copied it over. Given the
common provenance, I expect the copyright situation to be unambiguous.
With the two patches things look almost universally better in
regress/libm. I attached both logs from amd64.
Anybody has ideas for other tests that make sense to do? Maybe people
can help me run regress on less common platforms?
Thanks
Greg
>From a0b065bd3f5d48786f77f654dfb53cbf2617b0b3 Mon Sep 17 00:00:00 2001
From: Greg Steuck
Date: Mon, 10 Jan 2022 20:22:07 -0800
Subject: [PATCH 1/2] Copy cos(3) software implementation from FreeBSD-13
The result passes more tests from msun suite. In particular,
testacc(cos, -0.34061437849088045332L, 0.94254960031831729956L,
ALL_STD_EXCEPT, FE_INEXACT);
matches instead of being 1e-16 off.
---
lib/libm/src/k_cos.c | 45 ++--
lib/libm/src/s_cos.c | 6 +-
2 files changed, 23 insertions(+), 28 deletions(-)
diff --git a/lib/libm/src/k_cos.c b/lib/libm/src/k_cos.c
index 8f3882b6a00..0839243e90c 100644
--- a/lib/libm/src/k_cos.c
+++ b/lib/libm/src/k_cos.c
@@ -36,13 +36,17 @@
* ~ cos(x) - x*y,
*a correction term is necessary in cos(x) and hence
* cos(x+y) = 1 - (x*x/2 - (r - x*y))
- *For better accuracy when x > 0.3, let qx = |x|/4 with
- *the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
- *Then
- * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
- *Note that 1-qx and (x*x/2-qx) is EXACT here, and the
- *magnitude of the latter is at least a quarter of x*x/2,
- *thus, reducing the rounding error in the subtraction.
+ *For better accuracy, rearrange to
+ * cos(x+y) ~ w + (tmp + (r-x*y))
+ *where w = 1 - x*x/2 and tmp is a tiny correction term
+ *(1 - x*x/2 == w + tmp exactly in infinite precision).
+ *The exactness of w + tmp in infinite precision depends on w
+ *and tmp having the same precision as x. If they have extra
+ *precision due to compiler bugs, then the extra precision is
+ *only good provided it is retained in all terms of the final
+ *expression for cos(). Retention happens in all cases tested
+ *under FreeBSD, so don't pessimize things by forcibly clipping
+ *any extra precision in w.
*/
#include "math.h"
@@ -60,25 +64,12 @@ C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9,
0xBE8838D4 */
double
__kernel_cos(double x, double y)
{
- double a,hz,z,r,qx;
- int32_t ix;
- GET_HIGH_WORD(ix,x);
- ix &= 0x7fff; /* ix = |x|'s high word*/
- if(ix<0x3e40) { /* if x < 2**27 */
- if(((int)x)==0) return one; /* generate inexact */
- }
+ double hz,z,r,w;
+
z = x*x;
- r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6);
- if(ix < 0x3FD3) /* if |x| < 0.3 */
- return one - (0.5*z - (z*r - x*y));
- else {
- if(ix > 0x3fe9) { /* x > 0.78125 */
- qx = 0.28125;
- } else {
- INSERT_WORDS(qx,ix-0x0020,0); /* x/4 */
- }
- hz = 0.5*z-qx;
- a = one-qx;
- return a - (hz - (z*r-x*y));
- }
+ w = z*z;
+ r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
+ hz = 0.5*z;
+ w = one-hz;
+ return w + (((one-w)-hz) + (z*r-x*y));
}
diff --git a/lib/libm/src/s_cos.c b/lib/libm/src/s_cos.c
index 8b923d5fe61..1406504e9ab 100644
--- a/lib/libm/src/s_cos.c
+++ b/lib/libm/src/s_cos.c
@@ -57,7 +57,11 @@ cos(double x)
/* |x| ~< pi/4 */
ix &= 0x7fff;
- if(ix <= 0x3fe921fb) return __kernel_cos(x,z);
+ if(ix <= 0x3fe921fb) {
+ if(ix<0x3e46a09e) /* if x < 2**-27 * sqrt(2) */
+ if(((int)x)==0) return 1.0; /* generate inexact */
+ return __kernel_cos(x,z);
+ }
/* cos(Inf or NaN) is NaN */
else if (ix>=0x7ff0) return x-x;
--
2.34.1
>From 7461948f0e7ff41e5b3b9288d98edf13cb296742 Mon Sep 17 00:00:00 2001
From: Greg Steuck
Date: Sun, 9 Jan 2022 13:45:51 -0800
Subject: [PATCH 2/2] Unplug assembly implementations of
Re: tanf returns NaN for large inputs
With numpy/i386 we don't fix any of the currently broken tests and
surprisingly one new regression is introduced:
-13 failed, 10900 passed, 88 skipped, 108 deselected, 19 xfailed, 2
xpassed, 5 warnings in 206.14 seconds
+14 failed, 10899 passed, 88 skipped, 108 deselected, 19 xfailed, 2
xpassed, 5 warnings in 205.49 seconds
The new breakage is due to a test that checks for errors that <= 2
ulps using sin on 32-bit floats vs sin on 64-bit floats (shown below).
Still, it may be worth committing the proposed change and then
separately seeing if we can improve the C versions of these functions.
___ TestAVXFloat32Transcendental.test_sincos_float32 ___
self =
def test_sincos_float32(self):
np.random.seed(42)
N = 100
M = np.int_(N/20)
index = np.random.randint(low=0, high=N, size=M)
x_f32 = np.float32(np.random.uniform(low=-100.,high=100.,size=N))
# test coverage for elements > 117435.992f for which glibc is used
x_f32[index] = np.float32(10E+10*np.random.rand(M))
x_f64 = np.float64(x_f32)
> assert_array_max_ulp(np.sin(x_f32), np.float32(np.sin(x_f64)), maxulp=2)
E AssertionError: Arrays are not almost equal up to 2 ULP (max
difference is 65 ULP)
M = 5
N = 100
index = array([121958, 671155, 131932, ..., 738271, 310195, 233966])
self =
x_f32 = array([-10.577719, -35.353283, -97.29114 , ..., -80.99214
, -42.875526,
-87.8052 ], dtype=float32)
x_f64 = array([-10.57771873, -35.35328293, -97.2911377 , ..., -80.99214172,
-42.87552643, -87.80519867])
On Sun, Jan 9, 2022 at 4:57 PM Greg Steuck wrote:
>
> Daniel Dickman writes:
>
> > Here's the link to the commit Mark referenced:
> > https://github.com/NetBSD/src/commit/4f9e11b0dddf04640fe0553a9133a471af613627
> >
> > And then the actual implementations were removed in this commit:
> > https://github.com/NetBSD/src/commit/870f792ccadb412e522f37caec6028b0076a871b
> >
> > So I guess this is the list of functions to remove, Mark?
> >
> > I'm testing this on i386 with numpy to see if regress tests improve.
> >
> > s_cos.S
> > s_cosf.S
> > s_sin.S
> > s_sinf.S
> > s_tan.S
> > s_tanf.S
>
> This improves one of the regress tests mbuhl@ added. This one now
> passes:
>
> modified regress/lib/libm/msun/Makefile
> @@ -58,7 +58,6 @@ FAILING+= run-ctrig_test-1
> FAILING+= run-exponential_test-1
> FAILING+= run-invtrig_test-7
> FAILING+= run-next_test-{1,2,4}
> -FAILING+= run-trig_test-3
> . elif ${MACHINE} == arm64
> FAILING+= run-cexp_test-{1,7}
> FAILING+= run-ctrig_test-{1,5}
>
> But one new test is reported broken:
>
> run-trig_test-2
> 2 tests the accuracy of these functions over the primary range
> ./trig_test -r 2
>
> trig_test.c:257: 'fpequal_cs(_x, _y, 1)' evaluated to false
> *** Error 1 in . (Makefile:143 'run-trig_test-2')
> FAILED
>
> So, some attention is needed in this area.
>
> Thanks
> Greg
Re: tanf returns NaN for large inputs
Daniel Dickman writes:
> Here's the link to the commit Mark referenced:
> https://github.com/NetBSD/src/commit/4f9e11b0dddf04640fe0553a9133a471af613627
>
> And then the actual implementations were removed in this commit:
> https://github.com/NetBSD/src/commit/870f792ccadb412e522f37caec6028b0076a871b
>
> So I guess this is the list of functions to remove, Mark?
>
> I'm testing this on i386 with numpy to see if regress tests improve.
>
> s_cos.S
> s_cosf.S
> s_sin.S
> s_sinf.S
> s_tan.S
> s_tanf.S
This improves one of the regress tests mbuhl@ added. This one now
passes:
modified regress/lib/libm/msun/Makefile
@@ -58,7 +58,6 @@ FAILING+= run-ctrig_test-1
FAILING+= run-exponential_test-1
FAILING+= run-invtrig_test-7
FAILING+= run-next_test-{1,2,4}
-FAILING+= run-trig_test-3
. elif ${MACHINE} == arm64
FAILING+= run-cexp_test-{1,7}
FAILING+= run-ctrig_test-{1,5}
But one new test is reported broken:
run-trig_test-2
2 tests the accuracy of these functions over the primary range
./trig_test -r 2
trig_test.c:257: 'fpequal_cs(_x, _y, 1)' evaluated to false
*** Error 1 in . (Makefile:143 'run-trig_test-2')
FAILED
So, some attention is needed in this area.
Thanks
Greg
Re: tanf returns NaN for large inputs
> From: Daniel Dickman > Date: Sun, 9 Jan 2022 16:36:33 -0500 > > On Sun, Jan 9, 2022 at 4:18 PM Mark Kettenis wrote: > > > > > From: Greg Steuck > > > Date: Sun, 09 Jan 2022 12:47:14 -0800 > > > > > > Greg Steuck writes: > > > > > > > This was reduced from a ghc test. The results of the program differ > > > > between OpenBSD 7.0-current-amd64 and a couple of other systems: > > > > > > Thanks to phessler@ for testing on arm64 where the bug doesn't happen. > > > This patch makes OpenBSD-amd64 work the rest of the systems. I added > > > i386 as it was also similarly broken (but didn't test the change yet). > > > > As I said on icb, NetBSD removed all of the x86 assembly sin/cos/tan > > implementations because: > > > > "The x87 hardware uses a bad approximation to pi for argument reduction" > > > > So I think we should use the software fallbacks for all of these > > functions (and remove the broken assembly implementations). > > > > I don't think it makes sense to just remove tanf() and leave the > > others in place. > > > > > > Here's the link to the commit Mark referenced: > https://github.com/NetBSD/src/commit/4f9e11b0dddf04640fe0553a9133a471af613627 > > And then the actual implementations were removed in this commit: > https://github.com/NetBSD/src/commit/870f792ccadb412e522f37caec6028b0076a871b > > So I guess this is the list of functions to remove, Mark? Yes. > I'm testing this on i386 with numpy to see if regress tests improve. > > s_cos.S > s_cosf.S > s_sin.S > s_sinf.S > s_tan.S > s_tanf.S >
Re: tanf returns NaN for large inputs
On Sun, Jan 9, 2022 at 4:18 PM Mark Kettenis wrote: > > > From: Greg Steuck > > Date: Sun, 09 Jan 2022 12:47:14 -0800 > > > > Greg Steuck writes: > > > > > This was reduced from a ghc test. The results of the program differ > > > between OpenBSD 7.0-current-amd64 and a couple of other systems: > > > > Thanks to phessler@ for testing on arm64 where the bug doesn't happen. > > This patch makes OpenBSD-amd64 work the rest of the systems. I added > > i386 as it was also similarly broken (but didn't test the change yet). > > As I said on icb, NetBSD removed all of the x86 assembly sin/cos/tan > implementations because: > > "The x87 hardware uses a bad approximation to pi for argument reduction" > > So I think we should use the software fallbacks for all of these > functions (and remove the broken assembly implementations). > > I don't think it makes sense to just remove tanf() and leave the > others in place. > > Here's the link to the commit Mark referenced: https://github.com/NetBSD/src/commit/4f9e11b0dddf04640fe0553a9133a471af613627 And then the actual implementations were removed in this commit: https://github.com/NetBSD/src/commit/870f792ccadb412e522f37caec6028b0076a871b So I guess this is the list of functions to remove, Mark? I'm testing this on i386 with numpy to see if regress tests improve. s_cos.S s_cosf.S s_sin.S s_sinf.S s_tan.S s_tanf.S
Re: tanf returns NaN for large inputs
> From: Greg Steuck
> Date: Sun, 09 Jan 2022 12:47:14 -0800
>
> Greg Steuck writes:
>
> > This was reduced from a ghc test. The results of the program differ
> > between OpenBSD 7.0-current-amd64 and a couple of other systems:
>
> Thanks to phessler@ for testing on arm64 where the bug doesn't happen.
> This patch makes OpenBSD-amd64 work the rest of the systems. I added
> i386 as it was also similarly broken (but didn't test the change yet).
As I said on icb, NetBSD removed all of the x86 assembly sin/cos/tan
implementations because:
"The x87 hardware uses a bad approximation to pi for argument reduction"
So I think we should use the software fallbacks for all of these
functions (and remove the broken assembly implementations).
I don't think it makes sense to just remove tanf() and leave the
others in place.
> diff --git a/lib/libm/Makefile b/lib/libm/Makefile
> index 47cd94cac06..552e97ea0d3 100644
> --- a/lib/libm/Makefile
> +++ b/lib/libm/Makefile
> @@ -26,7 +26,7 @@ ARCH_SRCS = e_acos.S e_asin.S e_atan2.S e_exp.S e_fmod.S
> e_log.S e_log10.S \
> s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
> s_llrint.S s_llrintf.S s_lrint.S s_lrintf.S s_rint.S s_rintf.S\
> s_scalbnf.S s_significand.S s_significandf.S \
> - s_sin.S s_sinf.S s_tan.S s_tanf.S
> + s_sin.S s_sinf.S s_tan.S
> .elif (${MACHINE_ARCH} == "amd64")
> .PATH: ${.CURDIR}/arch/amd64
> CPPFLAGS+=-I${.CURDIR}/arch/amd64
> @@ -39,7 +39,7 @@ ARCH_SRCS = e_acos.S e_asin.S e_atan2.S e_exp.S e_fmod.S
> e_log.S e_log10.S \
> s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
> s_llrint.S s_llrintf.S s_lrint.S s_lrintf.S \
> s_rint.S s_rintf.S s_scalbnf.S s_significand.S \
> - s_significandf.S s_sin.S s_sinf.S s_tan.S s_tanf.S
> + s_significandf.S s_sin.S s_sinf.S s_tan.S
> .elif (${MACHINE_ARCH} == "hppa")
> .PATH: ${.CURDIR}/arch/hppa
> ARCH_SRCS = e_sqrt.c e_sqrtf.c e_remainder.c e_remainderf.c \
>
>
Re: tanf returns NaN for large inputs
Greg Steuck writes:
> This was reduced from a ghc test. The results of the program differ
> between OpenBSD 7.0-current-amd64 and a couple of other systems:
Thanks to phessler@ for testing on arm64 where the bug doesn't happen.
This patch makes OpenBSD-amd64 work the rest of the systems. I added
i386 as it was also similarly broken (but didn't test the change yet).
diff --git a/lib/libm/Makefile b/lib/libm/Makefile
index 47cd94cac06..552e97ea0d3 100644
--- a/lib/libm/Makefile
+++ b/lib/libm/Makefile
@@ -26,7 +26,7 @@ ARCH_SRCS = e_acos.S e_asin.S e_atan2.S e_exp.S e_fmod.S
e_log.S e_log10.S \
s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
s_llrint.S s_llrintf.S s_lrint.S s_lrintf.S s_rint.S s_rintf.S\
s_scalbnf.S s_significand.S s_significandf.S \
- s_sin.S s_sinf.S s_tan.S s_tanf.S
+ s_sin.S s_sinf.S s_tan.S
.elif (${MACHINE_ARCH} == "amd64")
.PATH: ${.CURDIR}/arch/amd64
CPPFLAGS+=-I${.CURDIR}/arch/amd64
@@ -39,7 +39,7 @@ ARCH_SRCS = e_acos.S e_asin.S e_atan2.S e_exp.S e_fmod.S
e_log.S e_log10.S \
s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
s_llrint.S s_llrintf.S s_lrint.S s_lrintf.S \
s_rint.S s_rintf.S s_scalbnf.S s_significand.S \
- s_significandf.S s_sin.S s_sinf.S s_tan.S s_tanf.S
+ s_significandf.S s_sin.S s_sinf.S s_tan.S
.elif (${MACHINE_ARCH} == "hppa")
.PATH: ${.CURDIR}/arch/hppa
ARCH_SRCS = e_sqrt.c e_sqrtf.c e_remainder.c e_remainderf.c \
tanf returns NaN for large inputs
This was reduced from a ghc test. The results of the program differ
between OpenBSD 7.0-current-amd64 and a couple of other systems:
% cat tanf.c
#include
#include
int main(int a, char**b) {
float x = 1e18;
printf("tanf(%f)=%f\n", x, tanf(x));
float y = 1e19;
printf("tanf(%f)=%f\n", y, tanf(y));
}
% uname -sr; cc tanf.c -o tanf -lm && ./tanf
FreeBSD 13.0-STABLE
tanf(99984306749440.00)=-0.222015
tanf(80506447872.00)=0.708482
% uname -sr; cc tanf.c -o tanf -lm; ./tanf
Linux 5.11.0-38-generic
tanf(99984306749440.00)=-0.222015
tanf(80506447872.00)=0.708482
% uname -sr; cc tanf.c -o tanf -lm && ./tanf
OpenBSD 7.0
tanf(99984306749440.00)=-0.220665
tanf(80506447872.00)=-nan
Notice also the precision loss starting at 1e18.
This behavior has likely been broken for a long time as I remember the
original ghc test to fail last year too.
Thanks
Greg
