ref: 3665783a1d1fe03db70cca452c615b439310a329
parent: 0f956d6bb55e0455326f4627024c1fa8c8adbff8
parent: 37ba2ec87069c9ca9bfeac71ec4dba660fcadf60
author: Ori Bernstein <[email protected]>
date: Wed May 9 15:21:16 EDT 2018
Merge remote-tracking branch 'npnth/libmath'
--- a/6/isel.c
+++ b/6/isel.c
@@ -939,8 +939,8 @@
/* These operators should never show up in the reduced trees,
* since they should have been replaced with more primitive
* expressions by now */
- case Obad: case Opreinc: case Opostinc: case Opredec: case Otern:
- case Opostdec: case Olor: case Oland: case Oaddeq:
+ case Obad: case Oauto: case Opreinc: case Opostinc: case Opredec:
+ case Otern: case Opostdec: case Olor: case Oland: case Oaddeq:
case Osubeq: case Omuleq: case Odiveq: case Omodeq: case Oboreq:
case Obandeq: case Obxoreq: case Obsleq: case Obsreq: case Omemb:
case Oslbase: case Osllen: case Ocast: case Outag: case Oudata:
--- a/lib/bld.sub
+++ b/lib/bld.sub
@@ -8,6 +8,7 @@
inifile
iter
json
+ math
regex
std
sys
--- /dev/null
+++ b/lib/math/bld.sub
@@ -1,0 +1,36 @@
+lib math =
+ fpmath.myr
+
+ # exp
+ exp-impl.myr
+
+ # rounding (to actual integers)
+ round-impl+posixy-x64-sse4.s
+ round-impl.myr
+
+ # fused-multiply-add
+ fma-impl+posixy-x64-fma.s
+ fma-impl.myr
+
+ # polynomial evaluation methods
+ poly-impl.myr
+
+ # scalb (multiply x by 2^m)
+ scale2-impl.myr
+
+ # sqrt
+ sqrt-impl+posixy-x64-sse2.s
+ sqrt-impl.myr
+
+ # summation
+ sum-impl.myr
+
+ # trunc, floor, ceil
+ trunc-impl+posixy-x64-sse4.s
+ trunc-impl.myr
+
+ # util
+ util.myr
+
+ lib ../std:std
+;;
--- /dev/null
+++ b/lib/math/exp-impl.myr
@@ -1,0 +1,396 @@
+use std
+
+use "fpmath"
+use "util"
+
+/*
+ See [Mul16] (6.2.1), [Tan89], and [Tan92]. While the exp and
+ expm1 functions are similar enough to be in the same file (Tang's
+ algorithms use the same S constants), they are not quite similar
+ enough to be in the same function.
+ */
+pkg math =
+ pkglocal const exp32 : (f : flt32 -> flt32)
+ pkglocal const exp64 : (f : flt64 -> flt64)
+
+ pkglocal const expm132 : (f : flt32 -> flt32)
+ pkglocal const expm164 : (f : flt64 -> flt64)
+;;
+
+extern const fma32 : (x : flt32, y : flt32, z : flt32 -> flt32)
+extern const fma64 : (x : flt64, y : flt64, z : flt64 -> flt64)
+extern const horner_polyu32 : (f : flt32, a : uint32[:] -> flt32)
+extern const horner_polyu64 : (f : flt64, a : uint64[:] -> flt64)
+
+type fltdesc(@f, @u, @i) = struct
+ explode : (f : @f -> (bool, @i, @u))
+ assem : (n : bool, e : @i, s : @u -> @f)
+ fma : (x : @f, y : @f, z : @f -> @f)
+ horner : (f : @f, a : @u[:] -> @f)
+ sgnmask : @u
+ tobits : (f : @f -> @u)
+ frombits : (u : @u -> @f)
+ inf : @u
+ nan : @u
+
+ /* For exp */
+ thresh_1_min : @u
+ thresh_1_max : @u
+ thresh_2 : @u
+ Ai : @u[:]
+
+ /* For expm1 */
+ thresh_tiny : @u
+ thresh_huge_neg : @u
+ Log3_4 : @u
+ Log5_4 : @u
+ Bi : @u[:]
+ precision : @u
+
+ /* For both */
+ nabs : @u
+ inv_L : @u
+ L1 : @u
+ L2 : @u
+ S : (@u, @u)[32]
+
+;;
+
+const desc32 : fltdesc(flt32, uint32, int32) = [
+ .explode = std.flt32explode,
+ .assem = std.flt32assem,
+ .fma = fma32,
+ .horner = horner_polyu32,
+ .sgnmask = (1 << 31),
+ .tobits = std.flt32bits,
+ .frombits = std.flt32frombits,
+ .inf = 0x7f800000,
+ .nan = 0x7fc00000,
+ .thresh_1_min = 0xc2cff1b4, /* ln(2^(-127 - 23)) ~= -103.9720770839917 */
+ .thresh_1_max = 0x42b17218, /* ln([2 - 2^(-24)] * 2^127) ~= 88.72283905206 */
+ .inv_L = 0x4238aa3b, /* L = log(2) / 32, this is 1/L in working precision */
+ .L1 = 0x3cb17200, /* L1 and L2 sum to give L in high precision, */
+ .L2 = 0x333fbe8e, /* and L1 has some trailing zeroes. */
+ .nabs = 9, /* L1 has 9 trailing zeroes in binary */
+ .Ai = [ /* Coefficients for approximating expm1 in a tiny range */
+ 0x3f000044,
+ 0x3e2aaaec,
+ ][:],
+ .Log3_4 = 0xbe934b11, /* log(1 - 1/4) < x < log(1 + 1/4) */
+ .Log5_4 = 0x3e647fbf, /* triggers special expm1 case */
+ .thresh_tiny = 0x33000000, /* similar to thresh_1_{min,max}, but for expm1 */+ .thresh_huge_neg = 0xc18aa122,
+ .Bi = [ /* Coefficients for approximating expm1 between log(3/4) and log(5/4) */
+ 0x3e2aaaaa,
+ 0x3d2aaaa0,
+ 0x3c0889ff,
+ 0x3ab64de5,
+ 0x394ab327,
+ ][:],
+ .S = [ /* Double-precise expansions of 2^(J/32) */
+ (0x3f800000, 0x00000000),
+ (0x3f82cd80, 0x35531585),
+ (0x3f85aac0, 0x34d9f312),
+ (0x3f889800, 0x35e8092e),
+ (0x3f8b95c0, 0x3471f546),
+ (0x3f8ea400, 0x36e62d17),
+ (0x3f91c3c0, 0x361b9d59),
+ (0x3f94f4c0, 0x36bea3fc),
+ (0x3f9837c0, 0x36c14637),
+ (0x3f9b8d00, 0x36e6e755),
+ (0x3f9ef500, 0x36c98247),
+ (0x3fa27040, 0x34c0c312),
+ (0x3fa5fec0, 0x36354d8b),
+ (0x3fa9a140, 0x3655a754),
+ (0x3fad5800, 0x36fba90b),
+ (0x3fb123c0, 0x36d6074b),
+ (0x3fb504c0, 0x36cccfe7),
+ (0x3fb8fb80, 0x36bd1d8c),
+ (0x3fbd0880, 0x368e7d60),
+ (0x3fc12c40, 0x35cca667),
+ (0x3fc56700, 0x36a84554),
+ (0x3fc9b980, 0x36f619b9),
+ (0x3fce2480, 0x35c151f8),
+ (0x3fd2a800, 0x366c8f89),
+ (0x3fd744c0, 0x36f32b5a),
+ (0x3fdbfb80, 0x36de5f6c),
+ (0x3fe0ccc0, 0x36776155),
+ (0x3fe5b900, 0x355cef90),
+ (0x3feac0c0, 0x355cfba5),
+ (0x3fefe480, 0x36e66f73),
+ (0x3ff52540, 0x36f45492),
+ (0x3ffa8380, 0x36cb6dc9),
+ ],
+ .precision = 24, /* threshold to prevent underflow in a S[high] + 2^-m */
+]
+
+const desc64 : fltdesc(flt64, uint64, int64) = [
+ .explode = std.flt64explode,
+ .assem = std.flt64assem,
+ .fma = fma64,
+ .horner = horner_polyu64,
+ .sgnmask = (1 << 63),
+ .tobits = std.flt64bits,
+ .frombits = std.flt64frombits,
+ .inf = 0x7ff0000000000000,
+ .nan = 0x7ff8000000000000,
+ .thresh_1_min = 0xc0874910d52d3052, /* ln(2^(-1023 - 52)) ~= -745.1332191019 */
+ .thresh_1_max = 0x40862e42fefa39ef, /* ln([2 - 2^(-53)] * 2^1023) ~= 709.78271289 */
+ .inv_L = 0x40471547652b82fe, /* below as in exp32 */
+ .L1 = 0x3f962e42fef00000,
+ .L2 = 0x3d8473de6af278ed,
+ .nabs = 20,
+ .Ai = [
+ 0x3fe0000000000000,
+ 0x3fc5555555548f7c,
+ 0x3fa5555555545d4e,
+ 0x3f811115b7aa905e,
+ 0x3f56c1728d739765,
+ ][:],
+ .Log3_4 = 0xbfd269621134db93,
+ .Log5_4 = 0x3fcc8ff7c79a9a22,
+ .thresh_tiny = 0x3c90000000000000,
+ .thresh_huge_neg = 0xc042b708872320e1,
+ .Bi = [
+ 0x3fc5555555555549,
+ 0x3fa55555555554b6,
+ 0x3f8111111111a9f3,
+ 0x3f56c16c16ce14c6,
+ 0x3f2a01a01159dd2d,
+ 0x3efa019f635825c4,
+ 0x3ec71e14bfe3db59,
+ 0x3e928295484734ea,
+ 0x3e5a2836aa646b96,
+ ][:],
+ .S = [
+ (0x3ff0000000000000, 0x0000000000000000),
+ (0x3ff059b0d3158540, 0x3d0a1d73e2a475b4),
+ (0x3ff0b5586cf98900, 0x3ceec5317256e308),
+ (0x3ff11301d0125b40, 0x3cf0a4ebbf1aed93),
+ (0x3ff172b83c7d5140, 0x3d0d6e6fbe462876),
+ (0x3ff1d4873168b980, 0x3d053c02dc0144c8),
+ (0x3ff2387a6e756200, 0x3d0c3360fd6d8e0b),
+ (0x3ff29e9df51fdec0, 0x3d009612e8afad12),
+ (0x3ff306fe0a31b700, 0x3cf52de8d5a46306),
+ (0x3ff371a7373aa9c0, 0x3ce54e28aa05e8a9),
+ (0x3ff3dea64c123400, 0x3d011ada0911f09f),
+ (0x3ff44e0860618900, 0x3d068189b7a04ef8),
+ (0x3ff4bfdad5362a00, 0x3d038ea1cbd7f621),
+ (0x3ff5342b569d4f80, 0x3cbdf0a83c49d86a),
+ (0x3ff5ab07dd485400, 0x3d04ac64980a8c8f),
+ (0x3ff6247eb03a5580, 0x3cd2c7c3e81bf4b7),
+ (0x3ff6a09e667f3bc0, 0x3ce921165f626cdd),
+ (0x3ff71f75e8ec5f40, 0x3d09ee91b8797785),
+ (0x3ff7a11473eb0180, 0x3cdb5f54408fdb37),
+ (0x3ff82589994cce00, 0x3cf28acf88afab35),
+ (0x3ff8ace5422aa0c0, 0x3cfb5ba7c55a192d),
+ (0x3ff93737b0cdc5c0, 0x3d027a280e1f92a0),
+ (0x3ff9c49182a3f080, 0x3cf01c7c46b071f3),
+ (0x3ffa5503b23e2540, 0x3cfc8b424491caf8),
+ (0x3ffae89f995ad380, 0x3d06af439a68bb99),
+ (0x3ffb7f76f2fb5e40, 0x3cdbaa9ec206ad4f),
+ (0x3ffc199bdd855280, 0x3cfc2220cb12a092),
+ (0x3ffcb720dcef9040, 0x3d048a81e5e8f4a5),
+ (0x3ffd5818dcfba480, 0x3cdc976816bad9b8),
+ (0x3ffdfc97337b9b40, 0x3cfeb968cac39ed3),
+ (0x3ffea4afa2a490c0, 0x3cf9858f73a18f5e),
+ (0x3fff50765b6e4540, 0x3c99d3e12dd8a18b),
+ ],
+ .precision = 53,
+]
+
+const exp32 = {f : flt32+ -> expgen(f, desc32)
+}
+
+const exp64 = {f : flt64+ -> expgen(f, desc64)
+}
+
+generic expgen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i+ var b = d.tobits(f)
+ var n, e, s
+ (n, e, s) = d.explode(f)
+
+ /*
+ Detect if exp(f) would round to outside representability.
+ We don't do bias adjustment, so Tang's thresholds can
+ be tightened.
+ */
+ if !n && b >= d.thresh_1_max
+ -> d.frombits(d.inf)
+ elif n && b > d.thresh_1_min
+ -> (0.0 : @f)
+ ;;
+
+ /* Detect if exp(f) is close to 1. */
+ if (b & ~d.sgnmask) <= d.thresh_2
+ -> (1.0 : @f)
+ ;;
+
+ /* Argument reduction to [ -ln(2)/64, ln(2)/64 ] */
+ var inv_L = d.frombits(d.inv_L)
+
+ var N = rn(f * inv_L)
+ var N2 = N % (32 : @i)
+ if N2 < 0
+ N2 += (32 : @i)
+ ;;
+ var N1 = N - N2
+
+ var R1 : @f, R2 : @f
+ var Nf = (N : @f)
+
+ /*
+ There are few enough significant bits that these are all
+ exact, without need for an FMA. R1 + R2 will approximate
+ (very well) f reduced into [ -ln(2)/64, ln(2)/64 ]
+ */
+ if std.abs(N) >= (1 << d.nabs)
+ R1 = (f - (N1 : @f) * d.frombits(d.L1)) - ((N2 : @f) * d.frombits(d.L1))
+ else
+ R1 = f - (N : @f) * d.frombits(d.L1)
+ ;;
+ R2 = -1.0 * (N : @f) * d.frombits(d.L2)
+
+ var M = (N1 : @i) / 32
+ var J = (N2 : std.size)
+
+ /* Compute a polynomial approximation of exp1m */
+ var R = R1 + R2
+ var Q = R * R * d.horner(R, d.Ai)
+ var P = R1 + (R2 + Q)
+
+ /* Expand out from reduced range */
+ var Su_hi, Su_lo
+ (Su_hi, Su_lo) = d.S[J]
+ var S_hi = d.frombits(Su_hi)
+ var S_lo = d.frombits(Su_lo)
+
+ var S = S_hi + S_lo
+ var unscaled = S_hi + (S_lo + (P * S))
+ var nu, eu, su
+ (nu, eu, su) = d.explode(unscaled)
+ var exp = d.assem(nu, eu + M, su)
+ if (d.tobits(exp) == 0)
+ /* Make sure we don't quite return 0 */
+ -> d.frombits(1)
+ ;;
+
+ -> exp
+}
+
+const expm132 = {f : flt32+ -> expm1gen(f, desc32)
+}
+
+const expm164 = {f : flt64+ -> expm1gen(f, desc64)
+}
+
+generic expm1gen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i+ var b = d.tobits(f)
+ var n, e, s
+ (n, e, s) = d.explode(f)
+
+ /* Special cases: +/- 0, inf, NaN, tiny, and huge */
+ if (b & ~d.sgnmask == 0)
+ -> f
+ elif n && (b & ~d.sgnmask == d.inf)
+ -> (-1.0 : @f)
+ elif (b & ~d.sgnmask == d.inf)
+ -> f
+ elif std.isnan(f)
+ -> d.frombits(d.nan)
+ elif (b & ~d.sgnmask) <= d.thresh_tiny
+ var two_to_large = d.assem(false, 100, 0)
+ var two_to_small = d.assem(false, -100, 0)
+ var abs_f = d.assem(false, e, s)
+ -> (two_to_large * f + abs_f) * two_to_small
+ elif !n && b >= d.thresh_1_max /* exp(x) = oo <=> expm1(x) = oo, as it turns out */
+ -> d.frombits(d.inf)
+ elif n && b >= d.thresh_huge_neg
+ -> (-1.0 : @f)
+ ;;
+
+ if ((n && b < d.Log3_4) || (!n && b < d.Log5_4))
+ /* Procedure 2 */
+
+ /* compute x^2 / 2 with extra precision */
+ var u = round(f, d)
+ var v = f - u
+ var y = u * u * (0.5 : @f)
+ var z = v * (f + u) * (0.5 : @f)
+ var q = f * f * f * d.horner(f, d.Bi)
+
+ var yn, ye, ys
+ (yn, ye, ys) = d.explode(y)
+ if (ye >= -7)
+ -> (u + y) + (q + (v + z))
+ else
+ -> f + (y + (q + z))
+ ;;
+ ;;
+
+ /* Procedure 1 */
+ var inv_L = d.frombits(d.inv_L)
+
+ var N = rn(f * inv_L)
+ var N2 = N % (32 : @i)
+ if N2 < 0
+ N2 += (32 : @i)
+ ;;
+ var N1 = N - N2
+
+ var R1 : @f, R2 : @f
+
+ /*
+ As in the exp case, R1 + R2 very well approximates f
+ reduced into [ -ln(2)/64, ln(2)/64 ]
+ */
+ if std.abs(N) >= (1 << d.nabs)
+ R1 = (f - (N1 : @f) * d.frombits(d.L1)) - ((N2 : @f) * d.frombits(d.L1))
+ else
+ R1 = f - (N : @f) * d.frombits(d.L1)
+ ;;
+ R2 = -1.0 * (N : @f) * d.frombits(d.L2)
+
+ var M = (N1 : @i) / 32
+ var J = (N2 : std.size)
+
+ /* Compute a polynomial approximation of exp1m */
+ var R = R1 + R2
+ var Q = R * R * d.horner(R, d.Ai)
+ var P = R1 + (R2 + Q)
+
+ /* Expand out */
+ var Su_hi, Su_lo
+ (Su_hi, Su_lo) = d.S[J]
+ var S_hi = d.frombits(Su_hi)
+ var S_lo = d.frombits(Su_lo)
+ var S = S_hi + S_lo
+
+ /*
+ Note that [Tan92] includes cases to avoid tripping the
+ underflow flag when not technically required. We currently
+ do not care about such flags, so that logic is omitted.
+ */
+ if M >= d.precision
+ -> scale2(S_hi + (S * P + (S_lo - scale2(1.0, -M))), M)
+ elif M <= -8
+ -> scale2(S_hi + (S * P + S_lo), M) - (1.0 : @f)
+ ;;
+
+ -> scale2((S_hi - scale2(1.0, -M)) + (S_hi * P + S_lo * (P + (1.0 : @f))), M)
+}
+
+generic round = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i+ var b = d.tobits(f)
+ var n, e, s
+ (n, e, s) = d.explode(f)
+ var mask = ~((1 << d.nabs) - 1)
+ if need_round_away(0, ((s & mask) : uint64), (d.nabs : int64) + 1)
+ -> d.assem(n, e, 1 + s & mask)
+ ;;
+ -> d.assem(n, e, s & mask)
+}
--- /dev/null
+++ b/lib/math/fma-impl+posixy-x64-fma.s
@@ -1,0 +1,13 @@
+.globl math$fma32
+.globl _math$fma32
+math$fma32:
+_math$fma32:
+ vfmadd132ss %xmm1, %xmm2, %xmm0
+ ret
+
+.globl math$fma64
+.globl _math$fma64
+math$fma64:
+_math$fma64:
+ vfmadd132sd %xmm1, %xmm2, %xmm0
+ ret
--- /dev/null
+++ b/lib/math/fma-impl.myr
@@ -1,0 +1,446 @@
+use std
+
+use "util"
+
+pkg math =
+ pkglocal const fma32 : (x : flt32, y : flt32, z : flt32 -> flt32)
+ pkglocal const fma64 : (x : flt64, y : flt64, z : flt64 -> flt64)
+;;
+
+const exp_mask32 : uint32 = 0xff << 23
+const exp_mask64 : uint64 = 0x7ff << 52
+
+pkglocal const fma32 = {x : flt32, y : flt32, z : flt32+ var xn, yn
+ (xn, _, _) = std.flt32explode(x)
+ (yn, _, _) = std.flt32explode(y)
+ var xd : flt64 = flt64fromflt32(x)
+ var yd : flt64 = flt64fromflt32(y)
+ var zd : flt64 = flt64fromflt32(z)
+ var prod : flt64 = xd * yd
+ var pn, pe, ps
+ (pn, pe, ps) = std.flt64explode(prod)
+ if pe == -1023
+ pe = -1022
+ ;;
+ if pn != (xn != yn)
+ /* In case of NaNs, sign might not have been preserved */
+ pn = (xn != yn)
+ prod = std.flt64assem(pn, pe, ps)
+ ;;
+
+ var r : flt64 = prod + zd
+ var rn, re, rs
+ (rn, re, rs) = std.flt64explode(r)
+
+ /*
+ At this point, r is probably the correct answer. The
+ only issue is the rounding.
+
+ Ex 1: If x*y > 0 and z is a tiny, negative number, then
+ adding z probably does no rounding. However, if
+ truncating to 23 bits of precision would cause round-to-even,
+ and that round would be upwards, then we need to remember
+ those trailing bits of z and cancel the rounding.
+
+ Ex 2: If x, y, z > 0, and z is small, with
+ last bit in flt64 |
+ last bit in flt32 v v
+ x * y = ...............101011..11
+ z = 10000...,
+ then x * y + z will be rounded to
+ ...............101100..00,
+ and then as a flt32 it will become
+ ...............110,
+ Even though, looking at the original bits, it doesn't
+ "deserve" the final rounding.
+
+ These can only happen if r is non-inf, non-NaN, and the
+ lower 29 bits correspond to "exactly halfway".
+ */
+ if re == 1024 || rs & 0x1fffffff != 0x10000000
+ -> flt32fromflt64(r)
+ ;;
+
+ /*
+ At this point, a rounding is about to happen. We need
+ to know what direction that rounding is, so that we can
+ tell if it's wrong. +1 means "away from 0", -1 means
+ "towards 0".
+ */
+ var zn, ze, zs
+ (zn, ze, zs) = std.flt64explode(zd)
+ var round_direction = 0
+ if rs & 0x20000000 == 0
+ round_direction = -1
+ else
+ round_direction = 1
+ ;;
+
+ var smaller, larger, smaller_e, larger_e
+ if pe > ze || (pe == ze && ps > zs)
+ (smaller, larger, smaller_e, larger_e) = (zs, ps, ze, pe)
+ else
+ (smaller, larger, smaller_e, larger_e) = (ps, zs, pe, ze)
+ ;;
+ var mask = shr((-1 : uint64), 64 - std.min(64, larger_e - smaller_e))
+ var prevent_rounding = false
+ if (round_direction > 0 && pn != zn) || (round_direction < 0 && pn == zn)
+ /*
+ The prospective rounding disagrees with the
+ signage. We are potentially in the case of Ex
+ 1.
+
+ Look at the bits (of the smaller flt64) that are
+ outside the range of r. If there are any such
+ bits, we need to cancel the rounding.
+
+ We certainly need to consider bits very far to
+ the right, but there's an awkwardness concerning
+ the bit just outside the flt64 range: it governed
+ round-to-even, so it might have had an effect.
+ We only care about bits which did not have an
+ effect. Therefore, we perform the subtraction
+ using only the bits from smaller that lie in
+ larger's range, then check whether the result
+ is susceptible to round-to-even.
+
+ (Since we only care about the last bit, and the
+ base is 2, subtraction or addition are equally
+ useful.)
+ */
+ if (larger ^ shr(smaller, larger_e - smaller_e)) & 0x1 == 0
+ prevent_rounding = smaller & mask != 0
+ ;;
+ else
+ /*
+ The prospective rounding agrees with the signage.
+ We are potentially in the case of Ex 2.
+
+ We just need to check if r was obtained by
+ rounding in the addition step. In this case, we
+ still check the smaller/larger, and we only
+ care about round-to-even. Any
+ rounding that happened previously is enough
+ reason to disqualify this next rounding.
+ */
+ prevent_rounding = (larger ^ shr(smaller, larger_e - smaller_e)) & 0x1 != 0
+ ;;
+
+ if prevent_rounding
+ if round_direction > 0
+ rs--
+ else
+ rs++
+ ;;
+ ;;
+
+ -> flt32fromflt64(std.flt64assem(rn, re, rs))
+}
+
+pkglocal const fma64 = {x : flt64, y : flt64, z : flt64+ var xn : bool, yn : bool, zn : bool
+ var xe : int64, ye : int64, ze : int64
+ var xs : uint64, ys : uint64, zs : uint64
+
+ var xb : uint64 = std.flt64bits(x)
+ var yb : uint64 = std.flt64bits(y)
+ var zb : uint64 = std.flt64bits(z)
+
+ /* check for both NaNs and infinities */
+ if xb & exp_mask64 == exp_mask64 || \
+ yb & exp_mask64 == exp_mask64
+ -> x * y + z
+ elif z == 0.0 || z == -0.0 || x * y == 0.0 || x * y == -0.0
+ -> x * y + z
+ elif zb & exp_mask64 == exp_mask64
+ -> z
+ ;;
+
+ (xn, xe, xs) = std.flt64explode(x)
+ (yn, ye, ys) = std.flt64explode(y)
+ (zn, ze, zs) = std.flt64explode(z)
+ if xe == -1023
+ xe = -1022
+ ;;
+ if ye == -1023
+ ye = -1022
+ ;;
+ if ze == -1023
+ ze = -1022
+ ;;
+
+ /* Keep product in high/low uint64s */
+ var xs_h : uint64 = xs >> 32
+ var ys_h : uint64 = ys >> 32
+ var xs_l : uint64 = xs & 0xffffffff
+ var ys_l : uint64 = ys & 0xffffffff
+
+ var t_l : uint64 = xs_l * ys_l
+ var t_m : uint64 = xs_l * ys_h + xs_h * ys_l
+ var t_h : uint64 = xs_h * ys_h
+
+ var prod_l : uint64 = t_l + (t_m << 32)
+ var prod_h : uint64 = t_h + (t_m >> 32)
+ if t_l > prod_l
+ prod_h++
+ ;;
+
+ var prod_n = xn != yn
+ var prod_lastbit_e = (xe - 52) + (ye - 52)
+ var prod_first1 = find_first1_64_hl(prod_h, prod_l, 105)
+ var prod_firstbit_e = prod_lastbit_e + prod_first1
+
+ var z_firstbit_e = ze
+ var z_lastbit_e = ze - 52
+ var z_first1 = 52
+
+ /* subnormals could throw firstbit_e calculations out of whack */
+ if (zb & exp_mask64 == 0)
+ z_first1 = find_first1_64(zs, z_first1)
+ z_firstbit_e = z_lastbit_e + z_first1
+ ;;
+
+ var res_n
+ var res_h = 0
+ var res_l = 0
+ var res_first1
+ var res_lastbit_e
+ var res_firstbit_e
+
+ if prod_n == zn
+ res_n = prod_n
+
+ /*
+ Align prod and z so that the top bit of the
+ result is either 53 or 54, then add.
+ */
+ if prod_firstbit_e >= z_firstbit_e
+ /*
+ [ prod_h ][ prod_l ]
+ [ z...
+ */
+ res_lastbit_e = prod_lastbit_e
+ (res_h, res_l) = (prod_h, prod_l)
+ (res_h, res_l) = add_shifted(res_h, res_l, zs, z_lastbit_e - prod_lastbit_e)
+ else
+ /*
+ [ prod_h ][ prod_l ]
+ [ z...
+ */
+ res_lastbit_e = z_lastbit_e - 64
+ res_h = zs
+ res_l = 0
+ if prod_lastbit_e >= res_lastbit_e + 64
+ /* In this situation, prod must be extremely subnormal */
+ res_h += shl(prod_l, prod_lastbit_e - res_lastbit_e - 64)
+ elif prod_lastbit_e >= res_lastbit_e
+ res_h += shl(prod_h, prod_lastbit_e - res_lastbit_e)
+ res_h += shr(prod_l, res_lastbit_e + 64 - prod_lastbit_e)
+ res_l += shl(prod_l, prod_lastbit_e - res_lastbit_e)
+ elif prod_lastbit_e + 64 >= res_lastbit_e
+ res_h += shr(prod_h, res_lastbit_e - prod_lastbit_e)
+ var l1 = shl(prod_h, prod_lastbit_e + 64 - res_lastbit_e)
+ var l2 = shr(prod_l, res_lastbit_e - prod_lastbit_e)
+ res_l = l1 + l2
+ if res_l < l1
+ res_h++
+ ;;
+ elif prod_lastbit_e + 128 >= res_lastbit_e
+ res_l += shr(prod_h, res_lastbit_e - prod_lastbit_e - 64)
+ ;;
+ ;;
+ else
+ match compare_hl_z(prod_h, prod_l, prod_firstbit_e, prod_lastbit_e, zs, z_firstbit_e, z_lastbit_e)
+ | `std.Equal: -> 0.0
+ | `std.Before:
+ /* prod > z */
+ res_n = prod_n
+ res_lastbit_e = prod_lastbit_e
+ (res_h, res_l) = sub_shifted(prod_h, prod_l, zs, z_lastbit_e - prod_lastbit_e)
+ | `std.After:
+ /* z > prod */
+ res_n = zn
+ res_lastbit_e = z_lastbit_e - 64
+ (res_h, res_l) = sub_shifted(zs, 0, prod_h, prod_lastbit_e + 64 - (z_lastbit_e - 64))
+ (res_h, res_l) = sub_shifted(res_h, res_l, prod_l, prod_lastbit_e - (z_lastbit_e - 64))
+ ;;
+ ;;
+
+ res_first1 = 64 + find_first1_64(res_h, 55)
+ if res_first1 == 63
+ res_first1 = find_first1_64(res_l, 63)
+ ;;
+ res_firstbit_e = res_first1 + res_lastbit_e
+
+ /*
+ Finally, res_h and res_l are the high and low bits of
+ the result. They now need to be assembled into a flt64.
+ Subnormals and infinities could be a problem.
+ */
+ var res_s = 0
+ if res_firstbit_e <= -1023
+ /* Subnormal case */
+ if res_lastbit_e + 128 < 12 - 1022
+ res_s = shr(res_h, 12 - 1022 - (res_lastbit_e + 128))
+ res_s |= shr(res_l, 12 - 1022 - (res_lastbit_e + 64))
+ elif res_lastbit_e + 64 < 12 - 1022
+ res_s = shl(res_h, -12 + (res_lastbit_e + 128) - (-1022))
+ res_s |= shr(res_l, 12 - 1022 - (res_lastbit_e + 64))
+ else
+ res_s = shl(res_h, -12 + (res_lastbit_e + 128) - (-1022))
+ res_s |= shl(res_l, -12 + (res_lastbit_e + 64) - (-1022))
+ ;;
+
+ if need_round_away(res_h, res_l, res_first1 + (-1074 - res_firstbit_e))
+ res_s++
+ ;;
+
+ /* No need for exponents, they are all zero */
+ var res = res_s
+ if res_n
+ res |= (1 << 63)
+ ;;
+ -> std.flt64frombits(res)
+ ;;
+
+ if res_firstbit_e >= 1024
+ /* Infinity case */
+ if res_n
+ -> std.flt64frombits(0xfff0000000000000)
+ else
+ -> std.flt64frombits(0x7ff0000000000000)
+ ;;
+ ;;
+
+ if res_first1 - 52 >= 64
+ res_s = shr(res_h, (res_first1 : int64) - 64 - 52)
+ if need_round_away(res_h, res_l, res_first1 - 52)
+ res_s++
+ ;;
+ elif res_first1 - 52 >= 0
+ res_s = shl(res_h, 64 - (res_first1 - 52))
+ res_s |= shr(res_l, res_first1 - 52)
+ if need_round_away(res_h, res_l, res_first1 - 52)
+ res_s++
+ ;;
+ else
+ res_s = shl(res_h, res_first1 - 52)
+ ;;
+
+ /* The res_s++s might have messed everything up */
+ if res_s & (1 << 53) != 0
+ res_s >= 1
+ res_firstbit_e++
+ if res_firstbit_e >= 1024
+ if res_n
+ -> std.flt64frombits(0xfff0000000000000)
+ else
+ -> std.flt64frombits(0x7ff0000000000000)
+ ;;
+ ;;
+ ;;
+
+ -> std.flt64assem(res_n, res_firstbit_e, res_s)
+}
+
+/*
+ Add (a << s) to [ h ][ l ], where if s < 0 then a corresponding
+ right-shift is used. This is aligned such that if s == 0, then
+ the result is [ h ][ l + a ]
+ */
+const add_shifted = {h : uint64, l : uint64, a : uint64, s : int64+ if s >= 64
+ -> (h + shl(a, s - 64), l)
+ elif s >= 0
+ var new_h = h + shr(a, 64 - s)
+ var sa = shl(a, s)
+ var new_l = l + sa
+ if new_l < l
+ new_h++
+ ;;
+ -> (new_h, new_l)
+ else
+ var new_h = h
+ var sa = shr(a, -s)
+ var new_l = l + sa
+ if new_l < l
+ new_h++
+ ;;
+ -> (new_h, new_l)
+ ;;
+}
+
+/* As above, but subtract (a << s) */
+const sub_shifted = {h : uint64, l : uint64, a : uint64, s : int64+ if s >= 64
+ -> (h - shl(a, s - 64), l)
+ elif s >= 0
+ var new_h = h - shr(a, 64 - s)
+ var sa = shl(a, s)
+ var new_l = l - sa
+ if sa > l
+ new_h--
+ ;;
+ -> (new_h, new_l)
+ else
+ var new_h = h
+ var sa = shr(a, -s)
+ var new_l = l - sa
+ if sa > l
+ new_h--
+ ;;
+ -> (new_h, new_l)
+ ;;
+}
+
+const compare_hl_z = {h : uint64, l : uint64, hl_firstbit_e : int64, hl_lastbit_e : int64, z : uint64, z_firstbit_e : int64, z_lastbit_e : int64+ if hl_firstbit_e > z_firstbit_e
+ -> `std.Before
+ elif hl_firstbit_e < z_firstbit_e
+ -> `std.After
+ ;;
+
+ var h_k : int64 = (hl_firstbit_e - hl_lastbit_e - 64)
+ var z_k : int64 = (z_firstbit_e - z_lastbit_e)
+ while h_k >= 0 && z_k >= 0
+ var h1 = h & shl(1, h_k) != 0
+ var z1 = z & shl(1, z_k) != 0
+ if h1 && !z1
+ -> `std.Before
+ elif !h1 && z1
+ -> `std.After
+ ;;
+ h_k--
+ z_k--
+ ;;
+
+ if z_k < 0
+ if (h & shr((-1 : uint64), 64 - h_k) != 0) || (l != 0)
+ -> `std.Before
+ else
+ -> `std.Equal
+ ;;
+ ;;
+
+ var l_k : int64 = 63
+ while l_k >= 0 && z_k >= 0
+ var l1 = l & shl(1, l_k) != 0
+ var z1 = z & shl(1, z_k) != 0
+ if l1 && !z1
+ -> `std.Before
+ elif !l1 && z1
+ -> `std.After
+ ;;
+ l_k--
+ z_k--
+ ;;
+
+ if (z_k < 0) && (l & shr((-1 : uint64), 64 - l_k) != 0)
+ -> `std.Before
+ elif (l_k < 0) && (z & shr((-1 : uint64), 64 - z_k) != 0)
+ -> `std.After
+ ;;
+
+ -> `std.Equal
+}
--- /dev/null
+++ b/lib/math/fpmath.myr
@@ -1,0 +1,150 @@
+use std
+
+pkg math =
+ trait fpmath @f =
+
+ /* exp-impl */
+ exp : (f : @f -> @f)
+ expm1 : (f : @f -> @f)
+
+ /* fma-impl */
+ fma : (x : @f, y : @f, z : @f -> @f)
+
+ /* poly-impl */
+ horner_poly : (x : @f, a : @f[:] -> @f)
+ horner_polyu : (x : @f, a : @u[:] -> @f)
+
+ /* scale2-impl */
+ scale2 : (f : @f, m : @i -> @f)
+
+ /* sqrt-impl */
+ sqrt : (f : @f -> @f)
+
+ /* sum-impl */
+ kahan_sum : (a : @f[:] -> @f)
+ priest_sum : (a : @f[:] -> @f)
+
+ /* trunc-impl */
+ trunc : (f : @f -> @f)
+ ceil : (f : @f -> @f)
+ floor : (f : @f -> @f)
+ ;;
+
+ trait roundable @f -> @i =
+ /* round-impl */
+ rn : (f : @f -> @i)
+ ;;
+
+ impl std.equatable flt32
+ impl std.equatable flt64
+ impl roundable flt64 -> int64
+ impl roundable flt32 -> int32
+ impl fpmath flt32
+ impl fpmath flt64
+;;
+
+/*
+ We consider two floating-point numbers equal if their bits are
+ equal. This does not treat NaNs specially: two distinct NaNs may
+ compare equal, or they may compare distinct (if they arise from
+ different bit patterns).
+
+ Additionally, +0.0 and -0.0 compare differently.
+ */
+impl std.equatable flt32 =
+ eq = {a : flt32, b : flt32; -> std.flt32bits(a) == std.flt32bits(b)}+;;
+
+impl std.equatable flt64 =
+ eq = {a : flt64, b : flt64; -> std.flt64bits(a) == std.flt64bits(b)}+;;
+
+impl roundable flt32 -> int32 =
+ rn = {f : flt32; -> rn32(f) }+;;
+
+impl roundable flt64 -> int64 =
+ rn = {f : flt64; -> rn64(f) }+;;
+
+impl fpmath flt32 =
+ fma = {x, y, z; -> fma32(x, y, z)}+
+ exp = {f; -> exp32(f)}+ expm1 = {f; -> expm132(f)}+
+ horner_poly = {f, a; -> horner_poly32(f, a)}+ horner_polyu = {f, a; -> horner_polyu32(f, a)}+
+ scale2 = {f, m; -> scale232(f, m)}+
+ sqrt = {f; -> sqrt32(f)}+
+ kahan_sum = {l; -> kahan_sum32(l) }+ priest_sum = {l; -> priest_sum32(l) }+
+ trunc = {f; -> trunc32(f)}+ floor = {f; -> floor32(f)}+ ceil = {f; -> ceil32(f)}+
+;;
+
+impl fpmath flt64 =
+ fma = {x, y, z; -> fma64(x, y, z)}+
+ exp = {f; -> exp64(f)}+ expm1 = {f; -> expm164(f)}+
+ horner_poly = {f, a; -> horner_poly64(f, a)}+ horner_polyu = {f, a; -> horner_polyu64(f, a)}+
+ scale2 = {f, m; -> scale264(f, m)}+
+ sqrt = {f; -> sqrt64(f)}+
+ kahan_sum = {l; -> kahan_sum64(l) }+ priest_sum = {l; -> priest_sum64(l) }+
+ trunc = {f; -> trunc64(f)}+ floor = {f; -> floor64(f)}+ ceil = {f; -> ceil64(f)}+;;
+
+extern const rn32 : (f : flt32 -> int32)
+extern const rn64 : (f : flt64 -> int64)
+
+extern const exp32 : (x : flt32 -> flt32)
+extern const exp64 : (x : flt64 -> flt64)
+
+extern const expm132 : (x : flt32 -> flt32)
+extern const expm164 : (x : flt64 -> flt64)
+
+extern const fma32 : (x : flt32, y : flt32, z : flt32 -> flt32)
+extern const fma64 : (x : flt64, y : flt64, z : flt64 -> flt64)
+
+extern const horner_poly32 : (f : flt32, a : flt32[:] -> flt32)
+extern const horner_poly64 : (f : flt64, a : flt64[:] -> flt64)
+
+extern const horner_polyu32 : (f : flt32, a : uint32[:] -> flt32)
+extern const horner_polyu64 : (f : flt64, a : uint64[:] -> flt64)
+
+extern const scale232 : (f : flt32, m : int32 -> flt32)
+extern const scale264 : (f : flt64, m : int64 -> flt64)
+
+extern const sqrt32 : (x : flt32 -> flt32)
+extern const sqrt64 : (x : flt64 -> flt64)
+
+extern const kahan_sum32 : (l : flt32[:] -> flt32)
+extern const kahan_sum64 : (l : flt64[:] -> flt64)
+
+extern const priest_sum32 : (l : flt32[:] -> flt32)
+extern const priest_sum64 : (l : flt64[:] -> flt64)
+
+extern const trunc32 : (f : flt32 -> flt32)
+extern const trunc64 : (f : flt64 -> flt64)
+
+extern const floor32 : (f : flt32 -> flt32)
+extern const floor64 : (f : flt64 -> flt64)
+
+extern const ceil32 : (f : flt32 -> flt32)
+extern const ceil64 : (f : flt64 -> flt64)
--- /dev/null
+++ b/lib/math/poly-impl.myr
@@ -1,0 +1,47 @@
+use std
+
+/* See [Mul16], section 5.1 */
+pkg math =
+ pkglocal const horner_poly32 : (f : flt32, a : flt32[:] -> flt32)
+ pkglocal const horner_poly64 : (f : flt64, a : flt64[:] -> flt64)
+
+ pkglocal const horner_polyu32 : (f : flt32, a : uint32[:] -> flt32)
+ pkglocal const horner_polyu64 : (f : flt64, a : uint64[:] -> flt64)
+;;
+
+extern const fma32 : (x : flt32, y : flt32, z : flt32 -> flt32)
+extern const fma64 : (x : flt64, y : flt64, z : flt64 -> flt64)
+
+const horner_poly32 = {f : flt32, a : flt32[:]+ var r : flt32 = 0.0
+ for var j = a.len - 1; j >= 0; j--
+ r = fma32(r, f, a[j])
+ ;;
+ -> r
+}
+
+const horner_poly64 = {f : flt64, a : flt64[:]+ var r : flt64 = 0.0
+ for var j = a.len - 1; j >= 0; j--
+ r = fma64(r, f, a[j])
+ ;;
+ -> r
+}
+
+const horner_polyu32 = {f : flt32, a : uint32[:]+ var r : flt32 = 0.0
+ for var j = a.len - 1; j >= 0; j--
+ r = fma32(r, f, std.flt32frombits(a[j]))
+ ;;
+ -> r
+}
+
+const horner_polyu64 = {f : flt64, a : uint64[:]+ var r : flt64 = 0.0
+ for var j = a.len - 1; j >= 0; j--
+ r = fma64(r, f, std.flt64frombits(a[j]))
+ ;;
+ -> r
+}
+
+/* TODO: add Knuth's method */
--- /dev/null
+++ b/lib/math/references
@@ -1,0 +1,27 @@
+References
+
+[KM06]
+Peter Kornerup and Jean-Michel Muller. “Choosing starting values
+for certain Newton–Raphson iterations”. In: Theoretical Computer
+Science 351 (1 2006), pp. 101–110. doi:
+https://doi.org/10.1016/j.tcs.2005.09.056.
+
+[Mul+10]
+Jean-Michel Muller et al. Handbook of floating-point arithmetic.
+Boston: Birkhäuser, 2010. isbn: 9780817647049.
+
+[Mul16]
+J. M. Muller. Elementary functions : algorithms and implementation.
+Third edition. New York: Birkhäuser, 2016. isbn: 9781489979810.
+
+[Tan89]
+Ping-Tak Peter Tang. “Table-driven Implementation of the Exponential
+Function in IEEE Floating-point Arithmetic”. In: ACM Trans. Math.
+Softw. 15.2 (June 1989), pp. 144–157. issn: 0098-3500. doi:
+10.1145/63522.214389. url: http://doi.acm.org/10.1145/63522.214389.
+
+[Tan92]
+Ping Tak Peter Tang. “Table-driven Implementation of the Expm1
+Function in IEEE Floating- point Arithmetic”. In: ACM Trans. Math.
+Softw. 18.2 (June 1992), pp. 211–222. issn: 0098-3500. doi:
+10.1145/146847.146928. url: http://doi.acm.org/10.1145/146847.146928.
--- /dev/null
+++ b/lib/math/round-impl+posixy-x64-sse4.s
@@ -1,0 +1,32 @@
+.globl math$rn32
+.globl _math$rn32
+math$rn32:
+_math$rn32:
+ pushq %rbp
+ movq %rsp, %rbp
+ subq $16, %rsp
+ fnstcw (%rsp)
+ fnstcw 8(%rsp)
+ andl $0xf3ff, 8(%rsp)
+ fldcw 8(%rsp)
+ cvtss2si %xmm0, %eax
+ fldcw (%rsp)
+ leave
+ ret
+
+.globl math$rn64
+.globl _math$rn64
+math$rn64:
+_math$rn64:
+ pushq %rbp
+ movq %rsp, %rbp
+ subq $16, %rsp
+ fnstcw (%rsp)
+ fnstcw 8(%rsp)
+ andl $0xf3ff, 8(%rsp)
+ fldcw 8(%rsp)
+ cvtsd2si %xmm0, %rax
+ fldcw (%rsp)
+ leave
+ ret
+
--- /dev/null
+++ b/lib/math/round-impl.myr
@@ -1,0 +1,70 @@
+use std
+
+use "util"
+
+pkg math =
+ const rn64 : (f : flt64 -> int64)
+ const rn32 : (f : flt32 -> int32)
+;;
+
+const rn64 = {f : flt64+ var n : bool, e : int64, s : uint64
+
+ (n, e, s) = std.flt64explode(f)
+
+ if e >= 63
+ -> -9223372036854775808
+ elif e >= 52
+ var shifted : int64 = (( s << (e - 52 : uint64)) : int64)
+ if n
+ -> -1 * shifted
+ else
+ -> shifted
+ ;;
+ elif e < -1
+ -> 0
+ ;;
+
+ var integral_s = (s >> (52 - e : uint64) : int64)
+
+ if need_round_away(0, s, 52 - e)
+ integral_s++
+ ;;
+
+ if n
+ -> -integral_s
+ else
+ -> integral_s
+ ;;
+}
+
+const rn32 = {f : flt32+ var n : bool, e : int32, s : uint32
+
+ (n, e, s) = std.flt32explode(f)
+
+ if e >= 31
+ -> -2147483648
+ elif e >= 23
+ var shifted : int32 = (( s << (e - 23 : uint32)) : int32)
+ if n
+ -> -1 * shifted
+ else
+ -> shifted
+ ;;
+ elif e < -1
+ -> 0
+ ;;
+
+ var integral_s = (s >> (23 - e : uint32) : int32)
+
+ if need_round_away(0, (s : uint64), (23 - e : int64))
+ integral_s++
+ ;;
+
+ if n
+ -> -integral_s
+ else
+ -> integral_s
+ ;;
+}
--- /dev/null
+++ b/lib/math/scale2-impl.myr
@@ -1,0 +1,73 @@
+use std
+
+use "fpmath"
+use "util"
+
+/*
+ The scaleB function recommended by the 2008 revision of IEEE
+ 754. Since only integer exponents are considered, the naive
+ approach works quite well.
+ */
+pkg math =
+ const scale232 : (f : flt32, m : int32 -> flt32)
+ const scale264 : (f : flt64, m : int64 -> flt64)
+;;
+
+const scale232 = {f : flt32, m : int32+ var n, e, s
+ (n, e, s) = std.flt32explode(f)
+ (n, e, s) = scale2gen(n, e, s, -126, 127, 24, m)
+ -> std.flt32assem(n, e, s)
+}
+
+const scale264 = {f : flt64, m : int64+ var n, e, s
+ (n, e, s) = std.flt64explode(f)
+ (n, e, s) = scale2gen(n, e, s, -1022, 1023, 53, m)
+ -> std.flt64assem(n, e, s)
+}
+
+generic scale2gen = {n : bool, e : @i, s : @u, emin : @i, emax : @i, p : @u, m : @i :: numeric, integral @i, numeric, integral @u+ if e < emin && s == 0
+ -> (n, e, s)
+ elif m == 0
+ -> (n, e, s)
+ elif m < 0
+ var sprime = s
+ var eprime = e
+ if e < emin
+ sprime = s >> (-m)
+ if need_round_away(0, (s : uint64), (-m : int64))
+ sprime++
+ ;;
+ eprime = emin - 1
+ elif e + m < emin
+ sprime = s >> (emin - m - e)
+ if need_round_away(0, (s : uint64), ((emin - m - e) : int64))
+ sprime++
+ ;;
+ eprime = emin - 1
+ else
+ eprime = e + m
+ ;;
+ -> (n, eprime, sprime)
+ ;;
+
+ if e < emin
+ var first_1 = find_first1_64((s : uint64), (p : int64))
+ var shift = p - (first_1 : @u) - 1
+ if m >= (shift : @i)
+ s = s << shift
+ m -= (shift : @i)
+ e = emin
+ else
+ -> (n, e, s << m)
+ ;;
+ ;;
+
+ if e + m > emax && s != 0
+ -> (n, emax + 1, 1 << (p - 1))
+ ;;
+
+ -> (n, e + m, s)
+}
--- /dev/null
+++ b/lib/math/sqrt-impl+posixy-x64-sse2.s
@@ -1,0 +1,13 @@
+.globl math$sqrt32
+.globl _math$sqrt32
+math$sqrt32:
+_math$sqrt32:
+ sqrtss %xmm0, %xmm0
+ ret
+
+.globl math$sqrt64
+.globl _math$sqrt64
+math$sqrt64:
+_math$sqrt64:
+ sqrtsd %xmm0, %xmm0
+ ret
--- /dev/null
+++ b/lib/math/sqrt-impl.myr
@@ -1,0 +1,189 @@
+use std
+
+use "fpmath"
+
+/* See [Mul+10], sections 5.4 and 8.7 */
+pkg math =
+ pkglocal const sqrt32 : (f : flt32 -> flt32)
+ pkglocal const sqrt64 : (f : flt64 -> flt64)
+;;
+
+extern const fma32 : (x : flt32, y : flt32, z : flt32 -> flt32)
+extern const fma64 : (x : flt64, y : flt64, z : flt64 -> flt64)
+
+type fltdesc(@f, @u, @i) = struct
+ explode : (f : @f -> (bool, @i, @u))
+ assem : (n : bool, e : @i, s : @u -> @f)
+ fma : (x : @f, y : @f, z : @f -> @f)
+ tobits : (f : @f -> @u)
+ frombits : (u : @u -> @f)
+ nan : @u
+ emin : @i
+ emax : @i
+ normmask : @u
+ sgnmask : @u
+ ab : (@u, @u)[:]
+ iterlim : int
+;;
+
+/*
+ The starting point of the N-R iteration of 1/sqrt, after significand
+ has been normalized to [1, 4).
+
+ See [KM06] for the construction and notation. Case p = -2. The
+ dividers (left values) are chosen roughly so that maximal error
+ of N-R, after 3 iterations, starting with the right value, is
+ less than an ulp (of the result). If g falls in [a_i, a_{i+1}),+ N-R should start with b_{i+1}.+
+ In the flt64 case, we need only one more iteration.
+ */
+const ab32 : (uint32, uint32)[:] = [
+ (0x3f800000, 0x3f800000), /* Nothing should ever get normalized to < 1.0 */
+ (0x3fa66666, 0x3f6f30ae), /* [1.0, 1.3 ) -> 0.9343365431 */
+ (0x3fd9999a, 0x3f5173ca), /* [1.3, 1.7 ) -> 0.8181730509 */
+ (0x40100000, 0x3f3691d3), /* [1.7, 2.25) -> 0.713162601 */
+ (0x40333333, 0x3f215342), /* [2.25, 2.8 ) -> 0.6301766634 */
+ (0x4059999a, 0x3f118e0e), /* [2.8, 3.4 ) -> 0.5685738325 */
+ (0x40800000, 0x3f053049), /* [3.4, 4.0 ) -> 0.520268023 */
+][:]
+
+const ab64 : (uint64, uint64)[:] = [
+ (0x3ff0000000000000, 0x3ff0000000000000), /* < 1.0 */
+ (0x3ff3333333333333, 0x3fee892ce1608cbc), /* [1.0, 1.2) -> 0.954245033445111356940060431953 */
+ (0x3ff6666666666666, 0x3fec1513a2184094), /* [1.2, 1.4) -> 0.877572838393478438234751592972 */
+ (0x3ffc000000000000, 0x3fe9878eb3e9ba20), /* [1.4, 1.75) -> 0.797797538178034670863780775107 */
+ (0x400199999999999a, 0x3fe6ccb14eeb238d), /* [1.75, 2.2) -> 0.712486890924184046447464879748 */
+ (0x400599999999999a, 0x3fe47717c17cd34f), /* [2.2, 2.7) -> 0.639537694840876969060161627567 */
+ (0x400b333333333333, 0x3fe258df212a8e9a), /* [2.7, 3.4) -> 0.573348583963212421465982515656 */
+ (0x4010000000000000, 0x3fe0a5989f2dc59a), /* [3.4, 4.0) -> 0.520214377304159869552790951275 */
+][:]
+
+const sqrt32 = {f : flt32+ const d : fltdesc(flt32, uint32, int32) = [
+ .explode = std.flt32explode,
+ .assem = std.flt32assem,
+ .fma = fma32,
+ .tobits = std.flt32bits,
+ .frombits = std.flt32frombits,
+ .nan = 0x7fc00000,
+ .emin = -127,
+ .emax = 128,
+ .normmask = 1 << 23,
+ .sgnmask = 1 << 31,
+ .ab = ab32,
+ .iterlim = 3,
+ ]
+ -> sqrtgen(f, d)
+}
+
+const sqrt64 = {f : flt64+ const d : fltdesc(flt64, uint64, int64) = [
+ .explode = std.flt64explode,
+ .assem = std.flt64assem,
+ .fma = fma64,
+ .tobits = std.flt64bits,
+ .frombits = std.flt64frombits,
+ .nan = 0x7ff8000000000000,
+ .emin = -1023,
+ .emax = 1024,
+ .normmask = 1 << 52,
+ .sgnmask = 1 << 63,
+ .ab = ab64,
+ .iterlim = 4,
+ ]
+ -> sqrtgen(f, d)
+}
+
+generic sqrtgen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i+ var n : bool, e : @i, s : @u, e2 : @i
+ (n, e, s) = d.explode(f)
+
+ /* Special cases: +/- 0.0, negative, NaN, and +inf */
+ if e == d.emin && s == 0
+ -> f
+ elif n || std.isnan(f)
+ /* Make sure to return a quiet NaN */
+ -> d.frombits(d.nan)
+ elif e == d.emax
+ -> f
+ ;;
+
+ /*
+ Remove a factor of 2^(even) to normalize significand.
+ */
+ if e == d.emin
+ e = d.emin + 1
+ while s & d.normmask == 0
+ s <<= 1
+ e--
+ ;;
+ ;;
+ if e % 2 != 0
+ e2 = e - 1
+ e = 1
+ else
+ e2 = e
+ e = 0
+ ;;
+
+ var a : @f = d.assem(false, e, s)
+ var au : @u = d.tobits(a)
+
+ /*
+ We shall perform iterated Newton-Raphson in order to
+ compute 1/sqrt(g), then multiply by g to obtain sqrt(g).
+ This is faster than calculating sqrt(g) directly because
+ it avoids division. (The multiplication by g is built
+ into Markstein's r, g, n variables.)
+ */
+ var xn : @f = d.frombits(0)
+ for (ai, beta) : d.ab
+ if au <= ai
+ xn = d.frombits(beta)
+ break
+ ;;
+ ;;
+
+ /* split up "x_{n+1} = x_n (3 - ax_n^2)/2" */+ var epsn = d.fma(-1.0 * a, xn * xn, 1.0)
+ var rn = 0.5 * epsn
+ var gn = a * xn
+ var hn = 0.5 * xn
+ for var j = 0; j < d.iterlim; ++j
+ rn = d.fma(-1.0 * gn, hn, 0.5)
+ gn = d.fma(gn, rn, gn)
+ hn = d.fma(hn, rn, hn)
+ ;;
+
+ /*
+ gn is almost what we want, except that we might want to
+ adjust by an ulp in one direction or the other. This is
+ the Tuckerman test.
+
+ Exhaustive testing has shown that we need only 3 adjustments
+ in the flt32 case (and it should be 4 in the flt64 case).
+ */
+ (_, e, s) = d.explode(gn)
+ e += (e2 / 2)
+ var r : @f = d.assem(false, e, s)
+
+ for var j = 0; j < d.iterlim; ++j
+ var r_plus_ulp : @f = d.frombits(d.tobits(r) + 1)
+ var r_minus_ulp : @f = d.frombits(d.tobits(r) - 1)
+
+ var delta_1 = d.fma(r, r_minus_ulp, -1.0 * f)
+ if d.tobits(delta_1) & d.sgnmask == 0
+ r = r_minus_ulp
+ else
+ var delta_2 = d.fma(r, r_plus_ulp, -1.0 * f)
+ if d.tobits(delta_2) & d.sgnmask != 0
+ r = r_plus_ulp
+ else
+ -> r
+ ;;
+ ;;
+ ;;
+
+ -> r
+}
--- /dev/null
+++ b/lib/math/sum-impl.myr
@@ -1,0 +1,105 @@
+use std
+
+/* For references, see [Mul+10] section 6.3 */
+pkg math =
+ pkglocal const kahan_sum32 : (l : flt32[:] -> flt32)
+ pkglocal const priest_sum32 : (l : flt32[:] -> flt32)
+
+ pkglocal const kahan_sum64: (l : flt64[:] -> flt64)
+ pkglocal const priest_sum64 : (l : flt64[:] -> flt64)
+;;
+
+type doomed_flt32_arr = flt32[:]
+type doomed_flt64_arr = flt64[:]
+
+impl disposable doomed_flt32_arr =
+ __dispose__ = {a : doomed_flt32_arr; std.slfree((a : flt32[:])) }+;;
+
+impl disposable doomed_flt64_arr =
+ __dispose__ = {a : doomed_flt64_arr; std.slfree((a : flt64[:])) }+;;
+
+/*
+ Kahan's compensated summation. Fast and reasonably accurate,
+ although cancellation can cause relative error blowup. For
+ something slower, but more accurate, use something like Priest's
+ doubly compensated sums.
+ */
+pkglocal const kahan_sum32 = {l; -> kahan_sum_gen(l, (0.0 : flt32))}+pkglocal const kahan_sum64 = {l; -> kahan_sum_gen(l, (0.0 : flt64))}+
+generic kahan_sum_gen = {l : @f[:], zero : @f :: numeric,floating @f+ if l.len == 0
+ -> zero
+ ;;
+
+ var s = zero
+ var c = zero
+ var y = zero
+ var t = zero
+
+ for x : l
+ y = x - c
+ t = s + y
+ c = (t - s) - y
+ s = t
+ ;;
+
+ -> s
+}
+
+/*
+ Priest's doubly compensated summation. Extremely accurate, but
+ relatively slow. For situations in which cancellation is not
+ expected, something like Kahan's compensated summation may be
+ more useful.
+ */
+pkglocal const priest_sum32 = {l : flt32[:]+ var l2 = std.sldup(l)
+ std.sort(l2, mag_cmp32)
+ auto (l2 : doomed_flt32_arr)
+ -> priest_sum_gen(l2, (0.0 : flt32))
+}
+
+const mag_cmp32 = {f : flt32, g : flt32+ var u = std.flt32bits(f) & ~(1 << 31)
+ var v = std.flt32bits(g) & ~(1 << 31)
+ -> std.numcmp(v, u)
+}
+
+pkglocal const priest_sum64 = {l : flt64[:]+ var l2 = std.sldup(l)
+ std.sort(l, mag_cmp64)
+ auto (l2 : doomed_flt64_arr)
+ -> priest_sum_gen(l2, (0.0 : flt64))
+}
+
+const mag_cmp64 = {f : flt64, g : flt64+ var u = std.flt64bits(f) & ~(1 << 63)
+ var v = std.flt64bits(g) & ~(1 << 63)
+ -> std.numcmp(v, u)
+}
+
+generic priest_sum_gen = {l : @f[:], zero : @f :: numeric,floating @f+ /* l should be sorted in descending order */
+ if l.len == 0
+ -> zero
+ ;;
+
+ var s = zero
+ var c = zero
+
+ for x : l
+ var y = c + x
+ var u = x - (y - c)
+ var t = (y + s)
+ var v = (y - (t - s))
+ var z = u + v
+ s = t + z
+ c = z - (s - t)
+ ;;
+
+ -> s
+}
+
--- /dev/null
+++ b/lib/math/test/exp-impl.myr
@@ -1,0 +1,285 @@
+use std
+use math
+use testr
+
+/*
+ Note: a major part of the algorithms are the S constants. They
+ are tested extensively in expm101 and expm102.
+ */
+const main = {+ testr.run([
+ [.name="exp-01", .fn = exp01],
+ [.name="exp-02", .fn = exp02],
+ [.name="exp-03", .fn = exp03],
+ [.name="exp-04", .fn = exp04],
+ [.name="expm1-01", .fn = expm101],
+ [.name="expm1-02", .fn = expm102],
+ [.name="expm1-03", .fn = expm103],
+ [.name="expm1-04", .fn = expm104],
+ ][:])
+}
+
+const exp01 = {c+ var inputs : (uint32, uint32)[:] = [
+ (0x00000000, 0x3f800000),
+ (0x34000000, 0x3f800001),
+ (0x3c000000, 0x3f810101),
+ (0x42000000, 0x568fa1fe),
+ (0xc2b00000, 0x0041edc4),
+ (0xc2b20000, 0x001840fc),
+ (0x7f7fffff, 0x7f800000),
+ (0x7f800000, 0x7f800000),
+ (0x7f800001, 0x7fc00000),
+ (0xc2cff1b3, 0x00000001),
+ (0xc2cff1b4, 0x00000001),
+ (0xc2cff1b5, 0000000000),
+ (0x42b17217, 0x7f7fff84),
+ (0x42b17218, 0x7f800000),
+ (0x42b17219, 0x7f800000),
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var yf : flt32 = std.flt32frombits(y)
+ var rf = math.exp(xf)
+ testr.check(c, rf == yf,
+ "exp(0x{b=16,w=8,p=0}) should be 0x{b=16,w=8,p=0}, was 0x{b=16,w=8,p=0}",+ x, y, std.flt32bits(rf))
+ ;;
+}
+
+const exp02 = {c+ var inputs : (uint64, uint64)[:] = [
+ (0x0000000000000000, 0x3ff0000000000000),
+ (0x3e50000000000000, 0x3ff0000004000000),
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var yf : flt64 = std.flt64frombits(y)
+ var rf = math.exp(xf)
+ testr.check(c, rf == yf,
+ "exp(0x{b=16,w=16,p=0}) should be 0x{b=16,w=16,p=0}, was 0x{b=16,w=16,p=0}",+ x, y, std.flt64bits(rf))
+ ;;
+}
+
+const exp03 = {c+ /*
+ Tang's algorithm has an error of up to 0.77 ulps. This
+ is not terrible (musl appears to follow it, for example).
+ Here we quarantine off some known-bad results.
+ */
+
+ var inputs : (uint32, uint32, uint32)[:] = [
+ (0x42020000, 0x56eccf79, 0x56eccf78),
+ (0x3ec40600, 0x3fbbb54b, 0x3fbbb54c),
+ ][:]
+
+ for (x, y_perfect, y_acceptable) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var ypf : flt32 = std.flt32frombits(y_perfect)
+ var yaf : flt32 = std.flt32frombits(y_acceptable)
+ var rf = math.exp(xf)
+ if rf != ypf && rf != yaf
+ testr.fail(c, "exp(0x{b=16,w=8,p=0}) was 0x{b=16,w=8,p=0}. It should have been 0x{b=16,w=8,p=0}, although we will also accept 0x{b=16,w=8,p=0}",+ x, std.flt32bits(rf), y_perfect, y_acceptable)
+ ;;
+ ;;
+}
+
+const exp04 = {c+ /*
+ Tang's algorithm has an error of up to 0.77 ulps. This
+ is not terrible (musl appears to follow it, for example).
+ Here we quarantine off some known-bad results.
+ */
+
+ var inputs : (uint64, uint64, uint64)[:] = [
+ (0x3cda000000000000, 0x3ff0000000000006, 0x3ff0000000000007),
+ (0x3d57020000000000, 0x3ff00000000005c0, 0x3ff00000000005c1),
+ (0x3d58020000000000, 0x3ff0000000000600, 0x3ff0000000000601),
+ (0xc087030000000000, 0x0000000000000c6d, 0x0000000000000c6e),
+ (0xc011070000000000, 0x3f8d039e34c59187, 0x3f8d039e34c59186),
+ (0xbd50070000000000, 0x3feffffffffff7fc, 0x3feffffffffff7fd),
+ (0xbd430e0000000000, 0x3feffffffffffb3c, 0x3feffffffffffb3d),
+ ][:]
+
+ for (x, y_perfect, y_acceptable) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var ypf : flt64 = std.flt64frombits(y_perfect)
+ var yaf : flt64 = std.flt64frombits(y_acceptable)
+ var rf = math.exp(xf)
+ if rf != ypf && rf != yaf
+ testr.fail(c, "exp(0x{b=16,w=16,p=0}) was 0x{b=16,w=16,p=0}. It should have been 0x{b=16,w=16,p=0}, although we will also accept 0x{b=16,w=16,p=0}",+ x, std.flt64bits(rf), y_perfect, y_acceptable)
+ ;;
+ ;;
+}
+
+const expm101 = {c+ var inputs : (uint32, uint32)[:] = [
+ (0x00000000, 0x00000000),
+ (0x80000000, 0x80000000),
+ (0x3f000000, 0x3f261299),
+ (0x3c000000, 0x3c008056),
+ (0x42000000, 0x568fa1fe),
+ (0xc2b00000, 0xbf800000),
+ (0xc2b20000, 0xbf800000),
+ (0x01000000, 0x01000000),
+ (0x40000000, 0x40cc7326),
+ (0x42b17200, 0x7f7ff404),
+ (0x415a3cf2, 0x494cd0e3),
+ (0x7f800000, 0x7f800000),
+ (0xff800000, 0xbf800000),
+ (0x7a2028b1, 0x7f800000),
+ (0xa201a23a, 0xa201a23a),
+ (0xc0000000, 0xbf5d5aab),
+ (0xbe934b10, 0xbe7fffff),
+ (0xbe934b11, 0xbe800000),
+ (0xbe934b12, 0xbe800001),
+ (0x3e647fbe, 0x3e800000),
+ (0x3e647fbf, 0x3e800000),
+ (0x3e647fc0, 0x3e800001),
+ (0xc0f744f5, 0xbf7fe31e),
+ (0x4210297a, 0x597f31f5), /* J = 0 */
+ (0x3f34c3cd, 0x3f83573d), /* J = 1 */
+ (0x3f3a52b6, 0x3f89087b), /* J = 2 */
+ (0xbf20e72b, 0xbeeee940), /* ... */
+ (0x41f4bd2a, 0x558c999f),
+ (0xc02a0418, 0xbf6e07cd),
+ (0xc0293a2a, 0xbf6dcec1),
+ (0x40b62779, 0x4393ca4b),
+ (0x3fc680ac, 0x406dc6a4),
+ (0x3fc9d2c6, 0x4075b516),
+ (0xbfedd645, 0xbf581273),
+ (0x3e70e5d1, 0x3e87cbdb),
+ (0xbeddcacc, 0xbeb3ffd9),
+ (0x3e8beb21, 0x3ea0e776),
+ (0x3e9ded31, 0x3eb8fe1f),
+ (0x40e8503c, 0x44b19ed8),
+ (0x40d265cb, 0x4432f91f),
+ (0xbea0f9bc, 0xbe8a2036),
+ (0x3ec42672, 0x3eef04c4),
+ (0xc140e8cc, 0xbf7fff9f),
+ (0x4117320e, 0x46467e73),
+ (0x3ee8b75d, 0x3f134ef0),
+ (0xc03f51f1, 0xbf731e4e),
+ (0x42733615, 0x6b52d3d4),
+ (0x3f02f2b5, 0x3f2af617),
+ (0xbf5ac925, 0xbf131660),
+ (0x40813277, 0x425eb7ac),
+ (0x41842e94, 0x4b64b1dd),
+ (0x41b0ba81, 0x4f6a0cc1),
+ (0xc061d7c2, 0xbf787d28),
+ (0xc0611682, 0xbf786657),
+ (0x40dcd7e9, 0x447827c5), /* J = 31 */
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var yf : flt32 = std.flt32frombits(y)
+ var rf = math.expm1(xf)
+ testr.check(c, rf == yf,
+ "expm1(0x{b=16,w=8,p=0}) should be 0x{b=16,w=8,p=0}, was 0x{b=16,w=8,p=0}",+ x, y, std.flt32bits(rf))
+ ;;
+}
+
+const expm102 = {c+ var inputs : (uint64, uint64)[:] = [
+ (0x0000000000000000, 0x0000000000000000),
+ (0x404ef04831cb65ed, 0x45834ac37c44b3d3),
+ (0x7ff0000000000000, 0x7ff0000000000000),
+ (0xfff0000000000000, 0xbff0000000000000),
+ (0x80318a89f290021a, 0x80318a89f290021a),
+ (0xc0180881a9e73af6, 0xbfefebdcaf24d5fe),
+ (0xc020cedaedb028c9, 0xbfeffe2a4ee5ba79),
+ (0xbfe62812ff80cb9e, 0xbfdff9cf6758cc6a), /* J = 0 */
+ (0xbfe526dab7e5054c, 0xbfdef44fe4876d02), /* J = 1 */
+ (0x3ff6ea5c51cbf0f2, 0x400980f836e2c055), /* J = 2 */
+ (0x403f4315360b2cdc, 0x42c12ad5f692ffff), /* ... */
+ (0x3fe93f778a9dc013, 0x3ff3381168aca01e),
+ (0x405023f373e01862, 0x45c1aa771de3ba9c),
+ (0xbfff56366f3b92de, 0xbfeb7c693f8bdbb8),
+ (0xc0159bcd07244a34, 0xbfefdb1461286124),
+ (0x40079bbc23a67f90, 0x40322039bbbb5e2e),
+ (0xbffe1933d2c0ed70, 0xbfeb1f6c166582bb),
+ (0xc009ea7d09d84479, 0xbfeebf01fa92741c),
+ (0x405d84f7563a2612, 0x4a946476fb27817e),
+ (0x4045f770a45f8e7f, 0x43e4da111dae0dfb),
+ (0x4070e9a83b352180, 0x585516d4e37bd9f6),
+ (0x3fefcdb2a528c6e8, 0x3ffb39ec926d8a50),
+ (0x3fd52b5e7995f00c, 0x3fd91739183464e0),
+ (0xbfd625845bbaf98f, 0xbfd2b893d222b8f2),
+ (0xbfd43c0407d114d5, 0xbfd15909a4bea38e),
+ (0x3fd93f8288681821, 0x3fdef406a9f7c4b1),
+ (0x3fdaad696fbfea76, 0x3fe08c8035acbd0b),
+ (0x405e9b6f4b4bc7f3, 0x4af8b6ad079946a7),
+ (0x3fdc85cf6b85bfaf, 0x3fe1f811193e5cf5),
+ (0xbfed55f9b317d7eb, 0xbfe334b0378703a0),
+ (0x406ee7f396a46c9b, 0x563a11333c75a10f),
+ (0x4031266dc880810e, 0x417ac458c1525e64),
+ (0xc019905c018b6d96, 0xbfeff243dffe774d),
+ (0x404dc17e83d7ee6b, 0x454cfbf6572d627f),
+ (0xc013dcd6fae06405, 0xbfefc6dfe34e5408),
+ (0x402762d99fa5cfda, 0x40fd3b9a2004f68b),
+ (0xbfe89074f132e353, 0xbfe12602fb3c9806),
+ (0x4077c7ea24627ae7, 0x623ea61f88281bd5),
+ (0xbff6a873237d8072, 0xbfe83c311800b6ee), /* J = 31 */
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var yf : flt64 = std.flt64frombits(y)
+ var rf = math.expm1(xf)
+ testr.check(c, rf == yf,
+ "expm1(0x{b=16,w=16,p=0}) should be 0x{b=16,w=16,p=0}, was 0x{b=16,w=16,p=0}",+ x, y, std.flt64bits(rf))
+ ;;
+}
+
+const expm103 = {c+ /*
+ As with exp, there is some accepted error in expm1.
+ */
+
+ var inputs : (uint32, uint32, uint32)[:] = [
+ (0x34000000, 0x34000001, 0x34000000),
+ (0xbe651dea, 0xbe4d4b4d, 0xbe4d4b4c),
+ ][:]
+
+ for (x, y_perfect, y_acceptable) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var ypf : flt32 = std.flt32frombits(y_perfect)
+ var yaf : flt32 = std.flt32frombits(y_acceptable)
+ var rf = math.expm1(xf)
+ if rf != ypf && rf != yaf
+ testr.fail(c, "expm1(0x{b=16,w=8,p=0}) was 0x{b=16,w=8,p=0}. It should have been 0x{b=16,w=8,p=0}, although we will also accept 0x{b=16,w=8,p=0}",+ x, std.flt32bits(rf), y_perfect, y_acceptable)
+ ;;
+ ;;
+}
+
+const expm104 = {c+ /*
+ As with exp, there is some accepted error in expm1.
+ */
+
+ var inputs : (uint64, uint64, uint64)[:] = [
+ (0xbf9d0b5aadc4d0ac, 0xbf9ca2e5b7bfa859, 0xbf9ca2e5b7bfa85a),
+ (0x3fc2dbb101fe0392, 0x3fc451731cc0e358, 0x3fc451731cc0e359),
+ (0x3fc8a39bc9c32fec, 0x3fcb2b988c3e0b2f, 0x3fcb2b988c3e0b30),
+ ][:]
+
+ for (x, y_perfect, y_acceptable) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var ypf : flt64 = std.flt64frombits(y_perfect)
+ var yaf : flt64 = std.flt64frombits(y_acceptable)
+ var rf = math.expm1(xf)
+ if rf != ypf && rf != yaf
+ testr.fail(c, "expm1(0x{b=16,w=16,p=0}) was 0x{b=16,w=16,p=0}. It should have been 0x{b=16,w=16,p=0}, although we will also accept 0x{b=16,w=16,p=0}",+ x, std.flt64bits(rf), y_perfect, y_acceptable)
+ ;;
+ ;;
+}
--- /dev/null
+++ b/lib/math/test/fma-impl.myr
@@ -1,0 +1,104 @@
+use std
+use math
+use testr
+
+const main = {+ testr.run([
+ [.name="fma-01", .fn = fma01],
+ [.name="fma-02", .fn = fma02],
+ ][:])
+}
+
+const fma01 = {c+ var inputs : (uint32, uint32, uint32, uint32)[:] = [
+ /*
+ These are mostly obtained by running fpmath-consensus
+ with seed 1234. Each (mostly) covers a different
+ corner case.
+ */
+ (0x000009a4, 0x00000000, 0x00000002, 0x00000002),
+ (0x69000000, 0x90008002, 0x68348026, 0x68348026),
+ (0x334802ab, 0x49113e8d, 0x90aea62e, 0x3ce2f4c3),
+ (0x5c35d8c1, 0x12dcb6e2, 0x6c1a8cc2, 0x6c1a8cc2),
+ (0xf6266d83, 0x2b3e04e8, 0x62f99bda, 0x62bbd79e),
+ (0x7278e907, 0x75f6c0f1, 0xf6f9b8e0, 0x7f800000),
+ (0xd7748eeb, 0x6737b23e, 0x68e3bbc7, 0xff2f7c71),
+ (0x7f373de4, 0x3dcf90f0, 0xd22ac17c, 0x7d9492ca),
+ (0xb50fce04, 0x00cd486d, 0x03800000, 0x03800000),
+ (0xbb600000, 0x43b7161a, 0x8684d442, 0xbfa03357),
+ (0xf26f8a00, 0x4bfac000, 0xc74ba9fc, 0xfeeaa06c),
+ (0x55d1fa60, 0x32f20000, 0x1b1fea3d, 0x49467eaf),
+ (0x29e26c00, 0x62352000, 0xa0e845af, 0x4ca032a9),
+ (0x287650f8, 0x7cd00000, 0x94e85d5e, 0x65c821c9),
+ (0x7689f580, 0x91418000, 0xaa2822ae, 0xc8508e21),
+ (0xbd813cc0, 0x421f0000, 0x9f098e17, 0xc0208977),
+ (0x3745461a, 0x4db9b736, 0xb6d7deff, 0x458f1cd8),
+ (0xa3ccfd37, 0x7f800000, 0xed328e70, 0xff800000),
+ (0xa3790205, 0x5033a3e6, 0xa001fd11, 0xb42ebbd5),
+ (0x83dd6ede, 0x31ddf8e6, 0x01fea4c8, 0x01fea4c7),
+ (0xa4988128, 0x099a41ad, 0x00800000, 0x00800000),
+ (0x1e0479cd, 0x91d5fcb4, 0x00800000, 0x00800000),
+ (0x2f413021, 0x0a3f5a4e, 0x80800483, 0x80800000),
+ (0x144dcd10, 0x12f4aba0, 0x80800000, 0x80800000),
+ (0x0d580b86, 0x435768a8, 0x966c8d6f, 0x966c5ffd),
+ (0xa19e9a6f, 0xb49af3e3, 0xa2468b59, 0xa2468b57),
+ (0xd119e996, 0x8e5ad0e3, 0x247e0028, 0x247e83b7),
+ (0x381adbc6, 0x00ee4f61, 0x005f2aeb, 0x005f2d2c),
+ (0x7008233c, 0x2a9613fb, 0x46affd02, 0x5b1f9e8a),
+ (0xe85018a1, 0x2cbd53ed, 0x3fcffab8, 0xd599e668),
+
+ /* These ones are especially tricky */
+ (0x65dbf098, 0xd5beb8b4, 0x7c23db61, 0x73027654),
+ (0xa4932927, 0xc565bc34, 0x316887af, 0x31688bcf),
+ (0xb080a420, 0x09e2e5ca, 0x807ff1bf, 0x80800000),
+ ][:]
+
+ for (x, y, z, r) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var yf : flt32 = std.flt32frombits(y)
+ var zf : flt32 = std.flt32frombits(z)
+ var rf = math.fma(xf, yf, zf)
+ testr.check(c, rf == std.flt32frombits(r),
+ "0x{b=16,w=8,p=0} * 0x{b=16,w=8,p=0} + 0x{b=16,w=8,p=0} should be 0x{b=16,w=8,p=0}, was 0x{b=16,w=8,p=0}",+ x, y, z, r, std.flt32bits(rf))
+ ;;
+}
+
+const fma02 = {c+ var inputs : (uint64, uint64, uint64, uint64)[:] = [
+ /*
+ These are mostly obtained by running fpmath-consensus
+ with seed 1234. Each (mostly) covers a different
+ corner case.
+ */
+ (0x0000000000000000, 0x0000000000000000, 0x0100000000000000, 0x0100000000000000),
+ (0x0000000000000000, 0x0000000000000000, 0x0200000000000000, 0x0200000000000000),
+ (0x00000000000009a4, 0x6900000000000002, 0x6834802690008002, 0x6834802690008002),
+ (0x49113e8d334802ab, 0x5c35d8c190aea62e, 0x6c1a8cc212dcb6e2, 0x6c1a8cc212dcb6e2),
+ (0x2b3e04e8f6266d83, 0xae84e20f62f99bda, 0xc9115a1ccea6ce27, 0xc9115a1ccea6ce27),
+ (0xa03ea9e9b09d932c, 0xded7bc19edcbf0c7, 0xbbc4c1f83b3f8f2e, 0x3f26be5f0c7b48e3),
+ (0xa5ec2141c1e6f339, 0xa2d80fc217f57b61, 0x00b3484b473ef1b8, 0x08d526cb86ee748d),
+ (0xccc6600ee88bb67c, 0xc692eeec9b51cf0f, 0xbf5f1ae3486401b0, 0x536a7a30857129db),
+ (0x5f9b9e449db17602, 0xbef22ae5b6a2b1c5, 0x6133e925e6bf8a12, 0x6133e925e6bf823b),
+ (0x7f851249841b6278, 0x3773388e53a375f4, 0x761c27fc2ffa57be, 0x7709506b0e99dc30),
+ (0x7c7cb20f3ca8af93, 0x800fd7f5cfd5baae, 0x14e4c09c9bb1e17e, 0xbc9c6a3fd0e58682),
+ (0xb5e8db2107f4463f, 0x614af740c0d7eb3b, 0xd7e3d25c4daa81e0, 0xd7e3d798d3ccdffb),
+ (0xae62c8be4cb45168, 0x90cc5236f3516c90, 0x0007f8b14f684558, 0x0007f9364eb1a815),
+ (0x5809f53e32a7e1ba, 0xcc647611ccaa5bf4, 0xdfbdb5c345ce7a56, 0xe480990da5526103),
+ (0xbb889d7f826438e1, 0x03bdaff82129696d, 0x000000dacab276ae, 0x8000009296c962f8),
+ (0x003d95525e2b057a, 0xbef738ea5717d89a, 0x800000089763d88c, 0x800000b456ed1a9c),
+ (0x0be868cb5a7180c8, 0x3357a30707ed947c, 0x80000050d6b86ac6, 0x000000cfa41cb229),
+ (0xbe535f4f8a7498af, 0x00d24adee12217b8, 0x0000005729e93fb0, 0x800000016d975af3),
+ (0x39d1968eb883f088, 0x856f286e3b268f0e, 0x800000d7cdd0ed70, 0x800001e9cf01a0ae),
+ ][:]
+
+ for (x, y, z, r) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var yf : flt64 = std.flt64frombits(y)
+ var zf : flt64 = std.flt64frombits(z)
+ var rf = math.fma(xf, yf, zf)
+ testr.check(c, rf == std.flt64frombits(r),
+ "0x{b=16,w=16,p=0} * 0x{b=16,w=16,p=0} + 0x{b=16,w=16,p=0} should be 0x{b=16,w=16,p=0}, was 0x{b=16,w=16,p=0}",+ x, y, z, r, std.flt64bits(rf))
+ ;;
+}
--- /dev/null
+++ b/lib/math/test/poly-impl.myr
@@ -1,0 +1,66 @@
+use std
+use math
+use testr
+
+const main = {+ testr.run([
+ [.name="horner-01", .fn = horner01],
+ [.name="horner-02", .fn = horner02],
+ ][:])
+}
+
+const horner01 = {c+ var inputs : (uint32, uint32[:], uint32)[:] = [
+ (0x00000000, [][:], 0x00000000),
+ (0xbf800000, [0x3f715545, 0x3f715546, 0x3f715544, 0x3f715545][:], 0xb4000000),
+ (0xbf800001, [0x3f715545, 0x3f715546, 0x3f715544, 0x3f715545][:], 0xb4bc5552),
+ ][:]
+
+ for (f, a, r) : inputs
+ var f2 = std.flt32frombits(f)
+ var a2 = [][:]
+ for aa : a
+ std.slpush(&a2, std.flt32frombits(aa))
+ ;;
+ var sf = math.horner_poly(f2, a2)
+ var s : uint32 = std.flt32bits(sf)
+ testr.eq(c, s, r)
+ ;;
+}
+
+const horner02 = {c+ var inputs : (uint64, uint64[:], uint64)[:] = [
+ (0x0000000000000000, [][:], 0x0000000000000000),
+ (0xc0003d3368ee1111, [
+ 0x3ff0000000000000,
+ 0x3fe0000000000000,
+ 0x3fc5555555555555,
+ 0x3fa5555555555555,
+ 0x3f81a7b9611a7b96,
+ 0x3f59c2d14ee4a102,
+ 0x3f3136f054ff42a4,
+ 0x3f05902ed525b6ee
+ ][:], 0x3fdb64730ab8fa29),
+ (0x40003d3368ee1111, [
+ 0x3ff0000000000000,
+ 0x3fe0000000000000,
+ 0x3fc5555555555555,
+ 0x3fa5555555555555,
+ 0x3f81a7b9611a7b96,
+ 0x3f59c2d14ee4a102,
+ 0x3f3136f054ff42a4,
+ 0x3f05902ed525b6ee
+ ][:], 0x400a331575ee40db),
+ ][:]
+
+ for (f, a, r) : inputs
+ var f2 = std.flt64frombits(f)
+ var a2 = [][:]
+ for aa : a
+ std.slpush(&a2, std.flt64frombits(aa))
+ ;;
+ var sf = math.horner_poly(f2, a2)
+ var s : uint64 = std.flt64bits(sf)
+ testr.eq(c, s, r)
+ ;;
+}
--- /dev/null
+++ b/lib/math/test/round-impl.myr
@@ -1,0 +1,82 @@
+use std
+use math
+use testr
+
+const main = {+ testr.run([
+ [.name = "round-01", .fn = round01],
+ [.name = "round-02", .fn = round02],
+ ][:])
+}
+
+const round01 = {c+ var inputs : (flt32, int32)[:] = [
+ (123.4, 123),
+ (0.0, 0),
+ (-0.0, 0),
+ (1.0, 1),
+ (1.1, 1),
+ (0.9, 1),
+ (15.3, 15),
+ (15.5, 16),
+ (15.7, 16),
+ (16.3, 16),
+ (16.5, 16),
+ (16.7, 17),
+ (-102.1, -102),
+ (-102.5, -102),
+ (-102.7, -103),
+ (-103.1, -103),
+ (-103.5, -104),
+ (-103.7, -104),
+ (2147483641.5, -2147483648),
+ (2147483646.5, -2147483648),
+ (2147483647.5, -2147483648),
+ (2147483649.0, -2147483648),
+ (-2147483641.5, -2147483648),
+ (-2147483646.5, -2147483648),
+ (-2147483647.5, -2147483648),
+ (-2147483649.0, -2147483648),
+ ][:]
+
+ for (f, g) : inputs
+ testr.eq(c, math.rn(f), g)
+ ;;
+}
+
+const round02 = {c+ var inputs : (flt64, int64)[:] = [
+ (123.4, 123),
+ (0.0, 0),
+ (-0.0, 0),
+ (1.0, 1),
+ (1.1, 1),
+ (0.9, 1),
+ (15.3, 15),
+ (15.5, 16),
+ (15.7, 16),
+ (16.3, 16),
+ (16.5, 16),
+ (16.7, 17),
+ (-102.1, -102),
+ (-102.5, -102),
+ (-102.7, -103),
+ (-103.1, -103),
+ (-103.5, -104),
+ (-103.7, -104),
+ (2147483641.5, 2147483642),
+ (2147483646.5, 2147483646),
+ (2147483647.5, 2147483648),
+ (2147483649.0, 2147483649),
+ (-2147483641.5, -2147483642),
+ (-2147483646.5, -2147483646),
+ (-2147483647.5, -2147483648),
+ (-2147483649.0, -2147483649),
+ (9223372036854775806.1, -9223372036854775808),
+ (-9223372036854775806.1, -9223372036854775808),
+ ][:]
+
+ for (f, g) : inputs
+ testr.eq(c, math.rn(f), g)
+ ;;
+}
--- /dev/null
+++ b/lib/math/test/scale2-impl.myr
@@ -1,0 +1,96 @@
+use std
+use math
+use testr
+
+const main = {+ testr.run([
+ [.name = "scale2-01", .fn = scale201],
+ [.name = "scale2-02", .fn = scale202],
+ [.name = "scale2-03", .fn = scale203],
+ [.name = "scale2-04", .fn = scale204],
+ ][:])
+}
+
+const scale201 = {c+ var inputsf : (flt32, int32, flt32)[:] = [
+ (0.0, 1, 0.0),
+ (-0.0, 2, -0.0),
+ (1.0, 3, 8.0),
+ (1.0, -3, 0.125),
+ (23.0, 2, 92.0),
+ (184.2, 10, 188620.8),
+ (0.00000234, 15, 0.07667712),
+ (1834.2, -31, 0.0000008541159331798554),
+ (4321.22341, 0, 4321.22341),
+ ][:]
+
+ for (f, m, g) : inputsf
+ testr.eq(c, math.scale2(f, m), g)
+ ;;
+}
+
+const scale202 = {c+ var inputsf : (flt64, int64, flt64)[:] = [
+ (0.0, 1, 0.0),
+ (-0.0, 2, -0.0),
+ (1.0, 3, 8.0),
+ (1.0, -3, 0.125),
+ (23.0, 2, 92.0),
+ (184.2, 10, 188620.8),
+ (0.00000234, 15, 0.07667712),
+ (1834.2, -31, 0.0000008541159331798554),
+ (4321.22341, 0, 4321.22341),
+ ][:]
+
+ for (f, m, g) : inputsf
+ testr.eq(c, math.scale2(f, m), g)
+ ;;
+}
+
+const scale203 = {c+ var inputsb : (uint32, int32, uint32)[:] = [
+ (0x00000000, 1, 0x00000000),
+ (0x7f38aa32, 0, 0x7f38aa32),
+ (0xaaaaaaaa, 0, 0xaaaaaaaa),
+ (0x43000000, -3, 0x41800000),
+ (0x000030a0, -8, 0x00000031),
+ (0x002f3030, -8, 0x00002f30),
+ (0x032f3030, -20, 0x0000015e),
+ (0x032aafff, -8, 0x00155600),
+ (0x002aafff, 8, 0x03aabffc),
+ (0x0000af31, 2, 0x0002bcc4),
+ (0x0000af31, 260, 0x7eaf3100),
+ (0x0000af31, 266, 0x7f800000),
+ (0x3f7ff404, 128, 0x7f7ff404),
+ ][:]
+
+ for (u, m, v) : inputsb
+ var f = std.flt32frombits(u)
+ var g = math.scale2(f, m)
+ var w = std.flt32bits(g)
+ testr.check(c, v == w, "scale2(0x{w=8,b=16,p=0}, {}) should be 0x{w=8,b=16,p=0}, was 0x{w=8,b=16,p=0}", u, m, v, w)+ ;;
+}
+
+const scale204 = {c+ var inputsb : (uint64, int64, uint64)[:] = [
+ (0x0000000000000000, 1, 0x0000000000000000),
+ (0x7f83785551aa873c, 0, 0x7f83785551aa873c),
+ (0xc2b00000aabbccdd, -1080, 0x800000400002aaef),
+ (0xc644fa802f33cfbd, -1, 0xc634fa802f33cfbd),
+ (0x8004fa802f33cfbd, -1, 0x80027d401799e7de),
+ (0x8004fa8fffffffff, -1, 0x80027d4800000000),
+ (0x0082aaffffffffff, 8, 0x0102aaffffffffff),
+ (0x000000ffffffffff, 1, 0x000001fffffffffe),
+ (0x000000ffffffffff, 1000, 0x3dcfffffffffe000),
+ (0x000000ffffffffff, 2000, 0x7c4fffffffffe000),
+ (0x000000ffffffffff, 2400, 0x7ff0000000000000),
+ ][:]
+
+ for (u, m, v) : inputsb
+ var f = std.flt64frombits(u)
+ var g = math.scale2(f, m)
+ var w = std.flt64bits(g)
+ testr.check(c, v == w, "scale2(0x{w=16,b=16,p=0}, {}) should be 0x{w=16,b=16,p=0}, was 0x{w=16,b=16,p=0}", u, m, v, w)+ ;;
+}
--- /dev/null
+++ b/lib/math/test/sqrt-impl.myr
@@ -1,0 +1,83 @@
+use std
+use math
+use testr
+
+const main = {+ testr.run([
+ [.name="sqrt-01", .fn = sqrt01],
+ [.name="sqrt-02", .fn = sqrt02],
+ ][:])
+}
+
+const sqrt01 = {c+ var inputs : (uint32, uint32)[:] = [
+ (0x00000000, 0x00000000),
+ (0x80000000, 0x80000000),
+ (0x80000001, 0x7ff80000),
+ (0x8aaaaaaa, 0x7ff80000),
+ (0x3f800000, 0x3f800000),
+ (0x40800000, 0x40000000),
+ (0x41100000, 0x40400000),
+ (0x3e800000, 0x3f000000),
+ (0x3a3a0000, 0x3cda35fe),
+ (0x017a1000, 0x207d038b),
+ (0x00fc0500, 0x20339b45),
+ (0x160b0000, 0x2abca321),
+ (0x00000800, 0x1d000000),
+ (0x7f690a00, 0x5f743ff8),
+ (0x7f5c0e00, 0x5f6d590c),
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var yf : flt32 = std.flt32frombits(y)
+ var rf = math.sqrt(xf)
+ testr.check(c, rf == yf,
+ "sqrt(0x{b=16,w=8,p=0}) should be 0x{b=16,w=8,p=0}, was 0x{b=16,w=8,p=0}",+ x, y, std.flt32bits(rf))
+ ;;
+}
+
+const sqrt02 = {c+ var inputs : (uint64, uint64)[:] = [
+ (0x0000000000000000ul, 0x0000000000000000ul),
+ (0x8000000000000000ul, 0x8000000000000000ul),
+ (0x8000000000000001ul, 0x7ff8000000000000ul),
+ (0x8aaaaaaaaaaaaaaaul, 0x7ff8000000000000ul),
+ (0x0606e437817acd16ul, 0x22fb10b36e4ae795ul),
+ (0x3ff0000000000000ul, 0x3ff0000000000000ul),
+ (0x4010000000000000ul, 0x4000000000000000ul),
+ (0x3fd0000000000000ul, 0x3fe0000000000000ul),
+ (0x1bbffa831c8f220eul, 0x2dd69ead9d353d6cul),
+ (0x3f0e0f7339499f0bul, 0x3f7f03d8229c8b81ul),
+ (0x3ca510f548e0f3ecul, 0x3e49f6bcadd1e806ul),
+ (0x044ef24a3cca214bul, 0x221f780430319d58ul),
+ (0x7ab034357a1e0474ul, 0x5d501a0593fd8d49ul),
+ (0x216b2df113b38de7ul, 0x30ad7dcc6f26285aul),
+ (0x2e2de34118496c06ul, 0x370eed0301fdade1ul),
+ (0x155bf26b4fb0b2c8ul, 0x2aa5255cf9bd799cul),
+ (0x4b8004df0ac137aaul, 0x45b6a40fee232f2aul),
+ (0x1acaf23d7b0bf80cul, 0x2d5d5d56beda3392ul),
+ (0x3f97ea4c6399a8e6ul, 0x3fc38fb000d55805ul),
+ (0x78f36ea1656dec48ul, 0x5c71a1fce3f370e4ul),
+ (0x409636d438489edbul, 0x4042da4eeac985aaul),
+ (0x72dfd27869ffd768ul, 0x5966907fc9668f57ul),
+ (0x1f483c585e4f03dcul, 0x2f9bd93c3bd1f884ul),
+ (0x7ade25ea6bb6464eul, 0x5d65f681bdbcdf4eul),
+ (0x24ffe5593b0836dbul, 0x3276973038d3bbddul),
+ (0x03e92ac739ec355eul, 0x21ec60eea1d102e8ul),
+ (0x76b656a961a4f64eul, 0x5b52e7cc1d30f55bul),
+ (0x5bc2fac208381d11ul, 0x4dd8a4f5203ab3d2ul),
+ (0x000578e105ac27aaul, 0x1ff2b6d3204e206eul),
+ (0x00057e1016b7c1edul, 0x1ff2bfae3a8e21bbul),
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var yf : flt64 = std.flt64frombits(y)
+ var rf = math.sqrt(xf)
+ testr.check(c, rf == yf,
+ "sqrt(0x{b=16,w=16,p=0}) should be 0x{b=16,w=16,p=0}, was 0x{b=16,w=16,p=0}",+ x, y, std.flt64bits(rf))
+ ;;
+}
--- /dev/null
+++ b/lib/math/test/sum-impl.myr
@@ -1,0 +1,36 @@
+use std
+use math
+use testr
+
+const main = {+ testr.run([
+ [.name = "sums-kahan-01", .fn = sums_kahan_01],
+ [.name = "sums-priest-01", .fn = sums_priest_01],
+ ][:])
+}
+
+const sums_kahan_01 = {c+ var sums : (flt32[:], flt32)[:] = [
+ ([1.0, 2.0, 3.0][:], 6.0),
+
+ /* Naive summing gives 1.0, actual answer is 2.0 */
+ ([33554432.0, 33554430.0, -16777215.0, -16777215.0, -16777215.0, -16777215.0][:], 3.0)
+ ][:]
+
+ for (a, r) : sums
+ var s = math.kahan_sum(a)
+ testr.eq(c, r, s)
+ ;;
+}
+
+const sums_priest_01 = {c+ var sums : (flt32[:], flt32)[:] = [
+ ([1.0, 2.0, 3.0][:], 6.0),
+ ([33554432.0, 33554430.0, -16777215.0, -16777215.0, -16777215.0, -16777215.0][:], 2.0)
+ ][:]
+
+ for (a, r) : sums
+ var s = math.priest_sum(a)
+ testr.eq(c, r, s)
+ ;;
+}
--- /dev/null
+++ b/lib/math/test/trunc-impl.myr
@@ -1,0 +1,124 @@
+use std
+use math
+use testr
+
+const main = {+ testr.run([
+ [.name = "trunc-01", .fn = trunc01],
+ [.name = "trunc-02", .fn = trunc02],
+ [.name = "floor-01", .fn = floor01],
+ [.name = "floor-02", .fn = floor02],
+ [.name = "ceil-01", .fn = ceil01],
+ [.name = "ceil-02", .fn = ceil02],
+ ][:])
+}
+
+const trunc01 = {c+ var flt32s : (flt32, flt32)[:] = [
+ (123.4, 123.0),
+ (0.0, 0.0),
+ (-0.0, -0.0),
+ (1.0, 1.0),
+ (1.1, 1.0),
+ (0.9, 0.0),
+ (10664524000000000000.0, 10664524000000000000.0),
+ (-3.5, -3.0),
+ (101.999, 101.0),
+ (std.flt32nan(), std.flt32nan()),
+ ][:]
+
+ for (f, g) : flt32s
+ testr.eq(c, math.trunc(f), g)
+ ;;
+}
+
+const trunc02 = {c+ var flt64s : (flt64, flt64)[:] = [
+ (0.0, 0.0),
+ (-0.0, -0.0),
+ (1.0, 1.0),
+ (1.1, 1.0),
+ (0.9, 0.0),
+ (13809453812721350000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.0, 13809453812721350000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.0),
+ (-3.5, -3.0),
+ (101.999, 101.0),
+ (std.flt64nan(), std.flt64nan()),
+ ][:]
+
+ for (f, g) : flt64s
+ testr.eq(c, math.trunc(f), g)
+ ;;
+}
+
+const floor01 = {c+ var flt32s : (flt32, flt32)[:] = [
+ (0.0, 0.0),
+ (-0.0, -0.0),
+ (0.5, 0.0),
+ (1.1, 1.0),
+ (10664524000000000000.0, 10664524000000000000.0),
+ (-3.5, -4.0),
+ (-101.999, -102.0),
+ (-126.999, -127.0),
+ (-127.999, -128.0),
+ (-128.999, -129.0),
+ (std.flt32nan(), std.flt32nan()),
+ ][:]
+
+ for (f, g) : flt32s
+ testr.eq(c, math.floor(f), g)
+ ;;
+}
+
+const floor02 = {c+ var flt64s : (flt64, flt64)[:] = [
+ (0.0, 0.0),
+ (-0.0, -0.0),
+ (0.5, 0.0),
+ (1.1, 1.0),
+ (13809453812721350000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.0, 13809453812721350000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.0),
+ (-3.5, -4.0),
+ (-101.999, -102.0),
+ (std.flt64nan(), std.flt64nan()),
+ ][:]
+
+ for (f, g) : flt64s
+ testr.eq(c, math.floor(f), g)
+ ;;
+}
+
+const ceil01 = {c+ var flt32s : (flt32, flt32)[:] = [
+ (0.0, 0.0),
+ (-0.0, -0.0),
+ (0.5, 1.0),
+ (-0.1, -0.0),
+ (1.1, 2.0),
+ (10664524000000000000.0, 10664524000000000000.0),
+ (-3.5, -3.0),
+ (-101.999, -101.0),
+ (std.flt32nan(), std.flt32nan()),
+ ][:]
+
+ for (f, g) : flt32s
+ testr.eq(c, math.ceil(f), g)
+ ;;
+}
+
+const ceil02 = {c+ var flt64s : (flt64, flt64)[:] = [
+ (0.0, 0.0),
+ (-0.0, -0.0),
+ (0.5, 1.0),
+ (-0.1, -0.0),
+ (1.1, 2.0),
+ (13809453812721350000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.0, 13809453812721350000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.0),
+ (-3.5, -3.0),
+ (-101.999, -101.0),
+ (std.flt64nan(), std.flt64nan()),
+ ][:]
+
+ for (f, g) : flt64s
+ testr.eq(c, math.ceil(f), g)
+ ;;
+}
--- /dev/null
+++ b/lib/math/trunc-impl+posixy-x64-sse4.s
@@ -1,0 +1,41 @@
+.globl math$trunc32
+.globl _math$trunc32
+math$trunc32:
+_math$trunc32:
+ roundss $0x03, %xmm0, %xmm0
+ ret
+
+.globl math$floor32
+.globl _math$floor32
+math$floor32:
+_math$floor32:
+ roundss $0x01, %xmm0, %xmm0
+ ret
+
+.globl math$ceil32
+.globl _math$ceil32
+math$ceil32:
+_math$ceil32:
+ roundss $0x02, %xmm0, %xmm0
+ ret
+
+.globl math$trunc64
+.globl _math$trunc64
+math$trunc64:
+_math$trunc64:
+ roundsd $0x03, %xmm0, %xmm0
+ ret
+
+.globl math$floor64
+.globl _math$floor64
+math$floor64:
+_math$floor64:
+ roundsd $0x01, %xmm0, %xmm0
+ ret
+
+.globl math$ceil64
+.globl _math$ceil64
+math$ceil64:
+_math$ceil64:
+ roundsd $0x02, %xmm0, %xmm0
+ ret
--- /dev/null
+++ b/lib/math/trunc-impl.myr
@@ -1,0 +1,103 @@
+use std
+
+pkg math =
+ pkglocal const trunc32 : (f : flt32 -> flt32)
+ pkglocal const floor32 : (f : flt32 -> flt32)
+ pkglocal const ceil32 : (f : flt32 -> flt32)
+ pkglocal const trunc64 : (f : flt64 -> flt64)
+ pkglocal const floor64 : (f : flt64 -> flt64)
+ pkglocal const ceil64 : (f : flt64 -> flt64)
+;;
+
+const Flt32NegMask : uint32 = (1 << 31)
+const Flt32SigMask : uint32 = (1 << 23) - 1
+
+const Flt64NegMask : uint64 = (1 << 63)
+const Flt64SigMask : uint64 = (1 << 52) - 1
+
+pkglocal const floor32 = {f : flt32+ var n, e, s
+ (n, e, s) = std.flt32explode(f)
+
+ /* Many special cases */
+ if e >= 23 || f == -0.0
+ -> f
+ elif e < 0
+ if n
+ -> -1.0
+ else
+ -> 0.0
+ ;;
+ ;;
+
+ if n
+ var fractional_mask = Flt32SigMask >> (e : uint32)
+ if s & fractional_mask == 0
+ -> f
+ else
+ /* Turns out the packing of exp and sig is useful */
+ var u : uint32 = std.flt32bits(f) & ~fractional_mask
+ u += ((1 << 23) >> (e : uint32))
+ -> std.flt32frombits(u)
+ ;;
+ ;;
+
+ var m : uint32 = (Flt32SigMask >> (e : uint32))
+ -> std.flt32assem(n, e, s & ~m)
+}
+
+pkglocal const trunc32 = {f : flt32+ if std.flt32bits(f) & Flt32NegMask != 0
+ -> -floor32(-f)
+ else
+ -> floor32(f)
+ ;;
+}
+
+pkglocal const ceil32 = {f : flt32+ -> -floor32(-f)
+}
+
+pkglocal const floor64 = {f : flt64+ var n, e, s
+ (n, e, s) = std.flt64explode(f)
+
+ /* Many special cases */
+ if e >= 52 || f == -0.0
+ -> f
+ elif e < 0
+ if n
+ -> -1.0
+ else
+ -> 0.0
+ ;;
+ ;;
+
+ if n
+ var fractional_mask = Flt64SigMask >> (e : uint64)
+ if s & fractional_mask == 0
+ -> f
+ else
+ /* Turns out the packing of exp and sig is useful */
+ var u : uint64 = std.flt64bits(f) & ~fractional_mask
+ u += ((1 << 52) >> (e : uint64))
+ -> std.flt64frombits(u)
+ ;;
+ ;;
+
+ var m : uint64 = (Flt64SigMask >> (e : uint64))
+ -> std.flt64assem(n, e, s & ~m)
+}
+
+pkglocal const trunc64 = {f : flt64+ if std.flt64bits(f) & Flt64NegMask != 0
+ -> -floor64(-f)
+ else
+ -> floor64(f)
+ ;;
+}
+
+pkglocal const ceil64 = {f : flt64+ -> -floor64(-f)
+}
+
--- /dev/null
+++ b/lib/math/util.myr
@@ -1,0 +1,174 @@
+use std
+
+pkg math =
+ const flt32fromflt64 : (f : flt64 -> flt32)
+ const flt64fromflt32 : (x : flt32 -> flt64)
+
+ /* For use in various normalizations */
+ const find_first1_64 : (b : uint64, start : int64 -> int64)
+ const find_first1_64_hl : (h : uint64, l : uint64, start : int64 -> int64)
+
+ /* >> and <<, but without wrapping when the shift is >= 64 */
+ const shr : (u : uint64, s : int64 -> uint64)
+ const shl : (u : uint64, s : int64 -> uint64)
+
+ /* Whether RN() requires incrementing after truncating */
+ const need_round_away : (h : uint64, l : uint64, bitpos_last : int64 -> bool)
+;;
+
+const flt64fromflt32 = {f : flt32+ var n, e, s
+ (n, e, s) = std.flt32explode(f)
+ var xs : uint64 = (s : uint64)
+ var xe : int64 = (e : int64)
+
+ if e == 128
+ -> std.flt64assem(n, 1024, xs)
+ elif e == -127
+ /*
+ All subnormals in single precision (except 0.0s)
+ can be upgraded to double precision, since the
+ exponent range is so much wider.
+ */
+ var first1 = find_first1_64(xs, 23)
+ if first1 < 0
+ -> std.flt64assem(n, -1023, 0)
+ ;;
+ xs = xs << (52 - (first1 : uint64))
+ xe = -126 - (23 - first1)
+ -> std.flt64assem(n, xe, xs)
+ ;;
+
+ -> std.flt64assem(n, xe, xs << (52 - 23))
+}
+
+const flt32fromflt64 = {f : flt64+ var n : bool, e : int64, s : uint64
+ (n, e, s) = std.flt64explode(f)
+ var ts : uint32
+ var te : int32 = (e : int32)
+
+ if e >= 128
+ if e == 1023 && s != 0
+ /* NaN */
+ -> std.flt32assem(n, 128, 1)
+ else
+ /* infinity */
+ -> std.flt32assem(n, 128, 0)
+ ;;
+ ;;
+
+ if e >= -127
+ /* normal */
+ ts = ((s >> (52 - 23)) : uint32)
+ if need_round_away(0, s, 52 - 23)
+ ts++
+ if ts & (1 << 24) != 0
+ ts >>= 1
+ te++
+ ;;
+ ;;
+ if te >= -126
+ -> std.flt32assem(n, te, ts)
+ ;;
+ ;;
+
+ /* subnormal already, will have to go to 0 */
+ if e == -1023
+ -> std.flt32assem(n, -127, 0)
+ ;;
+
+ /* subnormal (at least, it will be) */
+ te = -127
+ var shift : int64 = (52 - 23) + (-126 - e)
+ var ts1 = shr(s, shift)
+ ts = (ts1 : uint32)
+ if need_round_away(0, s, shift)
+ ts++
+ if ts & (1 << 23) != 0
+ /* false alarm, it's normal again */
+ te++
+ ;;
+ ;;
+ -> std.flt32assem(n, te, ts)
+}
+
+/* >> and <<, but without wrapping when the shift is >= 64 */
+const shr = {u : uint64, s : int64+ if (s : uint64) >= 64
+ -> 0
+ else
+ -> u >> (s : uint64)
+ ;;
+}
+
+const shl = {u : uint64, s : int64+ if (s : uint64) >= 64
+ -> 0
+ else
+ -> u << (s : uint64)
+ ;;
+}
+
+/* Find the first 1 bit in a bitstring */
+const find_first1_64 = {b : uint64, start : int64+ for var j = start; j >= 0; --j
+ var m = shl(1, j)
+ if b & m != 0
+ -> j
+ ;;
+ ;;
+
+ -> -1
+}
+
+const find_first1_64_hl = {h, l, start+ var first1_h = find_first1_64(h, start - 64)
+ if first1_h >= 0
+ -> first1_h + 64
+ ;;
+
+ -> find_first1_64(l, 63)
+}
+
+/*
+ For [ h ][ l ], where bitpos_last is the position of the last
+ bit that was included in the truncated result (l's last bit has
+ position 0), decide whether rounding up/away is needed. This is
+ true if
+
+ - following bitpos_last is a 1, then a non-zero sequence, or
+
+ - following bitpos_last is a 1, then a zero sequence, and the
+ round would be to even
+ */
+const need_round_away = {h : uint64, l : uint64, bitpos_last : int64+ var first_omitted_is_1 = false
+ var nonzero_beyond = false
+ if bitpos_last > 64
+ first_omitted_is_1 = h & shl(1, bitpos_last - 1 - 64) != 0
+ nonzero_beyond = nonzero_beyond || h & shr((-1 : uint64), 2 + 64 - (bitpos_last - 64)) != 0
+ nonzero_beyond = nonzero_beyond || (l != 0)
+ else
+ first_omitted_is_1 = l & shl(1, bitpos_last - 1) != 0
+ nonzero_beyond = nonzero_beyond || l & shr((-1 : uint64), 1 + 64 - bitpos_last) != 0
+ ;;
+
+ if !first_omitted_is_1
+ -> false
+ ;;
+
+ if nonzero_beyond
+ -> true
+ ;;
+
+ var hl_is_odd = false
+
+ if bitpos_last >= 64
+ hl_is_odd = h & shl(1, bitpos_last - 64) != 0
+ else
+ hl_is_odd = l & shl(1, bitpos_last) != 0
+ ;;
+
+ -> hl_is_odd
+}
--- a/lib/std/fltbits.myr
+++ b/lib/std/fltbits.myr
@@ -9,10 +9,10 @@
generic isnan : (f : @a -> bool) ::floating @a
const flt64frombits : (bits : uint64 -> flt64)
const flt32frombits : (bits : uint32 -> flt32)
- const flt64explode : (flt : flt64 -> (bool, uint64, int64))
- const flt32explode : (flt : flt32 -> (bool, uint32, int32))
- const flt64assem : (sign : bool, mant : uint64, exp : int64 -> flt64)
- const flt32assem : (sign : bool, mant : uint32, exp : int32 -> flt32)
+ const flt64explode : (flt : flt64 -> (bool, int64, uint64))
+ const flt32explode : (flt : flt32 -> (bool, int32, uint32))
+ const flt64assem : (sign : bool, exp : int64, mant : uint64 -> flt64)
+ const flt32assem : (sign : bool, exp : int32, mant : uint32 -> flt32)
;;
const flt64bits = {flt; -> (&flt : uint64#)#}@@ -20,6 +20,9 @@
const flt64frombits = {bits; -> (&bits : flt64#)#} const flt32frombits = {bits; -> (&bits : flt32#)#}+const Dblbias = 1023
+const Fltbias = 127
+
const flt64explode = {fltvar bits, isneg, mant, uexp, exp
@@ -31,17 +34,11 @@
/* add back the implicit bit if this is not a denormal */
if uexp != 0
mant |= 1ul << 52
- exp = (uexp : int64)
- else
- exp = 1
;;
- /*
- adjust for exponent bias. nb: because we are
- treating the mantissa as m.0 instead of 0.m,
- our exponent bias needs to be offset by the
- size of m
- */
- -> (isneg, mant, exp)
+
+ /* adjust for exponent bias */
+ exp = (uexp : int64) - Dblbias
+ -> (isneg, exp, mant)
}
const flt32explode = {flt@@ -55,33 +52,61 @@
/* add back the implicit bit if this is not a denormal */
if uexp != 0
mant |= 1 << 23
- exp = (uexp : int32)
- else
- exp = 1
;;
- /*
- adjust for exponent bias. nb: because we are
- treating the mantissa as m.0 instead of 0.m,
- our exponent bias needs to be offset by the
- size of m
- */
- -> (isneg, mant, exp)
+
+ /* adjust for exponent bias */
+ exp = (uexp : int32) - Fltbias
+ -> (isneg, exp, mant)
}
-const flt64assem = {sign, mant, exp+const flt64assem = {sign, exp, mantvar s, m, e
+ if exp <= -Dblbias && (mant & (1ul << 52) != 0)
+ var roundup = false
+ var shift : uint64 = ((1 - Dblbias - exp) : uint64)
+ var firstcut = mant & (1 << shift)
+ var restcut = mant & ((1 << shift) - 1)
+ var lastkept = mant & (1 << (shift + 1))
+ roundup = firstcut != 0 && (lastkept != 0 || restcut != 0)
+ mant >>= shift
+ exp = -Dblbias
+ if roundup
+ mant++
+ if (mant & (1ul << 52) != 0)
+ exp++
+ ;;
+ ;;
+ ;;
+
s = (sign : uint64)
- e = (exp : uint64) & 0x7ff
+ e = (exp + Dblbias : uint64) & 0x7ff
m = (mant : uint64) & ((1ul<<52) - 1)
-> std.flt64frombits((s << 63) | (e << 52) | m)
}
-const flt32assem = {sign, mant, exp+const flt32assem = {sign, exp, mantvar s, m, e
+ if exp <= -Fltbias && (mant & (1 << 23) != 0)
+ var roundup = false
+ var shift : uint32 = ((1 - Fltbias - exp) : uint32)
+ var firstcut = mant & (1 << shift)
+ var restcut = mant & ((1 << shift) - 1)
+ var lastkept = mant & (1 << (shift + 1))
+ roundup = firstcut != 0 && (lastkept != 0 || restcut != 0)
+ mant >>= shift
+ exp = -Fltbias
+ if roundup
+ mant++
+ if (mant & (1 << 23) != 0)
+ exp++
+ ;;
+ ;;
+ ;;
+
s = (sign : uint32)
- e = (exp : uint32) & 0xff
+ e = (exp + Fltbias : uint32) & 0xff
m = (mant : uint32) & ((1<<23) - 1)
-> std.flt32frombits(s << 31 | e << 23 | m)
--- a/lib/std/fltfmt.myr
+++ b/lib/std/fltfmt.myr
@@ -24,14 +24,16 @@
const flt64bfmt = {sb, val, mode, precisionvar isneg, exp, mant
- (isneg, mant, exp) = flt64explode(val)
- dragon4(sb, isneg, mant, (exp - 52 : int64), Dblbias, mode, precision)
+ (isneg, exp, mant) = flt64explode(val)
+ exp = max(exp, 1 - Dblbias)
+ dragon4(sb, isneg, mant, exp - 52, Dblbias, mode, precision)
}
const flt32bfmt = {sb, val, mode, precisionvar isneg, exp, mant
- (isneg, mant, exp) = flt32explode(val)
+ (isneg, exp, mant) = flt32explode(val)
+ exp = (max((exp : int64), 1 - Fltbias) : int32)
dragon4(sb, isneg, (mant : uint64), (exp - 23 : int64), Fltbias, mode, precision)
}
@@ -64,9 +66,9 @@
/* initialize */
roundup = false
u = mkbigint(0)
- r = bigshli(mkbigint(f), max(e - p, 0))
- s = bigshli(mkbigint(1), max(0, -(e - p)))
- mm = bigshli(mkbigint(1), max((e - p), 0))
+ r = bigshli(mkbigint(f), max(e, 0))
+ s = bigshli(mkbigint(1), max(0, -e))
+ mm = bigshli(mkbigint(1), max(e, 0))
mp = bigdup(mm)
/* fixup: unequal gaps */
--- a/lib/std/fltparse.myr
+++ b/lib/std/fltparse.myr
@@ -264,7 +264,7 @@
var sign, mant, exp
var za
- (sign, mant, exp) = std.flt64explode(z)
+ (sign, exp, mant) = std.flt64explode(z)
if std.abs((mant : int64) - (1l << 52) - 1) < (lim.nextinc : int64)
mant = 0
exp++
@@ -271,7 +271,7 @@
else
mant += lim.nextinc
;;
- za = std.flt64assem(sign, mant, exp)
+ za = std.flt64assem(sign, exp, mant)
-> za
}
--- a/lib/std/test/fltbits.myr
+++ b/lib/std/test/fltbits.myr
@@ -19,6 +19,8 @@
[.name = "bits-roundtrip-64", .fn = bitsround64],
[.name = "flt32bits", .fn = flt32bits],
[.name = "flt64bits", .fn = flt64bits],
+ [.name = "explode-roundtrip-32", .fn = exploderound32],
+ [.name = "explode-roundtrip-64", .fn = exploderound64],
][:])
}
@@ -95,6 +97,7 @@
(1.0, 0x3f800000),
(0.0000123, 0x374e5c19),
(-993.83, 0xc478751f),
+ (0.000000000000000000000000000000000000006054601, 0x0041edc4),
][:]
var uprime = std.flt32bits(f)
testr.check(c, u == uprime, "flt32bits wrong for {}: 0x{x} != 0x{x}", f, u, uprime)@@ -112,3 +115,65 @@
testr.check(c, u == uprime, "flt64bits wrong for {}: 0x{x} != 0x{x}", f, u, uprime);;
}
+
+const exploderound32 = {c+ for f : [1.0, 0.00001, 123.45, 1111111111111111.2, -1.9, -0.0001, 0.000000000000000000000000000000000000006054601, std.flt32nan()][:]
+ var n, e, s
+ (n, e, s) = std.flt32explode(f)
+ var g = std.flt32assem(n, e, s)
+ testr.check(c, f == g, "assem o explode non-identity: {} != {}", f, g)+ ;;
+
+ /*
+ The exponents and significands need to be rather specific
+ in order for flt32assem to work as expected
+ */
+ for (n, e, s) : [
+ (false, -127, 0),
+ (true, -127, 0),
+ (false, -127, 0x399),
+ (true, -127, 0x23),
+ (false, 45, (1 << 23) | 0x23),
+ (true, -12, (1 << 23) | 0x3a2),
+ (true, -126, (1 << 23) | 0x3a1),
+ (false, -127, 4320708),
+ ][:]
+ var m, f, t
+ (m, f, t) = std.flt32explode(std.flt32assem(n, e, s))
+ testr.check(c, n == m, "explode o assem non-identity: {} != {}", (n, e, s), (m, f, t))+ testr.check(c, e == f, "explode o assem non-identity: {} != {}", (n, e, s), (m, f, t))+ testr.check(c, s == t, "explode o assem non-identity: {} != {}", (n, e, s), (m, f, t))+ ;;
+}
+
+const exploderound64 = {c+ for f : [1.0, 0.00001, 123.45, 1111111111111e+309, -1.9, -0.0001, std.flt64nan()][:]
+ var n, e, s
+ (n, e, s) = std.flt64explode(f)
+ var g = std.flt64assem(n, e, s)
+ testr.check(c, f == g, "assem o explode non-identity: {} != {}", f, g)+ ;;
+
+ /*
+ The exponents and significands need to be rather specific
+ in order for flt32assem to work as expected
+ */
+ for (n, e, s) : [
+ (false, -1023, 0),
+ (true, -1023, 0),
+ (false, -1023, 0x399),
+ (true, -1023, 0x23),
+ (false, 45, (1 << 52) | 0xa33bc),
+ (true, -12, (1 << 52) | 0x3),
+ (true, -200, (1 << 52) | 0x11aabbcc),
+ (true, 543, (1 << 52) | 0x3a1),
+ (true, 1001, (1 << 52) | 0x99aa228),
+ ][:]
+ var m, f, t
+ (m, f, t) = std.flt64explode(std.flt64assem(n, e, s))
+ testr.check(c, n == m, "explode o assem non-identity: {} != {}", (n, e, s), (m, f, t))+ testr.check(c, e == f, "explode o assem non-identity: {} != {}", (n, e, s), (m, f, t))+ testr.check(c, s == t, "explode o assem non-identity: {} != {}", (n, e, s), (m, f, t))+ ;;
+}
+
--- a/lib/std/test/fmt.myr
+++ b/lib/std/test/fmt.myr
@@ -66,6 +66,7 @@
check("7b", "{x}", 123) check("0x7b", "0x{x}", 123) check("0.0", "{}", 0.0)+ check("-0.0", "{}", -0.0) check("0.3", "{}", 0.3) check("0.3", "{}", (0.3 : flt32)) check("1.0", "{}", 1.0)--- a/mbld/opts.myr
+++ b/mbld/opts.myr
@@ -34,7 +34,15 @@
const parseversion : (v : byte[:] -> (int, int, int))
/* not exactly portable, but good enough for now */
+ const CpuidSSE2 : uint64= 0x400000000000000
const CpuidSSE4 : uint64= 0x180000
+
+ /*
+ Intel manuals (vol 1, 14.5.3) say AVX, OSXSAVE also
+ needed. For full portability, XGETBV also needs to be
+ checked, though it isn't right now.
+ */
+ const CpuidFMA : uint64= 0x18001000
extern const cpufeatures : (-> uint64)
;;
--- a/mbld/syssel.myr
+++ b/mbld/syssel.myr
@@ -162,8 +162,14 @@
match opt_arch
| "x64":
tag(b, "x64")
+ if opt_cpufeatures & CpuidSSE2 == CpuidSSE2
+ tag(b, "sse2")
+ ;;
if opt_cpufeatures & CpuidSSE4 == CpuidSSE4
tag(b, "sse4")
+ ;;
+ if opt_cpufeatures & CpuidFMA == CpuidFMA
+ tag(b, "fma")
;;
| unknown:
std.fatal("unknown architecture {}\n", unknown)--- a/mi/flatten.c
+++ b/mi/flatten.c
@@ -114,10 +114,10 @@
static Node *
temp(Flattenctx *flatten, Node *e)
{- Node *t, *dcl;
+ Node *t;
assert(e->type == Nexpr);
- t = gentemp(e->loc, e->expr.type, &dcl);
+ t = gentemp(e->loc, e->expr.type, NULL);
return t;
}
@@ -225,23 +225,20 @@
}
static void
-dispose(Flattenctx *s, Stab *st)
+dispose(Flattenctx *s, Stab *st, size_t n)
{- Node *d, *call, *func, *val;
+ Node *e, *call, *func;
Trait *tr;
Type *ty;
size_t i;
tr = traittab[Tcdisp];
- /* dispose in reverse order of declaration */
- for (i = st->nautodcl; i-- > 0;) {- d = st->autodcl[i];
- ty = decltype(d);
- val = mkexpr(Zloc, Ovar, d->decl.name, NULL);
- val->expr.type = ty;
- val->expr.did = d->decl.did;
+ /* dispose in reverse order of appearance */
+ for (i = st->nautotmp; i-- > n;) {+ e = st->autotmp[i];
+ ty = exprtype(e);
func = traitfn(Zloc, tr, "__dispose__", ty);
- call = mkexpr(Zloc, Ocall, func, val, NULL);
+ call = mkexpr(Zloc, Ocall, func, e, NULL);
call->expr.type = mktype(Zloc, Tyvoid);
flatten(s, call);
}
@@ -460,18 +457,23 @@
/* returns 1 when the exit jump needs to be emitted */
static int
-exitscope(Flattenctx *s, Stab *stop, Srcloc loc, int x)
+exitscope(Flattenctx *s, Stab *stop, Srcloc loc, Exit x)
{+ Node *exit;
Stab *st;
for (st = s->curst;; st = st->super) {- if (st->exit[x]) {- jmp(s, st->exit[x]);
+ exit = st->exit[x];
+ if (st->ndisposed[x] < st->nautotmp) {+ st->exit[x] = genlbl(loc);
+ flatten(s, st->exit[x]);
+ dispose(s, st, st->ndisposed[x]);
+ st->ndisposed[x] = st->nautotmp;
+ }
+ if (exit) {+ jmp(s, exit);
return 0;
}
- st->exit[x] = genlbl(loc);
- flatten(s, st->exit[x]);
- dispose(s, st);
if ((!stop && st->isfunc) || st == stop) {return 1;
}
@@ -503,6 +505,12 @@
r = NULL;
args = n->expr.args;
switch (exprop(n)) {+ case Oauto:
+ r = rval(s, n->expr.args[0]);
+ t = temp(s, r);
+ r = asn(t, r);
+ lappend(&s->curst->autotmp, &s->curst->nautotmp, t);
+ break;
case Osize:
r = n; /* don't touch subexprs; they're a pseudo decl */
break;
@@ -696,7 +704,10 @@
flatten(s, n->block.stmts[i]);
}
assert(s->curst == n->block.scope);
- dispose(s, s->curst);
+ if (st->isfunc)
+ exitscope(s, NULL, Zloc, Xret);
+ else
+ dispose(s, s->curst, 0);
s->curst = st;
}
--- a/parse/dump.c
+++ b/parse/dump.c
@@ -186,7 +186,6 @@
findentf(fd, depth + 1, "isimport=%d\n", n->decl.isimport);
findentf(fd, depth + 1, "isnoret=%d\n", n->decl.isnoret);
findentf(fd, depth + 1, "isexportinit=%d\n", n->decl.isexportinit);
- findentf(fd, depth + 1, "isauto=%d\n", n->decl.isauto);
findentf(fd, depth, ")\n");
outsym(n, fd, depth + 1);
outnode(n->decl.init, fd, depth + 1);
--- a/parse/gram.y
+++ b/parse/gram.y
@@ -148,7 +148,7 @@
%type<node> littok literal lorexpr landexpr borexpr strlit bandexpr
%type<node> cmpexpr addexpr mulexpr shiftexpr prefixexpr ternexpr
%type<node> postfixexpr funclit seqlit tuplit name block stmt label
-%type<node> use fnparam declbody declcore typedeclcore autodecl structent
+%type<node> use fnparam declbody declcore typedeclcore structent
%type<node> arrayelt structelt tuphead ifstmt forstmt whilestmt
%type<node> matchstmt elifs optexprln loopcond optexpr match
@@ -420,8 +420,8 @@
}
;
-declbody: autodecl Tasn expr {$$ = $1; $1->decl.init = $3;}- | autodecl
+declbody: declcore Tasn expr {$$ = $1; $1->decl.init = $3;}+ | declcore
;
declcore: name {$$ = mkdecl($1->loc, $1, mktyvar($1->loc));}@@ -432,10 +432,6 @@
: name Tcolon type {$$ = mkdecl($1->loc, $1, $3);};
-autodecl: Tauto declcore {$$ = $2; $$->decl.isauto = 1;}- | declcore
- ;
-
name : Tident {$$ = mkname($1->loc, $1->id);} | Tident Tdot Tident {$$ = mknsname($3->loc, $1->id, $3->id);
@@ -771,7 +767,8 @@
shiftop : Tbsl | Tbsr;
prefixexpr
- : Tinc prefixexpr {$$ = mkexpr($1->loc, Opreinc, $2, NULL);}+ : Tauto prefixexpr {$$ = mkexpr($1->loc, Oauto, $2, NULL);}+ | Tinc prefixexpr {$$ = mkexpr($1->loc, Opreinc, $2, NULL);} | Tdec prefixexpr {$$ = mkexpr($1->loc, Opredec, $2, NULL);} | Tband prefixexpr {$$ = mkexpr($1->loc, Oaddr, $2, NULL);} | Tlnot prefixexpr {$$ = mkexpr($1->loc, Olnot, $2, NULL);}@@ -930,7 +927,7 @@
| /* empty */ {$$.nl = NULL; $$.nn = 0;};
-fnparam : autodecl {$$ = $1;}+fnparam : declcore {$$ = $1;} | Tgap { $$ = mkpseudodecl($1->loc, mktyvar($1->loc)); } | Tgap Tcolon type { $$ = mkpseudodecl($1->loc, $3); };
--- a/parse/infer.c
+++ b/parse/infer.c
@@ -263,7 +263,7 @@
Type *ty;
tr = traittab[Tcdisp];
- ty = decltype(n);
+ ty = exprtype(n);
assert(tr->nproto == 1);
if (hthas(tr->proto[0]->decl.impls, ty))
return;
@@ -883,7 +883,7 @@
if (update)
bsput(ty->trneed, tr->uid);
return 1;
- }
+ }
if (bshas(tm->traits, tr->uid))
return 1;
if (tm->name && ty->type == Tyname) {@@ -1641,6 +1641,13 @@
infernode(&n->expr.idx, NULL, NULL);
n = checkns(n, np);
switch (exprop(n)) {+ case Oauto: /* @a -> @a */
+ infersub(n, ret, sawret, &isconst);
+ t = type(args[0]);
+ constrain(n, t, traittab[Tcdisp]);
+ n->expr.isconst = isconst;
+ settype(n, t);
+ break;
/* all operands are same type */
case Oadd: /* @a + @a -> @a */
case Osub: /* @a - @a -> @a */
@@ -2058,8 +2065,6 @@
inferdecl(n);
if (hasparams(type(n)) && !ingeneric)
fatal(n, "generic type in non-generic near %s", ctxstr(n));
- if (n->decl.isauto)
- constrain(n, type(n), traittab[Tcdisp]);
popenv(n->decl.env);
indentdepth--;
if (n->decl.isgeneric)
@@ -2618,8 +2623,6 @@
if (streq(declname(n), "__init__"))
if (!initcompatible(tybase(decltype(n))))
fatal(n, "__init__ must be (->void), got %s", tystr(decltype(n)));
- if (n->decl.isauto)
- adddispspecialization(n, curstab());
popenv(n->decl.env);
break;
case Nblock:
@@ -2666,6 +2669,8 @@
settype(n->expr.args[0], exprtype(n));
settype(n->expr.args[0]->expr.args[0], exprtype(n));
}
+ if (exprop(n) == Oauto)
+ adddispspecialization(n, curstab());
for (i = 0; i < n->expr.nargs; i++)
typesub(n->expr.args[i], noerr);
if (!noerr)
@@ -2731,7 +2736,15 @@
for (i = 0; i < nspecializations; i++) {pushstab(specializationscope[i]);
n = specializations[i];
- if (n->type == Nexpr) {+ if (n->type == Nexpr && exprop(n) == Oauto) {+ tr = traittab[Tcdisp];
+ assert(tr->nproto == 1);
+ ty = exprtype(n);
+ dt = mktyfunc(n->loc, NULL, 0, mktype(n->loc, Tyvoid));
+ lappend(&dt->sub, &dt->nsub, ty);
+ d = specializedcl(tr->proto[0], ty, dt, &name);
+ htput(tr->proto[0]->decl.impls, ty, d);
+ } else if (n->type == Nexpr && exprop(n) == Ovar) {d = specializedcl(genericdecls[i], n->expr.param, n->expr.type, &name);
n->expr.args[0] = name;
n->expr.did = d->decl.did;
@@ -2753,14 +2766,6 @@
it = itertype(n->iterstmt.seq, mktype(n->loc, Tyvoid));
d = specializedcl(tr->proto[1], ty, it, &name);
htput(tr->proto[1]->decl.impls, ty, d);
- } else if (n->type == Ndecl && n->decl.isauto) {- tr = traittab[Tcdisp];
- assert(tr->nproto == 1);
- ty = decltype(n);
- dt = mktyfunc(n->loc, NULL, 0, mktype(n->loc, Tyvoid));
- lappend(&dt->sub, &dt->nsub, ty);
- d = specializedcl(tr->proto[0], ty, dt, &name);
- htput(tr->proto[0]->decl.impls, ty, d);
} else { die("unknown node for specialization\n");}
--- a/parse/ops.def
+++ b/parse/ops.def
@@ -1,5 +1,6 @@
/* operator name, is it pure, pretty name */
O(Obad, 1, OTmisc, "BAD")
+O(Oauto, 1, OTpre, "auto")
O(Oadd, 1, OTbin, "+")
O(Osub, 1, OTbin, "-")
O(Omul, 1, OTbin, "*")
@@ -105,4 +106,3 @@
O(Ouge, 1, OTmisc, NULL)
O(Oult, 1, OTmisc, NULL)
O(Oule, 1, OTmisc, NULL)
-
--- a/parse/parse.h
+++ b/parse/parse.h
@@ -108,10 +108,12 @@
Htab *lbl; /* labels */
Htab *impl; /* trait implementations: really a set of implemented traits. */
- Node **autodcl; /* declarations in dcl marked 'auto' */
- size_t nautodcl;
+ /* See mi/flatten.c for the following. */
+ Node **autotmp; /* temporaries for 'auto' expressions */
+ size_t nautotmp;
Node *exit[Nexits];
+ size_t ndisposed[Nexits];
};
struct Tyenv {@@ -331,7 +333,6 @@
char isnoret;
char isexportinit;
char isinit;
- char isauto;
} decl;
struct {--- a/parse/stab.c
+++ b/parse/stab.c
@@ -417,10 +417,6 @@
st = findstab(st, s->decl.name);
old = htget(st->dcl, s->decl.name);
- if (s->decl.isauto) {- assert(!old);
- lappend(&st->autodcl, &st->nautodcl, s);
- }
if (!old)
forcedcl(st, s);
else if (!mergedecl(old, s))
--- a/test/pkgtrait.myr
+++ b/test/pkgtrait.myr
@@ -8,7 +8,6 @@
;;
const main = {- var auto r : regex.regex#
- r = std.try(regex.compile(".*"))+ auto std.try(regex.compile(".*"))std.exit(42)
}