diff options
author | S. Gilles <sgilles@math.umd.edu> | 2019-06-07 23:58:10 -0400 |
---|---|---|
committer | S. Gilles <sgilles@math.umd.edu> | 2019-06-07 23:58:10 -0400 |
commit | e5833209e4cb397f5d17be79d81454b48002882a (patch) | |
tree | 36abb5af5acf54232d9ea860ddc37715aa871280 | |
parent | eff0ab25c1d542934e1a5d9df45da403777dd3c6 (diff) | |
download | mc-e5833209e4cb397f5d17be79d81454b48002882a.tar.gz |
Move hi/lo multiplication and addition routines to util.
-rw-r--r-- | lib/math/log-overkill.myr | 44 | ||||
-rw-r--r-- | lib/math/util.myr | 50 |
2 files changed, 50 insertions, 44 deletions
diff --git a/lib/math/log-overkill.myr b/lib/math/log-overkill.myr index 72e52c7..ecac77d 100644 --- a/lib/math/log-overkill.myr +++ b/lib/math/log-overkill.myr @@ -388,47 +388,3 @@ const logoverkill64 = {x : flt64 -> (lx_hi, lx_lo) } - -const hl_mult = {a_h : flt64, a_l : flt64, b_h : flt64, b_l : flt64 - /* - [ a_h ][ a_l ] * [ b_h ][ b_l ] - = - (A) [ a_h*b_h ] - (B) + [ a_h*b_l ] - (C) + [ a_l*b_h ] - (D) + [ a_l*b_l ] - - We therefore want to keep all of A, and the top halves of the two - smaller products B and C. - - To be pedantic, *_l could be much smaller than pictured above; *_h and - *_l need not butt up right against each other. But that won't cause - any problems; there's no way we'll ignore important information. - */ - var Ah, Al - (Ah, Al) = two_by_two64(a_h, b_h) - var Bh = a_h * b_l - var Ch = a_l * b_h - var resh, resl, t1, t2 - (resh, resl) = slow2sum(Bh, Ch) - (resl, t1) = fast2sum(Al, resl) - (resh, t2) = slow2sum(resh, resl) - (resh, resl) = slow2sum(Ah, resh) - -> (resh, resl + (t1 + t2)) -} - -const hl_add = {a_h : flt64, a_l : flt64, b_h : flt64, b_l : flt64 - /* - Not terribly clever, we just chain a couple of 2sums together. We are - free to impose the requirement that |a_h| > |b_h|, because we'll only - be using this for a = 1/5, 1/3, and the log(Fi)s from the C tables. - However, we can't guarantee that a_l > b_l. For example, compare C1[10] - to C2[18]. - */ - - var resh, resl, t1, t2 - (t1, t2) = slow2sum(a_l, b_l) - (resl, t1) = slow2sum(b_h, t1) - (resh, resl) = fast2sum(a_h, resl) - -> (resh, resl + (t1 + t2)) -} diff --git a/lib/math/util.myr b/lib/math/util.myr index 8223d13..b36715a 100644 --- a/lib/math/util.myr +++ b/lib/math/util.myr @@ -19,6 +19,12 @@ pkg math = const two_by_two64 : (x : flt64, y : flt64 -> (flt64, flt64)) const two_by_two32 : (x : flt32, y : flt32 -> (flt32, flt32)) + /* Multiply (a_hi + a_lo) * (b_hi + b_lo) to (z_hi + z_lo) */ + const hl_mult : (a_h : flt64, a_l : flt64, b_h : flt64, b_l : flt64 -> (flt64, flt64)) + + /* Add (a_hi + a_lo) * (b_hi + b_lo) to (z_hi + z_lo). Must have |a| > |b|. */ + const hl_add : (a_h : flt64, a_l : flt64, b_h : flt64, b_l : flt64 -> (flt64, flt64)) + /* Compare by magnitude */ const mag_cmp32 : (f : flt32, g : flt32 -> std.order) const mag_cmp64 : (f : flt64, g : flt64 -> std.order) @@ -251,6 +257,50 @@ const two_by_two32 = {x : flt32, y : flt32 -> (s, c) } +const hl_mult = {a_h : flt64, a_l : flt64, b_h : flt64, b_l : flt64 + /* + [ a_h ][ a_l ] * [ b_h ][ b_l ] + = + (A) [ a_h*b_h ] + (B) + [ a_h*b_l ] + (C) + [ a_l*b_h ] + (D) + [ a_l*b_l ] + + We therefore want to keep all of A, and the top halves of the two + smaller products B and C. + + To be pedantic, *_l could be much smaller than pictured above; *_h and + *_l need not butt up right against each other. But that won't cause + any problems; there's no way we'll ignore important information. + */ + var Ah, Al + (Ah, Al) = two_by_two64(a_h, b_h) + var Bh = a_h * b_l + var Ch = a_l * b_h + var resh, resl, t1, t2 + (resh, resl) = slow2sum(Bh, Ch) + (resl, t1) = fast2sum(Al, resl) + (resh, t2) = slow2sum(resh, resl) + (resh, resl) = slow2sum(Ah, resh) + -> (resh, resl + (t1 + t2)) +} + +const hl_add = {a_h : flt64, a_l : flt64, b_h : flt64, b_l : flt64 + /* + Not terribly clever, we just chain a couple of 2sums together. We are + free to impose the requirement that |a_h| > |b_h|, because we'll only + be using this for a = 1/5, 1/3, and the log(Fi)s from the C tables. + However, we can't guarantee that a_l > b_l. For example, compare C1[10] + to C2[18]. + */ + + var resh, resl, t1, t2 + (t1, t2) = slow2sum(a_l, b_l) + (resl, t1) = slow2sum(b_h, t1) + (resh, resl) = fast2sum(a_h, resl) + -> (resh, resl + (t1 + t2)) +} + const mag_cmp32 = {f : flt32, g : flt32 var u = std.flt32bits(f) & ~(1 << 31) var v = std.flt32bits(g) & ~(1 << 31) |