shithub: mc

ref: 0d7d70e5b1a4b5c7ae618d6d04fb30848b18c08f
dir: /lib/math/log-impl.myr/

View raw version
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)
}