| ... | ... | @@ -18,11 +18,14 @@ const debug_safety = false; |
| 18 | 18 | /// Returns the number of limbs needed to store `scalar`, which must be a |
| 19 | 19 | /// primitive integer value. |
| 20 | 20 | pub fn calcLimbLen(scalar: anytype) usize { |
| 21 | | if (scalar == 0) { |
| 22 | | return 1; |
| 23 | | } |
| 21 | const T = @TypeOf(scalar); |
| 22 | const max_scalar = switch (@typeInfo(T)) { |
| 23 | .Int => maxInt(T), |
| 24 | .ComptimeInt => scalar, |
| 25 | else => @compileError("parameter must be a primitive integer type"), |
| 26 | }; |
| 24 | 27 | |
| 25 | | const w_value = std.math.absCast(scalar); |
| 28 | const w_value = std.math.absCast(max_scalar); |
| 26 | 29 | return @divFloor(@intCast(Limb, math.log2(w_value)), limb_bits) + 1; |
| 27 | 30 | } |
| 28 | 31 | |
| ... | ... | @@ -33,7 +36,7 @@ pub fn calcToStringLimbsBufferLen(a_len: usize, base: u8) usize { |
| 33 | 36 | } |
| 34 | 37 | |
| 35 | 38 | pub fn calcDivLimbsBufferLen(a_len: usize, b_len: usize) usize { |
| 36 | | return calcMulLimbsBufferLen(a_len, b_len, 2) * 4; |
| 39 | return a_len + b_len + 4; |
| 37 | 40 | } |
| 38 | 41 | |
| 39 | 42 | pub fn calcMulLimbsBufferLen(a_len: usize, b_len: usize, aliases: usize) usize { |
| ... | ... | @@ -760,8 +763,8 @@ pub const Mutable = struct { |
| 760 | 763 | /// q may alias with a or b. |
| 761 | 764 | /// |
| 762 | 765 | /// Asserts there is enough memory to store q and r. |
| 763 | | /// The upper bound for r limb count is a.limbs.len. |
| 764 | | /// The upper bound for q limb count is given by `a.limbs.len + b.limbs.len + 1`. |
| 766 | /// The upper bound for r limb count is `b.limbs.len`. |
| 767 | /// The upper bound for q limb count is given by `a.limbs`. |
| 765 | 768 | /// |
| 766 | 769 | /// If `allocator` is provided, it will be used for temporary storage to improve |
| 767 | 770 | /// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm. |
| ... | ... | @@ -773,20 +776,115 @@ pub const Mutable = struct { |
| 773 | 776 | a: Const, |
| 774 | 777 | b: Const, |
| 775 | 778 | limbs_buffer: []Limb, |
| 776 | | allocator: ?*Allocator, |
| 777 | 779 | ) void { |
| 778 | | div(q, r, a, b, limbs_buffer, allocator); |
| 780 | const sep = a.limbs.len + 2; |
| 781 | var x = a.toMutable(limbs_buffer[0..sep]); |
| 782 | var y = b.toMutable(limbs_buffer[sep..]); |
| 783 | |
| 784 | div(q, r, &x, &y); |
| 785 | |
| 786 | // Note, `div` performs truncating division, which satisfies |
| 787 | // @divTrunc(a, b) * b + @rem(a, b) = a |
| 788 | // so r = a - @divTrunc(a, b) * b |
| 789 | // Note, @rem(a, -b) = @rem(-b, a) = -@rem(a, b) = -@rem(-a, -b) |
| 790 | // For divTrunc, we want to perform |
| 791 | // @divFloor(a, b) * b + @mod(a, b) = a |
| 792 | // Note: |
| 793 | // @divFloor(-a, b) |
| 794 | // = @divFloor(a, -b) |
| 795 | // = -@divCeil(a, b) |
| 796 | // = -@divFloor(a + b - 1, b) |
| 797 | // = -@divTrunc(a + b - 1, b) |
| 798 | |
| 799 | // Note (1): |
| 800 | // @divTrunc(a + b - 1, b) * b + @rem(a + b - 1, b) = a + b - 1 |
| 801 | // = @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) = a + b - 1 |
| 802 | // = @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) - b + 1 = a |
| 803 | |
| 804 | if (a.positive and b.positive) { |
| 805 | // Positive-positive case, don't need to do anything. |
| 806 | } else if (a.positive and !b.positive) { |
| 807 | // a/-b -> q is negative, and so we need to fix flooring. |
| 808 | // Subtract one to make the division flooring. |
| 809 | |
| 810 | // @divFloor(a, -b) * -b + @mod(a, -b) = a |
| 811 | // If b divides a exactly, we have @divFloor(a, -b) * -b = a |
| 812 | // Else, we have @divFloor(a, -b) * -b > a, so @mod(a, -b) becomes negative |
| 813 | |
| 814 | // We have: |
| 815 | // @divFloor(a, -b) * -b + @mod(a, -b) = a |
| 816 | // = -@divTrunc(a + b - 1, b) * -b + @mod(a, -b) = a |
| 817 | // = @divTrunc(a + b - 1, b) * b + @mod(a, -b) = a |
| 818 | |
| 819 | // Substitute a for (1): |
| 820 | // @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) - b + 1 = @divTrunc(a + b - 1, b) * b + @mod(a, -b) |
| 821 | // Yields: |
| 822 | // @mod(a, -b) = @rem(a - 1, b) - b + 1 |
| 823 | // Note that `r` holds @rem(a, b) at this point. |
| 824 | // |
| 825 | // If @rem(a, b) is not 0: |
| 826 | // @rem(a - 1, b) = @rem(a, b) - 1 |
| 827 | // => @mod(a, -b) = @rem(a, b) - 1 - b + 1 = @rem(a, b) - b |
| 828 | // Else: |
| 829 | // @rem(a - 1, b) = @rem(a + b - 1, b) = @rem(b - 1, b) = b - 1 |
| 830 | // => @mod(a, -b) = b - 1 - b + 1 = 0 |
| 831 | if (!r.eqZero()) { |
| 832 | q.addScalar(q.toConst(), -1); |
| 833 | r.positive = true; |
| 834 | r.sub(r.toConst(), y.toConst().abs()); |
| 835 | } |
| 836 | } else if (!a.positive and b.positive) { |
| 837 | // -a/b -> q is negative, and so we need to fix flooring. |
| 838 | // Subtract one to make the division flooring. |
| 839 | |
| 840 | // @divFloor(-a, b) * b + @mod(-a, b) = a |
| 841 | // If b divides a exactly, we have @divFloor(-a, b) * b = -a |
| 842 | // Else, we have @divFloor(-a, b) * b < -a, so @mod(-a, b) becomes positive |
| 843 | |
| 844 | // We have: |
| 845 | // @divFloor(-a, b) * b + @mod(-a, b) = -a |
| 846 | // = -@divTrunc(a + b - 1, b) * b + @mod(-a, b) = -a |
| 847 | // = @divTrunc(a + b - 1, b) * b - @mod(-a, b) = a |
| 848 | |
| 849 | // Substitute a for (1): |
| 850 | // @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) - b + 1 = @divTrunc(a + b - 1, b) * b - @mod(-a, b) |
| 851 | // Yields: |
| 852 | // @rem(a - 1, b) - b + 1 = -@mod(-a, b) |
| 853 | // => -@mod(-a, b) = @rem(a - 1, b) - b + 1 |
| 854 | // => @mod(-a, b) = -(@rem(a - 1, b) - b + 1) = -@rem(a - 1, b) + b - 1 |
| 855 | // |
| 856 | // If @rem(a, b) is not 0: |
| 857 | // @rem(a - 1, b) = @rem(a, b) - 1 |
| 858 | // => @mod(-a, b) = -(@rem(a, b) - 1) + b - 1 = -@rem(a, b) + 1 + b - 1 = -@rem(a, b) + b |
| 859 | // Else : |
| 860 | // @rem(a - 1, b) = b - 1 |
| 861 | // => @mod(-a, b) = -(b - 1) + b - 1 = 0 |
| 862 | if (!r.eqZero()) { |
| 863 | q.addScalar(q.toConst(), -1); |
| 864 | r.positive = false; |
| 865 | r.add(r.toConst(), y.toConst().abs()); |
| 866 | } |
| 867 | } else if (!a.positive and !b.positive) { |
| 868 | // a/b -> q is positive, don't need to do anything to fix flooring. |
| 779 | 869 | |
| 780 | | // Trunc -> Floor. |
| 781 | | if (a.positive and b.positive) return; |
| 870 | // @divFloor(-a, -b) * -b + @mod(-a, -b) = -a |
| 871 | // If b divides a exactly, we have @divFloor(-a, -b) * -b = -a |
| 872 | // Else, we have @divFloor(-a, -b) * -b > -a, so @mod(-a, -b) becomes negative |
| 782 | 873 | |
| 783 | | if ((!q.positive or q.eqZero()) and !r.eqZero()) { |
| 784 | | const one: Const = .{ .limbs = &[_]Limb{1}, .positive = true }; |
| 785 | | q.sub(q.toConst(), one); |
| 786 | | } |
| 874 | // We have: |
| 875 | // @divFloor(-a, -b) * -b + @mod(-a, -b) = -a |
| 876 | // = @divTrunc(a, b) * -b + @mod(-a, -b) = -a |
| 877 | // = @divTrunc(a, b) * b - @mod(-a, -b) = a |
| 878 | |
| 879 | // We also have: |
| 880 | // @divTrunc(a, b) * b + @rem(a, b) = a |
| 787 | 881 | |
| 788 | | r.mulNoAlias(q.toConst(), b, allocator); |
| 789 | | r.sub(a, r.toConst()); |
| 882 | // Substitute a: |
| 883 | // @divTrunc(a, b) * b + @rem(a, b) = @divTrunc(a, b) * b - @mod(-a, -b) |
| 884 | // => @rem(a, b) = -@mod(-a, -b) |
| 885 | // => @mod(-a, -b) = -@rem(a, b) |
| 886 | r.positive = false; |
| 887 | } |
| 790 | 888 | } |
| 791 | 889 | |
| 792 | 890 | /// q = a / b (rem r) |
| ... | ... | @@ -795,9 +893,8 @@ pub const Mutable = struct { |
| 795 | 893 | /// q may alias with a or b. |
| 796 | 894 | /// |
| 797 | 895 | /// Asserts there is enough memory to store q and r. |
| 798 | | /// The upper bound for r limb count is a.limbs.len. |
| 799 | | /// The upper bound for q limb count is given by `calcQuotientLimbLen`. This accounts |
| 800 | | /// for temporary space used by the division algorithm. |
| 896 | /// The upper bound for r limb count is `b.limbs.len`. |
| 897 | /// The upper bound for q limb count is given by `a.limbs.len`. |
| 801 | 898 | /// |
| 802 | 899 | /// If `allocator` is provided, it will be used for temporary storage to improve |
| 803 | 900 | /// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm. |
| ... | ... | @@ -809,10 +906,12 @@ pub const Mutable = struct { |
| 809 | 906 | a: Const, |
| 810 | 907 | b: Const, |
| 811 | 908 | limbs_buffer: []Limb, |
| 812 | | allocator: ?*Allocator, |
| 813 | 909 | ) void { |
| 814 | | div(q, r, a, b, limbs_buffer, allocator); |
| 815 | | r.positive = a.positive; |
| 910 | const sep = a.limbs.len + 2; |
| 911 | var x = a.toMutable(limbs_buffer[0..sep]); |
| 912 | var y = b.toMutable(limbs_buffer[sep..]); |
| 913 | |
| 914 | div(q, r, &x, &y); |
| 816 | 915 | } |
| 817 | 916 | |
| 818 | 917 | /// r = a << shift, in other words, r = a * 2^shift |
| ... | ... | @@ -1176,181 +1275,214 @@ pub const Mutable = struct { |
| 1176 | 1275 | result.copy(x.toConst()); |
| 1177 | 1276 | } |
| 1178 | 1277 | |
| 1179 | | /// Truncates by default. |
| 1180 | | fn div(quo: *Mutable, rem: *Mutable, a: Const, b: Const, limbs_buffer: []Limb, allocator: ?*Allocator) void { |
| 1181 | | assert(!b.eqZero()); // division by zero |
| 1182 | | assert(quo != rem); // illegal aliasing |
| 1278 | // Truncates by default. |
| 1279 | fn div(q: *Mutable, r: *Mutable, x: *Mutable, y: *Mutable) void { |
| 1280 | assert(!y.eqZero()); // division by zero |
| 1281 | assert(q != r); // illegal aliasing |
| 1282 | |
| 1283 | const q_positive = (x.positive == y.positive); |
| 1284 | const r_positive = x.positive; |
| 1183 | 1285 | |
| 1184 | | if (a.orderAbs(b) == .lt) { |
| 1185 | | // quo may alias a so handle rem first |
| 1186 | | rem.copy(a); |
| 1187 | | rem.positive = a.positive == b.positive; |
| 1286 | if (x.toConst().orderAbs(y.toConst()) == .lt) { |
| 1287 | // q may alias x so handle r first. |
| 1288 | r.copy(x.toConst()); |
| 1289 | r.positive = r_positive; |
| 1188 | 1290 | |
| 1189 | | quo.positive = true; |
| 1190 | | quo.len = 1; |
| 1191 | | quo.limbs[0] = 0; |
| 1291 | q.set(0); |
| 1192 | 1292 | return; |
| 1193 | 1293 | } |
| 1194 | 1294 | |
| 1195 | 1295 | // Handle trailing zero-words of divisor/dividend. These are not handled in the following |
| 1196 | 1296 | // algorithms. |
| 1197 | | const a_zero_limb_count = blk: { |
| 1198 | | var i: usize = 0; |
| 1199 | | while (i < a.limbs.len) : (i += 1) { |
| 1200 | | if (a.limbs[i] != 0) break; |
| 1201 | | } |
| 1202 | | break :blk i; |
| 1203 | | }; |
| 1204 | | const b_zero_limb_count = blk: { |
| 1205 | | var i: usize = 0; |
| 1206 | | while (i < b.limbs.len) : (i += 1) { |
| 1207 | | if (b.limbs[i] != 0) break; |
| 1208 | | } |
| 1209 | | break :blk i; |
| 1210 | | }; |
| 1297 | // Note, there must be a non-zero limb for either. |
| 1298 | // const x_trailing = std.mem.indexOfScalar(Limb, x.limbs[0..x.len], 0).?; |
| 1299 | // const y_trailing = std.mem.indexOfScalar(Limb, y.limbs[0..y.len], 0).?; |
| 1211 | 1300 | |
| 1212 | | const ab_zero_limb_count = math.min(a_zero_limb_count, b_zero_limb_count); |
| 1301 | const x_trailing = for (x.limbs[0..x.len]) |xi, i| { |
| 1302 | if (xi != 0) break i; |
| 1303 | } else unreachable; |
| 1213 | 1304 | |
| 1214 | | if (b.limbs.len - ab_zero_limb_count == 1) { |
| 1215 | | lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[ab_zero_limb_count..a.limbs.len], b.limbs[b.limbs.len - 1]); |
| 1216 | | quo.normalize(a.limbs.len - ab_zero_limb_count); |
| 1217 | | quo.positive = (a.positive == b.positive); |
| 1305 | const y_trailing = for (y.limbs[0..y.len]) |yi, i| { |
| 1306 | if (yi != 0) break i; |
| 1307 | } else unreachable; |
| 1218 | 1308 | |
| 1219 | | rem.len = 1; |
| 1220 | | rem.positive = true; |
| 1309 | const xy_trailing = math.min(x_trailing, y_trailing); |
| 1310 | |
| 1311 | if (y.len - xy_trailing == 1) { |
| 1312 | lldiv1(q.limbs, &r.limbs[0], x.limbs[xy_trailing..x.len], y.limbs[y.len - 1]); |
| 1313 | q.normalize(x.len - xy_trailing); |
| 1314 | q.positive = q_positive; |
| 1315 | |
| 1316 | r.len = 1; |
| 1317 | r.positive = r_positive; |
| 1221 | 1318 | } else { |
| 1222 | | // x and y are modified during division |
| 1223 | | const sep_len = calcMulLimbsBufferLen(a.limbs.len, b.limbs.len, 2); |
| 1224 | | const x_limbs = limbs_buffer[0 * sep_len ..][0..sep_len]; |
| 1225 | | const y_limbs = limbs_buffer[1 * sep_len ..][0..sep_len]; |
| 1226 | | const t_limbs = limbs_buffer[2 * sep_len ..][0..sep_len]; |
| 1227 | | const mul_limbs_buf = limbs_buffer[3 * sep_len ..][0..sep_len]; |
| 1228 | | |
| 1229 | | var x: Mutable = .{ |
| 1230 | | .limbs = x_limbs, |
| 1319 | // Shrink x, y such that the trailing zero limbs shared between are removed. |
| 1320 | var x0 = Mutable{ |
| 1321 | .limbs = x.limbs[xy_trailing..], |
| 1322 | .len = x.len - xy_trailing, |
| 1231 | 1323 | .positive = true, |
| 1232 | | .len = a.limbs.len - ab_zero_limb_count, |
| 1233 | 1324 | }; |
| 1234 | | var y: Mutable = .{ |
| 1235 | | .limbs = y_limbs, |
| 1325 | |
| 1326 | var y0 = Mutable{ |
| 1327 | .limbs = y.limbs[xy_trailing..], |
| 1328 | .len = y.len - xy_trailing, |
| 1236 | 1329 | .positive = true, |
| 1237 | | .len = b.limbs.len - ab_zero_limb_count, |
| 1238 | 1330 | }; |
| 1239 | 1331 | |
| 1240 | | // Shrink x, y such that the trailing zero limbs shared between are removed. |
| 1241 | | mem.copy(Limb, x.limbs, a.limbs[ab_zero_limb_count..a.limbs.len]); |
| 1242 | | mem.copy(Limb, y.limbs, b.limbs[ab_zero_limb_count..b.limbs.len]); |
| 1332 | divmod(q, r, &x0, &y0); |
| 1333 | q.positive = q_positive; |
| 1243 | 1334 | |
| 1244 | | divN(quo, rem, &x, &y, t_limbs, mul_limbs_buf, allocator); |
| 1245 | | quo.positive = (a.positive == b.positive); |
| 1335 | r.positive = r_positive; |
| 1246 | 1336 | } |
| 1247 | 1337 | |
| 1248 | | if (ab_zero_limb_count != 0) { |
| 1249 | | rem.shiftLeft(rem.toConst(), ab_zero_limb_count * limb_bits); |
| 1338 | if (xy_trailing != 0) { |
| 1339 | // Manually shift here since we know its limb aligned. |
| 1340 | mem.copyBackwards(Limb, r.limbs[xy_trailing..], r.limbs[0..r.len]); |
| 1341 | mem.set(Limb, r.limbs[0..xy_trailing], 0); |
| 1342 | r.len += xy_trailing; |
| 1250 | 1343 | } |
| 1251 | 1344 | } |
| 1252 | 1345 | |
| 1253 | 1346 | /// Handbook of Applied Cryptography, 14.20 |
| 1254 | 1347 | /// |
| 1255 | 1348 | /// x = qy + r where 0 <= r < y |
| 1256 | | fn divN( |
| 1349 | /// y is modified but returned intact. |
| 1350 | fn divmod( |
| 1257 | 1351 | q: *Mutable, |
| 1258 | 1352 | r: *Mutable, |
| 1259 | 1353 | x: *Mutable, |
| 1260 | 1354 | y: *Mutable, |
| 1261 | | tmp_limbs: []Limb, |
| 1262 | | mul_limb_buf: []Limb, |
| 1263 | | allocator: ?*Allocator, |
| 1264 | 1355 | ) void { |
| 1265 | | assert(y.len >= 2); |
| 1266 | | assert(x.len >= y.len); |
| 1267 | | assert(q.limbs.len >= x.len + y.len - 1); |
| 1268 | | |
| 1269 | | // See 3.2 |
| 1270 | | var backup_tmp_limbs: [3]Limb = undefined; |
| 1271 | | const t_limbs = if (tmp_limbs.len < 3) &backup_tmp_limbs else tmp_limbs; |
| 1272 | | |
| 1273 | | var tmp: Mutable = .{ |
| 1274 | | .limbs = t_limbs, |
| 1275 | | .len = 1, |
| 1276 | | .positive = true, |
| 1277 | | }; |
| 1278 | | tmp.limbs[0] = 0; |
| 1356 | // 0. |
| 1357 | // Normalize so that y[t] > b/2 |
| 1358 | const lz = @clz(Limb, y.limbs[y.len - 1]); |
| 1359 | const norm_shift = if (lz == 0 and y.toConst().isOdd()) |
| 1360 | limb_bits // Force an extra limb so that y is even. |
| 1361 | else |
| 1362 | lz; |
| 1279 | 1363 | |
| 1280 | | // Normalize so y > limb_bits / 2 (i.e. leading bit is set) and even |
| 1281 | | var norm_shift = @clz(Limb, y.limbs[y.len - 1]); |
| 1282 | | if (norm_shift == 0 and y.toConst().isOdd()) { |
| 1283 | | norm_shift = limb_bits; |
| 1284 | | } |
| 1285 | 1364 | x.shiftLeft(x.toConst(), norm_shift); |
| 1286 | 1365 | y.shiftLeft(y.toConst(), norm_shift); |
| 1287 | 1366 | |
| 1288 | 1367 | const n = x.len - 1; |
| 1289 | 1368 | const t = y.len - 1; |
| 1369 | const shift = n - t; |
| 1290 | 1370 | |
| 1291 | 1371 | // 1. |
| 1292 | | q.len = n - t + 1; |
| 1372 | // for 0 <= j <= n - t, set q[j] to 0 |
| 1373 | q.len = shift + 1; |
| 1293 | 1374 | q.positive = true; |
| 1294 | 1375 | mem.set(Limb, q.limbs[0..q.len], 0); |
| 1295 | 1376 | |
| 1296 | 1377 | // 2. |
| 1297 | | tmp.shiftLeft(y.toConst(), limb_bits * (n - t)); |
| 1298 | | while (x.toConst().order(tmp.toConst()) != .lt) { |
| 1299 | | q.limbs[n - t] += 1; |
| 1300 | | x.sub(x.toConst(), tmp.toConst()); |
| 1378 | // while x >= y * b^(n - t): |
| 1379 | // x -= y * b^(n - t) |
| 1380 | // q[n - t] += 1 |
| 1381 | // Note, this algorithm is performed only once if y[t] > radix/2 and y is even, which we |
| 1382 | // enforced in step 0. This means we can replace the while with an if. |
| 1383 | // Note, multiplication by b^(n - t) comes down to shifting to the right by n - t limbs. |
| 1384 | // We can also replace x >= y * b^(n - t) by x/b^(n - t) >= y, and use shifts for that. |
| 1385 | { |
| 1386 | // x >= y * b^(n - t) can be replaced by x/b^(n - t) >= y. |
| 1387 | |
| 1388 | // 'divide' x by b^(n - t) |
| 1389 | var tmp = Mutable{ |
| 1390 | .limbs = x.limbs[shift..], |
| 1391 | .len = x.len - shift, |
| 1392 | .positive = true, |
| 1393 | }; |
| 1394 | |
| 1395 | if (tmp.toConst().order(y.toConst()) != .lt) { |
| 1396 | // Perform x -= y * b^(n - t) |
| 1397 | // Note, we can subtract y from x[n - t..] and get the result without shifting. |
| 1398 | // We can also re-use tmp which already contains the relevant part of x. Note that |
| 1399 | // this also edits x. |
| 1400 | // Due to the check above, this cannot underflow. |
| 1401 | tmp.sub(tmp.toConst(), y.toConst()); |
| 1402 | |
| 1403 | // tmp.sub normalized tmp, but we need to normalize x now. |
| 1404 | x.limbs.len = tmp.limbs.len + shift; |
| 1405 | |
| 1406 | q.limbs[shift] += 1; |
| 1407 | } |
| 1301 | 1408 | } |
| 1302 | 1409 | |
| 1303 | 1410 | // 3. |
| 1411 | // for i from n down to t + 1, do |
| 1304 | 1412 | var i = n; |
| 1305 | | while (i > t) : (i -= 1) { |
| 1306 | | // 3.1 |
| 1413 | while (i >= t + 1) : (i -= 1) { |
| 1414 | const k = i - t - 1; |
| 1415 | // 3.1. |
| 1416 | // if x_i == y_t: |
| 1417 | // q[i - t - 1] = b - 1 |
| 1418 | // else: |
| 1419 | // q[i - t - 1] = (x[i] * b + x[i - 1]) / y[t] |
| 1307 | 1420 | if (x.limbs[i] == y.limbs[t]) { |
| 1308 | | q.limbs[i - t - 1] = maxInt(Limb); |
| 1421 | q.limbs[k] = maxInt(Limb); |
| 1309 | 1422 | } else { |
| 1310 | | const num = (@as(DoubleLimb, x.limbs[i]) << limb_bits) | @as(DoubleLimb, x.limbs[i - 1]); |
| 1311 | | const z = @intCast(Limb, num / @as(DoubleLimb, y.limbs[t])); |
| 1312 | | q.limbs[i - t - 1] = if (z > maxInt(Limb)) maxInt(Limb) else @as(Limb, z); |
| 1423 | const q0 = (@as(DoubleLimb, x.limbs[i]) << limb_bits) | @as(DoubleLimb, x.limbs[i - 1]); |
| 1424 | const n0 = @as(DoubleLimb, y.limbs[t]); |
| 1425 | q.limbs[k] = @intCast(Limb, q0 / n0); |
| 1313 | 1426 | } |
| 1314 | 1427 | |
| 1315 | 1428 | // 3.2 |
| 1316 | | tmp.limbs[0] = if (i >= 2) x.limbs[i - 2] else 0; |
| 1317 | | tmp.limbs[1] = if (i >= 1) x.limbs[i - 1] else 0; |
| 1318 | | tmp.limbs[2] = x.limbs[i]; |
| 1319 | | tmp.normalize(3); |
| 1429 | // while q[i - t - 1] * (y[t] * b + y[t - 1] > x[i] * b * b + x[i - 1] + x[i - 2]: |
| 1430 | // q[i - t - 1] -= 1 |
| 1431 | // Note, if y[t] > b / 2 this part is repeated no more than twice. |
| 1432 | |
| 1433 | // Extract from y. |
| 1434 | const y0 = if (t > 0) y.limbs[t - 1] else 0; |
| 1435 | const y1 = y.limbs[t]; |
| 1436 | |
| 1437 | // Extract from x. |
| 1438 | // Note, big endian. |
| 1439 | const tmp0 = [_]Limb{ |
| 1440 | x.limbs[i], |
| 1441 | if (i >= 1) x.limbs[i - 1] else 0, |
| 1442 | if (i >= 2) x.limbs[i - 2] else 0, |
| 1443 | }; |
| 1320 | 1444 | |
| 1321 | 1445 | while (true) { |
| 1322 | | // 2x1 limb multiplication unrolled against single-limb q[i-t-1] |
| 1323 | | var carry: Limb = 0; |
| 1324 | | r.limbs[0] = addMulLimbWithCarry(0, if (t >= 1) y.limbs[t - 1] else 0, q.limbs[i - t - 1], &carry); |
| 1325 | | r.limbs[1] = addMulLimbWithCarry(0, y.limbs[t], q.limbs[i - t - 1], &carry); |
| 1326 | | r.limbs[2] = carry; |
| 1327 | | r.normalize(3); |
| 1328 | | |
| 1329 | | if (r.toConst().orderAbs(tmp.toConst()) != .gt) { |
| 1446 | // Ad-hoc 2x1 multiplication with q[i - t - 1]. |
| 1447 | // Note, big endian. |
| 1448 | var tmp1 = [_]Limb{ 0, undefined, undefined }; |
| 1449 | tmp1[2] = addMulLimbWithCarry(0, y0, q.limbs[k], &tmp1[0]); |
| 1450 | tmp1[1] = addMulLimbWithCarry(0, y1, q.limbs[k], &tmp1[0]); |
| 1451 | |
| 1452 | // Big-endian compare |
| 1453 | if (mem.order(Limb, &tmp1, &tmp0) != .gt) |
| 1330 | 1454 | break; |
| 1331 | | } |
| 1332 | 1455 | |
| 1333 | | q.limbs[i - t - 1] -= 1; |
| 1456 | q.limbs[k] -= 1; |
| 1334 | 1457 | } |
| 1335 | 1458 | |
| 1336 | | // 3.3 |
| 1337 | | tmp.set(q.limbs[i - t - 1]); |
| 1338 | | tmp.mul(tmp.toConst(), y.toConst(), mul_limb_buf, allocator); |
| 1339 | | tmp.shiftLeft(tmp.toConst(), limb_bits * (i - t - 1)); |
| 1340 | | x.sub(x.toConst(), tmp.toConst()); |
| 1341 | | |
| 1342 | | if (!x.positive) { |
| 1343 | | tmp.shiftLeft(y.toConst(), limb_bits * (i - t - 1)); |
| 1344 | | x.add(x.toConst(), tmp.toConst()); |
| 1345 | | q.limbs[i - t - 1] -= 1; |
| 1459 | // 3.3. |
| 1460 | // x -= q[i - t - 1] * y * b^(i - t - 1) |
| 1461 | // Note, we multiply by a single limb here. |
| 1462 | // The shift doesn't need to be performed if we add the result of the first multiplication |
| 1463 | // to x[i - t - 1]. |
| 1464 | // mem.set(Limb, x.limbs, 0); |
| 1465 | const underflow = llmulLimb(.sub, x.limbs[k..x.len], y.limbs[0..y.len], q.limbs[k]); |
| 1466 | |
| 1467 | // 3.4. |
| 1468 | // if x < 0: |
| 1469 | // x += y * b^(i - t - 1) |
| 1470 | // q[i - t - 1] -= 1 |
| 1471 | // Note, we check for x < 0 using the underflow flag from the previous operation. |
| 1472 | if (underflow) { |
| 1473 | // While we didn't properly set the signedness of x, this operation should 'flow' it back to positive. |
| 1474 | llaccum(.add, x.limbs[k..x.len], y.limbs[0..y.len]); |
| 1475 | q.limbs[k] -= 1; |
| 1346 | 1476 | } |
| 1477 | |
| 1478 | x.normalize(x.len); |
| 1347 | 1479 | } |
| 1348 | 1480 | |
| 1349 | | // Denormalize |
| 1350 | 1481 | q.normalize(q.len); |
| 1351 | 1482 | |
| 1483 | // De-normalize r and y. |
| 1352 | 1484 | r.shiftRight(x.toConst(), norm_shift); |
| 1353 | | r.normalize(r.len); |
| 1485 | y.shiftRight(y.toConst(), norm_shift); |
| 1354 | 1486 | } |
| 1355 | 1487 | |
| 1356 | 1488 | /// Truncate an integer to a number of bits, following 2s-complement semantics. |
| ... | ... | @@ -1808,7 +1940,7 @@ pub const Const = struct { |
| 1808 | 1940 | while (q.len >= 2) { |
| 1809 | 1941 | // Passing an allocator here would not be helpful since this division is destroying |
| 1810 | 1942 | // information, not creating it. [TODO citation needed] |
| 1811 | | q.divTrunc(&r, q.toConst(), b, rest_of_the_limbs_buf, null); |
| 1943 | q.divTrunc(&r, q.toConst(), b, rest_of_the_limbs_buf); |
| 1812 | 1944 | |
| 1813 | 1945 | var r_word = r.limbs[0]; |
| 1814 | 1946 | var i: usize = 0; |
| ... | ... | @@ -2435,16 +2567,14 @@ pub const Managed = struct { |
| 2435 | 2567 | /// a / b are floored (rounded towards 0). |
| 2436 | 2568 | /// |
| 2437 | 2569 | /// Returns an error if memory could not be allocated. |
| 2438 | | /// |
| 2439 | | /// q's allocator is used for temporary storage to speed up the multiplication. |
| 2440 | 2570 | pub fn divFloor(q: *Managed, r: *Managed, a: Const, b: Const) !void { |
| 2441 | | try q.ensureCapacity(a.limbs.len + b.limbs.len + 1); |
| 2442 | | try r.ensureCapacity(a.limbs.len); |
| 2571 | try q.ensureCapacity(a.limbs.len); |
| 2572 | try r.ensureCapacity(b.limbs.len); |
| 2443 | 2573 | var mq = q.toMutable(); |
| 2444 | 2574 | var mr = r.toMutable(); |
| 2445 | 2575 | const limbs_buffer = try q.allocator.alloc(Limb, calcDivLimbsBufferLen(a.limbs.len, b.limbs.len)); |
| 2446 | 2576 | defer q.allocator.free(limbs_buffer); |
| 2447 | | mq.divFloor(&mr, a, b, limbs_buffer, q.allocator); |
| 2577 | mq.divFloor(&mr, a, b, limbs_buffer); |
| 2448 | 2578 | q.setMetadata(mq.positive, mq.len); |
| 2449 | 2579 | r.setMetadata(mr.positive, mr.len); |
| 2450 | 2580 | } |
| ... | ... | @@ -2454,16 +2584,14 @@ pub const Managed = struct { |
| 2454 | 2584 | /// a / b are truncated (rounded towards -inf). |
| 2455 | 2585 | /// |
| 2456 | 2586 | /// Returns an error if memory could not be allocated. |
| 2457 | | /// |
| 2458 | | /// q's allocator is used for temporary storage to speed up the multiplication. |
| 2459 | 2587 | pub fn divTrunc(q: *Managed, r: *Managed, a: Const, b: Const) !void { |
| 2460 | | try q.ensureCapacity(a.limbs.len + b.limbs.len + 1); |
| 2461 | | try r.ensureCapacity(a.limbs.len); |
| 2588 | try q.ensureCapacity(a.limbs.len); |
| 2589 | try r.ensureCapacity(b.limbs.len); |
| 2462 | 2590 | var mq = q.toMutable(); |
| 2463 | 2591 | var mr = r.toMutable(); |
| 2464 | 2592 | const limbs_buffer = try q.allocator.alloc(Limb, calcDivLimbsBufferLen(a.limbs.len, b.limbs.len)); |
| 2465 | 2593 | defer q.allocator.free(limbs_buffer); |
| 2466 | | mq.divTrunc(&mr, a, b, limbs_buffer, q.allocator); |
| 2594 | mq.divTrunc(&mr, a, b, limbs_buffer); |
| 2467 | 2595 | q.setMetadata(mq.positive, mq.len); |
| 2468 | 2596 | r.setMetadata(mr.positive, mr.len); |
| 2469 | 2597 | } |
| ... | ... | @@ -2893,20 +3021,22 @@ fn llmulaccLong(comptime op: AccOp, r: []Limb, a: []const Limb, b: []const Limb) |
| 2893 | 3021 | |
| 2894 | 3022 | var i: usize = 0; |
| 2895 | 3023 | while (i < b.len) : (i += 1) { |
| 2896 | | llmulLimb(op, r[i..], a, b[i]); |
| 3024 | _ = llmulLimb(op, r[i..], a, b[i]); |
| 2897 | 3025 | } |
| 2898 | 3026 | } |
| 2899 | 3027 | |
| 2900 | 3028 | /// r = r (op) y * xi |
| 2901 | 3029 | /// The result is computed modulo `r.len`. |
| 2902 | | fn llmulLimb(comptime op: AccOp, acc: []Limb, y: []const Limb, xi: Limb) void { |
| 3030 | /// Returns whether the operation overflowed. |
| 3031 | fn llmulLimb(comptime op: AccOp, acc: []Limb, y: []const Limb, xi: Limb) bool { |
| 2903 | 3032 | @setRuntimeSafety(debug_safety); |
| 2904 | 3033 | if (xi == 0) { |
| 2905 | | return; |
| 3034 | return false; |
| 2906 | 3035 | } |
| 2907 | 3036 | |
| 2908 | | var a_lo = acc[0..y.len]; |
| 2909 | | var a_hi = acc[y.len..]; |
| 3037 | const split = std.math.min(y.len, acc.len); |
| 3038 | var a_lo = acc[0..split]; |
| 3039 | var a_hi = acc[split..]; |
| 2910 | 3040 | |
| 2911 | 3041 | switch (op) { |
| 2912 | 3042 | .add => { |
| ... | ... | @@ -2920,6 +3050,8 @@ fn llmulLimb(comptime op: AccOp, acc: []Limb, y: []const Limb, xi: Limb) void { |
| 2920 | 3050 | while ((carry != 0) and (j < a_hi.len)) : (j += 1) { |
| 2921 | 3051 | carry = @boolToInt(@addWithOverflow(Limb, a_hi[j], carry, &a_hi[j])); |
| 2922 | 3052 | } |
| 3053 | |
| 3054 | return carry != 0; |
| 2923 | 3055 | }, |
| 2924 | 3056 | .sub => { |
| 2925 | 3057 | var borrow: Limb = 0; |
| ... | ... | @@ -2932,6 +3064,8 @@ fn llmulLimb(comptime op: AccOp, acc: []Limb, y: []const Limb, xi: Limb) void { |
| 2932 | 3064 | while ((borrow != 0) and (j < a_hi.len)) : (j += 1) { |
| 2933 | 3065 | borrow = @boolToInt(@subWithOverflow(Limb, a_hi[j], borrow, &a_hi[j])); |
| 2934 | 3066 | } |
| 3067 | |
| 3068 | return borrow != 0; |
| 2935 | 3069 | }, |
| 2936 | 3070 | } |
| 2937 | 3071 | } |
| ... | ... | @@ -3424,7 +3558,8 @@ fn llsquareBasecase(r: []Limb, x: []const Limb) void { |
| 3424 | 3558 | |
| 3425 | 3559 | for (x_norm) |v, i| { |
| 3426 | 3560 | // Accumulate all the x[i]*x[j] (with x!=j) products |
| 3427 | | llmulLimb(.add, r[2 * i + 1 ..], x_norm[i + 1 ..], v); |
| 3561 | const overflow = llmulLimb(.add, r[2 * i + 1 ..], x_norm[i + 1 ..], v); |
| 3562 | assert(!overflow); |
| 3428 | 3563 | } |
| 3429 | 3564 | |
| 3430 | 3565 | // Each product appears twice, multiply by 2 |
| ... | ... | @@ -3432,7 +3567,8 @@ fn llsquareBasecase(r: []Limb, x: []const Limb) void { |
| 3432 | 3567 | |
| 3433 | 3568 | for (x_norm) |v, i| { |
| 3434 | 3569 | // Compute and add the squares |
| 3435 | | llmulLimb(.add, r[2 * i ..], x[i .. i + 1], v); |
| 3570 | const overflow = llmulLimb(.add, r[2 * i ..], x[i .. i + 1], v); |
| 3571 | assert(!overflow); |
| 3436 | 3572 | } |
| 3437 | 3573 | } |
| 3438 | 3574 | |
Andrew Kelley <andrew@ziglang.org> 2021-10-23 22:49:33-04:00
GitHub <noreply@github.com> 2021-10-23 22:49:33-04:00
Signed by PGP key 4AEE18F83AFDEB23
lib