summaryrefslogtreecommitdiff log msg author committer range
path: root/lib/math
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
author committer S. Gilles 2018-06-30 17:36:51 -0400 S. Gilles 2018-06-30 19:36:32 -0400 9afaec2ee7ed8864f9f28c39533e2f2094301fba (patch) 3caa10552c66002f2acf60b2a717f369f3cfdbc5 /lib/math 09e7cfb56388b1ea79ee7e100e52cb5b34e6017c (diff) mc-9afaec2ee7ed8864f9f28c39533e2f2094301fba.tar.gz
Eliminate loop sorts in sin/cos reduction.
Diffstat (limited to 'lib/math')
-rw-r--r--lib/math/sin-impl.myr29
1 files changed, 18 insertions, 11 deletions
 diff --git a/lib/math/sin-impl.myr b/lib/math/sin-impl.myrindex 7812b7d..db9d47c 100644--- a/lib/math/sin-impl.myr+++ b/lib/math/sin-impl.myr@@ -624,23 +624,32 @@ const reduce = {x : flt64 /* Compute initial reduction -- this might not be sufficient. */ total_N = rn(x1 * std.flt64frombits(two_over_pi)) Nf = (-total_N : flt64)- (q[0], q[ 1], q[2]) = (x1, x2, x3)- (q[3], q[ 4]) = two_by_two(Nf, pi_o_2_0)- (q[5], q[ 6]) = two_by_two(Nf, pi_o_2_1)+ (q[0], q[ 4], q[5]) = (x1, x2, x3)+ (q[1], q[ 3]) = two_by_two(Nf, pi_o_2_0)+ (q[2], q[ 6]) = two_by_two(Nf, pi_o_2_1) (q[7], q[ 8]) = two_by_two(Nf, pi_o_2_2) (q[9], q[10]) = two_by_two(Nf, pi_o_2_3)- std.sort(q[:], fltabscmp)++ /*+ Sorting is very slow, but it's only the top five or so+ that are in question+ */+ std.sort(q[0:5], fltabscmp) (x1, x2) = double_compensated_sum(q[:]) while !(le_22(x1, x2, pi_o_4_0, pi_o_4_1) && le_22(-pi_o_4_0, -pi_o_4_1, x1, x2)) N = rn(x1 * std.flt64frombits(two_over_pi)) Nf = (-N : flt64)- (q[0], q[1]) = (x1, x2)- (q[2], q[3]) = two_by_two(Nf, pi_o_2_0)++ /*+ Sorting is slow. We know that x1 is roughly+ cancelled by Nf * pi_o_2_0, so line those up.+ */+ (q[0], q[2]) = (x1, x2)+ (q[1], q[3]) = two_by_two(Nf, pi_o_2_0) (q[4], q[5]) = two_by_two(Nf, pi_o_2_1) (q[6], q[7]) = two_by_two(Nf, pi_o_2_2) (q[8], q[9]) = two_by_two(Nf, pi_o_2_3)- std.sort(q[0:10], fltabscmp) (x1, x2) = double_compensated_sum(q[0:10]) total_N += (N % 4) ;;@@ -737,11 +746,10 @@ const huge_reduce_2pi = {x : flt64 (q[14], q[15]) = two_by_two(xc, c2) (q[16], q[17]) = two_by_two(xc, c3) - -> triple_compensate_sum(q[:])+ -> triple_compensated_sum(q[:]) } -/* TODO: move to sum-impl.myr, figure out a good API to allow optional sorting */-const triple_compensate_sum = {q : flt64[:]+const triple_compensated_sum = {q : flt64[:] /* TODO: verify, with GAPPA or something, that this is correct. */ std.sort(q, fltabscmp) var s1 : flt64, s2 : flt64, s3@@ -788,7 +796,6 @@ const two_by_two = {x : flt64, y : flt64 var a4 : flt64 = xl * yl /* By-hand compensated summation */- /* TODO: convert to normal compensated summation, move to sum-impl.myr */ var yy, u, t, v, z, s, c if a2 < a3 std.swap(&a3, &a2)