1#![doc(html_root_url = "https://docs.rs/num-integer/0.1")]
18#![no_std]
19
20use core::mem;
21use core::ops::Add;
22
23use num_traits::{Num, Signed, Zero};
24
25mod roots;
26pub use crate::roots::Roots;
27pub use crate::roots::{cbrt, nth_root, sqrt};
28
29mod average;
30pub use crate::average::Average;
31pub use crate::average::{average_ceil, average_floor};
32
33pub trait Integer: Sized + Num + PartialOrd + Ord + Eq {
34 fn div_floor(&self, other: &Self) -> Self;
51
52 fn mod_floor(&self, other: &Self) -> Self;
75
76 fn div_ceil(&self, other: &Self) -> Self {
93 let (q, r) = self.div_mod_floor(other);
94 if r.is_zero() {
95 q
96 } else {
97 q + Self::one()
98 }
99 }
100
101 fn gcd(&self, other: &Self) -> Self;
111
112 fn lcm(&self, other: &Self) -> Self;
123
124 #[inline]
138 fn gcd_lcm(&self, other: &Self) -> (Self, Self) {
139 (self.gcd(other), self.lcm(other))
140 }
141
142 #[inline]
159 fn extended_gcd(&self, other: &Self) -> ExtendedGcd<Self>
160 where
161 Self: Clone,
162 {
163 let mut s = (Self::zero(), Self::one());
164 let mut t = (Self::one(), Self::zero());
165 let mut r = (other.clone(), self.clone());
166
167 while !r.0.is_zero() {
168 let q = r.1.clone() / r.0.clone();
169 let f = |mut r: (Self, Self)| {
170 mem::swap(&mut r.0, &mut r.1);
171 r.0 = r.0 - q.clone() * r.1.clone();
172 r
173 };
174 r = f(r);
175 s = f(s);
176 t = f(t);
177 }
178
179 if r.1 >= Self::zero() {
180 ExtendedGcd {
181 gcd: r.1,
182 x: s.1,
183 y: t.1,
184 }
185 } else {
186 ExtendedGcd {
187 gcd: Self::zero() - r.1,
188 x: Self::zero() - s.1,
189 y: Self::zero() - t.1,
190 }
191 }
192 }
193
194 #[inline]
196 fn extended_gcd_lcm(&self, other: &Self) -> (ExtendedGcd<Self>, Self)
197 where
198 Self: Clone + Signed,
199 {
200 (self.extended_gcd(other), self.lcm(other))
201 }
202
203 #[deprecated(note = "Please use is_multiple_of instead")]
205 #[inline]
206 fn divides(&self, other: &Self) -> bool {
207 self.is_multiple_of(other)
208 }
209
210 fn is_multiple_of(&self, other: &Self) -> bool;
220
221 fn is_even(&self) -> bool;
231
232 fn is_odd(&self) -> bool;
242
243 fn div_rem(&self, other: &Self) -> (Self, Self);
261
262 fn div_mod_floor(&self, other: &Self) -> (Self, Self) {
280 (self.div_floor(other), self.mod_floor(other))
281 }
282
283 #[inline]
303 fn next_multiple_of(&self, other: &Self) -> Self
304 where
305 Self: Clone,
306 {
307 let m = self.mod_floor(other);
308 self.clone()
309 + if m.is_zero() {
310 Self::zero()
311 } else {
312 other.clone() - m
313 }
314 }
315
316 #[inline]
336 fn prev_multiple_of(&self, other: &Self) -> Self
337 where
338 Self: Clone,
339 {
340 self.clone() - self.mod_floor(other)
341 }
342
343 fn dec(&mut self)
354 where
355 Self: Clone,
356 {
357 *self = self.clone() - Self::one()
358 }
359
360 fn inc(&mut self)
371 where
372 Self: Clone,
373 {
374 *self = self.clone() + Self::one()
375 }
376}
377
378#[derive(Debug, Clone, Copy, PartialEq, Eq)]
385pub struct ExtendedGcd<A> {
386 pub gcd: A,
387 pub x: A,
388 pub y: A,
389}
390
391#[inline]
393pub fn div_rem<T: Integer>(x: T, y: T) -> (T, T) {
394 x.div_rem(&y)
395}
396#[inline]
398pub fn div_floor<T: Integer>(x: T, y: T) -> T {
399 x.div_floor(&y)
400}
401#[inline]
403pub fn mod_floor<T: Integer>(x: T, y: T) -> T {
404 x.mod_floor(&y)
405}
406#[inline]
408pub fn div_mod_floor<T: Integer>(x: T, y: T) -> (T, T) {
409 x.div_mod_floor(&y)
410}
411#[inline]
413pub fn div_ceil<T: Integer>(x: T, y: T) -> T {
414 x.div_ceil(&y)
415}
416
417#[inline(always)]
420pub fn gcd<T: Integer>(x: T, y: T) -> T {
421 x.gcd(&y)
422}
423#[inline(always)]
425pub fn lcm<T: Integer>(x: T, y: T) -> T {
426 x.lcm(&y)
427}
428
429#[inline(always)]
432pub fn gcd_lcm<T: Integer>(x: T, y: T) -> (T, T) {
433 x.gcd_lcm(&y)
434}
435
436macro_rules! impl_integer_for_isize {
437 ($T:ty, $test_mod:ident) => {
438 impl Integer for $T {
439 #[inline]
441 fn div_floor(&self, other: &Self) -> Self {
442 let (d, r) = self.div_rem(other);
445 if (r > 0 && *other < 0) || (r < 0 && *other > 0) {
446 d - 1
447 } else {
448 d
449 }
450 }
451
452 #[inline]
454 fn mod_floor(&self, other: &Self) -> Self {
455 let r = *self % *other;
458 if (r > 0 && *other < 0) || (r < 0 && *other > 0) {
459 r + *other
460 } else {
461 r
462 }
463 }
464
465 #[inline]
467 fn div_mod_floor(&self, other: &Self) -> (Self, Self) {
468 let (d, r) = self.div_rem(other);
471 if (r > 0 && *other < 0) || (r < 0 && *other > 0) {
472 (d - 1, r + *other)
473 } else {
474 (d, r)
475 }
476 }
477
478 #[inline]
479 fn div_ceil(&self, other: &Self) -> Self {
480 let (d, r) = self.div_rem(other);
481 if (r > 0 && *other > 0) || (r < 0 && *other < 0) {
482 d + 1
483 } else {
484 d
485 }
486 }
487
488 #[inline]
491 fn gcd(&self, other: &Self) -> Self {
492 let mut m = *self;
494 let mut n = *other;
495 if m == 0 || n == 0 {
496 return (m | n).abs();
497 }
498
499 let shift = (m | n).trailing_zeros();
501
502 if m == Self::min_value() || n == Self::min_value() {
511 return (1 << shift).abs();
512 }
513
514 m = m.abs();
516 n = n.abs();
517
518 m >>= m.trailing_zeros();
520 n >>= n.trailing_zeros();
521
522 while m != n {
523 if m > n {
524 m -= n;
525 m >>= m.trailing_zeros();
526 } else {
527 n -= m;
528 n >>= n.trailing_zeros();
529 }
530 }
531 m << shift
532 }
533
534 #[inline]
535 fn extended_gcd_lcm(&self, other: &Self) -> (ExtendedGcd<Self>, Self) {
536 let egcd = self.extended_gcd(other);
537 let lcm = if egcd.gcd.is_zero() {
539 Self::zero()
540 } else {
541 (*self * (*other / egcd.gcd)).abs()
542 };
543 (egcd, lcm)
544 }
545
546 #[inline]
549 fn lcm(&self, other: &Self) -> Self {
550 self.gcd_lcm(other).1
551 }
552
553 #[inline]
556 fn gcd_lcm(&self, other: &Self) -> (Self, Self) {
557 if self.is_zero() && other.is_zero() {
558 return (Self::zero(), Self::zero());
559 }
560 let gcd = self.gcd(other);
561 let lcm = (*self * (*other / gcd)).abs();
563 (gcd, lcm)
564 }
565
566 #[inline]
568 fn is_multiple_of(&self, other: &Self) -> bool {
569 if other.is_zero() {
570 return self.is_zero();
571 }
572 *self % *other == 0
573 }
574
575 #[inline]
577 fn is_even(&self) -> bool {
578 (*self) & 1 == 0
579 }
580
581 #[inline]
583 fn is_odd(&self) -> bool {
584 !self.is_even()
585 }
586
587 #[inline]
589 fn div_rem(&self, other: &Self) -> (Self, Self) {
590 (*self / *other, *self % *other)
591 }
592
593 #[inline]
595 fn next_multiple_of(&self, other: &Self) -> Self {
596 if *other == -1 {
598 return *self;
599 }
600
601 let m = Integer::mod_floor(self, other);
602 *self + if m == 0 { 0 } else { other - m }
603 }
604
605 #[inline]
607 fn prev_multiple_of(&self, other: &Self) -> Self {
608 if *other == -1 {
610 return *self;
611 }
612
613 *self - Integer::mod_floor(self, other)
614 }
615 }
616
617 #[cfg(test)]
618 mod $test_mod {
619 use crate::Integer;
620 use core::mem;
621
622 #[cfg(test)]
628 fn test_division_rule((n, d): ($T, $T), (q, r): ($T, $T)) {
629 assert_eq!(d * q + r, n);
630 }
631
632 #[test]
633 fn test_div_rem() {
634 fn test_nd_dr(nd: ($T, $T), qr: ($T, $T)) {
635 let (n, d) = nd;
636 let separate_div_rem = (n / d, n % d);
637 let combined_div_rem = n.div_rem(&d);
638
639 test_division_rule(nd, qr);
640
641 assert_eq!(separate_div_rem, qr);
642 assert_eq!(combined_div_rem, qr);
643 }
644
645 test_nd_dr((8, 3), (2, 2));
646 test_nd_dr((8, -3), (-2, 2));
647 test_nd_dr((-8, 3), (-2, -2));
648 test_nd_dr((-8, -3), (2, -2));
649
650 test_nd_dr((1, 2), (0, 1));
651 test_nd_dr((1, -2), (0, 1));
652 test_nd_dr((-1, 2), (0, -1));
653 test_nd_dr((-1, -2), (0, -1));
654 }
655
656 #[test]
657 fn test_div_mod_floor() {
658 fn test_nd_dm(nd: ($T, $T), dm: ($T, $T)) {
659 let (n, d) = nd;
660 let separate_div_mod_floor =
661 (Integer::div_floor(&n, &d), Integer::mod_floor(&n, &d));
662 let combined_div_mod_floor = Integer::div_mod_floor(&n, &d);
663
664 test_division_rule(nd, dm);
665
666 assert_eq!(separate_div_mod_floor, dm);
667 assert_eq!(combined_div_mod_floor, dm);
668 }
669
670 test_nd_dm((8, 3), (2, 2));
671 test_nd_dm((8, -3), (-3, -1));
672 test_nd_dm((-8, 3), (-3, 1));
673 test_nd_dm((-8, -3), (2, -2));
674
675 test_nd_dm((1, 2), (0, 1));
676 test_nd_dm((1, -2), (-1, -1));
677 test_nd_dm((-1, 2), (-1, 1));
678 test_nd_dm((-1, -2), (0, -1));
679 }
680
681 #[test]
682 fn test_gcd() {
683 assert_eq!((10 as $T).gcd(&2), 2 as $T);
684 assert_eq!((10 as $T).gcd(&3), 1 as $T);
685 assert_eq!((0 as $T).gcd(&3), 3 as $T);
686 assert_eq!((3 as $T).gcd(&3), 3 as $T);
687 assert_eq!((56 as $T).gcd(&42), 14 as $T);
688 assert_eq!((3 as $T).gcd(&-3), 3 as $T);
689 assert_eq!((-6 as $T).gcd(&3), 3 as $T);
690 assert_eq!((-4 as $T).gcd(&-2), 2 as $T);
691 }
692
693 #[test]
694 fn test_gcd_cmp_with_euclidean() {
695 fn euclidean_gcd(mut m: $T, mut n: $T) -> $T {
696 while m != 0 {
697 mem::swap(&mut m, &mut n);
698 m %= n;
699 }
700
701 n.abs()
702 }
703
704 for i in -127..=127 {
707 for j in -127..=127 {
708 assert_eq!(euclidean_gcd(i, j), i.gcd(&j));
709 }
710 }
711 }
712
713 #[test]
714 fn test_gcd_min_val() {
715 let min = <$T>::min_value();
716 let max = <$T>::max_value();
717 let max_pow2 = max / 2 + 1;
718 assert_eq!(min.gcd(&max), 1 as $T);
719 assert_eq!(max.gcd(&min), 1 as $T);
720 assert_eq!(min.gcd(&max_pow2), max_pow2);
721 assert_eq!(max_pow2.gcd(&min), max_pow2);
722 assert_eq!(min.gcd(&42), 2 as $T);
723 assert_eq!((42 as $T).gcd(&min), 2 as $T);
724 }
725
726 #[test]
727 #[should_panic]
728 fn test_gcd_min_val_min_val() {
729 let min = <$T>::min_value();
730 assert!(min.gcd(&min) >= 0);
731 }
732
733 #[test]
734 #[should_panic]
735 fn test_gcd_min_val_0() {
736 let min = <$T>::min_value();
737 assert!(min.gcd(&0) >= 0);
738 }
739
740 #[test]
741 #[should_panic]
742 fn test_gcd_0_min_val() {
743 let min = <$T>::min_value();
744 assert!((0 as $T).gcd(&min) >= 0);
745 }
746
747 #[test]
748 fn test_lcm() {
749 assert_eq!((1 as $T).lcm(&0), 0 as $T);
750 assert_eq!((0 as $T).lcm(&1), 0 as $T);
751 assert_eq!((1 as $T).lcm(&1), 1 as $T);
752 assert_eq!((-1 as $T).lcm(&1), 1 as $T);
753 assert_eq!((1 as $T).lcm(&-1), 1 as $T);
754 assert_eq!((-1 as $T).lcm(&-1), 1 as $T);
755 assert_eq!((8 as $T).lcm(&9), 72 as $T);
756 assert_eq!((11 as $T).lcm(&5), 55 as $T);
757 }
758
759 #[test]
760 fn test_gcd_lcm() {
761 use core::iter::once;
762 for i in once(0)
763 .chain((1..).take(127).flat_map(|a| once(a).chain(once(-a))))
764 .chain(once(-128))
765 {
766 for j in once(0)
767 .chain((1..).take(127).flat_map(|a| once(a).chain(once(-a))))
768 .chain(once(-128))
769 {
770 assert_eq!(i.gcd_lcm(&j), (i.gcd(&j), i.lcm(&j)));
771 }
772 }
773 }
774
775 #[test]
776 fn test_extended_gcd_lcm() {
777 use crate::ExtendedGcd;
778 use core::fmt::Debug;
779 use num_traits::NumAssign;
780
781 fn check<A: Copy + Debug + Integer + NumAssign>(a: A, b: A) {
782 let ExtendedGcd { gcd, x, y, .. } = a.extended_gcd(&b);
783 assert_eq!(gcd, x * a + y * b);
784 }
785
786 use core::iter::once;
787 for i in once(0)
788 .chain((1..).take(127).flat_map(|a| once(a).chain(once(-a))))
789 .chain(once(-128))
790 {
791 for j in once(0)
792 .chain((1..).take(127).flat_map(|a| once(a).chain(once(-a))))
793 .chain(once(-128))
794 {
795 check(i, j);
796 let (ExtendedGcd { gcd, .. }, lcm) = i.extended_gcd_lcm(&j);
797 assert_eq!((gcd, lcm), (i.gcd(&j), i.lcm(&j)));
798 }
799 }
800 }
801
802 #[test]
803 fn test_even() {
804 assert_eq!((-4 as $T).is_even(), true);
805 assert_eq!((-3 as $T).is_even(), false);
806 assert_eq!((-2 as $T).is_even(), true);
807 assert_eq!((-1 as $T).is_even(), false);
808 assert_eq!((0 as $T).is_even(), true);
809 assert_eq!((1 as $T).is_even(), false);
810 assert_eq!((2 as $T).is_even(), true);
811 assert_eq!((3 as $T).is_even(), false);
812 assert_eq!((4 as $T).is_even(), true);
813 }
814
815 #[test]
816 fn test_odd() {
817 assert_eq!((-4 as $T).is_odd(), false);
818 assert_eq!((-3 as $T).is_odd(), true);
819 assert_eq!((-2 as $T).is_odd(), false);
820 assert_eq!((-1 as $T).is_odd(), true);
821 assert_eq!((0 as $T).is_odd(), false);
822 assert_eq!((1 as $T).is_odd(), true);
823 assert_eq!((2 as $T).is_odd(), false);
824 assert_eq!((3 as $T).is_odd(), true);
825 assert_eq!((4 as $T).is_odd(), false);
826 }
827
828 #[test]
829 fn test_multiple_of_one_limits() {
830 for x in &[<$T>::min_value(), <$T>::max_value()] {
831 for one in &[1, -1] {
832 assert_eq!(Integer::next_multiple_of(x, one), *x);
833 assert_eq!(Integer::prev_multiple_of(x, one), *x);
834 }
835 }
836 }
837 }
838 };
839}
840
841impl_integer_for_isize!(i8, test_integer_i8);
842impl_integer_for_isize!(i16, test_integer_i16);
843impl_integer_for_isize!(i32, test_integer_i32);
844impl_integer_for_isize!(i64, test_integer_i64);
845impl_integer_for_isize!(i128, test_integer_i128);
846impl_integer_for_isize!(isize, test_integer_isize);
847
848macro_rules! impl_integer_for_usize {
849 ($T:ty, $test_mod:ident) => {
850 impl Integer for $T {
851 #[inline]
853 fn div_floor(&self, other: &Self) -> Self {
854 *self / *other
855 }
856
857 #[inline]
859 fn mod_floor(&self, other: &Self) -> Self {
860 *self % *other
861 }
862
863 #[inline]
864 fn div_ceil(&self, other: &Self) -> Self {
865 *self / *other + (0 != *self % *other) as Self
866 }
867
868 #[inline]
870 fn gcd(&self, other: &Self) -> Self {
871 let mut m = *self;
873 let mut n = *other;
874 if m == 0 || n == 0 {
875 return m | n;
876 }
877
878 let shift = (m | n).trailing_zeros();
880
881 m >>= m.trailing_zeros();
883 n >>= n.trailing_zeros();
884
885 while m != n {
886 if m > n {
887 m -= n;
888 m >>= m.trailing_zeros();
889 } else {
890 n -= m;
891 n >>= n.trailing_zeros();
892 }
893 }
894 m << shift
895 }
896
897 #[inline]
898 fn extended_gcd_lcm(&self, other: &Self) -> (ExtendedGcd<Self>, Self) {
899 let egcd = self.extended_gcd(other);
900 let lcm = if egcd.gcd.is_zero() {
902 Self::zero()
903 } else {
904 *self * (*other / egcd.gcd)
905 };
906 (egcd, lcm)
907 }
908
909 #[inline]
911 fn lcm(&self, other: &Self) -> Self {
912 self.gcd_lcm(other).1
913 }
914
915 #[inline]
918 fn gcd_lcm(&self, other: &Self) -> (Self, Self) {
919 if self.is_zero() && other.is_zero() {
920 return (Self::zero(), Self::zero());
921 }
922 let gcd = self.gcd(other);
923 let lcm = *self * (*other / gcd);
924 (gcd, lcm)
925 }
926
927 #[inline]
929 fn is_multiple_of(&self, other: &Self) -> bool {
930 if other.is_zero() {
931 return self.is_zero();
932 }
933 *self % *other == 0
934 }
935
936 #[inline]
938 fn is_even(&self) -> bool {
939 *self % 2 == 0
940 }
941
942 #[inline]
944 fn is_odd(&self) -> bool {
945 !self.is_even()
946 }
947
948 #[inline]
950 fn div_rem(&self, other: &Self) -> (Self, Self) {
951 (*self / *other, *self % *other)
952 }
953 }
954
955 #[cfg(test)]
956 mod $test_mod {
957 use crate::Integer;
958 use core::mem;
959
960 #[test]
961 fn test_div_mod_floor() {
962 assert_eq!(<$T as Integer>::div_floor(&10, &3), 3 as $T);
963 assert_eq!(<$T as Integer>::mod_floor(&10, &3), 1 as $T);
964 assert_eq!(<$T as Integer>::div_mod_floor(&10, &3), (3 as $T, 1 as $T));
965 assert_eq!(<$T as Integer>::div_floor(&5, &5), 1 as $T);
966 assert_eq!(<$T as Integer>::mod_floor(&5, &5), 0 as $T);
967 assert_eq!(<$T as Integer>::div_mod_floor(&5, &5), (1 as $T, 0 as $T));
968 assert_eq!(<$T as Integer>::div_floor(&3, &7), 0 as $T);
969 assert_eq!(<$T as Integer>::div_floor(&3, &7), 0 as $T);
970 assert_eq!(<$T as Integer>::mod_floor(&3, &7), 3 as $T);
971 assert_eq!(<$T as Integer>::div_mod_floor(&3, &7), (0 as $T, 3 as $T));
972 }
973
974 #[test]
975 fn test_gcd() {
976 assert_eq!((10 as $T).gcd(&2), 2 as $T);
977 assert_eq!((10 as $T).gcd(&3), 1 as $T);
978 assert_eq!((0 as $T).gcd(&3), 3 as $T);
979 assert_eq!((3 as $T).gcd(&3), 3 as $T);
980 assert_eq!((56 as $T).gcd(&42), 14 as $T);
981 }
982
983 #[test]
984 fn test_gcd_cmp_with_euclidean() {
985 fn euclidean_gcd(mut m: $T, mut n: $T) -> $T {
986 while m != 0 {
987 mem::swap(&mut m, &mut n);
988 m %= n;
989 }
990 n
991 }
992
993 for i in 0..=255 {
994 for j in 0..=255 {
995 assert_eq!(euclidean_gcd(i, j), i.gcd(&j));
996 }
997 }
998 }
999
1000 #[test]
1001 fn test_lcm() {
1002 assert_eq!((1 as $T).lcm(&0), 0 as $T);
1003 assert_eq!((0 as $T).lcm(&1), 0 as $T);
1004 assert_eq!((1 as $T).lcm(&1), 1 as $T);
1005 assert_eq!((8 as $T).lcm(&9), 72 as $T);
1006 assert_eq!((11 as $T).lcm(&5), 55 as $T);
1007 assert_eq!((15 as $T).lcm(&17), 255 as $T);
1008 }
1009
1010 #[test]
1011 fn test_gcd_lcm() {
1012 for i in (0..).take(256) {
1013 for j in (0..).take(256) {
1014 assert_eq!(i.gcd_lcm(&j), (i.gcd(&j), i.lcm(&j)));
1015 }
1016 }
1017 }
1018
1019 #[test]
1020 fn test_is_multiple_of() {
1021 assert!((0 as $T).is_multiple_of(&(0 as $T)));
1022 assert!((6 as $T).is_multiple_of(&(6 as $T)));
1023 assert!((6 as $T).is_multiple_of(&(3 as $T)));
1024 assert!((6 as $T).is_multiple_of(&(1 as $T)));
1025
1026 assert!(!(42 as $T).is_multiple_of(&(5 as $T)));
1027 assert!(!(5 as $T).is_multiple_of(&(3 as $T)));
1028 assert!(!(42 as $T).is_multiple_of(&(0 as $T)));
1029 }
1030
1031 #[test]
1032 fn test_even() {
1033 assert_eq!((0 as $T).is_even(), true);
1034 assert_eq!((1 as $T).is_even(), false);
1035 assert_eq!((2 as $T).is_even(), true);
1036 assert_eq!((3 as $T).is_even(), false);
1037 assert_eq!((4 as $T).is_even(), true);
1038 }
1039
1040 #[test]
1041 fn test_odd() {
1042 assert_eq!((0 as $T).is_odd(), false);
1043 assert_eq!((1 as $T).is_odd(), true);
1044 assert_eq!((2 as $T).is_odd(), false);
1045 assert_eq!((3 as $T).is_odd(), true);
1046 assert_eq!((4 as $T).is_odd(), false);
1047 }
1048 }
1049 };
1050}
1051
1052impl_integer_for_usize!(u8, test_integer_u8);
1053impl_integer_for_usize!(u16, test_integer_u16);
1054impl_integer_for_usize!(u32, test_integer_u32);
1055impl_integer_for_usize!(u64, test_integer_u64);
1056impl_integer_for_usize!(u128, test_integer_u128);
1057impl_integer_for_usize!(usize, test_integer_usize);
1058
1059pub struct IterBinomial<T> {
1061 a: T,
1062 n: T,
1063 k: T,
1064}
1065
1066impl<T> IterBinomial<T>
1067where
1068 T: Integer,
1069{
1070 pub fn new(n: T) -> IterBinomial<T> {
1089 IterBinomial {
1090 k: T::zero(),
1091 a: T::one(),
1092 n,
1093 }
1094 }
1095}
1096
1097impl<T> Iterator for IterBinomial<T>
1098where
1099 T: Integer + Clone,
1100{
1101 type Item = T;
1102
1103 fn next(&mut self) -> Option<T> {
1104 if self.k > self.n {
1105 return None;
1106 }
1107 self.a = if !self.k.is_zero() {
1108 multiply_and_divide(
1109 self.a.clone(),
1110 self.n.clone() - self.k.clone() + T::one(),
1111 self.k.clone(),
1112 )
1113 } else {
1114 T::one()
1115 };
1116 self.k = self.k.clone() + T::one();
1117 Some(self.a.clone())
1118 }
1119}
1120
1121fn multiply_and_divide<T: Integer + Clone>(r: T, a: T, b: T) -> T {
1125 let g = gcd(r.clone(), b.clone());
1127 r / g.clone() * (a / (b / g))
1128}
1129
1130pub fn binomial<T: Integer + Clone>(mut n: T, k: T) -> T {
1149 if k > n {
1151 return T::zero();
1152 }
1153 if k > n.clone() - k.clone() {
1154 return binomial(n.clone(), n - k);
1155 }
1156 let mut r = T::one();
1157 let mut d = T::one();
1158 loop {
1159 if d > k {
1160 break;
1161 }
1162 r = multiply_and_divide(r, n.clone(), d.clone());
1163 n = n - T::one();
1164 d = d + T::one();
1165 }
1166 r
1167}
1168
1169pub fn multinomial<T: Integer + Clone>(k: &[T]) -> T
1171where
1172 for<'a> T: Add<&'a T, Output = T>,
1173{
1174 let mut r = T::one();
1175 let mut p = T::zero();
1176 for i in k {
1177 p = p + i;
1178 r = r * binomial(p.clone(), i.clone());
1179 }
1180 r
1181}
1182
1183#[test]
1184fn test_lcm_overflow() {
1185 macro_rules! check {
1186 ($t:ty, $x:expr, $y:expr, $r:expr) => {{
1187 let x: $t = $x;
1188 let y: $t = $y;
1189 let o = x.checked_mul(y);
1190 assert!(
1191 o.is_none(),
1192 "sanity checking that {} input {} * {} overflows",
1193 stringify!($t),
1194 x,
1195 y
1196 );
1197 assert_eq!(x.lcm(&y), $r);
1198 assert_eq!(y.lcm(&x), $r);
1199 }};
1200 }
1201
1202 check!(i64, 46656000000000000, 600, 46656000000000000);
1204
1205 check!(i8, 0x40, 0x04, 0x40);
1206 check!(u8, 0x80, 0x02, 0x80);
1207 check!(i16, 0x40_00, 0x04, 0x40_00);
1208 check!(u16, 0x80_00, 0x02, 0x80_00);
1209 check!(i32, 0x4000_0000, 0x04, 0x4000_0000);
1210 check!(u32, 0x8000_0000, 0x02, 0x8000_0000);
1211 check!(i64, 0x4000_0000_0000_0000, 0x04, 0x4000_0000_0000_0000);
1212 check!(u64, 0x8000_0000_0000_0000, 0x02, 0x8000_0000_0000_0000);
1213}
1214
1215#[test]
1216fn test_iter_binomial() {
1217 macro_rules! check_simple {
1218 ($t:ty) => {{
1219 let n: $t = 3;
1220 let expected = [1, 3, 3, 1];
1221 for (b, &e) in IterBinomial::new(n).zip(&expected) {
1222 assert_eq!(b, e);
1223 }
1224 }};
1225 }
1226
1227 check_simple!(u8);
1228 check_simple!(i8);
1229 check_simple!(u16);
1230 check_simple!(i16);
1231 check_simple!(u32);
1232 check_simple!(i32);
1233 check_simple!(u64);
1234 check_simple!(i64);
1235
1236 macro_rules! check_binomial {
1237 ($t:ty, $n:expr) => {{
1238 let n: $t = $n;
1239 let mut k: $t = 0;
1240 for b in IterBinomial::new(n) {
1241 assert_eq!(b, binomial(n, k));
1242 k += 1;
1243 }
1244 }};
1245 }
1246
1247 check_binomial!(u8, 10);
1249 check_binomial!(i8, 9);
1250 check_binomial!(u16, 18);
1251 check_binomial!(i16, 17);
1252 check_binomial!(u32, 34);
1253 check_binomial!(i32, 33);
1254 check_binomial!(u64, 67);
1255 check_binomial!(i64, 66);
1256}
1257
1258#[test]
1259fn test_binomial() {
1260 macro_rules! check {
1261 ($t:ty, $x:expr, $y:expr, $r:expr) => {{
1262 let x: $t = $x;
1263 let y: $t = $y;
1264 let expected: $t = $r;
1265 assert_eq!(binomial(x, y), expected);
1266 if y <= x {
1267 assert_eq!(binomial(x, x - y), expected);
1268 }
1269 }};
1270 }
1271 check!(u8, 9, 4, 126);
1272 check!(u8, 0, 0, 1);
1273 check!(u8, 2, 3, 0);
1274
1275 check!(i8, 9, 4, 126);
1276 check!(i8, 0, 0, 1);
1277 check!(i8, 2, 3, 0);
1278
1279 check!(u16, 100, 2, 4950);
1280 check!(u16, 14, 4, 1001);
1281 check!(u16, 0, 0, 1);
1282 check!(u16, 2, 3, 0);
1283
1284 check!(i16, 100, 2, 4950);
1285 check!(i16, 14, 4, 1001);
1286 check!(i16, 0, 0, 1);
1287 check!(i16, 2, 3, 0);
1288
1289 check!(u32, 100, 2, 4950);
1290 check!(u32, 35, 11, 417225900);
1291 check!(u32, 14, 4, 1001);
1292 check!(u32, 0, 0, 1);
1293 check!(u32, 2, 3, 0);
1294
1295 check!(i32, 100, 2, 4950);
1296 check!(i32, 35, 11, 417225900);
1297 check!(i32, 14, 4, 1001);
1298 check!(i32, 0, 0, 1);
1299 check!(i32, 2, 3, 0);
1300
1301 check!(u64, 100, 2, 4950);
1302 check!(u64, 35, 11, 417225900);
1303 check!(u64, 14, 4, 1001);
1304 check!(u64, 0, 0, 1);
1305 check!(u64, 2, 3, 0);
1306
1307 check!(i64, 100, 2, 4950);
1308 check!(i64, 35, 11, 417225900);
1309 check!(i64, 14, 4, 1001);
1310 check!(i64, 0, 0, 1);
1311 check!(i64, 2, 3, 0);
1312}
1313
1314#[test]
1315fn test_multinomial() {
1316 macro_rules! check_binomial {
1317 ($t:ty, $k:expr) => {{
1318 let n: $t = $k.iter().fold(0, |acc, &x| acc + x);
1319 let k: &[$t] = $k;
1320 assert_eq!(k.len(), 2);
1321 assert_eq!(multinomial(k), binomial(n, k[0]));
1322 }};
1323 }
1324
1325 check_binomial!(u8, &[4, 5]);
1326
1327 check_binomial!(i8, &[4, 5]);
1328
1329 check_binomial!(u16, &[2, 98]);
1330 check_binomial!(u16, &[4, 10]);
1331
1332 check_binomial!(i16, &[2, 98]);
1333 check_binomial!(i16, &[4, 10]);
1334
1335 check_binomial!(u32, &[2, 98]);
1336 check_binomial!(u32, &[11, 24]);
1337 check_binomial!(u32, &[4, 10]);
1338
1339 check_binomial!(i32, &[2, 98]);
1340 check_binomial!(i32, &[11, 24]);
1341 check_binomial!(i32, &[4, 10]);
1342
1343 check_binomial!(u64, &[2, 98]);
1344 check_binomial!(u64, &[11, 24]);
1345 check_binomial!(u64, &[4, 10]);
1346
1347 check_binomial!(i64, &[2, 98]);
1348 check_binomial!(i64, &[11, 24]);
1349 check_binomial!(i64, &[4, 10]);
1350
1351 macro_rules! check_multinomial {
1352 ($t:ty, $k:expr, $r:expr) => {{
1353 let k: &[$t] = $k;
1354 let expected: $t = $r;
1355 assert_eq!(multinomial(k), expected);
1356 }};
1357 }
1358
1359 check_multinomial!(u8, &[2, 1, 2], 30);
1360 check_multinomial!(u8, &[2, 3, 0], 10);
1361
1362 check_multinomial!(i8, &[2, 1, 2], 30);
1363 check_multinomial!(i8, &[2, 3, 0], 10);
1364
1365 check_multinomial!(u16, &[2, 1, 2], 30);
1366 check_multinomial!(u16, &[2, 3, 0], 10);
1367
1368 check_multinomial!(i16, &[2, 1, 2], 30);
1369 check_multinomial!(i16, &[2, 3, 0], 10);
1370
1371 check_multinomial!(u32, &[2, 1, 2], 30);
1372 check_multinomial!(u32, &[2, 3, 0], 10);
1373
1374 check_multinomial!(i32, &[2, 1, 2], 30);
1375 check_multinomial!(i32, &[2, 3, 0], 10);
1376
1377 check_multinomial!(u64, &[2, 1, 2], 30);
1378 check_multinomial!(u64, &[2, 3, 0], 10);
1379
1380 check_multinomial!(i64, &[2, 1, 2], 30);
1381 check_multinomial!(i64, &[2, 3, 0], 10);
1382
1383 check_multinomial!(u64, &[], 1);
1384 check_multinomial!(u64, &[0], 1);
1385 check_multinomial!(u64, &[12345], 1);
1386}