| ... | @@ -182,7 +182,7 @@ fn divwide(comptime T: type, _u1: T, _u0: T, v: T, r: *T) T { | ... | @@ -182,7 +182,7 @@ fn divwide(comptime T: type, _u1: T, _u0: T, v: T, r: *T) T { |
| 182 | pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { | 182 | pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { |
| 183 | @setRuntimeSafety(compiler_rt.test_safety); | 183 | @setRuntimeSafety(compiler_rt.test_safety); |
| 184 | const HalfT = HalveInt(T, false).HalfT; | 184 | const HalfT = HalveInt(T, false).HalfT; |
| 185 | const SignedT = std.meta.Int(.signed, @bitSizeOf(T)); | 185 | const half_bits = @bitSizeOf(HalfT); |
| 186 | | 186 | |
| 187 | if (b_ > a_) { | 187 | if (b_ > a_) { |
| 188 | if (maybe_rem) |rem| { | 188 | if (maybe_rem) |rem| { |
| ... | @@ -214,26 +214,85 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { | ... | @@ -214,26 +214,85 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { |
| 214 | return @bitCast(q); | 214 | return @bitCast(q); |
| 215 | } | 215 | } |
| 216 | | 216 | |
| 217 | // 0 <= shift <= 63 | 217 | // Large-divisor case: b[hi] != 0, so the quotient fits in one HalfT word. |
| 218 | const shift: Log2Int(T) = @clz(b[hi]) - @clz(a[hi]); | 218 | // |
| 219 | var af: T = @bitCast(a); | 219 | // Trial quotient via divwide (Knuth Vol 2, Section 4.3.1): |
| 220 | var bf = @as(T, @bitCast(b)) << shift; | 220 | // Normalize the divisor so its high half has the MSB set, then use divwide |
| 221 | q = @bitCast(@as(T, 0)); | 221 | // on the top bits to get a trial quotient that is at most 1 too large. |
| | 222 | // This replaces the O(shift) bit-by-bit loop with O(1) operations. |
| | 223 | const s: Log2Int(HalfT) = @intCast(@clz(b[hi])); |
| | 224 | |
| | 225 | if (s == 0) { |
| | 226 | // b[hi] already has its MSB set, so b >= 2^(T_bits - 1). Since a >= b |
| | 227 | // (we passed the b_ > a_ check), a >= 2^(T_bits - 1) too, meaning |
| | 228 | // a[hi] also has its MSB set. Therefore a / b < 2, and the quotient |
| | 229 | // is exactly 1. |
| | 230 | q = @bitCast(@as(T, 0)); |
| | 231 | q[lo] = 1; |
| | 232 | if (maybe_rem) |rem| { |
| | 233 | rem.* = a_ - b_; |
| | 234 | } |
| | 235 | return @bitCast(q); |
| | 236 | } |
| 222 | | 237 | |
| 223 | for (0..shift + 1) |_| { | 238 | // Normalize b: shift left by s so bn_hi has its MSB set. |
| 224 | q[lo] <<= 1; | 239 | const sr: Log2Int(HalfT) = @intCast(half_bits - @as( |
| 225 | // Branchless version of: | 240 | std.math.IntFittingRange(0, half_bits), |
| 226 | // if (af >= bf) { | 241 | @intCast(s), |
| 227 | // af -= bf; | 242 | )); |
| 228 | // q[lo] |= 1; | 243 | const bn_hi: HalfT = (b[hi] << s) | (b[lo] >> sr); |
| 229 | // } | 244 | |
| 230 | const s = @as(SignedT, @bitCast(bf -% af -% 1)) >> (@bitSizeOf(T) - 1); | 245 | // Trial numerator: the top (half_bits + s) bits of (a << s), as [a2:a1]. |
| 231 | q[lo] |= @intCast(s & 1); | 246 | // a2 < bn_hi is guaranteed since a2 < 2^s and bn_hi >= 2^(half_bits - 1). |
| 232 | af -= bf & @as(T, @bitCast(s)); | 247 | const a2: HalfT = a[hi] >> sr; |
| 233 | bf >>= 1; | 248 | const a1: HalfT = (a[hi] << s) | (a[lo] >> sr); |
| | 249 | |
| | 250 | // Trial quotient via divwide: q_hat = floor([a2:a1] / bn_hi). |
| | 251 | // By Knuth's theorem (normalized divisor), q <= q_hat <= q + 1. |
| | 252 | var r_tmp: HalfT = undefined; |
| | 253 | var q_hat: HalfT = divwide(HalfT, a2, a1, bn_hi, &r_tmp); |
| | 254 | |
| | 255 | // Verify: q_hat * b must not exceed a. |
| | 256 | // Compute the product using HalfT * HalfT -> T widening multiplications, |
| | 257 | // which are native single-instruction ops when HalfT fits in a register |
| | 258 | // (e.g. u64 * u64 -> u128 via mulq on x86_64, mul on aarch64). |
| | 259 | // product = q_hat * [b[hi]:b[lo]] = [p_top : p_mid : p_lo] (3 half-words) |
| | 260 | const prod_lo: T = @as(T, q_hat) * @as(T, b[lo]); |
| | 261 | const prod_hi: T = @as(T, q_hat) * @as(T, b[hi]); |
| | 262 | |
| | 263 | const prod_lo_parts: [2]HalfT = @bitCast(prod_lo); |
| | 264 | const prod_hi_parts: [2]HalfT = @bitCast(prod_hi); |
| | 265 | |
| | 266 | const mid_add = @addWithOverflow(prod_hi_parts[lo], prod_lo_parts[hi]); |
| | 267 | var p_mid: HalfT = mid_add[0]; |
| | 268 | const p_top: HalfT = prod_hi_parts[hi] +% @as(HalfT, mid_add[1]); |
| | 269 | var p_lo: HalfT = prod_lo_parts[lo]; |
| | 270 | |
| | 271 | // If product > a, decrement q_hat (at most once, guaranteed by Knuth). |
| | 272 | if (p_top > 0 or p_mid > a[hi] or (p_mid == a[hi] and p_lo > a[lo])) { |
| | 273 | q_hat -= 1; |
| | 274 | // Subtract b from the product for correct remainder computation. |
| | 275 | // After correction, (q_hat * b) fits in T bits, so borrows into |
| | 276 | // p_top cancel it to zero -- we only need [p_mid:p_lo]. |
| | 277 | const sub_lo = @subWithOverflow(p_lo, b[lo]); |
| | 278 | p_lo = sub_lo[0]; |
| | 279 | const sub_mid = @subWithOverflow(p_mid, b[hi]); |
| | 280 | const sub_mid2 = @subWithOverflow(sub_mid[0], @as(HalfT, sub_lo[1])); |
| | 281 | p_mid = sub_mid2[0]; |
| 234 | } | 282 | } |
| | 283 | |
| | 284 | q = @bitCast(@as(T, 0)); |
| | 285 | q[lo] = q_hat; |
| | 286 | |
| 235 | if (maybe_rem) |rem| { | 287 | if (maybe_rem) |rem| { |
| 236 | rem.* = @bitCast(af); | 288 | // remainder = a - q_hat * b = [a[hi]:a[lo]] - [p_mid:p_lo] |
| | 289 | // This subtraction is non-negative since q_hat <= true quotient. |
| | 290 | const rem_lo = @subWithOverflow(a[lo], p_lo); |
| | 291 | r[lo] = rem_lo[0]; |
| | 292 | const rem_hi = @subWithOverflow(a[hi], p_mid); |
| | 293 | const rem_hi2 = @subWithOverflow(rem_hi[0], @as(HalfT, rem_lo[1])); |
| | 294 | r[hi] = rem_hi2[0]; |
| | 295 | rem.* = @bitCast(r); |
| 237 | } | 296 | } |
| 238 | return @bitCast(q); | 297 | return @bitCast(q); |
| 239 | } | 298 | } |
Koko Bhadra <amlandeep1912@gmail.com> 2026-03-03 10:15:26-05:00
Andrew Kelley <andrew@ziglang.org> 2026-03-05 20:22:19+01:00
lib