1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760
//! //!Traits for structures with both an addition and multiplication operation //! //!While ring-like structures can technically apply to any pair of operations, for ease of use //!and to stick with mathematical convention, these operations are taken to be specifically addition //!and multiplication where multiplication distributes over addition. //! //!# Implementation //! //!This module build upon and behaves in a way similar to the [group-like](crate::algebra::group_like) //!structures. However, in addition to the basic group traits, the following properties are also considered: //! * Distributivity: //! * Multiplication can act on summed elements independently, //! ie`z*(x+y)=z*x+z*y` and `(x+y)*z=x*z+y*z` for all `x`, `y`, and `z`. //! * Represented by the [`Distributive`] trait //! * Note that unlike other marker traits, this one takes a type parameter so that //! multiplication between types can also be marked //! (such as scalar-vector or matrix-vector multiplication) //! * The Exponentiation and Logarthm operations: //! * Unary operations that map between addition and multiplication //! * Represented with the [`Exponential`] trait //! * Properties regarding partial divisibility and factoring //! * These traits _essentially_ measure how "integer-like" a particular ring is by abstracting //! many properties and operations commonly associated with integers and natural numbers //! * These include, but are not limited to: //! * [Divisibility] //! * [Primality] //! * [GCD] //! * [Euclidean Division](EuclideanDiv) //! * [Factorizability](Factorizable) //! * For more information, see each individual trait's documentation //! //!# Usage //! //!Similarly to the [group-like](crate::algebra::group_like) module, there are a number of //!trait aliases corresponding to a system of ring-like algebraic structures that form a heirarchy //!as represented in the following diagram: //! //!```text //! //! Add Abelian Group Mul Semigroup //! | | //! -------Ring-------- //! | //! ----------Unital Ring------------------------- //! | | | | //! | Commutative Ring Domain | //! | | | | //! | --------------- | //! | | | //! | Integral Domain | //! | | | //!Division Ring ---GCD Domain--- | //! | | | | //! | UFD Bezout Domain | //! | | | | //! | ------PID------- | //! | | | //! | Euclidean Domain | //! | | | //! -------------Field----------- Exponential Ring //! | | //! ------Exponential Field-------- //! //!``` //! //!Note also that in addition to the above system, there is a complementary system of "Semirings" //!that follows almost the same exact structure but doesn't require a subtraction operation. //! //!For more information, see each structure's documentation //! //!# Naming //! //!It is worth noting that in mathematical literature, there is some leniency in what certain //!structures are called, so some trait aliases may be named differently than what some are used to. //! //!In particular: //! * Some mathematicians and literature use the term "Ring" to refer to what is called a "Unital Ring" //! here, instead preferring the term "Rng" for Rings without a multiplicative identity. In fact, //! this convention is _probably_ the majority considering how little use there is for non-unital rings. //! However, the choice was made against using the "Rng" system for the simple reason that it is //! obviously highly prone to typesetting errors and is arguably pretty unintuitive for those learning //! the terminology. //! * The term "Semiring" doesn't really have an established meaning in mathematical convention, so //! the use here reflects a particular design choice more than a particular standard. //! use crate::algebra::*; use core::iter::Iterator; ///A marker trait for stucts whose multiplication operation preserves addition, ///ie `z*(x+y)=z*x+z*y` and `(x+y)*z=x*z+y*z` for all `x`, `y`, and `z`. pub trait Distributive<T=Self> {} /// ///Common methods regarding multiplicative inverses in a ring or semiring /// ///Do note that while for some rings these methods are relatively quick (like the integers), ///for others, such as polynomials, this test might actually be relatively expensive pub trait Divisibility: Sized { ///Determines if there exists an element `x` such that `self*x` = `rhs` fn divides(self, rhs: Self) -> bool; ///Finds an element `x` such that `self*x` = `rhs` if such an element exists fn divide(self, rhs: Self) -> Option<Self>; ///Determines if this element has a multiplicative inverse fn unit(&self) -> bool; ///Finds this element's multiplicative inverse if it exists fn inverse(self) -> Option<Self>; } /// ///Methods for testing irreduciblity and primality /// ///Do note that for obvious reasons these methods should be assumed to be expensive by default. ///Furthermore, given the difficulty of checking irreduciblity or primality for general rings, ///it is entirely possible that there is simply no good algorithm for doing so, and as such, this ///trait is _not_ a requirement for any of the ring categorization traits. /// ///Furthermore, it is important to note that in general, while all primes are irreducible, ///the notions are **not** synonymous. In particular, in Z[sqrt(5)], 3 is irreducible ///but is not prime. However, if all pairs of elements have a [GCD], ///then all irreducibles are prime as well. pub trait Primality: Divisibility { ///Determines if this element cannot be produced from the product of two other non-invertible elements fn irreducible(&self) -> bool; /// ///Determines if this element is prime /// ///An element `x` is prime if when `y` is a multiple of `x`, `a*b=y` implies that `x` divides ///_either_ `a` _or_ `b`. Do note that this implies irreduciblity. /// fn prime(&self) -> bool; } /// ///A marker trait for [semirings](Semiring) where there are no nonzero elements that multiply to zero /// ///ie. for all `x` and `y` where `x!=0` and `y!=0`, `x*y!=0` /// pub trait NoZeroDivisors {} ///A trait for finding the greatest common divisor and least common multiple of two elements pub trait GCD: Sized { /// ///Finds the greatest common divisor of two elements /// ///An element `d` is a GCD of `x` and `y` if for all `z` where `z` divides both `x` and `y`, ///`z` also divides `d`. /// ///Do note that this element will often not be _entirely_ unique. However, it _is_ the case that ///any two GCDs in a [GCDDomain] will always differ by only multiplication by an invertible element. ///For example, while _both_ -4 and 4 are GCDs of 12 and 8 (as Integers), -4 can be ///transformed into 4 and vice versa by multiplying by -1 /// fn gcd(self, rhs: Self) -> Self; /// ///Finds the least common multiple of two elements /// ///An element `m` is a LCM of `x` and `y` if for all `z` where `x` and `y` both divide `z`, ///`m` also divides `z`. This element is also not always unique in the same way as the gcd. /// ///Additionally, the GCD and LCM can always be computed from each other ///by LCM(x,y) = (x*y)/GCD(x,y), so the existence of one always implies the existence of ///the other /// fn lcm(self, rhs: Self) -> Self; } ///A trait for finding the Bezout coefficients of a pair of elements pub trait Bezout: GCD { /// ///Bezout coefficients of a pair of elements /// ///The Bezout coefficients for `x` and `y` are a pair of elements `a` and `b` such that ///`a*x + b*y = GCD(x,y)`. Note that like with the [GCD], these are only unique up to units. ///However, the coefficients returned by this function must satisfy the defining identity for ///the _particular_ GCD as returned by [GCD::gcd()] /// #[inline] fn bezout_coefficients(self, rhs: Self) -> (Self, Self) { let (a,b,_) = self.bezout_with_gcd(rhs); (a,b) } /// ///Computes the [GCD](GCD::gcd()) and the [bezout coefficients](Bezout::bezout_coefficients) ///in one function call /// fn bezout_with_gcd(self, rhs: Self) -> (Self, Self, Self); } /// ///A marker trait for [semirings](Semiring) where each element's set of irreducible divisors is unique /// ///Note that this trait is independent of [`Factorizable`] and doesn't contain its own ///`factors()` method since there are a number of notable examples of rings ///where unique factorization is provable, but no known algorithm to find the factors is known. ///This is the case for integer polynomials for example. /// pub trait UniquelyFactorizable: Sized {} /// ///For [semirings](Semiring) that have a _known_ algorithm for decomposing elements into a set ///of irreducible factors /// ///# Notes on "Irreduciblity" /// ///An irreducible factor as defined here is any element that cannot be written as a product of two ///or more non-invertible elements. Note that this is technically, but not usually, different than ///a "prime" element. For more information, see [`Divisibility`]. /// ///# Uniqueness and Ordering /// ///This trait on its own makes **no guarantees** that the returned factors are unique or ordered /// in **any** way. The only requirements are that: ///* The product of the factors must be the original number ///* No factors of `1` are included ///* `0` returns `[0]` ///* Only *one* invertible element can be included and must be the first factor /// * For example, `-6` might return `[-1, 2, 3]` or `[-1, 3, 2]`, but not `[2, -1, 3]` /// ///Importantly, this means does that some structures **can** return different sets of factors ///at different times. For example, in ℤ[√5], `6 = 2*3 = (1+√5)*(1-√5)` are both valid irreducible ///factorizations. /// ///_However_, if the [`UniquelyFactorizable`] trait is _also_ implemented, this is restricted somewhat, ///but only up to ordering and units. /// ///In particular, given two lists of factors for an element: ///* There is still no guaranteed ordering, so the lists could be permuted somehow ///* Some factors in the second might be multiplied by some invertible element from the first /// or return an extra invertible element in the list. /// * For instance, `-10` might return `[-5, 2]` or `[5,-2]` or `[2,-5]` or even `[-1,5,2]` /// ///Furthermore, there is no requirement for any sort of ordering, as the wide range of possible ///implementing structures gives no obvious standard as to what that would be. /// ///_However_, implementors are highly encouraged to guarantee some sort of ordering and uniqueness ///themselves using their particular properties. For instance, the implementation on the primitive integers ///return all non-negative factors in order of increasing absolute magnitude with a `-1` at the beginning ///for negative numbers, or a polynomial could factor as monic polynomials of increasing degree. /// pub trait Factorizable: Sized { type Factors: Iterator<Item=Self>; fn factors(self) -> Self::Factors; } ///A trait for performing division with with remainder pub trait EuclideanDiv: Sized { ///The specific type used for the Euclidean Norm type Naturals: Natural; /// ///A measure of the "degree" of a given element. /// ///This function must satisfy the following restrictions: /// * `a.euclid_norm() <= (a*b).euclid_norm()` for all `a`,`b` /// * for `(q,r) = div_alg()`, `r.euclid_norm() < q.euclid_norm()` /// ///Beyond these, there are no other restrictions. In particular, the Euclidean Norm need not be unique. ///For example, the Euclidean norm of a polynomial is usually taken as the non-unique degree. /// fn euclid_norm(&self) -> Self::Naturals; ///The quotient from [Euclidean division](EuclideanDiv::div_alg) fn div_euc(self, rhs: Self) -> Self; ///The remainder from the [Euclidean division](EuclideanDiv::div_alg) fn rem_euc(self, rhs: Self) -> Self; /// ///Divides a pair of elements with remainder /// ///Given inputs `a` and `b`, this produces a quotient-remainder pair `(q,r)` such that: /// * `a = q*b + r` /// * `r.euclid_norm() < q.euclid_norm()` /// ///Do note that in general, this result is **not** unique, the most striking example ///simply being the Integers, for which **every** division with a remainder has the choice ///of a positive or negative remainder. (Take `13 = 4*3 + 1 = 5*3 - 2` for example) Furthermore, ///other rings, like the Guassian Integers, can have even more options, and others, like polynomials, ///can have _infinite_ possible options. /// ///As such, it is up to the implementor to decide what the canonical form of this result is ///and to communicate it as such, and this trait makes no guarantees that this has happened /// fn div_alg(self, rhs: Self) -> (Self, Self); } /// ///A unary operation mapping addition into multiplication /// ///Rigourously, and exponential operation is a mapping `E:R->R` such that: /// * `E(x+y) = E(x)*E(y)` whenever `x*y = y*x` /// * `E(x) != 0` /// ///In addition to these base properties, this trait also stipulates that this function _not_ map ///every element to 1, so as to rule out the trivial case. /// ///Furthmore, to clear up ambiguity, it is important to note that different variations on this ///definition exists. For instance: /// * As already mentioned, some may allow the mapping to be non-trivial /// * Some may also allow `E(x) = 0` /// * For the [Real](crate::analysis::Real) and [Complex](crate::analysis::Complex) exponentials, /// there are a [multitude][1] of equivalent definitions that have no true effect on the mapping /// ///More importantly though, most authors specify that the `E(x+y) = E(x)*E(y)` for _all_ `x` and `y` ///(such as on [Wikipedia][2]) ///However, doing so disallows both the Matrix exponential and Quaternion exponential as well ///as the Clifford Algebra exponential, all of which are, frankly, the only reason to make the exponential ///its own separate trait. /// /// # Effects on Ring structure /// ///It is worth noting that _any_ ring that has a non-trivial exponential operation must automatically ///have a characteristic of 0 (that is, `1+1+1+...+1` will never equal zero) and hence, has an ///embedding of the Integers within it. /// ///This fact is easily proven as follows: /// * assume `char(R) = n != 0` /// * then, for any `x` in `R`, `nx = x+...+x = 0`, and, /// the [Frobenious automorphism][3] gives that `(x+y)ⁿ = xⁿ + yⁿ` /// * Hence, `(E(x) - 1)ⁿ = E(x)ⁿ - 1 = E(nx) - 1 = E(0) - 1 = 0` /// * Thus, we have that that `E(x) = 1` contradicting our assumtions /// ///Additionally, any element that is the output of the expoenential _must_ be an invertible element, ///since <br> `1 = E(0) = E(x-x) = E(x)*E(-x)` implying `E(x) = E(x)⁻¹` /// /// # Uniqueness /// ///In general, this characterization of the exponential function is *not* unique. However, in the ///vast majority of cases, there is a canonical version that all others derive from _or_ there is only ///one non-trivial case. /// ///For example, most real-algebras have infinitely many exponentials, but we get a canonical form ///stipulating that the function satisfy the classic differential equation `E'(x) = E(x)` or some ///variant. /// /// [1]: https://en.wikipedia.org/wiki/Characterizations_of_the_exponential_function /// [2]: https://en.wikipedia.org/wiki/Exponential_field /// [3]: https://en.wikipedia.org/wiki/Frobenius_endomorphism /// pub trait Exponential: Sized { /// ///The exponential of this ring element /// ///Here, `exp(x)` is defined such that: /// * `exp(x+y) = exp(x)*exp(y)` for all `x` and `y` where `x*y=y*x` /// * `exp(x) != 0` /// * `exp(x)` is continuous (if applicable) /// * `d/dx exp(x)|ₓ₌₁ = 1` (if applicable) /// ///For most structures, this function is equivalent to the infinite series Σ x<sup>n</sup>/n! /// fn exp(self) -> Self; /// ///An inverse of [exp(x)](Exponential::exp) where `ln(1) = 0` /// ///This returns a `None` value whenever the inverse does not exist for the given input. /// /// ## Uniqueness and Continuity /// ///Do note, however, that for almost all non-[Real](crate::analysis::Real) structures, this function ///is not unique and can **never** be continuous. Of course, some of this ambiguity is resolved by ///stipulating that `ln(1) = 0`, but even so, some remains, ///and so, it is entirely up to the implementor to any specific canonical form if applicable. /// ///For example, the [Complex](crate::analysis::Complex) numbers, the natural logarithm *must* be discontinuous somewhere, ///and there are infinitely many choices as to where that is. However, usually, this ambiguity ///is removed by taking the imaginary component of the result between -π and π and setting ///the discontinuity to be on the negative real axis /// fn try_ln(self) -> Option<Self>; } ///A commutative and additive monoid with a distributive and associative multiplication operation pub trait Semiring = Distributive + AddMonoid + AddCommutative + MulSemigroup; ///A semiring with an identity element pub trait UnitalSemiring = Semiring + MulMonoid; ///A unital semiring where multiplication is commutative pub trait CommutativeSemiring = UnitalSemiring + MulCommutative; ///A semiring with a multiplicative inverse pub trait DivisionSemiring = UnitalSemiring + MulGroup; ///An additive abelian group with a distributive and associative multiplication operation pub trait Ring = Distributive + AddAbelianGroup + MulSemigroup; ///A ring with an identity element pub trait UnitalRing = Ring + MulMonoid; ///A unital ring where multiplication is commutative pub trait CommutativeRing = UnitalRing + MulCommutative; ///A ring with a multiplicative inverse pub trait DivisionRing = UnitalRing + MulGroup; ///A ring with an exponential operation pub trait ExponentialRing = UnitalRing + Exponential; ///A unital semiring with no pairs of nonzero elements that multiply to zero pub trait Semidomain = UnitalSemiring + Divisibility + NoZeroDivisors; ///A semidomain that is commutative pub trait IntegralSemidomain = Semidomain + MulCommutative; ///An integral semidomain where every pair of elements has a greatest common divisor pub trait GCDSemidomain = IntegralSemidomain + GCD; ///A GCD semidomain where every pair of elements is uniquely factorizable into irreducible elements (up to units) pub trait UFSemidomain = GCDSemidomain + UniquelyFactorizable; ///A UF semidomain with a division algorithm for dividing with a remainder pub trait EuclideanSemidomain = UFSemidomain + EuclideanDiv; ///A unital ring with no pairs of nonzero elements that multiply to zero pub trait Domain = UnitalRing + Divisibility + NoZeroDivisors; ///A domain that is commutative pub trait IntegralDomain = Domain + MulCommutative; ///A commutative ring where every pair of elements has a greatest common divisor pub trait GCDDomain = IntegralDomain + GCD; ///A commutative ring where every pair of elements has a weighted sum to their GCD pub trait BezoutDomain = GCDDomain + Bezout; ///A commutative ring that is uniquely factorizable into irreducible (up to units) pub trait UFD = GCDDomain + UniquelyFactorizable; /// ///An integral domain where every ideal is generated by one element /// ///ie. a UFD that is Bezout /// pub trait PID = UFD + BezoutDomain; ///A commutative ring with a division algorithm for dividing with a remainder pub trait EuclideanDomain = PID + EuclideanDiv; ///A set that is both an additive and multiplicative abelian group where multiplication distributes pub trait Field = CommutativeRing + MulGroup; ///A field with an exponential operation pub trait ExponentialField = Field + Exponential; // //Implementation for primitives // macro_rules! impl_dist { ($($t:ty)*) => { $( impl Distributive for $t{} impl Distributive for ::core::num::Wrapping<$t>{} )* }; } impl_dist!(usize u8 u16 u32 u64 u128 isize i8 i16 i32 i64 i128 f32 f64); macro_rules! impl_for_field { ($($t:ty)*) => { $( impl Divisibility for $t { #[inline(always)] fn divides(self, _rhs: Self) -> bool {true} #[inline(always)] fn divide(self, rhs: Self) -> Option<Self> {Some(rhs / self)} #[inline(always)] fn unit(&self) -> bool {true} #[inline(always)] fn inverse(self) -> Option<Self> {Some(self.inv())} } impl NoZeroDivisors for $t {} impl UniquelyFactorizable for $t {} )* } } impl_for_field!(f32 f64); ///Uses the [Euclidean Algorithm](https://en.wikipedia.org/wiki/Euclidean_algorithm) ///to find the [GCD] of two ring elements using [division with remainder](EuclideanDiv) pub fn euclidean<T>(lhs: T, rhs: T) -> T where T:CommutativeSemiring+EuclideanDiv { fn euclid<U>(a: U, b: U) -> U where U:CommutativeSemiring+EuclideanDiv { let r = a.rem_euc(b.clone()); if r.is_zero() { b }else{ euclid(b, r) } } if lhs.is_zero() || rhs.is_zero() {return T::zero()} if lhs.euclid_norm() >= rhs.euclid_norm() { euclid(lhs, rhs) } else { euclid(rhs, lhs) } } /// ///Uses the [Extended Euclidean Algorithm](https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm) ///to find the [GCD] of two ring elements _and_ their [bezout coefficients](Bezout) using ///[division with remainder](EuclideanDiv) /// pub fn extended_euclidean<T>(lhs: T, rhs: T) -> (T, T, T) where T:CommutativeRing+EuclideanDiv { fn euclid<U>(a: U, b: U) -> (U, U, U) where U:CommutativeRing+EuclideanDiv { let (q, r) = a.div_alg(b.clone()); if r.is_zero() { (U::zero(), U::one(), b) }else{ let (x1, y1, g) = euclid(b, r); (y1.clone(), x1 - q * y1, g) } } if lhs.is_zero() || rhs.is_zero() { return (T::zero(), T::zero(), T::zero()) } if lhs.euclid_norm() >= rhs.euclid_norm() { euclid(lhs, rhs) } else { let (y, x, g) = euclid(rhs, lhs); (x, y, g) } } /// ///Determines if a given Natural number is [prime](Primality) using the ///[Miller-Rabin primality test](https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test) /// ///The algorithm works in essence by searching for counter-examples to Fermat's Little theorem, ie, ///that `a^(p-1) = 1` for all `a` in `Z/pZ` when `p` is prime. In general, this tactic only gives a ///way to prove a number is _not_ prime, but with a few modifications to the check and by picking ///the right examples, it is possible to turn this into a deterministic test. \* /// ///Furthermore, this particular algorithm has the surprising benefit of having a runtime that is ///polynomial in the number of bits in the input. Of course, this does not guarantee that this function ///is particularly "fast" per se, but in testing, the algorithm runs reasonable fast for all primitive ///integer types. /// ///\* It is important to note that the particular method used to achieve a Deterministic Miller Robin ///algorithm in polynomial time, _technically_ depends on the Generalized Riemann Hypothesis. So I ///guess that means that for super huge numbers, this technically _could_ give a false positive... ¯\\\_(ツ)\_/¯ ///But hey, what _else_ is there? The AKS Primality Test? /// pub fn miller_rabin<Z:Natural>(n:Z) -> bool { //trivial cases if n <= Z::one() { return false; } if n == Z::two() { return true; } if n.even() { return false; } //decompose n-1 (ie. -1 in Z/nZ) into a product of the form d*2^s let mut d = n.clone() - Z::one(); let mut s = Z::zero(); while d.even() { s = s + Z::one(); d = d.div_two(); } #[inline] fn witness<N:Natural>(witness:u128, d:N, s:N, n:N) -> bool { _witness(N::try_from(witness).unwrap_or(N::zero()), d, s, n) } fn _witness<N:Natural>(mut a:N, mut d:N, mut s:N, n:N) -> bool { let mut r = a.clone(); d = d - N::one(); while d > N::zero() { if d.even() { d = d.div_two(); a = (a.clone() * a) % n.clone(); } else { d = d - N::one(); r = (r * a.clone()) % n.clone(); } } if r.is_one() { true } else { loop { if r == n.clone() - N::one() { return true; } if s.is_zero() { break; } else { s = s - N::one(); r = (r.clone() * r) % n.clone(); } } false } } //the breakpoints for each set of sufficient witnesses let b1 = Z::from_u32(2047u32); let b2 = Z::from_u32(1373653u32); let b3 = Z::from_u32(9080191u32); let b4 = Z::from_u32(25326001u32); let b5 = Z::from_u64(3215031751u64); let b6 = Z::from_u64(4759123141u64); let b7 = Z::from_u64(1122004669633u64); let b8 = Z::from_u64(2152302898747u64); let b9 = Z::from_u64(3474749660383u64); let b10 = Z::from_u64(341550071728321u64); let b11 = Z::from_u64(3825123056546413051u64); let b12 = Z::from_u128(318665857834031151167461u128); let b13 = Z::from_u128(3317044064679887385961981u128); if b1.map_or(true, |x| n < x) { witness(2, d.clone(), s.clone(), n.clone()) } else if b2.map_or(true, |x| n < x) { witness(2, d.clone(), s.clone(), n.clone()) && witness(3, d.clone(), s.clone(), n.clone()) } else if b3.map_or(true, |x| n < x) { witness(31, d.clone(), s.clone(), n.clone()) && witness(73, d.clone(), s.clone(), n.clone()) } else if b4.map_or(true, |x| n < x) { witness(2, d.clone(), s.clone(), n.clone()) && witness(3, d.clone(), s.clone(), n.clone()) && witness(5, d.clone(), s.clone(), n.clone()) } else if b5.map_or(true, |x| n < x) { witness(2, d.clone(), s.clone(), n.clone()) && witness(3, d.clone(), s.clone(), n.clone()) && witness(5, d.clone(), s.clone(), n.clone()) && witness(7, d.clone(), s.clone(), n.clone()) } else if b6.map_or(true, |x| n < x) { witness(2, d.clone(), s.clone(), n.clone()) && witness(7, d.clone(), s.clone(), n.clone()) && witness(61, d.clone(), s.clone(), n.clone()) } else if b7.map_or(true, |x| n < x) { witness(2, d.clone(), s.clone(), n.clone()) && witness(13, d.clone(), s.clone(), n.clone()) && witness(23, d.clone(), s.clone(), n.clone()) && witness(1662803, d.clone(), s.clone(), n.clone()) } else if b13.map_or(true, |x| n < x) { witness(2, d.clone(), s.clone(), n.clone()) && witness(3, d.clone(), s.clone(), n.clone()) && witness(5, d.clone(), s.clone(), n.clone()) && witness(7, d.clone(), s.clone(), n.clone()) && witness(11, d.clone(), s.clone(), n.clone()) && if b8.map_or(false, |x| n >= x) { witness(13, d.clone(), s.clone(), n.clone()) } else {true} && if b9.map_or(false, |x| n >= x) { witness(17, d.clone(), s.clone(), n.clone()) } else {true} && if b10.map_or(false, |x| n >= x) { witness(19, d.clone(), s.clone(), n.clone()) && witness(23, d.clone(), s.clone(), n.clone()) } else {true} && if b11.map_or(false, |x| n >= x) { witness(29, d.clone(), s.clone(), n.clone()) && witness(31, d.clone(), s.clone(), n.clone()) && witness(37, d.clone(), s.clone(), n.clone()) } else {true} && if b12.map_or(false, |x| n >= x) { witness(41, d.clone(), s.clone(), n.clone()) } else {true} } else { //in general, we need to check every witness below 2*ln(n)^2 let mut a = Z::two(); while Z::try_from(3.pow_n(a.clone().div_two())).unwrap_or(n.clone()) < n.clone().mul_two() { if !_witness(a.clone(), d.clone(), s.clone(), n.clone()) {return false} a = a + Z::one(); } true } } #[cfg(test)] mod tests { use crate::algebra::*; //TODO: add tests for euclidean and extended_euclidean #[test] fn primality() { assert!(18446744073709551557u64.prime()); assert!(!18446744073709551559u64.prime()); assert!(!18446744073709551555u64.prime()); assert!(2147483647u32.prime()); assert!(!2147483649u32.prime()); assert!(!2147483645u32.prime()); assert!(65521u16.prime()); assert!(65519u16.prime()); assert!(!65523u16.prime()); assert!(251u8.prime()); assert!(!253u8.prime()); assert!(!249u8.prime()); assert!(13u8.prime()); assert!(!15u8.prime()); } #[test] fn factor() { #[cfg(feature = "std")] { assert_eq!(91u8.factors().collect::<Vec<_>>(), vec![7,13] ); assert_eq!((-91i8).factors().collect::<Vec<_>>(), vec![-1,7,13] ); assert_eq!(360u16.factors().collect::<Vec<_>>(), vec![2,2,2,3,3,5] ); assert_eq!((-360i16).factors().collect::<Vec<_>>(), vec![-1,2,2,2,3,3,5] ); assert_eq!(1971813u32.factors().collect::<Vec<_>>(), vec![3,17,23,41,41] ); assert_eq!((-1971813i32).factors().collect::<Vec<_>>(), vec![-1,3,17,23,41,41] ); assert_eq!(1971813u32.factors().collect::<Vec<_>>(), vec![3,17,23,41,41] ); assert_eq!((-1971813i32).factors().collect::<Vec<_>>(), vec![-1,3,17,23,41,41] ); assert_eq!(0x344CAF57AB24A9u64.factors().collect::<Vec<_>>(), vec![101,101,103,103,107,107,109,109]); assert_eq!((-0x344CAF57AB24A9i64).factors().collect::<Vec<_>>(), vec![-1,101,101,103,103,107,107,109,109]); } fn factors_slice<Z:IntegerSubset+Factorizable>(x:Z, dest: &mut [Z]) -> usize { let mut i = 0; for f in x.factors() { if i<dest.len() { dest[i] = f; i += 1; } else { return i; } } return i; } let mut factors = [0xFF; 3]; //test 0 returns 0 assert_eq!(factors_slice(0u8, &mut factors), 1); assert_eq!(&factors, &[0,0xFF,0xFF]); //test 1 returns 1 assert_eq!(factors_slice(1u8, &mut factors), 0); assert_eq!(&factors, &[0,0xFF,0xFF]); //test the algorithm stopping at the end of the array assert_eq!(factors_slice(210u8, &mut factors), 3); assert_eq!(&factors, &[2,3,5]);//skips 7 let mut factors = [0;10]; //test -1 giving -1 assert_eq!(factors_slice(-1i64, &mut factors), 1); assert_eq!(&factors, &[-1,0,0,0,0,0,0,0,0,0]); assert_eq!(factors_slice(-0x344CAF57AB24A9i64, &mut factors), 9); assert_eq!(&factors, &[-1,101,101,103,103,107,107,109,109,0]); assert_eq!(factors_slice(0x344CAF57AB24A9i64, &mut factors), 8); assert_eq!(&factors, &[101,101,103,103,107,107,109,109,109,0]); } }