HAXY

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

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

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

1 const builtin = @import("builtin"); 2 3 const std = @import("std"); 4 const Log2Int = std.math.Log2Int; 5 6 const compiler_rt = @import("../compiler_rt.zig"); 7 const symbol = compiler_rt.symbol; 8 const HalveInt = compiler_rt.HalveInt; 9 10 comptime { 11 symbol(&__umodti3, "__umodti3"); 12 symbol(&__modti3, "__modti3"); 13 symbol(&__udivti3, "__udivti3"); 14 symbol(&__divti3, "__divti3"); 15 symbol(&__udivmodti4, "__udivmodti4"); 16 } 17 18 pub fn __udivmodti4(a: u128, b: u128, maybe_rem: ?*u128) callconv(.c) u128 { 19 return udivmod(u128, a, b, maybe_rem); 20 } 21 22 pub fn __divti3(a: i128, b: i128) callconv(.c) i128 { 23 return div(a, b); 24 } 25 26 inline fn div(a: i128, b: i128) i128 { 27 const s_a = a >> (128 - 1); 28 const s_b = b >> (128 - 1); 29 30 const an = (a ^ s_a) -% s_a; 31 const bn = (b ^ s_b) -% s_b; 32 33 const r = udivmod(u128, @bitCast(an), @bitCast(bn), null); 34 const s = s_a ^ s_b; 35 return (@as(i128, @bitCast(r)) ^ s) -% s; 36 } 37 38 pub fn __udivti3(a: u128, b: u128) callconv(.c) u128 { 39 return udivmod(u128, a, b, null); 40 } 41 42 pub fn __modti3(a: i128, b: i128) callconv(.c) i128 { 43 return mod(a, b); 44 } 45 46 inline fn mod(a: i128, b: i128) i128 { 47 const s_a = a >> (128 - 1); // s = a < 0 ? -1 : 0 48 const s_b = b >> (128 - 1); // s = b < 0 ? -1 : 0 49 50 const an = (a ^ s_a) -% s_a; // negate if s == -1 51 const bn = (b ^ s_b) -% s_b; // negate if s == -1 52 53 var r: u128 = undefined; 54 _ = udivmod(u128, @as(u128, @bitCast(an)), @as(u128, @bitCast(bn)), &r); 55 return (@as(i128, @bitCast(r)) ^ s_a) -% s_a; // negate if s == -1 56 } 57 58 pub fn __umodti3(a: u128, b: u128) callconv(.c) u128 { 59 var r: u128 = undefined; 60 _ = udivmod(u128, a, b, &r); 61 return r; 62 } 63 64 // Let _u1 and _u0 be the high and low limbs of U respectively. 65 // Returns U / v_ and sets r = U % v_. 66 fn divwide_generic(comptime T: type, _u1: T, _u0: T, v_: T, r: *T) T { 67 const HalfT = HalveInt(T, false).HalfT; 68 @setRuntimeSafety(compiler_rt.test_safety); 69 var v = v_; 70 71 const b = @as(T, 1) << (@bitSizeOf(T) / 2); 72 var un64: T = undefined; 73 var un10: T = undefined; 74 75 const s: Log2Int(T) = @intCast(@clz(v)); 76 if (s > 0) { 77 // Normalize divisor 78 v <<= s; 79 un64 = (_u1 << s) | (_u0 >> @intCast((@bitSizeOf(T) - @as(T, @intCast(s))))); 80 un10 = _u0 << s; 81 } else { 82 // Avoid undefined behavior of (u0 >> @bitSizeOf(T)) 83 un64 = _u1; 84 un10 = _u0; 85 } 86 87 // Break divisor up into two 32-bit digits 88 const vn1 = v >> (@bitSizeOf(T) / 2); 89 const vn0 = v & std.math.maxInt(HalfT); 90 91 // Break right half of dividend into two digits 92 const un1 = un10 >> (@bitSizeOf(T) / 2); 93 const un0 = un10 & std.math.maxInt(HalfT); 94 95 // Compute the first quotient digit, q1 96 var q1 = un64 / vn1; 97 var rhat = un64 -% q1 *% vn1; 98 99 // q1 has at most error 2. No more than 2 iterations 100 while (q1 >= b or q1 * vn0 > b * rhat + un1) { 101 q1 -= 1; 102 rhat += vn1; 103 if (rhat >= b) break; 104 } 105 106 const un21 = un64 *% b +% un1 -% q1 *% v; 107 108 // Compute the second quotient digit 109 var q0 = un21 / vn1; 110 rhat = un21 -% q0 *% vn1; 111 112 // q0 has at most error 2. No more than 2 iterations. 113 while (q0 >= b or q0 * vn0 > b * rhat + un0) { 114 q0 -= 1; 115 rhat += vn1; 116 if (rhat >= b) break; 117 } 118 119 r.* = (un21 *% b +% un0 -% q0 *% v) >> s; 120 return q1 *% b +% q0; 121 } 122 123 fn divwide(comptime T: type, _u1: T, _u0: T, v: T, r: *T) T { 124 @setRuntimeSafety(compiler_rt.test_safety); 125 if (T == u64 and builtin.target.cpu.arch == .x86_64 and builtin.target.os.tag != .windows) { 126 var rem: T = undefined; 127 const quo = asm ( 128 \\divq %[v] 129 : [_] "={rax}" (-> T), 130 [_] "={rdx}" (rem), 131 : [v] "r" (v), 132 [_] "{rax}" (_u0), 133 [_] "{rdx}" (_u1), 134 ); 135 r.* = rem; 136 return quo; 137 } else { 138 return divwide_generic(T, _u1, _u0, v, r); 139 } 140 } 141 142 // Returns a_ / b_ and sets maybe_rem = a_ % b. 143 pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { 144 @setRuntimeSafety(compiler_rt.test_safety); 145 const HalfT = HalveInt(T, false).HalfT; 146 const half_bits = @bitSizeOf(HalfT); 147 148 if (b_ > a_) { 149 if (maybe_rem) |rem| { 150 rem.* = a_; 151 } 152 return 0; 153 } 154 155 const a: [2]HalfT = @bitCast(a_); // [0] is low bits, [1] is high bits 156 const b: [2]HalfT = @bitCast(b_); // [0] is low bits, [1] is high bits 157 var q: [2]HalfT = undefined; 158 var r: [2]HalfT = undefined; 159 160 // When the divisor fits in 64 bits, we can use an optimized path 161 if (b[1] == 0) { 162 r[1] = 0; 163 if (a[1] < b[0]) { 164 // The result fits in 64 bits 165 q[1] = 0; 166 q[0] = divwide(HalfT, a[1], a[0], b[0], &r[0]); 167 } else { 168 // First, divide with the high part to get the remainder. After that a_hi < b_lo. 169 q[1] = a[1] / b[0]; 170 q[0] = divwide(HalfT, a[1] % b[0], a[0], b[0], &r[0]); 171 } 172 if (maybe_rem) |rem| { 173 rem.* = @bitCast(r); 174 } 175 return @bitCast(q); 176 } 177 178 // Large-divisor case: b[1] != 0, so the quotient fits in one HalfT word. 179 // 180 // Trial quotient via divwide (Knuth Vol 2, Section 4.3.1): 181 // Normalize the divisor so its high half has the MSB set, then use divwide 182 // on the top bits to get a trial quotient that is at most 1 too large. 183 // This replaces the O(shift) bit-by-bit loop with O(1) operations. 184 const s: Log2Int(HalfT) = @intCast(@clz(b[1])); 185 186 if (s == 0) { 187 // b[1] already has its MSB set, so b >= 2^(T_bits - 1). Since a >= b 188 // (we passed the b_ > a_ check), a >= 2^(T_bits - 1) too, meaning 189 // a[1] also has its MSB set. Therefore a / b < 2, and the quotient 190 // is exactly 1. 191 q = @bitCast(@as(T, 0)); 192 q[0] = 1; 193 if (maybe_rem) |rem| { 194 rem.* = a_ - b_; 195 } 196 return @bitCast(q); 197 } 198 199 // Normalize b: shift left by s so bn_hi has its MSB set. 200 const sr: Log2Int(HalfT) = @intCast(half_bits - @as( 201 std.math.IntFittingRange(0, half_bits), 202 @intCast(s), 203 )); 204 const bn_hi: HalfT = (b[1] << s) | (b[0] >> sr); 205 206 // Trial numerator: the top (half_bits + s) bits of (a << s), as [a2:a1]. 207 // a2 < bn_hi is guaranteed since a2 < 2^s and bn_hi >= 2^(half_bits - 1). 208 const a2: HalfT = a[1] >> sr; 209 const a1: HalfT = (a[1] << s) | (a[0] >> sr); 210 211 // Trial quotient via divwide: q_hat = floor([a2:a1] / bn_hi). 212 // By Knuth's theorem (normalized divisor), q <= q_hat <= q + 1. 213 var r_tmp: HalfT = undefined; 214 var q_hat: HalfT = divwide(HalfT, a2, a1, bn_hi, &r_tmp); 215 216 // Verify: q_hat * b must not exceed a. 217 // Compute the product using HalfT * HalfT -> T widening multiplications, 218 // which are native single-instruction ops when HalfT fits in a register 219 // (e.g. u64 * u64 -> u128 via mulq on x86_64, mul on aarch64). 220 // product = q_hat * [b[1]:b[0]] = [p_top : p_mid : p_lo] (3 half-words) 221 const prod_lo: T = @as(T, q_hat) * @as(T, b[0]); 222 const prod_hi: T = @as(T, q_hat) * @as(T, b[1]); 223 224 const prod_lo_parts: [2]HalfT = @bitCast(prod_lo); 225 const prod_hi_parts: [2]HalfT = @bitCast(prod_hi); 226 227 const mid_add = @addWithOverflow(prod_hi_parts[0], prod_lo_parts[1]); 228 var p_mid: HalfT = mid_add[0]; 229 const p_top: HalfT = prod_hi_parts[1] +% @as(HalfT, mid_add[1]); 230 var p_lo: HalfT = prod_lo_parts[0]; 231 232 // If product > a, decrement q_hat (at most once, guaranteed by Knuth). 233 if (p_top > 0 or p_mid > a[1] or (p_mid == a[1] and p_lo > a[0])) { 234 q_hat -= 1; 235 // Subtract b from the product for correct remainder computation. 236 // After correction, (q_hat * b) fits in T bits, so borrows into 237 // p_top cancel it to zero -- we only need [p_mid:p_lo]. 238 const sub_lo = @subWithOverflow(p_lo, b[0]); 239 p_lo = sub_lo[0]; 240 const sub_mid = @subWithOverflow(p_mid, b[1]); 241 const sub_mid2 = @subWithOverflow(sub_mid[0], @as(HalfT, sub_lo[1])); 242 p_mid = sub_mid2[0]; 243 } 244 245 q = @bitCast(@as(T, 0)); 246 q[0] = q_hat; 247 248 if (maybe_rem) |rem| { 249 // remainder = a - q_hat * b = [a[1]:a[0]] - [p_mid:p_lo] 250 // This subtraction is non-negative since q_hat <= true quotient. 251 const rem_lo = @subWithOverflow(a[0], p_lo); 252 r[0] = rem_lo[0]; 253 const rem_hi = @subWithOverflow(a[1], p_mid); 254 const rem_hi2 = @subWithOverflow(rem_hi[0], @as(HalfT, rem_lo[1])); 255 r[1] = rem_hi2[0]; 256 rem.* = @bitCast(r); 257 } 258 return @bitCast(q); 259 } 260 261 test { 262 _ = @import("modti3_test.zig"); 263 _ = @import("divti3_test.zig"); 264 _ = @import("udivmodti4_test.zig"); 265 } 266