HAXY

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

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

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

1 // Ported from musl, which is licensed under the MIT license: 2 // https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT 3 // 4 // https://git.musl-libc.org/cgit/musl/tree/src/math/__rem_pio2l.c 5 6 const std = @import("std"); 7 const math = std.math; 8 const ld = math.long_double; 9 10 const rem_pio2_large = @import("rem_pio2_large.zig").rem_pio2_large; 11 12 pub fn rem_pio2l(comptime T: type, x: T, y: *[2]T) i32 { 13 const impl = switch (T) { 14 f80 => struct { 15 const round1: i8 = 22; 16 const round2: i8 = 61; 17 const nx: i8 = 3; 18 const ny: i8 = 2; 19 20 const pio4: T = 0x1.921fb54442d1846ap-1; 21 // 64 bits of 2/pi 22 const invpio2: T = 6.36619772367581343076e-01; // 0xa2f9836e4e44152a.0p-64 23 // first 39 bits of pi/2 24 const pio2_1: f64 = 1.57079632679597125389e+00; // 0x3FF921FB, 0x54444000 25 // pi/2 - pio2_1 26 const pio2_1t: T = -1.07463465549719416346e-12; // -0x973dcb3b399d747f.0p-103 27 // second 39 bits of pi/2 28 const pio2_2: f64 = -1.07463465549783099519e-12; // -0x12e7b967674000.0p-92 29 // pi/2 - (pio2_1+pio2_2) 30 const pio2_2t: T = 6.36831716351095013979e-25; // 0xc51701b839a25205.0p-144 31 // pi/2 - (pio2_1+pio2_2+pio2_3) 32 const pio2_3t: T = -2.75299651904407171810e-37; // -0xbb5bf6c7ddd660ce.0p-185 33 // third 39 bits of pi/2 34 const pio2_3: f64 = 6.36831716351370313614e-25; // 0x18a2e037074000.0p-133 35 36 fn small(x_val: T) bool { 37 const se = ld.signExponent(x_val); 38 const top = ld.mantissaTop(x_val); 39 const lhs = (@as(u32, se & 0x7fff) << 16) | top; 40 const rhs: u32 = ((0x3fff + 25) << 16) | 0x921f >> 1 | 0x8000; 41 return lhs < rhs; 42 } 43 44 fn quobits(v: T) i32 { 45 const q: i32 = @intFromFloat(v); 46 return @intCast(@as(u32, @bitCast(q)) & 0x7fffffff); 47 } 48 }, 49 f128 => struct { 50 const round1: i8 = 51; 51 const round2: i8 = 119; 52 const nx: i8 = 5; 53 const ny: i8 = 3; 54 55 const pio4: T = 0x1.921fb54442d18469898cc51701b8p-1; 56 const invpio2: T = 6.3661977236758134307553505349005747e-01; 57 const pio2_1: T = 1.5707963267948966192292994253909555e+00; 58 const pio2_1t: T = 2.0222662487959507323996846200947577e-21; 59 const pio2_2: T = 2.0222662487959507323994779168837751e-21; 60 const pio2_2t: T = 2.0670321098263988236496903051604844e-43; 61 const pio2_3: T = 2.0670321098263988236499468110329591e-43; 62 const pio2_3t: T = -2.5650587247459238361625433492959285e-65; 63 64 fn small(x_val: T) bool { 65 const se = ld.signExponent(x_val); 66 const top = ld.mantissaTop(x_val); 67 const lhs = (@as(u32, se & 0x7fff) << 16) | top; 68 const rhs: u32 = ((0x3fff + 45) << 16) | 0x921f; 69 return lhs < rhs; 70 } 71 72 fn quobits(fn_val: T) i32 { 73 const q: i64 = @intFromFloat(fn_val); 74 return @intCast(@as(u64, @bitCast(q)) & 0x7fffffff); 75 } 76 }, 77 else => @compileError("rem_pio2l supports only f80 and f128, got: " ++ @typeName(T)), 78 }; 79 80 const x_se = ld.signExponent(x); 81 const ex: i32 = @intCast(x_se & 0x7fff); 82 83 if (impl.small(x)) { 84 // rint(x/(pi/2)) 85 const toint: T = 1.5 / math.floatEps(T); 86 var fn_ = x * impl.invpio2 + toint - toint; 87 var n = impl.quobits(fn_); 88 var r = x - fn_ * @as(T, impl.pio2_1); 89 var w = fn_ * impl.pio2_1t; // 1st round good to 102/180 bits 90 91 // Matters with directed rounding. 92 if (r - w < -impl.pio4) { 93 @branchHint(.unlikely); 94 n -= 1; 95 fn_ -= 1; 96 r = x - fn_ * @as(T, impl.pio2_1); 97 w = fn_ * impl.pio2_1t; 98 } else if (r - w > impl.pio4) { 99 @branchHint(.unlikely); 100 n += 1; 101 fn_ += 1; 102 r = x - fn_ * @as(T, impl.pio2_1); 103 w = fn_ * impl.pio2_1t; 104 } 105 106 y[0] = r - w; 107 108 const ey: i32 = @intCast(ld.signExponent(y[0]) & 0x7fff); 109 if (ex - ey > impl.round1) { 110 var t = r; 111 w = fn_ * impl.pio2_2; 112 r = t - w; 113 w = fn_ * impl.pio2_2t - ((t - r) - w); 114 y[0] = r - w; 115 const ey2: i32 = @intCast(ld.signExponent(y[0]) & 0x7fff); 116 if (ex - ey2 > impl.round2) { 117 t = r; 118 w = fn_ * impl.pio2_3; 119 r = t - w; 120 w = fn_ * impl.pio2_3t - ((t - r) - w); 121 y[0] = r - w; 122 } 123 } 124 y[1] = (r - y[0]) - w; 125 return n; 126 } 127 128 // all other (large) arguments 129 if (ex == 0x7fff) { // x is inf or NaN 130 y[0] = x - x; 131 y[1] = y[0]; 132 return 0; 133 } 134 135 var z: T = math.scalbn(@abs(x), -math.ilogb(x) + 23); 136 var tx: [impl.nx]f64 = undefined; 137 var ty: [impl.ny]f64 = undefined; 138 var i: usize = 0; 139 140 while (i < impl.nx - 1) : (i += 1) { 141 tx[i] = @floatFromInt(@as(i32, @intFromFloat(z))); 142 z = (z - @as(T, tx[i])) * 0x1p24; 143 } 144 145 tx[i] = @floatCast(z); 146 while (tx[i] == 0.0) { 147 i -= 1; 148 } 149 150 const n = rem_pio2_large( 151 tx[0..(i + 1)], 152 ty[0..impl.ny], 153 ex - 0x3fff - 23, 154 @intCast(i + 1), 155 impl.ny, 156 ); 157 var w: f64 = ty[1]; 158 if (impl.ny == 3) { 159 w += ty[2]; 160 } 161 const r = ty[0] + w; 162 w -= r - ty[0]; 163 164 if (x_se >> 15 != 0) { 165 y[0] = -@as(T, r); 166 y[1] = -@as(T, w); 167 return -n; 168 } 169 170 y[0] = @as(T, r); 171 y[1] = @as(T, w); 172 return n; 173 } 174