| ... | ... | @@ -774,17 +774,114 @@ pub const Mutable = struct { |
| 774 | 774 | b: Const, |
| 775 | 775 | limbs_buffer: []Limb, |
| 776 | 776 | ) void { |
| 777 | | div(q, r, a, b, limbs_buffer); |
| 777 | const sep = a.limbs.len + 2; |
| 778 | var x = a.toMutable(limbs_buffer[0..sep]); |
| 779 | var y = b.toMutable(limbs_buffer[sep..]); |
| 780 | |
| 781 | div(q, r, &x, &y); |
| 782 | |
| 783 | // Note, `div` performs truncating division, which satisfies |
| 784 | // @divTrunc(a, b) * b + @rem(a, b) = a |
| 785 | // so r = a - @divTrunc(a, b) * b |
| 786 | // Note, @rem(a, -b) = @rem(-b, a) = -@rem(a, b) = -@rem(-a, -b) |
| 787 | // For divTrunc, we want to perform |
| 788 | // @divFloor(a, b) * b + @mod(a, b) = a |
| 789 | // Note: |
| 790 | // @divFloor(-a, b) |
| 791 | // = @divFloor(a, -b) |
| 792 | // = -@divCeil(a, b) |
| 793 | // = -@divFloor(a + b - 1, b) |
| 794 | // = -@divTrunc(a + b - 1, b) |
| 795 | |
| 796 | // Note (1): |
| 797 | // @divTrunc(a + b - 1, b) * b + @rem(a + b - 1, b) = a + b - 1 |
| 798 | // = @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) = a + b - 1 |
| 799 | // = @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) - b + 1 = a |
| 800 | |
| 801 | if (a.positive and b.positive) { |
| 802 | // Positive-positive case, don't need to do anything. |
| 803 | } else if (a.positive and !b.positive) { |
| 804 | // a/-b -> q is negative, and so we need to fix flooring. |
| 805 | // Subtract one to make the division flooring. |
| 806 | |
| 807 | // @divFloor(a, -b) * -b + @mod(a, -b) = a |
| 808 | // If b divides a exactly, we have @divFloor(a, -b) * -b = a |
| 809 | // Else, we have @divFloor(a, -b) * -b > a, so @mod(a, -b) becomes negative |
| 810 | |
| 811 | // We have: |
| 812 | // @divFloor(a, -b) * -b + @mod(a, -b) = a |
| 813 | // = -@divTrunc(a + b - 1, b) * -b + @mod(a, -b) = a |
| 814 | // = @divTrunc(a + b - 1, b) * b + @mod(a, -b) = a |
| 815 | |
| 816 | // Substitute a for (1): |
| 817 | // @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) - b + 1 = @divTrunc(a + b - 1, b) * b + @mod(a, -b) |
| 818 | // Yields: |
| 819 | // @mod(a, -b) = @rem(a - 1, b) - b + 1 |
| 820 | // Note that `r` holds @rem(a, b) at this point. |
| 821 | // |
| 822 | // If @rem(a, b) is not 0: |
| 823 | // @rem(a - 1, b) = @rem(a, b) - 1 |
| 824 | // => @mod(a, -b) = @rem(a, b) - 1 - b + 1 = @rem(a, b) - b |
| 825 | // Else: |
| 826 | // @rem(a - 1, b) = @rem(a + b - 1, b) = @rem(b - 1, b) = b - 1 |
| 827 | // => @mod(a, -b) = b - 1 - b + 1 = 0 |
| 828 | if (!r.eqZero()) { |
| 829 | q.addScalar(q.toConst(), -1); |
| 830 | r.positive = true; |
| 831 | r.sub(r.toConst(), y.toConst().abs()); |
| 832 | } |
| 833 | } else if (!a.positive and b.positive) { |
| 834 | // -a/b -> q is negative, and so we need to fix flooring. |
| 835 | // Subtract one to make the division flooring. |
| 836 | |
| 837 | // @divFloor(-a, b) * b + @mod(-a, b) = a |
| 838 | // If b divides a exactly, we have @divFloor(-a, b) * b = -a |
| 839 | // Else, we have @divFloor(-a, b) * b < -a, so @mod(-a, b) becomes positive |
| 840 | |
| 841 | // We have: |
| 842 | // @divFloor(-a, b) * b + @mod(-a, b) = -a |
| 843 | // = -@divTrunc(a + b - 1, b) * b + @mod(-a, b) = -a |
| 844 | // = @divTrunc(a + b - 1, b) * b - @mod(-a, b) = a |
| 845 | |
| 846 | // Substitute a for (1): |
| 847 | // @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) - b + 1 = @divTrunc(a + b - 1, b) * b - @mod(-a, b) |
| 848 | // Yields: |
| 849 | // @rem(a - 1, b) - b + 1 = -@mod(-a, b) |
| 850 | // => -@mod(-a, b) = @rem(a - 1, b) - b + 1 |
| 851 | // => @mod(-a, b) = -(@rem(a - 1, b) - b + 1) = -@rem(a - 1, b) + b - 1 |
| 852 | // |
| 853 | // If @rem(a, b) is not 0: |
| 854 | // @rem(a - 1, b) = @rem(a, b) - 1 |
| 855 | // => @mod(-a, b) = -(@rem(a, b) - 1) + b - 1 = -@rem(a, b) + 1 + b - 1 = -@rem(a, b) + b |
| 856 | // Else : |
| 857 | // @rem(a - 1, b) = b - 1 |
| 858 | // => @mod(-a, b) = -(b - 1) + b - 1 = 0 |
| 859 | if (!r.eqZero()) { |
| 860 | q.addScalar(q.toConst(), -1); |
| 861 | r.positive = false; |
| 862 | r.add(r.toConst(), y.toConst().abs()); |
| 863 | } |
| 864 | } else if (!a.positive and !b.positive) { |
| 865 | // a/b -> q is positive, don't need to do anything to fix flooring. |
| 778 | 866 | |
| 779 | | // Trunc -> Floor. |
| 780 | | if (a.positive and b.positive) return; |
| 867 | // @divFloor(-a, -b) * -b + @mod(-a, -b) = -a |
| 868 | // If b divides a exactly, we have @divFloor(-a, -b) * -b = -a |
| 869 | // Else, we have @divFloor(-a, -b) * -b > -a, so @mod(-a, -b) becomes negative |
| 781 | 870 | |
| 782 | | if ((!q.positive or q.eqZero()) and !r.eqZero()) { |
| 783 | | q.addScalar(q.toConst(), -1); |
| 784 | | } |
| 871 | // We have: |
| 872 | // @divFloor(-a, -b) * -b + @mod(-a, -b) = -a |
| 873 | // = @divTrunc(a, b) * -b + @mod(-a, -b) = -a |
| 874 | // = @divTrunc(a, b) * b - @mod(-a, -b) = a |
| 875 | |
| 876 | // We also have: |
| 877 | // @divTrunc(a, b) * b + @rem(a, b) = a |
| 785 | 878 | |
| 786 | | r.mulNoAlias(q.toConst(), b, null); |
| 787 | | r.sub(a, r.toConst()); |
| 879 | // Substitute a: |
| 880 | // @divTrunc(a, b) * b + @rem(a, b) = @divTrunc(a, b) * b - @mod(-a, -b) |
| 881 | // => @rem(a, b) = -@mod(-a, -b) |
| 882 | // => @mod(-a, -b) = -@rem(a, b) |
| 883 | r.positive = false; |
| 884 | } |
| 788 | 885 | } |
| 789 | 886 | |
| 790 | 887 | /// q = a / b (rem r) |
| ... | ... | @@ -808,8 +905,11 @@ pub const Mutable = struct { |
| 808 | 905 | b: Const, |
| 809 | 906 | limbs_buffer: []Limb, |
| 810 | 907 | ) void { |
| 811 | | div(q, r, a, b, limbs_buffer); |
| 812 | | r.positive = a.positive; |
| 908 | const sep = a.limbs.len + 2; |
| 909 | var x = a.toMutable(limbs_buffer[0..sep]); |
| 910 | var y = b.toMutable(limbs_buffer[sep..]); |
| 911 | |
| 912 | div(q, r, &x, &y); |
| 813 | 913 | } |
| 814 | 914 | |
| 815 | 915 | /// r = a << shift, in other words, r = a * 2^shift |
| ... | ... | @@ -1173,84 +1273,78 @@ pub const Mutable = struct { |
| 1173 | 1273 | result.copy(x.toConst()); |
| 1174 | 1274 | } |
| 1175 | 1275 | |
| 1176 | | /// Truncates by default. |
| 1177 | | fn div(quo: *Mutable, rem: *Mutable, a: Const, b: Const, limbs_buffer: []Limb) void { |
| 1178 | | assert(!b.eqZero()); // division by zero |
| 1179 | | assert(quo != rem); // illegal aliasing |
| 1276 | // Truncates by default. |
| 1277 | fn div(q: *Mutable, r: *Mutable, x: *Mutable, y: *Mutable) void { |
| 1278 | assert(!y.eqZero()); // division by zero |
| 1279 | assert(q != r); // illegal aliasing |
| 1180 | 1280 | |
| 1181 | | if (a.orderAbs(b) == .lt) { |
| 1182 | | // quo may alias a so handle rem first |
| 1183 | | rem.copy(a); |
| 1184 | | rem.positive = a.positive == b.positive; |
| 1281 | const q_positive = (x.positive == y.positive); |
| 1282 | const r_positive = x.positive; |
| 1185 | 1283 | |
| 1186 | | quo.positive = true; |
| 1187 | | quo.len = 1; |
| 1188 | | quo.limbs[0] = 0; |
| 1284 | if (x.toConst().orderAbs(y.toConst()) == .lt) { |
| 1285 | // q may alias x so handle r first. |
| 1286 | r.copy(x.toConst()); |
| 1287 | r.positive = r_positive; |
| 1288 | |
| 1289 | q.set(0); |
| 1189 | 1290 | return; |
| 1190 | 1291 | } |
| 1191 | 1292 | |
| 1192 | 1293 | // Handle trailing zero-words of divisor/dividend. These are not handled in the following |
| 1193 | 1294 | // algorithms. |
| 1194 | | const a_zero_limb_count = blk: { |
| 1195 | | var i: usize = 0; |
| 1196 | | while (i < a.limbs.len) : (i += 1) { |
| 1197 | | if (a.limbs[i] != 0) break; |
| 1198 | | } |
| 1199 | | break :blk i; |
| 1200 | | }; |
| 1201 | | const b_zero_limb_count = blk: { |
| 1202 | | var i: usize = 0; |
| 1203 | | while (i < b.limbs.len) : (i += 1) { |
| 1204 | | if (b.limbs[i] != 0) break; |
| 1205 | | } |
| 1206 | | break :blk i; |
| 1207 | | }; |
| 1295 | // Note, there must be a non-zero limb for either. |
| 1296 | // const x_trailing = std.mem.indexOfScalar(Limb, x.limbs[0..x.len], 0).?; |
| 1297 | // const y_trailing = std.mem.indexOfScalar(Limb, y.limbs[0..y.len], 0).?; |
| 1208 | 1298 | |
| 1209 | | const ab_zero_limb_count = math.min(a_zero_limb_count, b_zero_limb_count); |
| 1299 | const x_trailing = for (x.limbs[0..x.len]) |xi, i| { |
| 1300 | if (xi != 0) break i; |
| 1301 | } else unreachable; |
| 1210 | 1302 | |
| 1211 | | if (b.limbs.len - ab_zero_limb_count == 1) { |
| 1212 | | lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[ab_zero_limb_count..a.limbs.len], b.limbs[b.limbs.len - 1]); |
| 1213 | | quo.normalize(a.limbs.len - ab_zero_limb_count); |
| 1214 | | quo.positive = (a.positive == b.positive); |
| 1303 | const y_trailing = for (y.limbs[0..y.len]) |yi, i| { |
| 1304 | if (yi != 0) break i; |
| 1305 | } else unreachable; |
| 1215 | 1306 | |
| 1216 | | rem.len = 1; |
| 1217 | | rem.positive = true; |
| 1218 | | } else { |
| 1219 | | // x and y are modified during division |
| 1220 | | const sep_len = a.limbs.len + 2; |
| 1221 | | const x_limbs = limbs_buffer[0 .. sep_len]; |
| 1222 | | const y_limbs = limbs_buffer[sep_len..]; |
| 1307 | const xy_trailing = math.min(x_trailing, y_trailing); |
| 1308 | |
| 1309 | if (y.len - xy_trailing == 1) { |
| 1310 | lldiv1(q.limbs, &r.limbs[0], x.limbs[xy_trailing..x.len], y.limbs[y.len - 1]); |
| 1311 | q.normalize(x.len - xy_trailing); |
| 1312 | q.positive = q_positive; |
| 1223 | 1313 | |
| 1224 | | var x: Mutable = .{ |
| 1225 | | .limbs = x_limbs, |
| 1314 | r.len = 1; |
| 1315 | r.positive = r_positive; |
| 1316 | } else { |
| 1317 | // Shrink x, y such that the trailing zero limbs shared between are removed. |
| 1318 | var x0 = Mutable{ |
| 1319 | .limbs = x.limbs[xy_trailing..], |
| 1320 | .len = x.len - xy_trailing, |
| 1226 | 1321 | .positive = true, |
| 1227 | | .len = a.limbs.len - ab_zero_limb_count, |
| 1228 | 1322 | }; |
| 1229 | | var y: Mutable = .{ |
| 1230 | | .limbs = y_limbs, |
| 1323 | |
| 1324 | var y0 = Mutable{ |
| 1325 | .limbs = y.limbs[xy_trailing..], |
| 1326 | .len = y.len - xy_trailing, |
| 1231 | 1327 | .positive = true, |
| 1232 | | .len = b.limbs.len - ab_zero_limb_count, |
| 1233 | 1328 | }; |
| 1234 | 1329 | |
| 1235 | | // Shrink x, y such that the trailing zero limbs shared between are removed. |
| 1236 | | mem.copy(Limb, x.limbs, a.limbs[ab_zero_limb_count..]); |
| 1237 | | mem.copy(Limb, y.limbs, b.limbs[ab_zero_limb_count..]); |
| 1330 | divmod(q, r, &x0, &y0); |
| 1331 | q.positive = q_positive; |
| 1238 | 1332 | |
| 1239 | | divmod(quo, rem, &x, &y); |
| 1240 | | quo.positive = (a.positive == b.positive); |
| 1333 | r.positive = r_positive; |
| 1241 | 1334 | } |
| 1242 | 1335 | |
| 1243 | | if (ab_zero_limb_count != 0) { |
| 1336 | if (xy_trailing != 0) { |
| 1244 | 1337 | // Manually shift here since we know its limb aligned. |
| 1245 | | mem.copyBackwards(Limb, rem.limbs[ab_zero_limb_count..], rem.limbs[0..rem.len]); |
| 1246 | | mem.set(Limb, rem.limbs[0..ab_zero_limb_count], 0); |
| 1247 | | rem.len += ab_zero_limb_count; |
| 1338 | mem.copyBackwards(Limb, r.limbs[xy_trailing..], r.limbs[0..r.len]); |
| 1339 | mem.set(Limb, r.limbs[0..xy_trailing], 0); |
| 1340 | r.len += xy_trailing; |
| 1248 | 1341 | } |
| 1249 | 1342 | } |
| 1250 | 1343 | |
| 1251 | 1344 | /// Handbook of Applied Cryptography, 14.20 |
| 1252 | 1345 | /// |
| 1253 | 1346 | /// x = qy + r where 0 <= r < y |
| 1347 | /// y is modified but returned intact. |
| 1254 | 1348 | fn divmod( |
| 1255 | 1349 | q: *Mutable, |
| 1256 | 1350 | r: *Mutable, |
| ... | ... | @@ -1349,7 +1443,7 @@ pub const Mutable = struct { |
| 1349 | 1443 | while (true) { |
| 1350 | 1444 | // Ad-hoc 2x1 multiplication with q[i - t - 1]. |
| 1351 | 1445 | // Note, big endian. |
| 1352 | | var tmp1 = [_]Limb{0, undefined, undefined}; |
| 1446 | var tmp1 = [_]Limb{ 0, undefined, undefined }; |
| 1353 | 1447 | tmp1[2] = addMulLimbWithCarry(0, y0, q.limbs[k], &tmp1[0]); |
| 1354 | 1448 | tmp1[1] = addMulLimbWithCarry(0, y1, q.limbs[k], &tmp1[0]); |
| 1355 | 1449 | |
| ... | ... | @@ -1366,7 +1460,7 @@ pub const Mutable = struct { |
| 1366 | 1460 | // The shift doesn't need to be performed if we add the result of the first multiplication |
| 1367 | 1461 | // to x[i - t - 1]. |
| 1368 | 1462 | // mem.set(Limb, x.limbs, 0); |
| 1369 | | const underflow = llmulLimb(.sub, x.limbs[k .. x.len], y.limbs[0 .. y.len], q.limbs[k]); |
| 1463 | const underflow = llmulLimb(.sub, x.limbs[k..x.len], y.limbs[0..y.len], q.limbs[k]); |
| 1370 | 1464 | |
| 1371 | 1465 | // 3.4. |
| 1372 | 1466 | // if x < 0: |
| ... | ... | @@ -1375,7 +1469,7 @@ pub const Mutable = struct { |
| 1375 | 1469 | // Note, we check for x < 0 using the underflow flag from the previous operation. |
| 1376 | 1470 | if (underflow) { |
| 1377 | 1471 | // While we didn't properly set the signedness of x, this operation should 'flow' it back to positive. |
| 1378 | | llaccum(.add, x.limbs[k .. x.len], y.limbs[0 .. y.len]); |
| 1472 | llaccum(.add, x.limbs[k..x.len], y.limbs[0..y.len]); |
| 1379 | 1473 | q.limbs[k] -= 1; |
| 1380 | 1474 | } |
| 1381 | 1475 | |
| ... | ... | @@ -1384,8 +1478,9 @@ pub const Mutable = struct { |
| 1384 | 1478 | |
| 1385 | 1479 | q.normalize(q.len); |
| 1386 | 1480 | |
| 1387 | | // De-normalize r. |
| 1481 | // De-normalize r and y. |
| 1388 | 1482 | r.shiftRight(x.toConst(), norm_shift); |
| 1483 | y.shiftRight(y.toConst(), norm_shift); |
| 1389 | 1484 | } |
| 1390 | 1485 | |
| 1391 | 1486 | /// Truncate an integer to a number of bits, following 2s-complement semantics. |
Robin Voetter <robin@voetter.nl> 2021-10-24 02:39:56+02:00