ref: e5833209e4cb397f5d17be79d81454b48002882a
parent: eff0ab25c1d542934e1a5d9df45da403777dd3c6
author: S. Gilles <sgilles@math.umd.edu>
date: Fri Jun 7 19:58:10 EDT 2019
Move hi/lo multiplication and addition routines to util.
--- a/lib/math/log-overkill.myr
+++ b/lib/math/log-overkill.myr
@@ -388,47 +388,3 @@
-> (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))
-}
--- a/lib/math/util.myr
+++ b/lib/math/util.myr
@@ -19,6 +19,12 @@
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)
@@ -249,6 +255,50 @@
var c : flt32 = (cL : 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