shithub: mc

Download patch

ref: d7de78142e5a68bce8a951a4567e0621671aab42
parent: c2cead0d7e100e3c6f808af1bfaea25e151c53c0
author: S. Gilles <sgilles@math.umd.edu>
date: Fri Jun 29 22:01:05 EDT 2018 

Be a little smarter about downscaling flt64s to flt32s in sin/cos.


--- a/lib/math/sin-impl.myr
+++ b/lib/math/sin-impl.myr
@@ -403,36 +403,39 @@
 ]
 
 const sin32 = {x : flt32
-	var r = sin64((x : flt64))
-	-> (r : flt32)
+	var s, s2
+	(s, s2, _, _) = w64((x : flt64), true, false)
+	-> round_down(s, s2)
 }
 
 const sin64 = {x : flt64
 	var s
-	(s, _) = w64(x, true, false)
+	(s, _, _, _) = w64(x, true, false)
 	-> s
 }
 
 const cos32 = {x : flt32
-	var r : flt64 = cos64((x : flt64))
-	-> (r : flt32)
+	var c, c2
+	(_, _, c, c2) = w64((x : flt64), false, true)
+	-> round_down(c, c2)
 }
 
 const cos64 = {x : flt64
 	var c
-	(_, c) = w64(x, false, true)
+	(_, _, c, _) = w64(x, false, true)
 	-> c
 }
 
 const sincos32 = {x : flt32
-	var r1 : flt64
-	var r2 : flt64
-	(r1, r2) = sincos64((x : flt64))
-	-> ((r1 : flt32), (r2 : flt32))
+	var s : flt64, s2 : flt64, c : flt64, c2 : flt64
+	(s, s2, c, c2) = w64((x : flt64), true, true)
+	-> (round_down(s, s2), round_down(c, c2))
 }
 
 const sincos64 = {x : flt64
-	-> w64(x, true, true)
+	var s, c
+	(s, _, c, _) = w64(x, true, true)
+	-> (s, c)
 }
 
 /* Calculate sin and/or cos */
@@ -439,12 +442,14 @@
 const w64 = {x : flt64, want_sin : bool, want_cos : bool
 	var sin_ret : flt64 = 0.0
 	var cos_ret : flt64 = 0.0
+	var sin_ret2 : flt64 = 0.0
+	var cos_ret2 : flt64 = 0.0
 
 	var e : int64
 	(_, e, _) = std.flt64explode(x)
 
 	if e == 1024
-		-> (std.flt64nan(), std.flt64nan())
+		-> (std.flt64nan(), 0.0, std.flt64nan(), 0.0)
 	;;
 
 	var N : int64
@@ -542,7 +547,7 @@
 		std.sort(q[:], fltabscmp)
 
 		/* TODO: use double-compensation instead */
-		(sin_ret, _, _) = triple_compensate_sum(q[:])
+		(sin_ret, sin_ret2, _) = triple_compensate_sum(q[:])
 	;;
 
 	if (want_cos && !swap_sin_cos) || (want_sin && swap_sin_cos)
@@ -553,26 +558,30 @@
 		std.sort(q[:], fltabscmp)
 
 		/* TODO: use double-compensation instead */
-		(cos_ret, _, _) = triple_compensate_sum(q[:])
+		(cos_ret, cos_ret2, _) = triple_compensate_sum(q[:])
 	;;
 
 	if first_negate_sin
 		sin_ret = -sin_ret
+		sin_ret2 = -sin_ret2
 	;;
 
 	if swap_sin_cos
 		std.swap(&sin_ret, &cos_ret)
+		std.swap(&sin_ret2, &cos_ret2)
 	;;
 
 	if then_negate_sin
 		sin_ret = -sin_ret
+		sin_ret2 = -sin_ret2
 	;;
 
 	if then_negate_cos
 		cos_ret = -cos_ret
+		cos_ret2 = -cos_ret2
 	;;
 
-	-> (sin_ret, cos_ret)
+	-> (sin_ret, sin_ret2, cos_ret, cos_ret2)
 }
 
 /* Reduce x to N*(pi/4) + x', with x' in [-pi/4, pi/4] */
@@ -834,4 +843,22 @@
 	var t : flt64 = b - z
 	-> (s, t)
 	
+}
+
+/*
+   Round a + b to a flt32. Only notable if round(a) is a rounding
+   tie, and b is non-zero
+ */
+const round_down = {a : flt64, b : flt64
+	var au : uint64 = std.flt64bits(a)
+	if au & 0x0000000070000000 == 0x0000000070000000
+		if b > 0.0
+			au++
+		elif b < 0.0
+			au--
+		;;
+		-> (std.flt64frombits(au) : flt32)
+	;;
+
+	-> (a : flt32)
 }