authorgravatar for christ@ophe.netChristophe Delage <christ@ophe.net> 2026-05-29 13:53:34+02:00
committergravatar for christ@ophe.netChristophe Delage <christ@ophe.net> 2026-05-29 13:53:34+02:00
log4cd396fc32cdc4488ecd9a36f2067b61142d727e
tree443edd74e71c199382fc0e4ac950645736850692
parent5e8a9ed9e87f9152e40c7493b57b0622ceb20f29

Finish the code sharing between logq, log2q and log10q


4 files changed, 95 insertions(+), 192 deletions(-)

lib/compiler_rt/log.zig+2-16
...@@ -452,21 +452,8 @@ pub fn __logx(a: f80) callconv(.c) f80 {...@@ -452,21 +452,8 @@ pub fn __logx(a: f80) callconv(.c) f80 {
452pub fn logq(x: f128) callconv(.c) f128 {452pub fn logq(x: f128) callconv(.c) f128 {
453 const impl = @import("log_f128.zig");453 const impl = @import("log_f128.zig");
454454
455 if (!math.isFinite(x)) {455 if (impl.specialCases(x)) |y|
456 if (math.isNan(x)) {456 return y;
457 if (math.isSignalNan(x)) math.raiseInvalid();
458 return math.nan(f128);
459 }
460 if (math.isPositiveInf(x)) return x;
461 }
462 if (x <= 0.0) {
463 if (x >= 0.0) {
464 math.raiseDivByZero();
465 return -math.inf(f128);
466 }
467 math.raiseInvalid();
468 return math.nan(f128);
469 }
470457
471 if (impl.Proc2.lo < x and x < impl.Proc2.hi) {458 if (impl.Proc2.lo < x and x < impl.Proc2.hi) {
472 // Polynomial approximation of log((1 + u / 2) / (1 - u / 2))459 // Polynomial approximation of log((1 + u / 2) / (1 - u / 2))
...@@ -634,7 +621,6 @@ pub fn logq(x: f128) callconv(.c) f128 {...@@ -634,7 +621,6 @@ pub fn logq(x: f128) callconv(.c) f128 {
634 .{ .hi = 0x1.60e32f44788d8ca7c895a0b5p-1, .lo = -0x1.3557995d063914a66aa81ead3fdbp-101 },621 .{ .hi = 0x1.60e32f44788d8ca7c895a0b5p-1, .lo = -0x1.3557995d063914a66aa81ead3fdbp-101 },
635 .{ .hi = 0x1.62e42fefa39ef35793c7673p-1, .lo = 0x1.f97b57a079a193394c5b16c5068cp-103 },622 .{ .hi = 0x1.62e42fefa39ef35793c7673p-1, .lo = 0x1.f97b57a079a193394c5b16c5068cp-103 },
636 };623 };
637
638 return impl.proc1(.{ .poly = poly, .tab = tab }, x);624 return impl.proc1(.{ .poly = poly, .tab = tab }, x);
639}625}
640626
lib/compiler_rt/log10.zig+29-87
...@@ -183,101 +183,47 @@ pub fn __log10x(a: f80) callconv(.c) f80 {...@@ -183,101 +183,47 @@ pub fn __log10x(a: f80) callconv(.c) f80 {
183/// Accuracy on 10 million random numbers near x = 1 (testing the proc2 case):183/// Accuracy on 10 million random numbers near x = 1 (testing the proc2 case):
184/// <= 0.5 ulp: 99.96%, worst case <= 0.565 ulp184/// <= 0.5 ulp: 99.96%, worst case <= 0.565 ulp
185pub fn log10q(x: f128) callconv(.c) f128 {185pub fn log10q(x: f128) callconv(.c) f128 {
186 if (!math.isFinite(x)) {186 const impl = @import("log_f128.zig");
187 if (math.isNan(x)) {187
188 if (math.isSignalNan(x)) math.raiseInvalid();188 if (impl.specialCases(x)) |y|
189 return math.nan(f128);189 return y;
190 }
191 if (math.isPositiveInf(x)) return x;
192 }
193 if (x <= 0.0) {
194 if (x >= 0.0) {
195 math.raiseDivByZero();
196 return -math.inf(f128);
197 }
198 math.raiseInvalid();
199 return math.nan(f128);
200 }
201 const bitsize = 7;
202 const size = 1 << bitsize;
203
204 // exp10(-1 / 16) rounded down
205 const proc2_lo: f128 = 0.939413062813475786119710824622305;
206 // exp10(1 / 16) rounded up
207 const proc2_hi: f128 = 1.0644944589178594295633905946428897;
208 if (proc2_lo < x and x < proc2_hi) {
209 const f = x - 1.0;
210 const g = 1 / (2 + f);
211 const u = 2 * f * g;
212 const v = u * u;
213 const uv = u * v;
214 const v64: f64 = @floatCast(v);
215190
191 if (impl.Proc2.lo < x and x < impl.Proc2.hi) {
216 // Polynomial approximation of log10((1 + u / 2) / (1 - u / 2))192 // Polynomial approximation of log10((1 + u / 2) / (1 - u / 2))
217 // in [2 * a / (2 + a), 2 * b / (2 + b)]193 // in [2 * a / (2 + a), 2 * b / (2 + b)]
218 // where a = exp(-1 / 16) - 1 and b = exp(1 / 16) - 1194 // where a = exp(-1 / 16) - 1 and b = exp(1 / 16) - 1
219 const p19 = 8.757839894876785986064901881670424e-8;195 const poly: impl.Proc2.Poly = .{
220 const p17 = 3.898090687230025454479255130305971e-7 + v64 * p19;196 .b1_hi = 0x1.bcb7b1526e50ep-2,
221 const p15 = 1.7671487851436387503930903139882346e-6 + v64 * p17;197 .b1_lo = 0x1.95355baaafad33dc323ee3460246p-57,
222 const p13 = 8.156071249646061672592451832809541e-6 + v * p15;198 .b3 = 3.619120682527098563759407657638483e-2,
223 const p11 = 3.855597318033137224139238937796866e-5 + v * p13;199 .b5 = 5.428681023790647845639111475407614e-3,
224 const p9 = 1.8849586888161971679466375197398104e-4 + v * p11;200 .b7 = 9.694073256769014010070232880860484e-4,
225 const p7 = 9.694073256769014010070232880860484e-4 + v * p9;201 .b9 = 1.8849586888161971679466375197398104e-4,
226 const p5 = 5.428681023790647845639111475407614e-3 + v * p7;202 .b11 = 3.855597318033137224139238937796866e-5,
227 const p3 = 3.619120682527098563759407657638483e-2;203 .b13 = 8.156071249646061672592451832809541e-6,
228204 .b15 = 1.7671487851436387503930903139882346e-6,
229 const q_hi = uv * p3;205 .b17 = 3.898090687230025454479255130305971e-7,
230 const q_lo = uv * v * p5;206 .b19 = 8.757839894876785986064901881670424e-8,
231207 };
232 const f_hi: f128 = @as(f64, @floatCast(f));208 return impl.proc2(.{ .poly = poly }, x);
233 const f_lo: f128 = f - f_hi;
234
235 const u_hi: f128 = @as(f64, @floatCast(u));
236 const u_lo: f128 = ((2 * (f - u_hi) - u_hi * f_hi) - u_hi * f_lo) * g;
237
238 const log10e_hi: f128 = 0x1.bcb7b1526e50ep-2;
239 const log10e_lo: f128 = 0x1.95355baaafad33dc323ee3460246p-57;
240 // t = u / log(10)
241 const t_hi = u_hi * log10e_hi;
242 const t_lo = u_lo * log10e_hi + u * log10e_lo;
243
244 // y = t + q
245 const y_hi = t_hi + q_hi;
246 const y_lo = t_lo + (t_hi - y_hi + q_hi) + q_lo;
247
248 return y_hi + y_lo;
249 }209 }
250210
251 const ym = @import("log.zig").frexp2(x);
252 const y = ym.significand;
253 const m = ym.exponent;
254
255 const F0 = @round(math.ldexp(y, bitsize));
256 const j0: usize = @intFromFloat(F0);
257 const j = j0 - size;
258 const F = math.ldexp(F0, -bitsize);
259 const f = y - F;
260
261 const u = (f + f) / (y + F);
262 const v = u * u;
263 const v64: f64 = @floatCast(v);
264
265 // Polynomial approximation of log10(1 + 2 * u / (2 - u))211 // Polynomial approximation of log10(1 + 2 * u / (2 - u))
266 // in [-(2 * fmax) / (2 + fmax), (2 * fmax) / (2 - fmax)]212 // in [-(2 * fmax) / (2 + fmax), (2 * fmax) / (2 - fmax)]
267 // where fmax = 0.5 / size213 // where fmax = 0.5 / size
268 const p11 = 3.8556341504143175800053507804546873e-5;214 const poly: impl.Proc1.Poly = .{
269 const p9 = 1.884958688754320118955531917460363e-4 + v64 * p11;215 .a1 = 0.4342944819032518276511289189166051,
270 const p7 = 9.694073256769014481942040422515466e-4 + v * p9;216 .a3 = 3.619120682527098563759407657655348e-2,
271 const p5 = 5.428681023790647845638954444458386e-3 + v * p7;217 .a5 = 5.428681023790647845638954444458386e-3,
272 const p3 = 3.619120682527098563759407657655348e-2 + v * p5;218 .a7 = 9.694073256769014481942040422515466e-4,
273 const p1 = 0.4342944819032518276511289189166051;219 .a9 = 1.884958688754320118955531917460363e-4,
274220 .a11 = 3.8556341504143175800053507804546873e-5,
275 const q = u * v * p3;221 };
276222
277 // log1p_tab[j].hi = 2^-n * round-to-integer(2^n * l)223 // log1p_tab[j].hi = 2^-n * round-to-integer(2^n * l)
278 // log1p_tab[j].lo = round-to-nearest-f128(l - log1p_tab[j].hi)224 // log1p_tab[j].lo = round-to-nearest-f128(l - log1p_tab[j].hi)
279 // where n = 97 and l = log10(1 + j / size)225 // where n = 97 and l = log10(1 + j / size)
280 const log1p_tab = [size + 1]struct { hi: f128, lo: f128 }{226 const tab = [impl.size + 1]impl.Proc1.HiLo{
281 .{ .hi = 0, .lo = 0 },227 .{ .hi = 0, .lo = 0 },
282 .{ .hi = 0x1.bafd47221ed2665c1ba949p-9, .lo = -0x1.eb6f20a90ad48515635f3b8a1d22p-104 },228 .{ .hi = 0x1.bafd47221ed2665c1ba949p-9, .lo = -0x1.eb6f20a90ad48515635f3b8a1d22p-104 },
283 .{ .hi = 0x1.b9476a4fcd10ed89b5a417p-8, .lo = 0x1.0b153c94bfd2527c3dce31e5e226p-100 },229 .{ .hi = 0x1.b9476a4fcd10ed89b5a417p-8, .lo = 0x1.0b153c94bfd2527c3dce31e5e226p-100 },
...@@ -408,11 +354,7 @@ pub fn log10q(x: f128) callconv(.c) f128 {...@@ -408,11 +354,7 @@ pub fn log10q(x: f128) callconv(.c) f128 {
408 .{ .hi = 0x1.32839e681fc6236e91f3dacap-2, .lo = -0x1.451bc31fd56e57af018a8d364cb8p-99 },354 .{ .hi = 0x1.32839e681fc6236e91f3dacap-2, .lo = -0x1.451bc31fd56e57af018a8d364cb8p-99 },
409 .{ .hi = 0x1.34413509f79fef311f12b358p-2, .lo = 0x1.6f922f04d5a618a87a3e69314bcep-102 },355 .{ .hi = 0x1.34413509f79fef311f12b358p-2, .lo = 0x1.6f922f04d5a618a87a3e69314bcep-102 },
410 };356 };
411 const xm: f128 = @floatFromInt(m);357 return impl.proc1(.{ .poly = poly, .tab = tab }, x);
412 const l_hi = xm * log1p_tab[128].hi + log1p_tab[j].hi;
413 const l_lo = xm * log1p_tab[128].lo + log1p_tab[j].lo;
414
415 return l_hi + (u * p1 + (q + l_lo));
416}358}
417359
418pub fn log10l(x: c_longdouble) callconv(.c) c_longdouble {360pub fn log10l(x: c_longdouble) callconv(.c) c_longdouble {
lib/compiler_rt/log2.zig+30-89
...@@ -176,101 +176,46 @@ pub fn __log2x(a: f80) callconv(.c) f80 {...@@ -176,101 +176,46 @@ pub fn __log2x(a: f80) callconv(.c) f80 {
176/// Accuracy on 10 million random numbers near x = 1 (testing the proc2 case):176/// Accuracy on 10 million random numbers near x = 1 (testing the proc2 case):
177/// <= 0.5 ulp: 99.86%, worst case <= 0.546 ulp177/// <= 0.5 ulp: 99.86%, worst case <= 0.546 ulp
178pub fn log2q(x: f128) callconv(.c) f128 {178pub fn log2q(x: f128) callconv(.c) f128 {
179 const bitsize = 7;179 const impl = @import("log_f128.zig");
180 const size = 1 << bitsize;
181 if (!math.isFinite(x)) {
182 if (math.isNan(x)) {
183 if (math.isSignalNan(x)) math.raiseInvalid();
184 return math.nan(f128);
185 }
186 if (math.isPositiveInf(x)) return x;
187 }
188 if (x <= 0.0) {
189 if (x >= 0.0) {
190 math.raiseDivByZero();
191 return -math.inf(f128);
192 }
193 math.raiseInvalid();
194 return math.nan(f128);
195 }
196180
197 // exp(-1 / 16) rounded down181 if (impl.specialCases(x)) |y|
198 const proc2_lo: f128 = 0.939413062813475786119710824622305;182 return y;
199 // exp(1 / 16) rounded up
200 const proc2_hi: f128 = 1.0644944589178594295633905946428897;
201 if (proc2_lo < x and x < proc2_hi) {
202 const f = x - 1.0;
203 const g = 1 / (2 + f);
204 const u = 2 * f * g;
205 const v = u * u;
206 const uv = u * v;
207 const v64: f64 = @floatCast(v);
208183
184 if (impl.Proc2.lo < x and x < impl.Proc2.hi) {
209 // Polynomial approximation of log2((1 + u / 2) / (1 - u / 2))185 // Polynomial approximation of log2((1 + u / 2) / (1 - u / 2))
210 // in [2 * a / (2 + a), 2 * b / (2 + b)]186 // in [2 * a / (2 + a), 2 * b / (2 + b)]
211 // where a = exp(-1 / 16) - 1 and b = exp(1 / 16) - 1187 // where a = exp(-1 / 16) - 1 and b = exp(1 / 16) - 1
212 const p19 = 2.909291439731657940692470637735429e-7;188 const poly: impl.Proc2.Poly = .{
213 const p17 = 1.294917697032820750200161813672143e-6 + v64 * p19;189 .b1_hi = 0x1.71547652b82fep0,
214 const p15 = 5.870341197214724685339102193694838e-6 + v64 * p17;190 .b1_lo = 0x1.777d0ffda0d23a7d11d6aef551bbp-56,
215 const p13 = 2.7093882228102330360125035037716968e-5 + v * p15;191 .b3 = 0.12022458674074695061332705675016125,
216 const p11 = 1.280801705334664325639770412440281e-4 + v * p13;192 .b5 = 1.8033688011112042591999058475816515e-2,
217 const p9 = 6.261697226080570342191671010883619e-4 + v * p11;193 .b7 = 3.2203014305557218914285331735164364e-3,
218 const p7 = 3.2203014305557218914285331735164364e-3 + v * p9;194 .b9 = 6.261697226080570342191671010883619e-4,
219 const p5 = 1.8033688011112042591999058475816515e-2 + v * p7;195 .b11 = 1.280801705334664325639770412440281e-4,
220 const p3 = 0.12022458674074695061332705675016125;196 .b13 = 2.7093882228102330360125035037716968e-5,
221197 .b15 = 5.870341197214724685339102193694838e-6,
222 const q_hi = uv * p3;198 .b17 = 1.294917697032820750200161813672143e-6,
223 const q_lo = uv * v * p5;199 .b19 = 2.909291439731657940692470637735429e-7,
224200 };
225 const f_hi: f128 = @as(f64, @floatCast(f));201 return impl.proc2(.{ .poly = poly }, x);
226 const f_lo: f128 = f - f_hi;
227
228 const u_hi: f128 = @as(f64, @floatCast(u));
229 const u_lo: f128 = ((2 * (f - u_hi) - u_hi * f_hi) - u_hi * f_lo) * g;
230
231 // t = u / log(2)
232 const log2e_hi: f128 = 0x1.71547652b82fep0;
233 const log2e_lo: f128 = 0x1.777d0ffda0d23a7d11d6aef551bbp-56;
234 const t_hi = u_hi * log2e_hi;
235 const t_lo = u_lo * log2e_hi + u * log2e_lo;
236
237 // y = t + q
238 const y_hi = t_hi + q_hi;
239 const y_lo = t_lo + (t_hi - y_hi + q_hi) + q_lo;
240
241 return y_hi + y_lo;
242 }202 }
243203
244 const ym = @import("log.zig").frexp2(x);
245 const y = ym.significand;
246 const m = ym.exponent;
247
248 const F0 = @round(math.ldexp(y, bitsize));
249 const j0: usize = @intFromFloat(F0);
250 const j = j0 - size;
251 const F = math.ldexp(F0, -bitsize);
252 const f = y - F;
253
254 const u = (f + f) / (y + F);
255 const v = u * u;
256 const v64: f64 = @floatCast(v);
257
258 // Polynomial approximation of log2(1 + 2 * u / (2 - u))204 // Polynomial approximation of log2(1 + 2 * u / (2 - u))
259 // in [-(2 * fmax) / (2 + fmax), (2 * fmax) / (2 - fmax)]205 // in [-(2 * fmax) / (2 + fmax), (2 * fmax) / (2 - fmax)]
260 // where fmax = 0.5 / size206 // where fmax = 0.5 / size
261 const p11 = 1.280813940786848788109850061222256e-4;207 const poly: impl.Proc1.Poly = .{
262 const p9 = 6.261697225875019234719395591078697e-4 + v64 * p11;208 .a1 = 1.442695040888963407359924681001892,
263 const p7 = 3.22030143055572204818095463930704e-3 + v * p9;209 .a3 = 0.12022458674074695061332705675072149,
264 const p5 = 1.8033688011112042591998536830294507e-2 + v * p7;210 .a5 = 1.8033688011112042591998536830294507e-2,
265 const p3 = 0.12022458674074695061332705675072149 + v * p5;211 .a7 = 3.22030143055572204818095463930704e-3,
266 const p1 = 1.442695040888963407359924681001892;212 .a9 = 6.261697225875019234719395591078697e-4,
267213 .a11 = 1.280813940786848788109850061222256e-4,
268 const q = u * v * p3;214 };
269215 // tab[j].hi = 2^-n * round-to-integer(2^n * l)
270 // log1p_tab[j].hi = 2^-n * round-to-integer(2^n * l)216 // tab[j].lo = round-to-nearest-f128(l - tab[j].hi)
271 // log1p_tab[j].lo = round-to-nearest-f128(l - log1p_tab[j].hi)
272 // where n = 97 and l = log2(1 + j / size)217 // where n = 97 and l = log2(1 + j / size)
273 const log1p_tab = [size + 1]struct { hi: f128, lo: f128 }{218 const tab = [impl.size + 1]impl.Proc1.HiLo{
274 .{ .hi = 0, .lo = 0 },219 .{ .hi = 0, .lo = 0 },
275 .{ .hi = 0x1.6fe50b6ef08517f8e37bp-7, .lo = 0x1.794f4441ccdf648f265a41e57d75p-99 },220 .{ .hi = 0x1.6fe50b6ef08517f8e37bp-7, .lo = 0x1.794f4441ccdf648f265a41e57d75p-99 },
276 .{ .hi = 0x1.6e79685c2d2298a6e27e212p-6, .lo = -0x1.fbd41ae7d5a2434912ad3fe21cfbp-100 },221 .{ .hi = 0x1.6e79685c2d2298a6e27e212p-6, .lo = -0x1.fbd41ae7d5a2434912ad3fe21cfbp-100 },
...@@ -401,11 +346,7 @@ pub fn log2q(x: f128) callconv(.c) f128 {...@@ -401,11 +346,7 @@ pub fn log2q(x: f128) callconv(.c) f128 {
401 .{ .hi = 0x1.fd1be4c7f2af942b221ce0d1p-1, .lo = 0x1.a275c854f5bb9732fae5130be48bp-104 },346 .{ .hi = 0x1.fd1be4c7f2af942b221ce0d1p-1, .lo = 0x1.a275c854f5bb9732fae5130be48bp-104 },
402 .{ .hi = 0x1p0, .lo = 0 },347 .{ .hi = 0x1p0, .lo = 0 },
403 };348 };
404 const xm: f128 = @floatFromInt(m);349 return impl.proc1(.{ .poly = poly, .tab = tab }, x);
405 const l_hi = xm * log1p_tab[128].hi + log1p_tab[j].hi;
406 const l_lo = xm * log1p_tab[128].lo + log1p_tab[j].lo;
407
408 return l_hi + (u * p1 + (q + l_lo));
409}350}
410351
411pub fn log2l(x: c_longdouble) callconv(.c) c_longdouble {352pub fn log2l(x: c_longdouble) callconv(.c) c_longdouble {
lib/compiler_rt/log_f128.zig+34
...@@ -1,9 +1,43 @@...@@ -1,9 +1,43 @@
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
1const std = @import("std");11const std = @import("std");
2const math = std.math;12const math = std.math;
313
4pub const log2size = 7;14pub const log2size = 7;
5pub const size = 1 << log2size;15pub const size = 1 << log2size;
616
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
7pub const Proc1 = struct {41pub const Proc1 = struct {
8 pub const Poly = struct {42 pub const Poly = struct {
9 a1: f128,43 a1: f128,