mod.rs (77498B)
1 #![allow(dead_code)] 2 3 use std::cmp::Ordering; 4 5 use num_bigint::{BigInt, Sign as BigSign}; 6 use num_traits::Signed; 7 8 type LimbVec = Vec<u64>; 9 10 #[derive(Clone, Copy, Debug, Eq, PartialEq)] 11 enum Sign { 12 Negative, 13 Zero, 14 Positive, 15 } 16 17 #[derive(Clone, Debug, Eq, PartialEq)] 18 pub(super) struct LimbInt { 19 sign: Sign, 20 limbs: LimbVec, 21 } 22 23 impl Default for LimbInt { 24 fn default() -> Self { 25 Self::zero() 26 } 27 } 28 29 pub(super) struct ExtendedGcd { 30 pub(super) x: LimbInt, 31 pub(super) y: LimbInt, 32 pub(super) gcd: LimbInt, 33 } 34 35 pub(super) struct LeftExtendedGcd { 36 pub(super) x: LimbInt, 37 pub(super) gcd: LimbInt, 38 } 39 40 #[derive(Default)] 41 pub(super) struct LimbScratch { 42 division: DivisionScratch, 43 linear_left: LimbVec, 44 linear_right: LimbVec, 45 } 46 47 #[derive(Default)] 48 struct DivisionScratch { 49 normalized_numerator: LimbVec, 50 normalized_denominator: LimbVec, 51 quotient: LimbVec, 52 remainder: LimbVec, 53 } 54 55 impl LimbInt { 56 pub(super) fn zero() -> Self { 57 Self { 58 sign: Sign::Zero, 59 limbs: LimbVec::new(), 60 } 61 } 62 63 pub(super) fn one() -> Self { 64 Self { 65 sign: Sign::Positive, 66 limbs: vec![1], 67 } 68 } 69 70 pub(super) fn from_u64(value: u64) -> Self { 71 if value == 0 { 72 Self::zero() 73 } else { 74 Self { 75 sign: Sign::Positive, 76 limbs: vec![value], 77 } 78 } 79 } 80 81 pub(super) fn from_i128(value: i128) -> Self { 82 if value == 0 { 83 return Self::zero(); 84 } 85 86 let sign = if value < 0 { 87 Sign::Negative 88 } else { 89 Sign::Positive 90 }; 91 let mut magnitude = value.unsigned_abs(); 92 let mut limbs = LimbVec::new(); 93 while magnitude > 0 { 94 limbs.push(magnitude as u64); 95 magnitude >>= 64; 96 } 97 Self { sign, limbs } 98 } 99 100 pub(super) fn from_bigint(value: &BigInt) -> Self { 101 if value == &BigInt::from(0) { 102 return Self::zero(); 103 } 104 105 let sign = if value.is_negative() { 106 Sign::Negative 107 } else { 108 Sign::Positive 109 }; 110 let magnitude = value 111 .abs() 112 .to_biguint() 113 .expect("absolute BigInt is non-negative"); 114 let mut limbs: LimbVec = magnitude.iter_u64_digits().collect(); 115 trim_leading_zero_limbs(&mut limbs); 116 Self { sign, limbs } 117 } 118 119 pub(super) fn to_bigint(&self) -> BigInt { 120 let sign = match self.sign { 121 Sign::Negative => BigSign::Minus, 122 Sign::Zero => BigSign::NoSign, 123 Sign::Positive => BigSign::Plus, 124 }; 125 BigInt::from_biguint(sign, self.abs_biguint()) 126 } 127 128 #[inline(always)] 129 pub(super) fn is_zero(&self) -> bool { 130 self.sign == Sign::Zero 131 } 132 133 #[inline(always)] 134 pub(super) fn is_negative(&self) -> bool { 135 self.sign == Sign::Negative 136 } 137 138 #[inline(always)] 139 pub(super) fn is_one(&self) -> bool { 140 self.sign == Sign::Positive && self.limbs.as_slice() == [1] 141 } 142 143 #[inline(always)] 144 pub(super) fn is_minus_one(&self) -> bool { 145 self.sign == Sign::Negative && self.limbs.as_slice() == [1] 146 } 147 148 #[inline(always)] 149 pub(super) fn bit_len(&self) -> u64 { 150 let Some(last) = self.limbs.last() else { 151 return 0; 152 }; 153 let top_bits = u64::BITS - last.leading_zeros(); 154 ((self.limbs.len() as u64 - 1) * 64) + u64::from(top_bits) 155 } 156 157 #[inline(always)] 158 pub(super) fn abs_cmp(&self, other: &Self) -> Ordering { 159 cmp_abs_limbs(&self.limbs, &other.limbs) 160 } 161 162 #[inline(always)] 163 pub(super) fn abs_cmp_double(&self, other: &Self) -> Ordering { 164 cmp_abs_to_double(&self.limbs, &other.limbs) 165 } 166 167 pub(super) fn cmp(&self, other: &Self) -> Ordering { 168 match (self.sign, other.sign) { 169 (Sign::Negative, Sign::Negative) => cmp_abs_limbs(&other.limbs, &self.limbs), 170 (Sign::Negative, _) => Ordering::Less, 171 (_, Sign::Negative) => Ordering::Greater, 172 (Sign::Zero, Sign::Zero) => Ordering::Equal, 173 (Sign::Zero, Sign::Positive) => Ordering::Less, 174 (Sign::Positive, Sign::Zero) => Ordering::Greater, 175 (Sign::Positive, Sign::Positive) => cmp_abs_limbs(&self.limbs, &other.limbs), 176 } 177 } 178 179 pub(super) fn negated(mut self) -> Self { 180 self.sign = match self.sign { 181 Sign::Negative => Sign::Positive, 182 Sign::Zero => Sign::Zero, 183 Sign::Positive => Sign::Negative, 184 }; 185 self 186 } 187 188 pub(super) fn negate_assign(&mut self) { 189 self.sign = match self.sign { 190 Sign::Negative => Sign::Positive, 191 Sign::Zero => Sign::Zero, 192 Sign::Positive => Sign::Negative, 193 }; 194 } 195 196 pub(super) fn add(&self, other: &Self) -> Self { 197 match (self.sign, other.sign) { 198 (Sign::Zero, _) => other.clone(), 199 (_, Sign::Zero) => self.clone(), 200 (Sign::Positive, Sign::Positive) => { 201 Self::from_parts(Sign::Positive, add_abs_limbs(&self.limbs, &other.limbs)) 202 } 203 (Sign::Negative, Sign::Negative) => { 204 Self::from_parts(Sign::Negative, add_abs_limbs(&self.limbs, &other.limbs)) 205 } 206 (Sign::Positive, Sign::Negative) => { 207 subtract_signed_abs(&self.limbs, Sign::Positive, &other.limbs, Sign::Negative) 208 } 209 (Sign::Negative, Sign::Positive) => { 210 subtract_signed_abs(&other.limbs, Sign::Positive, &self.limbs, Sign::Negative) 211 } 212 } 213 } 214 215 pub(super) fn add_owned(mut self, other: &Self) -> Self { 216 self.add_assign(other); 217 self 218 } 219 220 #[inline(always)] 221 pub(super) fn add_into(&self, other: &Self, output: &mut Self) { 222 match (self.sign, other.sign) { 223 (Sign::Zero, _) => { 224 *output = other.clone(); 225 } 226 (_, Sign::Zero) => { 227 *output = self.clone(); 228 } 229 (Sign::Positive, Sign::Positive) | (Sign::Negative, Sign::Negative) => { 230 add_abs_limbs_into(&self.limbs, &other.limbs, &mut output.limbs); 231 output.sign = if output.limbs.is_empty() { 232 Sign::Zero 233 } else { 234 self.sign 235 }; 236 } 237 (Sign::Positive, Sign::Negative) | (Sign::Negative, Sign::Positive) => { 238 let sign = combine_signed_abs_limbs_into( 239 self.sign, 240 &self.limbs, 241 other.sign, 242 &other.limbs, 243 &mut output.limbs, 244 ); 245 output.sign = if output.limbs.is_empty() { 246 Sign::Zero 247 } else { 248 sign 249 }; 250 } 251 } 252 } 253 254 pub(super) fn add_assign(&mut self, other: &Self) { 255 match (self.sign, other.sign) { 256 (Sign::Zero, _) => { 257 *self = other.clone(); 258 } 259 (_, Sign::Zero) => {} 260 (Sign::Positive, Sign::Positive) | (Sign::Negative, Sign::Negative) => { 261 add_abs_limbs_assign(&mut self.limbs, &other.limbs); 262 } 263 (Sign::Positive, Sign::Negative) => match cmp_abs_limbs(&self.limbs, &other.limbs) { 264 Ordering::Greater => sub_abs_limbs_assign(&mut self.limbs, &other.limbs), 265 Ordering::Less => { 266 self.limbs = sub_abs_limbs(&other.limbs, &self.limbs); 267 self.sign = Sign::Negative; 268 } 269 Ordering::Equal => *self = Self::zero(), 270 }, 271 (Sign::Negative, Sign::Positive) => match cmp_abs_limbs(&self.limbs, &other.limbs) { 272 Ordering::Greater => sub_abs_limbs_assign(&mut self.limbs, &other.limbs), 273 Ordering::Less => { 274 self.limbs = sub_abs_limbs(&other.limbs, &self.limbs); 275 self.sign = Sign::Positive; 276 } 277 Ordering::Equal => *self = Self::zero(), 278 }, 279 } 280 } 281 282 pub(super) fn sub(&self, other: &Self) -> Self { 283 let mut result = self.clone(); 284 result.sub_assign(other); 285 result 286 } 287 288 pub(super) fn sub_owned(mut self, other: &Self) -> Self { 289 self.sub_assign(other); 290 self 291 } 292 293 #[inline(always)] 294 pub(super) fn sub_into(&self, other: &Self, output: &mut Self) { 295 match (self.sign, other.sign) { 296 (_, Sign::Zero) => { 297 *output = self.clone(); 298 } 299 (Sign::Zero, _) => { 300 *output = other.clone().negated(); 301 } 302 (Sign::Positive, Sign::Positive) | (Sign::Negative, Sign::Negative) => { 303 match cmp_abs_limbs(&self.limbs, &other.limbs) { 304 Ordering::Greater => { 305 sub_abs_limbs_into(&self.limbs, &other.limbs, &mut output.limbs); 306 output.sign = self.sign; 307 } 308 Ordering::Less => { 309 sub_abs_limbs_into(&other.limbs, &self.limbs, &mut output.limbs); 310 output.sign = match self.sign { 311 Sign::Positive => Sign::Negative, 312 Sign::Negative => Sign::Positive, 313 Sign::Zero => Sign::Zero, 314 }; 315 } 316 Ordering::Equal => { 317 *output = Self::zero(); 318 } 319 } 320 } 321 (Sign::Positive, Sign::Negative) | (Sign::Negative, Sign::Positive) => { 322 add_abs_limbs_into(&self.limbs, &other.limbs, &mut output.limbs); 323 output.sign = if output.limbs.is_empty() { 324 Sign::Zero 325 } else { 326 self.sign 327 }; 328 } 329 } 330 } 331 332 pub(super) fn sub_assign(&mut self, other: &Self) { 333 match (self.sign, other.sign) { 334 (_, Sign::Zero) => {} 335 (Sign::Zero, _) => { 336 *self = other.clone().negated(); 337 } 338 (Sign::Positive, Sign::Positive) => match cmp_abs_limbs(&self.limbs, &other.limbs) { 339 Ordering::Greater => sub_abs_limbs_assign(&mut self.limbs, &other.limbs), 340 Ordering::Less => { 341 self.limbs = sub_abs_limbs(&other.limbs, &self.limbs); 342 self.sign = Sign::Negative; 343 } 344 Ordering::Equal => *self = Self::zero(), 345 }, 346 (Sign::Negative, Sign::Negative) => match cmp_abs_limbs(&self.limbs, &other.limbs) { 347 Ordering::Greater => sub_abs_limbs_assign(&mut self.limbs, &other.limbs), 348 Ordering::Less => { 349 self.limbs = sub_abs_limbs(&other.limbs, &self.limbs); 350 self.sign = Sign::Positive; 351 } 352 Ordering::Equal => *self = Self::zero(), 353 }, 354 (Sign::Positive, Sign::Negative) | (Sign::Negative, Sign::Positive) => { 355 add_abs_limbs_assign(&mut self.limbs, &other.limbs); 356 } 357 } 358 } 359 360 pub(super) fn mul(&self, other: &Self) -> Self { 361 if self.is_zero() || other.is_zero() { 362 return Self::zero(); 363 } 364 365 let sign = if self.sign == other.sign { 366 Sign::Positive 367 } else { 368 Sign::Negative 369 }; 370 let limbs = match (self.limbs.as_slice(), other.limbs.as_slice()) { 371 ([scalar], value) => mul_abs_one_limb(value, *scalar), 372 (value, [scalar]) => mul_abs_one_limb(value, *scalar), 373 _ => mul_abs_limbs(&self.limbs, &other.limbs), 374 }; 375 Self::from_parts(sign, limbs) 376 } 377 378 pub(super) fn mul_into(&self, other: &Self, output: &mut Self) { 379 if self.is_zero() || other.is_zero() { 380 *output = Self::zero(); 381 return; 382 } 383 384 output.sign = if self.sign == other.sign { 385 Sign::Positive 386 } else { 387 Sign::Negative 388 }; 389 match (self.limbs.as_slice(), other.limbs.as_slice()) { 390 ([scalar], value) => mul_abs_one_limb_into(value, *scalar, &mut output.limbs), 391 (value, [scalar]) => mul_abs_one_limb_into(value, *scalar, &mut output.limbs), 392 _ => mul_abs_limbs_into(&self.limbs, &other.limbs, &mut output.limbs), 393 } 394 if output.limbs.is_empty() { 395 output.sign = Sign::Zero; 396 } 397 } 398 399 pub(super) fn mul_i128(&self, scalar: i128) -> Self { 400 if self.is_zero() || scalar == 0 { 401 return Self::zero(); 402 } 403 if scalar == 1 { 404 return self.clone(); 405 } 406 if scalar == -1 { 407 return self.clone().negated(); 408 } 409 410 let sign = if scalar < 0 { 411 match self.sign { 412 Sign::Negative => Sign::Positive, 413 Sign::Positive => Sign::Negative, 414 Sign::Zero => Sign::Zero, 415 } 416 } else { 417 self.sign 418 }; 419 Self::from_parts(sign, mul_abs_small(&self.limbs, scalar.unsigned_abs())) 420 } 421 422 pub(super) fn linear_combination_i128_with_scratch( 423 left: &Self, 424 left_scalar: i128, 425 right: &Self, 426 right_scalar: i128, 427 scratch: &mut LimbScratch, 428 ) -> Self { 429 let mut output = Self::zero(); 430 Self::linear_combination_i128_into( 431 left, 432 left_scalar, 433 right, 434 right_scalar, 435 &mut output, 436 scratch, 437 ); 438 output 439 } 440 441 #[inline(always)] 442 fn linear_combination_i128_into( 443 left: &Self, 444 left_scalar: i128, 445 right: &Self, 446 right_scalar: i128, 447 output: &mut Self, 448 scratch: &mut LimbScratch, 449 ) { 450 if left_scalar == 0 && right_scalar == 0 { 451 *output = Self::zero(); 452 return; 453 } 454 if left_scalar == 0 { 455 right.mul_i128_into(right_scalar, output); 456 return; 457 } 458 if right_scalar == 0 { 459 left.mul_i128_into(left_scalar, output); 460 return; 461 } 462 463 let left_abs = left_scalar.unsigned_abs(); 464 let right_abs = right_scalar.unsigned_abs(); 465 let left_sign = signed_scalar_sign(left.sign, left_scalar); 466 let right_sign = signed_scalar_sign(right.sign, right_scalar); 467 468 if left_abs == 1 && right_abs == 1 { 469 let sign = combine_signed_abs_limbs_into( 470 left_sign, 471 &left.limbs, 472 right_sign, 473 &right.limbs, 474 &mut output.limbs, 475 ); 476 output.sign = if output.limbs.is_empty() { 477 Sign::Zero 478 } else { 479 sign 480 }; 481 return; 482 } 483 484 if left_abs == 1 { 485 mul_abs_small_into(&right.limbs, right_abs, &mut scratch.linear_right); 486 let sign = combine_signed_abs_limbs_into( 487 left_sign, 488 &left.limbs, 489 right_sign, 490 &scratch.linear_right, 491 &mut output.limbs, 492 ); 493 output.sign = if output.limbs.is_empty() { 494 Sign::Zero 495 } else { 496 sign 497 }; 498 return; 499 } 500 501 if right_abs == 1 { 502 mul_abs_small_into(&left.limbs, left_abs, &mut scratch.linear_left); 503 let sign = combine_signed_abs_limbs_into( 504 left_sign, 505 &scratch.linear_left, 506 right_sign, 507 &right.limbs, 508 &mut output.limbs, 509 ); 510 output.sign = if output.limbs.is_empty() { 511 Sign::Zero 512 } else { 513 sign 514 }; 515 return; 516 } 517 518 mul_abs_small_into(&left.limbs, left_abs, &mut scratch.linear_left); 519 mul_abs_small_into(&right.limbs, right_abs, &mut scratch.linear_right); 520 521 let sign = combine_signed_abs_limbs_into( 522 left_sign, 523 &scratch.linear_left, 524 right_sign, 525 &scratch.linear_right, 526 &mut output.limbs, 527 ); 528 529 output.sign = if output.limbs.is_empty() { 530 Sign::Zero 531 } else { 532 sign 533 }; 534 } 535 536 #[inline(always)] 537 fn mul_i128_into(&self, scalar: i128, output: &mut Self) { 538 if self.is_zero() || scalar == 0 { 539 *output = Self::zero(); 540 return; 541 } 542 543 let sign = if scalar < 0 { 544 match self.sign { 545 Sign::Negative => Sign::Positive, 546 Sign::Positive => Sign::Negative, 547 Sign::Zero => Sign::Zero, 548 } 549 } else { 550 self.sign 551 }; 552 output.sign = sign; 553 output.limbs.clear(); 554 mul_abs_small_into(&self.limbs, scalar.unsigned_abs(), &mut output.limbs); 555 if output.limbs.is_empty() { 556 output.sign = Sign::Zero; 557 } 558 } 559 560 pub(super) fn square(&self) -> Self { 561 if self.is_zero() { 562 return Self::zero(); 563 } 564 Self::from_parts(Sign::Positive, square_abs_limbs(&self.limbs)) 565 } 566 567 pub(super) fn square_into(&self, output: &mut Self) { 568 if self.is_zero() { 569 *output = Self::zero(); 570 return; 571 } 572 output.sign = Sign::Positive; 573 mul_abs_limbs_into(&self.limbs, &self.limbs, &mut output.limbs); 574 if output.limbs.is_empty() { 575 output.sign = Sign::Zero; 576 } 577 } 578 579 pub(super) fn shl_bits(&self, bits: usize) -> Self { 580 if self.is_zero() || bits == 0 { 581 return self.clone(); 582 } 583 Self::from_parts(self.sign, shl_abs_limbs(&self.limbs, bits)) 584 } 585 586 pub(super) fn shl_bits_into(&self, bits: usize, output: &mut Self) { 587 if self.is_zero() { 588 *output = Self::zero(); 589 return; 590 } 591 if bits == 0 { 592 *output = self.clone(); 593 return; 594 } 595 output.sign = self.sign; 596 shl_abs_limbs_into(&self.limbs, bits, &mut output.limbs); 597 if output.limbs.is_empty() { 598 output.sign = Sign::Zero; 599 } 600 } 601 602 pub(super) fn shl_bits_owned(mut self, bits: usize) -> Self { 603 if self.is_zero() || bits == 0 { 604 return self; 605 } 606 shl_abs_limbs_assign(&mut self.limbs, bits); 607 self 608 } 609 610 pub(super) fn shl_bits_assign(&mut self, bits: usize) { 611 if self.is_zero() || bits == 0 { 612 return; 613 } 614 shl_abs_limbs_assign(&mut self.limbs, bits); 615 } 616 617 pub(super) fn shr_abs_bits(&self, bits: usize) -> Self { 618 if self.is_zero() || bits == 0 { 619 return self.abs(); 620 } 621 Self::from_parts(Sign::Positive, shr_abs_limbs(&self.limbs, bits)) 622 } 623 624 pub(super) fn div_rem(&self, divisor: &Self) -> (Self, Self) { 625 let mut scratch = LimbScratch::default(); 626 self.div_rem_with_scratch(divisor, &mut scratch) 627 } 628 629 pub(super) fn div_rem_with_scratch( 630 &self, 631 divisor: &Self, 632 scratch: &mut LimbScratch, 633 ) -> (Self, Self) { 634 assert!(!divisor.is_zero(), "division by zero"); 635 if self.is_zero() { 636 return (Self::zero(), Self::zero()); 637 } 638 639 let (quotient_limbs, remainder_limbs) = 640 div_rem_abs_limbs_with_scratch(&self.limbs, &divisor.limbs, &mut scratch.division); 641 let quotient_sign = if quotient_limbs.is_empty() { 642 Sign::Zero 643 } else if self.sign == divisor.sign { 644 Sign::Positive 645 } else { 646 Sign::Negative 647 }; 648 let remainder_sign = if remainder_limbs.is_empty() { 649 Sign::Zero 650 } else { 651 self.sign 652 }; 653 654 ( 655 Self::from_parts(quotient_sign, quotient_limbs), 656 Self::from_parts(remainder_sign, remainder_limbs), 657 ) 658 } 659 660 pub(super) fn rem(&self, divisor: &Self) -> Self { 661 self.div_rem(divisor).1 662 } 663 664 pub(super) fn rem_with_scratch(&self, divisor: &Self, scratch: &mut LimbScratch) -> Self { 665 assert!(!divisor.is_zero(), "division by zero"); 666 if self.is_zero() { 667 return Self::zero(); 668 } 669 670 let remainder_limbs = 671 rem_abs_limbs_with_scratch(&self.limbs, &divisor.limbs, &mut scratch.division); 672 let remainder_sign = if remainder_limbs.is_empty() { 673 Sign::Zero 674 } else { 675 self.sign 676 }; 677 Self::from_parts(remainder_sign, remainder_limbs) 678 } 679 680 pub(super) fn div(&self, divisor: &Self) -> Self { 681 self.div_rem(divisor).0 682 } 683 684 pub(super) fn div_with_scratch(&self, divisor: &Self, scratch: &mut LimbScratch) -> Self { 685 assert!(!divisor.is_zero(), "division by zero"); 686 if self.is_zero() { 687 return Self::zero(); 688 } 689 match self.abs_cmp(divisor) { 690 Ordering::Less => return Self::zero(), 691 Ordering::Equal => { 692 return if self.sign == divisor.sign { 693 Self::one() 694 } else { 695 Self::from_i128(-1) 696 }; 697 } 698 Ordering::Greater => { 699 if cmp_abs_to_double(&self.limbs, &divisor.limbs) == Ordering::Less { 700 return if self.sign == divisor.sign { 701 Self::one() 702 } else { 703 Self::from_i128(-1) 704 }; 705 } 706 } 707 } 708 let quotient_limbs = 709 div_abs_limbs_with_scratch(&self.limbs, &divisor.limbs, &mut scratch.division); 710 let quotient_sign = if quotient_limbs.is_empty() { 711 Sign::Zero 712 } else if self.sign == divisor.sign { 713 Sign::Positive 714 } else { 715 Sign::Negative 716 }; 717 Self::from_parts(quotient_sign, quotient_limbs) 718 } 719 720 pub(super) fn div2_exact(&self) -> Self { 721 debug_assert!(self.limbs.first().copied().unwrap_or(0) & 1 == 0); 722 if self.is_zero() { 723 return Self::zero(); 724 } 725 Self::from_parts(self.sign, shr_abs_limbs(&self.limbs, 1)) 726 } 727 728 pub(super) fn div_floor(&self, divisor: &Self) -> Self { 729 let (quotient, remainder) = self.div_rem(divisor); 730 self.finish_div_floor(divisor, quotient, remainder) 731 } 732 733 pub(super) fn div_floor_with_scratch(&self, divisor: &Self, scratch: &mut LimbScratch) -> Self { 734 assert!(!divisor.is_zero(), "division by zero"); 735 if self.is_zero() { 736 return Self::zero(); 737 } 738 739 match self.abs_cmp(divisor) { 740 Ordering::Less => { 741 return if self.sign == divisor.sign { 742 Self::zero() 743 } else { 744 Self::from_i128(-1) 745 }; 746 } 747 Ordering::Equal => { 748 return if self.sign == divisor.sign { 749 Self::one() 750 } else { 751 Self::from_i128(-1) 752 }; 753 } 754 Ordering::Greater => { 755 if cmp_abs_to_double(&self.limbs, &divisor.limbs) == Ordering::Less { 756 return if self.sign == divisor.sign { 757 Self::one() 758 } else { 759 Self::from_i128(-2) 760 }; 761 } 762 } 763 } 764 765 if self.sign == divisor.sign { 766 let quotient_limbs = 767 div_abs_limbs_with_scratch(&self.limbs, &divisor.limbs, &mut scratch.division); 768 let quotient_sign = if quotient_limbs.is_empty() { 769 Sign::Zero 770 } else { 771 Sign::Positive 772 }; 773 return Self::from_parts(quotient_sign, quotient_limbs); 774 } 775 776 let (quotient, remainder) = self.div_rem_with_scratch(divisor, scratch); 777 self.finish_div_floor(divisor, quotient, remainder) 778 } 779 780 fn finish_div_floor(&self, divisor: &Self, quotient: Self, remainder: Self) -> Self { 781 if remainder.is_zero() || self.sign == divisor.sign { 782 quotient 783 } else { 784 quotient.sub(&Self::one()) 785 } 786 } 787 788 pub(super) fn mod_positive(&self, modulus: &Self) -> Self { 789 let mut scratch = LimbScratch::default(); 790 self.mod_positive_with_scratch(modulus, &mut scratch) 791 } 792 793 pub(super) fn mod_positive_with_scratch( 794 &self, 795 modulus: &Self, 796 scratch: &mut LimbScratch, 797 ) -> Self { 798 assert!(!modulus.is_zero(), "division by zero"); 799 if self.is_zero() { 800 return Self::zero(); 801 } 802 match self.abs_cmp(modulus) { 803 Ordering::Less => { 804 return if self.is_negative() { 805 self.add(modulus) 806 } else { 807 self.clone() 808 }; 809 } 810 Ordering::Equal | Ordering::Greater => {} 811 } 812 813 let remainder_limbs = 814 rem_abs_limbs_with_scratch(&self.limbs, &modulus.limbs, &mut scratch.division); 815 let remainder_sign = if remainder_limbs.is_empty() { 816 Sign::Zero 817 } else { 818 self.sign 819 }; 820 let mut remainder = Self::from_parts(remainder_sign, remainder_limbs); 821 if remainder.is_negative() { 822 remainder = remainder.add(modulus); 823 } 824 remainder 825 } 826 827 pub(super) fn abs(&self) -> Self { 828 if self.is_zero() { 829 Self::zero() 830 } else { 831 Self { 832 sign: Sign::Positive, 833 limbs: self.limbs.clone(), 834 } 835 } 836 } 837 838 #[inline(always)] 839 pub(super) fn shifted_low_word(&self, shift_bits: u64) -> u64 { 840 debug_assert!(!self.is_negative()); 841 let digit_index = usize::try_from(shift_bits / 64).unwrap_or(usize::MAX); 842 let offset = (shift_bits % 64) as u32; 843 let low = self.limbs.get(digit_index).copied().unwrap_or(0); 844 if offset == 0 { 845 return low; 846 } 847 let high = self.limbs.get(digit_index + 1).copied().unwrap_or(0); 848 (low >> offset) | (high << (64 - offset)) 849 } 850 851 pub(super) fn extended_gcd(&self, other: &Self) -> ExtendedGcd { 852 let mut scratch = LimbScratch::default(); 853 self.extended_gcd_with_scratch(other, &mut scratch) 854 } 855 856 pub(super) fn extended_gcd_with_scratch( 857 &self, 858 other: &Self, 859 scratch: &mut LimbScratch, 860 ) -> ExtendedGcd { 861 let mut old_r = self.clone(); 862 let mut r = other.clone(); 863 let mut old_s = Self::one(); 864 let mut s = Self::zero(); 865 let mut old_t = Self::zero(); 866 let mut t = Self::one(); 867 let mut next_old_r = Self::zero(); 868 let mut next_r = Self::zero(); 869 let mut next_old_s = Self::zero(); 870 let mut next_s = Self::zero(); 871 let mut next_old_t = Self::zero(); 872 let mut next_t = Self::zero(); 873 874 while !r.is_zero() { 875 if old_r.is_one() { 876 break; 877 } 878 if r.is_one() { 879 old_r = r; 880 old_s = s; 881 old_t = t; 882 break; 883 } 884 885 if !old_r.is_negative() && !r.is_negative() && old_r.abs_cmp(&r) == Ordering::Less { 886 std::mem::swap(&mut old_r, &mut r); 887 std::mem::swap(&mut old_s, &mut s); 888 std::mem::swap(&mut old_t, &mut t); 889 continue; 890 } 891 892 if !old_r.is_negative() && !r.is_negative() { 893 let bits = old_r.bit_len().saturating_sub(63); 894 let mut rr2 = old_r.shifted_low_word(bits); 895 let mut rr1 = r.shifted_low_word(bits); 896 897 let mut aa2 = 0_i128; 898 let mut aa1 = 1_i128; 899 let mut bb2 = 1_i128; 900 let mut bb1 = 0_i128; 901 let mut steps = 0_u32; 902 903 while rr1 != 0 { 904 let q = rr2 / rr1; 905 if q == 0 { 906 break; 907 } 908 let next_r = rr2 - q * rr1; 909 let q = i128::from(q); 910 let next_a = aa2 - q * aa1; 911 let next_b = bb2 - q * bb1; 912 let next_r_i = i128::from(next_r); 913 let rr1_minus_next_r = i128::from(rr1 - next_r); 914 915 if steps & 1 == 1 { 916 if next_r_i < -next_b || rr1_minus_next_r < next_a - aa1 { 917 break; 918 } 919 } else if next_r_i < -next_a || rr1_minus_next_r < next_b - bb1 { 920 break; 921 } 922 923 rr2 = rr1; 924 rr1 = next_r; 925 aa2 = aa1; 926 aa1 = next_a; 927 bb2 = bb1; 928 bb1 = next_b; 929 steps += 1; 930 } 931 932 if steps != 0 { 933 Self::linear_combination_i128_into( 934 &old_r, 935 bb2, 936 &r, 937 aa2, 938 &mut next_old_r, 939 scratch, 940 ); 941 Self::linear_combination_i128_into(&r, aa1, &old_r, bb1, &mut next_r, scratch); 942 Self::linear_combination_i128_into( 943 &old_s, 944 bb2, 945 &s, 946 aa2, 947 &mut next_old_s, 948 scratch, 949 ); 950 Self::linear_combination_i128_into(&s, aa1, &old_s, bb1, &mut next_s, scratch); 951 Self::linear_combination_i128_into( 952 &old_t, 953 bb2, 954 &t, 955 aa2, 956 &mut next_old_t, 957 scratch, 958 ); 959 Self::linear_combination_i128_into(&t, aa1, &old_t, bb1, &mut next_t, scratch); 960 961 std::mem::swap(&mut old_r, &mut next_old_r); 962 std::mem::swap(&mut r, &mut next_r); 963 std::mem::swap(&mut old_s, &mut next_old_s); 964 std::mem::swap(&mut s, &mut next_s); 965 std::mem::swap(&mut old_t, &mut next_old_t); 966 std::mem::swap(&mut t, &mut next_t); 967 968 if old_r.is_negative() { 969 old_r = old_r.negated(); 970 old_s = old_s.negated(); 971 old_t = old_t.negated(); 972 } 973 if r.is_negative() { 974 r = r.negated(); 975 s = s.negated(); 976 t = t.negated(); 977 } 978 continue; 979 } 980 } 981 982 let (quotient, next_r) = old_r.div_rem_with_scratch(&r, scratch); 983 old_r = r; 984 r = next_r; 985 986 quotient.mul_into(&s, &mut next_s); 987 let mut fallback_next_s = std::mem::take(&mut old_s); 988 fallback_next_s.sub_assign(&next_s); 989 old_s = s; 990 s = fallback_next_s; 991 992 quotient.mul_into(&t, &mut next_t); 993 let mut fallback_next_t = std::mem::take(&mut old_t); 994 fallback_next_t.sub_assign(&next_t); 995 old_t = t; 996 t = fallback_next_t; 997 } 998 999 if old_r.is_negative() { 1000 old_r = old_r.negated(); 1001 old_s = old_s.negated(); 1002 old_t = old_t.negated(); 1003 } 1004 1005 ExtendedGcd { 1006 x: old_s, 1007 y: old_t, 1008 gcd: old_r, 1009 } 1010 } 1011 1012 pub(super) fn left_extended_gcd_with_scratch( 1013 &self, 1014 other: &Self, 1015 scratch: &mut LimbScratch, 1016 ) -> LeftExtendedGcd { 1017 let mut old_r = self.clone(); 1018 let mut r = other.clone(); 1019 let mut old_s = Self::one(); 1020 let mut s = Self::zero(); 1021 let mut next_old_r = Self::zero(); 1022 let mut next_r = Self::zero(); 1023 let mut next_old_s = Self::zero(); 1024 let mut next_s = Self::zero(); 1025 1026 while !r.is_zero() { 1027 if old_r.is_one() { 1028 break; 1029 } 1030 if r.is_one() { 1031 old_r = r; 1032 old_s = s; 1033 break; 1034 } 1035 1036 if !old_r.is_negative() && !r.is_negative() && old_r.abs_cmp(&r) == Ordering::Less { 1037 std::mem::swap(&mut old_r, &mut r); 1038 std::mem::swap(&mut old_s, &mut s); 1039 continue; 1040 } 1041 1042 if !old_r.is_negative() && !r.is_negative() { 1043 let bits = old_r.bit_len().saturating_sub(63); 1044 let mut rr2 = old_r.shifted_low_word(bits); 1045 let mut rr1 = r.shifted_low_word(bits); 1046 1047 let mut aa2 = 0_i128; 1048 let mut aa1 = 1_i128; 1049 let mut bb2 = 1_i128; 1050 let mut bb1 = 0_i128; 1051 let mut steps = 0_u32; 1052 1053 while rr1 != 0 { 1054 let q = rr2 / rr1; 1055 if q == 0 { 1056 break; 1057 } 1058 let next_r = rr2 - q * rr1; 1059 let q = i128::from(q); 1060 let next_a = aa2 - q * aa1; 1061 let next_b = bb2 - q * bb1; 1062 let next_r_i = i128::from(next_r); 1063 let rr1_minus_next_r = i128::from(rr1 - next_r); 1064 1065 if steps & 1 == 1 { 1066 if next_r_i < -next_b || rr1_minus_next_r < next_a - aa1 { 1067 break; 1068 } 1069 } else if next_r_i < -next_a || rr1_minus_next_r < next_b - bb1 { 1070 break; 1071 } 1072 1073 rr2 = rr1; 1074 rr1 = next_r; 1075 aa2 = aa1; 1076 aa1 = next_a; 1077 bb2 = bb1; 1078 bb1 = next_b; 1079 steps += 1; 1080 } 1081 1082 if steps != 0 { 1083 Self::linear_combination_i128_into( 1084 &old_r, 1085 bb2, 1086 &r, 1087 aa2, 1088 &mut next_old_r, 1089 scratch, 1090 ); 1091 Self::linear_combination_i128_into(&r, aa1, &old_r, bb1, &mut next_r, scratch); 1092 Self::linear_combination_i128_into( 1093 &old_s, 1094 bb2, 1095 &s, 1096 aa2, 1097 &mut next_old_s, 1098 scratch, 1099 ); 1100 Self::linear_combination_i128_into(&s, aa1, &old_s, bb1, &mut next_s, scratch); 1101 1102 std::mem::swap(&mut old_r, &mut next_old_r); 1103 std::mem::swap(&mut r, &mut next_r); 1104 std::mem::swap(&mut old_s, &mut next_old_s); 1105 std::mem::swap(&mut s, &mut next_s); 1106 1107 if old_r.is_negative() { 1108 old_r = old_r.negated(); 1109 old_s = old_s.negated(); 1110 } 1111 if r.is_negative() { 1112 r = r.negated(); 1113 s = s.negated(); 1114 } 1115 continue; 1116 } 1117 } 1118 1119 let (quotient, next_r) = old_r.div_rem_with_scratch(&r, scratch); 1120 old_r = r; 1121 r = next_r; 1122 1123 quotient.mul_into(&s, &mut next_s); 1124 let mut fallback_next_s = std::mem::take(&mut old_s); 1125 fallback_next_s.sub_assign(&next_s); 1126 old_s = s; 1127 s = fallback_next_s; 1128 } 1129 1130 if old_r.is_negative() { 1131 old_r = old_r.negated(); 1132 old_s = old_s.negated(); 1133 } 1134 1135 LeftExtendedGcd { 1136 x: old_s, 1137 gcd: old_r, 1138 } 1139 } 1140 1141 pub(super) fn left_extended_gcd_positive_with_scratch( 1142 &self, 1143 other: &Self, 1144 scratch: &mut LimbScratch, 1145 ) -> LeftExtendedGcd { 1146 debug_assert!(!self.is_negative()); 1147 debug_assert!(!other.is_negative()); 1148 1149 let mut old_r; 1150 let mut r; 1151 let mut old_s; 1152 let mut s; 1153 if self.abs_cmp(other) == Ordering::Less { 1154 old_r = other.clone(); 1155 r = self.clone(); 1156 old_s = Self::zero(); 1157 s = Self::one(); 1158 } else { 1159 old_r = self.clone(); 1160 r = other.clone(); 1161 old_s = Self::one(); 1162 s = Self::zero(); 1163 } 1164 let mut next_old_r = Self::zero(); 1165 let mut next_r = Self::zero(); 1166 let mut next_old_s = Self::zero(); 1167 let mut next_s = Self::zero(); 1168 1169 while !r.is_zero() { 1170 if old_r.is_one() { 1171 break; 1172 } 1173 if r.is_one() { 1174 old_r = r; 1175 old_s = s; 1176 break; 1177 } 1178 1179 let bits = old_r.bit_len().saturating_sub(63); 1180 let mut rr2 = old_r.shifted_low_word(bits); 1181 let mut rr1 = r.shifted_low_word(bits); 1182 1183 let mut aa2 = 0_i128; 1184 let mut aa1 = 1_i128; 1185 let mut bb2 = 1_i128; 1186 let mut bb1 = 0_i128; 1187 let mut steps = 0_u32; 1188 1189 while rr1 != 0 { 1190 let q = rr2 / rr1; 1191 if q == 0 { 1192 break; 1193 } 1194 let next_r = rr2 - q * rr1; 1195 let q = i128::from(q); 1196 let next_a = aa2 - q * aa1; 1197 let next_b = bb2 - q * bb1; 1198 let next_r_i = i128::from(next_r); 1199 let rr1_minus_next_r = i128::from(rr1 - next_r); 1200 1201 if steps & 1 == 1 { 1202 if next_r_i < -next_b || rr1_minus_next_r < next_a - aa1 { 1203 break; 1204 } 1205 } else if next_r_i < -next_a || rr1_minus_next_r < next_b - bb1 { 1206 break; 1207 } 1208 1209 rr2 = rr1; 1210 rr1 = next_r; 1211 aa2 = aa1; 1212 aa1 = next_a; 1213 bb2 = bb1; 1214 bb1 = next_b; 1215 steps += 1; 1216 } 1217 1218 if steps != 0 { 1219 Self::linear_combination_i128_into(&old_r, bb2, &r, aa2, &mut next_old_r, scratch); 1220 Self::linear_combination_i128_into(&r, aa1, &old_r, bb1, &mut next_r, scratch); 1221 Self::linear_combination_i128_into(&old_s, bb2, &s, aa2, &mut next_old_s, scratch); 1222 Self::linear_combination_i128_into(&s, aa1, &old_s, bb1, &mut next_s, scratch); 1223 1224 std::mem::swap(&mut old_r, &mut next_old_r); 1225 std::mem::swap(&mut r, &mut next_r); 1226 std::mem::swap(&mut old_s, &mut next_old_s); 1227 std::mem::swap(&mut s, &mut next_s); 1228 1229 if old_r.is_negative() { 1230 old_r = old_r.negated(); 1231 old_s = old_s.negated(); 1232 } 1233 if r.is_negative() { 1234 r = r.negated(); 1235 s = s.negated(); 1236 } 1237 continue; 1238 } 1239 1240 let (quotient, next_r) = old_r.div_rem_with_scratch(&r, scratch); 1241 old_r = r; 1242 r = next_r; 1243 1244 quotient.mul_into(&s, &mut next_s); 1245 let mut fallback_next_s = std::mem::take(&mut old_s); 1246 fallback_next_s.sub_assign(&next_s); 1247 old_s = s; 1248 s = fallback_next_s; 1249 } 1250 1251 if old_r.is_negative() { 1252 old_r = old_r.negated(); 1253 old_s = old_s.negated(); 1254 } 1255 1256 LeftExtendedGcd { 1257 x: old_s, 1258 gcd: old_r, 1259 } 1260 } 1261 1262 fn from_parts(sign: Sign, mut limbs: LimbVec) -> Self { 1263 trim_leading_zero_limbs(&mut limbs); 1264 let sign = if limbs.is_empty() { Sign::Zero } else { sign }; 1265 Self { sign, limbs } 1266 } 1267 1268 fn abs_biguint(&self) -> num_bigint::BigUint { 1269 let mut bytes = Vec::with_capacity(self.limbs.len() * 8); 1270 for limb in &self.limbs { 1271 bytes.extend_from_slice(&limb.to_le_bytes()); 1272 } 1273 num_bigint::BigUint::from_bytes_le(&bytes) 1274 } 1275 } 1276 1277 pub(super) fn xgcd_partial( 1278 r2: &LimbInt, 1279 r1: &LimbInt, 1280 threshold: &LimbInt, 1281 ) -> (LimbInt, LimbInt, LimbInt, LimbInt) { 1282 let mut scratch = LimbScratch::default(); 1283 xgcd_partial_with_scratch(r2, r1, threshold, &mut scratch) 1284 } 1285 1286 pub(super) fn xgcd_partial_with_scratch( 1287 r2: &LimbInt, 1288 r1: &LimbInt, 1289 threshold: &LimbInt, 1290 scratch: &mut LimbScratch, 1291 ) -> (LimbInt, LimbInt, LimbInt, LimbInt) { 1292 let mut r2 = r2.clone(); 1293 let mut r1 = r1.clone(); 1294 let mut co2 = LimbInt::zero(); 1295 let mut co1 = LimbInt::from_i128(-1); 1296 let mut next_r2 = LimbInt::zero(); 1297 let mut next_r1 = LimbInt::zero(); 1298 let mut next_co2 = LimbInt::zero(); 1299 let mut next_co1 = LimbInt::zero(); 1300 1301 while !r1.is_zero() && r1.cmp(threshold) == Ordering::Greater { 1302 let bits = r2.bit_len().saturating_sub(63); 1303 let mut rr2 = r2.shifted_low_word(bits); 1304 let mut rr1 = r1.shifted_low_word(bits); 1305 let threshold_word = threshold.shifted_low_word(bits); 1306 1307 let mut aa2 = 0_i128; 1308 let mut aa1 = 1_i128; 1309 let mut bb2 = 1_i128; 1310 let mut bb1 = 0_i128; 1311 let mut steps = 0_u32; 1312 1313 while rr1 != 0 && rr1 > threshold_word { 1314 let q = rr2 / rr1; 1315 if q == 0 { 1316 break; 1317 } 1318 let next_r = rr2 - q * rr1; 1319 let q = i128::from(q); 1320 let next_a = aa2 - q * aa1; 1321 let next_b = bb2 - q * bb1; 1322 let next_r_i = i128::from(next_r); 1323 let rr1_minus_next_r = i128::from(rr1 - next_r); 1324 1325 if steps & 1 == 1 { 1326 if next_r_i < -next_b || rr1_minus_next_r < next_a - aa1 { 1327 break; 1328 } 1329 } else if next_r_i < -next_a || rr1_minus_next_r < next_b - bb1 { 1330 break; 1331 } 1332 1333 rr2 = rr1; 1334 rr1 = next_r; 1335 aa2 = aa1; 1336 aa1 = next_a; 1337 bb2 = bb1; 1338 bb1 = next_b; 1339 steps += 1; 1340 } 1341 1342 if steps == 0 { 1343 let (q, next_r) = r2.div_rem_with_scratch(&r1, scratch); 1344 q.mul_into(&co1, &mut next_co1); 1345 let mut next_co = std::mem::take(&mut co2); 1346 next_co.sub_assign(&next_co1); 1347 r2 = r1; 1348 r1 = next_r; 1349 co2 = co1; 1350 co1 = next_co; 1351 } else { 1352 LimbInt::linear_combination_i128_into(&r2, bb2, &r1, aa2, &mut next_r2, scratch); 1353 LimbInt::linear_combination_i128_into(&r1, aa1, &r2, bb1, &mut next_r1, scratch); 1354 LimbInt::linear_combination_i128_into(&co2, bb2, &co1, aa2, &mut next_co2, scratch); 1355 LimbInt::linear_combination_i128_into(&co1, aa1, &co2, bb1, &mut next_co1, scratch); 1356 1357 std::mem::swap(&mut r2, &mut next_r2); 1358 std::mem::swap(&mut r1, &mut next_r1); 1359 std::mem::swap(&mut co2, &mut next_co2); 1360 std::mem::swap(&mut co1, &mut next_co1); 1361 1362 if r1.is_negative() { 1363 r1 = r1.negated(); 1364 co1 = co1.negated(); 1365 } 1366 if r2.is_negative() { 1367 r2 = r2.negated(); 1368 co2 = co2.negated(); 1369 } 1370 } 1371 } 1372 1373 if r2.is_negative() { 1374 r2 = r2.negated(); 1375 co2 = co2.negated(); 1376 co1 = co1.negated(); 1377 } 1378 1379 (co2, co1, r2, r1) 1380 } 1381 1382 fn subtract_signed_abs( 1383 positive_limbs: &[u64], 1384 positive_sign: Sign, 1385 negative_limbs: &[u64], 1386 negative_sign: Sign, 1387 ) -> LimbInt { 1388 match cmp_abs_limbs(positive_limbs, negative_limbs) { 1389 Ordering::Greater => { 1390 LimbInt::from_parts(positive_sign, sub_abs_limbs(positive_limbs, negative_limbs)) 1391 } 1392 Ordering::Less => { 1393 LimbInt::from_parts(negative_sign, sub_abs_limbs(negative_limbs, positive_limbs)) 1394 } 1395 Ordering::Equal => LimbInt::zero(), 1396 } 1397 } 1398 1399 fn signed_scalar_sign(value_sign: Sign, scalar: i128) -> Sign { 1400 if scalar == 0 || value_sign == Sign::Zero { 1401 Sign::Zero 1402 } else if scalar < 0 { 1403 match value_sign { 1404 Sign::Negative => Sign::Positive, 1405 Sign::Zero => Sign::Zero, 1406 Sign::Positive => Sign::Negative, 1407 } 1408 } else { 1409 value_sign 1410 } 1411 } 1412 1413 #[inline(always)] 1414 fn combine_signed_abs_limbs_into( 1415 left_sign: Sign, 1416 left: &[u64], 1417 right_sign: Sign, 1418 right: &[u64], 1419 output: &mut LimbVec, 1420 ) -> Sign { 1421 output.clear(); 1422 match (left_sign, right_sign) { 1423 (Sign::Zero, Sign::Zero) => Sign::Zero, 1424 (Sign::Zero, _) => { 1425 output.extend_from_slice(right); 1426 right_sign 1427 } 1428 (_, Sign::Zero) => { 1429 output.extend_from_slice(left); 1430 left_sign 1431 } 1432 (Sign::Positive, Sign::Positive) | (Sign::Negative, Sign::Negative) => { 1433 add_abs_limbs_into(left, right, output); 1434 left_sign 1435 } 1436 (Sign::Positive, Sign::Negative) | (Sign::Negative, Sign::Positive) => { 1437 match cmp_abs_limbs(left, right) { 1438 Ordering::Greater => { 1439 sub_abs_limbs_into(left, right, output); 1440 left_sign 1441 } 1442 Ordering::Less => { 1443 sub_abs_limbs_into(right, left, output); 1444 right_sign 1445 } 1446 Ordering::Equal => Sign::Zero, 1447 } 1448 } 1449 } 1450 } 1451 1452 #[inline(always)] 1453 fn cmp_abs_limbs(left: &[u64], right: &[u64]) -> Ordering { 1454 match left.len().cmp(&right.len()) { 1455 Ordering::Equal => left.iter().rev().cmp(right.iter().rev()), 1456 other => other, 1457 } 1458 } 1459 1460 fn cmp_abs_to_double(left: &[u64], right: &[u64]) -> Ordering { 1461 debug_assert!(!right.is_empty()); 1462 let doubled_len = right.len() + usize::from(right.last().copied().unwrap_or(0) >> 63 != 0); 1463 match left.len().cmp(&doubled_len) { 1464 Ordering::Equal => {} 1465 other => return other, 1466 } 1467 1468 for index in (0..doubled_len).rev() { 1469 let doubled_limb = shifted_left_one_limb_unchecked(right, index); 1470 match left[index].cmp(&doubled_limb) { 1471 Ordering::Equal => {} 1472 other => return other, 1473 } 1474 } 1475 Ordering::Equal 1476 } 1477 1478 fn shifted_left_one_limb_unchecked(limbs: &[u64], index: usize) -> u64 { 1479 debug_assert!(index <= limbs.len()); 1480 let low = if index < limbs.len() { 1481 limbs[index] << 1 1482 } else { 1483 0 1484 }; 1485 let carry = if index == 0 { 1486 0 1487 } else { 1488 limbs[index - 1] >> 63 1489 }; 1490 low | carry 1491 } 1492 1493 fn add_abs_limbs(left: &[u64], right: &[u64]) -> LimbVec { 1494 let mut output = LimbVec::new(); 1495 add_abs_limbs_into(left, right, &mut output); 1496 output 1497 } 1498 1499 #[inline(always)] 1500 #[allow(clippy::uninit_vec)] 1501 fn add_abs_limbs_into(left: &[u64], right: &[u64], output: &mut LimbVec) { 1502 output.clear(); 1503 let len = left.len().max(right.len()); 1504 output.reserve(len + 1); 1505 // SAFETY: u64 has no drop glue, reserve guarantees room for len + 1 limbs, 1506 // and the loops initialize every slot in 0..len before the vector is read. 1507 unsafe { 1508 output.set_len(len); 1509 } 1510 let mut carry = 0_u64; 1511 let shared_len = left.len().min(right.len()); 1512 1513 for index in 0..shared_len { 1514 let (sum, carry_a) = left[index].overflowing_add(right[index]); 1515 let (sum, carry_b) = sum.overflowing_add(carry); 1516 output[index] = sum; 1517 carry = u64::from(carry_a || carry_b); 1518 } 1519 1520 let remaining = if left.len() > right.len() { 1521 &left[shared_len..] 1522 } else { 1523 &right[shared_len..] 1524 }; 1525 for (offset, limb) in remaining.iter().copied().enumerate() { 1526 let (sum, next_carry) = limb.overflowing_add(carry); 1527 output[shared_len + offset] = sum; 1528 carry = u64::from(next_carry); 1529 } 1530 if carry != 0 { 1531 output.push(carry); 1532 } 1533 } 1534 1535 fn add_abs_limbs_assign(left: &mut LimbVec, right: &[u64]) { 1536 let len = left.len().max(right.len()); 1537 let original_left_len = left.len(); 1538 left.resize(len, 0); 1539 let mut carry = 0_u64; 1540 let shared_len = original_left_len.min(right.len()); 1541 1542 for index in 0..shared_len { 1543 let (sum, carry_a) = left[index].overflowing_add(right[index]); 1544 let (sum, carry_b) = sum.overflowing_add(carry); 1545 left[index] = sum; 1546 carry = u64::from(carry_a || carry_b); 1547 } 1548 1549 if right.len() > original_left_len { 1550 for index in shared_len..right.len() { 1551 let (sum, next_carry) = right[index].overflowing_add(carry); 1552 left[index] = sum; 1553 carry = u64::from(next_carry); 1554 } 1555 } else { 1556 for left_limb in left.iter_mut().take(original_left_len).skip(shared_len) { 1557 let (sum, next_carry) = left_limb.overflowing_add(carry); 1558 *left_limb = sum; 1559 carry = u64::from(next_carry); 1560 } 1561 } 1562 1563 if carry != 0 { 1564 left.push(carry); 1565 } 1566 } 1567 1568 fn sub_abs_limbs(left: &[u64], right: &[u64]) -> LimbVec { 1569 let mut output = LimbVec::new(); 1570 sub_abs_limbs_into(left, right, &mut output); 1571 output 1572 } 1573 1574 #[inline(always)] 1575 #[allow(clippy::uninit_vec)] 1576 fn sub_abs_limbs_into(left: &[u64], right: &[u64], output: &mut LimbVec) { 1577 debug_assert!(cmp_abs_limbs(left, right) != Ordering::Less); 1578 output.clear(); 1579 output.reserve(left.len()); 1580 // SAFETY: u64 has no drop glue, reserve guarantees room for left.len() 1581 // limbs, and the loops initialize every slot before trimming reads them. 1582 unsafe { 1583 output.set_len(left.len()); 1584 } 1585 let mut borrow = 0_u64; 1586 let shared_len = right.len(); 1587 1588 for index in 0..shared_len { 1589 let (difference, borrow_a) = left[index].overflowing_sub(right[index]); 1590 let (difference, borrow_b) = difference.overflowing_sub(borrow); 1591 output[index] = difference; 1592 borrow = u64::from(borrow_a || borrow_b); 1593 } 1594 1595 for (offset, left_limb) in left[shared_len..].iter().copied().enumerate() { 1596 let (difference, next_borrow) = left_limb.overflowing_sub(borrow); 1597 output[shared_len + offset] = difference; 1598 borrow = u64::from(next_borrow); 1599 } 1600 1601 debug_assert_eq!(borrow, 0); 1602 trim_leading_zero_limbs(output); 1603 } 1604 1605 fn sub_abs_limbs_assign(left: &mut LimbVec, right: &[u64]) { 1606 debug_assert!(cmp_abs_limbs(left, right) != Ordering::Less); 1607 let mut borrow = 0_u64; 1608 let shared_len = right.len(); 1609 1610 for index in 0..shared_len { 1611 let (difference, borrow_a) = left[index].overflowing_sub(right[index]); 1612 let (difference, borrow_b) = difference.overflowing_sub(borrow); 1613 left[index] = difference; 1614 borrow = u64::from(borrow_a || borrow_b); 1615 } 1616 1617 let mut index = shared_len; 1618 while borrow != 0 && index < left.len() { 1619 let (difference, next_borrow) = left[index].overflowing_sub(borrow); 1620 left[index] = difference; 1621 borrow = u64::from(next_borrow); 1622 index += 1; 1623 } 1624 1625 debug_assert_eq!(borrow, 0); 1626 trim_leading_zero_limbs(left); 1627 } 1628 1629 fn mul_abs_limbs(left: &[u64], right: &[u64]) -> LimbVec { 1630 let mut output = LimbVec::new(); 1631 mul_abs_limbs_into(left, right, &mut output); 1632 output 1633 } 1634 1635 fn mul_abs_limbs_into(left: &[u64], right: &[u64], output: &mut LimbVec) { 1636 output.clear(); 1637 if left.is_empty() || right.is_empty() { 1638 return; 1639 } 1640 1641 let (outer, inner) = if left.len() <= right.len() { 1642 (left, right) 1643 } else { 1644 (right, left) 1645 }; 1646 1647 output.resize(outer.len() + inner.len(), 0); 1648 for (left_index, left_limb) in outer.iter().copied().enumerate() { 1649 let mut carry = 0_u128; 1650 for (right_index, right_limb) in inner.iter().copied().enumerate() { 1651 let output_index = left_index + right_index; 1652 let product = u128::from(left_limb) * u128::from(right_limb) 1653 + u128::from(output[output_index]) 1654 + carry; 1655 output[output_index] = product as u64; 1656 carry = product >> 64; 1657 } 1658 1659 let mut output_index = left_index + inner.len(); 1660 while carry != 0 { 1661 let sum = u128::from(output[output_index]) + carry; 1662 output[output_index] = sum as u64; 1663 carry = sum >> 64; 1664 output_index += 1; 1665 } 1666 } 1667 1668 trim_leading_zero_limbs(output); 1669 } 1670 1671 fn mul_abs_one_limb(value: &[u64], scalar: u64) -> LimbVec { 1672 let mut output = LimbVec::new(); 1673 mul_abs_one_limb_into(value, scalar, &mut output); 1674 output 1675 } 1676 1677 #[inline(always)] 1678 #[allow(clippy::uninit_vec)] 1679 fn mul_abs_one_limb_into(value: &[u64], scalar: u64, output: &mut LimbVec) { 1680 output.clear(); 1681 if value.is_empty() || scalar == 0 { 1682 return; 1683 } 1684 if scalar == 1 { 1685 output.extend_from_slice(value); 1686 return; 1687 } 1688 1689 output.reserve(value.len() + 1); 1690 let mut carry = 0_u128; 1691 // SAFETY: u64 has no drop glue, reserve guarantees room for value.len() + 1 1692 // limbs, and the loop initializes every slot before the vector is read. 1693 unsafe { 1694 output.set_len(value.len() + 1); 1695 } 1696 for (index, limb) in value.iter().copied().enumerate() { 1697 let product = u128::from(limb) * u128::from(scalar) + carry; 1698 output[index] = product as u64; 1699 carry = product >> 64; 1700 } 1701 if carry != 0 { 1702 output[value.len()] = carry as u64; 1703 } else { 1704 output.truncate(value.len()); 1705 } 1706 } 1707 1708 fn square_abs_limbs(value: &[u64]) -> LimbVec { 1709 mul_abs_limbs(value, value) 1710 } 1711 1712 fn mul_abs_small(value: &[u64], scalar: u128) -> LimbVec { 1713 let mut output = LimbVec::new(); 1714 mul_abs_small_into(value, scalar, &mut output); 1715 output 1716 } 1717 1718 #[inline(always)] 1719 fn mul_abs_small_into(value: &[u64], scalar: u128, output: &mut LimbVec) { 1720 output.clear(); 1721 if value.is_empty() || scalar == 0 { 1722 return; 1723 } 1724 if scalar == 1 { 1725 output.extend_from_slice(value); 1726 return; 1727 } 1728 let scalar_low = scalar as u64; 1729 let scalar_high = (scalar >> 64) as u64; 1730 if scalar_high == 0 { 1731 mul_abs_one_limb_into(value, scalar_low, output); 1732 return; 1733 } 1734 1735 output.resize(value.len() + 2, 0); 1736 1737 if scalar_low != 0 { 1738 for (index, limb) in value.iter().copied().enumerate() { 1739 add_u128_at(output, index, u128::from(limb) * u128::from(scalar_low)); 1740 } 1741 } 1742 if scalar_high != 0 { 1743 for (index, limb) in value.iter().copied().enumerate() { 1744 add_u128_at( 1745 output, 1746 index + 1, 1747 u128::from(limb) * u128::from(scalar_high), 1748 ); 1749 } 1750 } 1751 1752 trim_leading_zero_limbs(output); 1753 } 1754 1755 fn add_u128_at(output: &mut LimbVec, index: usize, value: u128) { 1756 let low = value as u64; 1757 let high = (value >> 64) as u64; 1758 ensure_len(output, index + 2); 1759 1760 let (sum_low, carry_low) = output[index].overflowing_add(low); 1761 output[index] = sum_low; 1762 1763 let (sum_high, carry_high_a) = output[index + 1].overflowing_add(high); 1764 let (sum_high, carry_high_b) = sum_high.overflowing_add(u64::from(carry_low)); 1765 output[index + 1] = sum_high; 1766 1767 let mut carry = u64::from(carry_high_a) + u64::from(carry_high_b); 1768 let mut carry_index = index + 2; 1769 while carry != 0 { 1770 ensure_len(output, carry_index + 1); 1771 let (sum, overflowed) = output[carry_index].overflowing_add(carry); 1772 output[carry_index] = sum; 1773 carry = u64::from(overflowed); 1774 carry_index += 1; 1775 } 1776 } 1777 1778 fn ensure_len(output: &mut LimbVec, len: usize) { 1779 if output.len() < len { 1780 output.resize(len, 0); 1781 } 1782 } 1783 1784 fn div_rem_abs_limbs(numerator: &[u64], denominator: &[u64]) -> (LimbVec, LimbVec) { 1785 let mut scratch = DivisionScratch::default(); 1786 div_rem_abs_limbs_with_scratch(numerator, denominator, &mut scratch) 1787 } 1788 1789 fn div_rem_abs_limbs_with_scratch( 1790 numerator: &[u64], 1791 denominator: &[u64], 1792 scratch: &mut DivisionScratch, 1793 ) -> (LimbVec, LimbVec) { 1794 div_rem_abs_limbs_into_scratch(numerator, denominator, scratch, true, true); 1795 (scratch.quotient.clone(), scratch.remainder.clone()) 1796 } 1797 1798 fn div_abs_limbs_with_scratch( 1799 numerator: &[u64], 1800 denominator: &[u64], 1801 scratch: &mut DivisionScratch, 1802 ) -> LimbVec { 1803 div_rem_abs_limbs_into_scratch(numerator, denominator, scratch, true, false); 1804 scratch.quotient.clone() 1805 } 1806 1807 fn rem_abs_limbs_with_scratch( 1808 numerator: &[u64], 1809 denominator: &[u64], 1810 scratch: &mut DivisionScratch, 1811 ) -> LimbVec { 1812 div_rem_abs_limbs_into_scratch(numerator, denominator, scratch, false, true); 1813 scratch.remainder.clone() 1814 } 1815 1816 #[inline(always)] 1817 fn div_rem_abs_limbs_into_scratch( 1818 numerator: &[u64], 1819 denominator: &[u64], 1820 scratch: &mut DivisionScratch, 1821 keep_quotient: bool, 1822 keep_remainder: bool, 1823 ) { 1824 debug_assert!(!denominator.is_empty()); 1825 scratch.quotient.clear(); 1826 scratch.remainder.clear(); 1827 1828 if numerator.is_empty() { 1829 return; 1830 } 1831 if cmp_abs_limbs(numerator, denominator) == Ordering::Less { 1832 if keep_remainder { 1833 scratch.remainder.extend_from_slice(numerator); 1834 } 1835 return; 1836 } 1837 if denominator == [1] { 1838 if keep_quotient { 1839 scratch.quotient.extend_from_slice(numerator); 1840 } 1841 return; 1842 } 1843 if denominator.len() == 1 { 1844 div_rem_abs_one_limb_into( 1845 numerator, 1846 denominator[0], 1847 &mut scratch.quotient, 1848 &mut scratch.remainder, 1849 keep_quotient, 1850 keep_remainder, 1851 ); 1852 return; 1853 } 1854 1855 let shift = denominator.last().copied().unwrap_or(0).leading_zeros() as usize; 1856 shl_abs_limbs_into(numerator, shift, &mut scratch.normalized_numerator); 1857 shl_abs_limbs_into(denominator, shift, &mut scratch.normalized_denominator); 1858 scratch.normalized_numerator.push(0); 1859 1860 let denominator_len = scratch.normalized_denominator.len(); 1861 let quotient_len = scratch.normalized_numerator.len() - denominator_len; 1862 if keep_quotient { 1863 scratch.quotient.resize(quotient_len, 0); 1864 } 1865 1866 for quotient_index in (0..quotient_len).rev() { 1867 let (mut qhat, mut rhat, mut rhat_overflowed) = estimate_quotient( 1868 scratch.normalized_numerator[quotient_index + denominator_len], 1869 scratch.normalized_numerator[quotient_index + denominator_len - 1], 1870 scratch.normalized_denominator[denominator_len - 1], 1871 ); 1872 1873 if denominator_len > 1 { 1874 let next_denominator_limb = scratch.normalized_denominator[denominator_len - 2]; 1875 let next_numerator_limb = 1876 scratch.normalized_numerator[quotient_index + denominator_len - 2]; 1877 while !rhat_overflowed 1878 && quotient_too_large(qhat, rhat, next_denominator_limb, next_numerator_limb) 1879 { 1880 qhat -= 1; 1881 let (next_rhat, overflowed) = 1882 rhat.overflowing_add(scratch.normalized_denominator[denominator_len - 1]); 1883 rhat = next_rhat; 1884 rhat_overflowed = overflowed; 1885 } 1886 } 1887 1888 if sub_mul_at( 1889 &mut scratch.normalized_numerator, 1890 &scratch.normalized_denominator, 1891 qhat, 1892 quotient_index, 1893 ) { 1894 qhat -= 1; 1895 add_at( 1896 &mut scratch.normalized_numerator, 1897 &scratch.normalized_denominator, 1898 quotient_index, 1899 ); 1900 } 1901 if keep_quotient { 1902 scratch.quotient[quotient_index] = qhat; 1903 } 1904 } 1905 1906 if keep_remainder { 1907 shr_abs_limbs_into( 1908 &scratch.normalized_numerator[..denominator_len], 1909 shift, 1910 &mut scratch.remainder, 1911 ); 1912 trim_leading_zero_limbs(&mut scratch.remainder); 1913 } 1914 if keep_quotient { 1915 trim_leading_zero_limbs(&mut scratch.quotient); 1916 } 1917 } 1918 1919 fn div_rem_abs_one_limb(numerator: &[u64], denominator: u64) -> (LimbVec, LimbVec) { 1920 let mut quotient = LimbVec::new(); 1921 let mut remainder = LimbVec::new(); 1922 div_rem_abs_one_limb_into( 1923 numerator, 1924 denominator, 1925 &mut quotient, 1926 &mut remainder, 1927 true, 1928 true, 1929 ); 1930 (quotient, remainder) 1931 } 1932 1933 fn div_rem_abs_one_limb_into( 1934 numerator: &[u64], 1935 denominator: u64, 1936 quotient: &mut LimbVec, 1937 remainder_output: &mut LimbVec, 1938 keep_quotient: bool, 1939 keep_remainder: bool, 1940 ) { 1941 debug_assert_ne!(denominator, 0); 1942 quotient.clear(); 1943 remainder_output.clear(); 1944 if keep_quotient { 1945 quotient.resize(numerator.len(), 0); 1946 } 1947 let mut remainder = 0_u128; 1948 1949 for (index, limb) in numerator.iter().copied().enumerate().rev() { 1950 let value = (remainder << 64) | u128::from(limb); 1951 if keep_quotient { 1952 quotient[index] = (value / u128::from(denominator)) as u64; 1953 } 1954 remainder = value % u128::from(denominator); 1955 } 1956 1957 if keep_quotient { 1958 trim_leading_zero_limbs(quotient); 1959 } 1960 if keep_remainder && remainder != 0 { 1961 remainder_output.push(remainder as u64); 1962 } 1963 } 1964 1965 #[inline(always)] 1966 fn estimate_quotient(high: u64, low: u64, denominator_high: u64) -> (u64, u64, bool) { 1967 if high == denominator_high { 1968 let (remainder, overflowed) = low.overflowing_add(denominator_high); 1969 (u64::MAX, remainder, overflowed) 1970 } else { 1971 let numerator = (u128::from(high) << 64) | u128::from(low); 1972 ( 1973 (numerator / u128::from(denominator_high)) as u64, 1974 (numerator % u128::from(denominator_high)) as u64, 1975 false, 1976 ) 1977 } 1978 } 1979 1980 #[inline(always)] 1981 fn quotient_too_large(qhat: u64, rhat: u64, denominator_next: u64, numerator_next: u64) -> bool { 1982 let left = u128::from(qhat) * u128::from(denominator_next); 1983 let right = (u128::from(rhat) << 64) | u128::from(numerator_next); 1984 left > right 1985 } 1986 1987 #[inline(always)] 1988 fn sub_mul_at(target: &mut [u64], value: &[u64], multiplier: u64, offset: usize) -> bool { 1989 if multiplier == 0 { 1990 return false; 1991 } 1992 1993 let mut carry = 0_u128; 1994 for (index, value_limb) in value.iter().copied().enumerate() { 1995 let product = u128::from(multiplier) * u128::from(value_limb) + carry; 1996 let product_low = product as u64; 1997 carry = product >> 64; 1998 1999 let target_index = offset + index; 2000 let (difference, borrowed) = target[target_index].overflowing_sub(product_low); 2001 target[target_index] = difference; 2002 carry += u128::from(borrowed); 2003 } 2004 2005 let target_index = offset + value.len(); 2006 let carry_low = carry as u64; 2007 let carry_high = carry >> 64; 2008 let (difference, borrowed) = target[target_index].overflowing_sub(carry_low); 2009 target[target_index] = difference; 2010 borrowed || carry_high != 0 2011 } 2012 2013 #[inline(always)] 2014 fn add_at(target: &mut [u64], value: &[u64], offset: usize) { 2015 let mut carry = 0_u128; 2016 for (index, value_limb) in value.iter().copied().enumerate() { 2017 let target_index = offset + index; 2018 let sum = u128::from(target[target_index]) + u128::from(value_limb) + carry; 2019 target[target_index] = sum as u64; 2020 carry = sum >> 64; 2021 } 2022 2023 let mut target_index = offset + value.len(); 2024 while carry != 0 { 2025 let sum = u128::from(target[target_index]) + carry; 2026 target[target_index] = sum as u64; 2027 carry = sum >> 64; 2028 target_index += 1; 2029 } 2030 } 2031 2032 fn shl_abs_limbs(value: &[u64], bits: usize) -> LimbVec { 2033 let mut output = LimbVec::new(); 2034 shl_abs_limbs_into(value, bits, &mut output); 2035 output 2036 } 2037 2038 fn shl_abs_limbs_into(value: &[u64], bits: usize, output: &mut LimbVec) { 2039 output.clear(); 2040 if value.is_empty() { 2041 return; 2042 } 2043 if bits == 0 { 2044 output.extend_from_slice(value); 2045 return; 2046 } 2047 2048 let limb_shift = bits / 64; 2049 let bit_shift = (bits % 64) as u32; 2050 if limb_shift == 0 && bit_shift != 0 { 2051 shl_abs_limbs_small_into(value, bit_shift, output); 2052 return; 2053 } 2054 output.resize(limb_shift + value.len() + usize::from(bit_shift != 0), 0); 2055 let mut carry = 0_u64; 2056 2057 for (index, limb) in value.iter().copied().enumerate() { 2058 output[index + limb_shift] = if bit_shift == 0 { 2059 limb 2060 } else { 2061 (limb << bit_shift) | carry 2062 }; 2063 carry = if bit_shift == 0 { 2064 0 2065 } else { 2066 limb >> (64 - bit_shift) 2067 }; 2068 } 2069 if bit_shift != 0 { 2070 output[limb_shift + value.len()] = carry; 2071 } 2072 2073 trim_leading_zero_limbs(output); 2074 } 2075 2076 #[allow(clippy::uninit_vec)] 2077 fn shl_abs_limbs_small_into(value: &[u64], bit_shift: u32, output: &mut LimbVec) { 2078 debug_assert!((1..64).contains(&bit_shift)); 2079 output.reserve(value.len() + 1); 2080 // SAFETY: u64 has no drop glue, reserve guarantees room for value.len() + 1 2081 // limbs, and the loop initializes every slot before the vector is read. 2082 unsafe { 2083 output.set_len(value.len() + 1); 2084 } 2085 let mut carry = 0_u64; 2086 for (index, limb) in value.iter().copied().enumerate() { 2087 output[index] = (limb << bit_shift) | carry; 2088 carry = limb >> (64 - bit_shift); 2089 } 2090 if carry != 0 { 2091 output[value.len()] = carry; 2092 } else { 2093 output.truncate(value.len()); 2094 } 2095 } 2096 2097 fn shl_abs_limbs_assign(value: &mut LimbVec, bits: usize) { 2098 if value.is_empty() || bits == 0 { 2099 return; 2100 } 2101 2102 let limb_shift = bits / 64; 2103 let bit_shift = (bits % 64) as u32; 2104 if limb_shift == 0 && bit_shift != 0 { 2105 shl_abs_limbs_small_assign(value, bit_shift); 2106 return; 2107 } 2108 let original_len = value.len(); 2109 value.resize(original_len + limb_shift + usize::from(bit_shift != 0), 0); 2110 2111 for index in (0..original_len).rev() { 2112 let limb = value[index]; 2113 let output_index = index + limb_shift; 2114 value[output_index] = limb << bit_shift; 2115 if bit_shift != 0 { 2116 value[output_index + 1] |= limb >> (64 - bit_shift); 2117 } 2118 } 2119 2120 value[..limb_shift].fill(0); 2121 trim_leading_zero_limbs(value); 2122 } 2123 2124 fn shl_abs_limbs_small_assign(value: &mut LimbVec, bit_shift: u32) { 2125 debug_assert!((1..64).contains(&bit_shift)); 2126 let mut carry = 0_u64; 2127 for limb in value.iter_mut() { 2128 let current = *limb; 2129 *limb = (current << bit_shift) | carry; 2130 carry = current >> (64 - bit_shift); 2131 } 2132 if carry != 0 { 2133 value.push(carry); 2134 } 2135 } 2136 2137 fn shr_abs_limbs(value: &[u64], bits: usize) -> LimbVec { 2138 let mut output = LimbVec::new(); 2139 shr_abs_limbs_into(value, bits, &mut output); 2140 output 2141 } 2142 2143 fn shr_abs_limbs_into(value: &[u64], bits: usize, output: &mut LimbVec) { 2144 output.clear(); 2145 let limb_shift = bits / 64; 2146 if limb_shift >= value.len() { 2147 return; 2148 } 2149 2150 let bit_shift = (bits % 64) as u32; 2151 if limb_shift == 0 && bit_shift != 0 { 2152 shr_abs_limbs_small_into(value, bit_shift, output); 2153 return; 2154 } 2155 output.reserve(value.len() - limb_shift); 2156 for index in limb_shift..value.len() { 2157 let mut limb = value[index] >> bit_shift; 2158 if bit_shift != 0 { 2159 limb |= value.get(index + 1).copied().unwrap_or(0) << (64 - bit_shift); 2160 } 2161 output.push(limb); 2162 } 2163 2164 trim_leading_zero_limbs(output); 2165 } 2166 2167 fn shr_abs_limbs_small_into(value: &[u64], bit_shift: u32, output: &mut LimbVec) { 2168 debug_assert!((1..64).contains(&bit_shift)); 2169 output.reserve(value.len()); 2170 for index in 0..value.len() { 2171 let limb = (value[index] >> bit_shift) 2172 | (value.get(index + 1).copied().unwrap_or(0) << (64 - bit_shift)); 2173 output.push(limb); 2174 } 2175 trim_leading_zero_limbs(output); 2176 } 2177 2178 fn trim_leading_zero_limbs(limbs: &mut LimbVec) { 2179 while limbs.last() == Some(&0) { 2180 limbs.pop(); 2181 } 2182 } 2183 2184 fn abs_bit_len(limbs: &[u64]) -> usize { 2185 let Some(last) = limbs.last() else { 2186 return 0; 2187 }; 2188 ((limbs.len() - 1) * 64) + (u64::BITS - last.leading_zeros()) as usize 2189 } 2190 2191 fn get_abs_bit(limbs: &[u64], bit: usize) -> bool { 2192 let limb = bit / 64; 2193 let offset = bit % 64; 2194 limbs 2195 .get(limb) 2196 .map(|value| (value & (1_u64 << offset)) != 0) 2197 .unwrap_or(false) 2198 } 2199 2200 fn set_abs_bit(limbs: &mut LimbVec, bit: usize) { 2201 let limb = bit / 64; 2202 let offset = bit % 64; 2203 if limbs.len() <= limb { 2204 limbs.resize(limb + 1, 0); 2205 } 2206 limbs[limb] |= 1_u64 << offset; 2207 } 2208 2209 #[cfg(test)] 2210 mod tests { 2211 use num_bigint::BigInt; 2212 use num_integer::Integer; 2213 use num_traits::Signed; 2214 use proptest::prelude::*; 2215 2216 use super::{LimbInt, LimbScratch}; 2217 2218 fn arb_bigint() -> impl Strategy<Value = BigInt> { 2219 proptest::collection::vec(any::<u64>(), 0..=24).prop_flat_map(|limbs| { 2220 any::<bool>().prop_map(move |negative| { 2221 let mut bytes = Vec::with_capacity(limbs.len() * 8); 2222 for limb in &limbs { 2223 bytes.extend_from_slice(&limb.to_le_bytes()); 2224 } 2225 let magnitude = num_bigint::BigUint::from_bytes_le(&bytes); 2226 let sign = if magnitude == num_bigint::BigUint::from(0_u8) { 2227 num_bigint::Sign::NoSign 2228 } else if negative { 2229 num_bigint::Sign::Minus 2230 } else { 2231 num_bigint::Sign::Plus 2232 }; 2233 BigInt::from_biguint(sign, magnitude) 2234 }) 2235 }) 2236 } 2237 2238 proptest! { 2239 #![proptest_config(ProptestConfig { cases: 512, .. ProptestConfig::default() })] 2240 2241 #[test] 2242 fn roundtrips_bigint(value in arb_bigint()) { 2243 prop_assert_eq!(LimbInt::from_bigint(&value).to_bigint(), value); 2244 } 2245 2246 #[test] 2247 fn add_matches_bigint(left in arb_bigint(), right in arb_bigint()) { 2248 let actual = LimbInt::from_bigint(&left).add(&LimbInt::from_bigint(&right)).to_bigint(); 2249 prop_assert_eq!(actual, left + right); 2250 } 2251 2252 #[test] 2253 fn sub_matches_bigint(left in arb_bigint(), right in arb_bigint()) { 2254 let actual = LimbInt::from_bigint(&left).sub(&LimbInt::from_bigint(&right)).to_bigint(); 2255 prop_assert_eq!(actual, left - right); 2256 } 2257 2258 #[test] 2259 fn mul_matches_bigint(left in arb_bigint(), right in arb_bigint()) { 2260 let actual = LimbInt::from_bigint(&left).mul(&LimbInt::from_bigint(&right)).to_bigint(); 2261 prop_assert_eq!(actual, left * right); 2262 } 2263 2264 #[test] 2265 fn square_matches_bigint(value in arb_bigint()) { 2266 let actual = LimbInt::from_bigint(&value).square().to_bigint(); 2267 prop_assert_eq!(actual, &value * &value); 2268 } 2269 2270 #[test] 2271 fn shifts_match_positive_bigint(value in arb_bigint(), bits in 0_usize..512) { 2272 let value = value.abs(); 2273 let limbs = LimbInt::from_bigint(&value); 2274 prop_assert_eq!(limbs.shl_bits(bits).to_bigint(), &value << bits); 2275 prop_assert_eq!(limbs.shr_abs_bits(bits).to_bigint(), &value >> bits); 2276 2277 let mut owned_shift = limbs.clone(); 2278 owned_shift = owned_shift.shl_bits_owned(bits); 2279 prop_assert_eq!(owned_shift.to_bigint(), &value << bits); 2280 } 2281 2282 #[test] 2283 fn shifted_low_word_matches_num_bigint(value in arb_bigint(), shift in 0_u64..1536) { 2284 let value = value.abs(); 2285 let limbs = LimbInt::from_bigint(&value); 2286 let mut digits = value.iter_u64_digits(); 2287 let digit_index = usize::try_from(shift / 64).unwrap_or(usize::MAX); 2288 let offset = (shift % 64) as u32; 2289 let low = digits.nth(digit_index).unwrap_or(0); 2290 let expected = if offset == 0 { 2291 low 2292 } else { 2293 let high = digits.next().unwrap_or(0); 2294 (low >> offset) | (high << (64 - offset)) 2295 }; 2296 2297 prop_assert_eq!(limbs.shifted_low_word(shift), expected); 2298 } 2299 2300 #[test] 2301 fn div_rem_matches_bigint(left in arb_bigint(), right in arb_bigint().prop_filter("non-zero divisor", |value| value != &BigInt::from(0))) { 2302 let left_limbs = LimbInt::from_bigint(&left); 2303 let right_limbs = LimbInt::from_bigint(&right); 2304 let (actual_q, actual_r) = left_limbs.div_rem(&right_limbs); 2305 let mut scratch = LimbScratch::default(); 2306 let actual_div = left_limbs.div_with_scratch(&right_limbs, &mut scratch); 2307 let actual_floor = left_limbs.div_floor_with_scratch(&right_limbs, &mut scratch); 2308 let (expected_q, expected_r) = left.div_rem(&right); 2309 2310 prop_assert_eq!(actual_q.to_bigint(), expected_q); 2311 prop_assert_eq!(actual_r.to_bigint(), expected_r); 2312 prop_assert_eq!(actual_div.to_bigint(), actual_q.to_bigint()); 2313 prop_assert_eq!(actual_floor.to_bigint(), left.div_floor(&right)); 2314 } 2315 2316 #[test] 2317 fn extended_gcd_matches_bezout( 2318 left in arb_bigint(), 2319 right in arb_bigint().prop_filter("not both zero", |right| right != &BigInt::from(0)), 2320 ) { 2321 let left_limbs = LimbInt::from_bigint(&left); 2322 let right_limbs = LimbInt::from_bigint(&right); 2323 let actual = left_limbs.extended_gcd(&right_limbs); 2324 let actual_x = actual.x.to_bigint(); 2325 let actual_y = actual.y.to_bigint(); 2326 let actual_gcd = actual.gcd.to_bigint(); 2327 2328 prop_assert_eq!(&left * actual_x + &right * actual_y, actual_gcd.clone()); 2329 prop_assert_eq!(actual_gcd.clone(), left.gcd(&right)); 2330 } 2331 2332 #[test] 2333 fn left_extended_gcd_matches_modular_bezout( 2334 left in arb_bigint(), 2335 right in arb_bigint().prop_filter("non-zero modulus", |value| value != &BigInt::from(0)), 2336 ) { 2337 let left_limbs = LimbInt::from_bigint(&left); 2338 let right_limbs = LimbInt::from_bigint(&right); 2339 let mut scratch = LimbScratch::default(); 2340 let actual = left_limbs.left_extended_gcd_with_scratch(&right_limbs, &mut scratch); 2341 let actual_x = actual.x.to_bigint(); 2342 let actual_gcd = actual.gcd.to_bigint(); 2343 2344 prop_assert_eq!(actual_gcd.clone(), left.gcd(&right)); 2345 prop_assert_eq!((&left * actual_x - actual_gcd) % &right, BigInt::from(0)); 2346 } 2347 2348 #[test] 2349 fn mod_positive_matches_corrected_bigint(value in arb_bigint(), modulus in arb_bigint().prop_filter("positive modulus", |value| value > &BigInt::from(0))) { 2350 let value_limbs = LimbInt::from_bigint(&value); 2351 let modulus_limbs = LimbInt::from_bigint(&modulus); 2352 let mut expected = value % &modulus; 2353 if expected.is_negative() { 2354 expected += &modulus; 2355 } 2356 2357 prop_assert_eq!(value_limbs.mod_positive(&modulus_limbs).to_bigint(), expected); 2358 } 2359 } 2360 2361 #[test] 2362 fn from_i128_handles_i128_min() { 2363 assert_eq!( 2364 LimbInt::from_i128(i128::MIN).to_bigint(), 2365 BigInt::from(i128::MIN) 2366 ); 2367 } 2368 2369 #[test] 2370 fn xgcd_partial_matches_bigint_reference() { 2371 let mut state = 0x811c_9dc5_0123_4567_u64; 2372 2373 for case in 0..128 { 2374 let mut r2 = next_positive_bigint(&mut state, 6); 2375 let mut r1 = next_positive_bigint(&mut state, 6); 2376 if r2 < r1 { 2377 std::mem::swap(&mut r2, &mut r1); 2378 } 2379 2380 let threshold_seed = next_positive_bigint(&mut state, 3); 2381 let threshold = (threshold_seed % &r1).max(BigInt::from(1)); 2382 2383 let expected = crate::domain::vdf::arithmetic::xgcd_partial(&r2, &r1, &threshold); 2384 let actual = super::xgcd_partial( 2385 &LimbInt::from_bigint(&r2), 2386 &LimbInt::from_bigint(&r1), 2387 &LimbInt::from_bigint(&threshold), 2388 ); 2389 2390 assert_eq!(actual.0.to_bigint(), expected.0, "case={case} co2"); 2391 assert_eq!(actual.1.to_bigint(), expected.1, "case={case} co1"); 2392 assert_eq!(actual.2.to_bigint(), expected.2, "case={case} r2"); 2393 assert_eq!(actual.3.to_bigint(), expected.3, "case={case} r1"); 2394 } 2395 } 2396 2397 fn next_positive_bigint(state: &mut u64, max_limbs: usize) -> BigInt { 2398 let limb_count = (next_u64(state) as usize % max_limbs) + 1; 2399 let mut bytes = Vec::with_capacity(limb_count * 8); 2400 for _ in 0..limb_count { 2401 bytes.extend_from_slice(&next_u64(state).to_le_bytes()); 2402 } 2403 let value = num_bigint::BigUint::from_bytes_le(&bytes); 2404 BigInt::from(value.max(num_bigint::BigUint::from(1_u8))) 2405 } 2406 2407 fn next_u64(state: &mut u64) -> u64 { 2408 *state = state 2409 .wrapping_mul(6_364_136_223_846_793_005) 2410 .wrapping_add(1); 2411 *state 2412 } 2413 }