ref: 34adf4b62333b2bfd3496d3b61d62bc7a00f3140
parent: 70b4fc9127a3e8863f128642d62d2fbcbf7121be
author: S. Gilles <sgilles@math.umd.edu>
date: Thu Jun 28 10:51:24 EDT 2018
Fix local reduction from [-pi/2, pi/2] to [-pi/4, pi/4]
--- a/lib/math/ancillary/pi-constants.c
+++ b/lib/math/ancillary/pi-constants.c
@@ -60,6 +60,19 @@
}
printf("\n");+ printf("1000 bits of pi/4:\n");+ mpfr_const_pi(pi, MPFR_RNDN);
+ mpfr_div_si(pi, pi, 4, MPFR_RNDN);
+
+ for (int bits_obtained = 0; bits_obtained < 1000; bits_obtained += 53) {+ d = mpfr_get_d(pi, MPFR_RNDN);
+ u = FLT64_TO_UINT64(d);
+ printf("%#018lx\n", u);+ mpfr_set_d(t, d, MPFR_RNDN);
+ mpfr_sub(pi, pi, t, MPFR_RNDN);
+ }
+
+ printf("\n"); printf("Pre-computed 2/pi:\n");mpfr_const_pi(pi, MPFR_RNDN);
mpfr_set_si(t, 2, MPFR_RNDN);
--- a/lib/math/sin-impl.myr
+++ b/lib/math/sin-impl.myr
@@ -58,6 +58,14 @@
0x35b4cf98e804177d,
]
+/* Pi/4 in lots of detail */
+const pi_over_4 : uint64[4] = [
+ 0x3fe921fb54442d18,
+ 0x3c81a62633145c07,
+ 0xb90f1976b7ed8fbc,
+ 0x35a4cf98e804177d,
+]
+
/* Pre-computed inverses */
const two_over_pi : uint64 = 0x3fe45f306dc9c883 /* 1/(pi/2) */
const oneohtwofour_over_pi : uint64 = 0x40745f306dc9c883 /* 1/(pi/(4 * 256)) */
@@ -584,6 +592,10 @@
var pi_o_2_1 : flt64 = std.flt64frombits(pi_over_2[1])
var pi_o_2_2 : flt64 = std.flt64frombits(pi_over_2[2])
var pi_o_2_3 : flt64 = std.flt64frombits(pi_over_2[3])
+ var pi_o_4_0 : flt64 = std.flt64frombits(pi_over_4[0])
+ var pi_o_4_1 : flt64 = std.flt64frombits(pi_over_4[1])
+ var pi_o_4_2 : flt64 = std.flt64frombits(pi_over_4[2])
+ var pi_o_4_3 : flt64 = std.flt64frombits(pi_over_4[3])
var total_N = 0
var q : flt64[11]
@@ -598,8 +610,8 @@
(x1, x2, x3) = triple_compensate_sum(q[:])
total_N += (N % 4)
- if le_33(x1, x2, x3, pi_o_2_0, pi_o_2_1, pi_o_2_2) && \
- le_33(-pi_o_2_0, -pi_o_2_1, -pi_o_2_2, x1, x2, x3)
+ if le_33(x1, x2, x3, pi_o_4_0, pi_o_4_1, pi_o_4_2) && \
+ le_33(-pi_o_4_0, -pi_o_4_1, -pi_o_4_2, x1, x2, x3)
break
;;
;;
--- a/lib/math/test/sin-impl.myr
+++ b/lib/math/test/sin-impl.myr
@@ -56,7 +56,6 @@
var rsu2 = std.flt32bits(rsf2)
var rcu2 = std.flt32bits(rcf2)
-
testr.check(c, rsf1 == rsf2 && rcf1 == rcf2,
"sincos(0x{b=16,w=8,p=0}) is (0x{b=16,w=8,p=0}, 0x{b=16,w=8,p=0}), individual results (0x{b=16,w=8,p=0}, 0x{b=16,w=8,p=0})",x, rsu1, rcu1, rsu2, rcu2)
@@ -73,7 +72,7 @@
(0x4100000000000000, 0xbfeff8bd7b10d6b0, 0x3fa58ced65ec8b50),
(0x4b01000000000000, 0xbfe3e9527dc75f12, 0x3fe90cf80997c963),
(0x4b11000000000000, 0xbfef2cb48ed49aa6, 0x3fcce246843789ad),
- //(0x020400000a0c0000, 0x020400000a0c0000, 0x3ff0000000000000),
+ (0x020400000a0c0000, 0x020400000a0c0000, 0x3ff0000000000000),
][:]
for (x, ys, yc) : inputs
@@ -88,7 +87,6 @@
var rsu2 = std.flt64bits(rsf2)
var rcu2 = std.flt64bits(rcf2)
-
testr.check(c, rsf1 == rsf2 && rcf1 == rcf2,
"sincos(0x{b=16,w=16,p=0}) is (0x{b=16,w=16,p=0}, 0x{b=16,w=16,p=0}), individual results (0x{b=16,w=16,p=0}, 0x{b=16,w=16,p=0})",x, rsu1, rcu1, rsu2, rcu2)
@@ -100,6 +98,31 @@
}
const sincos03 = {c+ var inputs : (uint64, uint64, uint64, uint64, uint64)[:] = [
+ (0x5101000000000000, 0x3fe9706123d509f1, 0xbfe369af9695aba1, 0x3fe9706123d509f0, 0xbfe369af9695aba0),
+ ][:]
+
+ for (x, ys_perfect, yc_perfect, ys_acceptable, yc_acceptable) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var rsf1, rcf1, rsf2, rcf2
+ (rsf1, rcf1) = math.sincos(xf)
+ rsf2 = math.sin(xf)
+ rcf2 = math.cos(xf)
+
+ var rsu1 = std.flt64bits(rsf1)
+ var rcu1 = std.flt64bits(rcf1)
+ var rsu2 = std.flt64bits(rsf2)
+ var rcu2 = std.flt64bits(rcf2)
+
+ testr.check(c, rsf1 == rsf2 && rcf1 == rcf2,
+ "sincos(0x{b=16,w=16,p=0}) is (0x{b=16,w=16,p=0}, 0x{b=16,w=16,p=0}), individual results (0x{b=16,w=16,p=0}, 0x{b=16,w=16,p=0})",+ x, rsu1, rcu1, rsu2, rcu2)
+
+ testr.check(c, (same64(rsu1, ys_perfect) || same64(rsu1, ys_acceptable)) && \
+ (same64(rcu1, yc_perfect) || same64(rcu1, yc_acceptable)),
+ "sincos(0x{b=16,w=16,p=0}) should be (0x{b=16,w=16,p=0}, 0x{b=16,w=16,p=0}), will also accept (0x{b=16,w=16,p=0}, 0x{b=16,w=16,p=0}), was (0x{b=16,w=16,p=0}, 0x{b=16,w=16,p=0})",+ x, ys_perfect, yc_perfect, ys_acceptable, yc_acceptable, rsu1, rcu1)
+ ;;
}
const sincos04 = {c