1/// Implementation of "Table-driven implementation of the logarithm function in IEEE floating-point arithmetic"
2/// by PTP Tang in ACM Transactions on Mathematical Software (TOMS), 1990
3///
4/// https://dl.acm.org/doi/pdf/10.1145/98267.98294
5///
6/// Adapted to work for f128 and bases 2 and 10 by Christophe Delage.
7///
8/// This file contains the code shared between logq, log2q and log10q.
9const log_f128 = @This();
10
11const std = @import("std");
12const math = std.math;
13
14pub const log2size = 7;
15pub const size = 1 << log2size;
16
17/// Filter out special cases for log in bases {e,2,10}.
18///
19/// If x is finite and positive, returns null.
20/// Returns the appropriate NaN or inf otherwise.
21pub fn specialCases(x: f128) ?f128 {
22 if (!math.isFinite(x)) {
23 if (math.isNan(x)) {
24 if (math.isSignalNan(x)) math.raiseInvalid();
25 return math.nan(f128);
26 }
27 if (math.isPositiveInf(x)) return x;
28 }
29 if (x <= 0.0) {
30 if (x >= 0.0) {
31 math.raiseDivByZero();
32 return -math.inf(f128);
33 }
34 math.raiseInvalid();
35 return math.nan(f128);
36 }
37
38 return null;
39}
40
41pub const Proc1 = struct {
42 pub const Poly = struct {
43 a1: f128,
44 a3: f128,
45 a5: f128,
46 a7: f128,
47 a9: f64,
48 a11: f64,
49 };
50 pub const HiLo = struct { hi: f128, lo: f128 };
51 poly: Poly,
52 tab: [size + 1]HiLo,
53};
54
55pub fn proc1(comptime p: Proc1, x: f128) f128 {
56 const ym = frexp2(x);
57 const y = ym.significand;
58 const m = ym.exponent;
59
60 const F0 = @round(math.ldexp(y, log2size));
61 const j0: usize = @intFromFloat(F0);
62 const j = j0 - size;
63 const F = math.ldexp(F0, -log2size);
64 const f = y - F;
65
66 const u = (f + f) / (y + F);
67 const v = u * u;
68 const v64: f64 = @floatCast(v);
69
70 const p9 = p.poly.a9 + v64 * p.poly.a11;
71 const p7 = p.poly.a7 + v * p9;
72 const p5 = p.poly.a5 + v * p7;
73 const p3 = p.poly.a3 + v * p5;
74
75 const q = u * v * p3;
76
77 const xm: f128 = @floatFromInt(m);
78 const l_hi = xm * p.tab[128].hi + p.tab[j].hi;
79 const l_lo = xm * p.tab[128].lo + p.tab[j].lo;
80
81 if (comptime p.poly.a1 == 1.0)
82 return l_hi + (u + (q + l_lo))
83 else
84 return l_hi + (u * p.poly.a1 + (q + l_lo));
85}
86
87pub const Proc2 = struct {
88 // exp(-1 / 16) rounded down
89 pub const lo: f128 = 0.939413062813475786119710824622305;
90 // exp(1 / 16) rounded up
91 pub const hi: f128 = 1.0644944589178594295633905946428897;
92
93 pub const Poly = struct {
94 b1_hi: f128,
95 b1_lo: f128,
96 b3: f128,
97 b5: f128,
98 b7: f128,
99 b9: f128,
100 b11: f128,
101 b13: f128,
102 b15: f64,
103 b17: f64,
104 b19: f64,
105 };
106
107 poly: Poly,
108};
109
110pub fn proc2(comptime p: Proc2, x: f128) f128 {
111 std.debug.assert(Proc2.lo < x and x < Proc2.hi);
112
113 const f = x - 1.0;
114 const g = 1 / (2 + f);
115 const u = 2 * f * g;
116 const v = u * u;
117 const uv = u * v;
118 const v64: f64 = @floatCast(v);
119
120 const p17 = p.poly.b17 + v64 * p.poly.b19;
121 const p15 = p.poly.b15 + v64 * p17;
122 const p13 = p.poly.b13 + v * p15;
123 const p11 = p.poly.b11 + v * p13;
124 const p9 = p.poly.b9 + v * p11;
125 const p7 = p.poly.b7 + v * p9;
126 const p5 = p.poly.b5 + v * p7;
127
128 const q_hi = uv * p.poly.b3;
129 const q_lo = uv * v * p5;
130
131 const f_hi: f128 = @as(f64, @floatCast(f));
132 const f_lo = f - f_hi;
133
134 const u_hi: f128 = @as(f64, @floatCast(u));
135 const u_lo = ((2 * (f - u_hi) - u_hi * f_hi) - u_hi * f_lo) * g;
136
137 if (comptime p.poly.b1_hi == 1.0 and p.poly.b1_lo == 0.0)
138 return u_hi + (u_lo + (q_hi + q_lo));
139
140 // t = u * p.poly.b1
141 const t_hi = u_hi * p.poly.b1_hi;
142 const t_lo = u_lo * p.poly.b1_hi + u * p.poly.b1_lo;
143
144 // y = t + q
145 const y_hi = t_hi + q_hi;
146 const y_lo = t_lo + (t_hi - y_hi + q_hi) + q_lo;
147
148 return y_hi + y_lo;
149}
150
151/// Returns (f, k) such that x = f * 2^k and f in [1,2).
152/// Asserts that x is finite and positive.
153pub fn frexp2(x: f128) math.Frexp(f128) {
154 std.debug.assert(math.isFinite(x));
155 std.debug.assert(x > 0.0);
156
157 const bits: u128 = @bitCast(x);
158 const uexp: i32 = @intCast(bits >> 112);
159
160 std.debug.assert(uexp >= 0);
161
162 if (uexp == 0) {
163 const shift: u7 = @intCast(@clz(bits) - 15);
164
165 const exp = -@as(i32, shift) - 0x3ffe;
166 const frac: f128 = @bitCast((bits << shift) | (0x3fff << 112));
167 return .{ .significand = frac, .exponent = exp };
168 }
169
170 const exp = uexp - 0x3fff;
171 const frac: f128 = @bitCast((0x3fff << 112) | ((bits << 16) >> 16));
172 return .{ .significand = frac, .exponent = exp };
173}