ref: 690ab7aa8e6c56a4ae8a2bd5bce58d1df681234e
dir: /lib/math/log-overkill.myr/
use std
use "fpmath"
use "impls"
use "util"
/*
This is an implementation of log(x) using the following idea, based on [Tan90]
First, reduce to 2^e · xs, with xs ∈ [1, 2).
xs = F1 + f1, with
F1 = 1 + j1/2^5
j1 ∈ {1, 2, …, 2^5 - 1}
f1 ∈ [0, 2^-5)
log(xs) = log(F1) + log(1 + f1/F1)
1 + f1/F1 = F2 + f2, with
F2 = 1 + j1/2^10
j2 ∈ {1, 2, …, 2^5 - 1}
f2 ∈ [0, 2^-10)
log(xs) = log(F1) + log(F2) + log(1 + f2/F2)
…
log(xs) = log(F1) + log(F2) + log(F3) + log(F4) + log(1 + f4/F4)
And f4/F4 < 2^-20, so we can get 100 bits of precision using a
degree 5 polynomial.
The specific choice of using 4 tables, each with 2^5 entries, may
be improvable. It's a trade-off between storage for the tables and
the number of floating point operations to chain the results
together.
*/
pkg math =
pkglocal const logoverkill32 : (x : flt32 -> (flt32, flt32))
pkglocal const logoverkill64 : (x : flt64 -> (flt64, flt64))
;;
/*
Ci is a table such that, for Ci[j] = (L1, L2, I1, I2),
L1, L2 are log(1 + j·2^-(5i))
I1, I2 are 1/(1 + j·2^-(5i))
*/
const C1 : (uint64, uint64, uint64, uint64)[32] = [
(000000000000000000, 000000000000000000, 0x3ff0000000000000, 000000000000000000), /* j = 0 */
(0x3f9f829b0e783300, 0x3c333e3f04f1ef23, 0x3fef07c1f07c1f08, 0xbc7f07c1f07c1f08), /* j = 1 */
(0x3faf0a30c01162a6, 0x3c485f325c5bbacd, 0x3fee1e1e1e1e1e1e, 0x3c6e1e1e1e1e1e1e), /* j = 2 */
(0x3fb6f0d28ae56b4c, 0xbc5906d99184b992, 0x3fed41d41d41d41d, 0x3c80750750750750), /* j = 3 */
(0x3fbe27076e2af2e6, 0xbc361578001e0162, 0x3fec71c71c71c71c, 0x3c8c71c71c71c71c), /* j = 4 */
(0x3fc29552f81ff523, 0x3c6301771c407dbf, 0x3febacf914c1bad0, 0xbc8bacf914c1bad0), /* j = 5 */
(0x3fc5ff3070a793d4, 0xbc5bc60efafc6f6e, 0x3feaf286bca1af28, 0x3c8af286bca1af28), /* j = 6 */
(0x3fc9525a9cf456b4, 0x3c6d904c1d4e2e26, 0x3fea41a41a41a41a, 0x3c80690690690690), /* j = 7 */
(0x3fcc8ff7c79a9a22, 0xbc64f689f8434012, 0x3fe999999999999a, 0xbc8999999999999a), /* j = 8 */
(0x3fcfb9186d5e3e2b, 0xbc6caaae64f21acb, 0x3fe8f9c18f9c18fa, 0xbc7f3831f3831f38), /* j = 9 */
(0x3fd1675cababa60e, 0x3c2ce63eab883717, 0x3fe8618618618618, 0x3c88618618618618), /* j = 10 */
(0x3fd2e8e2bae11d31, 0xbc78f4cdb95ebdf9, 0x3fe7d05f417d05f4, 0x3c67d05f417d05f4), /* j = 11 */
(0x3fd4618bc21c5ec2, 0x3c7f42decdeccf1d, 0x3fe745d1745d1746, 0xbc7745d1745d1746), /* j = 12 */
(0x3fd5d1bdbf5809ca, 0x3c74236383dc7fe1, 0x3fe6c16c16c16c17, 0xbc7f49f49f49f49f), /* j = 13 */
(0x3fd739d7f6bbd007, 0xbc78c76ceb014b04, 0x3fe642c8590b2164, 0x3c7642c8590b2164), /* j = 14 */
(0x3fd89a3386c1425b, 0xbc729639dfbbf0fb, 0x3fe5c9882b931057, 0x3c7310572620ae4c), /* j = 15 */
(0x3fd9f323ecbf984c, 0xbc4a92e513217f5c, 0x3fe5555555555555, 0x3c85555555555555), /* j = 16 */
(0x3fdb44f77bcc8f63, 0xbc7cd04495459c78, 0x3fe4e5e0a72f0539, 0x3c8e0a72f0539783), /* j = 17 */
(0x3fdc8ff7c79a9a22, 0xbc74f689f8434012, 0x3fe47ae147ae147b, 0xbc6eb851eb851eb8), /* j = 18 */
(0x3fddd46a04c1c4a1, 0xbc70467656d8b892, 0x3fe4141414141414, 0x3c64141414141414), /* j = 19 */
(0x3fdf128f5faf06ed, 0xbc7328df13bb38c3, 0x3fe3b13b13b13b14, 0xbc83b13b13b13b14), /* j = 20 */
(0x3fe02552a5a5d0ff, 0xbc7cb1cb51408c00, 0x3fe3521cfb2b78c1, 0x3c7a90e7d95bc60a), /* j = 21 */
(0x3fe0be72e4252a83, 0xbc8259da11330801, 0x3fe2f684bda12f68, 0x3c82f684bda12f68), /* j = 22 */
(0x3fe154c3d2f4d5ea, 0xbc859c33171a6876, 0x3fe29e4129e4129e, 0x3c804a7904a7904a), /* j = 23 */
(0x3fe1e85f5e7040d0, 0x3c7ef62cd2f9f1e3, 0x3fe2492492492492, 0x3c82492492492492), /* j = 24 */
(0x3fe2795e1289b11b, 0xbc6487c0c246978e, 0x3fe1f7047dc11f70, 0x3c81f7047dc11f70), /* j = 25 */
(0x3fe307d7334f10be, 0x3c6fb590a1f566da, 0x3fe1a7b9611a7b96, 0x3c61a7b9611a7b96), /* j = 26 */
(0x3fe393e0d3562a1a, 0xbc858eef67f2483a, 0x3fe15b1e5f75270d, 0x3c415b1e5f75270d), /* j = 27 */
(0x3fe41d8fe84672ae, 0x3c89192f30bd1806, 0x3fe1111111111111, 0x3c61111111111111), /* j = 28 */
(0x3fe4a4f85db03ebb, 0x3c313dfa3d3761b6, 0x3fe0c9714fbcda3b, 0xbc7f79b47582192e), /* j = 29 */
(0x3fe52a2d265bc5ab, 0xbc61883750ea4d0a, 0x3fe0842108421084, 0x3c70842108421084), /* j = 30 */
(0x3fe5ad404c359f2d, 0xbc435955683f7196, 0x3fe0410410410410, 0x3c80410410410410), /* j = 31 */
]
const C2 : (uint64, uint64, uint64, uint64)[32] = [
(000000000000000000, 000000000000000000, 0x3ff0000000000000, 000000000000000000), /* j = 0 */
(0x3f4ffc00aa8ab110, 0xbbe0fecbeb9b6cdb, 0x3feff801ff801ff8, 0x3c2ff801ff801ff8), /* j = 1 */
(0x3f5ff802a9ab10e6, 0x3bfe29e3a153e3b2, 0x3feff007fc01ff00, 0x3c8ff007fc01ff00), /* j = 2 */
(0x3f67f7047d7983da, 0x3c0a275a19204e80, 0x3fefe811f28a186e, 0xbc849093915301bf), /* j = 3 */
(0x3f6ff00aa2b10bc0, 0x3c02821ad5a6d353, 0x3fefe01fe01fe020, 0xbc6fe01fe01fe020), /* j = 4 */
(0x3f73f38a60f06489, 0x3c1693c937494046, 0x3fefd831c1cdbed1, 0x3c8e89d3b75ace7e), /* j = 5 */
(0x3f77ee11ebd82e94, 0xbc161e96e2fc5d90, 0x3fefd04794a10e6a, 0x3c881bd63ea20ced), /* j = 6 */
(0x3f7be79c70058ec9, 0xbbd964fefef02b62, 0x3fefc86155aa1659, 0xbc6b8fc468497f61), /* j = 7 */
(0x3f7fe02a6b106789, 0xbbce44b7e3711ebf, 0x3fefc07f01fc07f0, 0x3c6fc07f01fc07f0), /* j = 8 */
(0x3f81ebde2d1997e6, 0xbbfffa46e1b2ec81, 0x3fefb8a096acfacc, 0xbc82962e18495af3), /* j = 9 */
(0x3f83e7295d25a7d9, 0xbbeff29a11443a06, 0x3fefb0c610d5e939, 0xbc5cb8337f41db5c), /* j = 10 */
(0x3f85e1f703ecbe50, 0x3c23eb0bb43693b9, 0x3fefa8ef6d92aca5, 0x3c7cd0c1eaba7f22), /* j = 11 */
(0x3f87dc475f810a77, 0xbc116d7687d3df21, 0x3fefa11caa01fa12, 0xbc7aaff02f71aaff), /* j = 12 */
(0x3f89d61aadc6bd8d, 0xbc239b097b525947, 0x3fef994dc3455e8d, 0xbc82546d9bc5bd59), /* j = 13 */
(0x3f8bcf712c74384c, 0xbc1f6842688f499a, 0x3fef9182b6813baf, 0x3c6b210c54d70f4a), /* j = 14 */
(0x3f8dc84b19123815, 0xbc25f0e2d267d821, 0x3fef89bb80dcc421, 0xbc8e7da8c7156f9d), /* j = 15 */
(0x3f8fc0a8b0fc03e4, 0xbc183092c59642a1, 0x3fef81f81f81f820, 0xbc8f81f81f81f820), /* j = 16 */
(0x3f90dc4518afcc88, 0xbc379c0189fdfe78, 0x3fef7a388f9da20f, 0x3c7f99b2c82d3fb1), /* j = 17 */
(0x3f91d7f7eb9eebe7, 0xbc2d41fe63d2dbf9, 0x3fef727cce5f530a, 0x3c84643cedd1cfd9), /* j = 18 */
(0x3f92d36cefb557c3, 0xbc132fa3e4f20cf7, 0x3fef6ac4d8f95f7a, 0x3c8e8ed9770a8dde), /* j = 19 */
(0x3f93cea44346a575, 0xbc10cb5a902b3a1c, 0x3fef6310aca0dbb5, 0x3c8d2e19807d8c43), /* j = 20 */
(0x3f94c99e04901ded, 0xbc2cd10505ada0d6, 0x3fef5b60468d989f, 0xbc805a26b4cad717), /* j = 21 */
(0x3f95c45a51b8d389, 0xbc3b10b6c3ec21b4, 0x3fef53b3a3fa204e, 0x3c8450467c5430f3), /* j = 22 */
(0x3f96bed948d1b7d1, 0xbc1058290fde6de1, 0x3fef4c0ac223b2bc, 0x3c815c2df7afcd24), /* j = 23 */
(0x3f97b91b07d5b11b, 0xbc35b602ace3a510, 0x3fef44659e4a4271, 0x3c85fc17734c36b8), /* j = 24 */
(0x3f98b31faca9b00e, 0x3c1d5a46da6f6772, 0x3fef3cc435b0713c, 0x3c81d0a7e69ea094), /* j = 25 */
(0x3f99ace7551cc514, 0x3c33409c1df8167f, 0x3fef3526859b8cec, 0x3c2f3526859b8cec), /* j = 26 */
(0x3f9aa6721ee835aa, 0xbc24a3a50b6c5621, 0x3fef2d8c8b538c0f, 0xbc88a987ac35964a), /* j = 27 */
(0x3f9b9fc027af9198, 0xbbf0ae69229dc868, 0x3fef25f644230ab5, 0x3c594ed8175c78b3), /* j = 28 */
(0x3f9c98d18d00c814, 0xbc250589df0f25bf, 0x3fef1e63ad57473c, 0xbc8bf54d8dbc6a00), /* j = 29 */
(0x3f9d91a66c543cc4, 0xbc1d34e608cbdaab, 0x3fef16d4c4401f17, 0xbc759ddff074959e), /* j = 30 */
(0x3f9e8a3ee30cdcac, 0x3c07086b1c00b395, 0x3fef0f4986300ba6, 0xbc811b6b7ee8766a), /* j = 31 */
]
const C3 : (uint64, uint64, uint64, uint64)[32] = [
(000000000000000000, 000000000000000000, 0x3ff0000000000000, 000000000000000000), /* j = 0 */
(0x3effffe0002aaa6b, 0xbb953bbbe6661d42, 0x3fefffc0007fff00, 0x3c2fffc0007fff00), /* j = 1 */
(0x3f0fffc000aaa8ab, 0xbba3bbc110fec82c, 0x3fefff8001fff800, 0x3c6fff8001fff800), /* j = 2 */
(0x3f17ffb8011ffaf0, 0x3b984c534f3d9b6a, 0x3fefff40047fe501, 0xbc8780f2fa4e222b), /* j = 3 */
(0x3f1fff8002aa9aab, 0x3b910e6678af0afc, 0x3fefff0007ffc002, 0xbbdfff0007ffc002), /* j = 4 */
(0x3f23ff9c029a9723, 0x3bc1b965303b23b1, 0x3feffec00c7f8305, 0xbc6e30d217cb1211), /* j = 5 */
(0x3f27ff70047fd782, 0xbbced098a5c0aff0, 0x3feffe8011ff280a, 0x3c6f8685b1bbab34), /* j = 6 */
(0x3f2bff3c07250a51, 0xbbc89dd6d6bad8c1, 0x3feffe40187ea913, 0xbc7f8346d2208239), /* j = 7 */
(0x3f2fff000aaa2ab1, 0x3ba0bbc04dc4e3dc, 0x3feffe001ffe0020, 0xbc2ffe001ffe0020), /* j = 8 */
(0x3f31ff5e07979982, 0xbbce0e704817ebcd, 0x3feffdc0287d2733, 0x3c7f32ce6d7c4d43), /* j = 9 */
(0x3f33ff380a6a0e74, 0x3bdb81fcb95bc1fe, 0x3feffd8031fc184e, 0x3c69e5fa087756ad), /* j = 10 */
(0x3f35ff0e0ddc70a1, 0x3bacf6f3d97a3c05, 0x3feffd403c7acd72, 0x3c860b1b0bacff22), /* j = 11 */
(0x3f37fee011febc18, 0x3bd2b9bcf5d3f323, 0x3feffd0047f940a2, 0xbc5e5d274451985a), /* j = 12 */
(0x3f39feae16e0ec8b, 0xbbb6137aceeb34b1, 0x3feffcc054776bdf, 0x3c56b1b1f3ed39e8), /* j = 13 */
(0x3f3bfe781c92fd4a, 0xbbc4ed10713cc126, 0x3feffc8061f5492c, 0xbc19fd284f974b74), /* j = 14 */
(0x3f3dfe3e2324e946, 0x3bc0916462dd5deb, 0x3feffc407072d28b, 0x3c84eb0c748a57ca), /* j = 15 */
(0x3f3ffe002aa6ab11, 0x3b999e2b62cc632d, 0x3feffc007ff00200, 0xbc7ffc007ff00200), /* j = 16 */
(0x3f40fedf19941e6e, 0xbbb194c2e0aa6338, 0x3feffbc0906cd18c, 0x3c75b11e79f3cd9f), /* j = 17 */
(0x3f41febc1e5ccc3c, 0x3bdc657d895d3592, 0x3feffb80a1e93b34, 0xbc84d11299626e29), /* j = 18 */
(0x3f42fe9723b55bac, 0x3bd6cfb73e538464, 0x3feffb40b46538fa, 0xbc8d42a81b0bfc39), /* j = 19 */
(0x3f43fe7029a5c947, 0xbbd4d578bf46e36a, 0x3feffb00c7e0c4e1, 0x3c7e673fde054f2c), /* j = 20 */
(0x3f44fe4730361165, 0x3be400d77e93f2fd, 0x3feffac0dc5bd8ee, 0x3c8a38b2b2aeaf57), /* j = 21 */
(0x3f45fe1c376e3031, 0xbbd524eb8a5ae7f6, 0x3feffa80f1d66f25, 0xbc6a5778f73582ce), /* j = 22 */
(0x3f46fdef3f5621a3, 0xbbdf09d734886d52, 0x3feffa4108508189, 0xbc81a45478d24a37), /* j = 23 */
(0x3f47fdc047f5e185, 0xbbebfa5c57d202d3, 0x3feffa011fca0a1e, 0x3c6a5b0eed338657), /* j = 24 */
(0x3f48fd8f51556b70, 0xbbd0f4f2e08fd201, 0x3feff9c1384302e9, 0x3c8b9a1be68cf877), /* j = 25 */
(0x3f49fd5c5b7cbace, 0x3beb6ed49f17d42d, 0x3feff98151bb65ef, 0x3c82d92be315df8f), /* j = 26 */
(0x3f4afd276673cada, 0x3bd3222545da594f, 0x3feff9416c332d34, 0x3c8dbbba66ae573a), /* j = 27 */
(0x3f4bfcf07242969d, 0x3bc5db4d2b3efe1c, 0x3feff90187aa52be, 0xbc698a69b8df8f19), /* j = 28 */
(0x3f4cfcb77ef118f1, 0x3becc55406f300fb, 0x3feff8c1a420d091, 0xbc8032d47bdbf02c), /* j = 29 */
(0x3f4dfc7c8c874c82, 0xbbe863e9d57a176f, 0x3feff881c196a0b2, 0x3c858cf2f70e18b2), /* j = 30 */
(0x3f4efc3f9b0d2bc8, 0x3bd1e8e8a5f5b8b7, 0x3feff841e00bbd28, 0x3c782227ba60dc8b), /* j = 31 */
]
const C4 : (uint64, uint64, uint64, uint64)[32] = [
(000000000000000000, 000000000000000000, 0x3ff0000000000000, 000000000000000000), /* j = 0 */
(0x3eafffff00000aab, 0xbb45755553bbbbd1, 0x3feffffe00002000, 0xbc2ffffe00002000), /* j = 1 */
(0x3ebffffe00002aab, 0xbb5655553bbbbe66, 0x3feffffc00008000, 0xbc5ffffc00008000), /* j = 2 */
(0x3ec7fffdc0004800, 0xbb343ffcf666dfe6, 0x3feffffa00012000, 0xbc7afffaf000f300), /* j = 3 */
(0x3ecffffc0000aaab, 0xbb6d5553bbbc1111, 0x3feffff800020000, 0xbc8ffff800020000), /* j = 4 */
(0x3ed3fffce000a6ab, 0xbb7f1952e455f818, 0x3feffff600031fff, 0x3c4801387f9e581f), /* j = 5 */
(0x3ed7fffb80012000, 0xbb743ff9ecceb2cc, 0x3feffff400047ffe, 0x3c8400287ff0d006), /* j = 6 */
(0x3edbfff9e001c955, 0xbb702e9d89490dc5, 0x3feffff200061ffd, 0x3c84804b07df2c8e), /* j = 7 */
(0x3edffff80002aaaa, 0xbb75553bbbc66662, 0x3feffff00007fffc, 0x3baffff00007fffc), /* j = 8 */
(0x3ee1fffaf001e5ff, 0x3b797c2e21b72cff, 0x3fefffee000a1ffa, 0x3c8380cd078cabc1), /* j = 9 */
(0x3ee3fff9c0029aa9, 0x3b8c8ad1ba965260, 0x3fefffec000c7ff8, 0x3c780270fe7960f4), /* j = 10 */
(0x3ee5fff870037754, 0xbb8d0c6bc1b51bdd, 0x3fefffea000f1ff6, 0xbc897e36793a8ca8), /* j = 11 */
(0x3ee7fff700047ffd, 0x3b8e006132f6735a, 0x3fefffe80011fff3, 0xbc8ffd7801e5fe94), /* j = 12 */
(0x3ee9fff57005b8a7, 0x3b771277672a8835, 0x3fefffe600151fef, 0xbc74f906f5aa5866), /* j = 13 */
(0x3eebfff3c0072551, 0xbb86c9d894dd7427, 0x3fefffe400187feb, 0xbc8bfb4f841a6c69), /* j = 14 */
(0x3eedfff1f008c9fa, 0xbb7701aebecf7ae4, 0x3fefffe2001c1fe6, 0xbc8779d1fdcb2212), /* j = 15 */
(0x3eeffff0000aaaa3, 0xbb8553bbbd110fec, 0x3fefffe0001fffe0, 0x3befffe0001fffe0), /* j = 16 */
(0x3ef0fff6f80665a6, 0xbb9b95400571451d, 0x3fefffde00241fda, 0xbc8875ce02d51cfe), /* j = 17 */
(0x3ef1fff5e00797fa, 0xbb9a0e8ef2f395cc, 0x3fefffdc00287fd2, 0x3c8c0cd071958038), /* j = 18 */
(0x3ef2fff4b808ee4d, 0x3b984638f069f71f, 0x3fefffda002d1fca, 0x3c8a8fe8751bf4ef), /* j = 19 */
(0x3ef3fff3800a6aa1, 0xbb794b915f81751a, 0x3fefffd80031ffc2, 0xbc8fec781869e17c), /* j = 20 */
(0x3ef4fff2380c0ef4, 0x3b80a43b5348b6b9, 0x3fefffd600371fb8, 0xbc8668429728999b), /* j = 21 */
(0x3ef5fff0e00ddd47, 0x3b624a1f136cc09c, 0x3fefffd4003c7fad, 0xbc77c6cf4ea2f3e0), /* j = 22 */
(0x3ef6ffef780fd79a, 0xbb9a716c4210cdd0, 0x3fefffd200421fa1, 0xbc5aeeb948d5a74d), /* j = 23 */
(0x3ef7ffee0011ffec, 0xbb8ff3d9a8c98613, 0x3fefffd00047ff94, 0x3c143fe1a02d8fbc), /* j = 24 */
(0x3ef8ffec7814583d, 0x3b9f7bc8a4e1db73, 0x3fefffce004e1f86, 0xbc6141450a04205a), /* j = 25 */
(0x3ef9ffeae016e28f, 0xbb8cb889999dd735, 0x3fefffcc00547f77, 0xbc83c837daa53cb3), /* j = 26 */
(0x3efaffe93819a0e0, 0xbb9be60d8a0b5ba8, 0x3fefffca005b1f66, 0x3c7d81be350f0677), /* j = 27 */
(0x3efbffe7801c9530, 0xbb873b12aed4646c, 0x3fefffc80061ff55, 0xbc8fb4f8834d1a39), /* j = 28 */
(0x3efcffe5b81fc17f, 0x3b9fe950a87b2785, 0x3fefffc600691f41, 0x3c8dd655eb844520), /* j = 29 */
(0x3efdffe3e02327cf, 0xbb9bfd7604f796f2, 0x3fefffc400707f2d, 0x3c618b7f1a71ae6b), /* j = 30 */
(0x3efeffe1f826ca1d, 0xbb60ae99630e566e, 0x3fefffc200781f17, 0x3c80f0bb2d9557af), /* j = 31 */
]
const logoverkill32 = {x : flt32
var x64 : flt64 = (x : flt64)
var l64 : flt64 = log64(x64)
var y1 : flt32 = (l64 : flt32)
var y2 : flt32 = ((l64 - (y1 : flt64)) : flt32)
-> (y1, y2)
}
const logoverkill64 = {x : flt64
var xn, xe, xs
(xn, xe, xs) = std.flt64explode(x)
var xsf = std.flt64assem(false, 0, xs)
var xb = std.flt64bits(x)
var tn, te, ts
var t1, t2
/* Special cases */
if std.isnan(x)
-> (std.flt64frombits(0x7ff8000000000000), std.flt64frombits(0x7ff8000000000000))
elif xe <= -1023 && xs == 0
/* log (+/- 0) is -infinity */
-> (std.flt64frombits(0xfff0000000000000), std.flt64frombits(0xfff0000000000000))
elif xn
/* log(-anything) is NaN */
-> (std.flt64frombits(0x7ff8000000000000), std.flt64frombits(0x7ff8000000000000))
;;
if xe <= -1023
/*
We depend on being able to pick bits out of xs as if it were normal, so
normalize any subnormals.
*/
xe++
var check = 1 << 52
while xs & check == 0
xs <<= 1
xe--
;;
xsf = std.flt64assem(false, 0, xs)
;;
var shift = 0
var non_trivial = 0
var then_invert = false
var j1, F1, f1, logF1_hi, logF1_lo, F1_inv_hi, F1_inv_lo, fF1_hi, fF1_lo
/* F1 */
if xe == -1 && xs > 0x001d1ff05e41cfba
/*
If we reduced to [1, 2) unconditionally, then values of x like 0.999… =
2^-1 · 1.999… would cause subtractive cancellation; we'd compute
log(1.999…), then subtract out log(2) at the end. They'd agree on the
first n bits, and we'd lose n bits of precision.
This is only a problem for exponent -1, and for xs large enough;
outside that, the numbers are so different that we won't lose precision
by cancelling. But here, we compute 1/x, proceed (with exponent 0), and
flip all the signs at the end.
*/
xe = 0
var xinv_hi = 1.0 / x
var xinv_lo = fma64(-1.0 * xinv_hi, x, 1.0) / x
(tn, te, ts) = std.flt64explode(xinv_hi)
non_trivial = ((47 >= te) : uint64) * ((47 - te < 64) : uint64)
shift = non_trivial * ((47 - te) : uint64)
j1 = non_trivial * ((ts >> shift) & 0x1f)
var F1m1 = scale2((j1 : flt64), -5)
F1 = 1.0 + F1m1
var f1_hi, f1_lo
(f1_hi, f1_lo) = fast2sum(xinv_hi - F1, xinv_lo)
(logF1_hi, logF1_lo, F1_inv_hi, F1_inv_lo) = C1[j1]
/* Compute 1 + f1/F1 */
(fF1_hi, fF1_lo) = two_by_two64(f1_hi, std.flt64frombits(F1_inv_hi))
(fF1_lo, t1) = slow2sum(fF1_lo, f1_lo * std.flt64frombits(F1_inv_hi))
(fF1_lo, t2) = slow2sum(fF1_lo, f1_hi * std.flt64frombits(F1_inv_lo))
(fF1_hi, fF1_lo) = fast2sum(fF1_hi, fF1_lo + (t1 + t2))
then_invert = true
else
j1 = (xs & 0x000f800000000000) >> 47
F1 = std.flt64assem(false, 0, xs & 0xffff800000000000)
f1 = xsf - F1
(logF1_hi, logF1_lo, F1_inv_hi, F1_inv_lo) = C1[j1]
/* Compute 1 + f1/F1 */
(fF1_hi, fF1_lo) = two_by_two64(f1, std.flt64frombits(F1_inv_hi))
(fF1_lo, t1) = slow2sum(fF1_lo, f1 * std.flt64frombits(F1_inv_lo))
(fF1_hi, fF1_lo) = fast2sum(fF1_hi, fF1_lo)
fF1_lo += t1
;;
/* F2 */
(tn, te, ts) = std.flt64explode(fF1_hi)
non_trivial = ((42 >= te) : uint64) * ((42 - te < 64) : uint64)
shift = non_trivial * ((42 - te) : uint64)
var j2 = non_trivial * ((ts >> shift) & 0x1f)
var F2m1 = scale2((j2 : flt64), -10)
var F2 = 1.0 + F2m1
var f2_hi, f2_lo
(f2_hi, f2_lo) = fast2sum(fF1_hi - F2m1, fF1_lo)
var logF2_hi, logF2_lo, F2_inv_hi, F2_inv_lo
(logF2_hi, logF2_lo, F2_inv_hi, F2_inv_lo) = C2[j2]
/* Compute 1 + f2/F2 */
var fF2_hi, fF2_lo
(fF2_hi, fF2_lo) = two_by_two64(f2_hi, std.flt64frombits(F2_inv_hi))
(fF2_lo, t1) = slow2sum(fF2_lo, f2_lo * std.flt64frombits(F2_inv_hi))
(fF2_lo, t2) = slow2sum(fF2_lo, f2_hi * std.flt64frombits(F2_inv_lo))
(fF2_hi, fF2_lo) = fast2sum(fF2_hi, fF2_lo + (t1 + t2))
/* F3 (just like F2) */
(tn, te, ts) = std.flt64explode(fF2_hi)
non_trivial = ((37 >= te) : uint64) * ((37 - te < 64) : uint64)
shift = non_trivial * ((37 - te) : uint64)
var j3 = non_trivial * ((ts >> shift) & 0x1f)
var F3m1 = scale2((j3 : flt64), -15)
var F3 = 1.0 + F3m1
var f3_hi, f3_lo
(f3_hi, f3_lo) = fast2sum(fF2_hi - F3m1, fF2_lo)
var logF3_hi, logF3_lo, F3_inv_hi, F3_inv_lo
(logF3_hi, logF3_lo, F3_inv_hi, F3_inv_lo) = C3[j3]
/* Compute 1 + f3/F3 */
var fF3_hi, fF3_lo
(fF3_hi, fF3_lo) = two_by_two64(f3_hi, std.flt64frombits(F3_inv_hi))
(fF3_lo, t1) = slow2sum(fF3_lo, f3_lo * std.flt64frombits(F3_inv_hi))
(fF3_lo, t2) = slow2sum(fF3_lo, f3_hi * std.flt64frombits(F3_inv_lo))
(fF3_hi, fF3_lo) = fast2sum(fF3_hi, fF3_lo + (t1 + t2))
/* F4 (just like F2) */
(tn, te, ts) = std.flt64explode(fF3_hi)
non_trivial = ((32 >= te) : uint64) * ((32 - te < 64) : uint64)
shift = non_trivial * ((32 - te) : uint64)
var j4 = non_trivial * ((ts >> shift) & 0x1f)
var F4m1 = scale2((j4 : flt64), -20)
var F4 = 1.0 + F4m1
var f4_hi, f4_lo
(f4_hi, f4_lo) = fast2sum(fF3_hi - F4m1, fF3_lo)
var logF4_hi, logF4_lo, F4_inv_hi, F4_inv_lo
(logF4_hi, logF4_lo, F4_inv_hi, F4_inv_lo) = C4[j4]
/* Compute 1 + f4/F4 */
var fF4_hi, fF4_lo
(fF4_hi, fF4_lo) = two_by_two64(f4_hi, std.flt64frombits(F4_inv_hi))
(fF4_lo, t1) = slow2sum(fF4_lo, f4_lo * std.flt64frombits(F4_inv_hi))
(fF4_lo, t2) = slow2sum(fF4_lo, f4_hi * std.flt64frombits(F4_inv_lo))
(fF4_hi, fF4_lo) = fast2sum(fF4_hi, fF4_lo + (t1 + t2))
/*
L = log(1 + f4/F4); we'd like to use horner_polyu, but since we have
_hi and _lo, it becomes more complicated.
*/
var L_hi, L_lo
/* r = (1/5) · x */
(L_hi, L_lo) = hl_mult(std.flt64frombits(0x3fc999999999999a), std.flt64frombits(0xbc6999999999999a), fF4_hi, fF4_lo)
/* r = r - 1/4 */
(t1, t2) = fast2sum(std.flt64frombits(0xbfd0000000000000), L_lo)
(L_hi, L_lo) = fast2sum(t1, L_hi)
L_lo += t2
/* r = r · x */
(L_hi, L_lo) = hl_mult(L_hi, L_lo, fF4_hi, fF4_lo)
/* r = r + 1/3 */
(L_hi, L_lo) = hl_add(std.flt64frombits(0x3fd5555555555555), std.flt64frombits(0x3c75555555555555), L_hi, L_lo)
/* r = r · x */
(L_hi, L_lo) = hl_mult(L_hi, L_lo, fF4_hi, fF4_lo)
/* r = r - 1/2 */
(t1, t2) = fast2sum(std.flt64frombits(0xbfe0000000000000), L_lo)
(L_hi, L_lo) = fast2sum(t1, L_hi)
L_lo += t2
/* r = r · x */
(L_hi, L_lo) = hl_mult(L_hi, L_lo, fF4_hi, fF4_lo)
/* r = r + 1 */
(t1, t2) = fast2sum(1.0, L_lo)
(L_hi, L_lo) = fast2sum(t1, L_hi)
L_lo += t2
/* r = r · x */
(L_hi, L_lo) = hl_mult(L_hi, L_lo, fF4_hi, fF4_lo)
/*
Finally, compute log(F1) + log(F2) + log(F3) + log(F4) + L. We may
assume F1 > F2 > F3 > F4 > F5, since the only way this is disrupted is
if some Fi == 1.0, in which case the log is 0 and the fast2sum works
out either way. We can also assume each F1,2,3 > L.
*/
var lsig_hi, lsig_lo
/* log(F4) + L, slow because they're the same order of magnitude */
(t1, t2) = slow2sum(std.flt64frombits(logF4_lo), L_lo)
(lsig_lo, t1) = slow2sum(L_hi, t1)
(lsig_hi, lsig_lo) = slow2sum(std.flt64frombits(logF4_hi), lsig_lo)
(lsig_hi, lsig_lo) = hl_add(std.flt64frombits(logF3_hi), std.flt64frombits(logF3_lo), lsig_hi, lsig_lo + (t1 + t2))
(lsig_hi, lsig_lo) = hl_add(std.flt64frombits(logF2_hi), std.flt64frombits(logF2_lo), lsig_hi, lsig_lo)
(lsig_hi, lsig_lo) = hl_add(std.flt64frombits(logF1_hi), std.flt64frombits(logF1_lo), lsig_hi, lsig_lo)
/* Oh yeah, and we need xe * log(2) */
var xel2_hi, xel2_lo, lx_hi, lx_lo
(xel2_hi, xel2_lo) = hl_mult((xe : flt64), 0.0, std.flt64frombits(0x3fe62e42fefa39ef), std.flt64frombits(0x3c7abc9e3b39803f))
(t1, t2) = slow2sum(xel2_lo, lsig_lo)
(lx_lo, t1) = slow2sum(lsig_hi, t1)
(lx_hi, lx_lo) = slow2sum(xel2_hi, lx_lo)
(lx_hi, lx_lo) = slow2sum(lx_hi, lx_lo + (t1 + t2))
if then_invert
-> (-1.0 * lx_hi, -1.0 * lx_lo)
;;
-> (lx_hi, lx_lo)
}