zig

fork of https://codeberg.org/ziglang/zig
Log | Files | Refs | README | LICENSE

commit 8bf2f80cf17cfab92315316882ba1aa3101ee379 (tree)
parent 3f1dead2fc5922b588fbfb108f421ca957d6934a
Author: Christophe Delage <christ@ophe.net>
Date:   Thu, 28 May 2026 17:01:16 +0200

Add f128 support for @log()

Diffstat:
Mlib/compiler_rt/log.zig | 294++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 291 insertions(+), 3 deletions(-)

diff --git a/lib/compiler_rt/log.zig b/lib/compiler_rt/log.zig @@ -8,6 +8,7 @@ const std = @import("std"); const math = std.math; const expect = std.testing.expect; const expectEqual = std.testing.expectEqual; +const expectApproxEqRel = std.testing.expectApproxEqRel; const compiler_rt = @import("../compiler_rt.zig"); const symbol = @import("../compiler_rt.zig").symbol; @@ -436,9 +437,234 @@ pub fn __logx(a: f80) callconv(.c) f80 { return @floatCast(logq(a)); } -pub fn logq(a: f128) callconv(.c) f128 { - // TODO: more correct implementation - return log(@floatCast(a)); +/// Implementation of "Table-driven implementation of the exponential function in IEEE floating-point arithmetic" +/// By PTP Tang in ACM Transactions on Mathematical Software (TOMS), 1989 +/// +/// https://dl.acm.org/doi/pdf/10.1145/63522.214389 +/// +/// Adapted to work for f128 by Christophe Delage. +/// +/// Accuracy on 100 million random numbers in [0, inf) (exponent uniformly random) +/// <= 0.5 ulp: 99.99%, worst case <= 0.530 ulp +/// +/// Accuracy on 10 million random numbers near x = 1 (testing the proc2 case): +/// <= 0.5 ulp: 99.85%, worst case <= 0.532 ulp +pub fn logq(x: f128) callconv(.c) f128 { + const bitsize = 7; + const size = 1 << bitsize; + + if (!math.isFinite(x)) { + if (math.isNan(x)) { + if (math.isSignalNan(x)) math.raiseInvalid(); + return math.nan(f128); + } + if (math.isPositiveInf(x)) return x; + } + if (x <= 0.0) { + if (x >= 0.0) { + math.raiseDivByZero(); + return -math.inf(f128); + } + math.raiseInvalid(); + return math.nan(f128); + } + + // exp(-1 / 16) rounded down + const proc2_lo: f128 = 0.939413062813475786119710824622305; + // exp(1 / 16) rounded up + const proc2_hi: f128 = 1.0644944589178594295633905946428897; + if (proc2_lo < x and x < proc2_hi) { + const f = x - 1.0; + + const g = 1 / (2 + f); + const u = 2 * f * g; + const v = u * u; + const uv = u * v; + + const v64: f64 = @floatCast(v); + + // use f64 for the last few coefficients to improve performance + const p19 = 2.0165671588771827537210411918018159e-7; + const p17 = 8.97568550755477160981619052649713e-7 + v64 * p19; + const p15 = 4.069010449774280288178309893970754e-6 + v64 * p17; + const p13 = 1.8780048076832339308077858301484127e-5 + v * p15; + const p11 = 8.87784090909092440759545734146088e-5 + v * p13; + const p9 = 4.340277777777777776216500817402857e-4 + v * p11; + const p7 = 2.2321428571428571428572328745789477e-3 + v * p9; + const p5 = 1.249999999999999999999999997455655e-2 + v * p7; + + const q_hi = uv * 8.333333333333333333333333333333581e-2; + const q_lo = uv * v * p5; + const q = q_hi + q_lo; + + const fa: f128 = @as(f64, @floatCast(f)); + const ua: f128 = @as(f64, @floatCast(u)); + + const fb: f128 = f - fa; + const ub: f128 = ((2 * (f - ua) - ua * fa) - ua * fb) * g; + + return ua + (ub + q); + } + + const ym = frexp2(x); + const y = ym.significand; + const m = ym.exponent; + + const F0 = @round(math.ldexp(y, bitsize)); + const j0: usize = @intFromFloat(F0); + const j = j0 - size; + const F = math.ldexp(F0, -bitsize); + const f = y - F; + + const u = (f + f) / (y + F); + const v = u * u; + + const v64: f64 = @floatCast(v); + // use f64 for the last coefficient to improve performance + const p11 = 8.877925718782769769445565656611838e-5; + const p9 = 4.340277777635300605611118803507141e-4 + v64 * p11; + const p7 = 2.2321428571428572515097318595359542e-3 + v * p9; + const p5 = 1.249999999999999999999963839743372e-2 + v * p7; + const p3 = 8.333333333333333333333333333372414e-2 + v * p5; + + const q = u * v * p3; + + const log1p_tab = [size + 1]struct { hi: f128, lo: f128 }{ + .{ .hi = 0, .lo = 0 }, + .{ .hi = 0x1.fe02a6b106788fc3769039p-8, .lo = 0x1.dc282d2b3db2c3ef9a073a876702p-100 }, + .{ .hi = 0x1.fc0a8b0fc03e3cf9eda74d4p-7, .lo = -0x1.0a8552414fc416fc223acca2ebfp-100 }, + .{ .hi = 0x1.7b91b07d5b11aa927f54c72p-6, .lo = -0x1.287fc46561dfab5bc5cceecb4882p-99 }, + .{ .hi = 0x1.f829b0e7833004cf8fc13c8p-6, .lo = -0x1.0dd605151051eb3220ca52e20939p-100 }, + .{ .hi = 0x1.39e87b9febd5fa9015b202bp-5, .lo = -0x1.1bac6e550a3c3dc859cfe2e178a8p-99 }, + .{ .hi = 0x1.77458f632dcfc4634f2a1eep-5, .lo = 0x1.2960b1e4dfb80d9544ec6583eb3ap-99 }, + .{ .hi = 0x1.b42dd711971bec28d14c7dap-5, .lo = -0x1.2645ad50c7672fc0eb08d862221dp-102 }, + .{ .hi = 0x1.f0a30c01162a6617cc9716fp-5, .lo = -0x1.4cd0ece597165991495b4d31cf5cp-101 }, + .{ .hi = 0x1.16536eea37ae0e8625c173ep-4, .lo = -0x1.66d0dc92deb7d2ccbd2caa9640cap-99 }, + .{ .hi = 0x1.341d7961bd1d092998376108p-4, .lo = -0x1.976457ef2f89ad243dcc3578cf7ep-99 }, + .{ .hi = 0x1.51b073f06183f69278e686ap-4, .lo = 0x1.7c8ac25e4e3f04de1f086f5cb4b1p-99 }, + .{ .hi = 0x1.6f0d28ae56b4b9be499b9edp-4, .lo = 0x1.9b640ce50c1ef65087fdf23812f4p-100 }, + .{ .hi = 0x1.8c345d6319b20f5acb42a66p-4, .lo = -0x1.254bca8fd9fc1bf283b3b4b8662dp-100 }, + .{ .hi = 0x1.a926d3a4ad563650bd22a9cp-4, .lo = 0x1.d5263cd4fb3f11769cc680ef5589p-99 }, + .{ .hi = 0x1.c5e548f5bc74315d617ef818p-4, .lo = -0x1.e4e896269950723c88d353ee9c18p-100 }, + .{ .hi = 0x1.e27076e2af2e5e9ea87ffe2p-4, .lo = -0x1.61eaa246b143bfe80906a822f768p-104 }, + .{ .hi = 0x1.fec9131dbeabaaa2e5199f9p-4, .lo = 0x1.9271dff48f15d409017630a93931p-99 }, + .{ .hi = 0x1.0d77e7cd08e596697717a40cp-3, .lo = 0x1.574712132d3f6340e183be2031c6p-102 }, + .{ .hi = 0x1.1b72ad52f67a029060468e58p-3, .lo = 0x1.ae73f3bc7ec84ca997609536a037p-99 }, + .{ .hi = 0x1.29552f81ff5234c05dc7102p-3, .lo = -0x1.20b2ef60436f8f081d60452c9fc1p-100 }, + .{ .hi = 0x1.371fc201e8f743bcd96c55e4p-3, .lo = -0x1.d80d17e0cd92558ad6fcd608bb1bp-100 }, + .{ .hi = 0x1.44d2b6ccb7d1e67d3d950f88p-3, .lo = -0x1.e1f3be9a83374584faad83fa4fecp-103 }, + .{ .hi = 0x1.526e5e3a1b437a2e401d6e3cp-3, .lo = 0x1.6334db798c76a888aa87317b14f8p-100 }, + .{ .hi = 0x1.5ff3070a793d3c873e20a074p-3, .lo = -0x1.edc45019551c65501060ce71fa98p-99 }, + .{ .hi = 0x1.6d60fe719d21c8d54765c4ccp-3, .lo = -0x1.790e412e6d3ed18e4a7c22362e49p-101 }, + .{ .hi = 0x1.7ab890210d9091be36b2d6ap-3, .lo = 0x1.820191ff8525362042cad5d8c597p-101 }, + .{ .hi = 0x1.87fa06520c910902009017dcp-3, .lo = 0x1.32ef5a55704b6b7eb4ebea28a6cdp-100 }, + .{ .hi = 0x1.9525a9cf456b47641307538cp-3, .lo = -0x1.da62766be8258611d132d71d84acp-101 }, + .{ .hi = 0x1.a23bc1fe2b563193711b07a8p-3, .lo = 0x1.98c27e3f1b66d8b32de61c04cf95p-99 }, + .{ .hi = 0x1.af3c94e80bff2d8ce601937cp-3, .lo = 0x1.9eb976769b8b9a5d50ca14a7622fp-100 }, + .{ .hi = 0x1.bc286742d8cd629f9ce890ep-3, .lo = 0x1.ea9e1e2c3dca46c83cd6d19e5e9bp-99 }, + .{ .hi = 0x1.c8ff7c79a9a21ac25d81ef3p-3, .lo = -0x1.1976d471342b17dca47c9d1e1d98p-105 }, + .{ .hi = 0x1.d5c216b4fbb915b910d65f94p-3, .lo = -0x1.5ff1e1c98c2ed4063968ad2332f8p-100 }, + .{ .hi = 0x1.e27076e2af2e5e9ea87ffe2p-3, .lo = -0x1.61eaa246b143bfe80906a822f768p-103 }, + .{ .hi = 0x1.ef0adcbdc59365218de5437p-3, .lo = 0x1.06429f5a5098683a47c3a3ca5835p-100 }, + .{ .hi = 0x1.fb9186d5e3e2a8d55466c378p-3, .lo = 0x1.4d2ca09202c222ad91640bd8b2fcp-99 }, + .{ .hi = 0x1.0402594b4d040dae27bd0b6p-2, .lo = -0x1.16a1bbb899f343f105ee37cafa25p-100 }, + .{ .hi = 0x1.0a324e27390e35f73f7a0188p-2, .lo = -0x1.fe78eb90fe52820fb4690e4910e3p-99 }, + .{ .hi = 0x1.1058bf9ae4ad5189fa0ab4ccp-2, .lo = -0x1.9c60f598d3a32076376a5960e735p-99 }, + .{ .hi = 0x1.1675cababa60e039cc7d571p-2, .lo = 0x1.b8b823f067d04a43c19f534c3c8ep-100 }, + .{ .hi = 0x1.1c898c16999fafbc68e75404p-2, .lo = -0x1.c443cc477d114a1350a26c9b5535p-100 }, + .{ .hi = 0x1.22941fbcf7965a242853da76p-2, .lo = -0x1.5e685a2caa590b0f13a31703b136p-101 }, + .{ .hi = 0x1.2895a13de86a35eb49304fc2p-2, .lo = -0x1.f8d3a52b8aa6834f50a903b31253p-99 }, + .{ .hi = 0x1.2e8e2bae11d309c2cc91a85p-2, .lo = 0x1.03679bdbbd6b7d91378a909a6793p-99 }, + .{ .hi = 0x1.347dd9a987d54d645674fedcp-2, .lo = 0x1.821ee510a580b3b30ef57f83ca28p-99 }, + .{ .hi = 0x1.3a64c556945e9c72f35cd74p-2, .lo = 0x1.a11beb7a3cee7e029e46e1334dfep-99 }, + .{ .hi = 0x1.404308686a7e3bd0c127df4cp-2, .lo = 0x1.92985641827d9d91a2da0d4f2207p-100 }, + .{ .hi = 0x1.4618bc21c5ec27d0b7b37b34p-2, .lo = -0x1.c65df511a65b67fe9778d8694229p-101 }, + .{ .hi = 0x1.4be5f957778a0db4c9949f7p-2, .lo = -0x1.3cdc28d5974f3185a1f581529353p-101 }, + .{ .hi = 0x1.51aad872df82d09c93d60cfap-2, .lo = 0x1.5e311d4f4f357cbfae095f10fcd3p-99 }, + .{ .hi = 0x1.5767717455a6c549ab6ca0dap-2, .lo = -0x1.f42ff0747cbcce6c0841fd7ceb87p-100 }, + .{ .hi = 0x1.5d1bdbf5809ca508d8e0f72p-2, .lo = -0x1.eea60c7f4b594bd65b44f6b20634p-104 }, + .{ .hi = 0x1.62c82f2b9c7952f6f5f22a6p-2, .lo = 0x1.ca2e7226c55dd257f44b5002c8cdp-102 }, + .{ .hi = 0x1.686c81e9b14aec442be1014ep-2, .lo = 0x1.c34b25bdbda38672b32ca47f535dp-101 }, + .{ .hi = 0x1.6e08eaa2ba1e38c139318d72p-2, .lo = -0x1.07a1f9d2a3058cd3f047b933d485p-99 }, + .{ .hi = 0x1.739d7f6bbd0069ce24c53faep-2, .lo = -0x1.8210d2a910f7918ad34542221b6cp-99 }, + .{ .hi = 0x1.792a55fdd47a27c15da47fa8p-2, .lo = -0x1.297ea603cd10e7c702842a1590aep-100 }, + .{ .hi = 0x1.7eaf83b82afc364b3a5e7b4ap-2, .lo = 0x1.50437160cbfcbf71ee8d4b3cd067p-100 }, + .{ .hi = 0x1.842d1da1e8b17493b1465e12p-2, .lo = -0x1.899770fb9eb8e0c7f06a12cd88c3p-100 }, + .{ .hi = 0x1.89a3386c1425ab5a71881104p-2, .lo = -0x1.f549800739afc97f5d4b9c15fc12p-101 }, + .{ .hi = 0x1.8f11e873662c77e1769d5698p-2, .lo = 0x1.a29979cbfcbc0e45410ee8ca0d17p-100 }, + .{ .hi = 0x1.947941c2116faba4cdd147d2p-2, .lo = -0x1.f22a2b6b19ed11af82f2c0e6730ep-99 }, + .{ .hi = 0x1.99d958117e08acba92eec478p-2, .lo = 0x1.8dfba4bb71c95f44a3225b5d8f8fp-101 }, + .{ .hi = 0x1.9f323ecbf984bf2b68d766f4p-2, .lo = 0x1.4886067d20ffb34547d7c2b38ad8p-104 }, + .{ .hi = 0x1.a484090e5bb0a2bfca6b70ecp-2, .lo = -0x1.62da6c6290af3949ca12eb97e6bbp-99 }, + .{ .hi = 0x1.a9cec9a9a08498d484ff52f2p-2, .lo = 0x1.50d7be11fc8ee26768c44e9f35aap-100 }, + .{ .hi = 0x1.af1293247786b1133844a15ep-2, .lo = -0x1.ebf9e3b2fb68378f9b7aa0ef6685p-101 }, + .{ .hi = 0x1.b44f77bcc8f628cbeedaae98p-2, .lo = 0x1.c3c779029d5a684301e7888d5449p-99 }, + .{ .hi = 0x1.b9858969310fb598fb14f88ep-2, .lo = 0x1.e1cf79039d5a31010282d8264afap-99 }, + .{ .hi = 0x1.beb4d9da71b7bf7861d37abcp-2, .lo = -0x1.f29b495a7c83dffb9899ad6da116p-101 }, + .{ .hi = 0x1.c3dd7a7cdad4d73b3c14b7aap-2, .lo = -0x1.649f44ae71a827ebf6601862d429p-100 }, + .{ .hi = 0x1.c8ff7c79a9a21ac25d81ef3p-2, .lo = -0x1.1976d471342b17dca47c9d1e1d98p-104 }, + .{ .hi = 0x1.ce1af0b85f3eb7b7d2bcaadp-2, .lo = 0x1.33a4e218c8c00b2f6f2b23998791p-99 }, + .{ .hi = 0x1.d32fe7e00ebd561dec8cbebep-2, .lo = 0x1.1b7e8c55a6ffaf7eb72c1ff893cdp-99 }, + .{ .hi = 0x1.d83e7258a2f3e50515ba2ecap-2, .lo = -0x1.7778ad7b599f3edc40a13fe57cp-99 }, + .{ .hi = 0x1.dd46a04c1c4a0bee626a49d2p-2, .lo = -0x1.23c02c15f2f66328a06054db5e01p-101 }, + .{ .hi = 0x1.e24881a7c6c261cbd8f45954p-2, .lo = 0x1.48c6c5403df1764e1ed219643a85p-99 }, + .{ .hi = 0x1.e744261d68787e37da36f3ccp-2, .lo = -0x1.2e44bfb5e71d55a11feaedf56a94p-100 }, + .{ .hi = 0x1.ec399d2468cc0175cee53f36p-2, .lo = -0x1.8d2027bb4681af8a138d77633327p-99 }, + .{ .hi = 0x1.f128f5faf06ecb35c83b1132p-2, .lo = -0x1.851461536ddd17ac271709c2b464p-101 }, + .{ .hi = 0x1.f6123fa7028ac61456c3cb6cp-2, .lo = 0x1.a0a432a2414cc0b049c0fb73b4dep-99 }, + .{ .hi = 0x1.faf588f78f31ed9afb3e4ea8p-2, .lo = 0x1.afec6d4cde2ef184dc7b6e634ba1p-100 }, + .{ .hi = 0x1.ffd2e0857f4985597d0364c8p-2, .lo = 0x1.af246e3379d3ea68e29866b3b38bp-100 }, + .{ .hi = 0x1.02552a5a5d0fec69c695d7eep-1, .lo = 0x1.fff1a3725c5c6261a75dc4f512adp-99 }, + .{ .hi = 0x1.04bdf9da926d265fcc1008b2p-1, .lo = 0x1.df6a6d08e4470f10c7053f04f1ep-99 }, + .{ .hi = 0x1.0723e5c1cdf404e579638911p-1, .lo = 0x1.baad95bd2eba9089087bae740509p-99 }, + .{ .hi = 0x1.0986f4f573520b91fda94ff4p-1, .lo = 0x1.fe038113f135a26389f13c4bde3ap-100 }, + .{ .hi = 0x1.0be72e4252a82b69897bb33ep-1, .lo = -0x1.9649bc990440ca2c12ee56f6c90bp-108 }, + .{ .hi = 0x1.0e44985d1cc8bf6eae5de969p-1, .lo = 0x1.8f8d1fbacf70a2da185e84859e36p-99 }, + .{ .hi = 0x1.109f39e2d4c96fde3ec9b0b8p-1, .lo = -0x1.b9d06df4e0c8337c0da8f63aefc7p-99 }, + .{ .hi = 0x1.12f719593efbc53012319c7dp-1, .lo = 0x1.a96893b2757f50123379aecd147cp-102 }, + .{ .hi = 0x1.154c3d2f4d5e9a98f33a3966p-1, .lo = -0x1.d7f56258b99e197c61c0a23b2403p-101 }, + .{ .hi = 0x1.179eabbd899a0bfc60e6fa08p-1, .lo = -0x1.68f2a71c8279b89eb8392b7a41ep-100 }, + .{ .hi = 0x1.19ee6b467c96ecc5cbdd7782p-1, .lo = -0x1.0c2a8ef8715f93d3c8e2c9016713p-100 }, + .{ .hi = 0x1.1c3b81f713c24bc94f8ecdfcp-1, .lo = -0x1.d103880b3dc864d107231c87eaa7p-100 }, + .{ .hi = 0x1.1e85f5e7040d03dec59a5f3ep-1, .lo = 0x1.e35f2e7ca4f4817696ad39f9a4b2p-100 }, + .{ .hi = 0x1.20cdcd192ab6d93503d0f75cp-1, .lo = -0x1.3db0bb5bf2b76be256c1a2a08a8p-103 }, + .{ .hi = 0x1.23130d7bebf4282de368722dp-1, .lo = -0x1.6dfadc7219019f05279c9cfd2b08p-99 }, + .{ .hi = 0x1.2555bce98f7cb3c043adad1p-1, .lo = 0x1.e71084750a06eb301ae8308feca1p-100 }, + .{ .hi = 0x1.2795e1289b11aeb783f3db97p-1, .lo = -0x1.e3801fe56c1467b5e622105c5e41p-99 }, + .{ .hi = 0x1.29d37fec2b08ac85cd6cba5dp-1, .lo = -0x1.824a2f303e0aa6995ef5598cc5b7p-100 }, + .{ .hi = 0x1.2c0e9ed448e8bb97a9c31ba3p-1, .lo = -0x1.8676adfad5c83dea45d7349693e3p-99 }, + .{ .hi = 0x1.2e47436e4026840542186922p-1, .lo = 0x1.ab1252cf29452c88ebc5ce2f36f2p-101 }, + .{ .hi = 0x1.307d7334f10be1fb590a1f56p-1, .lo = 0x1.b66f70c05b80999c08cdd48a021p-99 }, + .{ .hi = 0x1.32b1339121d71320556b67b2p-1, .lo = 0x1.6b06715ba133462ddf94bde6b3a5p-100 }, + .{ .hi = 0x1.34e289d9ce1d316eb92d885dp-1, .lo = -0x1.b151b59c44058fa92837dec71351p-101 }, + .{ .hi = 0x1.37117b54747b5c5dd024844ep-1, .lo = -0x1.037076a726793d7e93c9b2f3eef6p-100 }, + .{ .hi = 0x1.393e0d3562a19a9c4426036ep-1, .lo = -0x1.cf1c277a3a0d863fefb367ef8616p-102 }, + .{ .hi = 0x1.3b68449fffc22af8edec1859p-1, .lo = 0x1.b341d6de6d9b9591a54790d29adcp-100 }, + .{ .hi = 0x1.3d9026a7156faa404263d0adp-1, .lo = 0x1.3cf6b809d96954adf1ff9360bd04p-100 }, + .{ .hi = 0x1.3fb5b84d16f425b4e9d505cbp-1, .lo = -0x1.110c9bd8986629a5a46a5834774cp-99 }, + .{ .hi = 0x1.41d8fe84672ae6464bcc2f46p-1, .lo = 0x1.779538890dd44eadeb32e848f817p-105 }, + .{ .hi = 0x1.43f9fe2f9ce677a727b9b60fp-1, .lo = -0x1.e6927c891167a0306333dca3b8a4p-101 }, + .{ .hi = 0x1.4618bc21c5ec27d0b7b37b34p-1, .lo = -0x1.c65df511a65b67fe9778d8694229p-100 }, + .{ .hi = 0x1.48353d1ea88df73d5e8bb302p-1, .lo = -0x1.7b4f3e104187cc8aca358d9263f9p-104 }, + .{ .hi = 0x1.4a4f85db03ebb0227bf47a6fp-1, .lo = -0x1.e49ce91908e6e2750b818bba40e6p-100 }, + .{ .hi = 0x1.4c679afccee39b168ecdd318p-1, .lo = 0x1.bf619db0d889bbe38f559618dbd3p-99 }, + .{ .hi = 0x1.4e7d811b75bb09cb09856458p-1, .lo = 0x1.5770d0c5ebca2047bba2c9ee4f52p-99 }, + .{ .hi = 0x1.50913cc01686b4bcb3a5b0b5p-1, .lo = 0x1.b0f69791cf6c54c4e5d96f1d584fp-99 }, + .{ .hi = 0x1.52a2d265bc5aaee77c8af15bp-1, .lo = 0x1.7aea7bed0920f6a4c39b31ea388bp-100 }, + .{ .hi = 0x1.54b2467999497a915428b43ep-1, .lo = -0x1.f434bb5d154a84758a2a5033748cp-99 }, + .{ .hi = 0x1.56bf9d5b3f399411c6217364p-1, .lo = -0x1.a6323ea9ce40a3caf6baebad2c64p-104 }, + .{ .hi = 0x1.58cadb5cd79893092f25d931p-1, .lo = -0x1.d3d5761cff2a6a51ddc9c241881cp-99 }, + .{ .hi = 0x1.5ad404c359f2cfb29aaa5f02p-1, .lo = 0x1.cd40845839e04578161ccf77753bp-100 }, + .{ .hi = 0x1.5cdb1dc6c17648cf6e3c5d71p-1, .lo = -0x1.f3a84fed7a8e6b8a7923783c6bf7p-99 }, + .{ .hi = 0x1.5ee02a924167570d6095fd24p-1, .lo = 0x1.03ef4c629f04ae4e7a29309e6688p-100 }, + .{ .hi = 0x1.60e32f44788d8ca7c895a0b5p-1, .lo = -0x1.3557995d063914a66aa81ead3fdbp-101 }, + .{ .hi = 0x1.62e42fefa39ef35793c7673p-1, .lo = 0x1.f97b57a079a193394c5b16c5068cp-103 }, + }; + const xm: f128 = @floatFromInt(m); + const l_hi = xm * log1p_tab[128].hi + log1p_tab[j].hi; + const l_lo = xm * log1p_tab[128].lo + log1p_tab[j].lo; + + return l_hi + (u + (q + l_lo)); } pub fn logl(x: c_longdouble) callconv(.c) c_longdouble { @@ -450,6 +676,30 @@ pub fn logl(x: c_longdouble) callconv(.c) c_longdouble { } } +/// Returns (f, k) such that x = f * 2^k and f in [1,2). +/// Asserts that x is finite and positive. +pub fn frexp2(x: f128) math.Frexp(f128) { + std.debug.assert(math.isFinite(x)); + std.debug.assert(x > 0.0); + + const bits: u128 = @bitCast(x); + const uexp: i32 = @intCast(bits >> 112); + + std.debug.assert(uexp >= 0); + + if (uexp == 0) { + const shift: u7 = @intCast(@clz(bits) - 15); + + const exp = -@as(i32, shift) - 0x3ffe; + const frac: f128 = @bitCast((bits << shift) | (0x3fff << 112)); + return .{ .significand = frac, .exponent = exp }; + } + + const exp = uexp - 0x3fff; + const frac: f128 = @bitCast((0x3fff << 112) | ((bits << 16) >> 16)); + return .{ .significand = frac, .exponent = exp }; +} + test "logf() special" { try expectEqual(logf(0.0), -math.inf(f32)); try expectEqual(logf(-0.0), -math.inf(f32)); @@ -519,3 +769,41 @@ test "log() boundary" { try expectEqual(log(0x1p-1022), -0x1.6232bdd7abcd2p+9); // First subnormal try expect(math.isNan(log(-0x1p-1022))); // First negative subnormal } + +test "logq() special" { + try expectEqual(logq(0.0), -math.inf(f128)); + try expectEqual(logq(-0.0), -math.inf(f128)); + try expect(math.isPositiveZero(logq(1.0))); + // Sadly, the rounding gods decided that 0.9999999999999999999999999999999999 + // is the correctly rounded value of logq(math.e) + try expectApproxEqRel(logq(math.e), 1.0, math.floatEpsAt(f128, 1.0)); + try expectEqual(logq(math.inf(f128)), math.inf(f128)); + try expect(math.isNan(logq(-1.0))); + try expect(math.isNan(logq(-math.inf(f128)))); + try expect(math.isNan(logq(math.nan(f128)))); + try expect(math.isNan(logq(math.snan(f128)))); +} + +test "logq() boundary" { + try expectEqual(logq(0x1.ffffffffffffffffffffffffffffp16383), 0x1.62e42fefa39ef35793c7673007e6p13); // Max input value + try expectEqual(logq(0x1p-16494), -0x1.6546282207802c89d24d65e96274p13); // Min positive input value + try expect(math.isNan(logq(-0x1p-16494))); // Min negative input value + try expectEqual(logq(0x1.0000000000000000000000000001p0), 0x1.ffffffffffffffffffffffffffffp-113); // Last value before result reaches +0 + try expectEqual(logq(0x1.ffffffffffffffffffffffffffffp-1), -0x1p-113); // Last value before result reaches -0 + try expectEqual(logq(0x1p-16382), -0x1.62d918ce2421d65ff90ac8f4ce66p13); // First subnormal + try expect(math.isNan(logq(-0x1p-16382))); // First negative subnormal +} + +test "logq() sanity" { + try expectEqual(logq(4.151135979023751199079583784623537e-4), -7.7869583453055243113993340258295346e0); + try expectEqual(logq(9.614234245933828353176667689130293e-14), -2.9972946567656004014786271559909435e1); + try expectEqual(logq(1.012889803704721484375e13), 2.9946413646144315985379677542014356e1); + try expectEqual(logq(2.397741857206453154086912e24), 5.613656963346284538829358703465392e1); + try expectEqual(logq(3.442377567808290806386655232e27), 6.3405959896920645453203836625419693e1); + try expectEqual(logq(1.0689155158234028407981544637594257e-8), -1.835403614606774451014272772421113e1); + try expectEqual(logq(1.4813913545768791536741499811327596e-10), -2.263286917934202003739900705050399e1); + try expectEqual(logq(4.518948965781299591064453125e10), 2.453413036705097282892685629562292e1); + try expectEqual(logq(1.200355637363589375e14), 3.2418809179272977400408325788186897e1); + try expectEqual(logq(6.6145398293682003021240234375e9), 2.261253606737223221601998075023261e1); + try expectEqual(logq(5.16179116383965741056e20), 4.7692985503915646405875629300054525e1); +}