ref: 696f29e1650e45b27a949d513248a52f400531ee
parent: e08060746f8d0f1f42f8eeefb4fa6a699a6d331e
author: S. Gilles <sgilles@math.umd.edu>
date: Wed May 1 14:05:36 EDT 2019
Implement rootn.
--- a/lib/math/bld.sub
+++ b/lib/math/bld.sub
@@ -22,7 +22,7 @@
# polynomial evaluation methods
poly-impl.myr
- # x^n
+ # x^n and x^1/q
log-overkill.myr
pown-impl.myr
--- a/lib/math/fpmath.myr
+++ b/lib/math/fpmath.myr
@@ -25,6 +25,7 @@
/* pown-impl */
pown : (x : @f, n : @i -> @f)
+ rootn : (x : @f, n : @u -> @f)
/* powr-impl */
powr : (x : @f, y : @f -> @f)
@@ -107,6 +108,7 @@
horner_polyu = {x, a; -> horner_polyu32(x, a)} pown = {x, n; -> pown32(x, n)}+ rootn = {x, q; -> rootn32(x, q)} powr = {x, y; -> powr32(x, y)}@@ -146,6 +148,7 @@
horner_polyu = {x, a; -> horner_polyu64(x, a)} pown = {x, n; -> pown64(x, n)}+ rootn = {x, q; -> rootn64(x, q)} powr = {x, y; -> powr64(x, y)}--- a/lib/math/impls.myr
+++ b/lib/math/impls.myr
@@ -35,6 +35,9 @@
pkglocal extern const powr32 : (x : flt32, y : flt32 -> flt32)
pkglocal extern const powr64 : (x : flt64, y : flt64 -> flt64)
+ pkglocal extern const rootn32 : (x : flt32, q : uint32 -> flt32)
+ pkglocal extern const rootn64 : (x : flt64, q : uint64 -> flt64)
+
pkglocal extern const scale232 : (x : flt32, m : int32 -> flt32)
pkglocal extern const scale264 : (x : flt64, m : int64 -> flt64)
--- a/lib/math/pown-impl.myr
+++ b/lib/math/pown-impl.myr
@@ -9,11 +9,14 @@
/*
This is an implementation of pown: computing x^n where n is an
integer. We sort of follow [PEB04], but without their high-radix
- log_2.
+ log_2. Instead, we use log-overkill, which should be good enough.
*/
pkg math =
pkglocal const pown32 : (x : flt32, n : int32 -> flt32)
pkglocal const pown64 : (x : flt64, n : int64 -> flt64)
+
+ pkglocal const rootn32 : (x : flt32, q : uint32 -> flt32)
+ pkglocal const rootn64 : (x : flt64, q : uint64 -> flt64)
;;
type fltdesc(@f, @u, @i) = struct
@@ -105,7 +108,7 @@
elif std.isnan(x)
/* Propagate NaN (why doesn't this come first? Ask IEEE.) */
-> d.frombits(d.nan)
- elif (x == 0.0)
+ elif (x == 0.0 || x == -0.0)
if n < 0 && (n % 2 == 1) && xn
/* (+/- 0)^n = +/- oo */
-> d.frombits(d.neginf)
@@ -196,8 +199,7 @@
var G1, G2
(G1, G2) = double_compensated_sum(ls3[0:6])
- var base1 = exp(G1)
- var base2 = exp(G2)
+ var base = exp(G1) + G2
var pow_xen = xe * n
var pow = pow_xen + I
if pow_xen / n != xe || (I > 0 && d.imax - I < pow_xen) || (I < 0 && d.imin - I > pow_xen)
@@ -215,5 +217,99 @@
;;
;;
- -> ult_sgn * scale2(base1 * base2, pow)
+ -> ult_sgn * scale2(base, pow)
+}
+
+/*
+ Rootn is barely different enough from pown to justify being split out
+ into an entirely separate function.
+ */
+const rootn32 = {x : flt32, q : uint32+ -> rootngen(x, q, desc32)
+}
+
+const rootn64 = {x : flt64, q : uint64+ -> rootngen(x, q, desc64)
+}
+
+generic rootngen = {x : @f, q : @u, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i+ var xb
+ xb = d.tobits(x)
+
+ var xn : bool, xe : @i, xs : @u
+ (xn, xe, xs) = d.explode(x)
+
+ var qf : @f = (q : @f)
+
+ /*
+ Special cases. Note we do not follow IEEE exceptions.
+ */
+ if q == 0
+ /* "for any x (even a zero, quiet NaN, or infinity" */
+ -> 1.0
+ elif std.isnan(x)
+ -> d.frombits(d.nan)
+ elif (x == 0.0 || x == -0.0)
+ if xn && q % 2 == 1
+ /* (+/- 0)^1/q = +/- oo (q odd) */
+ -> d.assem(xn, d.emax, 0)
+ else
+ -> d.frombits(d.inf)
+ ;;
+ elif q == 1
+ /* Anything^1/1 is itself */
+ -> x
+ ;;
+
+ /* As in pown */
+ var ult_sgn = 1.0
+ if xn && (q % 2 == 1)
+ ult_sgn = -1.0
+ ;;
+
+ /* Similar to pown. Let e/q = E + psi, with E an integer.
+
+ x^(1/q) = e^(log(xs)/q) * 2^(e/q)
+ = e^(log(xs)/q) * 2^(psi) * 2^E
+ = e^(log(xs)/q) * e^(log(2) * psi) * 2^E
+ = e^( log(xs)/q + log(2) * psi ) * 2^E
+
+ I've opted to do things just in terms of natural base here because we
+ don't have an integer part, I, that we can slide over in infinite
+ precision.
+ */
+
+ /* Calculate 1/q in very high precision */
+ var r1 = 1.0 / qf
+ var r2 = -math.fma(r1, 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 E : @i
+ if q > std.abs(xe)
+ /* Don't cast q to @i unless we're sure it's in small range */
+ E = 0
+ else
+ E = xe / (q : @i)
+ ;;
+ var qpsi = xe - q * E
+ var psi_hi = (qpsi : @f) / qf
+ 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 G1, G2
+ (G1, G2) = double_compensated_sum(ls1[0:12])
+ /* G1 + G2 approximates log(xs)/q + log(2)*psi */
+
+ var base = exp(G1) + G2
+
+ -> ult_sgn * scale2(base, E)
}
--- a/lib/math/test/pown-impl.myr
+++ b/lib/math/test/pown-impl.myr
@@ -8,6 +8,9 @@
[.name="pown-01", .fn = pown01],
[.name="pown-02", .fn = pown02],
[.name="pown-03", .fn = pown03],
+ [.name="rootn-01", .fn = rootn01],
+ [.name="rootn-02", .fn = rootn02],
+ [.name="rootn-03", .fn = rootn03],
][:])
}
@@ -115,5 +118,351 @@
testr.check(c, rf == zf_perfect || rf == zf_accepted,
"pown(0x{b=16,w=8,p=0}, {}) should be 0x{b=16,w=8,p=0}, will also accept 0x{b=16,w=8,p=0}, was 0x{b=16,w=8,p=0}",x, yi, z_perfect, z_accepted, std.flt32bits(rf))
+ ;;
+}
+
+/* Mostly obtained from fpmath-consensus */
+const rootn01 = {c+ var inputs : (uint32, uint32, uint32)[:] = [
+ (0x69000000, 0x00008002, 0x3f803994),
+ (0x334802ab, 0x00003e8d, 0x3f7fbaf1),
+ (0x5c35d8c1, 0x0000b6e2, 0x3f801be9),
+ (0x9bda2b3e, 0x000062f9, 0xbf7f806b),
+ (0xae84e20f, 0x0000ce27, 0xbf7fe2ca),
+ (0x5a1ccea6, 0x0000c911, 0x3f801786),
+ (0xada99760, 0x0000d3a3, 0xbf7fe22a),
+ (0x7278e907, 0x0000c0f1, 0x3f802eeb),
+ (0xb8e075f6, 0x0000f6f9, 0xbf7ff686),
+ (0x175e379b, 0x00000180, 0x3f5d7ee7),
+ (0xceb6dd80, 0x00003675, 0xbf8031c1),
+ (0xd19204f3, 0x0000a605, 0xbf801359),
+ (0x61017aaf, 0x0000db1f, 0x3f801b25),
+ (0x55b9adba, 0x00007b21, 0x3f80201c),
+ (0x5894d094, 0x0000445b, 0x3f80413f),
+ (0x2a9468ea, 0x00008074, 0x3f7fc64d),
+ (0x4b19a1ed, 0x00004792, 0x3f801cda),
+ (0x2cfab2a7, 0x000077e2, 0x3f7fc936),
+ (0xc2eb245e, 0x000064b7, 0xbf80060f),
+ (0x15aee28d, 0x000081fe, 0x3f7f8e0d),
+ (0x39ac732d, 0x0000d980, 0x3f7ff690),
+ (0x3fb113e7, 0x00000f65, 0x3f8002b3),
+ (0x038d3590, 0x0000753a, 0x3f7f4ad2),
+ (0x035de442, 0x0000133d, 0x3f7bb498),
+ (0x40796a30, 0x00000bdc, 0x3f800eaf),
+ (0x609bff02, 0x0000a14e, 0x3f80247b),
+ (0xcba6522d, 0x00007635, 0xbf80124d),
+ (0x47785e73, 0x00005cc5, 0x3f800f44),
+ (0x2f36f143, 0x00004d5a, 0x3f7fb586),
+ (0x629fabc8, 0x0000169e, 0x3f811503),
+ (0x8378d70a, 0x0000c19b, 0xbf7f9212),
+ (0x0c372300, 0x0000b0ed, 0x3f7f994c),
+ (0x033701e7, 0x0000d215, 0x3f7f9a50),
+ (0x2322aae2, 0x0000c8ff, 0x3f7fce01),
+ (0x27863ba4, 0x00007a79, 0x3f7fba97),
+ (0xd46d669e, 0x0000eb4b, 0xbf800fcd),
+ (0x67a19577, 0x00005d29, 0x3f804c99),
+ (0x097e2e4f, 0x0000229d, 0x3f7dd89e),
+ (0x4c8143e3, 0x00006a31, 0x3f8015be),
+ (0x461238df, 0x00009496, 0x3f8007e1),
+ (0x99e77846, 0x0000bf63, 0xbf7fba5e),
+ (0xec95e5c2, 0x00007177, 0xbf8046a1),
+ (0x65142783, 0x0000a7e0, 0x3f8027c6),
+ (0x08967e92, 0x00004803, 0x3f7ef214),
+ (0xb97681b2, 0x0000bce3, 0xbf7ff4ad),
+ (0x71799aa6, 0x00008d7f, 0x3f803ebe),
+ (0x0af6b9f6, 0x0000db62, 0x3f7fab15),
+ (0xc1e6f339, 0x00002141, 0xbf800cf2),
+ (0x7b61a5ec, 0x000017f5, 0x3f81bec1),
+ (0x484b473e, 0x000000b3, 0x3f89103e),
+ (0x47f3e7ff, 0x000028d5, 0x3f8024cf),
+ (0x34fa7ba5, 0x00000ed1, 0x3f7f049b),
+ (0x912a6fe3, 0x00001fe1, 0xbf7dfea9),
+ (0x1c2132c7, 0x0000dd4c, 0x3f7fc75c),
+ (0x0f734df0, 0x00007463, 0x3f7f6db0),
+ (0x8a793e29, 0x0000ae57, 0xbf7f9429),
+ (0xe8019ac4, 0x00007e21, 0xbf80390a),
+ (0x40a490da, 0x00003d5f, 0x3f80036a),
+ (0x2d62ad71, 0x000003b3, 0x3f794f7e),
+ (0x7ebcc761, 0x0000ff26, 0x3f802c0a),
+ (0x47e20599, 0x0000bc08, 0x3f8007f0),
+ (0xf4cbf264, 0x0000ff79, 0xbf802511),
+ (0xe03f717e, 0x00001cd1, 0xbf80ca8a),
+ (0xb09d932c, 0x0000a9e9, 0xbf7fe0fd),
+ (0xf0c7a03e, 0x0000edcb, 0xbf8024d3),
+ (0x1c6b20bf, 0x0000f715, 0x3f7fcda9),
+ (0x3d601f07, 0x0000044a, 0x3f7f52ce),
+ (0x49d80fd8, 0x0000e98f, 0x3f8007e3),
+ (0x8b914821, 0x000032cf, 0xbf7e966d),
+ (0x48cb418c, 0x00001830, 0x3f80448c),
+ (0x62b61073, 0x0000a215, 0x3f80269e),
+ (0x58f31fb9, 0x000012e9, 0x3f80efce),
+ (0x60152416, 0x0000dbfa, 0x3f801a51),
+ (0x359ba65f, 0x0000e1bf, 0x3f7ff081),
+ (0xd99d5091, 0x0000035f, 0xbf857db8),
+ (0x3a9936f6, 0x0000a6bc, 0x3f7ff5a2),
+ (0x4ae113c3, 0x00001c93, 0x3f8046ea),
+ (0xe7536898, 0x00004fc7, 0xbf8058c9),
+ (0x3978f008, 0x0000f5c4, 0x3f7ff74f),
+ (0x4d34407a, 0x00007eff, 0x3f801337),
+ (0x7b0cf7a7, 0x00005243, 0x3f8080c0),
+ (0x495d7a43, 0x000034f1, 0x3f80212f),
+ (0xcf0fccc6, 0x00009b51, 0xbf8011cf),
+ (0x1ae34864, 0x0000bf5f, 0x3f7fbc30),
+ (0x3a6dea49, 0x00006fab, 0x3f7feff2),
+ (0xba17a03e, 0x0000a27d, 0xbf7ff441),
+ (0x0ef7f2a3, 0x00002d4d, 0x3f7e84f7),
+ (0x50fd68c9, 0x0000959c, 0x3f8014c1),
+ (0xca42a083, 0x00003dc1, 0xbf801f0e),
+ (0x543f6cce, 0x00004ffe, 0x3f802e27),
+ (0x5624f440, 0x0000166f, 0x3f80b3e9),
+ (0xa3383196, 0x0000f9a1, 0xbf7fd7de),
+ (0x333bad9e, 0x0000f204, 0x3f7fee14),
+ (0x99e751e6, 0x0000a0d7, 0xbf7fad26),
+ (0xa99d1aeb, 0x00006ee3, 0xbf7fba1a),
+ (0x13931fa4, 0x000062e8, 0x3f7f62ac),
+ (0x3adb5dff, 0x0000238d, 0x3f7fd1fb),
+ (0x37d89731, 0x00002480, 0x3f7fb5f2),
+ (0x79b51b50, 0x00005246, 0x3f807de0),
+ (0x58e3ea42, 0x0000411e, 0x3f804555),
+ (0xf629153b, 0x00003f79, 0xbf809948),
+ (0x05ee0cf9, 0x0000dd5a, 0x3f7fa3cb),
+ (0xdb9c056f, 0x000036f3, 0xbf805b02),
+ (0xa1727619, 0x0000bb6d, 0xbf7fc725),
+ (0x532b7eb2, 0x00009d85, 0x3f801636),
+ (0xdc287d13, 0x000075c5, 0xbf802b45),
+ (0x98d77f06, 0x00008d67, 0xbf7f9f22),
+ (0x091ffed7, 0x0000f434, 0x3f7fb114),
+ (0x2442ebc7, 0x0000161c, 0x3f7e4ce7),
+ (0x7b6b9715, 0x000003b6, 0x3f8bb33c),
+ (0xf9f57824, 0x00001077, 0xbf827c3e),
+ (0x7f44f304, 0x00006b7b, 0x3f806985),
+ (0x5e594f64, 0x00001c0c, 0x3f80c3f7),
+ (0x4555f61f, 0x000005aa, 0x3f80b86f),
+ (0x7ca236a2, 0x00009020, 0x3f804b66),
+ (0xeb146554, 0x0000e35d, 0xbf80220d),
+ (0x9134b38e, 0x0000a9f7, 0xbf7f9f7f),
+ (0x4f08632d, 0x0000001e, 0x400344f0),
+ (0xc47188f1, 0x0000da49, 0xbf800408),
+ (0x3e834b88, 0x0000b59d, 0x3f7ffe15),
+ (0x24c33b41, 0x00009f22, 0x3f7fc47e),
+ (0x606d0172, 0x00004f12, 0x3f804a04),
+ (0xff34f73f, 0x0000be43, 0xbf803b82),
+ (0x092bdf75, 0x0000c4af, 0x3f7f9e1e),
+ (0xf4f77091, 0x000084f5, 0xbf804772),
+ (0x4805d58d, 0x00003b0f, 0x3f8019a5),
+ (0x4575ae9f, 0x00000be8, 0x3f80591a),
+ (0x74017371, 0x0000f486, 0x3f802620),
+ (0x00fe0261, 0x000041ad, 0x3f7eaf1c),
+ (0xd3b5457f, 0x00006007, 0xbf802571),
+ (0xf2beb37f, 0x00005811, 0xbf806781),
+ (0xe79c1ece, 0x0000f353, 0xbf801d4a),
+ (0x491e82ab, 0x00000d65, 0x3f808025),
+ (0x52b24e0d, 0x0000db4e, 0x3f800f92),
+ (0x34aef047, 0x00003549, 0x3f7fb847),
+ (0x3c71cfbc, 0x0000c215, 0x3f7ffa70),
+ (0xd0760aba, 0x0000a127, 0xbf8012b1),
+ (0x51645923, 0x00007734, 0x3f801aaf),
+ (0x17626bc4, 0x00009713, 0x3f7fa1e5),
+ (0x996d0358, 0x0000a495, 0xbf7fadfe),
+ (0xf04e6f93, 0x0000f993, 0xbf8022bf),
+ (0xc22ed0a6, 0x000097b1, 0xbf800330),
+ (0x32e367fa, 0x00006d48, 0x3f7fd724),
+ (0x78d6dde1, 0x000030c3, 0x3f80d174),
+ (0xb1c91a14, 0x0000a597, 0xbf7fe2b3),
+ (0x0c65fbcc, 0x0000a15d, 0x3f7f8fc3),
+ (0x4d8b7f7d, 0x00004500, 0x3f80242f),
+ (0x5d11fcf5, 0x0000b0b6, 0x3f801dbb),
+ (0x9e6f453d, 0x0000abc1, 0xbf7fbbbf),
+ (0xbbc76737, 0x000068e3, 0xbf7ff38d),
+ ][:]
+
+ for (x, y, z) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var zf : flt32 = std.flt32frombits(z)
+ var rf = math.rootn(xf, y)
+ testr.check(c, rf == zf,
+ "rootn(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, std.flt32bits(rf))
+ ;;
+}
+
+/* Mostly obtained from fpmath-consensus */
+const rootn02 = {c+ var inputs : (uint64, uint64, uint64)[:] = [
+ (0x334802ab68348026, 0x0000000000003e8d, 0x3fefb88966d0fc73),
+ (0xd8c190aea62e4911, 0x0000000000005c35, 0xbff0300baeb0b7c1),
+ (0x6c1a8cc212dcb6e2, 0x0000000000006d83, 0x3ff04833427e3888),
+ (0x9bda2b3e04e8f626, 0x00000000000062f9, 0xbfef7fa3dab09e46),
+ (0x7278e907f500a3f8, 0x000000000000c0f1, 0x3ff02ebee17c0535),
+ (0x379bf6f9b8e075f6, 0x000000000000175e, 0x3fef828a0794d0f4),
+ (0x3675ceb6dd800180, 0x00000000000004f3, 0x3fed742f5217a3e3),
+ (0xadbadb1f61017aaf, 0x00000000000055b9, 0xbfefb4fbe7b380cb),
+ (0x445b5894d0947b21, 0x00000000000068ea, 0x3ff0077cfae0e59f),
+ (0x4b19a1ed80742a94, 0x0000000000004792, 0x3ff01bc8158a397e),
+ (0x245e77e22cfab2a7, 0x000000000000c2eb, 0x3fefcdf63591ded3),
+ (0x81fe15aee28d64b7, 0x000000000000732d, 0xbfef43574db9d30f),
+ (0x3fb113e7d98039ac, 0x0000000000000f65, 0x3feffa5fc962713e),
+ (0x35905a30d561f2bf, 0x000000000000038d, 0x3fec322a1c71a627),
+ (0x133d035de442753a, 0x0000000000006a30, 0x3fef6bf924df7782),
+ (0x609bff020bdc4079, 0x000000000000a14e, 0x3ff0241a72527be1),
+ (0xda117635cba6522d, 0x000000000000bf6d, 0xbff0184bcf841bdd),
+ (0x5cc547785e736a8a, 0x000000000000d5f6, 0x3ff017fc8e894de4),
+ (0x169e629fabc84d5a, 0x000000000000d70a, 0x3fefbc184b5a8ce0),
+ (0x0c372300c19b8378, 0x000000000000b0ed, 0x3fef98eaa8c53873),
+ (0x7a7927863ba4c8ff, 0x000000000000669e, 0x3ff0667d720dc43a),
+ (0x67a19577eb4bd46d, 0x0000000000005d29, 0x3ff04c500bc75bcd),
+ (0x54ada98488faddea, 0x000000000000e716, 0x3ff00ff58f16d301),
+ (0x6e9cea0e69033396, 0x0000000000002e4f, 0x3ff0b6d3444208b8),
+ (0x4c8143e3229d097e, 0x0000000000006a31, 0x3ff0150ead03b40c),
+ (0x651427837177ec95, 0x000000000000a7e0, 0x3ff02773cacfa8b5),
+ (0xc1e6f339db620af6, 0x0000000000002141, 0xbff00a86966000a8),
+ (0x0fc217f57b61a5ec, 0x000000000000a2d8, 0x3fef97ad65938ca6),
+ (0x00b3484b473ef1b8, 0x000000000000e7ff, 0x3fef9fd6d39a83d6),
+ (0xca7aa4908a53ef22, 0x000000000000fa71, 0xbff0077a9e7c016f),
+ (0x1fe1912a6fe31066, 0x00000000000032c7, 0x3fef23011b63aa48),
+ (0x0f734df0dd4c1c21, 0x0000000000007463, 0x3fef6d7d445473f3),
+ (0x3e297530d3b5006a, 0x0000000000008a79, 0x3feffb769ccd4055),
+ (0x7e21e8019ac4ae57, 0x000000000000d927, 0x3ff033242ee77746),
+ (0x40a490dafe3c8f94, 0x0000000000003d5f, 0x3ff0020dc3f0fddd),
+ (0x380703b32d62ad71, 0x000000000000e19e, 0x3feff393518626b3),
+ (0xff267ebcc7616b60, 0x0000000000000599, 0xbffa1908fc246f7e),
+ (0x717eff79f4cbf264, 0x000000000000e03f, 0x3ff02767a208cf56),
+ (0x20bfbbc4c1f83b3f, 0x0000000000001c6b, 0x3fee83a0699e6c36),
+ (0x7d42efec285cf715, 0x0000000000001f07, 0x3ff16e30212edc86),
+ (0x49d80fd8044a3d60, 0x000000000000e98f, 0x3ff0078992781047),
+ (0x482175028dd2b016, 0x0000000000008b91, 0x3ff00a6ed2799c1c),
+ (0x183048cb418c32cf, 0x0000000000001073, 0x3fecd1ca430e9cfe),
+ (0x58f31fb9a21562b6, 0x00000000000012e9, 0x3ff0f1992458de0a),
+ (0xe1bf359ba65ff3a2, 0x0000000000005091, 0xbff04b25b6f847e2),
+ (0x3a9936f6035fd99d, 0x000000000000a6bc, 0x3feff4a79126f616),
+ (0x68981c934ae113c3, 0x000000000000e753, 0x3ff01f4f70f11e3c),
+ (0x7b0cf7a77eff4d34, 0x0000000000005243, 0x3ff081862de5afb8),
+ (0xb67c34f1495d7a43, 0x000000000000e88b, 0xbfeff197aaa3b984),
+ (0x1ae3486401b0c692, 0x000000000000bf5f, 0x3fefbb9658a9354c),
+ (0x131831ae8a34c781, 0x00000000000085ba, 0x3fef89ddfaa2c88c),
+ (0x6fab3a6dea498fa4, 0x000000000000214f, 0x3ff1065cc2526e23),
+ (0xba17a03edd8492fc, 0x000000000000a27d, 0xbfeff3414c393510),
+ (0x0f162d4d0ef7f2a3, 0x00000000000096f5, 0x3fef8df7cd14a0f6),
+ (0x959c50fd68c97d24, 0x000000000000a083, 0xbfefa2f483a03c4a),
+ (0x543f6cce3dc1ca42, 0x0000000000004ffe, 0x3ff02d4dd65dc569),
+ (0x3196166f5624f440, 0x000000000000a338, 0x3fefe0ddd56997df),
+ (0xa99d1aeba0d799e7, 0x0000000000006ee3, 0xbfefb8df548c5ad3),
+ (0xf54c62e813931fa4, 0x000000000000a8e3, 0xbff0387441fa51b1),
+ (0x3e0692bf720b6590, 0x0000000000005dff, 0x3feff8ce1beeb9d7),
+ (0x37d89731238d3adb, 0x0000000000002480, 0x3fefb1c0ca8765cd),
+ (0xea42524679b51b50, 0x00000000000058e3, 0xbff0555ef0ed7477),
+ (0xa114a7c140f2411e, 0x000000000000153b, 0xbfee0c2c31c06409),
+ (0x05ee0cf93f79f629, 0x000000000000dd5a, 0x3fefa3869d2f33b5),
+ (0x32d536f3db9c056f, 0x000000000000f4c1, 0x3fefed07198c0efa),
+ (0xcf5aec86a4a4431a, 0x0000000000007619, 0xbff0173cd52e0e5b),
+ (0x532b7eb2bb6da172, 0x0000000000009d85, 0x3ff015ba24ef5b78),
+ (0x75c5dc287d139f4c, 0x0000000000008e09, 0x3ff043d1a1f65b17),
+ (0x98d77f06d0e8b970, 0x0000000000008d67, 0xbfef9e79eb59770a),
+ (0xfed77aa8d806eae9, 0x000000000000091f, 0xbff5925ef56fc014),
+ (0x161c2442ebc7f434, 0x0000000000009715, 0x3fef9e57929a6625),
+ (0xf9f5782403b67b6b, 0x0000000000001077, 0xbff2a3a17e416195),
+ (0x4f646b7b7f44f304, 0x0000000000005e59, 0x3ff01d2dd9a06afc),
+ (0x05aa4555f61f1c0c, 0x00000000000036a2, 0x3fee8e1b6e4eaa00),
+ (0xeb14655490207ca2, 0x000000000000e35d, 0xbff021d012fbe596),
+ (0x632da9f79134b38e, 0x0000000000004f08, 0x3ff04fe63c107761),
+ (0xb002ee74613c001e, 0x00000000000088f1, 0xbfefd6d71311d064),
+ (0x3e834b88da49c471, 0x000000000000b59d, 0x3feffd3974e40828),
+ (0xc0f99f2224c33b41, 0x000000000000cc4f, 0xbff000e7ce49b2bc),
+ (0x4f12606d01720bee, 0x000000000000f73f, 0x3ff00ae0e034c121),
+ (0x7091c4af092bdf75, 0x000000000000f4f7, 0x3ff02361a1c39170),
+ (0x9cd8a907541784f5, 0x000000000000d58d, 0xbfefc5e62708a0c8),
+ (0x4575ae9f3b0f4805, 0x0000000000000be8, 0x3ff0533646954400),
+ (0x6007d3b5457f41ad, 0x000000000000b37f, 0x3ff01fdadf837895),
+ (0xdb4e52b24e0d0d65, 0x0000000000006a95, 0xbff02dd27fee44ff),
+ (0x34aef04706ecaabf, 0x0000000000003549, 0x3fefb564daef5008),
+ (0x0abac2153c71cfbc, 0x000000000000d076, 0x3fefa5ec51fa7ab7),
+ (0x773451645923a127, 0x0000000000006bc4, 0x3ff05c07325ec692),
+ (0xbad3a495996d0358, 0x0000000000009bb1, 0xbfeff45e4413ba3c),
+ (0xb0b65d11fcf5d6e4, 0x000000000000453d, 0xbfefb25bdc2c0214),
+ (0x6737b23ed7748eeb, 0x000000000000bbc7, 0x3ff0254aea7c80df),
+ (0x5f6916142f797d9e, 0x0000000000009809, 0x3ff024e788316dbd),
+ (0x7a8f5038a7b92cdc, 0x0000000000008bce, 0x3ff04b162df278f5),
+ (0x1f4d0ff0979df02c, 0x0000000000001af7, 0x3fee5d908618bac2),
+ (0x6f03d2de48dee0fe, 0x000000000000e61b, 0x3ff02477f0034973),
+ (0x1e4afd5f6caf1807, 0x000000000000599a, 0x3fef7bd601619e3d),
+ (0xc25a19fe20abf149, 0x000000000000a5ab, 0xbff0029788b92d5a),
+ (0xc062122d42551e02, 0x00000000000028f9, 0xbff001f156db6c8e),
+ (0x2bf4524675375aae, 0x000000000000dfce, 0x3fefe061880887b6),
+ (0x490975d75209d47f, 0x000000000000d119, 0x3ff007bbc73f5f5b),
+ (0x7690977400ec7c06, 0x000000000000449f, 0x3ff08fb9e3610f41),
+ (0x4dce16e24ccb8569, 0x0000000000002163, 0x3ff04a61a70ca92c),
+ (0xbde36bee4aef46b9, 0x00000000000064c3, 0xbfeff8cce7bfe9f2),
+ (0x7b4abc0bc6b1f379, 0x000000000000fd63, 0x3ff029c7c5f425ce),
+ (0xc1be99a58a1bf978, 0x000000000000b6c9, 0xbff001c18a2b9a22),
+ (0x2bcf92298660f438, 0x0000000000005261, 0x3fefa9c0b3d6fb68),
+ (0x692d6449c6d32e79, 0x00000000000016d1, 0x3ff14da388c3aeed),
+ (0x6808089a9e5f3472, 0x000000000000c416, 0x3ff02472d12a9377),
+ (0xe855b67873f90e51, 0x000000000000f169, 0xbff01dce297308a2),
+ (0x7e24f8ae9c6f90d6, 0x000000000000a5fe, 0x3ff0430c1fe9dc3b),
+ (0xa469fabcdcd3543a, 0x000000000000a1bb, 0xbfefc3d622fa82ab),
+ (0x6d77f9a29af7be7c, 0x0000000000004770, 0x3ff072af80d7d6ac),
+ (0x75341f99e487dfbe, 0x000000000000bee6, 0x3ff031d131a50b02),
+ (0xdef4bbea0e26770d, 0x000000000000981f, 0xbff024592e6326e1),
+ (0x0c18acaf50a11a46, 0x000000000000b84b, 0x3fef9cc922213ab4),
+ (0x70980bf8070f97d2, 0x0000000000008817, 0x3ff03ff17f5b76c5),
+ (0xa5963748db6b8d18, 0x000000000000abdb, 0xbfefc9c686a8402e),
+ (0xc625b36e73f2b7ed, 0x000000000000e733, 0xbff004c5d5208b95),
+ (0x774bc3efcd8fcc44, 0x000000000000c45b, 0x3ff0325660227318),
+ (0x04b5e00c65d79618, 0x0000000000009c96, 0x3fef7adf163b5aef),
+ (0x6ae07ba20272ebac, 0x00000000000007fa, 0x3ff433d17cb4b601),
+ (0x50b80190af0a03e2, 0x0000000000009810, 0x3ff013a2b7a6b7b8),
+ (0x27aeff1996633c47, 0x000000000000330d, 0x3fef591fa7d02a66),
+ (0x1c8858751ddd6052, 0x000000000000f0cb, 0x3fefcbfe01ce8a83),
+ (0x5c0c2c03da9fb955, 0x000000000000a53e, 0x3ff01e4d2118b879),
+ (0x2e79bd50977593f9, 0x000000000000eb71, 0x3fefe5ba721d89fb),
+ (0x11b6e7e0eb456ae3, 0x000000000000046b, 0x3fe456c1ddc124da),
+ (0x3c63a15b31f146ff, 0x0000000000007e26, 0x3feff6091df20fa6),
+ (0xc2c11021660c7001, 0x00000000000008f1, 0xbff03850593a5874),
+ (0x0b5642870c5af8c8, 0x000000000000e169, 0x3fefad9c2124ffb0),
+ (0x1f1c6db084f8fdc6, 0x000000000000488e, 0x3fef610352128281),
+ (0x806821a4eebbb7b5, 0x000000000000d13d, 0xbfef94f592506f40),
+ (0x3c4feaf03e236d2f, 0x000000000000dfd5, 0x3feffa410deb8908),
+ (0x51d67587767a81f7, 0x0000000000003e94, 0x3ff03316bb566043),
+ (0x70fffac8f637436f, 0x0000000000006d7a, 0x3ff0504cd3f6b593),
+ (0x73a17306c0360b29, 0x000000000000db6c, 0x3ff02a0517906960),
+ (0x5657a1bdb19ae9f1, 0x000000000000c90d, 0x3ff013d3ac7a2a45),
+ (0xfa7b0c3536a1f012, 0x000000000000f6f5, 0xbff02a48e27b6f86),
+ (0x8bfaa51cf7de3bac, 0x000000000000a2d1, 0xbfef8f88eb007ed7),
+ (0x7e3f5f615830f3b2, 0x000000000000dee6, 0x3ff031e7f58070e6),
+ (0x4cc9285e52ab609c, 0x0000000000009ee2, 0x3ff00e615439065e),
+ (0x003defe634a70af2, 0x000000000000889f, 0x3fef5c3511c56e96),
+ (0xcc639a60b891b351, 0x000000000000daf7, 0xbff00a1b3fdf2a90),
+ (0x5b9933f008846f16, 0x0000000000008463, 0x3ff025402d098c4e),
+ (0xe0bc4cbf6bd74d8f, 0x000000000000bd8b, 0xbff01ed2c4e821fc),
+ (0x31d4f2baa91a9e8e, 0x000000000000244d, 0x3fef774c954b40bf),
+ (0x01283d1c679f5652, 0x0000000000008647, 0x3fef5bc18f5e292f),
+ ][:]
+
+ for (x, y, z) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var zf : flt64 = std.flt64frombits(z)
+ var rf = math.rootn(xf, y)
+ testr.check(c, rf == zf,
+ "rootn(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, std.flt64bits(rf))
+ ;;
+}
+
+/* Quarantine for off-by-one */
+const rootn03 = {c+ var inputs : (uint64, uint64, uint64, uint64)[:] = [
+ (0x0261f48674017371, 0x00000000000000fe, 0x3fb16b9c5a6b95b5, 0x3fb16b9c5a6b95b4),
+ (0x0ed134fa7ba5ed3e, 0x000000000000031e, 0x3fe02b4b4fd4e785, 0x3fe02b4b4fd4e784),
+ (0xfaddaa55205a6a8e, 0x0000000000003edf, 0xbff0a9bf57c870b8, 0xbff0a9bf57c870b7),
+ ][:]
+
+ for (x, q, y_perfect, y_acceptable) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var rf : flt64 = math.rootn(xf, q)
+ var ru : uint64 = std.flt64bits(rf)
+
+ testr.check(c, ru == y_perfect || ru == y_acceptable,
+ "rootn(0x{b=16,w=16,p=0}, {}) should be 0x{b=16,w=16,p=0}, will also accept 0x{b=16,w=16,p=0}, was 0x{b=16,w=16,p=0}",+ x, q, y_perfect, y_acceptable, ru)
;;
}