| ... | ... | @@ -44,6 +44,11 @@ pub fn calcMulLimbsBufferLen(a_len: usize, b_len: usize, aliases: usize) usize { |
| 44 | 44 | return aliases * math.max(a_len, b_len); |
| 45 | 45 | } |
| 46 | 46 | |
| 47 | pub fn calcMulWrapLimbsBufferLen(bit_count: usize, a_len: usize, b_len: usize, aliases: usize) usize { |
| 48 | const req_limbs = calcTwosCompLimbCount(bit_count); |
| 49 | return aliases * math.min(req_limbs, math.max(a_len, b_len)); |
| 50 | } |
| 51 | |
| 47 | 52 | pub fn calcSetStringLimbsBufferLen(base: u8, string_len: usize) usize { |
| 48 | 53 | const limb_count = calcSetStringLimbCount(base, string_len); |
| 49 | 54 | return calcMulLimbsBufferLen(limb_count, limb_count, 2); |
| ... | ... | @@ -58,6 +63,11 @@ pub fn calcPowLimbsBufferLen(a_bit_count: usize, y: usize) usize { |
| 58 | 63 | return 2 + (a_bit_count * y + (limb_bits - 1)) / limb_bits; |
| 59 | 64 | } |
| 60 | 65 | |
| 66 | // Compute the number of limbs required to store a 2s-complement number of `bit_count` bits. |
| 67 | pub fn calcTwosCompLimbCount(bit_count: usize) usize { |
| 68 | return std.math.divCeil(usize, bit_count, @bitSizeOf(Limb)) catch unreachable; |
| 69 | } |
| 70 | |
| 61 | 71 | /// a + b * c + *carry, sets carry to the overflow bits |
| 62 | 72 | pub fn addMulLimbWithCarry(a: Limb, b: Limb, c: Limb, carry: *Limb) Limb { |
| 63 | 73 | @setRuntimeSafety(debug_safety); |
| ... | ... | @@ -81,6 +91,33 @@ pub fn addMulLimbWithCarry(a: Limb, b: Limb, c: Limb, carry: *Limb) Limb { |
| 81 | 91 | return r1; |
| 82 | 92 | } |
| 83 | 93 | |
| 94 | /// a - b * c - *carry, sets carry to the overflow bits |
| 95 | fn subMulLimbWithBorrow(a: Limb, b: Limb, c: Limb, carry: *Limb) Limb { |
| 96 | // r1 = a - *carry |
| 97 | var r1: Limb = undefined; |
| 98 | const c1: Limb = @boolToInt(@subWithOverflow(Limb, a, carry.*, &r1)); |
| 99 | |
| 100 | // r2 = b * c |
| 101 | const bc = @as(DoubleLimb, std.math.mulWide(Limb, b, c)); |
| 102 | const r2 = @truncate(Limb, bc); |
| 103 | const c2 = @truncate(Limb, bc >> limb_bits); |
| 104 | |
| 105 | // r1 = r1 - r2 |
| 106 | const c3: Limb = @boolToInt(@subWithOverflow(Limb, r1, r2, &r1)); |
| 107 | carry.* = c1 + c2 + c3; |
| 108 | |
| 109 | return r1; |
| 110 | } |
| 111 | |
| 112 | /// Used to indicate either limit of a 2s-complement integer. |
| 113 | pub const TwosCompIntLimit = enum { |
| 114 | // The low limit, either 0x00 (unsigned) or (-)0x80 (signed) for an 8-bit integer. |
| 115 | min, |
| 116 | |
| 117 | // The high limit, either 0xFF (unsigned) or 0x7F (signed) for an 8-bit integer. |
| 118 | max, |
| 119 | }; |
| 120 | |
| 84 | 121 | /// A arbitrary-precision big integer, with a fixed set of mutable limbs. |
| 85 | 122 | pub const Mutable = struct { |
| 86 | 123 | /// Raw digits. These are: |
| ... | ... | @@ -282,6 +319,75 @@ pub const Mutable = struct { |
| 282 | 319 | self.positive = positive; |
| 283 | 320 | } |
| 284 | 321 | |
| 322 | /// Set self to either bound of a 2s-complement integer. |
| 323 | /// Note: The result is still sign-magnitude, not twos complement! In order to convert the |
| 324 | /// result to twos complement, it is sufficient to take the absolute value. |
| 325 | /// |
| 326 | /// Asserts the result fits in `r`. An upper bound on the number of limbs needed by |
| 327 | /// r is `calcTwosCompLimbCount(bit_count)`. |
| 328 | pub fn setTwosCompIntLimit( |
| 329 | r: *Mutable, |
| 330 | limit: TwosCompIntLimit, |
| 331 | signedness: std.builtin.Signedness, |
| 332 | bit_count: usize, |
| 333 | ) void { |
| 334 | // Handle zero-bit types. |
| 335 | if (bit_count == 0) { |
| 336 | r.set(0); |
| 337 | return; |
| 338 | } |
| 339 | |
| 340 | const req_limbs = calcTwosCompLimbCount(bit_count); |
| 341 | const bit = @truncate(Log2Limb, bit_count - 1); |
| 342 | const signmask = @as(Limb, 1) << bit; // 0b0..010..0 where 1 is the sign bit. |
| 343 | const mask = (signmask << 1) -% 1; // 0b0..011..1 where the leftmost 1 is the sign bit. |
| 344 | |
| 345 | r.positive = true; |
| 346 | |
| 347 | switch (signedness) { |
| 348 | .signed => switch (limit) { |
| 349 | .min => { |
| 350 | // Negative bound, signed = -0x80. |
| 351 | r.len = req_limbs; |
| 352 | mem.set(Limb, r.limbs[0 .. r.len - 1], 0); |
| 353 | r.limbs[r.len - 1] = signmask; |
| 354 | r.positive = false; |
| 355 | }, |
| 356 | .max => { |
| 357 | // Positive bound, signed = 0x7F |
| 358 | // Note, in this branch we need to normalize because the first bit is |
| 359 | // supposed to be 0. |
| 360 | |
| 361 | // Special case for 1-bit integers. |
| 362 | if (bit_count == 1) { |
| 363 | r.set(0); |
| 364 | } else { |
| 365 | const new_req_limbs = calcTwosCompLimbCount(bit_count - 1); |
| 366 | const msb = @truncate(Log2Limb, bit_count - 2); |
| 367 | const new_signmask = @as(Limb, 1) << msb; // 0b0..010..0 where 1 is the sign bit. |
| 368 | const new_mask = (new_signmask << 1) -% 1; // 0b0..001..1 where the rightmost 0 is the sign bit. |
| 369 | |
| 370 | r.len = new_req_limbs; |
| 371 | std.mem.set(Limb, r.limbs[0 .. r.len - 1], maxInt(Limb)); |
| 372 | r.limbs[r.len - 1] = new_mask; |
| 373 | } |
| 374 | }, |
| 375 | }, |
| 376 | .unsigned => switch (limit) { |
| 377 | .min => { |
| 378 | // Min bound, unsigned = 0x00 |
| 379 | r.set(0); |
| 380 | }, |
| 381 | .max => { |
| 382 | // Max bound, unsigned = 0xFF |
| 383 | r.len = req_limbs; |
| 384 | std.mem.set(Limb, r.limbs[0 .. r.len - 1], maxInt(Limb)); |
| 385 | r.limbs[r.len - 1] = mask; |
| 386 | }, |
| 387 | }, |
| 388 | } |
| 389 | } |
| 390 | |
| 285 | 391 | /// r = a + scalar |
| 286 | 392 | /// |
| 287 | 393 | /// r and a may be aliases. |
| ... | ... | @@ -295,102 +401,220 @@ pub const Mutable = struct { |
| 295 | 401 | return add(r, a, operand); |
| 296 | 402 | } |
| 297 | 403 | |
| 298 | | /// r = a + b |
| 299 | | /// |
| 404 | /// Base implementation for addition. Adds `max(a.limbs.len, b.limbs.len)` elements from a and b, |
| 405 | /// and returns whether any overflow occured. |
| 300 | 406 | /// r, a and b may be aliases. |
| 301 | 407 | /// |
| 302 | | /// Asserts the result fits in `r`. An upper bound on the number of limbs needed by |
| 303 | | /// r is `math.max(a.limbs.len, b.limbs.len) + 1`. |
| 304 | | pub fn add(r: *Mutable, a: Const, b: Const) void { |
| 408 | /// Asserts r has enough elements to hold the result. The upper bound is `max(a.limbs.len, b.limbs.len)`. |
| 409 | fn addCarry(r: *Mutable, a: Const, b: Const) bool { |
| 305 | 410 | if (a.eqZero()) { |
| 306 | 411 | r.copy(b); |
| 307 | | return; |
| 412 | return false; |
| 308 | 413 | } else if (b.eqZero()) { |
| 309 | 414 | r.copy(a); |
| 310 | | return; |
| 311 | | } |
| 312 | | |
| 313 | | if (a.limbs.len == 1 and b.limbs.len == 1 and a.positive == b.positive) { |
| 314 | | var o: Limb = undefined; |
| 315 | | if (!@addWithOverflow(Limb, a.limbs[0], b.limbs[0], &o)) { |
| 316 | | r.limbs[0] = o; |
| 317 | | r.len = 1; |
| 318 | | r.positive = a.positive; |
| 319 | | return; |
| 320 | | } |
| 321 | | } |
| 322 | | |
| 323 | | if (a.positive != b.positive) { |
| 415 | return false; |
| 416 | } else if (a.positive != b.positive) { |
| 324 | 417 | if (a.positive) { |
| 325 | 418 | // (a) + (-b) => a - b |
| 326 | | r.sub(a, b.abs()); |
| 419 | return r.subCarry(a, b.abs()); |
| 327 | 420 | } else { |
| 328 | 421 | // (-a) + (b) => b - a |
| 329 | | r.sub(b, a.abs()); |
| 422 | return r.subCarry(b, a.abs()); |
| 330 | 423 | } |
| 331 | 424 | } else { |
| 425 | r.positive = a.positive; |
| 332 | 426 | if (a.limbs.len >= b.limbs.len) { |
| 333 | | lladd(r.limbs[0..], a.limbs, b.limbs); |
| 334 | | r.normalize(a.limbs.len + 1); |
| 427 | const c = lladdcarry(r.limbs, a.limbs, b.limbs); |
| 428 | r.normalize(a.limbs.len); |
| 429 | return c != 0; |
| 335 | 430 | } else { |
| 336 | | lladd(r.limbs[0..], b.limbs, a.limbs); |
| 337 | | r.normalize(b.limbs.len + 1); |
| 431 | const c = lladdcarry(r.limbs, b.limbs, a.limbs); |
| 432 | r.normalize(b.limbs.len); |
| 433 | return c != 0; |
| 338 | 434 | } |
| 435 | } |
| 436 | } |
| 339 | 437 | |
| 340 | | r.positive = a.positive; |
| 438 | /// r = a + b |
| 439 | /// r, a and b may be aliases. |
| 440 | /// |
| 441 | /// Asserts the result fits in `r`. An upper bound on the number of limbs needed by |
| 442 | /// r is `math.max(a.limbs.len, b.limbs.len) + 1`. |
| 443 | pub fn add(r: *Mutable, a: Const, b: Const) void { |
| 444 | if (r.addCarry(a, b)) { |
| 445 | // Fix up the result. Note that addCarry normalizes by a.limbs.len or b.limbs.len, |
| 446 | // so we need to set the length here. |
| 447 | const msl = math.max(a.limbs.len, b.limbs.len); |
| 448 | // `[add|sub]Carry` normalizes by `msl`, so we need to fix up the result manually here. |
| 449 | // Note, the fact that it normalized means that the intermediary limbs are zero here. |
| 450 | r.len = msl + 1; |
| 451 | r.limbs[msl] = 1; // If this panics, there wasn't enough space in `r`. |
| 341 | 452 | } |
| 342 | 453 | } |
| 343 | 454 | |
| 344 | | /// r = a - b |
| 455 | /// r = a + b with 2s-complement wrapping semantics. |
| 456 | /// r, a and b may be aliases |
| 345 | 457 | /// |
| 458 | /// Asserts the result fits in `r`. An upper bound on the number of limbs needed by |
| 459 | /// r is `calcTwosCompLimbCount(bit_count)`. |
| 460 | pub fn addWrap(r: *Mutable, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) void { |
| 461 | const req_limbs = calcTwosCompLimbCount(bit_count); |
| 462 | |
| 463 | // Slice of the upper bits if they exist, these will be ignored and allows us to use addCarry to determine |
| 464 | // if an overflow occured. |
| 465 | const x = Const{ |
| 466 | .positive = a.positive, |
| 467 | .limbs = a.limbs[0..math.min(req_limbs, a.limbs.len)], |
| 468 | }; |
| 469 | |
| 470 | const y = Const{ |
| 471 | .positive = b.positive, |
| 472 | .limbs = b.limbs[0..math.min(req_limbs, b.limbs.len)], |
| 473 | }; |
| 474 | |
| 475 | if (r.addCarry(x, y)) { |
| 476 | // There are two possibilities here: |
| 477 | // - We overflowed req_limbs. In this case, the carry is ignored, as it would be removed by |
| 478 | // truncate anyway. |
| 479 | // - a and b had less elements than req_limbs, and those were overflowed. This case needs to be handled. |
| 480 | // Note: after this we still might need to wrap. |
| 481 | const msl = math.max(a.limbs.len, b.limbs.len); |
| 482 | if (msl < req_limbs) { |
| 483 | r.limbs[msl] = 1; |
| 484 | r.len = req_limbs; |
| 485 | } |
| 486 | } |
| 487 | |
| 488 | r.truncate(r.toConst(), signedness, bit_count); |
| 489 | } |
| 490 | |
| 491 | /// r = a + b with 2s-complement saturating semantics. |
| 346 | 492 | /// r, a and b may be aliases. |
| 347 | 493 | /// |
| 348 | | /// Asserts the result fits in `r`. An upper bound on the number of limbs needed by |
| 349 | | /// r is `math.max(a.limbs.len, b.limbs.len) + 1`. The +1 is not needed if both operands are positive. |
| 350 | | pub fn sub(r: *Mutable, a: Const, b: Const) void { |
| 351 | | if (a.positive != b.positive) { |
| 494 | /// Assets the result fits in `r`. Upper bound on the number of limbs needed by |
| 495 | /// r is `calcTwosCompLimbCount(bit_count)`. |
| 496 | pub fn addSat(r: *Mutable, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) void { |
| 497 | const req_limbs = calcTwosCompLimbCount(bit_count); |
| 498 | |
| 499 | // Slice of the upper bits if they exist, these will be ignored and allows us to use addCarry to determine |
| 500 | // if an overflow occured. |
| 501 | const x = Const{ |
| 502 | .positive = a.positive, |
| 503 | .limbs = a.limbs[0..math.min(req_limbs, a.limbs.len)], |
| 504 | }; |
| 505 | |
| 506 | const y = Const{ |
| 507 | .positive = b.positive, |
| 508 | .limbs = b.limbs[0..math.min(req_limbs, b.limbs.len)], |
| 509 | }; |
| 510 | |
| 511 | if (r.addCarry(x, y)) { |
| 512 | // There are two possibilities here: |
| 513 | // - We overflowed req_limbs, in which case we need to saturate. |
| 514 | // - a and b had less elements than req_limbs, and those were overflowed. |
| 515 | // Note: In this case, might _also_ need to saturate. |
| 516 | const msl = math.max(a.limbs.len, b.limbs.len); |
| 517 | if (msl < req_limbs) { |
| 518 | r.limbs[msl] = 1; |
| 519 | r.len = req_limbs; |
| 520 | // Note: Saturation may still be required if msl == req_limbs - 1 |
| 521 | } else { |
| 522 | // Overflowed req_limbs, definitely saturate. |
| 523 | r.setTwosCompIntLimit(if (r.positive) .max else .min, signedness, bit_count); |
| 524 | } |
| 525 | } |
| 526 | |
| 527 | // Saturate if the result didn't fit. |
| 528 | r.saturate(r.toConst(), signedness, bit_count); |
| 529 | } |
| 530 | |
| 531 | /// Base implementation for subtraction. Subtracts `max(a.limbs.len, b.limbs.len)` elements from a and b, |
| 532 | /// and returns whether any overflow occured. |
| 533 | /// r, a and b may be aliases. |
| 534 | /// |
| 535 | /// Asserts r has enough elements to hold the result. The upper bound is `max(a.limbs.len, b.limbs.len)`. |
| 536 | fn subCarry(r: *Mutable, a: Const, b: Const) bool { |
| 537 | if (a.eqZero()) { |
| 538 | r.copy(b); |
| 539 | r.positive = !b.positive; |
| 540 | return false; |
| 541 | } else if (b.eqZero()) { |
| 542 | r.copy(a); |
| 543 | return false; |
| 544 | } else if (a.positive != b.positive) { |
| 352 | 545 | if (a.positive) { |
| 353 | 546 | // (a) - (-b) => a + b |
| 354 | | r.add(a, b.abs()); |
| 547 | return r.addCarry(a, b.abs()); |
| 355 | 548 | } else { |
| 356 | | // (-a) - (b) => -(a + b) |
| 357 | | r.add(a.abs(), b); |
| 358 | | r.positive = false; |
| 549 | // (-a) - (b) => -a + -b |
| 550 | return r.addCarry(a, b.negate()); |
| 359 | 551 | } |
| 360 | | } else { |
| 361 | | if (a.positive) { |
| 552 | } else if (a.positive) { |
| 553 | if (a.order(b) != .lt) { |
| 362 | 554 | // (a) - (b) => a - b |
| 363 | | if (a.order(b) != .lt) { |
| 364 | | llsub(r.limbs[0..], a.limbs[0..a.limbs.len], b.limbs[0..b.limbs.len]); |
| 365 | | r.normalize(a.limbs.len); |
| 366 | | r.positive = true; |
| 367 | | } else { |
| 368 | | llsub(r.limbs[0..], b.limbs[0..b.limbs.len], a.limbs[0..a.limbs.len]); |
| 369 | | r.normalize(b.limbs.len); |
| 370 | | r.positive = false; |
| 371 | | } |
| 555 | const c = llsubcarry(r.limbs, a.limbs, b.limbs); |
| 556 | r.normalize(a.limbs.len); |
| 557 | r.positive = true; |
| 558 | return c != 0; |
| 372 | 559 | } else { |
| 560 | // (a) - (b) => -b + a => -(b - a) |
| 561 | const c = llsubcarry(r.limbs, b.limbs, a.limbs); |
| 562 | r.normalize(b.limbs.len); |
| 563 | r.positive = false; |
| 564 | return c != 0; |
| 565 | } |
| 566 | } else { |
| 567 | if (a.order(b) == .lt) { |
| 373 | 568 | // (-a) - (-b) => -(a - b) |
| 374 | | if (a.order(b) == .lt) { |
| 375 | | llsub(r.limbs[0..], a.limbs[0..a.limbs.len], b.limbs[0..b.limbs.len]); |
| 376 | | r.normalize(a.limbs.len); |
| 377 | | r.positive = false; |
| 378 | | } else { |
| 379 | | llsub(r.limbs[0..], b.limbs[0..b.limbs.len], a.limbs[0..a.limbs.len]); |
| 380 | | r.normalize(b.limbs.len); |
| 381 | | r.positive = true; |
| 382 | | } |
| 569 | const c = llsubcarry(r.limbs, a.limbs, b.limbs); |
| 570 | r.normalize(a.limbs.len); |
| 571 | r.positive = false; |
| 572 | return c != 0; |
| 573 | } else { |
| 574 | // (-a) - (-b) => --b + -a => b - a |
| 575 | const c = llsubcarry(r.limbs, b.limbs, a.limbs); |
| 576 | r.normalize(b.limbs.len); |
| 577 | r.positive = true; |
| 578 | return c != 0; |
| 383 | 579 | } |
| 384 | 580 | } |
| 385 | 581 | } |
| 386 | 582 | |
| 583 | /// r = a - b |
| 584 | /// |
| 585 | /// r, a and b may be aliases. |
| 586 | /// |
| 587 | /// Asserts the result fits in `r`. An upper bound on the number of limbs needed by |
| 588 | /// r is `math.max(a.limbs.len, b.limbs.len) + 1`. The +1 is not needed if both operands are positive. |
| 589 | pub fn sub(r: *Mutable, a: Const, b: Const) void { |
| 590 | r.add(a, b.negate()); |
| 591 | } |
| 592 | |
| 593 | /// r = a - b with 2s-complement wrapping semantics. |
| 594 | /// |
| 595 | /// r, a and b may be aliases |
| 596 | /// Asserts the result fits in `r`. An upper bound on the number of limbs needed by |
| 597 | /// r is `calcTwosCompLimbCount(bit_count)`. |
| 598 | pub fn subWrap(r: *Mutable, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) void { |
| 599 | r.addWrap(a, b.negate(), signedness, bit_count); |
| 600 | } |
| 601 | |
| 602 | /// r = a - b with 2s-complement saturating semantics. |
| 603 | /// r, a and b may be aliases. |
| 604 | /// |
| 605 | /// Assets the result fits in `r`. Upper bound on the number of limbs needed by |
| 606 | /// r is `calcTwosCompLimbCount(bit_count)`. |
| 607 | pub fn subSat(r: *Mutable, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) void { |
| 608 | r.addSat(a, b.negate(), signedness, bit_count); |
| 609 | } |
| 610 | |
| 387 | 611 | /// rma = a * b |
| 388 | 612 | /// |
| 389 | 613 | /// `rma` may alias with `a` or `b`. |
| 390 | 614 | /// `a` and `b` may alias with each other. |
| 391 | 615 | /// |
| 392 | 616 | /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by |
| 393 | | /// rma is given by `a.limbs.len + b.limbs.len + 1`. |
| 617 | /// rma is given by `a.limbs.len + b.limbs.len`. |
| 394 | 618 | /// |
| 395 | 619 | /// `limbs_buffer` is used for temporary storage. The amount required is given by `calcMulLimbsBufferLen`. |
| 396 | 620 | pub fn mul(rma: *Mutable, a: Const, b: Const, limbs_buffer: []Limb, allocator: ?*Allocator) void { |
| ... | ... | @@ -419,7 +643,7 @@ pub const Mutable = struct { |
| 419 | 643 | /// `a` and `b` may alias with each other. |
| 420 | 644 | /// |
| 421 | 645 | /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by |
| 422 | | /// rma is given by `a.limbs.len + b.limbs.len + 1`. |
| 646 | /// rma is given by `a.limbs.len + b.limbs.len`. |
| 423 | 647 | /// |
| 424 | 648 | /// If `allocator` is provided, it will be used for temporary storage to improve |
| 425 | 649 | /// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm. |
| ... | ... | @@ -435,14 +659,89 @@ pub const Mutable = struct { |
| 435 | 659 | } |
| 436 | 660 | } |
| 437 | 661 | |
| 438 | | mem.set(Limb, rma.limbs[0 .. a.limbs.len + b.limbs.len + 1], 0); |
| 662 | mem.set(Limb, rma.limbs[0 .. a.limbs.len + b.limbs.len], 0); |
| 439 | 663 | |
| 440 | | llmulacc(allocator, rma.limbs, a.limbs, b.limbs); |
| 664 | llmulacc(.add, allocator, rma.limbs, a.limbs, b.limbs); |
| 441 | 665 | |
| 442 | 666 | rma.normalize(a.limbs.len + b.limbs.len); |
| 443 | 667 | rma.positive = (a.positive == b.positive); |
| 444 | 668 | } |
| 445 | 669 | |
| 670 | /// rma = a * b with 2s-complement wrapping semantics. |
| 671 | /// |
| 672 | /// `rma` may alias with `a` or `b`. |
| 673 | /// `a` and `b` may alias with each other. |
| 674 | /// |
| 675 | /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by |
| 676 | /// rma is given by `a.limbs.len + b.limbs.len`. |
| 677 | /// |
| 678 | /// `limbs_buffer` is used for temporary storage. The amount required is given by `calcMulWrapLimbsBufferLen`. |
| 679 | pub fn mulWrap( |
| 680 | rma: *Mutable, |
| 681 | a: Const, |
| 682 | b: Const, |
| 683 | signedness: std.builtin.Signedness, |
| 684 | bit_count: usize, |
| 685 | limbs_buffer: []Limb, |
| 686 | allocator: ?*Allocator, |
| 687 | ) void { |
| 688 | var buf_index: usize = 0; |
| 689 | const req_limbs = calcTwosCompLimbCount(bit_count); |
| 690 | |
| 691 | const a_copy = if (rma.limbs.ptr == a.limbs.ptr) blk: { |
| 692 | const start = buf_index; |
| 693 | const a_len = math.min(req_limbs, a.limbs.len); |
| 694 | mem.copy(Limb, limbs_buffer[buf_index..], a.limbs[0..a_len]); |
| 695 | buf_index += a_len; |
| 696 | break :blk a.toMutable(limbs_buffer[start..buf_index]).toConst(); |
| 697 | } else a; |
| 698 | |
| 699 | const b_copy = if (rma.limbs.ptr == b.limbs.ptr) blk: { |
| 700 | const start = buf_index; |
| 701 | const b_len = math.min(req_limbs, b.limbs.len); |
| 702 | mem.copy(Limb, limbs_buffer[buf_index..], b.limbs[0..b_len]); |
| 703 | buf_index += b_len; |
| 704 | break :blk a.toMutable(limbs_buffer[start..buf_index]).toConst(); |
| 705 | } else b; |
| 706 | |
| 707 | return rma.mulWrapNoAlias(a_copy, b_copy, signedness, bit_count, allocator); |
| 708 | } |
| 709 | |
| 710 | /// rma = a * b with 2s-complement wrapping semantics. |
| 711 | /// |
| 712 | /// `rma` may not alias with `a` or `b`. |
| 713 | /// `a` and `b` may alias with each other. |
| 714 | /// |
| 715 | /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by |
| 716 | /// rma is given by `a.limbs.len + b.limbs.len`. |
| 717 | /// |
| 718 | /// If `allocator` is provided, it will be used for temporary storage to improve |
| 719 | /// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm. |
| 720 | pub fn mulWrapNoAlias( |
| 721 | rma: *Mutable, |
| 722 | a: Const, |
| 723 | b: Const, |
| 724 | signedness: std.builtin.Signedness, |
| 725 | bit_count: usize, |
| 726 | allocator: ?*Allocator, |
| 727 | ) void { |
| 728 | assert(rma.limbs.ptr != a.limbs.ptr); // illegal aliasing |
| 729 | assert(rma.limbs.ptr != b.limbs.ptr); // illegal aliasing |
| 730 | |
| 731 | const req_limbs = calcTwosCompLimbCount(bit_count); |
| 732 | |
| 733 | // We can ignore the upper bits here, those results will be discarded anyway. |
| 734 | const a_limbs = a.limbs[0..math.min(req_limbs, a.limbs.len)]; |
| 735 | const b_limbs = b.limbs[0..math.min(req_limbs, b.limbs.len)]; |
| 736 | |
| 737 | mem.set(Limb, rma.limbs[0..req_limbs], 0); |
| 738 | |
| 739 | llmulacc(.add, allocator, rma.limbs, a_limbs, b_limbs); |
| 740 | rma.normalize(math.min(req_limbs, a.limbs.len + b.limbs.len)); |
| 741 | rma.positive = (a.positive == b.positive); |
| 742 | rma.truncate(rma.toConst(), signedness, bit_count); |
| 743 | } |
| 744 | |
| 446 | 745 | /// rma = a * a |
| 447 | 746 | /// |
| 448 | 747 | /// `rma` may not alias with `a`. |
| ... | ... | @@ -458,7 +757,7 @@ pub const Mutable = struct { |
| 458 | 757 | |
| 459 | 758 | mem.set(Limb, rma.limbs, 0); |
| 460 | 759 | |
| 461 | | llsquare_basecase(rma.limbs, a.limbs); |
| 760 | llsquareBasecase(rma.limbs, a.limbs); |
| 462 | 761 | |
| 463 | 762 | rma.normalize(2 * a.limbs.len + 1); |
| 464 | 763 | rma.positive = true; |
| ... | ... | @@ -980,6 +1279,102 @@ pub const Mutable = struct { |
| 980 | 1279 | r.normalize(r.len); |
| 981 | 1280 | } |
| 982 | 1281 | |
| 1282 | /// Truncate an integer to a number of bits, following 2s-complement semantics. |
| 1283 | /// r may alias a. |
| 1284 | /// |
| 1285 | /// Asserts `r` has enough storage to store the result. |
| 1286 | /// The upper bound is `calcTwosCompLimbCount(a.len)`. |
| 1287 | pub fn truncate(r: *Mutable, a: Const, signedness: std.builtin.Signedness, bit_count: usize) void { |
| 1288 | const req_limbs = calcTwosCompLimbCount(bit_count); |
| 1289 | |
| 1290 | // Handle 0-bit integers. |
| 1291 | if (req_limbs == 0 or a.eqZero()) { |
| 1292 | r.set(0); |
| 1293 | return; |
| 1294 | } |
| 1295 | |
| 1296 | const bit = @truncate(Log2Limb, bit_count - 1); |
| 1297 | const signmask = @as(Limb, 1) << bit; // 0b0..010...0 where 1 is the sign bit. |
| 1298 | const mask = (signmask << 1) -% 1; // 0b0..01..1 where the leftmost 1 is the sign bit. |
| 1299 | |
| 1300 | if (!a.positive) { |
| 1301 | // Convert the integer from sign-magnitude into twos-complement. |
| 1302 | // -x = ~(x - 1) |
| 1303 | // Note, we simply take req_limbs * @bitSizeOf(Limb) as the |
| 1304 | // target bit count. |
| 1305 | |
| 1306 | r.addScalar(a.abs(), -1); |
| 1307 | |
| 1308 | // Zero-extend the result |
| 1309 | if (req_limbs > r.len) { |
| 1310 | mem.set(Limb, r.limbs[r.len..req_limbs], 0); |
| 1311 | } |
| 1312 | |
| 1313 | // Truncate to required number of limbs. |
| 1314 | assert(r.limbs.len >= req_limbs); |
| 1315 | r.len = req_limbs; |
| 1316 | |
| 1317 | // Without truncating, we can already peek at the sign bit of the result here. |
| 1318 | // Note that it will be 0 if the result is negative, as we did not apply the flip here. |
| 1319 | // If the result is negative, we have |
| 1320 | // -(-x & mask) |
| 1321 | // = ~(~(x - 1) & mask) + 1 |
| 1322 | // = ~(~((x - 1) | ~mask)) + 1 |
| 1323 | // = ((x - 1) | ~mask)) + 1 |
| 1324 | // Note, this is only valid for the target bits and not the upper bits |
| 1325 | // of the most significant limb. Those still need to be cleared. |
| 1326 | // Also note that `mask` is zero for all other bits, reducing to the identity. |
| 1327 | // This means that we still need to use & mask to clear off the upper bits. |
| 1328 | |
| 1329 | if (signedness == .signed and r.limbs[r.len - 1] & signmask == 0) { |
| 1330 | // Re-add the one and negate to get the result. |
| 1331 | r.limbs[r.len - 1] &= mask; |
| 1332 | // Note, addition cannot require extra limbs here as we did a subtraction before. |
| 1333 | r.addScalar(r.toConst(), 1); |
| 1334 | r.normalize(r.len); |
| 1335 | r.positive = false; |
| 1336 | } else { |
| 1337 | llnot(r.limbs[0..r.len]); |
| 1338 | r.limbs[r.len - 1] &= mask; |
| 1339 | r.normalize(r.len); |
| 1340 | } |
| 1341 | } else { |
| 1342 | r.copy(a); |
| 1343 | if (r.len < req_limbs) { |
| 1344 | // Integer fits within target bits, no wrapping required. |
| 1345 | return; |
| 1346 | } |
| 1347 | |
| 1348 | r.len = req_limbs; |
| 1349 | r.limbs[r.len - 1] &= mask; |
| 1350 | r.normalize(r.len); |
| 1351 | |
| 1352 | if (signedness == .signed and r.limbs[r.len - 1] & signmask != 0) { |
| 1353 | // Convert 2s-complement back to sign-magnitude. |
| 1354 | // Sign-extend the upper bits so that they are inverted correctly. |
| 1355 | r.limbs[r.len - 1] |= ~mask; |
| 1356 | llnot(r.limbs[0..r.len]); |
| 1357 | |
| 1358 | // Note, can only overflow if r holds 0xFFF...F which can only happen if |
| 1359 | // a holds 0. |
| 1360 | r.addScalar(r.toConst(), 1); |
| 1361 | |
| 1362 | r.positive = false; |
| 1363 | } |
| 1364 | } |
| 1365 | } |
| 1366 | |
| 1367 | /// Saturate an integer to a number of bits, following 2s-complement semantics. |
| 1368 | /// r may alias a. |
| 1369 | /// |
| 1370 | /// Asserts `r` has enough storage to store the result. |
| 1371 | /// The upper bound is `calcTwosCompLimbCount(a.len)`. |
| 1372 | pub fn saturate(r: *Mutable, a: Const, signedness: std.builtin.Signedness, bit_count: usize) void { |
| 1373 | if (!a.fitsInTwosComp(signedness, bit_count)) { |
| 1374 | r.setTwosCompIntLimit(if (r.positive) .max else .min, signedness, bit_count); |
| 1375 | } |
| 1376 | } |
| 1377 | |
| 983 | 1378 | /// Normalize a possible sequence of leading zeros. |
| 984 | 1379 | /// |
| 985 | 1380 | /// [1, 2, 3, 4, 0] -> [1, 2, 3, 4] |
| ... | ... | @@ -1040,6 +1435,13 @@ pub const Const = struct { |
| 1040 | 1435 | }; |
| 1041 | 1436 | } |
| 1042 | 1437 | |
| 1438 | pub fn negate(self: Const) Const { |
| 1439 | return .{ |
| 1440 | .limbs = self.limbs, |
| 1441 | .positive = !self.positive, |
| 1442 | }; |
| 1443 | } |
| 1444 | |
| 1043 | 1445 | pub fn isOdd(self: Const) bool { |
| 1044 | 1446 | return self.limbs[0] & 1 != 0; |
| 1045 | 1447 | } |
| ... | ... | @@ -1643,6 +2045,21 @@ pub const Managed = struct { |
| 1643 | 2045 | self.setMetadata(m.positive, m.len); |
| 1644 | 2046 | } |
| 1645 | 2047 | |
| 2048 | /// Set self to either bound of a 2s-complement integer. |
| 2049 | /// Note: The result is still sign-magnitude, not twos complement! In order to convert the |
| 2050 | /// result to twos complement, it is sufficient to take the absolute value. |
| 2051 | pub fn setTwosCompIntLimit( |
| 2052 | r: *Managed, |
| 2053 | limit: TwosCompIntLimit, |
| 2054 | signedness: std.builtin.Signedness, |
| 2055 | bit_count: usize, |
| 2056 | ) !void { |
| 2057 | try r.ensureCapacity(calcTwosCompLimbCount(bit_count)); |
| 2058 | var m = r.toMutable(); |
| 2059 | m.setTwosCompIntLimit(limit, signedness, bit_count); |
| 2060 | r.setMetadata(m.positive, m.len); |
| 2061 | } |
| 2062 | |
| 1646 | 2063 | /// Converts self to a string in the requested base. Memory is allocated from the provided |
| 1647 | 2064 | /// allocator and not the one present in self. |
| 1648 | 2065 | pub fn toString(self: Managed, allocator: *Allocator, base: u8, case: std.fmt.Case) ![]u8 { |
| ... | ... | @@ -1741,6 +2158,32 @@ pub const Managed = struct { |
| 1741 | 2158 | r.setMetadata(m.positive, m.len); |
| 1742 | 2159 | } |
| 1743 | 2160 | |
| 2161 | /// r = a + b with 2s-complement wrapping semantics. |
| 2162 | /// |
| 2163 | /// r, a and b may be aliases. If r aliases a or b, then caller must call |
| 2164 | /// `r.ensureTwosCompCapacity` prior to calling `add`. |
| 2165 | /// |
| 2166 | /// Returns an error if memory could not be allocated. |
| 2167 | pub fn addWrap(r: *Managed, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) Allocator.Error!void { |
| 2168 | try r.ensureTwosCompCapacity(bit_count); |
| 2169 | var m = r.toMutable(); |
| 2170 | m.addWrap(a, b, signedness, bit_count); |
| 2171 | r.setMetadata(m.positive, m.len); |
| 2172 | } |
| 2173 | |
| 2174 | /// r = a + b with 2s-complement saturating semantics. |
| 2175 | /// |
| 2176 | /// r, a and b may be aliases. If r aliases a or b, then caller must call |
| 2177 | /// `r.ensureTwosCompCapacity` prior to calling `add`. |
| 2178 | /// |
| 2179 | /// Returns an error if memory could not be allocated. |
| 2180 | pub fn addSat(r: *Managed, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) Allocator.Error!void { |
| 2181 | try r.ensureTwosCompCapacity(bit_count); |
| 2182 | var m = r.toMutable(); |
| 2183 | m.addSat(a, b, signedness, bit_count); |
| 2184 | r.setMetadata(m.positive, m.len); |
| 2185 | } |
| 2186 | |
| 1744 | 2187 | /// r = a - b |
| 1745 | 2188 | /// |
| 1746 | 2189 | /// r, a and b may be aliases. |
| ... | ... | @@ -1753,6 +2196,32 @@ pub const Managed = struct { |
| 1753 | 2196 | r.setMetadata(m.positive, m.len); |
| 1754 | 2197 | } |
| 1755 | 2198 | |
| 2199 | /// r = a - b with 2s-complement wrapping semantics. |
| 2200 | /// |
| 2201 | /// r, a and b may be aliases. If r aliases a or b, then caller must call |
| 2202 | /// `r.ensureTwosCompCapacity` prior to calling `add`. |
| 2203 | /// |
| 2204 | /// Returns an error if memory could not be allocated. |
| 2205 | pub fn subWrap(r: *Managed, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) Allocator.Error!void { |
| 2206 | try r.ensureTwosCompCapacity(bit_count); |
| 2207 | var m = r.toMutable(); |
| 2208 | m.subWrap(a, b, signedness, bit_count); |
| 2209 | r.setMetadata(m.positive, m.len); |
| 2210 | } |
| 2211 | |
| 2212 | /// r = a - b with 2s-complement saturating semantics. |
| 2213 | /// |
| 2214 | /// r, a and b may be aliases. If r aliases a or b, then caller must call |
| 2215 | /// `r.ensureTwosCompCapacity` prior to calling `add`. |
| 2216 | /// |
| 2217 | /// Returns an error if memory could not be allocated. |
| 2218 | pub fn subSat(r: *Managed, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) Allocator.Error!void { |
| 2219 | try r.ensureTwosCompCapacity(bit_count); |
| 2220 | var m = r.toMutable(); |
| 2221 | m.subSat(a, b, signedness, bit_count); |
| 2222 | r.setMetadata(m.positive, m.len); |
| 2223 | } |
| 2224 | |
| 1756 | 2225 | /// rma = a * b |
| 1757 | 2226 | /// |
| 1758 | 2227 | /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b. |
| ... | ... | @@ -1781,6 +2250,39 @@ pub const Managed = struct { |
| 1781 | 2250 | rma.setMetadata(m.positive, m.len); |
| 1782 | 2251 | } |
| 1783 | 2252 | |
| 2253 | /// rma = a * b with 2s-complement wrapping semantics. |
| 2254 | /// |
| 2255 | /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b. |
| 2256 | /// If rma aliases a or b, then caller must call `ensureTwosCompCapacity` |
| 2257 | /// prior to calling `mul`. |
| 2258 | /// |
| 2259 | /// Returns an error if memory could not be allocated. |
| 2260 | /// |
| 2261 | /// rma's allocator is used for temporary storage to speed up the multiplication. |
| 2262 | pub fn mulWrap(rma: *Managed, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) !void { |
| 2263 | var alias_count: usize = 0; |
| 2264 | if (rma.limbs.ptr == a.limbs.ptr) |
| 2265 | alias_count += 1; |
| 2266 | if (rma.limbs.ptr == b.limbs.ptr) |
| 2267 | alias_count += 1; |
| 2268 | |
| 2269 | try rma.ensureTwosCompCapacity(bit_count); |
| 2270 | var m = rma.toMutable(); |
| 2271 | if (alias_count == 0) { |
| 2272 | m.mulWrapNoAlias(a, b, signedness, bit_count, rma.allocator); |
| 2273 | } else { |
| 2274 | const limb_count = calcMulWrapLimbsBufferLen(bit_count, a.limbs.len, b.limbs.len, alias_count); |
| 2275 | const limbs_buffer = try rma.allocator.alloc(Limb, limb_count); |
| 2276 | defer rma.allocator.free(limbs_buffer); |
| 2277 | m.mulWrap(a, b, signedness, bit_count, limbs_buffer, rma.allocator); |
| 2278 | } |
| 2279 | rma.setMetadata(m.positive, m.len); |
| 2280 | } |
| 2281 | |
| 2282 | pub fn ensureTwosCompCapacity(r: *Managed, bit_count: usize) !void { |
| 2283 | try r.ensureCapacity(calcTwosCompLimbCount(bit_count)); |
| 2284 | } |
| 2285 | |
| 1784 | 2286 | pub fn ensureAddScalarCapacity(r: *Managed, a: Const, scalar: anytype) !void { |
| 1785 | 2287 | try r.ensureCapacity(math.max(a.limbs.len, calcLimbLen(scalar)) + 1); |
| 1786 | 2288 | } |
| ... | ... | @@ -1941,30 +2443,58 @@ pub const Managed = struct { |
| 1941 | 2443 | rma.setMetadata(rma_mut.positive, rma_mut.len); |
| 1942 | 2444 | } |
| 1943 | 2445 | } |
| 2446 | |
| 2447 | /// r = truncate(Int(signedness, bit_count), a) |
| 2448 | pub fn truncate(r: *Managed, a: Const, signedness: std.builtin.Signedness, bit_count: usize) !void { |
| 2449 | try r.ensureCapacity(calcTwosCompLimbCount(bit_count)); |
| 2450 | var m = r.toMutable(); |
| 2451 | m.truncate(a, signedness, bit_count); |
| 2452 | r.setMetadata(m.positive, m.len); |
| 2453 | } |
| 2454 | |
| 2455 | /// r = saturate(Int(signedness, bit_count), a) |
| 2456 | pub fn saturate(r: *Managed, a: Const, signedness: std.builtin.Signedness, bit_count: usize) !void { |
| 2457 | try r.ensureCapacity(calcTwosCompLimbCount(bit_count)); |
| 2458 | var m = r.toMutable(); |
| 2459 | m.saturate(a, signedness, bit_count); |
| 2460 | r.setMetadata(m.positive, m.len); |
| 2461 | } |
| 2462 | }; |
| 2463 | |
| 2464 | /// Different operators which can be used in accumulation style functions |
| 2465 | /// (llmulacc, llmulaccKaratsuba, llmulaccLong, llmulLimb). In all these functions, |
| 2466 | /// a computed value is accumulated with an existing result. |
| 2467 | const AccOp = enum { |
| 2468 | /// The computed value is added to the result. |
| 2469 | add, |
| 2470 | |
| 2471 | /// The computed value is subtracted from the result. |
| 2472 | sub, |
| 1944 | 2473 | }; |
| 1945 | 2474 | |
| 1946 | 2475 | /// Knuth 4.3.1, Algorithm M. |
| 1947 | 2476 | /// |
| 2477 | /// r = r (op) a * b |
| 1948 | 2478 | /// r MUST NOT alias any of a or b. |
| 1949 | | fn llmulacc(opt_allocator: ?*Allocator, r: []Limb, a: []const Limb, b: []const Limb) void { |
| 2479 | /// |
| 2480 | /// The result is computed modulo `r.len`. When `r.len >= a.len + b.len`, no overflow occurs. |
| 2481 | fn llmulacc(comptime op: AccOp, opt_allocator: ?*Allocator, r: []Limb, a: []const Limb, b: []const Limb) void { |
| 1950 | 2482 | @setRuntimeSafety(debug_safety); |
| 2483 | assert(r.len >= a.len); |
| 2484 | assert(r.len >= b.len); |
| 1951 | 2485 | |
| 1952 | | const a_norm = a[0..llnormalize(a)]; |
| 1953 | | const b_norm = b[0..llnormalize(b)]; |
| 1954 | | var x = a_norm; |
| 1955 | | var y = b_norm; |
| 1956 | | if (a_norm.len > b_norm.len) { |
| 1957 | | x = b_norm; |
| 1958 | | y = a_norm; |
| 2486 | // Order greatest first. |
| 2487 | var x = a; |
| 2488 | var y = b; |
| 2489 | if (a.len < b.len) { |
| 2490 | x = b; |
| 2491 | y = a; |
| 1959 | 2492 | } |
| 1960 | 2493 | |
| 1961 | | assert(r.len >= x.len + y.len + 1); |
| 1962 | | |
| 1963 | | // 48 is a pretty abitrary size chosen based on performance of a factorial program. |
| 1964 | 2494 | k_mul: { |
| 1965 | | if (x.len > 48) { |
| 2495 | if (y.len > 48) { |
| 1966 | 2496 | if (opt_allocator) |allocator| { |
| 1967 | | llmulacc_karatsuba(allocator, r, x, y) catch |err| switch (err) { |
| 2497 | llmulaccKaratsuba(op, allocator, r, x, y) catch |err| switch (err) { |
| 1968 | 2498 | error.OutOfMemory => break :k_mul, // handled below |
| 1969 | 2499 | }; |
| 1970 | 2500 | return; |
| ... | ... | @@ -1972,83 +2502,191 @@ fn llmulacc(opt_allocator: ?*Allocator, r: []Limb, a: []const Limb, b: []const L |
| 1972 | 2502 | } |
| 1973 | 2503 | } |
| 1974 | 2504 | |
| 1975 | | // Basecase multiplication |
| 1976 | | var i: usize = 0; |
| 1977 | | while (i < x.len) : (i += 1) { |
| 1978 | | llmulDigit(r[i..], y, x[i]); |
| 1979 | | } |
| 2505 | llmulaccLong(op, r, x, y); |
| 1980 | 2506 | } |
| 1981 | 2507 | |
| 1982 | 2508 | /// Knuth 4.3.1, Algorithm M. |
| 1983 | 2509 | /// |
| 2510 | /// r = r (op) a * b |
| 1984 | 2511 | /// r MUST NOT alias any of a or b. |
| 1985 | | fn llmulacc_karatsuba(allocator: *Allocator, r: []Limb, x: []const Limb, y: []const Limb) error{OutOfMemory}!void { |
| 2512 | /// |
| 2513 | /// The result is computed modulo `r.len`. When `r.len >= a.len + b.len`, no overflow occurs. |
| 2514 | fn llmulaccKaratsuba( |
| 2515 | comptime op: AccOp, |
| 2516 | allocator: *Allocator, |
| 2517 | r: []Limb, |
| 2518 | a: []const Limb, |
| 2519 | b: []const Limb, |
| 2520 | ) error{OutOfMemory}!void { |
| 1986 | 2521 | @setRuntimeSafety(debug_safety); |
| 2522 | assert(r.len >= a.len); |
| 2523 | assert(a.len >= b.len); |
| 1987 | 2524 | |
| 1988 | | assert(r.len >= x.len + y.len + 1); |
| 2525 | // Classical karatsuba algorithm: |
| 2526 | // a = a1 * B + a0 |
| 2527 | // b = b1 * B + b0 |
| 2528 | // Where a0, b0 < B |
| 2529 | // |
| 2530 | // We then have: |
| 2531 | // ab = a * b |
| 2532 | // = (a1 * B + a0) * (b1 * B + b0) |
| 2533 | // = a1 * b1 * B * B + a1 * B * b0 + a0 * b1 * B + a0 * b0 |
| 2534 | // = a1 * b1 * B * B + (a1 * b0 + a0 * b1) * B + a0 * b0 |
| 2535 | // |
| 2536 | // Note that: |
| 2537 | // a1 * b0 + a0 * b1 |
| 2538 | // = (a1 + a0)(b1 + b0) - a1 * b1 - a0 * b0 |
| 2539 | // = (a0 - a1)(b1 - b0) + a1 * b1 + a0 * b0 |
| 2540 | // |
| 2541 | // This yields: |
| 2542 | // ab = p2 * B^2 + (p0 + p1 + p2) * B + p0 |
| 2543 | // |
| 2544 | // Where: |
| 2545 | // p0 = a0 * b0 |
| 2546 | // p1 = (a0 - a1)(b1 - b0) |
| 2547 | // p2 = a1 * b1 |
| 2548 | // |
| 2549 | // Note, (a0 - a1) and (b1 - b0) produce values -B < x < B, and so we need to mind the sign here. |
| 2550 | // We also have: |
| 2551 | // 0 <= p0 <= 2B |
| 2552 | // -2B <= p1 <= 2B |
| 2553 | // |
| 2554 | // Note, when B is a multiple of the limb size, multiplies by B amount to shifts or |
| 2555 | // slices of a limbs array. |
| 2556 | // |
| 2557 | // This function computes the result of the multiplication modulo r.len. This means: |
| 2558 | // - p2 and p1 only need to be computed modulo r.len - B. |
| 2559 | // - In the case of p2, p2 * B^2 needs to be added modulo r.len - 2 * B. |
| 2560 | |
| 2561 | const split = b.len / 2; // B |
| 1989 | 2562 | |
| 1990 | | const split = @divFloor(x.len, 2); |
| 1991 | | var x0 = x[0..split]; |
| 1992 | | var x1 = x[split..x.len]; |
| 1993 | | var y0 = y[0..split]; |
| 1994 | | var y1 = y[split..y.len]; |
| 2563 | const limbs_after_split = r.len - split; // Limbs to compute for p1 and p2. |
| 2564 | const limbs_after_split2 = r.len - split * 2; // Limbs to add for p2 * B^2. |
| 2565 | |
| 2566 | // For a0 and b0 we need the full range. |
| 2567 | const a0 = a[0..llnormalize(a[0..split])]; |
| 2568 | const b0 = b[0..llnormalize(b[0..split])]; |
| 2569 | |
| 2570 | // For a1 and b1 we only need `limbs_after_split` limbs. |
| 2571 | const a1 = blk: { |
| 2572 | var a1 = a[split..]; |
| 2573 | a1.len = math.min(llnormalize(a1), limbs_after_split); |
| 2574 | break :blk a1; |
| 2575 | }; |
| 1995 | 2576 | |
| 1996 | | var tmp = try allocator.alloc(Limb, x1.len + y1.len + 1); |
| 2577 | const b1 = blk: { |
| 2578 | var b1 = b[split..]; |
| 2579 | b1.len = math.min(llnormalize(b1), limbs_after_split); |
| 2580 | break :blk b1; |
| 2581 | }; |
| 2582 | |
| 2583 | // Note that the above slices relative to `split` work because we have a.len > b.len. |
| 2584 | |
| 2585 | // We need some temporary memory to store intermediate results. |
| 2586 | // Note, we can reduce the amount of temporaries we need by reordering the computation here: |
| 2587 | // ab = p2 * B^2 + (p0 + p1 + p2) * B + p0 |
| 2588 | // = p2 * B^2 + (p0 * B + p1 * B + p2 * B) + p0 |
| 2589 | // = (p2 * B^2 + p2 * B) + (p0 * B + p0) + p1 * B |
| 2590 | |
| 2591 | // Allocate at least enough memory to be able to multiply the upper two segments of a and b, assuming |
| 2592 | // no overflow. |
| 2593 | const tmp = try allocator.alloc(Limb, a.len - split + b.len - split); |
| 1997 | 2594 | defer allocator.free(tmp); |
| 1998 | | mem.set(Limb, tmp, 0); |
| 1999 | 2595 | |
| 2000 | | llmulacc(allocator, tmp, x1, y1); |
| 2596 | // Compute p2. |
| 2597 | // Note, we don't need to compute all of p2, just enough limbs to satisfy r. |
| 2598 | const p2_limbs = math.min(limbs_after_split, a1.len + b1.len); |
| 2001 | 2599 | |
| 2002 | | var length = llnormalize(tmp); |
| 2003 | | _ = llaccum(r[split..], tmp[0..length]); |
| 2004 | | _ = llaccum(r[split * 2 ..], tmp[0..length]); |
| 2600 | mem.set(Limb, tmp[0..p2_limbs], 0); |
| 2601 | llmulacc(.add, allocator, tmp[0..p2_limbs], a1[0..math.min(a1.len, p2_limbs)], b1[0..math.min(b1.len, p2_limbs)]); |
| 2602 | const p2 = tmp[0..llnormalize(tmp[0..p2_limbs])]; |
| 2005 | 2603 | |
| 2006 | | mem.set(Limb, tmp[0..length], 0); |
| 2604 | // Add p2 * B to the result. |
| 2605 | llaccum(op, r[split..], p2); |
| 2007 | 2606 | |
| 2008 | | llmulacc(allocator, tmp, x0, y0); |
| 2607 | // Add p2 * B^2 to the result if required. |
| 2608 | if (limbs_after_split2 > 0) { |
| 2609 | llaccum(op, r[split * 2 ..], p2[0..math.min(p2.len, limbs_after_split2)]); |
| 2610 | } |
| 2611 | |
| 2612 | // Compute p0. |
| 2613 | // Since a0.len, b0.len <= split and r.len >= split * 2, the full width of p0 needs to be computed. |
| 2614 | const p0_limbs = a0.len + b0.len; |
| 2615 | mem.set(Limb, tmp[0..p0_limbs], 0); |
| 2616 | llmulacc(.add, allocator, tmp[0..p0_limbs], a0, b0); |
| 2617 | const p0 = tmp[0..llnormalize(tmp[0..p0_limbs])]; |
| 2618 | |
| 2619 | // Add p0 to the result. |
| 2620 | llaccum(op, r, p0); |
| 2009 | 2621 | |
| 2010 | | length = llnormalize(tmp); |
| 2011 | | _ = llaccum(r[0..], tmp[0..length]); |
| 2012 | | _ = llaccum(r[split..], tmp[0..length]); |
| 2622 | // Add p0 * B to the result. In this case, we may not need all of it. |
| 2623 | llaccum(op, r[split..], p0[0..math.min(limbs_after_split, p0.len)]); |
| 2013 | 2624 | |
| 2014 | | const x_cmp = llcmp(x1, x0); |
| 2015 | | const y_cmp = llcmp(y1, y0); |
| 2016 | | if (x_cmp * y_cmp == 0) { |
| 2625 | // Finally, compute and add p1. |
| 2626 | // From now on we only need `limbs_after_split` limbs for a0 and b0, since the result of the |
| 2627 | // following computation will be added * B. |
| 2628 | const a0x = a0[0..std.math.min(a0.len, limbs_after_split)]; |
| 2629 | const b0x = b0[0..std.math.min(b0.len, limbs_after_split)]; |
| 2630 | |
| 2631 | const j0_sign = llcmp(a0x, a1); |
| 2632 | const j1_sign = llcmp(b1, b0x); |
| 2633 | |
| 2634 | if (j0_sign * j1_sign == 0) { |
| 2635 | // p1 is zero, we don't need to do any computation at all. |
| 2017 | 2636 | return; |
| 2018 | 2637 | } |
| 2019 | | const x0_len = llnormalize(x0); |
| 2020 | | const x1_len = llnormalize(x1); |
| 2021 | | var j0 = try allocator.alloc(Limb, math.max(x0_len, x1_len)); |
| 2022 | | defer allocator.free(j0); |
| 2023 | | if (x_cmp == 1) { |
| 2024 | | llsub(j0, x1[0..x1_len], x0[0..x0_len]); |
| 2638 | |
| 2639 | mem.set(Limb, tmp, 0); |
| 2640 | |
| 2641 | // p1 is nonzero, so compute the intermediary terms j0 = a0 - a1 and j1 = b1 - b0. |
| 2642 | // Note that in this case, we again need some storage for intermediary results |
| 2643 | // j0 and j1. Since we have tmp.len >= 2B, we can store both |
| 2644 | // intermediaries in the already allocated array. |
| 2645 | const j0 = tmp[0 .. a.len - split]; |
| 2646 | const j1 = tmp[a.len - split ..]; |
| 2647 | |
| 2648 | // Ensure that no subtraction overflows. |
| 2649 | if (j0_sign == 1) { |
| 2650 | // a0 > a1. |
| 2651 | _ = llsubcarry(j0, a0x, a1); |
| 2025 | 2652 | } else { |
| 2026 | | llsub(j0, x0[0..x0_len], x1[0..x1_len]); |
| 2653 | // a0 < a1. |
| 2654 | _ = llsubcarry(j0, a1, a0x); |
| 2027 | 2655 | } |
| 2028 | 2656 | |
| 2029 | | const y0_len = llnormalize(y0); |
| 2030 | | const y1_len = llnormalize(y1); |
| 2031 | | var j1 = try allocator.alloc(Limb, math.max(y0_len, y1_len)); |
| 2032 | | defer allocator.free(j1); |
| 2033 | | if (y_cmp == 1) { |
| 2034 | | llsub(j1, y1[0..y1_len], y0[0..y0_len]); |
| 2657 | if (j1_sign == 1) { |
| 2658 | // b1 > b0. |
| 2659 | _ = llsubcarry(j1, b1, b0x); |
| 2035 | 2660 | } else { |
| 2036 | | llsub(j1, y0[0..y0_len], y1[0..y1_len]); |
| 2661 | // b1 > b0. |
| 2662 | _ = llsubcarry(j1, b0x, b1); |
| 2037 | 2663 | } |
| 2038 | | if (x_cmp == y_cmp) { |
| 2039 | | mem.set(Limb, tmp[0..length], 0); |
| 2040 | | llmulacc(allocator, tmp, j0, j1); |
| 2041 | 2664 | |
| 2042 | | length = llnormalize(tmp); |
| 2043 | | llsub(r[split..], r[split..], tmp[0..length]); |
| 2665 | if (j0_sign * j1_sign == 1) { |
| 2666 | // If j0 and j1 are both positive, we now have: |
| 2667 | // p1 = j0 * j1 |
| 2668 | // If j0 and j1 are both negative, we now have: |
| 2669 | // p1 = -j0 * -j1 = j0 * j1 |
| 2670 | // In this case we can add p1 to the result using llmulacc. |
| 2671 | llmulacc(op, allocator, r[split..], j0[0..llnormalize(j0)], j1[0..llnormalize(j1)]); |
| 2044 | 2672 | } else { |
| 2045 | | llmulacc(allocator, r[split..], j0, j1); |
| 2673 | // In this case either j0 or j1 is negative, an we have: |
| 2674 | // p1 = -(j0 * j1) |
| 2675 | // Now we need to subtract instead of accumulate. |
| 2676 | const inverted_op = if (op == .add) .sub else .add; |
| 2677 | llmulacc(inverted_op, allocator, r[split..], j0[0..llnormalize(j0)], j1[0..llnormalize(j1)]); |
| 2046 | 2678 | } |
| 2047 | 2679 | } |
| 2048 | 2680 | |
| 2049 | | // r = r + a |
| 2050 | | fn llaccum(r: []Limb, a: []const Limb) Limb { |
| 2681 | /// r = r (op) a. |
| 2682 | /// The result is computed modulo `r.len`. |
| 2683 | fn llaccum(comptime op: AccOp, r: []Limb, a: []const Limb) void { |
| 2051 | 2684 | @setRuntimeSafety(debug_safety); |
| 2685 | if (op == .sub) { |
| 2686 | _ = llsubcarry(r, r, a); |
| 2687 | return; |
| 2688 | } |
| 2689 | |
| 2052 | 2690 | assert(r.len != 0 and a.len != 0); |
| 2053 | 2691 | assert(r.len >= a.len); |
| 2054 | 2692 | |
| ... | ... | @@ -2065,8 +2703,6 @@ fn llaccum(r: []Limb, a: []const Limb) Limb { |
| 2065 | 2703 | while ((carry != 0) and i < r.len) : (i += 1) { |
| 2066 | 2704 | carry = @boolToInt(@addWithOverflow(Limb, r[i], carry, &r[i])); |
| 2067 | 2705 | } |
| 2068 | | |
| 2069 | | return carry; |
| 2070 | 2706 | } |
| 2071 | 2707 | |
| 2072 | 2708 | /// Returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively for limbs. |
| ... | ... | @@ -2097,24 +2733,55 @@ pub fn llcmp(a: []const Limb, b: []const Limb) i8 { |
| 2097 | 2733 | } |
| 2098 | 2734 | } |
| 2099 | 2735 | |
| 2100 | | fn llmulDigit(acc: []Limb, y: []const Limb, xi: Limb) void { |
| 2736 | /// r = r (op) y * xi |
| 2737 | /// The result is computed modulo `r.len`. When `r.len >= a.len + b.len`, no overflow occurs. |
| 2738 | fn llmulaccLong(comptime op: AccOp, r: []Limb, a: []const Limb, b: []const Limb) void { |
| 2739 | @setRuntimeSafety(debug_safety); |
| 2740 | assert(r.len >= a.len); |
| 2741 | assert(a.len >= b.len); |
| 2742 | |
| 2743 | var i: usize = 0; |
| 2744 | while (i < b.len) : (i += 1) { |
| 2745 | llmulLimb(op, r[i..], a, b[i]); |
| 2746 | } |
| 2747 | } |
| 2748 | |
| 2749 | /// r = r (op) y * xi |
| 2750 | /// The result is computed modulo `r.len`. |
| 2751 | fn llmulLimb(comptime op: AccOp, acc: []Limb, y: []const Limb, xi: Limb) void { |
| 2101 | 2752 | @setRuntimeSafety(debug_safety); |
| 2102 | 2753 | if (xi == 0) { |
| 2103 | 2754 | return; |
| 2104 | 2755 | } |
| 2105 | 2756 | |
| 2106 | | var carry: Limb = 0; |
| 2107 | 2757 | var a_lo = acc[0..y.len]; |
| 2108 | 2758 | var a_hi = acc[y.len..]; |
| 2109 | 2759 | |
| 2110 | | var j: usize = 0; |
| 2111 | | while (j < a_lo.len) : (j += 1) { |
| 2112 | | a_lo[j] = @call(.{ .modifier = .always_inline }, addMulLimbWithCarry, .{ a_lo[j], y[j], xi, &carry }); |
| 2113 | | } |
| 2760 | switch (op) { |
| 2761 | .add => { |
| 2762 | var carry: Limb = 0; |
| 2763 | var j: usize = 0; |
| 2764 | while (j < a_lo.len) : (j += 1) { |
| 2765 | a_lo[j] = addMulLimbWithCarry(a_lo[j], y[j], xi, &carry); |
| 2766 | } |
| 2114 | 2767 | |
| 2115 | | j = 0; |
| 2116 | | while ((carry != 0) and (j < a_hi.len)) : (j += 1) { |
| 2117 | | carry = @boolToInt(@addWithOverflow(Limb, a_hi[j], carry, &a_hi[j])); |
| 2768 | j = 0; |
| 2769 | while ((carry != 0) and (j < a_hi.len)) : (j += 1) { |
| 2770 | carry = @boolToInt(@addWithOverflow(Limb, a_hi[j], carry, &a_hi[j])); |
| 2771 | } |
| 2772 | }, |
| 2773 | .sub => { |
| 2774 | var borrow: Limb = 0; |
| 2775 | var j: usize = 0; |
| 2776 | while (j < a_lo.len) : (j += 1) { |
| 2777 | a_lo[j] = subMulLimbWithBorrow(a_lo[j], y[j], xi, &borrow); |
| 2778 | } |
| 2779 | |
| 2780 | j = 0; |
| 2781 | while ((borrow != 0) and (j < a_hi.len)) : (j += 1) { |
| 2782 | borrow = @boolToInt(@subWithOverflow(Limb, a_hi[j], borrow, &a_hi[j])); |
| 2783 | } |
| 2784 | }, |
| 2118 | 2785 | } |
| 2119 | 2786 | } |
| 2120 | 2787 | |
| ... | ... | @@ -2133,10 +2800,10 @@ fn llnormalize(a: []const Limb) usize { |
| 2133 | 2800 | } |
| 2134 | 2801 | |
| 2135 | 2802 | /// Knuth 4.3.1, Algorithm S. |
| 2136 | | fn llsub(r: []Limb, a: []const Limb, b: []const Limb) void { |
| 2803 | fn llsubcarry(r: []Limb, a: []const Limb, b: []const Limb) Limb { |
| 2137 | 2804 | @setRuntimeSafety(debug_safety); |
| 2138 | 2805 | assert(a.len != 0 and b.len != 0); |
| 2139 | | assert(a.len > b.len or (a.len == b.len and a[a.len - 1] >= b[b.len - 1])); |
| 2806 | assert(a.len >= b.len); |
| 2140 | 2807 | assert(r.len >= a.len); |
| 2141 | 2808 | |
| 2142 | 2809 | var i: usize = 0; |
| ... | ... | @@ -2153,15 +2820,21 @@ fn llsub(r: []Limb, a: []const Limb, b: []const Limb) void { |
| 2153 | 2820 | borrow = @boolToInt(@subWithOverflow(Limb, a[i], borrow, &r[i])); |
| 2154 | 2821 | } |
| 2155 | 2822 | |
| 2156 | | assert(borrow == 0); |
| 2823 | return borrow; |
| 2824 | } |
| 2825 | |
| 2826 | fn llsub(r: []Limb, a: []const Limb, b: []const Limb) void { |
| 2827 | @setRuntimeSafety(debug_safety); |
| 2828 | assert(a.len > b.len or (a.len == b.len and a[a.len - 1] >= b[b.len - 1])); |
| 2829 | assert(llsubcarry(r, a, b) == 0); |
| 2157 | 2830 | } |
| 2158 | 2831 | |
| 2159 | 2832 | /// Knuth 4.3.1, Algorithm A. |
| 2160 | | fn lladd(r: []Limb, a: []const Limb, b: []const Limb) void { |
| 2833 | fn lladdcarry(r: []Limb, a: []const Limb, b: []const Limb) Limb { |
| 2161 | 2834 | @setRuntimeSafety(debug_safety); |
| 2162 | 2835 | assert(a.len != 0 and b.len != 0); |
| 2163 | 2836 | assert(a.len >= b.len); |
| 2164 | | assert(r.len >= a.len + 1); |
| 2837 | assert(r.len >= a.len); |
| 2165 | 2838 | |
| 2166 | 2839 | var i: usize = 0; |
| 2167 | 2840 | var carry: Limb = 0; |
| ... | ... | @@ -2177,7 +2850,13 @@ fn lladd(r: []Limb, a: []const Limb, b: []const Limb) void { |
| 2177 | 2850 | carry = @boolToInt(@addWithOverflow(Limb, a[i], carry, &r[i])); |
| 2178 | 2851 | } |
| 2179 | 2852 | |
| 2180 | | r[i] = carry; |
| 2853 | return carry; |
| 2854 | } |
| 2855 | |
| 2856 | fn lladd(r: []Limb, a: []const Limb, b: []const Limb) void { |
| 2857 | @setRuntimeSafety(debug_safety); |
| 2858 | assert(r.len >= a.len + 1); |
| 2859 | r[a.len] = lladdcarry(r, a, b); |
| 2181 | 2860 | } |
| 2182 | 2861 | |
| 2183 | 2862 | /// Knuth 4.3.1, Exercise 16. |
| ... | ... | @@ -2258,6 +2937,15 @@ fn llshr(r: []Limb, a: []const Limb, shift: usize) void { |
| 2258 | 2937 | } |
| 2259 | 2938 | } |
| 2260 | 2939 | |
| 2940 | // r = ~r |
| 2941 | fn llnot(r: []Limb) void { |
| 2942 | @setRuntimeSafety(debug_safety); |
| 2943 | |
| 2944 | for (r) |*elem| { |
| 2945 | elem.* = ~elem.*; |
| 2946 | } |
| 2947 | } |
| 2948 | |
| 2261 | 2949 | // r = a | b with 2s complement semantics. |
| 2262 | 2950 | // r may alias. |
| 2263 | 2951 | // a and b must not be 0. |
| ... | ... | @@ -2554,7 +3242,7 @@ fn llsignedxor(r: []Limb, a: []const Limb, a_positive: bool, b: []const Limb, b_ |
| 2554 | 3242 | } |
| 2555 | 3243 | |
| 2556 | 3244 | /// r MUST NOT alias x. |
| 2557 | | fn llsquare_basecase(r: []Limb, x: []const Limb) void { |
| 3245 | fn llsquareBasecase(r: []Limb, x: []const Limb) void { |
| 2558 | 3246 | @setRuntimeSafety(debug_safety); |
| 2559 | 3247 | |
| 2560 | 3248 | const x_norm = x; |
| ... | ... | @@ -2577,7 +3265,7 @@ fn llsquare_basecase(r: []Limb, x: []const Limb) void { |
| 2577 | 3265 | |
| 2578 | 3266 | for (x_norm) |v, i| { |
| 2579 | 3267 | // Accumulate all the x[i]*x[j] (with x!=j) products |
| 2580 | | llmulDigit(r[2 * i + 1 ..], x_norm[i + 1 ..], v); |
| 3268 | llmulLimb(.add, r[2 * i + 1 ..], x_norm[i + 1 ..], v); |
| 2581 | 3269 | } |
| 2582 | 3270 | |
| 2583 | 3271 | // Each product appears twice, multiply by 2 |
| ... | ... | @@ -2585,7 +3273,7 @@ fn llsquare_basecase(r: []Limb, x: []const Limb) void { |
| 2585 | 3273 | |
| 2586 | 3274 | for (x_norm) |v, i| { |
| 2587 | 3275 | // Compute and add the squares |
| 2588 | | llmulDigit(r[2 * i ..], x[i .. i + 1], v); |
| 3276 | llmulLimb(.add, r[2 * i ..], x[i .. i + 1], v); |
| 2589 | 3277 | } |
| 2590 | 3278 | } |
| 2591 | 3279 | |
| ... | ... | @@ -2624,12 +3312,12 @@ fn llpow(r: []Limb, a: []const Limb, b: u32, tmp_limbs: []Limb) void { |
| 2624 | 3312 | while (i < exp_bits) : (i += 1) { |
| 2625 | 3313 | // Square |
| 2626 | 3314 | mem.set(Limb, tmp2, 0); |
| 2627 | | llsquare_basecase(tmp2, tmp1[0..llnormalize(tmp1)]); |
| 3315 | llsquareBasecase(tmp2, tmp1[0..llnormalize(tmp1)]); |
| 2628 | 3316 | mem.swap([]Limb, &tmp1, &tmp2); |
| 2629 | 3317 | // Multiply by a |
| 2630 | 3318 | if (@shlWithOverflow(u32, exp, 1, &exp)) { |
| 2631 | 3319 | mem.set(Limb, tmp2, 0); |
| 2632 | | llmulacc(null, tmp2, tmp1[0..llnormalize(tmp1)], a); |
| 3320 | llmulacc(.add, null, tmp2, tmp1[0..llnormalize(tmp1)], a); |
| 2633 | 3321 | mem.swap([]Limb, &tmp1, &tmp2); |
| 2634 | 3322 | } |
| 2635 | 3323 | } |
Andrew Kelley <andrew@ziglang.org> 2021-10-04 14:29:09-04:00
GitHub <noreply@github.com> 2021-10-04 14:29:09-04:00
Signed by PGP key 4AEE18F83AFDEB23