| 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.31 and greater. |
| 16 | |
| 17 | #![doc (html_root_url = "https://docs.rs/num-integer/0.1" )] |
| 18 | #![no_std ] |
| 19 | |
| 20 | use core::mem; |
| 21 | use core::ops::Add; |
| 22 | |
| 23 | use num_traits::{Num, Signed, Zero}; |
| 24 | |
| 25 | mod roots; |
| 26 | pub use crate::roots::Roots; |
| 27 | pub use crate::roots::{cbrt, nth_root, sqrt}; |
| 28 | |
| 29 | mod average; |
| 30 | pub use crate::average::Average; |
| 31 | pub use crate::average::{average_ceil, average_floor}; |
| 32 | |
| 33 | pub trait Integer: Sized + Num + PartialOrd + Ord + Eq { |
| 34 | /// Floored integer division. |
| 35 | /// |
| 36 | /// # Examples |
| 37 | /// |
| 38 | /// ~~~ |
| 39 | /// # use num_integer::Integer; |
| 40 | /// assert!(( 8).div_floor(& 3) == 2); |
| 41 | /// assert!(( 8).div_floor(&-3) == -3); |
| 42 | /// assert!((-8).div_floor(& 3) == -3); |
| 43 | /// assert!((-8).div_floor(&-3) == 2); |
| 44 | /// |
| 45 | /// assert!(( 1).div_floor(& 2) == 0); |
| 46 | /// assert!(( 1).div_floor(&-2) == -1); |
| 47 | /// assert!((-1).div_floor(& 2) == -1); |
| 48 | /// assert!((-1).div_floor(&-2) == 0); |
| 49 | /// ~~~ |
| 50 | fn div_floor(&self, other: &Self) -> Self; |
| 51 | |
| 52 | /// Floored integer modulo, satisfying: |
| 53 | /// |
| 54 | /// ~~~ |
| 55 | /// # use num_integer::Integer; |
| 56 | /// # let n = 1; let d = 1; |
| 57 | /// assert!(n.div_floor(&d) * d + n.mod_floor(&d) == n) |
| 58 | /// ~~~ |
| 59 | /// |
| 60 | /// # Examples |
| 61 | /// |
| 62 | /// ~~~ |
| 63 | /// # use num_integer::Integer; |
| 64 | /// assert!(( 8).mod_floor(& 3) == 2); |
| 65 | /// assert!(( 8).mod_floor(&-3) == -1); |
| 66 | /// assert!((-8).mod_floor(& 3) == 1); |
| 67 | /// assert!((-8).mod_floor(&-3) == -2); |
| 68 | /// |
| 69 | /// assert!(( 1).mod_floor(& 2) == 1); |
| 70 | /// assert!(( 1).mod_floor(&-2) == -1); |
| 71 | /// assert!((-1).mod_floor(& 2) == 1); |
| 72 | /// assert!((-1).mod_floor(&-2) == -1); |
| 73 | /// ~~~ |
| 74 | fn mod_floor(&self, other: &Self) -> Self; |
| 75 | |
| 76 | /// Ceiled integer division. |
| 77 | /// |
| 78 | /// # Examples |
| 79 | /// |
| 80 | /// ~~~ |
| 81 | /// # use num_integer::Integer; |
| 82 | /// assert_eq!(( 8).div_ceil( &3), 3); |
| 83 | /// assert_eq!(( 8).div_ceil(&-3), -2); |
| 84 | /// assert_eq!((-8).div_ceil( &3), -2); |
| 85 | /// assert_eq!((-8).div_ceil(&-3), 3); |
| 86 | /// |
| 87 | /// assert_eq!(( 1).div_ceil( &2), 1); |
| 88 | /// assert_eq!(( 1).div_ceil(&-2), 0); |
| 89 | /// assert_eq!((-1).div_ceil( &2), 0); |
| 90 | /// assert_eq!((-1).div_ceil(&-2), 1); |
| 91 | /// ~~~ |
| 92 | 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 | /// Greatest Common Divisor (GCD). |
| 102 | /// |
| 103 | /// # Examples |
| 104 | /// |
| 105 | /// ~~~ |
| 106 | /// # use num_integer::Integer; |
| 107 | /// assert_eq!(6.gcd(&8), 2); |
| 108 | /// assert_eq!(7.gcd(&3), 1); |
| 109 | /// ~~~ |
| 110 | fn gcd(&self, other: &Self) -> Self; |
| 111 | |
| 112 | /// Lowest Common Multiple (LCM). |
| 113 | /// |
| 114 | /// # Examples |
| 115 | /// |
| 116 | /// ~~~ |
| 117 | /// # use num_integer::Integer; |
| 118 | /// assert_eq!(7.lcm(&3), 21); |
| 119 | /// assert_eq!(2.lcm(&4), 4); |
| 120 | /// assert_eq!(0.lcm(&0), 0); |
| 121 | /// ~~~ |
| 122 | fn lcm(&self, other: &Self) -> Self; |
| 123 | |
| 124 | /// Greatest Common Divisor (GCD) and |
| 125 | /// Lowest Common Multiple (LCM) together. |
| 126 | /// |
| 127 | /// Potentially more efficient than calling `gcd` and `lcm` |
| 128 | /// individually for identical inputs. |
| 129 | /// |
| 130 | /// # Examples |
| 131 | /// |
| 132 | /// ~~~ |
| 133 | /// # use num_integer::Integer; |
| 134 | /// assert_eq!(10.gcd_lcm(&4), (2, 20)); |
| 135 | /// assert_eq!(8.gcd_lcm(&9), (1, 72)); |
| 136 | /// ~~~ |
| 137 | #[inline ] |
| 138 | fn gcd_lcm(&self, other: &Self) -> (Self, Self) { |
| 139 | (self.gcd(other), self.lcm(other)) |
| 140 | } |
| 141 | |
| 142 | /// Greatest common divisor and Bézout coefficients. |
| 143 | /// |
| 144 | /// # Examples |
| 145 | /// |
| 146 | /// ~~~ |
| 147 | /// # fn main() { |
| 148 | /// # use num_integer::{ExtendedGcd, Integer}; |
| 149 | /// # use num_traits::NumAssign; |
| 150 | /// fn check<A: Copy + Integer + NumAssign>(a: A, b: A) -> bool { |
| 151 | /// let ExtendedGcd { gcd, x, y, .. } = a.extended_gcd(&b); |
| 152 | /// gcd == x * a + y * b |
| 153 | /// } |
| 154 | /// assert!(check(10isize, 4isize)); |
| 155 | /// assert!(check(8isize, 9isize)); |
| 156 | /// # } |
| 157 | /// ~~~ |
| 158 | #[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 | /// Greatest common divisor, least common multiple, and Bézout coefficients. |
| 195 | #[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, use `is_multiple_of` instead. |
| 204 | #[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 | /// Returns `true` if `self` is a multiple of `other`. |
| 211 | /// |
| 212 | /// # Examples |
| 213 | /// |
| 214 | /// ~~~ |
| 215 | /// # use num_integer::Integer; |
| 216 | /// assert_eq!(9.is_multiple_of(&3), true); |
| 217 | /// assert_eq!(3.is_multiple_of(&9), false); |
| 218 | /// ~~~ |
| 219 | fn is_multiple_of(&self, other: &Self) -> bool; |
| 220 | |
| 221 | /// Returns `true` if the number is even. |
| 222 | /// |
| 223 | /// # Examples |
| 224 | /// |
| 225 | /// ~~~ |
| 226 | /// # use num_integer::Integer; |
| 227 | /// assert_eq!(3.is_even(), false); |
| 228 | /// assert_eq!(4.is_even(), true); |
| 229 | /// ~~~ |
| 230 | fn is_even(&self) -> bool; |
| 231 | |
| 232 | /// Returns `true` if the number is odd. |
| 233 | /// |
| 234 | /// # Examples |
| 235 | /// |
| 236 | /// ~~~ |
| 237 | /// # use num_integer::Integer; |
| 238 | /// assert_eq!(3.is_odd(), true); |
| 239 | /// assert_eq!(4.is_odd(), false); |
| 240 | /// ~~~ |
| 241 | fn is_odd(&self) -> bool; |
| 242 | |
| 243 | /// Simultaneous truncated integer division and modulus. |
| 244 | /// Returns `(quotient, remainder)`. |
| 245 | /// |
| 246 | /// # Examples |
| 247 | /// |
| 248 | /// ~~~ |
| 249 | /// # use num_integer::Integer; |
| 250 | /// assert_eq!(( 8).div_rem( &3), ( 2, 2)); |
| 251 | /// assert_eq!(( 8).div_rem(&-3), (-2, 2)); |
| 252 | /// assert_eq!((-8).div_rem( &3), (-2, -2)); |
| 253 | /// assert_eq!((-8).div_rem(&-3), ( 2, -2)); |
| 254 | /// |
| 255 | /// assert_eq!(( 1).div_rem( &2), ( 0, 1)); |
| 256 | /// assert_eq!(( 1).div_rem(&-2), ( 0, 1)); |
| 257 | /// assert_eq!((-1).div_rem( &2), ( 0, -1)); |
| 258 | /// assert_eq!((-1).div_rem(&-2), ( 0, -1)); |
| 259 | /// ~~~ |
| 260 | fn div_rem(&self, other: &Self) -> (Self, Self); |
| 261 | |
| 262 | /// Simultaneous floored integer division and modulus. |
| 263 | /// Returns `(quotient, remainder)`. |
| 264 | /// |
| 265 | /// # Examples |
| 266 | /// |
| 267 | /// ~~~ |
| 268 | /// # use num_integer::Integer; |
| 269 | /// assert_eq!(( 8).div_mod_floor( &3), ( 2, 2)); |
| 270 | /// assert_eq!(( 8).div_mod_floor(&-3), (-3, -1)); |
| 271 | /// assert_eq!((-8).div_mod_floor( &3), (-3, 1)); |
| 272 | /// assert_eq!((-8).div_mod_floor(&-3), ( 2, -2)); |
| 273 | /// |
| 274 | /// assert_eq!(( 1).div_mod_floor( &2), ( 0, 1)); |
| 275 | /// assert_eq!(( 1).div_mod_floor(&-2), (-1, -1)); |
| 276 | /// assert_eq!((-1).div_mod_floor( &2), (-1, 1)); |
| 277 | /// assert_eq!((-1).div_mod_floor(&-2), ( 0, -1)); |
| 278 | /// ~~~ |
| 279 | fn div_mod_floor(&self, other: &Self) -> (Self, Self) { |
| 280 | (self.div_floor(other), self.mod_floor(other)) |
| 281 | } |
| 282 | |
| 283 | /// Rounds up to nearest multiple of argument. |
| 284 | /// |
| 285 | /// # Notes |
| 286 | /// |
| 287 | /// For signed types, `a.next_multiple_of(b) = a.prev_multiple_of(b.neg())`. |
| 288 | /// |
| 289 | /// # Examples |
| 290 | /// |
| 291 | /// ~~~ |
| 292 | /// # use num_integer::Integer; |
| 293 | /// assert_eq!(( 16).next_multiple_of(& 8), 16); |
| 294 | /// assert_eq!(( 23).next_multiple_of(& 8), 24); |
| 295 | /// assert_eq!(( 16).next_multiple_of(&-8), 16); |
| 296 | /// assert_eq!(( 23).next_multiple_of(&-8), 16); |
| 297 | /// assert_eq!((-16).next_multiple_of(& 8), -16); |
| 298 | /// assert_eq!((-23).next_multiple_of(& 8), -16); |
| 299 | /// assert_eq!((-16).next_multiple_of(&-8), -16); |
| 300 | /// assert_eq!((-23).next_multiple_of(&-8), -24); |
| 301 | /// ~~~ |
| 302 | #[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 | /// Rounds down to nearest multiple of argument. |
| 317 | /// |
| 318 | /// # Notes |
| 319 | /// |
| 320 | /// For signed types, `a.prev_multiple_of(b) = a.next_multiple_of(b.neg())`. |
| 321 | /// |
| 322 | /// # Examples |
| 323 | /// |
| 324 | /// ~~~ |
| 325 | /// # use num_integer::Integer; |
| 326 | /// assert_eq!(( 16).prev_multiple_of(& 8), 16); |
| 327 | /// assert_eq!(( 23).prev_multiple_of(& 8), 16); |
| 328 | /// assert_eq!(( 16).prev_multiple_of(&-8), 16); |
| 329 | /// assert_eq!(( 23).prev_multiple_of(&-8), 24); |
| 330 | /// assert_eq!((-16).prev_multiple_of(& 8), -16); |
| 331 | /// assert_eq!((-23).prev_multiple_of(& 8), -24); |
| 332 | /// assert_eq!((-16).prev_multiple_of(&-8), -16); |
| 333 | /// assert_eq!((-23).prev_multiple_of(&-8), -16); |
| 334 | /// ~~~ |
| 335 | #[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 | /// Decrements self by one. |
| 344 | /// |
| 345 | /// # Examples |
| 346 | /// |
| 347 | /// ~~~ |
| 348 | /// # use num_integer::Integer; |
| 349 | /// let mut x: i32 = 43; |
| 350 | /// x.dec(); |
| 351 | /// assert_eq!(x, 42); |
| 352 | /// ~~~ |
| 353 | fn dec(&mut self) |
| 354 | where |
| 355 | Self: Clone, |
| 356 | { |
| 357 | *self = self.clone() - Self::one() |
| 358 | } |
| 359 | |
| 360 | /// Increments self by one. |
| 361 | /// |
| 362 | /// # Examples |
| 363 | /// |
| 364 | /// ~~~ |
| 365 | /// # use num_integer::Integer; |
| 366 | /// let mut x: i32 = 41; |
| 367 | /// x.inc(); |
| 368 | /// assert_eq!(x, 42); |
| 369 | /// ~~~ |
| 370 | fn inc(&mut self) |
| 371 | where |
| 372 | Self: Clone, |
| 373 | { |
| 374 | *self = self.clone() + Self::one() |
| 375 | } |
| 376 | } |
| 377 | |
| 378 | /// Greatest common divisor and Bézout coefficients |
| 379 | /// |
| 380 | /// ```no_build |
| 381 | /// let e = isize::extended_gcd(a, b); |
| 382 | /// assert_eq!(e.gcd, e.x*a + e.y*b); |
| 383 | /// ``` |
| 384 | #[derive (Debug, Clone, Copy, PartialEq, Eq)] |
| 385 | pub struct ExtendedGcd<A> { |
| 386 | pub gcd: A, |
| 387 | pub x: A, |
| 388 | pub y: A, |
| 389 | } |
| 390 | |
| 391 | /// Simultaneous integer division and modulus |
| 392 | #[inline ] |
| 393 | pub fn div_rem<T: Integer>(x: T, y: T) -> (T, T) { |
| 394 | x.div_rem(&y) |
| 395 | } |
| 396 | /// Floored integer division |
| 397 | #[inline ] |
| 398 | pub fn div_floor<T: Integer>(x: T, y: T) -> T { |
| 399 | x.div_floor(&y) |
| 400 | } |
| 401 | /// Floored integer modulus |
| 402 | #[inline ] |
| 403 | pub fn mod_floor<T: Integer>(x: T, y: T) -> T { |
| 404 | x.mod_floor(&y) |
| 405 | } |
| 406 | /// Simultaneous floored integer division and modulus |
| 407 | #[inline ] |
| 408 | pub fn div_mod_floor<T: Integer>(x: T, y: T) -> (T, T) { |
| 409 | x.div_mod_floor(&y) |
| 410 | } |
| 411 | /// Ceiled integer division |
| 412 | #[inline ] |
| 413 | pub fn div_ceil<T: Integer>(x: T, y: T) -> T { |
| 414 | x.div_ceil(&y) |
| 415 | } |
| 416 | |
| 417 | /// Calculates the Greatest Common Divisor (GCD) of the number and `other`. The |
| 418 | /// result is always non-negative. |
| 419 | #[inline (always)] |
| 420 | pub fn gcd<T: Integer>(x: T, y: T) -> T { |
| 421 | x.gcd(&y) |
| 422 | } |
| 423 | /// Calculates the Lowest Common Multiple (LCM) of the number and `other`. |
| 424 | #[inline (always)] |
| 425 | pub fn lcm<T: Integer>(x: T, y: T) -> T { |
| 426 | x.lcm(&y) |
| 427 | } |
| 428 | |
| 429 | /// Calculates the Greatest Common Divisor (GCD) and |
| 430 | /// Lowest Common Multiple (LCM) of the number and `other`. |
| 431 | #[inline (always)] |
| 432 | pub fn gcd_lcm<T: Integer>(x: T, y: T) -> (T, T) { |
| 433 | x.gcd_lcm(&y) |
| 434 | } |
| 435 | |
| 436 | macro_rules! impl_integer_for_isize { |
| 437 | ($T:ty, $test_mod:ident) => { |
| 438 | impl Integer for $T { |
| 439 | /// Floored integer division |
| 440 | #[inline] |
| 441 | fn div_floor(&self, other: &Self) -> Self { |
| 442 | // Algorithm from [Daan Leijen. _Division and Modulus for Computer Scientists_, |
| 443 | // December 2001](http://research.microsoft.com/pubs/151917/divmodnote-letter.pdf) |
| 444 | 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 | /// Floored integer modulo |
| 453 | #[inline] |
| 454 | fn mod_floor(&self, other: &Self) -> Self { |
| 455 | // Algorithm from [Daan Leijen. _Division and Modulus for Computer Scientists_, |
| 456 | // December 2001](http://research.microsoft.com/pubs/151917/divmodnote-letter.pdf) |
| 457 | 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 | /// Calculates `div_floor` and `mod_floor` simultaneously |
| 466 | #[inline] |
| 467 | fn div_mod_floor(&self, other: &Self) -> (Self, Self) { |
| 468 | // Algorithm from [Daan Leijen. _Division and Modulus for Computer Scientists_, |
| 469 | // December 2001](http://research.microsoft.com/pubs/151917/divmodnote-letter.pdf) |
| 470 | 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 | /// Calculates the Greatest Common Divisor (GCD) of the number and |
| 489 | /// `other`. The result is always non-negative. |
| 490 | #[inline] |
| 491 | fn gcd(&self, other: &Self) -> Self { |
| 492 | // Use Stein's algorithm |
| 493 | let mut m = *self; |
| 494 | let mut n = *other; |
| 495 | if m == 0 || n == 0 { |
| 496 | return (m | n).abs(); |
| 497 | } |
| 498 | |
| 499 | // find common factors of 2 |
| 500 | let shift = (m | n).trailing_zeros(); |
| 501 | |
| 502 | // The algorithm needs positive numbers, but the minimum value |
| 503 | // can't be represented as a positive one. |
| 504 | // It's also a power of two, so the gcd can be |
| 505 | // calculated by bitshifting in that case |
| 506 | |
| 507 | // Assuming two's complement, the number created by the shift |
| 508 | // is positive for all numbers except gcd = abs(min value) |
| 509 | // The call to .abs() causes a panic in debug mode |
| 510 | if m == Self::min_value() || n == Self::min_value() { |
| 511 | return (1 << shift).abs(); |
| 512 | } |
| 513 | |
| 514 | // guaranteed to be positive now, rest like unsigned algorithm |
| 515 | m = m.abs(); |
| 516 | n = n.abs(); |
| 517 | |
| 518 | // divide n and m by 2 until odd |
| 519 | 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 | // should not have to recalculate abs |
| 538 | 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 | /// Calculates the Lowest Common Multiple (LCM) of the number and |
| 547 | /// `other`. |
| 548 | #[inline] |
| 549 | fn lcm(&self, other: &Self) -> Self { |
| 550 | self.gcd_lcm(other).1 |
| 551 | } |
| 552 | |
| 553 | /// Calculates the Greatest Common Divisor (GCD) and |
| 554 | /// Lowest Common Multiple (LCM) of the number and `other`. |
| 555 | #[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 | // should not have to recalculate abs |
| 562 | let lcm = (*self * (*other / gcd)).abs(); |
| 563 | (gcd, lcm) |
| 564 | } |
| 565 | |
| 566 | /// Returns `true` if the number is a multiple of `other`. |
| 567 | #[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 | /// Returns `true` if the number is divisible by `2` |
| 576 | #[inline] |
| 577 | fn is_even(&self) -> bool { |
| 578 | (*self) & 1 == 0 |
| 579 | } |
| 580 | |
| 581 | /// Returns `true` if the number is not divisible by `2` |
| 582 | #[inline] |
| 583 | fn is_odd(&self) -> bool { |
| 584 | !self.is_even() |
| 585 | } |
| 586 | |
| 587 | /// Simultaneous truncated integer division and modulus. |
| 588 | #[inline] |
| 589 | fn div_rem(&self, other: &Self) -> (Self, Self) { |
| 590 | (*self / *other, *self % *other) |
| 591 | } |
| 592 | |
| 593 | /// Rounds up to nearest multiple of argument. |
| 594 | #[inline] |
| 595 | fn next_multiple_of(&self, other: &Self) -> Self { |
| 596 | // Avoid the overflow of `MIN % -1` |
| 597 | 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 | /// Rounds down to nearest multiple of argument. |
| 606 | #[inline] |
| 607 | fn prev_multiple_of(&self, other: &Self) -> Self { |
| 608 | // Avoid the overflow of `MIN % -1` |
| 609 | 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 | /// Checks that the division rule holds for: |
| 623 | /// |
| 624 | /// - `n`: numerator (dividend) |
| 625 | /// - `d`: denominator (divisor) |
| 626 | /// - `qr`: quotient and remainder |
| 627 | #[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 | // gcd(-128, b) = 128 is not representable as positive value |
| 705 | // for i8 |
| 706 | 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 | |
| 841 | impl_integer_for_isize!(i8, test_integer_i8); |
| 842 | impl_integer_for_isize!(i16, test_integer_i16); |
| 843 | impl_integer_for_isize!(i32, test_integer_i32); |
| 844 | impl_integer_for_isize!(i64, test_integer_i64); |
| 845 | impl_integer_for_isize!(i128, test_integer_i128); |
| 846 | impl_integer_for_isize!(isize, test_integer_isize); |
| 847 | |
| 848 | macro_rules! impl_integer_for_usize { |
| 849 | ($T:ty, $test_mod:ident) => { |
| 850 | impl Integer for $T { |
| 851 | /// Unsigned integer division. Returns the same result as `div` (`/`). |
| 852 | #[inline] |
| 853 | fn div_floor(&self, other: &Self) -> Self { |
| 854 | *self / *other |
| 855 | } |
| 856 | |
| 857 | /// Unsigned integer modulo operation. Returns the same result as `rem` (`%`). |
| 858 | #[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 | /// Calculates the Greatest Common Divisor (GCD) of the number and `other` |
| 869 | #[inline] |
| 870 | fn gcd(&self, other: &Self) -> Self { |
| 871 | // Use Stein's algorithm |
| 872 | let mut m = *self; |
| 873 | let mut n = *other; |
| 874 | if m == 0 || n == 0 { |
| 875 | return m | n; |
| 876 | } |
| 877 | |
| 878 | // find common factors of 2 |
| 879 | let shift = (m | n).trailing_zeros(); |
| 880 | |
| 881 | // divide n and m by 2 until odd |
| 882 | 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 | // should not have to recalculate abs |
| 901 | let lcm = if egcd.gcd.is_zero() { |
| 902 | Self::zero() |
| 903 | } else { |
| 904 | *self * (*other / egcd.gcd) |
| 905 | }; |
| 906 | (egcd, lcm) |
| 907 | } |
| 908 | |
| 909 | /// Calculates the Lowest Common Multiple (LCM) of the number and `other`. |
| 910 | #[inline] |
| 911 | fn lcm(&self, other: &Self) -> Self { |
| 912 | self.gcd_lcm(other).1 |
| 913 | } |
| 914 | |
| 915 | /// Calculates the Greatest Common Divisor (GCD) and |
| 916 | /// Lowest Common Multiple (LCM) of the number and `other`. |
| 917 | #[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 | /// Returns `true` if the number is a multiple of `other`. |
| 928 | #[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 | /// Returns `true` if the number is divisible by `2`. |
| 937 | #[inline] |
| 938 | fn is_even(&self) -> bool { |
| 939 | *self % 2 == 0 |
| 940 | } |
| 941 | |
| 942 | /// Returns `true` if the number is not divisible by `2`. |
| 943 | #[inline] |
| 944 | fn is_odd(&self) -> bool { |
| 945 | !self.is_even() |
| 946 | } |
| 947 | |
| 948 | /// Simultaneous truncated integer division and modulus. |
| 949 | #[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 | |
| 1052 | impl_integer_for_usize!(u8, test_integer_u8); |
| 1053 | impl_integer_for_usize!(u16, test_integer_u16); |
| 1054 | impl_integer_for_usize!(u32, test_integer_u32); |
| 1055 | impl_integer_for_usize!(u64, test_integer_u64); |
| 1056 | impl_integer_for_usize!(u128, test_integer_u128); |
| 1057 | impl_integer_for_usize!(usize, test_integer_usize); |
| 1058 | |
| 1059 | /// An iterator over binomial coefficients. |
| 1060 | pub struct IterBinomial<T> { |
| 1061 | a: T, |
| 1062 | n: T, |
| 1063 | k: T, |
| 1064 | } |
| 1065 | |
| 1066 | impl<T> IterBinomial<T> |
| 1067 | where |
| 1068 | T: Integer, |
| 1069 | { |
| 1070 | /// For a given n, iterate over all binomial coefficients binomial(n, k), for k=0...n. |
| 1071 | /// |
| 1072 | /// Note that this might overflow, depending on `T`. For the primitive |
| 1073 | /// integer types, the following n are the largest ones for which there will |
| 1074 | /// be no overflow: |
| 1075 | /// |
| 1076 | /// type | n |
| 1077 | /// -----|--- |
| 1078 | /// u8 | 10 |
| 1079 | /// i8 | 9 |
| 1080 | /// u16 | 18 |
| 1081 | /// i16 | 17 |
| 1082 | /// u32 | 34 |
| 1083 | /// i32 | 33 |
| 1084 | /// u64 | 67 |
| 1085 | /// i64 | 66 |
| 1086 | /// |
| 1087 | /// For larger n, `T` should be a bigint type. |
| 1088 | pub fn new(n: T) -> IterBinomial<T> { |
| 1089 | IterBinomial { |
| 1090 | k: T::zero(), |
| 1091 | a: T::one(), |
| 1092 | n, |
| 1093 | } |
| 1094 | } |
| 1095 | } |
| 1096 | |
| 1097 | impl<T> Iterator for IterBinomial<T> |
| 1098 | where |
| 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 | |
| 1121 | /// Calculate r * a / b, avoiding overflows and fractions. |
| 1122 | /// |
| 1123 | /// Assumes that b divides r * a evenly. |
| 1124 | fn multiply_and_divide<T: Integer + Clone>(r: T, a: T, b: T) -> T { |
| 1125 | // See http://blog.plover.com/math/choose-2.html for the idea. |
| 1126 | let g: T = gcd(x:r.clone(), y:b.clone()); |
| 1127 | r / g.clone() * (a / (b / g)) |
| 1128 | } |
| 1129 | |
| 1130 | /// Calculate the binomial coefficient. |
| 1131 | /// |
| 1132 | /// Note that this might overflow, depending on `T`. For the primitive integer |
| 1133 | /// types, the following n are the largest ones possible such that there will |
| 1134 | /// be no overflow for any k: |
| 1135 | /// |
| 1136 | /// type | n |
| 1137 | /// -----|--- |
| 1138 | /// u8 | 10 |
| 1139 | /// i8 | 9 |
| 1140 | /// u16 | 18 |
| 1141 | /// i16 | 17 |
| 1142 | /// u32 | 34 |
| 1143 | /// i32 | 33 |
| 1144 | /// u64 | 67 |
| 1145 | /// i64 | 66 |
| 1146 | /// |
| 1147 | /// For larger n, consider using a bigint type for `T`. |
| 1148 | pub fn binomial<T: Integer + Clone>(mut n: T, k: T) -> T { |
| 1149 | // See http://blog.plover.com/math/choose.html for the idea. |
| 1150 | if k > n { |
| 1151 | return T::zero(); |
| 1152 | } |
| 1153 | if k > n.clone() - k.clone() { |
| 1154 | return binomial(n.clone(), k:n - k); |
| 1155 | } |
| 1156 | let mut r: T = T::one(); |
| 1157 | let mut d: T = T::one(); |
| 1158 | loop { |
| 1159 | if d > k { |
| 1160 | break; |
| 1161 | } |
| 1162 | r = multiply_and_divide(r, a:n.clone(), b:d.clone()); |
| 1163 | n = n - T::one(); |
| 1164 | d = d + T::one(); |
| 1165 | } |
| 1166 | r |
| 1167 | } |
| 1168 | |
| 1169 | /// Calculate the multinomial coefficient. |
| 1170 | pub fn multinomial<T: Integer + Clone>(k: &[T]) -> T |
| 1171 | where |
| 1172 | for<'a> T: Add<&'a T, Output = T>, |
| 1173 | { |
| 1174 | let mut r: T = T::one(); |
| 1175 | let mut p: T = T::zero(); |
| 1176 | for i: &T in k { |
| 1177 | p = p + i; |
| 1178 | r = r * binomial(n:p.clone(), k:i.clone()); |
| 1179 | } |
| 1180 | r |
| 1181 | } |
| 1182 | |
| 1183 | #[test ] |
| 1184 | fn 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 | // Original bug (Issue #166) |
| 1203 | 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 ] |
| 1216 | fn 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 the largest n for which there is no overflow. |
| 1248 | 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 ] |
| 1259 | fn 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 ] |
| 1315 | fn 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 | } |
| 1387 | |