HAXY

▄▄ ▀▀▀▄ ▄▄ ▀▀▀▀▀▄ ▄▀▀▀▀ ▀▀▀▀▀▀ ▄▀▀▀▀▀▀ ▀▀▀▀▀▀ ▀ ▀ ▄▀▀▀▀▀▀ ▀▀▀▀▀▀ ▄▀ ▄▀ ▄▀▀▀▀▀▀▀▀ ▀▀▀▄▄ ▀▀▀▀▀▀ ▀▀ ▀▄▀▀ ▄▀ ▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▄ ▀▀▀▀▀▄ ▄▀▀▀▀▀▀▀ ▄ ▄▀▄▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▄ ▄▀▀▀▀▀▀▀▀▀▀▄ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▄ ▀ ▀▀▀▀▀▀ ▄▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀ ▄▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▄▀▀▀▀▀▀▀▀▀▀▀▀▀▀▄ ▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▄ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▄▄▄▀▀▀▀▄▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▄▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▄ ▀▀▀▀▀▀▀▀▀▀▀▀▀▄▀▀▀▀▀▀▀▀▀▀▄ ▄▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▄▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▄▄ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▄▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀ ▀▀▀▀▀▀▀▀▄ ▀▀▀ ▀▀▀▀▀▀▀▀ ▀▀▀ ▀▀▀ ▀▀▀▀▀ ▀▀ ▀▀▀▀▀▀ ▀▀▀▄▀▄ ▄▄▄▀▀ ▀▀▀▀▀▀ ▀▀▄▀▀▀▀▀▀▀ ▀▀▀▀▀▄▀▀▀▀▀ ▀▀▀▀▀▄ ▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▄▀ ▄▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀▀▀▀▀▀▀▀▀▀ ▀▀▀ ▀

╭ log_f128.zig ─────────────────────────────────────────────────╮ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ ╰───────────────────────────────────────────────────────────────╯

╭ find file ─ ─ ─ ─ ─ ─ ─ ─ ─ ─╮ ╭────╮ ╭ clone ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ╮ │ │ │http│ ssh │ │ ╰ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─╯ ╰────╯ ╰ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ╯

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. 9 const log_f128 = @This(); 10 11 const std = @import("std"); 12 const math = std.math; 13 14 pub const log2size = 7; 15 pub 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. 21 pub 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 41 pub 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 55 pub 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 87 pub 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 110 pub 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. 153 pub 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 } 174