ref: 6a113c23192c0877185aadddbc58b943ed68dd7b
dir: /lib/math/log-impl.myr/
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)
/* Constants from [Tan90], note that [128] contains accurate log(2) */
pkglocal const accurate_logs32 : (uint32, uint32)[129]
pkglocal const accurate_logs64 : (uint64, uint64)[129]
;;
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)[:]
Ai : @u[:]
/* For procedure 2 */
Bi : @u[:]
Mtruncmask : @u
;;
/* Accurate representations for log(1 + j/2^7), all j */
const accurate_logs32 = [
(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 */
]
const accurate_logs64 = [
(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),
]
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_logs32[:],
.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 = accurate_logs64[:],
.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)
}