| ... | @@ -61,12 +61,6 @@ pub fn __umodti3(a: u128, b: u128) callconv(.c) u128 { | ... | @@ -61,12 +61,6 @@ pub fn __umodti3(a: u128, b: u128) callconv(.c) u128 { |
| 61 | return r; | 61 | return r; |
| 62 | } | 62 | } |
| 63 | | 63 | |
| 64 | const lo = switch (builtin.cpu.arch.endian()) { | | |
| 65 | .big => 1, | | |
| 66 | .little => 0, | | |
| 67 | }; | | |
| 68 | const hi = 1 - lo; | | |
| 69 | | | |
| 70 | // Let _u1 and _u0 be the high and low limbs of U respectively. | 64 | // Let _u1 and _u0 be the high and low limbs of U respectively. |
| 71 | // Returns U / v_ and sets r = U % v_. | 65 | // Returns U / v_ and sets r = U % v_. |
| 72 | fn divwide_generic(comptime T: type, _u1: T, _u0: T, v_: T, r: *T) T { | 66 | fn divwide_generic(comptime T: type, _u1: T, _u0: T, v_: T, r: *T) T { |
| ... | @@ -158,22 +152,22 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { | ... | @@ -158,22 +152,22 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { |
| 158 | return 0; | 152 | return 0; |
| 159 | } | 153 | } |
| 160 | | 154 | |
| 161 | const a: [2]HalfT = @bitCast(a_); | 155 | const a: [2]HalfT = @bitCast(a_); // [0] is low bits, [1] is high bits |
| 162 | const b: [2]HalfT = @bitCast(b_); | 156 | const b: [2]HalfT = @bitCast(b_); // [0] is low bits, [1] is high bits |
| 163 | var q: [2]HalfT = undefined; | 157 | var q: [2]HalfT = undefined; |
| 164 | var r: [2]HalfT = undefined; | 158 | var r: [2]HalfT = undefined; |
| 165 | | 159 | |
| 166 | // When the divisor fits in 64 bits, we can use an optimized path | 160 | // When the divisor fits in 64 bits, we can use an optimized path |
| 167 | if (b[hi] == 0) { | 161 | if (b[1] == 0) { |
| 168 | r[hi] = 0; | 162 | r[1] = 0; |
| 169 | if (a[hi] < b[lo]) { | 163 | if (a[1] < b[0]) { |
| 170 | // The result fits in 64 bits | 164 | // The result fits in 64 bits |
| 171 | q[hi] = 0; | 165 | q[1] = 0; |
| 172 | q[lo] = divwide(HalfT, a[hi], a[lo], b[lo], &r[lo]); | 166 | q[0] = divwide(HalfT, a[1], a[0], b[0], &r[0]); |
| 173 | } else { | 167 | } else { |
| 174 | // First, divide with the high part to get the remainder. After that a_hi < b_lo. | 168 | // First, divide with the high part to get the remainder. After that a_hi < b_lo. |
| 175 | q[hi] = a[hi] / b[lo]; | 169 | q[1] = a[1] / b[0]; |
| 176 | q[lo] = divwide(HalfT, a[hi] % b[lo], a[lo], b[lo], &r[lo]); | 170 | q[0] = divwide(HalfT, a[1] % b[0], a[0], b[0], &r[0]); |
| 177 | } | 171 | } |
| 178 | if (maybe_rem) |rem| { | 172 | if (maybe_rem) |rem| { |
| 179 | rem.* = @bitCast(r); | 173 | rem.* = @bitCast(r); |
| ... | @@ -181,21 +175,21 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { | ... | @@ -181,21 +175,21 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { |
| 181 | return @bitCast(q); | 175 | return @bitCast(q); |
| 182 | } | 176 | } |
| 183 | | 177 | |
| 184 | // Large-divisor case: b[hi] != 0, so the quotient fits in one HalfT word. | 178 | // Large-divisor case: b[1] != 0, so the quotient fits in one HalfT word. |
| 185 | // | 179 | // |
| 186 | // Trial quotient via divwide (Knuth Vol 2, Section 4.3.1): | 180 | // Trial quotient via divwide (Knuth Vol 2, Section 4.3.1): |
| 187 | // Normalize the divisor so its high half has the MSB set, then use divwide | 181 | // Normalize the divisor so its high half has the MSB set, then use divwide |
| 188 | // on the top bits to get a trial quotient that is at most 1 too large. | 182 | // on the top bits to get a trial quotient that is at most 1 too large. |
| 189 | // This replaces the O(shift) bit-by-bit loop with O(1) operations. | 183 | // This replaces the O(shift) bit-by-bit loop with O(1) operations. |
| 190 | const s: Log2Int(HalfT) = @intCast(@clz(b[hi])); | 184 | const s: Log2Int(HalfT) = @intCast(@clz(b[1])); |
| 191 | | 185 | |
| 192 | if (s == 0) { | 186 | if (s == 0) { |
| 193 | // b[hi] already has its MSB set, so b >= 2^(T_bits - 1). Since a >= b | 187 | // b[1] already has its MSB set, so b >= 2^(T_bits - 1). Since a >= b |
| 194 | // (we passed the b_ > a_ check), a >= 2^(T_bits - 1) too, meaning | 188 | // (we passed the b_ > a_ check), a >= 2^(T_bits - 1) too, meaning |
| 195 | // a[hi] also has its MSB set. Therefore a / b < 2, and the quotient | 189 | // a[1] also has its MSB set. Therefore a / b < 2, and the quotient |
| 196 | // is exactly 1. | 190 | // is exactly 1. |
| 197 | q = @bitCast(@as(T, 0)); | 191 | q = @bitCast(@as(T, 0)); |
| 198 | q[lo] = 1; | 192 | q[0] = 1; |
| 199 | if (maybe_rem) |rem| { | 193 | if (maybe_rem) |rem| { |
| 200 | rem.* = a_ - b_; | 194 | rem.* = a_ - b_; |
| 201 | } | 195 | } |
| ... | @@ -207,12 +201,12 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { | ... | @@ -207,12 +201,12 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { |
| 207 | std.math.IntFittingRange(0, half_bits), | 201 | std.math.IntFittingRange(0, half_bits), |
| 208 | @intCast(s), | 202 | @intCast(s), |
| 209 | )); | 203 | )); |
| 210 | const bn_hi: HalfT = (b[hi] << s) | (b[lo] >> sr); | 204 | const bn_hi: HalfT = (b[1] << s) | (b[0] >> sr); |
| 211 | | 205 | |
| 212 | // Trial numerator: the top (half_bits + s) bits of (a << s), as [a2:a1]. | 206 | // Trial numerator: the top (half_bits + s) bits of (a << s), as [a2:a1]. |
| 213 | // a2 < bn_hi is guaranteed since a2 < 2^s and bn_hi >= 2^(half_bits - 1). | 207 | // a2 < bn_hi is guaranteed since a2 < 2^s and bn_hi >= 2^(half_bits - 1). |
| 214 | const a2: HalfT = a[hi] >> sr; | 208 | const a2: HalfT = a[1] >> sr; |
| 215 | const a1: HalfT = (a[hi] << s) | (a[lo] >> sr); | 209 | const a1: HalfT = (a[1] << s) | (a[0] >> sr); |
| 216 | | 210 | |
| 217 | // Trial quotient via divwide: q_hat = floor([a2:a1] / bn_hi). | 211 | // Trial quotient via divwide: q_hat = floor([a2:a1] / bn_hi). |
| 218 | // By Knuth's theorem (normalized divisor), q <= q_hat <= q + 1. | 212 | // By Knuth's theorem (normalized divisor), q <= q_hat <= q + 1. |
| ... | @@ -223,42 +217,42 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { | ... | @@ -223,42 +217,42 @@ pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { |
| 223 | // Compute the product using HalfT * HalfT -> T widening multiplications, | 217 | // Compute the product using HalfT * HalfT -> T widening multiplications, |
| 224 | // which are native single-instruction ops when HalfT fits in a register | 218 | // which are native single-instruction ops when HalfT fits in a register |
| 225 | // (e.g. u64 * u64 -> u128 via mulq on x86_64, mul on aarch64). | 219 | // (e.g. u64 * u64 -> u128 via mulq on x86_64, mul on aarch64). |
| 226 | // product = q_hat * [b[hi]:b[lo]] = [p_top : p_mid : p_lo] (3 half-words) | 220 | // product = q_hat * [b[1]:b[0]] = [p_top : p_mid : p_lo] (3 half-words) |
| 227 | const prod_lo: T = @as(T, q_hat) * @as(T, b[lo]); | 221 | const prod_lo: T = @as(T, q_hat) * @as(T, b[0]); |
| 228 | const prod_hi: T = @as(T, q_hat) * @as(T, b[hi]); | 222 | const prod_hi: T = @as(T, q_hat) * @as(T, b[1]); |
| 229 | | 223 | |
| 230 | const prod_lo_parts: [2]HalfT = @bitCast(prod_lo); | 224 | const prod_lo_parts: [2]HalfT = @bitCast(prod_lo); |
| 231 | const prod_hi_parts: [2]HalfT = @bitCast(prod_hi); | 225 | const prod_hi_parts: [2]HalfT = @bitCast(prod_hi); |
| 232 | | 226 | |
| 233 | const mid_add = @addWithOverflow(prod_hi_parts[lo], prod_lo_parts[hi]); | 227 | const mid_add = @addWithOverflow(prod_hi_parts[0], prod_lo_parts[1]); |
| 234 | var p_mid: HalfT = mid_add[0]; | 228 | var p_mid: HalfT = mid_add[0]; |
| 235 | const p_top: HalfT = prod_hi_parts[hi] +% @as(HalfT, mid_add[1]); | 229 | const p_top: HalfT = prod_hi_parts[1] +% @as(HalfT, mid_add[1]); |
| 236 | var p_lo: HalfT = prod_lo_parts[lo]; | 230 | var p_lo: HalfT = prod_lo_parts[0]; |
| 237 | | 231 | |
| 238 | // If product > a, decrement q_hat (at most once, guaranteed by Knuth). | 232 | // If product > a, decrement q_hat (at most once, guaranteed by Knuth). |
| 239 | if (p_top > 0 or p_mid > a[hi] or (p_mid == a[hi] and p_lo > a[lo])) { | 233 | if (p_top > 0 or p_mid > a[1] or (p_mid == a[1] and p_lo > a[0])) { |
| 240 | q_hat -= 1; | 234 | q_hat -= 1; |
| 241 | // Subtract b from the product for correct remainder computation. | 235 | // Subtract b from the product for correct remainder computation. |
| 242 | // After correction, (q_hat * b) fits in T bits, so borrows into | 236 | // After correction, (q_hat * b) fits in T bits, so borrows into |
| 243 | // p_top cancel it to zero -- we only need [p_mid:p_lo]. | 237 | // p_top cancel it to zero -- we only need [p_mid:p_lo]. |
| 244 | const sub_lo = @subWithOverflow(p_lo, b[lo]); | 238 | const sub_lo = @subWithOverflow(p_lo, b[0]); |
| 245 | p_lo = sub_lo[0]; | 239 | p_lo = sub_lo[0]; |
| 246 | const sub_mid = @subWithOverflow(p_mid, b[hi]); | 240 | const sub_mid = @subWithOverflow(p_mid, b[1]); |
| 247 | const sub_mid2 = @subWithOverflow(sub_mid[0], @as(HalfT, sub_lo[1])); | 241 | const sub_mid2 = @subWithOverflow(sub_mid[0], @as(HalfT, sub_lo[1])); |
| 248 | p_mid = sub_mid2[0]; | 242 | p_mid = sub_mid2[0]; |
| 249 | } | 243 | } |
| 250 | | 244 | |
| 251 | q = @bitCast(@as(T, 0)); | 245 | q = @bitCast(@as(T, 0)); |
| 252 | q[lo] = q_hat; | 246 | q[0] = q_hat; |
| 253 | | 247 | |
| 254 | if (maybe_rem) |rem| { | 248 | if (maybe_rem) |rem| { |
| 255 | // remainder = a - q_hat * b = [a[hi]:a[lo]] - [p_mid:p_lo] | 249 | // remainder = a - q_hat * b = [a[1]:a[0]] - [p_mid:p_lo] |
| 256 | // This subtraction is non-negative since q_hat <= true quotient. | 250 | // This subtraction is non-negative since q_hat <= true quotient. |
| 257 | const rem_lo = @subWithOverflow(a[lo], p_lo); | 251 | const rem_lo = @subWithOverflow(a[0], p_lo); |
| 258 | r[lo] = rem_lo[0]; | 252 | r[0] = rem_lo[0]; |
| 259 | const rem_hi = @subWithOverflow(a[hi], p_mid); | 253 | const rem_hi = @subWithOverflow(a[1], p_mid); |
| 260 | const rem_hi2 = @subWithOverflow(rem_hi[0], @as(HalfT, rem_lo[1])); | 254 | const rem_hi2 = @subWithOverflow(rem_hi[0], @as(HalfT, rem_lo[1])); |
| 261 | r[hi] = rem_hi2[0]; | 255 | r[1] = rem_hi2[0]; |
| 262 | rem.* = @bitCast(r); | 256 | rem.* = @bitCast(r); |
| 263 | } | 257 | } |
| 264 | return @bitCast(q); | 258 | return @bitCast(q); |