Yi Kong | c1a6315 | 2021-02-03 15:04:59 +0800 | [diff] [blame] | 1 | // Copyright 2013-2014 The Rust Project Developers. See the COPYRIGHT |
| 2 | // file at the top-level directory of this distribution and at |
| 3 | // http://rust-lang.org/COPYRIGHT. |
| 4 | // |
| 5 | // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or |
| 6 | // http://www.apache.org/licenses/LICENSE-2.0> or the MIT license |
| 7 | // <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your |
| 8 | // option. This file may not be copied, modified, or distributed |
| 9 | // except according to those terms. |
| 10 | |
| 11 | //! Integer trait and functions. |
| 12 | //! |
| 13 | //! ## Compatibility |
| 14 | //! |
| 15 | //! The `num-integer` crate is tested for rustc 1.8 and greater. |
| 16 | |
| 17 | #![doc(html_root_url = "https://docs.rs/num-integer/0.1")] |
| 18 | #![no_std] |
| 19 | #[cfg(feature = "std")] |
| 20 | extern crate std; |
| 21 | |
| 22 | extern crate num_traits as traits; |
| 23 | |
| 24 | use core::mem; |
| 25 | use core::ops::Add; |
| 26 | |
| 27 | use traits::{Num, Signed, Zero}; |
| 28 | |
| 29 | mod roots; |
| 30 | pub use roots::Roots; |
| 31 | pub use roots::{cbrt, nth_root, sqrt}; |
| 32 | |
| 33 | mod average; |
| 34 | pub use average::Average; |
| 35 | pub use average::{average_ceil, average_floor}; |
| 36 | |
| 37 | pub trait Integer: Sized + Num + PartialOrd + Ord + Eq { |
| 38 | /// Floored integer division. |
| 39 | /// |
| 40 | /// # Examples |
| 41 | /// |
| 42 | /// ~~~ |
| 43 | /// # use num_integer::Integer; |
| 44 | /// assert!(( 8).div_floor(& 3) == 2); |
| 45 | /// assert!(( 8).div_floor(&-3) == -3); |
| 46 | /// assert!((-8).div_floor(& 3) == -3); |
| 47 | /// assert!((-8).div_floor(&-3) == 2); |
| 48 | /// |
| 49 | /// assert!(( 1).div_floor(& 2) == 0); |
| 50 | /// assert!(( 1).div_floor(&-2) == -1); |
| 51 | /// assert!((-1).div_floor(& 2) == -1); |
| 52 | /// assert!((-1).div_floor(&-2) == 0); |
| 53 | /// ~~~ |
| 54 | fn div_floor(&self, other: &Self) -> Self; |
| 55 | |
| 56 | /// Floored integer modulo, satisfying: |
| 57 | /// |
| 58 | /// ~~~ |
| 59 | /// # use num_integer::Integer; |
| 60 | /// # let n = 1; let d = 1; |
| 61 | /// assert!(n.div_floor(&d) * d + n.mod_floor(&d) == n) |
| 62 | /// ~~~ |
| 63 | /// |
| 64 | /// # Examples |
| 65 | /// |
| 66 | /// ~~~ |
| 67 | /// # use num_integer::Integer; |
| 68 | /// assert!(( 8).mod_floor(& 3) == 2); |
| 69 | /// assert!(( 8).mod_floor(&-3) == -1); |
| 70 | /// assert!((-8).mod_floor(& 3) == 1); |
| 71 | /// assert!((-8).mod_floor(&-3) == -2); |
| 72 | /// |
| 73 | /// assert!(( 1).mod_floor(& 2) == 1); |
| 74 | /// assert!(( 1).mod_floor(&-2) == -1); |
| 75 | /// assert!((-1).mod_floor(& 2) == 1); |
| 76 | /// assert!((-1).mod_floor(&-2) == -1); |
| 77 | /// ~~~ |
| 78 | fn mod_floor(&self, other: &Self) -> Self; |
| 79 | |
| 80 | /// Ceiled integer division. |
| 81 | /// |
| 82 | /// # Examples |
| 83 | /// |
| 84 | /// ~~~ |
| 85 | /// # use num_integer::Integer; |
| 86 | /// assert_eq!(( 8).div_ceil( &3), 3); |
| 87 | /// assert_eq!(( 8).div_ceil(&-3), -2); |
| 88 | /// assert_eq!((-8).div_ceil( &3), -2); |
| 89 | /// assert_eq!((-8).div_ceil(&-3), 3); |
| 90 | /// |
| 91 | /// assert_eq!(( 1).div_ceil( &2), 1); |
| 92 | /// assert_eq!(( 1).div_ceil(&-2), 0); |
| 93 | /// assert_eq!((-1).div_ceil( &2), 0); |
| 94 | /// assert_eq!((-1).div_ceil(&-2), 1); |
| 95 | /// ~~~ |
| 96 | fn div_ceil(&self, other: &Self) -> Self { |
| 97 | let (q, r) = self.div_mod_floor(other); |
| 98 | if r.is_zero() { |
| 99 | q |
| 100 | } else { |
| 101 | q + Self::one() |
| 102 | } |
| 103 | } |
| 104 | |
| 105 | /// Greatest Common Divisor (GCD). |
| 106 | /// |
| 107 | /// # Examples |
| 108 | /// |
| 109 | /// ~~~ |
| 110 | /// # use num_integer::Integer; |
| 111 | /// assert_eq!(6.gcd(&8), 2); |
| 112 | /// assert_eq!(7.gcd(&3), 1); |
| 113 | /// ~~~ |
| 114 | fn gcd(&self, other: &Self) -> Self; |
| 115 | |
| 116 | /// Lowest Common Multiple (LCM). |
| 117 | /// |
| 118 | /// # Examples |
| 119 | /// |
| 120 | /// ~~~ |
| 121 | /// # use num_integer::Integer; |
| 122 | /// assert_eq!(7.lcm(&3), 21); |
| 123 | /// assert_eq!(2.lcm(&4), 4); |
| 124 | /// assert_eq!(0.lcm(&0), 0); |
| 125 | /// ~~~ |
| 126 | fn lcm(&self, other: &Self) -> Self; |
| 127 | |
| 128 | /// Greatest Common Divisor (GCD) and |
| 129 | /// Lowest Common Multiple (LCM) together. |
| 130 | /// |
| 131 | /// Potentially more efficient than calling `gcd` and `lcm` |
| 132 | /// individually for identical inputs. |
| 133 | /// |
| 134 | /// # Examples |
| 135 | /// |
| 136 | /// ~~~ |
| 137 | /// # use num_integer::Integer; |
| 138 | /// assert_eq!(10.gcd_lcm(&4), (2, 20)); |
| 139 | /// assert_eq!(8.gcd_lcm(&9), (1, 72)); |
| 140 | /// ~~~ |
| 141 | #[inline] |
| 142 | fn gcd_lcm(&self, other: &Self) -> (Self, Self) { |
| 143 | (self.gcd(other), self.lcm(other)) |
| 144 | } |
| 145 | |
| 146 | /// Greatest common divisor and Bézout coefficients. |
| 147 | /// |
| 148 | /// # Examples |
| 149 | /// |
| 150 | /// ~~~ |
| 151 | /// # extern crate num_integer; |
| 152 | /// # extern crate num_traits; |
| 153 | /// # fn main() { |
| 154 | /// # use num_integer::{ExtendedGcd, Integer}; |
| 155 | /// # use num_traits::NumAssign; |
| 156 | /// fn check<A: Copy + Integer + NumAssign>(a: A, b: A) -> bool { |
| 157 | /// let ExtendedGcd { gcd, x, y, .. } = a.extended_gcd(&b); |
| 158 | /// gcd == x * a + y * b |
| 159 | /// } |
| 160 | /// assert!(check(10isize, 4isize)); |
| 161 | /// assert!(check(8isize, 9isize)); |
| 162 | /// # } |
| 163 | /// ~~~ |
| 164 | #[inline] |
| 165 | fn extended_gcd(&self, other: &Self) -> ExtendedGcd<Self> |
| 166 | where |
| 167 | Self: Clone, |
| 168 | { |
| 169 | let mut s = (Self::zero(), Self::one()); |
| 170 | let mut t = (Self::one(), Self::zero()); |
| 171 | let mut r = (other.clone(), self.clone()); |
| 172 | |
| 173 | while !r.0.is_zero() { |
| 174 | let q = r.1.clone() / r.0.clone(); |
| 175 | let f = |mut r: (Self, Self)| { |
| 176 | mem::swap(&mut r.0, &mut r.1); |
| 177 | r.0 = r.0 - q.clone() * r.1.clone(); |
| 178 | r |
| 179 | }; |
| 180 | r = f(r); |
| 181 | s = f(s); |
| 182 | t = f(t); |
| 183 | } |
| 184 | |
| 185 | if r.1 >= Self::zero() { |
| 186 | ExtendedGcd { |
| 187 | gcd: r.1, |
| 188 | x: s.1, |
| 189 | y: t.1, |
| 190 | _hidden: (), |
| 191 | } |
| 192 | } else { |
| 193 | ExtendedGcd { |
| 194 | gcd: Self::zero() - r.1, |
| 195 | x: Self::zero() - s.1, |
| 196 | y: Self::zero() - t.1, |
| 197 | _hidden: (), |
| 198 | } |
| 199 | } |
| 200 | } |
| 201 | |
| 202 | /// Greatest common divisor, least common multiple, and Bézout coefficients. |
| 203 | #[inline] |
| 204 | fn extended_gcd_lcm(&self, other: &Self) -> (ExtendedGcd<Self>, Self) |
| 205 | where |
| 206 | Self: Clone + Signed, |
| 207 | { |
| 208 | (self.extended_gcd(other), self.lcm(other)) |
| 209 | } |
| 210 | |
| 211 | /// Deprecated, use `is_multiple_of` instead. |
| 212 | fn divides(&self, other: &Self) -> bool; |
| 213 | |
| 214 | /// Returns `true` if `self` is a multiple of `other`. |
| 215 | /// |
| 216 | /// # Examples |
| 217 | /// |
| 218 | /// ~~~ |
| 219 | /// # use num_integer::Integer; |
| 220 | /// assert_eq!(9.is_multiple_of(&3), true); |
| 221 | /// assert_eq!(3.is_multiple_of(&9), false); |
| 222 | /// ~~~ |
| 223 | fn is_multiple_of(&self, other: &Self) -> bool; |
| 224 | |
| 225 | /// Returns `true` if the number is even. |
| 226 | /// |
| 227 | /// # Examples |
| 228 | /// |
| 229 | /// ~~~ |
| 230 | /// # use num_integer::Integer; |
| 231 | /// assert_eq!(3.is_even(), false); |
| 232 | /// assert_eq!(4.is_even(), true); |
| 233 | /// ~~~ |
| 234 | fn is_even(&self) -> bool; |
| 235 | |
| 236 | /// Returns `true` if the number is odd. |
| 237 | /// |
| 238 | /// # Examples |
| 239 | /// |
| 240 | /// ~~~ |
| 241 | /// # use num_integer::Integer; |
| 242 | /// assert_eq!(3.is_odd(), true); |
| 243 | /// assert_eq!(4.is_odd(), false); |
| 244 | /// ~~~ |
| 245 | fn is_odd(&self) -> bool; |
| 246 | |
| 247 | /// Simultaneous truncated integer division and modulus. |
| 248 | /// Returns `(quotient, remainder)`. |
| 249 | /// |
| 250 | /// # Examples |
| 251 | /// |
| 252 | /// ~~~ |
| 253 | /// # use num_integer::Integer; |
| 254 | /// assert_eq!(( 8).div_rem( &3), ( 2, 2)); |
| 255 | /// assert_eq!(( 8).div_rem(&-3), (-2, 2)); |
| 256 | /// assert_eq!((-8).div_rem( &3), (-2, -2)); |
| 257 | /// assert_eq!((-8).div_rem(&-3), ( 2, -2)); |
| 258 | /// |
| 259 | /// assert_eq!(( 1).div_rem( &2), ( 0, 1)); |
| 260 | /// assert_eq!(( 1).div_rem(&-2), ( 0, 1)); |
| 261 | /// assert_eq!((-1).div_rem( &2), ( 0, -1)); |
| 262 | /// assert_eq!((-1).div_rem(&-2), ( 0, -1)); |
| 263 | /// ~~~ |
| 264 | fn div_rem(&self, other: &Self) -> (Self, Self); |
| 265 | |
| 266 | /// Simultaneous floored integer division and modulus. |
| 267 | /// Returns `(quotient, remainder)`. |
| 268 | /// |
| 269 | /// # Examples |
| 270 | /// |
| 271 | /// ~~~ |
| 272 | /// # use num_integer::Integer; |
| 273 | /// assert_eq!(( 8).div_mod_floor( &3), ( 2, 2)); |
| 274 | /// assert_eq!(( 8).div_mod_floor(&-3), (-3, -1)); |
| 275 | /// assert_eq!((-8).div_mod_floor( &3), (-3, 1)); |
| 276 | /// assert_eq!((-8).div_mod_floor(&-3), ( 2, -2)); |
| 277 | /// |
| 278 | /// assert_eq!(( 1).div_mod_floor( &2), ( 0, 1)); |
| 279 | /// assert_eq!(( 1).div_mod_floor(&-2), (-1, -1)); |
| 280 | /// assert_eq!((-1).div_mod_floor( &2), (-1, 1)); |
| 281 | /// assert_eq!((-1).div_mod_floor(&-2), ( 0, -1)); |
| 282 | /// ~~~ |
| 283 | fn div_mod_floor(&self, other: &Self) -> (Self, Self) { |
| 284 | (self.div_floor(other), self.mod_floor(other)) |
| 285 | } |
| 286 | |
| 287 | /// Rounds up to nearest multiple of argument. |
| 288 | /// |
| 289 | /// # Notes |
| 290 | /// |
| 291 | /// For signed types, `a.next_multiple_of(b) = a.prev_multiple_of(b.neg())`. |
| 292 | /// |
| 293 | /// # Examples |
| 294 | /// |
| 295 | /// ~~~ |
| 296 | /// # use num_integer::Integer; |
| 297 | /// assert_eq!(( 16).next_multiple_of(& 8), 16); |
| 298 | /// assert_eq!(( 23).next_multiple_of(& 8), 24); |
| 299 | /// assert_eq!(( 16).next_multiple_of(&-8), 16); |
| 300 | /// assert_eq!(( 23).next_multiple_of(&-8), 16); |
| 301 | /// assert_eq!((-16).next_multiple_of(& 8), -16); |
| 302 | /// assert_eq!((-23).next_multiple_of(& 8), -16); |
| 303 | /// assert_eq!((-16).next_multiple_of(&-8), -16); |
| 304 | /// assert_eq!((-23).next_multiple_of(&-8), -24); |
| 305 | /// ~~~ |
| 306 | #[inline] |
| 307 | fn next_multiple_of(&self, other: &Self) -> Self |
| 308 | where |
| 309 | Self: Clone, |
| 310 | { |
| 311 | let m = self.mod_floor(other); |
| 312 | self.clone() |
| 313 | + if m.is_zero() { |
| 314 | Self::zero() |
| 315 | } else { |
| 316 | other.clone() - m |
| 317 | } |
| 318 | } |
| 319 | |
| 320 | /// Rounds down to nearest multiple of argument. |
| 321 | /// |
| 322 | /// # Notes |
| 323 | /// |
| 324 | /// For signed types, `a.prev_multiple_of(b) = a.next_multiple_of(b.neg())`. |
| 325 | /// |
| 326 | /// # Examples |
| 327 | /// |
| 328 | /// ~~~ |
| 329 | /// # use num_integer::Integer; |
| 330 | /// assert_eq!(( 16).prev_multiple_of(& 8), 16); |
| 331 | /// assert_eq!(( 23).prev_multiple_of(& 8), 16); |
| 332 | /// assert_eq!(( 16).prev_multiple_of(&-8), 16); |
| 333 | /// assert_eq!(( 23).prev_multiple_of(&-8), 24); |
| 334 | /// assert_eq!((-16).prev_multiple_of(& 8), -16); |
| 335 | /// assert_eq!((-23).prev_multiple_of(& 8), -24); |
| 336 | /// assert_eq!((-16).prev_multiple_of(&-8), -16); |
| 337 | /// assert_eq!((-23).prev_multiple_of(&-8), -16); |
| 338 | /// ~~~ |
| 339 | #[inline] |
| 340 | fn prev_multiple_of(&self, other: &Self) -> Self |
| 341 | where |
| 342 | Self: Clone, |
| 343 | { |
| 344 | self.clone() - self.mod_floor(other) |
| 345 | } |
| 346 | } |
| 347 | |
| 348 | /// Greatest common divisor and Bézout coefficients |
| 349 | /// |
| 350 | /// ```no_build |
| 351 | /// let e = isize::extended_gcd(a, b); |
| 352 | /// assert_eq!(e.gcd, e.x*a + e.y*b); |
| 353 | /// ``` |
| 354 | #[derive(Debug, Clone, Copy, PartialEq, Eq)] |
| 355 | pub struct ExtendedGcd<A> { |
| 356 | pub gcd: A, |
| 357 | pub x: A, |
| 358 | pub y: A, |
| 359 | _hidden: (), |
| 360 | } |
| 361 | |
| 362 | /// Simultaneous integer division and modulus |
| 363 | #[inline] |
| 364 | pub fn div_rem<T: Integer>(x: T, y: T) -> (T, T) { |
| 365 | x.div_rem(&y) |
| 366 | } |
| 367 | /// Floored integer division |
| 368 | #[inline] |
| 369 | pub fn div_floor<T: Integer>(x: T, y: T) -> T { |
| 370 | x.div_floor(&y) |
| 371 | } |
| 372 | /// Floored integer modulus |
| 373 | #[inline] |
| 374 | pub fn mod_floor<T: Integer>(x: T, y: T) -> T { |
| 375 | x.mod_floor(&y) |
| 376 | } |
| 377 | /// Simultaneous floored integer division and modulus |
| 378 | #[inline] |
| 379 | pub fn div_mod_floor<T: Integer>(x: T, y: T) -> (T, T) { |
| 380 | x.div_mod_floor(&y) |
| 381 | } |
| 382 | /// Ceiled integer division |
| 383 | #[inline] |
| 384 | pub fn div_ceil<T: Integer>(x: T, y: T) -> T { |
| 385 | x.div_ceil(&y) |
| 386 | } |
| 387 | |
| 388 | /// Calculates the Greatest Common Divisor (GCD) of the number and `other`. The |
| 389 | /// result is always positive. |
| 390 | #[inline(always)] |
| 391 | pub fn gcd<T: Integer>(x: T, y: T) -> T { |
| 392 | x.gcd(&y) |
| 393 | } |
| 394 | /// Calculates the Lowest Common Multiple (LCM) of the number and `other`. |
| 395 | #[inline(always)] |
| 396 | pub fn lcm<T: Integer>(x: T, y: T) -> T { |
| 397 | x.lcm(&y) |
| 398 | } |
| 399 | |
| 400 | /// Calculates the Greatest Common Divisor (GCD) and |
| 401 | /// Lowest Common Multiple (LCM) of the number and `other`. |
| 402 | #[inline(always)] |
| 403 | pub fn gcd_lcm<T: Integer>(x: T, y: T) -> (T, T) { |
| 404 | x.gcd_lcm(&y) |
| 405 | } |
| 406 | |
| 407 | macro_rules! impl_integer_for_isize { |
| 408 | ($T:ty, $test_mod:ident) => { |
| 409 | impl Integer for $T { |
| 410 | /// Floored integer division |
| 411 | #[inline] |
| 412 | fn div_floor(&self, other: &Self) -> Self { |
| 413 | // Algorithm from [Daan Leijen. _Division and Modulus for Computer Scientists_, |
| 414 | // December 2001](http://research.microsoft.com/pubs/151917/divmodnote-letter.pdf) |
| 415 | let (d, r) = self.div_rem(other); |
| 416 | if (r > 0 && *other < 0) || (r < 0 && *other > 0) { |
| 417 | d - 1 |
| 418 | } else { |
| 419 | d |
| 420 | } |
| 421 | } |
| 422 | |
| 423 | /// Floored integer modulo |
| 424 | #[inline] |
| 425 | fn mod_floor(&self, other: &Self) -> Self { |
| 426 | // Algorithm from [Daan Leijen. _Division and Modulus for Computer Scientists_, |
| 427 | // December 2001](http://research.microsoft.com/pubs/151917/divmodnote-letter.pdf) |
| 428 | let r = *self % *other; |
| 429 | if (r > 0 && *other < 0) || (r < 0 && *other > 0) { |
| 430 | r + *other |
| 431 | } else { |
| 432 | r |
| 433 | } |
| 434 | } |
| 435 | |
| 436 | /// Calculates `div_floor` and `mod_floor` simultaneously |
| 437 | #[inline] |
| 438 | fn div_mod_floor(&self, other: &Self) -> (Self, Self) { |
| 439 | // Algorithm from [Daan Leijen. _Division and Modulus for Computer Scientists_, |
| 440 | // December 2001](http://research.microsoft.com/pubs/151917/divmodnote-letter.pdf) |
| 441 | let (d, r) = self.div_rem(other); |
| 442 | if (r > 0 && *other < 0) || (r < 0 && *other > 0) { |
| 443 | (d - 1, r + *other) |
| 444 | } else { |
| 445 | (d, r) |
| 446 | } |
| 447 | } |
| 448 | |
| 449 | #[inline] |
| 450 | fn div_ceil(&self, other: &Self) -> Self { |
| 451 | let (d, r) = self.div_rem(other); |
| 452 | if (r > 0 && *other > 0) || (r < 0 && *other < 0) { |
| 453 | d + 1 |
| 454 | } else { |
| 455 | d |
| 456 | } |
| 457 | } |
| 458 | |
| 459 | /// Calculates the Greatest Common Divisor (GCD) of the number and |
| 460 | /// `other`. The result is always positive. |
| 461 | #[inline] |
| 462 | fn gcd(&self, other: &Self) -> Self { |
| 463 | // Use Stein's algorithm |
| 464 | let mut m = *self; |
| 465 | let mut n = *other; |
| 466 | if m == 0 || n == 0 { |
| 467 | return (m | n).abs(); |
| 468 | } |
| 469 | |
| 470 | // find common factors of 2 |
| 471 | let shift = (m | n).trailing_zeros(); |
| 472 | |
| 473 | // The algorithm needs positive numbers, but the minimum value |
| 474 | // can't be represented as a positive one. |
| 475 | // It's also a power of two, so the gcd can be |
| 476 | // calculated by bitshifting in that case |
| 477 | |
| 478 | // Assuming two's complement, the number created by the shift |
| 479 | // is positive for all numbers except gcd = abs(min value) |
| 480 | // The call to .abs() causes a panic in debug mode |
| 481 | if m == Self::min_value() || n == Self::min_value() { |
| 482 | return (1 << shift).abs(); |
| 483 | } |
| 484 | |
| 485 | // guaranteed to be positive now, rest like unsigned algorithm |
| 486 | m = m.abs(); |
| 487 | n = n.abs(); |
| 488 | |
| 489 | // divide n and m by 2 until odd |
| 490 | m >>= m.trailing_zeros(); |
| 491 | n >>= n.trailing_zeros(); |
| 492 | |
| 493 | while m != n { |
| 494 | if m > n { |
| 495 | m -= n; |
| 496 | m >>= m.trailing_zeros(); |
| 497 | } else { |
| 498 | n -= m; |
| 499 | n >>= n.trailing_zeros(); |
| 500 | } |
| 501 | } |
| 502 | m << shift |
| 503 | } |
| 504 | |
| 505 | #[inline] |
| 506 | fn extended_gcd_lcm(&self, other: &Self) -> (ExtendedGcd<Self>, Self) { |
| 507 | let egcd = self.extended_gcd(other); |
| 508 | // should not have to recalculate abs |
| 509 | let lcm = if egcd.gcd.is_zero() { |
| 510 | Self::zero() |
| 511 | } else { |
| 512 | (*self * (*other / egcd.gcd)).abs() |
| 513 | }; |
| 514 | (egcd, lcm) |
| 515 | } |
| 516 | |
| 517 | /// Calculates the Lowest Common Multiple (LCM) of the number and |
| 518 | /// `other`. |
| 519 | #[inline] |
| 520 | fn lcm(&self, other: &Self) -> Self { |
| 521 | self.gcd_lcm(other).1 |
| 522 | } |
| 523 | |
| 524 | /// Calculates the Greatest Common Divisor (GCD) and |
| 525 | /// Lowest Common Multiple (LCM) of the number and `other`. |
| 526 | #[inline] |
| 527 | fn gcd_lcm(&self, other: &Self) -> (Self, Self) { |
| 528 | if self.is_zero() && other.is_zero() { |
| 529 | return (Self::zero(), Self::zero()); |
| 530 | } |
| 531 | let gcd = self.gcd(other); |
| 532 | // should not have to recalculate abs |
| 533 | let lcm = (*self * (*other / gcd)).abs(); |
| 534 | (gcd, lcm) |
| 535 | } |
| 536 | |
| 537 | /// Deprecated, use `is_multiple_of` instead. |
| 538 | #[inline] |
| 539 | fn divides(&self, other: &Self) -> bool { |
| 540 | self.is_multiple_of(other) |
| 541 | } |
| 542 | |
| 543 | /// Returns `true` if the number is a multiple of `other`. |
| 544 | #[inline] |
| 545 | fn is_multiple_of(&self, other: &Self) -> bool { |
| 546 | *self % *other == 0 |
| 547 | } |
| 548 | |
| 549 | /// Returns `true` if the number is divisible by `2` |
| 550 | #[inline] |
| 551 | fn is_even(&self) -> bool { |
| 552 | (*self) & 1 == 0 |
| 553 | } |
| 554 | |
| 555 | /// Returns `true` if the number is not divisible by `2` |
| 556 | #[inline] |
| 557 | fn is_odd(&self) -> bool { |
| 558 | !self.is_even() |
| 559 | } |
| 560 | |
| 561 | /// Simultaneous truncated integer division and modulus. |
| 562 | #[inline] |
| 563 | fn div_rem(&self, other: &Self) -> (Self, Self) { |
| 564 | (*self / *other, *self % *other) |
| 565 | } |
| 566 | } |
| 567 | |
| 568 | #[cfg(test)] |
| 569 | mod $test_mod { |
| 570 | use core::mem; |
| 571 | use Integer; |
| 572 | |
| 573 | /// Checks that the division rule holds for: |
| 574 | /// |
| 575 | /// - `n`: numerator (dividend) |
| 576 | /// - `d`: denominator (divisor) |
| 577 | /// - `qr`: quotient and remainder |
| 578 | #[cfg(test)] |
| 579 | fn test_division_rule((n, d): ($T, $T), (q, r): ($T, $T)) { |
| 580 | assert_eq!(d * q + r, n); |
| 581 | } |
| 582 | |
| 583 | #[test] |
| 584 | fn test_div_rem() { |
| 585 | fn test_nd_dr(nd: ($T, $T), qr: ($T, $T)) { |
| 586 | let (n, d) = nd; |
| 587 | let separate_div_rem = (n / d, n % d); |
| 588 | let combined_div_rem = n.div_rem(&d); |
| 589 | |
| 590 | assert_eq!(separate_div_rem, qr); |
| 591 | assert_eq!(combined_div_rem, qr); |
| 592 | |
| 593 | test_division_rule(nd, separate_div_rem); |
| 594 | test_division_rule(nd, combined_div_rem); |
| 595 | } |
| 596 | |
| 597 | test_nd_dr((8, 3), (2, 2)); |
| 598 | test_nd_dr((8, -3), (-2, 2)); |
| 599 | test_nd_dr((-8, 3), (-2, -2)); |
| 600 | test_nd_dr((-8, -3), (2, -2)); |
| 601 | |
| 602 | test_nd_dr((1, 2), (0, 1)); |
| 603 | test_nd_dr((1, -2), (0, 1)); |
| 604 | test_nd_dr((-1, 2), (0, -1)); |
| 605 | test_nd_dr((-1, -2), (0, -1)); |
| 606 | } |
| 607 | |
| 608 | #[test] |
| 609 | fn test_div_mod_floor() { |
| 610 | fn test_nd_dm(nd: ($T, $T), dm: ($T, $T)) { |
| 611 | let (n, d) = nd; |
| 612 | let separate_div_mod_floor = (n.div_floor(&d), n.mod_floor(&d)); |
| 613 | let combined_div_mod_floor = n.div_mod_floor(&d); |
| 614 | |
| 615 | assert_eq!(separate_div_mod_floor, dm); |
| 616 | assert_eq!(combined_div_mod_floor, dm); |
| 617 | |
| 618 | test_division_rule(nd, separate_div_mod_floor); |
| 619 | test_division_rule(nd, combined_div_mod_floor); |
| 620 | } |
| 621 | |
| 622 | test_nd_dm((8, 3), (2, 2)); |
| 623 | test_nd_dm((8, -3), (-3, -1)); |
| 624 | test_nd_dm((-8, 3), (-3, 1)); |
| 625 | test_nd_dm((-8, -3), (2, -2)); |
| 626 | |
| 627 | test_nd_dm((1, 2), (0, 1)); |
| 628 | test_nd_dm((1, -2), (-1, -1)); |
| 629 | test_nd_dm((-1, 2), (-1, 1)); |
| 630 | test_nd_dm((-1, -2), (0, -1)); |
| 631 | } |
| 632 | |
| 633 | #[test] |
| 634 | fn test_gcd() { |
| 635 | assert_eq!((10 as $T).gcd(&2), 2 as $T); |
| 636 | assert_eq!((10 as $T).gcd(&3), 1 as $T); |
| 637 | assert_eq!((0 as $T).gcd(&3), 3 as $T); |
| 638 | assert_eq!((3 as $T).gcd(&3), 3 as $T); |
| 639 | assert_eq!((56 as $T).gcd(&42), 14 as $T); |
| 640 | assert_eq!((3 as $T).gcd(&-3), 3 as $T); |
| 641 | assert_eq!((-6 as $T).gcd(&3), 3 as $T); |
| 642 | assert_eq!((-4 as $T).gcd(&-2), 2 as $T); |
| 643 | } |
| 644 | |
| 645 | #[test] |
| 646 | fn test_gcd_cmp_with_euclidean() { |
| 647 | fn euclidean_gcd(mut m: $T, mut n: $T) -> $T { |
| 648 | while m != 0 { |
| 649 | mem::swap(&mut m, &mut n); |
| 650 | m %= n; |
| 651 | } |
| 652 | |
| 653 | n.abs() |
| 654 | } |
| 655 | |
| 656 | // gcd(-128, b) = 128 is not representable as positive value |
| 657 | // for i8 |
| 658 | for i in -127..127 { |
| 659 | for j in -127..127 { |
| 660 | assert_eq!(euclidean_gcd(i, j), i.gcd(&j)); |
| 661 | } |
| 662 | } |
| 663 | |
| 664 | // last value |
| 665 | // FIXME: Use inclusive ranges for above loop when implemented |
| 666 | let i = 127; |
| 667 | for j in -127..127 { |
| 668 | assert_eq!(euclidean_gcd(i, j), i.gcd(&j)); |
| 669 | } |
| 670 | assert_eq!(127.gcd(&127), 127); |
| 671 | } |
| 672 | |
| 673 | #[test] |
| 674 | fn test_gcd_min_val() { |
| 675 | let min = <$T>::min_value(); |
| 676 | let max = <$T>::max_value(); |
| 677 | let max_pow2 = max / 2 + 1; |
| 678 | assert_eq!(min.gcd(&max), 1 as $T); |
| 679 | assert_eq!(max.gcd(&min), 1 as $T); |
| 680 | assert_eq!(min.gcd(&max_pow2), max_pow2); |
| 681 | assert_eq!(max_pow2.gcd(&min), max_pow2); |
| 682 | assert_eq!(min.gcd(&42), 2 as $T); |
| 683 | assert_eq!((42 as $T).gcd(&min), 2 as $T); |
| 684 | } |
| 685 | |
| 686 | #[test] |
| 687 | #[should_panic] |
| 688 | fn test_gcd_min_val_min_val() { |
| 689 | let min = <$T>::min_value(); |
| 690 | assert!(min.gcd(&min) >= 0); |
| 691 | } |
| 692 | |
| 693 | #[test] |
| 694 | #[should_panic] |
| 695 | fn test_gcd_min_val_0() { |
| 696 | let min = <$T>::min_value(); |
| 697 | assert!(min.gcd(&0) >= 0); |
| 698 | } |
| 699 | |
| 700 | #[test] |
| 701 | #[should_panic] |
| 702 | fn test_gcd_0_min_val() { |
| 703 | let min = <$T>::min_value(); |
| 704 | assert!((0 as $T).gcd(&min) >= 0); |
| 705 | } |
| 706 | |
| 707 | #[test] |
| 708 | fn test_lcm() { |
| 709 | assert_eq!((1 as $T).lcm(&0), 0 as $T); |
| 710 | assert_eq!((0 as $T).lcm(&1), 0 as $T); |
| 711 | assert_eq!((1 as $T).lcm(&1), 1 as $T); |
| 712 | assert_eq!((-1 as $T).lcm(&1), 1 as $T); |
| 713 | assert_eq!((1 as $T).lcm(&-1), 1 as $T); |
| 714 | assert_eq!((-1 as $T).lcm(&-1), 1 as $T); |
| 715 | assert_eq!((8 as $T).lcm(&9), 72 as $T); |
| 716 | assert_eq!((11 as $T).lcm(&5), 55 as $T); |
| 717 | } |
| 718 | |
| 719 | #[test] |
| 720 | fn test_gcd_lcm() { |
| 721 | use core::iter::once; |
| 722 | for i in once(0) |
| 723 | .chain((1..).take(127).flat_map(|a| once(a).chain(once(-a)))) |
| 724 | .chain(once(-128)) |
| 725 | { |
| 726 | for j in once(0) |
| 727 | .chain((1..).take(127).flat_map(|a| once(a).chain(once(-a)))) |
| 728 | .chain(once(-128)) |
| 729 | { |
| 730 | assert_eq!(i.gcd_lcm(&j), (i.gcd(&j), i.lcm(&j))); |
| 731 | } |
| 732 | } |
| 733 | } |
| 734 | |
| 735 | #[test] |
| 736 | fn test_extended_gcd_lcm() { |
| 737 | use core::fmt::Debug; |
| 738 | use traits::NumAssign; |
| 739 | use ExtendedGcd; |
| 740 | |
| 741 | fn check<A: Copy + Debug + Integer + NumAssign>(a: A, b: A) { |
| 742 | let ExtendedGcd { gcd, x, y, .. } = a.extended_gcd(&b); |
| 743 | assert_eq!(gcd, x * a + y * b); |
| 744 | } |
| 745 | |
| 746 | use core::iter::once; |
| 747 | for i in once(0) |
| 748 | .chain((1..).take(127).flat_map(|a| once(a).chain(once(-a)))) |
| 749 | .chain(once(-128)) |
| 750 | { |
| 751 | for j in once(0) |
| 752 | .chain((1..).take(127).flat_map(|a| once(a).chain(once(-a)))) |
| 753 | .chain(once(-128)) |
| 754 | { |
| 755 | check(i, j); |
| 756 | let (ExtendedGcd { gcd, .. }, lcm) = i.extended_gcd_lcm(&j); |
| 757 | assert_eq!((gcd, lcm), (i.gcd(&j), i.lcm(&j))); |
| 758 | } |
| 759 | } |
| 760 | } |
| 761 | |
| 762 | #[test] |
| 763 | fn test_even() { |
| 764 | assert_eq!((-4 as $T).is_even(), true); |
| 765 | assert_eq!((-3 as $T).is_even(), false); |
| 766 | assert_eq!((-2 as $T).is_even(), true); |
| 767 | assert_eq!((-1 as $T).is_even(), false); |
| 768 | assert_eq!((0 as $T).is_even(), true); |
| 769 | assert_eq!((1 as $T).is_even(), false); |
| 770 | assert_eq!((2 as $T).is_even(), true); |
| 771 | assert_eq!((3 as $T).is_even(), false); |
| 772 | assert_eq!((4 as $T).is_even(), true); |
| 773 | } |
| 774 | |
| 775 | #[test] |
| 776 | fn test_odd() { |
| 777 | assert_eq!((-4 as $T).is_odd(), false); |
| 778 | assert_eq!((-3 as $T).is_odd(), true); |
| 779 | assert_eq!((-2 as $T).is_odd(), false); |
| 780 | assert_eq!((-1 as $T).is_odd(), true); |
| 781 | assert_eq!((0 as $T).is_odd(), false); |
| 782 | assert_eq!((1 as $T).is_odd(), true); |
| 783 | assert_eq!((2 as $T).is_odd(), false); |
| 784 | assert_eq!((3 as $T).is_odd(), true); |
| 785 | assert_eq!((4 as $T).is_odd(), false); |
| 786 | } |
| 787 | } |
| 788 | }; |
| 789 | } |
| 790 | |
| 791 | impl_integer_for_isize!(i8, test_integer_i8); |
| 792 | impl_integer_for_isize!(i16, test_integer_i16); |
| 793 | impl_integer_for_isize!(i32, test_integer_i32); |
| 794 | impl_integer_for_isize!(i64, test_integer_i64); |
| 795 | impl_integer_for_isize!(isize, test_integer_isize); |
| 796 | #[cfg(has_i128)] |
| 797 | impl_integer_for_isize!(i128, test_integer_i128); |
| 798 | |
| 799 | macro_rules! impl_integer_for_usize { |
| 800 | ($T:ty, $test_mod:ident) => { |
| 801 | impl Integer for $T { |
| 802 | /// Unsigned integer division. Returns the same result as `div` (`/`). |
| 803 | #[inline] |
| 804 | fn div_floor(&self, other: &Self) -> Self { |
| 805 | *self / *other |
| 806 | } |
| 807 | |
| 808 | /// Unsigned integer modulo operation. Returns the same result as `rem` (`%`). |
| 809 | #[inline] |
| 810 | fn mod_floor(&self, other: &Self) -> Self { |
| 811 | *self % *other |
| 812 | } |
| 813 | |
| 814 | #[inline] |
| 815 | fn div_ceil(&self, other: &Self) -> Self { |
| 816 | *self / *other + (0 != *self % *other) as Self |
| 817 | } |
| 818 | |
| 819 | /// Calculates the Greatest Common Divisor (GCD) of the number and `other` |
| 820 | #[inline] |
| 821 | fn gcd(&self, other: &Self) -> Self { |
| 822 | // Use Stein's algorithm |
| 823 | let mut m = *self; |
| 824 | let mut n = *other; |
| 825 | if m == 0 || n == 0 { |
| 826 | return m | n; |
| 827 | } |
| 828 | |
| 829 | // find common factors of 2 |
| 830 | let shift = (m | n).trailing_zeros(); |
| 831 | |
| 832 | // divide n and m by 2 until odd |
| 833 | m >>= m.trailing_zeros(); |
| 834 | n >>= n.trailing_zeros(); |
| 835 | |
| 836 | while m != n { |
| 837 | if m > n { |
| 838 | m -= n; |
| 839 | m >>= m.trailing_zeros(); |
| 840 | } else { |
| 841 | n -= m; |
| 842 | n >>= n.trailing_zeros(); |
| 843 | } |
| 844 | } |
| 845 | m << shift |
| 846 | } |
| 847 | |
| 848 | #[inline] |
| 849 | fn extended_gcd_lcm(&self, other: &Self) -> (ExtendedGcd<Self>, Self) { |
| 850 | let egcd = self.extended_gcd(other); |
| 851 | // should not have to recalculate abs |
| 852 | let lcm = if egcd.gcd.is_zero() { |
| 853 | Self::zero() |
| 854 | } else { |
| 855 | *self * (*other / egcd.gcd) |
| 856 | }; |
| 857 | (egcd, lcm) |
| 858 | } |
| 859 | |
| 860 | /// Calculates the Lowest Common Multiple (LCM) of the number and `other`. |
| 861 | #[inline] |
| 862 | fn lcm(&self, other: &Self) -> Self { |
| 863 | self.gcd_lcm(other).1 |
| 864 | } |
| 865 | |
| 866 | /// Calculates the Greatest Common Divisor (GCD) and |
| 867 | /// Lowest Common Multiple (LCM) of the number and `other`. |
| 868 | #[inline] |
| 869 | fn gcd_lcm(&self, other: &Self) -> (Self, Self) { |
| 870 | if self.is_zero() && other.is_zero() { |
| 871 | return (Self::zero(), Self::zero()); |
| 872 | } |
| 873 | let gcd = self.gcd(other); |
| 874 | let lcm = *self * (*other / gcd); |
| 875 | (gcd, lcm) |
| 876 | } |
| 877 | |
| 878 | /// Deprecated, use `is_multiple_of` instead. |
| 879 | #[inline] |
| 880 | fn divides(&self, other: &Self) -> bool { |
| 881 | self.is_multiple_of(other) |
| 882 | } |
| 883 | |
| 884 | /// Returns `true` if the number is a multiple of `other`. |
| 885 | #[inline] |
| 886 | fn is_multiple_of(&self, other: &Self) -> bool { |
| 887 | *self % *other == 0 |
| 888 | } |
| 889 | |
| 890 | /// Returns `true` if the number is divisible by `2`. |
| 891 | #[inline] |
| 892 | fn is_even(&self) -> bool { |
| 893 | *self % 2 == 0 |
| 894 | } |
| 895 | |
| 896 | /// Returns `true` if the number is not divisible by `2`. |
| 897 | #[inline] |
| 898 | fn is_odd(&self) -> bool { |
| 899 | !self.is_even() |
| 900 | } |
| 901 | |
| 902 | /// Simultaneous truncated integer division and modulus. |
| 903 | #[inline] |
| 904 | fn div_rem(&self, other: &Self) -> (Self, Self) { |
| 905 | (*self / *other, *self % *other) |
| 906 | } |
| 907 | } |
| 908 | |
| 909 | #[cfg(test)] |
| 910 | mod $test_mod { |
| 911 | use core::mem; |
| 912 | use Integer; |
| 913 | |
| 914 | #[test] |
| 915 | fn test_div_mod_floor() { |
| 916 | assert_eq!((10 as $T).div_floor(&(3 as $T)), 3 as $T); |
| 917 | assert_eq!((10 as $T).mod_floor(&(3 as $T)), 1 as $T); |
| 918 | assert_eq!((10 as $T).div_mod_floor(&(3 as $T)), (3 as $T, 1 as $T)); |
| 919 | assert_eq!((5 as $T).div_floor(&(5 as $T)), 1 as $T); |
| 920 | assert_eq!((5 as $T).mod_floor(&(5 as $T)), 0 as $T); |
| 921 | assert_eq!((5 as $T).div_mod_floor(&(5 as $T)), (1 as $T, 0 as $T)); |
| 922 | assert_eq!((3 as $T).div_floor(&(7 as $T)), 0 as $T); |
| 923 | assert_eq!((3 as $T).mod_floor(&(7 as $T)), 3 as $T); |
| 924 | assert_eq!((3 as $T).div_mod_floor(&(7 as $T)), (0 as $T, 3 as $T)); |
| 925 | } |
| 926 | |
| 927 | #[test] |
| 928 | fn test_gcd() { |
| 929 | assert_eq!((10 as $T).gcd(&2), 2 as $T); |
| 930 | assert_eq!((10 as $T).gcd(&3), 1 as $T); |
| 931 | assert_eq!((0 as $T).gcd(&3), 3 as $T); |
| 932 | assert_eq!((3 as $T).gcd(&3), 3 as $T); |
| 933 | assert_eq!((56 as $T).gcd(&42), 14 as $T); |
| 934 | } |
| 935 | |
| 936 | #[test] |
| 937 | fn test_gcd_cmp_with_euclidean() { |
| 938 | fn euclidean_gcd(mut m: $T, mut n: $T) -> $T { |
| 939 | while m != 0 { |
| 940 | mem::swap(&mut m, &mut n); |
| 941 | m %= n; |
| 942 | } |
| 943 | n |
| 944 | } |
| 945 | |
| 946 | for i in 0..255 { |
| 947 | for j in 0..255 { |
| 948 | assert_eq!(euclidean_gcd(i, j), i.gcd(&j)); |
| 949 | } |
| 950 | } |
| 951 | |
| 952 | // last value |
| 953 | // FIXME: Use inclusive ranges for above loop when implemented |
| 954 | let i = 255; |
| 955 | for j in 0..255 { |
| 956 | assert_eq!(euclidean_gcd(i, j), i.gcd(&j)); |
| 957 | } |
| 958 | assert_eq!(255.gcd(&255), 255); |
| 959 | } |
| 960 | |
| 961 | #[test] |
| 962 | fn test_lcm() { |
| 963 | assert_eq!((1 as $T).lcm(&0), 0 as $T); |
| 964 | assert_eq!((0 as $T).lcm(&1), 0 as $T); |
| 965 | assert_eq!((1 as $T).lcm(&1), 1 as $T); |
| 966 | assert_eq!((8 as $T).lcm(&9), 72 as $T); |
| 967 | assert_eq!((11 as $T).lcm(&5), 55 as $T); |
| 968 | assert_eq!((15 as $T).lcm(&17), 255 as $T); |
| 969 | } |
| 970 | |
| 971 | #[test] |
| 972 | fn test_gcd_lcm() { |
| 973 | for i in (0..).take(256) { |
| 974 | for j in (0..).take(256) { |
| 975 | assert_eq!(i.gcd_lcm(&j), (i.gcd(&j), i.lcm(&j))); |
| 976 | } |
| 977 | } |
| 978 | } |
| 979 | |
| 980 | #[test] |
| 981 | fn test_is_multiple_of() { |
| 982 | assert!((6 as $T).is_multiple_of(&(6 as $T))); |
| 983 | assert!((6 as $T).is_multiple_of(&(3 as $T))); |
| 984 | assert!((6 as $T).is_multiple_of(&(1 as $T))); |
| 985 | } |
| 986 | |
| 987 | #[test] |
| 988 | fn test_even() { |
| 989 | assert_eq!((0 as $T).is_even(), true); |
| 990 | assert_eq!((1 as $T).is_even(), false); |
| 991 | assert_eq!((2 as $T).is_even(), true); |
| 992 | assert_eq!((3 as $T).is_even(), false); |
| 993 | assert_eq!((4 as $T).is_even(), true); |
| 994 | } |
| 995 | |
| 996 | #[test] |
| 997 | fn test_odd() { |
| 998 | assert_eq!((0 as $T).is_odd(), false); |
| 999 | assert_eq!((1 as $T).is_odd(), true); |
| 1000 | assert_eq!((2 as $T).is_odd(), false); |
| 1001 | assert_eq!((3 as $T).is_odd(), true); |
| 1002 | assert_eq!((4 as $T).is_odd(), false); |
| 1003 | } |
| 1004 | } |
| 1005 | }; |
| 1006 | } |
| 1007 | |
| 1008 | impl_integer_for_usize!(u8, test_integer_u8); |
| 1009 | impl_integer_for_usize!(u16, test_integer_u16); |
| 1010 | impl_integer_for_usize!(u32, test_integer_u32); |
| 1011 | impl_integer_for_usize!(u64, test_integer_u64); |
| 1012 | impl_integer_for_usize!(usize, test_integer_usize); |
| 1013 | #[cfg(has_i128)] |
| 1014 | impl_integer_for_usize!(u128, test_integer_u128); |
| 1015 | |
| 1016 | /// An iterator over binomial coefficients. |
| 1017 | pub struct IterBinomial<T> { |
| 1018 | a: T, |
| 1019 | n: T, |
| 1020 | k: T, |
| 1021 | } |
| 1022 | |
| 1023 | impl<T> IterBinomial<T> |
| 1024 | where |
| 1025 | T: Integer, |
| 1026 | { |
| 1027 | /// For a given n, iterate over all binomial coefficients binomial(n, k), for k=0...n. |
| 1028 | /// |
| 1029 | /// Note that this might overflow, depending on `T`. For the primitive |
| 1030 | /// integer types, the following n are the largest ones for which there will |
| 1031 | /// be no overflow: |
| 1032 | /// |
| 1033 | /// type | n |
| 1034 | /// -----|--- |
| 1035 | /// u8 | 10 |
| 1036 | /// i8 | 9 |
| 1037 | /// u16 | 18 |
| 1038 | /// i16 | 17 |
| 1039 | /// u32 | 34 |
| 1040 | /// i32 | 33 |
| 1041 | /// u64 | 67 |
| 1042 | /// i64 | 66 |
| 1043 | /// |
| 1044 | /// For larger n, `T` should be a bigint type. |
| 1045 | pub fn new(n: T) -> IterBinomial<T> { |
| 1046 | IterBinomial { |
| 1047 | k: T::zero(), |
| 1048 | a: T::one(), |
| 1049 | n: n, |
| 1050 | } |
| 1051 | } |
| 1052 | } |
| 1053 | |
| 1054 | impl<T> Iterator for IterBinomial<T> |
| 1055 | where |
| 1056 | T: Integer + Clone, |
| 1057 | { |
| 1058 | type Item = T; |
| 1059 | |
| 1060 | fn next(&mut self) -> Option<T> { |
| 1061 | if self.k > self.n { |
| 1062 | return None; |
| 1063 | } |
| 1064 | self.a = if !self.k.is_zero() { |
| 1065 | multiply_and_divide( |
| 1066 | self.a.clone(), |
| 1067 | self.n.clone() - self.k.clone() + T::one(), |
| 1068 | self.k.clone(), |
| 1069 | ) |
| 1070 | } else { |
| 1071 | T::one() |
| 1072 | }; |
| 1073 | self.k = self.k.clone() + T::one(); |
| 1074 | Some(self.a.clone()) |
| 1075 | } |
| 1076 | } |
| 1077 | |
| 1078 | /// Calculate r * a / b, avoiding overflows and fractions. |
| 1079 | /// |
| 1080 | /// Assumes that b divides r * a evenly. |
| 1081 | fn multiply_and_divide<T: Integer + Clone>(r: T, a: T, b: T) -> T { |
| 1082 | // See http://blog.plover.com/math/choose-2.html for the idea. |
| 1083 | let g = gcd(r.clone(), b.clone()); |
| 1084 | r / g.clone() * (a / (b / g)) |
| 1085 | } |
| 1086 | |
| 1087 | /// Calculate the binomial coefficient. |
| 1088 | /// |
| 1089 | /// Note that this might overflow, depending on `T`. For the primitive integer |
| 1090 | /// types, the following n are the largest ones possible such that there will |
| 1091 | /// be no overflow for any k: |
| 1092 | /// |
| 1093 | /// type | n |
| 1094 | /// -----|--- |
| 1095 | /// u8 | 10 |
| 1096 | /// i8 | 9 |
| 1097 | /// u16 | 18 |
| 1098 | /// i16 | 17 |
| 1099 | /// u32 | 34 |
| 1100 | /// i32 | 33 |
| 1101 | /// u64 | 67 |
| 1102 | /// i64 | 66 |
| 1103 | /// |
| 1104 | /// For larger n, consider using a bigint type for `T`. |
| 1105 | pub fn binomial<T: Integer + Clone>(mut n: T, k: T) -> T { |
| 1106 | // See http://blog.plover.com/math/choose.html for the idea. |
| 1107 | if k > n { |
| 1108 | return T::zero(); |
| 1109 | } |
| 1110 | if k > n.clone() - k.clone() { |
| 1111 | return binomial(n.clone(), n - k); |
| 1112 | } |
| 1113 | let mut r = T::one(); |
| 1114 | let mut d = T::one(); |
| 1115 | loop { |
| 1116 | if d > k { |
| 1117 | break; |
| 1118 | } |
| 1119 | r = multiply_and_divide(r, n.clone(), d.clone()); |
| 1120 | n = n - T::one(); |
| 1121 | d = d + T::one(); |
| 1122 | } |
| 1123 | r |
| 1124 | } |
| 1125 | |
| 1126 | /// Calculate the multinomial coefficient. |
| 1127 | pub fn multinomial<T: Integer + Clone>(k: &[T]) -> T |
| 1128 | where |
| 1129 | for<'a> T: Add<&'a T, Output = T>, |
| 1130 | { |
| 1131 | let mut r = T::one(); |
| 1132 | let mut p = T::zero(); |
| 1133 | for i in k { |
| 1134 | p = p + i; |
| 1135 | r = r * binomial(p.clone(), i.clone()); |
| 1136 | } |
| 1137 | r |
| 1138 | } |
| 1139 | |
| 1140 | #[test] |
| 1141 | fn test_lcm_overflow() { |
| 1142 | macro_rules! check { |
| 1143 | ($t:ty, $x:expr, $y:expr, $r:expr) => {{ |
| 1144 | let x: $t = $x; |
| 1145 | let y: $t = $y; |
| 1146 | let o = x.checked_mul(y); |
| 1147 | assert!( |
| 1148 | o.is_none(), |
| 1149 | "sanity checking that {} input {} * {} overflows", |
| 1150 | stringify!($t), |
| 1151 | x, |
| 1152 | y |
| 1153 | ); |
| 1154 | assert_eq!(x.lcm(&y), $r); |
| 1155 | assert_eq!(y.lcm(&x), $r); |
| 1156 | }}; |
| 1157 | } |
| 1158 | |
| 1159 | // Original bug (Issue #166) |
| 1160 | check!(i64, 46656000000000000, 600, 46656000000000000); |
| 1161 | |
| 1162 | check!(i8, 0x40, 0x04, 0x40); |
| 1163 | check!(u8, 0x80, 0x02, 0x80); |
| 1164 | check!(i16, 0x40_00, 0x04, 0x40_00); |
| 1165 | check!(u16, 0x80_00, 0x02, 0x80_00); |
| 1166 | check!(i32, 0x4000_0000, 0x04, 0x4000_0000); |
| 1167 | check!(u32, 0x8000_0000, 0x02, 0x8000_0000); |
| 1168 | check!(i64, 0x4000_0000_0000_0000, 0x04, 0x4000_0000_0000_0000); |
| 1169 | check!(u64, 0x8000_0000_0000_0000, 0x02, 0x8000_0000_0000_0000); |
| 1170 | } |
| 1171 | |
| 1172 | #[test] |
| 1173 | fn test_iter_binomial() { |
| 1174 | macro_rules! check_simple { |
| 1175 | ($t:ty) => {{ |
| 1176 | let n: $t = 3; |
| 1177 | let expected = [1, 3, 3, 1]; |
| 1178 | for (b, &e) in IterBinomial::new(n).zip(&expected) { |
| 1179 | assert_eq!(b, e); |
| 1180 | } |
| 1181 | }}; |
| 1182 | } |
| 1183 | |
| 1184 | check_simple!(u8); |
| 1185 | check_simple!(i8); |
| 1186 | check_simple!(u16); |
| 1187 | check_simple!(i16); |
| 1188 | check_simple!(u32); |
| 1189 | check_simple!(i32); |
| 1190 | check_simple!(u64); |
| 1191 | check_simple!(i64); |
| 1192 | |
| 1193 | macro_rules! check_binomial { |
| 1194 | ($t:ty, $n:expr) => {{ |
| 1195 | let n: $t = $n; |
| 1196 | let mut k: $t = 0; |
| 1197 | for b in IterBinomial::new(n) { |
| 1198 | assert_eq!(b, binomial(n, k)); |
| 1199 | k += 1; |
| 1200 | } |
| 1201 | }}; |
| 1202 | } |
| 1203 | |
| 1204 | // Check the largest n for which there is no overflow. |
| 1205 | check_binomial!(u8, 10); |
| 1206 | check_binomial!(i8, 9); |
| 1207 | check_binomial!(u16, 18); |
| 1208 | check_binomial!(i16, 17); |
| 1209 | check_binomial!(u32, 34); |
| 1210 | check_binomial!(i32, 33); |
| 1211 | check_binomial!(u64, 67); |
| 1212 | check_binomial!(i64, 66); |
| 1213 | } |
| 1214 | |
| 1215 | #[test] |
| 1216 | fn test_binomial() { |
| 1217 | macro_rules! check { |
| 1218 | ($t:ty, $x:expr, $y:expr, $r:expr) => {{ |
| 1219 | let x: $t = $x; |
| 1220 | let y: $t = $y; |
| 1221 | let expected: $t = $r; |
| 1222 | assert_eq!(binomial(x, y), expected); |
| 1223 | if y <= x { |
| 1224 | assert_eq!(binomial(x, x - y), expected); |
| 1225 | } |
| 1226 | }}; |
| 1227 | } |
| 1228 | check!(u8, 9, 4, 126); |
| 1229 | check!(u8, 0, 0, 1); |
| 1230 | check!(u8, 2, 3, 0); |
| 1231 | |
| 1232 | check!(i8, 9, 4, 126); |
| 1233 | check!(i8, 0, 0, 1); |
| 1234 | check!(i8, 2, 3, 0); |
| 1235 | |
| 1236 | check!(u16, 100, 2, 4950); |
| 1237 | check!(u16, 14, 4, 1001); |
| 1238 | check!(u16, 0, 0, 1); |
| 1239 | check!(u16, 2, 3, 0); |
| 1240 | |
| 1241 | check!(i16, 100, 2, 4950); |
| 1242 | check!(i16, 14, 4, 1001); |
| 1243 | check!(i16, 0, 0, 1); |
| 1244 | check!(i16, 2, 3, 0); |
| 1245 | |
| 1246 | check!(u32, 100, 2, 4950); |
| 1247 | check!(u32, 35, 11, 417225900); |
| 1248 | check!(u32, 14, 4, 1001); |
| 1249 | check!(u32, 0, 0, 1); |
| 1250 | check!(u32, 2, 3, 0); |
| 1251 | |
| 1252 | check!(i32, 100, 2, 4950); |
| 1253 | check!(i32, 35, 11, 417225900); |
| 1254 | check!(i32, 14, 4, 1001); |
| 1255 | check!(i32, 0, 0, 1); |
| 1256 | check!(i32, 2, 3, 0); |
| 1257 | |
| 1258 | check!(u64, 100, 2, 4950); |
| 1259 | check!(u64, 35, 11, 417225900); |
| 1260 | check!(u64, 14, 4, 1001); |
| 1261 | check!(u64, 0, 0, 1); |
| 1262 | check!(u64, 2, 3, 0); |
| 1263 | |
| 1264 | check!(i64, 100, 2, 4950); |
| 1265 | check!(i64, 35, 11, 417225900); |
| 1266 | check!(i64, 14, 4, 1001); |
| 1267 | check!(i64, 0, 0, 1); |
| 1268 | check!(i64, 2, 3, 0); |
| 1269 | } |
| 1270 | |
| 1271 | #[test] |
| 1272 | fn test_multinomial() { |
| 1273 | macro_rules! check_binomial { |
| 1274 | ($t:ty, $k:expr) => {{ |
| 1275 | let n: $t = $k.iter().fold(0, |acc, &x| acc + x); |
| 1276 | let k: &[$t] = $k; |
| 1277 | assert_eq!(k.len(), 2); |
| 1278 | assert_eq!(multinomial(k), binomial(n, k[0])); |
| 1279 | }}; |
| 1280 | } |
| 1281 | |
| 1282 | check_binomial!(u8, &[4, 5]); |
| 1283 | |
| 1284 | check_binomial!(i8, &[4, 5]); |
| 1285 | |
| 1286 | check_binomial!(u16, &[2, 98]); |
| 1287 | check_binomial!(u16, &[4, 10]); |
| 1288 | |
| 1289 | check_binomial!(i16, &[2, 98]); |
| 1290 | check_binomial!(i16, &[4, 10]); |
| 1291 | |
| 1292 | check_binomial!(u32, &[2, 98]); |
| 1293 | check_binomial!(u32, &[11, 24]); |
| 1294 | check_binomial!(u32, &[4, 10]); |
| 1295 | |
| 1296 | check_binomial!(i32, &[2, 98]); |
| 1297 | check_binomial!(i32, &[11, 24]); |
| 1298 | check_binomial!(i32, &[4, 10]); |
| 1299 | |
| 1300 | check_binomial!(u64, &[2, 98]); |
| 1301 | check_binomial!(u64, &[11, 24]); |
| 1302 | check_binomial!(u64, &[4, 10]); |
| 1303 | |
| 1304 | check_binomial!(i64, &[2, 98]); |
| 1305 | check_binomial!(i64, &[11, 24]); |
| 1306 | check_binomial!(i64, &[4, 10]); |
| 1307 | |
| 1308 | macro_rules! check_multinomial { |
| 1309 | ($t:ty, $k:expr, $r:expr) => {{ |
| 1310 | let k: &[$t] = $k; |
| 1311 | let expected: $t = $r; |
| 1312 | assert_eq!(multinomial(k), expected); |
| 1313 | }}; |
| 1314 | } |
| 1315 | |
| 1316 | check_multinomial!(u8, &[2, 1, 2], 30); |
| 1317 | check_multinomial!(u8, &[2, 3, 0], 10); |
| 1318 | |
| 1319 | check_multinomial!(i8, &[2, 1, 2], 30); |
| 1320 | check_multinomial!(i8, &[2, 3, 0], 10); |
| 1321 | |
| 1322 | check_multinomial!(u16, &[2, 1, 2], 30); |
| 1323 | check_multinomial!(u16, &[2, 3, 0], 10); |
| 1324 | |
| 1325 | check_multinomial!(i16, &[2, 1, 2], 30); |
| 1326 | check_multinomial!(i16, &[2, 3, 0], 10); |
| 1327 | |
| 1328 | check_multinomial!(u32, &[2, 1, 2], 30); |
| 1329 | check_multinomial!(u32, &[2, 3, 0], 10); |
| 1330 | |
| 1331 | check_multinomial!(i32, &[2, 1, 2], 30); |
| 1332 | check_multinomial!(i32, &[2, 3, 0], 10); |
| 1333 | |
| 1334 | check_multinomial!(u64, &[2, 1, 2], 30); |
| 1335 | check_multinomial!(u64, &[2, 3, 0], 10); |
| 1336 | |
| 1337 | check_multinomial!(i64, &[2, 1, 2], 30); |
| 1338 | check_multinomial!(i64, &[2, 3, 0], 10); |
| 1339 | |
| 1340 | check_multinomial!(u64, &[], 1); |
| 1341 | check_multinomial!(u64, &[0], 1); |
| 1342 | check_multinomial!(u64, &[12345], 1); |
| 1343 | } |