diff options
author | S. Gilles <sgilles@math.umd.edu> | 2019-06-08 15:30:55 -0400 |
---|---|---|
committer | S. Gilles <sgilles@math.umd.edu> | 2019-06-09 04:21:50 -0400 |
commit | 2002d2fa999d7efeb71aa64a32a1a87c2f4b8399 (patch) | |
tree | 01b1e71d3dd5ec94297cd8d8b877ba884e67c383 /lib | |
parent | 64c5e8f1227c1a98ab996e95425ab43791c0ee2e (diff) | |
download | mc-2002d2fa999d7efeb71aa64a32a1a87c2f4b8399.tar.gz |
Apply changes of pown to rootn. Faster, better edge handling.
Diffstat (limited to 'lib')
-rw-r--r-- | lib/math/ancillary/ulp.gp | 4 | ||||
-rw-r--r-- | lib/math/log-overkill.myr | 2 | ||||
-rw-r--r-- | lib/math/pown-impl.myr | 79 | ||||
-rw-r--r-- | lib/math/test/pown-impl.myr | 6 |
4 files changed, 74 insertions, 17 deletions
diff --git a/lib/math/ancillary/ulp.gp b/lib/math/ancillary/ulp.gp index 9a4b155..c4861b4 100644 --- a/lib/math/ancillary/ulp.gp +++ b/lib/math/ancillary/ulp.gp @@ -31,7 +31,7 @@ e = bitand(a, 0x7f800000) >> 23; s = bitand(a, 0x007fffff); - if(e != 0, s = bitor(s, 0x00800000),); + if(e != 0, s = bitor(s, 0x00800000), s = 2.0 * s); s = s * 2.0^(-23); e = e - 127; return((-1)^n * s * 2^(e)); @@ -42,7 +42,7 @@ e = bitand(a, 0x7ff0000000000000) >> 52; s = bitand(a, 0x000fffffffffffff); - if(e != 0, s = bitor(s, 0x0010000000000000),); + if(e != 0, s = bitor(s, 0x0010000000000000), s = 2.0 * s); s = s * 2.0^(-52); e = e - 1023; return((-1)^n * 2^(e) * s); diff --git a/lib/math/log-overkill.myr b/lib/math/log-overkill.myr index 6c74486..dcc2569 100644 --- a/lib/math/log-overkill.myr +++ b/lib/math/log-overkill.myr @@ -40,8 +40,6 @@ pkg math = pkglocal const logoverkill64 : (x : flt64 -> (flt64, flt64)) ;; - - /* Ci is a table such that, for Ci[j] = (L1, L2, I1, I2), L1, L2 are log(1 + j·2^-(5i)) diff --git a/lib/math/pown-impl.myr b/lib/math/pown-impl.myr index 7e6c5b8..6dfec19 100644 --- a/lib/math/pown-impl.myr +++ b/lib/math/pown-impl.myr @@ -38,6 +38,8 @@ type fltdesc(@f, @u, @i) = struct floor : (x : @f -> @f) log_overkill : (x : @f -> (@f, @f)) precision : @i + nosgn_mask : @u + implicit_bit : @u emin : @i emax : @i imax : @i @@ -62,6 +64,8 @@ const desc32 : fltdesc(flt32, uint32, int32) = [ .floor = floor32, .log_overkill = logoverkill32, .precision = 24, + .nosgn_mask = 0x7fffffff, + .implicit_bit = 23, .emin = -126, .emax = 127, .imax = 2147483647, /* For detecting overflow in final exponent */ @@ -86,6 +90,8 @@ const desc64 : fltdesc(flt64, uint64, int64) = [ .floor = floor64, .log_overkill = logoverkill64, .precision = 53, + .nosgn_mask = 0x7fffffffffffffff, + .implicit_bit = 52, .emin = -1022, .emax = 1023, .imax = 9223372036854775807, @@ -181,6 +187,9 @@ generic powngen = {x : @f, n : @i, d : fltdesc(@f, @u, @i) :: numeric,floating,s But first: do some rough calculations: if we can show n*log(xs) has the same sign as n*e, and n*e would cause overflow, then we might as well return right now. + + This also takes care of subnormals very nicely, so we don't have to do + any special handling to reconstitute xs "right", as we do in rootn. */ var exp_rough_estimate = n * xe if n > 0 && (exp_rough_estimate > d.emax + 1 || (exp_rough_estimate / n != xe)) @@ -286,6 +295,19 @@ generic rootngen = {x : @f, q : @u, d : fltdesc(@f, @u, @i) :: numeric,floating, elif q == 1 /* Anything^1/1 is itself */ -> x + elif xe < d.emin + /* + Subnormals are actually a problem. If we naively reconstitute xs, it + will be wildly wrong and won't match up with the exponent. So let's + pretend we have unbounded exponent range. We know the loop terminates + because we covered the +/-0.0 case above. + */ + xe++ + var check = 1 << d.implicit_bit + while xs & check == 0 + xs <<= 1 + xe-- + ;; ;; /* As in pown */ @@ -294,6 +316,33 @@ generic rootngen = {x : @f, q : @u, d : fltdesc(@f, @u, @i) :: numeric,floating, ult_sgn = -1.0 ;; + /* + If we're looking at (1 + 2^-h)^1/q, and the answer will be 1 + e, with + (1 + e)^q = 1 + 2^-h, then for q and h large enough, e might be below + the representable range. Specifically, + + (1 + e)^q ≅ 1 + qe + (q choose 2)e^2 + ... + + So (using single-precision as the example) + + (1 + 2^-23)^q ≅ 1 + q 2^-23 + (absolutely tiny terms) + + And anything in [1, 1 + q 2^-24) will just truncate to 1.0 when + calculated. + */ + if xe == 0 + var cutoff = scale2(qf, -1 * d.precision - 1) + 1.0 + if (xb & d.nosgn_mask) < d.tobits(cutoff) + -> 1.0 + ;; + elif xe == -1 + /* Something similar for (1 - e)^q */ + var cutoff = 1.0 - scale2(qf, -1 * d.precision - 1) + if (xb & d.nosgn_mask) > d.tobits(cutoff) + -> 1.0 + ;; + ;; + /* Similar to pown. Let e/q = E + psi, with E an integer. x^(1/q) = e^(log(xs)/q) * 2^(e/q) @@ -307,14 +356,13 @@ generic rootngen = {x : @f, q : @u, d : fltdesc(@f, @u, @i) :: numeric,floating, */ /* Calculate 1/q in very high precision */ - var r1 = 1.0 / qf - var r2 = -math.fma(r1, qf, -1.0) / qf + var qinv_hi = 1.0 / qf + var qinv_lo = -math.fma(qinv_hi, qf, -1.0) / qf var ln_xs_hi, ln_xs_lo (ln_xs_hi, ln_xs_lo) = d.log_overkill(d.assem(false, 0, xs)) - var ls1 : @f[12] - (ls1[0], ls1[1]) = d.two_by_two(ln_xs_hi, r1) - (ls1[2], ls1[3]) = d.two_by_two(ln_xs_hi, r2) - (ls1[4], ls1[5]) = d.two_by_two(ln_xs_lo, r1) + + var G1, G2 + (G1, G2) = d.split_mul(ln_xs_hi, ln_xs_lo, qinv_hi, qinv_lo) var E : @i if q > std.abs(xe) @@ -328,15 +376,20 @@ generic rootngen = {x : @f, q : @u, d : fltdesc(@f, @u, @i) :: numeric,floating, var psi_lo = -math.fma(psi_hi, qf, -(qpsi : @f)) / qf var log2_hi, log2_lo (log2_hi, log2_lo) = d.C[128] - (ls1[ 6], ls1[ 7]) = d.two_by_two(psi_hi, d.frombits(log2_hi)) - (ls1[ 8], ls1[ 9]) = d.two_by_two(psi_hi, d.frombits(log2_lo)) - (ls1[10], ls1[11]) = d.two_by_two(psi_lo, d.frombits(log2_hi)) + var H1, H2 + (H1, H2) = d.split_mul(psi_hi, psi_lo, d.frombits(log2_hi), d.frombits(log2_lo)) - var G1, G2 - (G1, G2) = double_compensated_sum(ls1[0:12]) - /* G1 + G2 approximates log(xs)/q + log(2)*psi */ + var J1, J2, t1, t2 + /* + We can't use split_add; we don't kow the relative magitudes of G and H + */ + (t1, t2) = slow2sum(G2, H2) + (J2, t1) = slow2sum(H1, t1) + (J1, J2) = slow2sum(G1, J2) + J2 = J2 + (t1 + t2) - var base = exp(G1) + G2 + /* J1 + J2 approximates log(xs)/q + log(2)*psi */ + var base = exp(J1) + J2 -> ult_sgn * scale2(base, E) } diff --git a/lib/math/test/pown-impl.myr b/lib/math/test/pown-impl.myr index d9d431f..7308d25 100644 --- a/lib/math/test/pown-impl.myr +++ b/lib/math/test/pown-impl.myr @@ -457,6 +457,12 @@ const rootn02 = {c (0xe0bc4cbf6bd74d8f, 0x000000000000bd8b, 0xbff01ed2c4e821fc), (0x31d4f2baa91a9e8e, 0x000000000000244d, 0x3fef774c954b40bf), (0x01283d1c679f5652, 0x0000000000008647, 0x3fef5bc18f5e292f), + (0x80003d8a341ee060, 0x0000000000009c71, 0xbfef6f873f76b7cd), + (0xbfecf0fc4dc97b93, 0x0000000000005f4b, 0xbfeffff75d0a25fe), + (0xbfe2fb84944a35ee, 0x000000000000e947, 0xbfefffeda94c6d07), + (0xbfe0ef0c05bd84ab, 0x0000000000007165, 0xbfefffd205e2a1c8), + (0xbfe7354e962bdcb3, 0x000000000000076b, 0xbfeffe9d4844aad6), + (0xbfea556dd1eb1e58, 0x00000000000095cb, 0xbfeffff557890356), ][:] for (x, y, z) : inputs |