ref: 601f0b7368ff5970f7e2cdeb88d339586c857cf8
parent: 66a472dc689c3a7b2b8f077fdf01c0dc2fdc1f7d
author: S. Gilles <[email protected]>
date: Fri May 11 10:46:02 EDT 2018
Implement log and log1p.
--- a/lib/math/bld.sub
+++ b/lib/math/bld.sub
@@ -1,7 +1,7 @@
lib math =
fpmath.myr
- # exp
+ # exp and expm1
exp-impl.myr
# rounding (to actual integers)
@@ -11,6 +11,9 @@
# fused-multiply-add
fma-impl+posixy-x64-fma.s
fma-impl.myr
+
+ # log and log1p
+ log-impl.myr
# polynomial evaluation methods
poly-impl.myr
--- a/lib/math/fpmath.myr
+++ b/lib/math/fpmath.myr
@@ -10,6 +10,10 @@
/* fma-impl */
fma : (x : @f, y : @f, z : @f -> @f)
+ /* log-impl */
+ log : (x : @f -> @f)
+ log1p : (x : @f -> @f)
+
/* poly-impl */
horner_poly : (x : @f, a : @f[:] -> @f)
horner_polyu : (x : @f, a : @u[:] -> @f)
@@ -130,6 +134,9 @@
extern const log32 : (x : flt32 -> flt32)
extern const log64 : (x : flt64 -> flt64)
+
+extern const log1p32 : (x : flt32 -> flt32)
+extern const log1p64 : (x : flt64 -> flt64)
extern const horner_poly32 : (x : flt32, a : flt32[:] -> flt32)
extern const horner_poly64 : (x : flt64, a : flt64[:] -> flt64)
--- /dev/null
+++ b/lib/math/log-impl.myr
@@ -1,0 +1,579 @@
+use std
+
+use "fpmath"
+use "util"
+
+/*
+ See [Mul16] (6.2.2) and [Tan90].
+ */
+pkg math =
+ pkglocal const log32 : (x : flt32 -> flt32)
+ pkglocal const log64 : (x : flt64 -> flt64)
+
+ pkglocal const log1p32 : (x : flt32 -> flt32)
+ pkglocal const log1p64 : (x : 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)
+ horner : (f : @f, a : @u[:] -> @f)
+ tobits : (f : @f -> @u)
+ frombits : (u : @u -> @f)
+ sgnmask : @u
+ sig8mask : @u
+ sig8last : @u
+ emin : @i
+ emax : @i
+ precision : @u
+ inf : @u
+ ninf : @u
+ nan : @u
+
+ /* For log */
+ logT1 : @u
+ logT2 : @u
+
+ /* For log1p */
+ log1pT1 : @u
+ log1pT2 : @u
+ T3exp : @u
+
+ /* For procedure 1 */
+ C : (@u, @u)[129]
+ Ai : @u[:]
+
+ /* For procedure 2 */
+ Bi : @u[:]
+ Mtruncmask : @u
+;;
+
+const desc32 : fltdesc(flt32, uint32, int32) = [
+ .explode = std.flt32explode,
+ .assem = std.flt32assem,
+ .horner = horner_polyu32,
+ .tobits = std.flt32bits,
+ .frombits = std.flt32frombits,
+ .sgnmask = (1 << 31),
+ .sig8mask = 0xffff0000, /* Mask to get 8 significant bits */
+ .sig8last = 16, /* Last bit kept when masking */
+ .emin = -126,
+ .emax = 127,
+ .precision = 24,
+ .inf = 0x7f800000,
+ .ninf = 0xff800000,
+ .nan = 0x7fc00000,
+ .logT1 = 0x3f707d5f, /* Just smaller than e^(-1/16) ~= 0.939413 */
+ .logT2 = 0x3f88415b, /* Just larger than e^(1/16) ~= 1.06449 */
+ .log1pT1 = 0xbd782a03, /* Just smaller than e^(-1/16) - 1 ~= -0.0605869 */
+ .log1pT2 = 0x3d8415ac, /* Just larger than e^(1/16) - 1 ~= 0.06449445 */
+ .T3exp = 26, /* Beyond 2^T3exp, 1 + x rounds to x */
+ .C = [ /* Accurate representations for log(1 + j/2^7), all j */
+ (0000000000, 0000000000),
+ (0x3bff0000, 0x3429ac42),
+ (0x3c7e0000, 0x35a8b0fc),
+ (0x3cbdc000, 0x368d83eb),
+ (0x3cfc2000, 0xb6b278c4),
+ (0x3d1cf000, 0x3687b9ff),
+ (0x3d3ba000, 0x3631ec66),
+ (0x3d5a1000, 0x36dd7119),
+ (0x3d785000, 0x35c30046),
+ (0x3d8b2800, 0x365bba8e),
+ (0x3d9a1000, 0xb621a791),
+ (0x3da8d800, 0x34e7e0c3),
+ (0x3db78800, 0xb635d46a),
+ (0x3dc61800, 0x368bac63),
+ (0x3dd49000, 0x36da7496),
+ (0x3de2f000, 0x36a91eb8),
+ (0x3df13800, 0x34edc55e),
+ (0x3dff6800, 0xb6dd9c48),
+ (0x3e06bc00, 0xb44197b9),
+ (0x3e0db800, 0x36ab54be),
+ (0x3e14ac00, 0xb6b41f80),
+ (0x3e1b9000, 0xb4f7f85c),
+ (0x3e226800, 0x36adb32e),
+ (0x3e293800, 0xb650e2f2),
+ (0x3e2ff800, 0x36c1c29e),
+ (0x3e36b000, 0x35fe719d),
+ (0x3e3d5c00, 0x3590210e),
+ (0x3e43fc00, 0x36819483),
+ (0x3e4a9400, 0xb6958c2f),
+ (0x3e511c00, 0x36f07f8b),
+ (0x3e57a000, 0xb6dac5fd),
+ (0x3e5e1400, 0x354e85b2),
+ (0x3e648000, 0xb5838656),
+ (0x3e6ae000, 0x3685ad3f),
+ (0x3e713800, 0x356dc55e),
+ (0x3e778400, 0x36b72f71),
+ (0x3e7dc800, 0x36436af2),
+ (0x3e820200, 0xb6d35a59),
+ (0x3e851a00, 0xb6d8ec63),
+ (0x3e882c00, 0x363f9ae5),
+ (0x3e8b3a00, 0x36e55d5d),
+ (0x3e8e4400, 0x36c60b4d),
+ (0x3e914a00, 0x34fde7bd),
+ (0x3e944a00, 0x36d09ef4),
+ (0x3e974800, 0xb6ea28f7),
+ (0x3e9a3e00, 0x36ecd4c4),
+ (0x3e9d3200, 0x36455694),
+ (0x3ea02200, 0xb6779796),
+ (0x3ea30c00, 0x363c21c6),
+ (0x3ea5f200, 0x36fcabbc),
+ (0x3ea8d600, 0xb693c690),
+ (0x3eabb400, 0xb60e8baa),
+ (0x3eae8e00, 0xb51029fe),
+ (0x3eb16400, 0x353cae72),
+ (0x3eb43600, 0x3601e9b1),
+ (0x3eb70400, 0x366aa2ba),
+ (0x3eb9ce00, 0x36bfb5df),
+ (0x3ebc9600, 0xb6d50116),
+ (0x3ebf5800, 0xb5f88faa),
+ (0x3ec21600, 0x368ed0f4),
+ (0x3ec4d200, 0xb64793ec),
+ (0x3ec78800, 0x36f439b3),
+ (0x3eca3c00, 0x36a0e109),
+ (0x3eccec00, 0x36ac08bf),
+ (0x3ecf9a00, 0xb6e09a03),
+ (0x3ed24200, 0x3410e5bb),
+ (0x3ed4e800, 0xb69b2b30),
+ (0x3ed78a00, 0xb6b66dc4),
+ (0x3eda2800, 0xb6084337),
+ (0x3edcc200, 0x36c4b499),
+ (0x3edf5a00, 0x3659da72),
+ (0x3ee1ee00, 0x36bd3e6d),
+ (0x3ee48000, 0xb6038656),
+ (0x3ee70e00, 0xb687a3d0),
+ (0x3ee99800, 0xb4c0ff8a),
+ (0x3eec2000, 0xb6c6d3af),
+ (0x3eeea400, 0xb6afd9f2),
+ (0x3ef12400, 0x3601a7c7),
+ (0x3ef3a200, 0x351875a2),
+ (0x3ef61c00, 0x36ce9234),
+ (0x3ef89400, 0x3675faf0),
+ (0x3efb0a00, 0xb6e02c7f),
+ (0x3efd7a00, 0x36c47bc8),
+ (0x3effea00, 0xb68fbd40),
+ (0x3f012b00, 0xb6d5a5a3),
+ (0x3f025f00, 0xb444adb2),
+ (0x3f039200, 0xb551f190),
+ (0x3f04c300, 0x36f4f573),
+ (0x3f05f400, 0xb6d1bdad),
+ (0x3f072200, 0x36985d1d),
+ (0x3f085000, 0xb6c61d2b),
+ (0x3f097c00, 0xb6e6a6c1),
+ (0x3f0aa600, 0x35f4bd35),
+ (0x3f0bcf00, 0x36abbd8a),
+ (0x3f0cf700, 0x36568cf9),
+ (0x3f0e1e00, 0xb67c11d8),
+ (0x3f0f4300, 0xb4a18fbf),
+ (0x3f106700, 0xb5cb9b55),
+ (0x3f118a00, 0xb6f28414),
+ (0x3f12ab00, 0xb6062ce1),
+ (0x3f13cb00, 0xb576bb27),
+ (0x3f14ea00, 0xb68013d5),
+ (0x3f160700, 0x369ed449),
+ (0x3f172400, 0xb6bc91c0),
+ (0x3f183f00, 0xb68ccb0f),
+ (0x3f195900, 0xb6cc6ede),
+ (0x3f1a7100, 0x3689d9ce),
+ (0x3f1b8900, 0xb684ab8c),
+ (0x3f1c9f00, 0x34d3562a),
+ (0x3f1db400, 0x36094000),
+ (0x3f1ec800, 0x359a9c56),
+ (0x3f1fdb00, 0xb60f65d2),
+ (0x3f20ec00, 0x36fe8467),
+ (0x3f21fd00, 0xb368318d),
+ (0x3f230c00, 0x36bc21c6),
+ (0x3f241b00, 0xb6c2e157),
+ (0x3f252800, 0xb67449f8),
+ (0x3f263400, 0xb64a0662),
+ (0x3f273f00, 0xb67dc915),
+ (0x3f284900, 0xb6c33fe9),
+ (0x3f295100, 0x36d265bc),
+ (0x3f2a5900, 0x360cf333),
+ (0x3f2b6000, 0xb6454982),
+ (0x3f2c6500, 0x36db5cd8),
+ (0x3f2d6a00, 0x34186b3e),
+ (0x3f2e6e00, 0xb6e2393f),
+ (0x3f2f7000, 0x35aa4906),
+ (0x3f307200, 0xb6d0bb87),
+ (0x3f317200, 0x35bfbe8e), /* Note C[128] is log2 */
+ ],
+ .Ai = [ 0x3daaaac2 ][:], /* Coefficients for log(1 + f/F) */
+ .Bi = [ /* Coefficients for log(1 + f) in terms of a = 2f/(2 + f) */
+ 0x3daaaaa9,
+ 0x3c4d0095,
+ ][:],
+ .Mtruncmask = 0xfffff000, /* Mask to get 12 significant bits */
+]
+
+const desc64 : fltdesc(flt64, uint64, int64) = [
+ .explode = std.flt64explode,
+ .assem = std.flt64assem,
+ .horner = horner_polyu64,
+ .tobits = std.flt64bits,
+ .frombits = std.flt64frombits,
+ .sgnmask = (1 << 63),
+ .sig8mask = 0xffffe00000000000, /* Mask to get 8 significant bits */
+ .sig8last = 45, /* Last bit kept when masking */
+ .emin = -1022,
+ .emax = 1023,
+ .precision = 53,
+ .inf = 0x7ff0000000000000,
+ .ninf = 0xfff0000000000000,
+ .nan = 0x7ff8000000000000,
+ .logT1 = 0x3fee0fabfbc702a3, /* Just smaller than e^(-1/16) ~= 0.939413 */
+ .logT2 = 0x3ff1082b577d34ee, /* Just larger than e^(1/16) ~= 1.06449 */
+ .log1pT1 = 0xbfaf0540428fd5c4, /* Just smaller than e^(-1/16) - 1 ~= -0.0605869 */
+ .log1pT2 = 0x3fb082b577d34ed8, /* Just larger than e^(1/16) - 1 ~= 0.06449445 */
+ .T3exp = 55, /* Beyond 2^T3exp, 1 + x rounds to x */
+ .C = [
+ (000000000000000000, 000000000000000000),
+ (0x3f7fe02a6b200000, 0xbd6f30ee07912df9),
+ (0x3f8fc0a8b1000000, 0xbd5fe0e183092c59),
+ (0x3f97b91b07d80000, 0xbd62772ab6c0559c),
+ (0x3f9f829b0e780000, 0x3d2980267c7e09e4),
+ (0x3fa39e87ba000000, 0xbd642a056fea4dfd),
+ (0x3fa77458f6340000, 0xbd62303b9cb0d5e1),
+ (0x3fab42dd71180000, 0x3d671bec28d14c7e),
+ (0x3faf0a30c0100000, 0x3d662a6617cc9717),
+ (0x3fb16536eea40000, 0xbd60a3e2f3b47d18),
+ (0x3fb341d7961c0000, 0xbd4717b6b33e44f8),
+ (0x3fb51b073f060000, 0x3d383f69278e686a),
+ (0x3fb6f0d28ae60000, 0xbd62968c836cc8c2),
+ (0x3fb8c345d6320000, 0xbd5937c294d2f567),
+ (0x3fba926d3a4a0000, 0x3d6aac6ca17a4554),
+ (0x3fbc5e548f5c0000, 0xbd4c5e7514f4083f),
+ (0x3fbe27076e2a0000, 0x3d6e5cbd3d50fffc),
+ (0x3fbfec9131dc0000, 0xbd354555d1ae6607),
+ (0x3fc0d77e7cd10000, 0xbd6c69a65a23a170),
+ (0x3fc1b72ad52f0000, 0x3d69e80a41811a39),
+ (0x3fc29552f8200000, 0xbd35b967f4471dfc),
+ (0x3fc371fc201f0000, 0xbd6c22f10c9a4ea8),
+ (0x3fc44d2b6ccb0000, 0x3d6f4799f4f6543e),
+ (0x3fc526e5e3a20000, 0xbd62f21746ff8a47),
+ (0x3fc5ff3070a80000, 0xbd6b0b0de3077d7e),
+ (0x3fc6d60fe71a0000, 0xbd56f1b955c4d1da),
+ (0x3fc7ab8902110000, 0xbd537b720e4a694b),
+ (0x3fc87fa065210000, 0xbd5b77b7effb7f41),
+ (0x3fc9525a9cf40000, 0x3d65ad1d904c1d4e),
+ (0x3fca23bc1fe30000, 0xbd62a739b23b93e1),
+ (0x3fcaf3c94e810000, 0xbd600349cc67f9b2),
+ (0x3fcbc286742e0000, 0xbd6cca75818c5dbc),
+ (0x3fcc8ff7c79b0000, 0xbd697794f689f843),
+ (0x3fcd5c216b500000, 0xbd611ba91bbca682),
+ (0x3fce27076e2b0000, 0xbd3a342c2af0003c),
+ (0x3fcef0adcbdc0000, 0x3d664d948637950e),
+ (0x3fcfb9186d5e0000, 0x3d5f1546aaa3361c),
+ (0x3fd0402594b50000, 0xbd67df928ec217a5),
+ (0x3fd0a324e2738000, 0x3d50e35f73f7a018),
+ (0x3fd1058bf9ae8000, 0xbd6a9573b02faa5a),
+ (0x3fd1675cabab8000, 0x3d630701ce63eab9),
+ (0x3fd1c898c1698000, 0x3d59fafbc68e7540),
+ (0x3fd22941fbcf8000, 0xbd3a6976f5eb0963),
+ (0x3fd2895a13de8000, 0x3d3a8d7ad24c13f0),
+ (0x3fd2e8e2bae10000, 0x3d5d309c2cc91a85),
+ (0x3fd347dd9a988000, 0xbd25594dd4c58092),
+ (0x3fd3a64c55698000, 0xbd6d0b1c68651946),
+ (0x3fd4043086868000, 0x3d63f1de86093efa),
+ (0x3fd4618bc21c8000, 0xbd609ec17a426426),
+ (0x3fd4be5f95778000, 0xbd3d7c92cd9ad824),
+ (0x3fd51aad872e0000, 0xbd3f4bd8db0a7cc1),
+ (0x3fd5767717458000, 0xbd62c9d5b2a49af9),
+ (0x3fd5d1bdbf580000, 0x3d4394a11b1c1ee4),
+ (0x3fd62c82f2ba0000, 0xbd6c356848506ead),
+ (0x3fd686c81e9b0000, 0x3d54aec442be1015),
+ (0x3fd6e08eaa2b8000, 0x3d60f1c609c98c6c),
+ (0x3fd739d7f6bc0000, 0xbd67fcb18ed9d603),
+ (0x3fd792a55fdd8000, 0xbd6c2ec1f512dc03),
+ (0x3fd7eaf83b828000, 0x3d67e1b259d2f3da),
+ (0x3fd842d1da1e8000, 0x3d462e927628cbc2),
+ (0x3fd89a3386c18000, 0xbd6ed2a52c73bf78),
+ (0x3fd8f11e87368000, 0xbd5d3881e8962a96),
+ (0x3fd947941c210000, 0x3d56faba4cdd147d),
+ (0x3fd99d9581180000, 0xbd5f753456d113b8),
+ (0x3fd9f323ecbf8000, 0x3d584bf2b68d766f),
+ (0x3fda484090e58000, 0x3d6d8515fe535b87),
+ (0x3fda9cec9a9a0000, 0x3d40931a909fea5e),
+ (0x3fdaf12932478000, 0xbd3e53bb31eed7a9),
+ (0x3fdb44f77bcc8000, 0x3d4ec5197ddb55d3),
+ (0x3fdb985896930000, 0x3d50fb598fb14f89),
+ (0x3fdbeb4d9da70000, 0x3d5b7bf7861d37ac),
+ (0x3fdc3dd7a7cd8000, 0x3d66a6b9d9e0a5bd),
+ (0x3fdc8ff7c79a8000, 0x3d5a21ac25d81ef3),
+ (0x3fdce1af0b860000, 0xbd48290905a86aa6),
+ (0x3fdd32fe7e010000, 0xbd542a9e21373414),
+ (0x3fdd83e7258a0000, 0x3d679f2828add176),
+ (0x3fddd46a04c20000, 0xbd6dafa08cecadb1),
+ (0x3fde24881a7c8000, 0xbd53d9e34270ba6b),
+ (0x3fde744261d68000, 0x3d3e1f8df68dbcf3),
+ (0x3fdec399d2468000, 0x3d49802eb9dca7e7),
+ (0x3fdf128f5faf0000, 0x3d3bb2cd720ec44c),
+ (0x3fdf6123fa700000, 0x3d645630a2b61e5b),
+ (0x3fdfaf588f790000, 0xbd49c24ca098362b),
+ (0x3fdffd2e08580000, 0xbd46cf54d05f9367),
+ (0x3fe02552a5a5c000, 0x3d60fec69c695d7f),
+ (0x3fe04bdf9da94000, 0xbd692d9a033eff75),
+ (0x3fe0723e5c1cc000, 0x3d6f404e57963891),
+ (0x3fe0986f4f574000, 0xbd55be8dc04ad601),
+ (0x3fe0be72e4254000, 0xbd657d49676844cc),
+ (0x3fe0e44985d1c000, 0x3d5917edd5cbbd2d),
+ (0x3fe109f39e2d4000, 0x3d592dfbc7d93617),
+ (0x3fe12f7195940000, 0xbd6043acfedce638),
+ (0x3fe154c3d2f4c000, 0x3d65e9a98f33a396),
+ (0x3fe179eabbd88000, 0x3d69a0bfc60e6fa0),
+ (0x3fe19ee6b467c000, 0x3d52dd98b97baef0),
+ (0x3fe1c3b81f714000, 0xbd3eda1b58389902),
+ (0x3fe1e85f5e704000, 0x3d1a07bd8b34be7c),
+ (0x3fe20cdcd192c000, 0xbd64926cafc2f08a),
+ (0x3fe23130d7bec000, 0xbd17afa4392f1ba7),
+ (0x3fe2555bce990000, 0xbd506987f78a4a5e),
+ (0x3fe2795e1289c000, 0xbd5dca290f81848d),
+ (0x3fe29d37fec2c000, 0xbd5eea6f465268b4),
+ (0x3fe2c0e9ed448000, 0x3d5d1772f5386374),
+ (0x3fe2e47436e40000, 0x3d334202a10c3491),
+ (0x3fe307d7334f0000, 0x3d60be1fb590a1f5),
+ (0x3fe32b1339120000, 0x3d6d71320556b67b),
+ (0x3fe34e289d9d0000, 0xbd6e2ce9146d277a),
+ (0x3fe37117b5474000, 0x3d4ed71774092113),
+ (0x3fe393e0d3564000, 0xbd65e6563bbd9fc9),
+ (0x3fe3b6844a000000, 0xbd3eea838909f3d3),
+ (0x3fe3d9026a714000, 0x3d66faa404263d0b),
+ (0x3fe3fb5b84d18000, 0xbd60bda4b162afa3),
+ (0x3fe41d8fe8468000, 0xbd5aa33736867a17),
+ (0x3fe43f9fe2f9c000, 0x3d5ccef4e4f736c2),
+ (0x3fe4618bc21c4000, 0x3d6ec27d0b7b37b3),
+ (0x3fe48353d1ea8000, 0x3d51bee7abd17660),
+ (0x3fe4a4f85db04000, 0xbd244fdd840b8591),
+ (0x3fe4c679afcd0000, 0xbd61c64e971322ce),
+ (0x3fe4e7d811b74000, 0x3d6bb09cb0985646),
+ (0x3fe50913cc018000, 0xbd6794b434c5a4f5),
+ (0x3fe52a2d265bc000, 0x3d46abb9df22bc57),
+ (0x3fe54b2467998000, 0x3d6497a915428b44),
+ (0x3fe56bf9d5b40000, 0xbd58cd7dc73bd194),
+ (0x3fe58cadb5cd8000, 0xbd49db3db43689b4),
+ (0x3fe5ad404c358000, 0x3d6f2cfb29aaa5f0),
+ (0x3fe5cdb1dc6c0000, 0x3d67648cf6e3c5d7),
+ (0x3fe5ee02a9240000, 0x3d667570d6095fd2),
+ (0x3fe60e32f4478000, 0x3d51b194f912b417),
+ (0x3fe62e42fefa4000, 0xbd48432a1b0e2634),
+ ],
+ .Ai = [
+ 0x3fb5555555550286,
+ 0x3f8999a0bc712416,
+ ][:],
+ .Bi = [
+ 0x3fb55555555554e6,
+ 0x3f89999999bac6d4,
+ 0x3f62492307f1519f,
+ 0x3f3c8034c85dfff0,
+ ][:],
+ .Mtruncmask = 0xfffffffff0000000, /* Mask to get 24 significant bits */
+]
+
+const log32 = {x : flt32+ -> loggen(x, desc32)
+}
+
+const log64 = {x : flt64+ -> loggen(x, desc64)
+}
+
+generic loggen = {x : @f, d : fltdesc(@f, @u, @i) :: numeric,floating @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i+ var b = d.tobits(x)
+ var n : bool, e : @i, s : @u
+ (n, e, s) = d.explode(x)
+
+ /*
+ Special cases for NaN, +/- 0, < 0, inf, and 1. There are
+ certain exceptions (inexact, division by 0, &c) that
+ should be flagged in these cases, which we do not honor
+ currently. See [Tan90].
+ */
+ if std.isnan(x)
+ -> d.frombits(d.nan)
+ elif (b & ~d.sgnmask == 0)
+ -> d.frombits(d.ninf)
+ elif n
+ -> d.frombits(d.nan)
+ elif (b == d.inf)
+ -> x
+ elif std.eq(x, (1.0 : @f))
+ -> (0.0 : @f)
+ ;;
+
+ /* If x is close to 1, polynomial log1p(x - 1) will be sufficient */
+ if (d.logT1 < b && b < d.logT2)
+ -> procedure_2(x - (1.0 : @f), d)
+ ;;
+
+ /*
+ Reduce x to 2^m * (F + f), with (F + f) in [1, 2), so
+ procedure_2's tables work. We also require that F have
+ only 8 significant bits.
+ */
+ var m : @i, Y : @f, F : @f, f : @f
+
+ if e < d.emin
+ /* Normalize significand */
+ var first_1 = find_first1_64((s : uint64), (d.precision : int64))
+ var offset = (d.precision : @u) - 1 - (first_1 : @u)
+ s = s << offset
+ e = d.emin - offset
+ ;;
+
+ m = e
+ Y = d.assem(false, 0, s)
+ if need_round_away(0, (s : uint64), (d.sig8last : int64))
+ F = d.frombits((d.tobits(Y) & d.sig8mask) + (1 << d.sig8last))
+ else
+ F = d.frombits(d.tobits(Y) & d.sig8mask)
+ ;;
+
+ f = Y - F
+
+ -> procedure_1(m, F, f, Y, d)
+}
+
+const log1p32 = {x : flt32+ -> log1pgen(x, desc32)
+}
+
+const log1p64 = {x : flt64+ -> log1pgen(x, desc64)
+}
+
+generic log1pgen = {x : @f, d : fltdesc(@f, @u, @i) :: numeric,floating @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i+ var b = d.tobits(x)
+ var n, e, s
+ (n, e, s) = d.explode(x)
+
+ /*
+ Special cases for NaN, +/- 0, < 0, inf, and 1. There are
+ certain exceptions (inexact, division by 0, &c) that
+ should be flagged in these cases, which we do not honor
+ currently. See [Tan90].
+ */
+ if std.isnan(x)
+ -> d.frombits(d.nan)
+ elif (b & ~d.sgnmask == 0)
+ -> x
+ elif std.eq(x, (-1.0 : @f))
+ -> d.frombits(d.nan | d.sgnmask)
+ elif x < (-1.0 : @f)
+ -> d.frombits(d.nan)
+ elif (b == d.inf)
+ -> x
+ ;;
+
+ /* If x is small enough that 1 + x rounds to 1, return x */
+ if e < (-d.precision : @i)
+ -> x
+ ;;
+
+ /* If x is close to 0, use polynomial */
+ if (n && b < d.log1pT1) || (!n && b < d.log1pT2)
+ -> procedure_2(x, d)
+ ;;
+
+ /*
+ Reduce x m, F, f as in log case. However, since we're
+ approximating 1 + x, more care has to be taken (for
+ example: 1 + x might be infinity).
+ */
+ var Y, m, F, f
+ if e > d.T3exp
+ Y = x
+ else
+ Y = (1.0 : @f) + x
+ ;;
+
+ /*
+ y must be normal, otherwise x would have been -1 +
+ (subnormal), but that would round to -1.
+ */
+ var ny, ey, sy
+ (ny, ey, sy) = d.explode(Y)
+ m = ey
+ Y = d.assem(ny, 0, sy)
+ if need_round_away(0, (sy : uint64), (d.sig8last : int64))
+ F = d.frombits((d.tobits(Y) & d.sig8mask) + (1 << d.sig8last))
+ else
+ F = d.frombits(d.tobits(Y) & d.sig8mask)
+ ;;
+
+ /*
+ f is trickier to compute than in the exp case, because
+ the scale of the 1 is unknown near x.
+ */
+ if m <= -2
+ f = Y - F
+ elif m <= d.precision - 1
+ f = (d.assem(false, -m, 0) - F) + scale2(x, -m)
+ else
+ f = (scale2(x, -m) - F) + d.assem(false, -m, 0)
+ ;;
+
+ -> procedure_1(m, F, f, Y, d)
+}
+
+/* Approximate log(2^m * (F + f)) by tables */
+generic procedure_1 = {m : @i, F : @f, f : @f, Y : @f, d : fltdesc(@f, @u, @i) :: numeric,floating @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i+ /*
+ We must compute log(2^m * (F + f)) = m log(2) + log(F)
+ + log(1 + f/F). Only this last term need be approximated,
+ since log(2) and log(F) may be precomputed.
+
+ For computing log(1 + f/F), [Tan90] gives two alternatives.
+ We choose step 3', which requires floating-point division,
+ but allows us to save approximately 2.5 KiB of precomputed
+ values.
+
+ F is some 1 + j2^(-7), so first we compute j. Note that
+ j could actually be 128 (Ex: x = 0x4effac00.)
+ */
+ var j
+ var nF, eF, sF
+ (nF, eF, sF) = d.explode(F)
+ if eF != 0
+ j = 128
+ else
+ j = 0x7f & (((d.sig8mask & sF) >> d.sig8last) - 0x80)
+ ;;
+
+ var Cu_hi, Cu_lo, log2_hi, log2_lo
+ (Cu_hi, Cu_lo) = d.C[j]
+ (log2_hi, log2_lo) = d.C[128]
+
+ var L_hi = (m : @f) * d.frombits(log2_hi) + d.frombits(Cu_hi)
+ var L_lo = (m : @f) * d.frombits(log2_lo) + d.frombits(Cu_lo)
+
+ var u = ((2.0 : @f) * f)/(Y + F)
+ var v = u * u
+ var q = u * v * d.horner(v, d.Ai)
+
+ -> L_hi + (u + (q + L_lo))
+}
+
+/* Approximate log1p by polynomial */
+generic procedure_2 = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i+ var g = (1.0 : @f)/((2.0 : @f) + f)
+ var u = (2.0 : @f) * f * g
+ var v = u * u
+ var q = u * v * d.horner(v, d.Bi)
+
+ /*
+ 1 / (2 + f) in working precision was good enough for the
+ polynomial evaluation, but to complete the approximation
+ we need to add 2f/(2 + f) with higher precision than
+ working. So we go back and compute better, split u.
+ */
+ var u1 = d.frombits(d.Mtruncmask & d.tobits(u))
+ var f1 = d.frombits(d.Mtruncmask & d.tobits(f))
+ var f2 = f - f1
+ var u2 = (((2.0 : @f) * (f - u1) - u1 * f1) - u1 * f2) * g
+ -> u1 + (u2 + q)
+}
--- a/lib/math/references
+++ b/lib/math/references
@@ -20,8 +20,14 @@
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.
+[Tan90]
+Ping-Tak Peter Tang. “Table-driven Implementation of the Logarithm
+Function in IEEE Floating-point Arithmetic”. In: ACM Trans. Math.
+Softw. 16.4 (Dec. 1990), pp. 378–400. issn: 0098-3500. doi:
+10.1145/98267.98294. url: http://doi.acm.org/10.1145/98267.98294.
+
[Tan92]
Ping Tak Peter Tang. “Table-driven Implementation of the Expm1
-Function in IEEE Floating- point Arithmetic”. In: ACM Trans. Math.
+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/test/log-impl.myr
@@ -1,0 +1,742 @@
+use std
+use math
+use testr
+
+/*
+ Note: a major part of the algorithms are the C constants. They
+ are tested extensively in log{,1p}0{1,2}+ */
+const main = {+ testr.run([
+ [.name="log-01", .fn = log01],
+ [.name="log-02", .fn = log02],
+ [.name="log-03", .fn = log03],
+ [.name="log-04", .fn = log04],
+ [.name="log1p-01", .fn = log1p01],
+ [.name="log1p-02", .fn = log1p02],
+ [.name="log1p-03", .fn = log1p03],
+ [.name="log1p-04", .fn = log1p04],
+ ][:])
+}
+
+const log01 = {c+ var inputs : (uint32, uint32)[:] = [
+ (0x00000000, 0xff800000),
+ (0x01000000, 0xc2ad496b),
+ (0x3f060000, 0xbf25b7eb),
+ (0x3f871b00, 0x3d5d49cd),
+ (0x3f710000, 0xbd77518e),
+ (0x4effac00, 0x41abe3e7),
+ (0x477fb7b6, 0x41316d93),
+ (0x41ff8653, 0x405db02b),
+ (0x7f800000, 0x7f800000),
+ (0x3f800000, 0x00000000),
+ (0x000009a4, 0xc2beef7f),
+ (0x00000002, 0xc2cd2bec), /* C[0] */
+ (0x61017aaf, 0x4239cf35), /* C[1] */
+ (0x5702333f, 0x4202613d), /* C[2] */
+ (0x27837177, 0xc204fa63), /* ... */
+ (0x3603beba, 0xc152415d),
+ (0x6c04905b, 0x4276e689),
+ (0x4805d58d, 0x413d3fca),
+ (0x3e0692bf, 0xc001e117),
+ (0x4f08632d, 0x41ac6883),
+ (0x4b88da49, 0x41859e88),
+ (0x0389bcca, 0xc2a6356c),
+ (0x720b6590, 0x428c2fb3),
+ (0x1c0c5e59, 0xc2447c1e),
+ (0x7b0cf7a7, 0x42a5297b),
+ (0x488e1f1c, 0x41494d03),
+ (0x7a8f5038, 0x42a3cf0a),
+ (0x35905a30, 0xc15be22b),
+ (0x49113e8d, 0x4154bd2b),
+ (0x5811f2be, 0x420861b9),
+ (0x73129cb1, 0x428f0f52),
+ (0x6f93b5c4, 0x42855ee6),
+ (0x5894d094, 0x420b3b6c),
+ (0x08967e92, 0xc2982b29),
+ (0x3396e716, 0xc183c476),
+ (0x68981c93, 0x42640ae9),
+ (0x6718afeb, 0x425bbd6e),
+ (0x000009a4, 0xc2beef7f),
+ (0x6c1a8cc2, 0x427783ac),
+ (0x609bff02, 0x4237c835),
+ (0x229d097e, 0xc21ffe0a),
+ (0x031e0ed1, 0xc2a751dc),
+ (0x491e82ab, 0x4156232c),
+ (0x629fabc8, 0x4242f72e),
+ (0x7b2155b9, 0x42a56e93),
+ (0x36a205aa, 0xc143daeb),
+ (0x2322aae2, 0xc21d142f),
+ (0x1fa46ee3, 0xc230719c),
+ (0x3ba4c8ff, 0xc0a95cb2),
+ (0x0e26770d, 0xc288b7b7),
+ (0x1d273361, 0xc23e3d6e),
+ (0x01283d1c, 0xc2acbd76),
+ (0x5d2967a1, 0x4224b42b),
+ (0x2c29e575, 0xc1d5ff25),
+ (0x12ab16d1, 0xc2785f53),
+ (0x39ac732d, 0xc10050a6),
+ (0x41ad00fe, 0x4044ba54),
+ (0x522da14e, 0x41cf9c59),
+ (0x6caf1807, 0x427ac941),
+ (0x6a30133d, 0x426cf210),
+ (0x4b309bc3, 0x41821d44),
+ (0x043197b1, 0xc2a45069),
+ (0x00b3484b, 0xc2adffcd),
+ (0x6d347bde, 0x427dae16),
+ (0x68348026, 0x4261f45a),
+ (0x5c35d8c1, 0x421f712d),
+ (0x0c372300, 0xc28e1269),
+ (0x6737b23e, 0x425c7ac2),
+ (0x5c38b8d6, 0x421f813d),
+ (0x753a038d, 0x429514c2),
+ (0x0abac215, 0xc292310f),
+ (0x333bad9e, 0xc187915f),
+ (0x623cc68a, 0x4240dcdc),
+ (0x2b3e04e8, 0xc1e03107),
+ (0x473ef1b8, 0x412cc12a),
+ (0x3b3f8f2e, 0xc0bab99c),
+ (0x7140973b, 0x428a0f6a),
+ (0x44c1d017, 0x40eb152c),
+ (0x2ec2d4dc, 0xc1b92cda),
+ (0x6bc47734, 0x4275b39e),
+ (0x5cc54778, 0x42228a5e),
+ (0x78463aba, 0x429d86ab),
+ (0x1a46981f, 0xc24e2fee),
+ (0x334802ab, 0xc1870f08),
+ (0x68c97d24, 0x42652ac6),
+ (0x464a3cc7, 0x41177e43),
+ (0x48cb418c, 0x414f067e),
+ (0x71cc037f, 0x428b8fcf),
+ (0x2d4d0ef7, 0xc1c966c5),
+ (0x6cce3dc1, 0x427b70e9),
+ (0x2e4f6e9c, 0xc1be3812),
+ (0x32cf8b91, 0xc18c4edd),
+ (0x0d515df5, 0xc28b0818),
+ (0x1cd1e03f, 0xc2401a70),
+ (0x75535300, 0x42955613),
+ (0x0cd3dd45, 0xc28c64ea),
+ (0x1cd566e5, 0xc2400960),
+ (0x57d628d5, 0x4207249d),
+ (0x75d75209, 0x4296c28e),
+ (0x0fd8044a, 0xc28409a1),
+ (0x0358fcb6, 0xc2a6af9e),
+ (0x4d5a2f36, 0x4199fc7c),
+ (0x3adb5dff, 0xc0cc9173),
+ (0x0bdc4079, 0xc28f16d1),
+ (0x12dcb6e2, 0xc2775a87),
+ (0x175e379b, 0xc25e5f89),
+ (0x38df6a31, 0xc1125a5c),
+ (0x3d601f07, 0xc039f502),
+ (0x4ae113c3, 0x417d04b7),
+ (0x77e22cfa, 0x429c674e),
+ (0x6fe31066, 0x42863b0d),
+ (0x486401b0, 0x4145c607),
+ (0x0f653fb1, 0xc2854e14),
+ (0x39659200, 0xc106d3e7),
+ (0x51e6f204, 0x41cc58fc),
+ (0x13e7d980, 0xc2719c90),
+ (0x5f691614, 0x42311213),
+ (0x68ea445b, 0x4265c51f),
+ (0x1c6b20bf, 0xc2426be1),
+ (0x7b6b9715, 0x42a6306d),
+ (0x606d0172, 0x4236aeb7),
+ (0x3a6dea49, 0xc0e026ca),
+ (0x0beecc4f, 0xc28eed6d),
+ (0x056fdd5a, 0xc2a0f0bb),
+ (0x4df0dd4c, 0x41a05296),
+ (0x3c71cfbc, 0xc086e8ac),
+ (0x5e736a8a, 0x422bb2ea),
+ (0x47f3e7ff, 0x413bc301),
+ (0x17f57b61, 0xc25b33cc),
+ (0x3675ceb6, 0xc14846c5),
+ (0x75f6c0f1, 0x42970853),
+ (0x1077f9f5, 0xc2826017),
+ (0x7278e907, 0x428d588a),
+ (0x62f99bda, 0x4244c0af),
+ (0x407af5c4, 0x3faee68b),
+ (0x067c2ef4, 0xc29e114e),
+ (0x7f7d456e, 0x42b16c9b),
+ (0x4ffe543f, 0x41b6f03f),
+ (0x4efebc08, 0x41abdc61),
+ (0x2f7ffc46, 0xc1b17236), /* C[128] */
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var yf : flt32 = std.flt32frombits(y)
+ var rf = math.log(xf)
+ testr.check(c, rf == yf,
+ "log(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))
+ ;;
+
+ for nan : [ 0x7fc00000, 0x7fd00000, 0x7fffffff, 0xc147a5e3 ][:]
+ testr.check(c, std.isnan(math.log(std.flt32frombits(nan))),
+ "log(0x{b=16,w=8,p=0}) should be NaN", nan)+ ;;
+}
+
+const log02 = {c+ var inputs : (uint64, uint64)[:] = [
+ (0x0000000000000000, 0xfff0000000000000),
+ (0x4000000000000000, 0x3fe62e42fefa39ef),
+ (0x7ff0000000000000, 0x7ff0000000000000),
+ (0x3ff0000000000000, 0x0000000000000000),
+ (0x3fee0fabffc702a3, 0xbfafffffbbdf52b4),
+ (0x6834802690008002, 0x407bea2785dd467d),
+ (0x00000000000009a4, 0xc08705080132de98),
+ (0x49113e8d334802ab, 0x4059518f80c8b520),
+ (0x70fffac8f637436f, 0x408100f58e19ab8f),
+ (0x6900000000000002, 0x407c765cf8301757), /* C[0] */
+ (0x6a30133d035de442, 0x407d4927a622132e), /* C[1] */
+ (0x248037d89731238d, 0xc0730472f99e1f4a), /* C[2] */
+ (0x0ab05ff44dd62adb, 0xc082744e51a6f0b5), /* ... */
+ (0x2b108aea88779bcb, 0xc06cef4a2f1c6596),
+ (0x4b309bc3b92e2dfc, 0x405f3371b7645f6c),
+ (0x1990c9fa6a4c0713, 0xc07a98b52fca339e),
+ (0x0740d2f4341bcac4, 0xc083a512fd9739f7),
+ (0x2c4100ab8cb78c63, 0xc06b48fa89cd7e6c),
+ (0x4ae113c3a6bc3a99, 0x405e576b1bac2f98),
+ (0x49113e8d334802ab, 0x4059518f80c8b520),
+ (0x54815a632f7ffc46, 0x406c840d22b38727),
+ (0x482175028dd2b016, 0x4056b8ec87c62e57),
+ (0x3f7197c085acb2f4, 0xc015cd158b41fb6a),
+ (0x4651b338cfcc7ad5, 0x4051b353db2b8859),
+ (0x1cd1e03f717eff79, 0xc07857016bce40ea),
+ (0x0261f48674017371, 0xc0855513d4bd24c2),
+ (0x7b922a8c7d7b8896, 0x4084ab1d75ac1cd5),
+ (0x3fb23d2c3bf35ec8, 0xc0052208742a9d29),
+ (0x5f92651d0e7b356e, 0x4075edf38e954d24),
+ (0x7f7285f6b3e07dfa, 0x4086031261f39110),
+ (0x2322aae2d2150337, 0xc073f62fbbbec476),
+ (0x7eb2bb6da1727619, 0x4085c09e8f1dadd9),
+ (0x71d2e2eb43d8256d, 0x40814a60e0b279cb),
+ (0x74630f734df0dd4c, 0x40822dcdd6310ba8),
+ (0x6fe31066fa71ca7a, 0x40809e8d868734a7),
+ (0x00000000000009a4, 0xc08705080132de98),
+ (0x5e736a8abf6dda11, 0x40752730815a6166),
+ (0x6903781bd44c88f2, 0x407c7980c8a083ab),
+ (0x28b3963c1648a637, 0xc0701a602dded579),
+ (0x6f93b5c49bb1bad3, 0x40808317f139e8e2),
+ (0x6133e925e6bf8a12, 0x40770f91495b0b5f),
+ (0x3c44077c04158248, 0xc04455e5edbac018),
+ (0x6e042c14d7b5903c, 0x407ff14c8c3b9be2),
+ (0x0dc44db612e4d131, 0xc08162df301e5813),
+ (0x4f646b7b7f44f304, 0x40656e70d305a121),
+ (0x6834802690008002, 0x407bea2785dd467d),
+ (0x43d4a2bb177f4778, 0x40459d6223a7d41b),
+ (0x6634cc7aab4228cc, 0x407a877e7a78169b),
+ (0x5894d0947b2155b9, 0x407115cf1b30a32d),
+ (0x5624f4404ffe543f, 0x406ecac8a9c25223),
+ (0x58751ddd6052330d, 0x4070ffdbd350379f),
+ (0x7f75348c0f5498fb, 0x408604275047d724),
+ (0x7265517fa66418ae, 0x40817d410871254b),
+ (0x08158bfe5ae12b55, 0xc0835b01ec0d209d),
+ (0x39659200fba6c321, 0xc0521ed4910ed7c0),
+ (0x7635cba6522da14e, 0x4082cfafdb79820f),
+ (0x5c35d8c190aea62e, 0x407399d2e3f6c83c),
+ (0x36f6035fd99d5091, 0xc058dfa002d358b7),
+ (0x57d628d547f3e7ff, 0x407091b9f622cd4e),
+ (0x78463abad2eed0c7, 0x408386d5e3383f6b),
+ (0x4b76524c54f71c43, 0x405ff7cf887ecb05),
+ (0x524675375aaefb86, 0x40696dcc32c58d99),
+ (0x1806909d6b4d344f, 0xc07ba93c617cc0e8),
+ (0x7306c0360b296d7a, 0x4081b539dff3fe2b),
+ (0x3396e71654ada984, 0xc0611c4da154c729),
+ (0x51e6f204333bad9e, 0x4068e9668d521be8),
+ (0x1af71f4d0ff0979d, 0xc0799f9949497468),
+ (0x1d273361f32be48e, 0xc0781b61d584de4f),
+ (0x05776dbf84a723dc, 0xc084433c2faca744),
+ (0x541784f5f4f77091, 0x406bf2841f114315),
+ (0x2577a13952a873f2, 0xc07258125f58596d),
+ (0x01e7b0ed0c372300, 0xc0857f389ab75fb2),
+ (0x13e7d98039ac732d, 0xc07e8450364cae3a),
+ (0x1077f9f5782403b6, 0xc08073195883aa08),
+ (0x62e813931fa46ee3, 0x40783e0bf5cda73d),
+ (0x131831ae8a34c781, 0xc07f14422da7ff1f),
+ (0x4f08632da9f79134, 0x4064ef09d28eeda0),
+ (0x1eb87975b31b171e, 0xc0770544a50ec830),
+ (0x11989fbebae2eb0f, 0xc0800f1295e3f21e),
+ (0x5c38b8d6d96e1ffe, 0x40739bcd56c393a8),
+ (0x7278e907f500a3f8, 0x4081840b7f2ad445),
+ (0x73f90e51c4166808, 0x4082092d01e05c94),
+ (0x7a7927863ba4c8ff, 0x408449e7d7ea5790),
+ (0x504943add8ae6c70, 0x4066abc8767dbc43),
+ (0x481953e382eee069, 0x4056a4615edb9861),
+ (0x3e297530d3b5006a, 0xc033a307aa45595f),
+ (0x71799aa6bce3b976, 0x40812b8ab6921ace),
+ (0x28f9c062122d4255, 0xc06fd345b377c762),
+ (0x2c29e575a79ea752, 0xc06b67e06908d1f4),
+ (0x753a038d35905a30, 0x4082786127c0c498),
+ (0x3d9a1e1a1f21c1bc, 0xc039d97d98b10ccb),
+ (0x156a3759cb620ff3, 0xc07d78a18c7ae709),
+ (0x1bea5f7f0998cd3b, 0xc078f72383135970),
+ (0x6c1a8cc212dcb6e2, 0x407e9de4bd07d597),
+ (0x089a9e5f347220f2, 0xc0832cf47cf5e0f7),
+ (0x5c1abf7963c10a0c, 0x407386e1b0c37460),
+ (0x0f1ae87bb321fb89, 0xc080ec2b87d49467),
+ (0x407af5c43978f008, 0x4018448cf475c732),
+ (0x12ab16d1692d6449, 0xc07f601524bb16d3),
+ (0x48cb418c32cf8b91, 0x4058910d56870021),
+ (0x720b6590a8e3f54c, 0x40815dfd6485181b),
+ (0x46eb81a434b9f620, 0x40535ecb7294cc14),
+ (0x7b6b9715161c2442, 0x40849dd29d738baf),
+ (0x4eebb8fcb686a19f, 0x4064c6c75b7c96b0),
+ (0x4cfbe7fe8242657f, 0x4062176353551bf1),
+ (0x609bff020bdc4079, 0x4076a61dec83f7c3),
+ (0x0f5c17c620a239e8, 0xc080d5a506dd57ea),
+ (0x4f1c3cfb0a12764c, 0x406509e91a2092fc),
+ (0x2fac53e44bfddb93, 0xc0668ae29cc4834f),
+ (0x371c7eaeb3110aa5, 0xc05876628c3c6820),
+ (0x56fcaa4913caa832, 0x406ff52903fd288b),
+ (0x623cc68a255eedf9, 0x4077c6e7cd2766f5),
+ (0x4bfcde0e60de47ad, 0x4060b5948be5f1a9),
+ (0x285cf7151c6b20bf, 0xc07056a876ed47dd),
+ (0x44bd10aeecadc14d, 0x404aa358796a4dca),
+ (0x7f7d456ec4df8ff2, 0x408606bb7f61f92f),
+ (0x395d6841a6821a8a, 0xc052375b513417c1),
+ (0x33ed717b73129cb1, 0xc060a55c64cbc4e9),
+ (0x667dad261da3f8f5, 0x407ab98af553a20c),
+ (0x42cdc8a5274640b7, 0x403fd020a35da739),
+ (0x3a6dea498fa485ba, 0xc04e883bb08ae507),
+ (0x2b3e04e8f6266d83, 0xc06cafdc166f9da3),
+ (0x488e1f1c6db084f8, 0x4057e601112df1ff),
+ (0x6cce3dc1ca42a083, 0x407f19f86830564f),
+ (0x697e599a1e4afd5f, 0x407cce3d2d79ec2e),
+ (0x491e82abf353e79c, 0x40597613f67eb37a),
+ (0x169eabc19e6f453d, 0xc07ca3673fa3d24d),
+ (0x4efebc0847e20599, 0x4064e04285d5e27e),
+ (0x544eeb58043197b1, 0x406c3dcfe12bf3bb),
+ (0x27aeff1996633c47, 0xc070cf914713d900),
+ (0x6caf1807e61b6f03, 0x407f043c0805dba1),
+ (0x3aef3641195d57ab, 0xc04bbd04d4b5e1f2),
+ (0x2e4f6e9cea0e6903, 0xc0686f887e310813),
+ (0x4f5f7177b7abc69e, 0x40656612da92283f),
+ (0x629fabc84d5a2f36, 0x40780afb4bbfbe4c),
+ (0x418fce8b9cdac427, 0x40320409991d5e65),
+ (0x056fdd5a05ee0cf9, 0xc0844651eb4324d0),
+ (0x70fffac8f637436f, 0x408100f58e19ab8f), /* C[128] */
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var yf : flt64 = std.flt64frombits(y)
+ var rf = math.log(xf)
+ testr.check(c, rf == yf,
+ "log(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))
+ ;;
+
+ for nan : [
+ 0x7ff8000000000000,
+ 0x7ff9000000000000,
+ 0x7fffffffffffffff,
+ 0xc147a5e354789328,
+ ][:]
+ testr.check(c, std.isnan(math.log(std.flt64frombits(nan))),
+ "log(0x{b=16,w=16,p=0}) should be NaN", nan)+ ;;
+}
+
+const log03 = {c+ /*
+ The [Tan90], steps 1-3' implementation have error bounds
+ of about 0.6 ulps, so we do not obtain last-bit accuracy.
+ Here are some known-bad results.
+ */
+
+ var inputs : (uint32, uint32, uint32)[:] = [
+ (0x3f610400, 0xbe041a91, 0xbe041a92),
+ (0x3fc70700, 0x3ee200bf, 0x3ee200c0),
+ (0x3f610400, 0xbe041a91, 0xbe041a92),
+ (0x3e360700, 0xbfdd18a7, 0xbfdd18a8),
+ ][:]
+
+ 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.log(xf)
+ if rf != ypf && rf != yaf
+ testr.fail(c, "log(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 log04 = {c+ /*
+ The [Tan90], steps 1-3' implementation have error bounds
+ of about 0.6 ulps, so we do not obtain last-bit accuracy.
+ Here are some known-bad results.
+ */
+
+ var inputs : (uint64, uint64, uint64)[:] = [
+ (0x3c71cfbc354934ae, 0xc0435ac0222f1703, 0xc0435ac0222f1704),
+ (0x35f0681e2059a1bb, 0xc05bb8387abe5fcf, 0xc05bb8387abe5fd0),
+ (0x40d268e6c4ad9588, 0x4023b04f15e91586, 0x4023b04f15e91585),
+ ][:]
+
+ 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.log(xf)
+ if rf != ypf && rf != yaf
+ testr.fail(c, "log(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 log1p01 = {c+ var inputs : (uint32, uint32)[:] = [
+ (0x00000000, 0x00000000),
+ (0xbf700700, 0xc0318e1e),
+ (0x69000000, 0x42661ff7), /* C[0] */
+ (0x61017aaf, 0x4239cf35), /* C[1] */
+ (0x5702333f, 0x4202613d), /* C[2] */
+ (0x6e030844, 0x4280f8e2), /* ... */
+ (0xbf5f1ae3, 0xc00351a2),
+ (0x6c04905b, 0x4276e689),
+ (0x4805d58d, 0x413d3fd2),
+ (0x7c06e112, 0x42a7d8a8),
+ (0x4f08632d, 0x41ac6883),
+ (0x4b88da49, 0x41859e88),
+ (0x6b898457, 0x42744651),
+ (0x720b6590, 0x428c2fb3),
+ (0x5c0c2c03, 0x421e66a2),
+ (0x7b0cf7a7, 0x42a5297b),
+ (0x488e1f1c, 0x41494d06),
+ (0x7a8f5038, 0x42a3cf0a),
+ (0x768fe527, 0x4298b9fb),
+ (0x49113e8d, 0x4154bd2d),
+ (0x5811f2be, 0x420861b9),
+ (0x73129cb1, 0x428f0f52),
+ (0x6f93b5c4, 0x42855ee6),
+ (0x5894d094, 0x420b3b6c),
+ (0x3facbad0, 0x3f5aaba7),
+ (0x6f16daf7, 0x428406cc),
+ (0x68981c93, 0x42640ae9),
+ (0x6718afeb, 0x425bbd6e),
+ (0x4b19a1ed, 0x4180ffd5),
+ (0x6c1a8cc2, 0x427783ac),
+ (0x609bff02, 0x4237c835),
+ (0x407af5c4, 0x3fcbf9dc),
+ (0x459de4ba, 0x41087217),
+ (0x491e82ab, 0x4156232d),
+ (0x629fabc8, 0x4242f72e),
+ (0x7b2155b9, 0x42a56e93),
+ (0x7ba20272, 0x42a6d39a),
+ (0x44229f49, 0x40cf561a),
+ (0x782403b6, 0x429d25a9),
+ (0x5624f440, 0x41fb8fe5),
+ (0x7e263c63, 0x42adcf3f),
+ (0x6026873d, 0x42354554),
+ (0x52a873f2, 0x41d4e9ec),
+ (0x5d2967a1, 0x4224b42b),
+ (0x41a1dbe2, 0x40438dc0),
+ (0x52ab609c, 0x41d50d2b),
+ (0x632ba15e, 0x424606ed),
+ (0x692d6449, 0x426756c6),
+ (0x522da14e, 0x41cf9c59),
+ (0x6caf1807, 0x427ac941),
+ (0x6a30133d, 0x426cf210),
+ (0x4b309bc3, 0x41821d44),
+ (0x47b19da7, 0x4136aff5),
+ (0x7eb2bb6d, 0x42af573f),
+ (0x6d347bde, 0x427dae16),
+ (0x68348026, 0x4261f45a),
+ (0x5c35d8c1, 0x421f712d),
+ (0x4fb68eb5, 0x41b44934),
+ (0x6737b23e, 0x425c7ac2),
+ (0x5c38b8d6, 0x421f813d),
+ (0x753a038d, 0x429514c2),
+ (0x733b43e1, 0x428f8ca0),
+ (0x613c001e, 0x423b4d15),
+ (0x623cc68a, 0x4240dcdc),
+ (0x6e3e4d6b, 0x4281b7f2),
+ (0x473ef1b8, 0x412cc13f),
+ (0x50c01f91, 0x41bfc8ef),
+ (0x7140973b, 0x428a0f6a),
+ (0x44c1d017, 0x40eb1a74),
+ (0x7d42efec, 0x42ab5b02),
+ (0x6bc47734, 0x4275b39e),
+ (0x5cc54778, 0x42228a5e),
+ (0x78463aba, 0x429d86ab),
+ (0x50473690, 0x41ba8795),
+ (0x4fc7e753, 0x41b5031a),
+ (0x68c97d24, 0x42652ac6),
+ (0x464a3cc7, 0x41177e94),
+ (0x48cb418c, 0x414f0681),
+ (0x71cc037f, 0x428b8fcf),
+ (0x6b4d344f, 0x42731a66),
+ (0x6cce3dc1, 0x427b70e9),
+ (0xbe43ff34, 0xbe598dc6),
+ (0x42cdc8a5, 0x40949654),
+ (0x5450d14b, 0x41e74489),
+ (0x6052330d, 0x423633cf),
+ (0x75535300, 0x42955613),
+ (0x67d3bac5, 0x425fd1fa),
+ (0x43d4a2bb, 0x40c1c32f),
+ (0x57d628d5, 0x4207249d),
+ (0x75d75209, 0x4296c28e),
+ (0x65d79618, 0x4254cd54),
+ (0x665934fc, 0x42579ac8),
+ (0x4d5a2f36, 0x4199fc7c),
+ (0x4e5b5517, 0x41a51e5d),
+ (0x51dba523, 0x41cbf23d),
+ (0x51dd3656, 0x41cc00cd),
+ (0x3f3b4d9e, 0x3f0c9047),
+ (0x745f3fc3, 0x4292ac65),
+ (0x43dfe6db, 0x40c36926),
+ (0x4ae113c3, 0x417d04b7),
+ (0x77e22cfa, 0x429c674e),
+ (0xbf6399e7, 0xc00cb9a1),
+ (0x486401b0, 0x4145c60b),
+ (0x76e50770, 0x4299a7f1),
+ (0x53661704, 0x41dcf415),
+ (0x51e6f204, 0x41cc58fc),
+ (0x62e81393, 0x4244761b),
+ (0x5f691614, 0x42311213),
+ (0x68ea445b, 0x4265c51f),
+ (0x4beb7a84, 0x4189f603),
+ (0x7b6b9715, 0x42a6306d),
+ (0x606d0172, 0x4236aeb7),
+ (0x6a6d8093, 0x426e2483),
+ (0x4aef46b9, 0x417dff4a),
+ (0x40cff51e, 0x4000f145),
+ (0x4df0dd4c, 0x41a05296),
+ (0x6771f0ee, 0x425d94c7),
+ (0x5e736a8a, 0x422bb2ea),
+ (0x47f3e7ff, 0x413bc30a),
+ (0x58751ddd, 0x420a74a6),
+ (0x6675ba2a, 0x4258191d),
+ (0x75f6c0f1, 0x42970853),
+ (0x6d77f9a2, 0x427ef365),
+ (0x7278e907, 0x428d588a),
+ (0x62f99bda, 0x4244c0af),
+ (0x5b7b5415, 0x421b30f9),
+ (0x4cfbe7fe, 0x41959740),
+ (0x3f79f629, 0x3f2e6894),
+ (0x4ffe543f, 0x41b6f03f),
+ (0x4efebc08, 0x41abdc61),
+ (0x70fffac8, 0x42893e34), /* C[128] */
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var yf : flt32 = std.flt32frombits(y)
+ var rf = math.log1p(xf)
+ testr.check(c, rf == yf,
+ "log1p(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 log1p02 = {c+ var inputs : (uint64, uint64)[:] = [
+ (0x0000000000000000, 0x0000000000000000),
+ (0x6900000000000002, 0x407c765cf8301757), /* C[0] */
+ (0x6a30133d035de442, 0x407d4927a622132e), /* C[1] */
+ (0x545031c8b889228e, 0x406c3f4c487ae367), /* C[2] */
+ (0xbfdf2f686343ab7f, 0xbfe56047e9fe5abe), /* ... */
+ (0x52b08a1d750f04ff, 0x4069ff462f9f1cbc),
+ (0x4b309bc3b92e2dfc, 0x405f3371b7645f6c),
+ (0x5580cc2d69e348b5, 0x406de5e6ca65b1cd),
+ (0x4fe0e879a2736200, 0x406619d8f80fa200),
+ (0x4ad0fbf9d10cd8a2, 0x405e2ab53092727b),
+ (0x4ae113c3a6bc3a99, 0x405e576b1bac2f98),
+ (0x49113e8d334802ab, 0x4059518f80c8b520),
+ (0x54815a632f7ffc46, 0x406c840d22b38727),
+ (0x482175028dd2b016, 0x4056b8ec87c62e57),
+ (0x46719e36e5f0da58, 0x40520bc0c4781e1f),
+ (0x4651b338cfcc7ad5, 0x4051b353db2b8859),
+ (0x41a1dbe267002305, 0x4032d32be3bb9cb9),
+ (0x5811f2beb37f6007, 0x4070bab72614d157),
+ (0x7b922a8c7d7b8896, 0x4084ab1d75ac1cd5),
+ (0x74a2309e5aa7b2e9, 0x4082439c5e5fec3b),
+ (0x40d268e6c4ad9588, 0x4023b05609c9a18b),
+ (0x7f7285f6b3e07dfa, 0x4086031261f39110),
+ (0x54f29bb41ec2d2b9, 0x406d218d0b1c6484),
+ (0x7eb2bb6da1727619, 0x4085c09e8f1dadd9),
+ (0x71d2e2eb43d8256d, 0x40814a60e0b279cb),
+ (0x74630f734df0dd4c, 0x40822dcdd6310ba8),
+ (0x6fe31066fa71ca7a, 0x40809e8d868734a7),
+ (0x50e345290de2dd6c, 0x406780ec610f0c26),
+ (0x5e736a8abf6dda11, 0x40752730815a6166),
+ (0x6903781bd44c88f2, 0x407c7980c8a083ab),
+ (0x78b39b662df2692e, 0x4083aca5c3801c1c),
+ (0x6f93b5c49bb1bad3, 0x40808317f139e8e2),
+ (0x6133e925e6bf8a12, 0x40770f91495b0b5f),
+ (0x5a63fa758c92398a, 0x407256c5e940a2db),
+ (0x6e042c14d7b5903c, 0x407ff14c8c3b9be2),
+ (0x5f7445e6b0eee22e, 0x4075d9537ce0389b),
+ (0x4f646b7b7f44f304, 0x40656e70d305a121),
+ (0x6834802690008002, 0x407bea2785dd467d),
+ (0x43d4a2bb177f4778, 0x40459d6223a7d41b),
+ (0x6634cc7aab4228cc, 0x407a877e7a78169b),
+ (0x5894d0947b2155b9, 0x407115cf1b30a32d),
+ (0x5624f4404ffe543f, 0x406ecac8a9c25223),
+ (0x58751ddd6052330d, 0x4070ffdbd350379f),
+ (0x7f75348c0f5498fb, 0x408604275047d724),
+ (0x7265517fa66418ae, 0x40817d410871254b),
+ (0x7515838bfa241f5a, 0x40826bc50abc4fa5),
+ (0x45459662ae3323f5, 0x404d9bc7c072e470),
+ (0x7635cba6522da14e, 0x4082cfafdb79820f),
+ (0x5c35d8c190aea62e, 0x407399d2e3f6c83c),
+ (0x4555f61f1c0c5e59, 0x404df6b397ade3ba),
+ (0x57d628d547f3e7ff, 0x407091b9f622cd4e),
+ (0x78463abad2eed0c7, 0x408386d5e3383f6b),
+ (0x4b76524c54f71c43, 0x405ff7cf887ecb05),
+ (0x524675375aaefb86, 0x40696dcc32c58d99),
+ (0x68569f0d424db12a, 0x407c01e8fcc54a3b),
+ (0x7306c0360b296d7a, 0x4081b539dff3fe2b),
+ (0x4ad6e0cb2febd13d, 0x405e3dc5d9b757f7),
+ (0x51e6f204333bad9e, 0x4068e9668d521be8),
+ (0x4097178f382d0b5a, 0x401d32395ff3dc65),
+ (0x7f373de46ebc9a27, 0x4085eeb4db53f788),
+ (0x66d76e22a34e58cd, 0x407af84dc23a7f99),
+ (0x541784f5f4f77091, 0x406bf2841f114315),
+ (0x627797a59dd57660, 0x4077f016d94fac5c),
+ (0x7ab7bb38889422cf, 0x40845f9ed64c4859),
+ (0x4fc7e75368981c93, 0x4065f890d0df7af1),
+ (0x75d8018df3774adb, 0x4082af304ec36d3f),
+ (0x62e813931fa46ee3, 0x40783e0bf5cda73d),
+ (0x75984e941fdcceb2, 0x4082991b8dfdda1c),
+ (0x4f08632da9f79134, 0x4064ef09d28eeda0),
+ (0x4ee87d4f36b08d11, 0x4064c2cf83f46cf4),
+ (0x46a8931cf6708831, 0x4052a622d2cf5019),
+ (0x5c38b8d6d96e1ffe, 0x40739bcd56c393a8),
+ (0x7278e907f500a3f8, 0x4081840b7f2ad445),
+ (0x73f90e51c4166808, 0x4082092d01e05c94),
+ (0x7a7927863ba4c8ff, 0x408449e7d7ea5790),
+ (0x504943add8ae6c70, 0x4066abc8767dbc43),
+ (0x481953e382eee069, 0x4056a4615edb9861),
+ (0x6b8984577920fba9, 0x407e397207c69605),
+ (0x71799aa6bce3b976, 0x40812b8ab6921ace),
+ (0x6019b68e88bc2e04, 0x40764c0873609927),
+ (0x4bf9d05bd54aa5c3, 0x4060b200ae465dfa),
+ (0x753a038d35905a30, 0x4082786127c0c498),
+ (0x617a151591b1b79f, 0x4077403fbb1aa8c9),
+ (0x503a33bad761d1c1, 0x406696c4be84ff41),
+ (0x4d2a6cc069c8de6f, 0x4062582f440da781),
+ (0x6c1a8cc212dcb6e2, 0x407e9de4bd07d597),
+ (0x6b5a9349dd94a791, 0x407e18d31a0f5866),
+ (0x5c1abf7963c10a0c, 0x407386e1b0c37460),
+ (0x4d2ad7c63e651a60, 0x406258afdaa5e19a),
+ (0x407af5c43978f008, 0x401846ebf755d962),
+ (0x73cb2d9a04724837, 0x4081f930d140b8b3),
+ (0x48cb418c32cf8b91, 0x4058910d56870021),
+ (0x720b6590a8e3f54c, 0x40815dfd6485181b),
+ (0x46eb81a434b9f620, 0x40535ecb7294cc14),
+ (0x7b6b9715161c2442, 0x40849dd29d738baf),
+ (0x4eebb8fcb686a19f, 0x4064c6c75b7c96b0),
+ (0x4cfbe7fe8242657f, 0x4062176353551bf1),
+ (0x609bff020bdc4079, 0x4076a61dec83f7c3),
+ (0x7b9c14acd61754cb, 0x4084ae996873a75c),
+ (0x4f1c3cfb0a12764c, 0x406509e91a2092fc),
+ (0x6cfc6da63a428918, 0x407f3a4094ee1891),
+ (0x74ec84ae3e4aeb42, 0x40825d639229d899),
+ (0x56fcaa4913caa832, 0x406ff52903fd288b),
+ (0x623cc68a255eedf9, 0x4077c6e7cd2766f5),
+ (0x4bfcde0e60de47ad, 0x4060b5948be5f1a9),
+ (0x7b0cf7a77eff4d34, 0x40847cf0fbe4ff93),
+ (0x44bd10aeecadc14d, 0x404aa358796a4dca),
+ (0x7f7d456ec4df8ff2, 0x408606bb7f61f92f),
+ (0x7b8d69195829f3f8, 0x4084a96c99b76f20),
+ (0x7f6d78ed8101a55a, 0x4086013df53b153f),
+ (0x667dad261da3f8f5, 0x407ab98af553a20c),
+ (0x42cdc8a5274640b7, 0x403fd020a35da73e),
+ (0x5a6de65a705a667d, 0x40725d396e1dd345),
+ (0x795dfa770c749b2c, 0x4083e77ef991f95a),
+ (0x488e1f1c6db084f8, 0x4057e601112df1ff),
+ (0x6cce3dc1ca42a083, 0x407f19f86830564f),
+ (0x697e599a1e4afd5f, 0x407cce3d2d79ec2e),
+ (0x491e82abf353e79c, 0x40597613f67eb37a),
+ (0x722ea4ce7f96d2ff, 0x408169f9e9352797),
+ (0x4efebc0847e20599, 0x4064e04285d5e27e),
+ (0x544eeb58043197b1, 0x406c3dcfe12bf3bb),
+ (0x46ff046b11b6e7e0, 0x405392d817dda2d2),
+ (0x6caf1807e61b6f03, 0x407f043c0805dba1),
+ (0x50ef387e0fc4a5e1, 0x4067905d34c36a51),
+ (0x4cbf6bd74d8f3edf, 0x4061c276224e5d71),
+ (0x4f5f7177b7abc69e, 0x40656612da92283f),
+ (0x629fabc84d5a2f36, 0x40780afb4bbfbe4c),
+ (0x418fce8b9cdac427, 0x40320409995dc1e7),
+ (0x46cfe9b12985df6e, 0x40530f94e5ddae5e),
+ (0x70fffac8f637436f, 0x408100f58e19ab8f), /* C[128] */
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt64 = std.flt64frombits(x)
+ var yf : flt64 = std.flt64frombits(y)
+ var rf = math.log1p(xf)
+ testr.check(c, rf == yf,
+ "log1p(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 log1p03 = {c+ /*
+ As with log, there is some accepted error in log1p.
+ */
+
+ var inputs : (uint32, uint32, uint32)[:] = [
+ (0x49c68d15, 0x4164d4d5, 0x4164d4d4),
+ (0x3d86912c, 0x3d8254a9, 0x3d8254a8),
+ (0x3dd7210e, 0x3dcc905b, 0x3dcc905c),
+ (0x3d986e71, 0x3d93067e, 0x3d93067f),
+ (0xbe1eefcb, 0xbe2cb799, 0xbe2cb798),
+ (0x3e057287, 0x3dfae18d, 0x3dfae18c),
+ (0x424d8fe0, 0x407d5bc1, 0x407d5bc2),
+ (0xb95cb5e9, 0xb95cbbdb, 0xb95cbbdc),
+ (0x3de66745, 0x3dda56fd, 0x3dda56fc),
+ ][:]
+
+ 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.log1p(xf)
+ if rf != ypf && rf != yaf
+ testr.fail(c, "log1p(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 log1p04 = {c+ /*
+ As with log, there is some accepted error in log1p.
+ */
+
+ var inputs : (uint64, uint64, uint64)[:] = [
+ (0xbf8d2fb5e91b21dc, 0xbf8d65764edb0cd6, 0xbf8d65764edb0cd5),
+ (0x3fc855690a4a67e1, 0x3fc64708ed6e9abb, 0x3fc64708ed6e9aba),
+ (0xbfafb59aa6bb5f14, 0xbfb05dee438595dd, 0xbfb05dee438595de),
+ (0x3f896e0154c1be37, 0x3f8945eb78442aa1, 0x3f8945eb78442aa0),
+ (0x3fb09ef0bcfe6932, 0x3fb01a8404c5051a, 0x3fb01a8404c50519),
+ (0x3fa071dec13893e8, 0x3fa02fad06dc3334, 0x3fa02fad06dc3335),
+ (0x4000d2445e953eb4, 0x3ff21dbfa8f28f5d, 0x3ff21dbfa8f28f5c),
+ (0xbfe37c5eda902f8d ,0xbfee0b40d5f061d7, 0xbfee0b40d5f061d6),
+ (0x400dd2fe516cced3, 0x3ff8db2a8f466eeb, 0x3ff8db2a8f466eea),
+ (0xbfb5d9612ba5b9bf, 0xbfb6d6962ad7508b, 0xbfb6d6962ad7508c),
+ (0x40c512345c72e7f9, 0x4022929892b71a96, 0x4022929892b71a95),
+ (0x47409b795894785f, 0x405448ab9f468935, 0x405448ab9f468936),
+ ][:]
+
+ 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.log1p(xf)
+ if rf != ypf && rf != yaf
+ testr.fail(c, "log1p(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)
+ ;;
+ ;;
+}