| ... | @@ -1,201 +1,149 @@ | ... | @@ -1,201 +1,149 @@ |
| | 1 | const std = @import("std"); |
| 1 | const builtin = @import("builtin"); | 2 | const builtin = @import("builtin"); |
| 2 | const is_test = builtin.is_test; | 3 | const is_test = builtin.is_test; |
| 3 | const native_endian = builtin.cpu.arch.endian(); | 4 | const Log2Int = std.math.Log2Int; |
| 4 | const std = @import("std"); | 5 | const HalveInt = @import("common.zig").HalveInt; |
| 5 | | 6 | |
| 6 | const low = switch (native_endian) { | 7 | const lo = switch (builtin.cpu.arch.endian()) { |
| 7 | .Big => 1, | 8 | .Big => 1, |
| 8 | .Little => 0, | 9 | .Little => 0, |
| 9 | }; | 10 | }; |
| 10 | const high = 1 - low; | 11 | const hi = 1 - lo; |
| 11 | | 12 | |
| 12 | pub fn udivmod(comptime DoubleInt: type, a: DoubleInt, b: DoubleInt, maybe_rem: ?*DoubleInt) DoubleInt { | 13 | // Let _u1 and _u0 be the high and low limbs of U respectively. |
| | 14 | // Returns U / v_ and sets r = U % v_. |
| | 15 | fn divwide_generic(comptime T: type, _u1: T, _u0: T, v_: T, r: *T) T { |
| | 16 | const HalfT = HalveInt(T, false).HalfT; |
| 13 | @setRuntimeSafety(is_test); | 17 | @setRuntimeSafety(is_test); |
| | 18 | var v = v_; |
| 14 | | 19 | |
| 15 | const double_int_bits = @typeInfo(DoubleInt).Int.bits; | 20 | const b = @as(T, 1) << (@bitSizeOf(T) / 2); |
| 16 | const single_int_bits = @divExact(double_int_bits, 2); | 21 | var un64: T = undefined; |
| 17 | const SingleInt = std.meta.Int(.unsigned, single_int_bits); | 22 | var un10: T = undefined; |
| 18 | const SignedDoubleInt = std.meta.Int(.signed, double_int_bits); | 23 | |
| 19 | const Log2SingleInt = std.math.Log2Int(SingleInt); | 24 | const s = @intCast(Log2Int(T), @clz(v)); |
| 20 | | 25 | if (s > 0) { |
| 21 | const n = @bitCast([2]SingleInt, a); | 26 | // Normalize divisor |
| 22 | const d = @bitCast([2]SingleInt, b); | 27 | v <<= s; |
| 23 | var q: [2]SingleInt = undefined; | 28 | un64 = (_u1 << s) | (_u0 >> @intCast(Log2Int(T), (@bitSizeOf(T) - @intCast(T, s)))); |
| 24 | var r: [2]SingleInt = undefined; | 29 | un10 = _u0 << s; |
| 25 | var sr: c_uint = undefined; | 30 | } else { |
| 26 | // special cases, X is unknown, K != 0 | 31 | // Avoid undefined behavior of (u0 >> @bitSizeOf(T)) |
| 27 | if (n[high] == 0) { | 32 | un64 = _u1; |
| 28 | if (d[high] == 0) { | 33 | un10 = _u0; |
| 29 | // 0 X | 34 | } |
| 30 | // --- | 35 | |
| 31 | // 0 X | 36 | // Break divisor up into two 32-bit digits |
| 32 | if (maybe_rem) |rem| { | 37 | const vn1 = v >> (@bitSizeOf(T) / 2); |
| 33 | rem.* = n[low] % d[low]; | 38 | const vn0 = v & std.math.maxInt(HalfT); |
| 34 | } | 39 | |
| 35 | return n[low] / d[low]; | 40 | // Break right half of dividend into two digits |
| 36 | } | 41 | const un1 = un10 >> (@bitSizeOf(T) / 2); |
| 37 | // 0 X | 42 | const un0 = un10 & std.math.maxInt(HalfT); |
| 38 | // --- | 43 | |
| 39 | // K X | 44 | // Compute the first quotient digit, q1 |
| | 45 | var q1 = un64 / vn1; |
| | 46 | var rhat = un64 -% q1 *% vn1; |
| | 47 | |
| | 48 | // q1 has at most error 2. No more than 2 iterations |
| | 49 | while (q1 >= b or q1 * vn0 > b * rhat + un1) { |
| | 50 | q1 -= 1; |
| | 51 | rhat += vn1; |
| | 52 | if (rhat >= b) break; |
| | 53 | } |
| | 54 | |
| | 55 | var un21 = un64 *% b +% un1 -% q1 *% v; |
| | 56 | |
| | 57 | // Compute the second quotient digit |
| | 58 | var q0 = un21 / vn1; |
| | 59 | rhat = un21 -% q0 *% vn1; |
| | 60 | |
| | 61 | // q0 has at most error 2. No more than 2 iterations. |
| | 62 | while (q0 >= b or q0 * vn0 > b * rhat + un0) { |
| | 63 | q0 -= 1; |
| | 64 | rhat += vn1; |
| | 65 | if (rhat >= b) break; |
| | 66 | } |
| | 67 | |
| | 68 | r.* = (un21 *% b +% un0 -% q0 *% v) >> s; |
| | 69 | return q1 *% b +% q0; |
| | 70 | } |
| | 71 | |
| | 72 | fn divwide(comptime T: type, _u1: T, _u0: T, v: T, r: *T) T { |
| | 73 | @setRuntimeSafety(is_test); |
| | 74 | if (T == u64 and builtin.target.cpu.arch == .x86_64 and builtin.target.os.tag != .windows) { |
| | 75 | var rem: T = undefined; |
| | 76 | const quo = asm ( |
| | 77 | \\divq %[v] |
| | 78 | : [_] "={rax}" (-> T), |
| | 79 | [_] "={rdx}" (rem), |
| | 80 | : [v] "r" (v), |
| | 81 | [_] "{rax}" (_u0), |
| | 82 | [_] "{rdx}" (_u1), |
| | 83 | ); |
| | 84 | r.* = rem; |
| | 85 | return quo; |
| | 86 | } else { |
| | 87 | return divwide_generic(T, _u1, _u0, v, r); |
| | 88 | } |
| | 89 | } |
| | 90 | |
| | 91 | // Returns a_ / b_ and sets maybe_rem = a_ % b. |
| | 92 | pub fn udivmod(comptime T: type, a_: T, b_: T, maybe_rem: ?*T) T { |
| | 93 | @setRuntimeSafety(is_test); |
| | 94 | const HalfT = HalveInt(T, false).HalfT; |
| | 95 | const SignedT = std.meta.Int(.signed, @bitSizeOf(T)); |
| | 96 | |
| | 97 | if (b_ > a_) { |
| 40 | if (maybe_rem) |rem| { | 98 | if (maybe_rem) |rem| { |
| 41 | rem.* = n[low]; | 99 | rem.* = a_; |
| 42 | } | 100 | } |
| 43 | return 0; | 101 | return 0; |
| 44 | } | 102 | } |
| 45 | // n[high] != 0 | 103 | |
| 46 | if (d[low] == 0) { | 104 | var a = @bitCast([2]HalfT, a_); |
| 47 | if (d[high] == 0) { | 105 | var b = @bitCast([2]HalfT, b_); |
| 48 | // K X | 106 | var q: [2]HalfT = undefined; |
| 49 | // --- | 107 | var r: [2]HalfT = undefined; |
| 50 | // 0 0 | 108 | |
| 51 | if (maybe_rem) |rem| { | 109 | // When the divisor fits in 64 bits, we can use an optimized path |
| 52 | rem.* = n[high] % d[low]; | 110 | if (b[hi] == 0) { |
| 53 | } | 111 | r[hi] = 0; |
| 54 | return n[high] / d[low]; | 112 | if (a[hi] < b[lo]) { |
| 55 | } | 113 | // The result fits in 64 bits |
| 56 | // d[high] != 0 | 114 | q[hi] = 0; |
| 57 | if (n[low] == 0) { | 115 | q[lo] = divwide(HalfT, a[hi], a[lo], b[lo], &r[lo]); |
| 58 | // K 0 | | |
| 59 | // --- | | |
| 60 | // K 0 | | |
| 61 | if (maybe_rem) |rem| { | | |
| 62 | r[high] = n[high] % d[high]; | | |
| 63 | r[low] = 0; | | |
| 64 | rem.* = @bitCast(DoubleInt, r); | | |
| 65 | } | | |
| 66 | return n[high] / d[high]; | | |
| 67 | } | | |
| 68 | // K K | | |
| 69 | // --- | | |
| 70 | // K 0 | | |
| 71 | if ((d[high] & (d[high] - 1)) == 0) { | | |
| 72 | // d is a power of 2 | | |
| 73 | if (maybe_rem) |rem| { | | |
| 74 | r[low] = n[low]; | | |
| 75 | r[high] = n[high] & (d[high] - 1); | | |
| 76 | rem.* = @bitCast(DoubleInt, r); | | |
| 77 | } | | |
| 78 | return n[high] >> @intCast(Log2SingleInt, @ctz(d[high])); | | |
| 79 | } | | |
| 80 | // K K | | |
| 81 | // --- | | |
| 82 | // K 0 | | |
| 83 | sr = @bitCast(c_uint, @as(c_int, @clz(d[high])) - @as(c_int, @clz(n[high]))); | | |
| 84 | // 0 <= sr <= single_int_bits - 2 or sr large | | |
| 85 | if (sr > single_int_bits - 2) { | | |
| 86 | if (maybe_rem) |rem| { | | |
| 87 | rem.* = a; | | |
| 88 | } | | |
| 89 | return 0; | | |
| 90 | } | | |
| 91 | sr += 1; | | |
| 92 | // 1 <= sr <= single_int_bits - 1 | | |
| 93 | // q.all = a << (double_int_bits - sr); | | |
| 94 | q[low] = 0; | | |
| 95 | q[high] = n[low] << @intCast(Log2SingleInt, single_int_bits - sr); | | |
| 96 | // r.all = a >> sr; | | |
| 97 | r[high] = n[high] >> @intCast(Log2SingleInt, sr); | | |
| 98 | r[low] = (n[high] << @intCast(Log2SingleInt, single_int_bits - sr)) | (n[low] >> @intCast(Log2SingleInt, sr)); | | |
| 99 | } else { | | |
| 100 | // d[low] != 0 | | |
| 101 | if (d[high] == 0) { | | |
| 102 | // K X | | |
| 103 | // --- | | |
| 104 | // 0 K | | |
| 105 | if ((d[low] & (d[low] - 1)) == 0) { | | |
| 106 | // d is a power of 2 | | |
| 107 | if (maybe_rem) |rem| { | | |
| 108 | rem.* = n[low] & (d[low] - 1); | | |
| 109 | } | | |
| 110 | if (d[low] == 1) { | | |
| 111 | return a; | | |
| 112 | } | | |
| 113 | sr = @ctz(d[low]); | | |
| 114 | q[high] = n[high] >> @intCast(Log2SingleInt, sr); | | |
| 115 | q[low] = (n[high] << @intCast(Log2SingleInt, single_int_bits - sr)) | (n[low] >> @intCast(Log2SingleInt, sr)); | | |
| 116 | return @bitCast(DoubleInt, q); | | |
| 117 | } | | |
| 118 | // K X | | |
| 119 | // --- | | |
| 120 | // 0 K | | |
| 121 | sr = 1 + single_int_bits + @as(c_uint, @clz(d[low])) - @as(c_uint, @clz(n[high])); | | |
| 122 | // 2 <= sr <= double_int_bits - 1 | | |
| 123 | // q.all = a << (double_int_bits - sr); | | |
| 124 | // r.all = a >> sr; | | |
| 125 | if (sr == single_int_bits) { | | |
| 126 | q[low] = 0; | | |
| 127 | q[high] = n[low]; | | |
| 128 | r[high] = 0; | | |
| 129 | r[low] = n[high]; | | |
| 130 | } else if (sr < single_int_bits) { | | |
| 131 | // 2 <= sr <= single_int_bits - 1 | | |
| 132 | q[low] = 0; | | |
| 133 | q[high] = n[low] << @intCast(Log2SingleInt, single_int_bits - sr); | | |
| 134 | r[high] = n[high] >> @intCast(Log2SingleInt, sr); | | |
| 135 | r[low] = (n[high] << @intCast(Log2SingleInt, single_int_bits - sr)) | (n[low] >> @intCast(Log2SingleInt, sr)); | | |
| 136 | } else { | | |
| 137 | // single_int_bits + 1 <= sr <= double_int_bits - 1 | | |
| 138 | q[low] = n[low] << @intCast(Log2SingleInt, double_int_bits - sr); | | |
| 139 | q[high] = (n[high] << @intCast(Log2SingleInt, double_int_bits - sr)) | (n[low] >> @intCast(Log2SingleInt, sr - single_int_bits)); | | |
| 140 | r[high] = 0; | | |
| 141 | r[low] = n[high] >> @intCast(Log2SingleInt, sr - single_int_bits); | | |
| 142 | } | | |
| 143 | } else { | 116 | } else { |
| 144 | // K X | 117 | // First, divide with the high part to get the remainder. After that a_hi < b_lo. |
| 145 | // --- | 118 | q[hi] = a[hi] / b[lo]; |
| 146 | // K K | 119 | q[lo] = divwide(HalfT, a[hi] % b[lo], a[lo], b[lo], &r[lo]); |
| 147 | sr = @bitCast(c_uint, @as(c_int, @clz(d[high])) - @as(c_int, @clz(n[high]))); | | |
| 148 | // 0 <= sr <= single_int_bits - 1 or sr large | | |
| 149 | if (sr > single_int_bits - 1) { | | |
| 150 | if (maybe_rem) |rem| { | | |
| 151 | rem.* = a; | | |
| 152 | } | | |
| 153 | return 0; | | |
| 154 | } | | |
| 155 | sr += 1; | | |
| 156 | // 1 <= sr <= single_int_bits | | |
| 157 | // q.all = a << (double_int_bits - sr); | | |
| 158 | // r.all = a >> sr; | | |
| 159 | q[low] = 0; | | |
| 160 | if (sr == single_int_bits) { | | |
| 161 | q[high] = n[low]; | | |
| 162 | r[high] = 0; | | |
| 163 | r[low] = n[high]; | | |
| 164 | } else { | | |
| 165 | r[high] = n[high] >> @intCast(Log2SingleInt, sr); | | |
| 166 | r[low] = (n[high] << @intCast(Log2SingleInt, single_int_bits - sr)) | (n[low] >> @intCast(Log2SingleInt, sr)); | | |
| 167 | q[high] = n[low] << @intCast(Log2SingleInt, single_int_bits - sr); | | |
| 168 | } | | |
| 169 | } | 120 | } |
| | 121 | if (maybe_rem) |rem| { |
| | 122 | rem.* = @bitCast(T, r); |
| | 123 | } |
| | 124 | return @bitCast(T, q); |
| 170 | } | 125 | } |
| 171 | // Not a special case | 126 | |
| 172 | // q and r are initialized with: | 127 | // 0 <= shift <= 63 |
| 173 | // q.all = a << (double_int_bits - sr); | 128 | var shift: Log2Int(T) = @clz(b[hi]) - @clz(a[hi]); |
| 174 | // r.all = a >> sr; | 129 | var af = @bitCast(T, a); |
| 175 | // 1 <= sr <= double_int_bits - 1 | 130 | var bf = @bitCast(T, b) << shift; |
| 176 | var carry: u32 = 0; | 131 | q = @bitCast([2]HalfT, @as(T, 0)); |
| 177 | var r_all: DoubleInt = undefined; | 132 | |
| 178 | while (sr > 0) : (sr -= 1) { | 133 | for (0..shift + 1) |_| { |
| 179 | // r:q = ((r:q) << 1) | carry | 134 | q[lo] <<= 1; |
| 180 | r[high] = (r[high] << 1) | (r[low] >> (single_int_bits - 1)); | 135 | // Branchless version of: |
| 181 | r[low] = (r[low] << 1) | (q[high] >> (single_int_bits - 1)); | 136 | // if (af >= bf) { |
| 182 | q[high] = (q[high] << 1) | (q[low] >> (single_int_bits - 1)); | 137 | // af -= bf; |
| 183 | q[low] = (q[low] << 1) | carry; | 138 | // q[lo] |= 1; |
| 184 | // carry = 0; | | |
| 185 | // if (r.all >= b) | | |
| 186 | // { | | |
| 187 | // r.all -= b; | | |
| 188 | // carry = 1; | | |
| 189 | // } | 139 | // } |
| 190 | r_all = @bitCast(DoubleInt, r); | 140 | const s = @bitCast(SignedT, bf -% af -% 1) >> (@bitSizeOf(T) - 1); |
| 191 | const s: SignedDoubleInt = @bitCast(SignedDoubleInt, b -% r_all -% 1) >> (double_int_bits - 1); | 141 | q[lo] |= @intCast(HalfT, s & 1); |
| 192 | carry = @intCast(u32, s & 1); | 142 | af -= bf & @bitCast(T, s); |
| 193 | r_all -= b & @bitCast(DoubleInt, s); | 143 | bf >>= 1; |
| 194 | r = @bitCast([2]SingleInt, r_all); | | |
| 195 | } | 144 | } |
| 196 | const q_all = (@bitCast(DoubleInt, q) << 1) | carry; | | |
| 197 | if (maybe_rem) |rem| { | 145 | if (maybe_rem) |rem| { |
| 198 | rem.* = r_all; | 146 | rem.* = @bitCast(T, af); |
| 199 | } | 147 | } |
| 200 | return q_all; | 148 | return @bitCast(T, q); |
| 201 | } | 149 | } |