gauss.rs
24.7 kB · rust · 692 lines
1use crate::prime::is_prime;2use crate::spiral::flags;34const ROOT3: f64 = 1.732_050_807_568_877_2;56/// The two rings of whole numbers in the plane, each a pair (a, b) on its own lattice.7#[derive(Clone, Copy, Debug, PartialEq, Eq)]8pub enum Ring {9 /// a + b i on the square lattice: norm a^2 + b^2, four units, the window a square.10 Gaussian,11 /// a + b omega on the hexagonal lattice, omega a cube root of one: norm a^2 - a b + b^2, six units, the window a hexagon.12 Eisenstein,13}1415impl Ring {16 /// Reads a ring from its name.17 pub fn named(name: &str) -> Option<Ring> {18 match name {19 "gaussian" => Some(Ring::Gaussian),20 "eisenstein" => Some(Ring::Eisenstein),21 _ => None,22 }23 }24 /// Returns the norm of a point: its squared length.25 pub fn norm(self, a: i64, b: i64) -> u64 {26 match self {27 Ring::Gaussian => (a * a + b * b) as u64,28 Ring::Eisenstein => (a * a - a * b + b * b) as u64,29 }30 }31 /// Returns the product of two points.32 pub fn mul(self, (a, b): (i64, i64), (c, d): (i64, i64)) -> (i64, i64) {33 match self {34 Ring::Gaussian => (a * c - b * d, a * d + b * c),35 Ring::Eisenstein => (a * c - b * d, a * d + b * c - b * d),36 }37 }38 /// Returns the point turned anticlockwise by one unit: a quarter turn or a sixth.39 pub fn turn(self, a: i64, b: i64) -> (i64, i64) {40 match self {41 Ring::Gaussian => (-b, a),42 Ring::Eisenstein => (a - b, a),43 }44 }45 /// Returns the count of units: 4 or 6.46 pub fn units(self) -> usize {47 match self {48 Ring::Gaussian => 4,49 Ring::Eisenstein => 6,50 }51 }52 /// Returns the order of the symmetry of the picture, the units and the mirror: 8 or 12.53 pub fn symmetry(self) -> usize {54 2 * self.units()55 }56 /// Returns the unit multiples of a point, the point first, turning anticlockwise.57 pub fn associates(self, a: i64, b: i64) -> Vec<(i64, i64)> {58 let mut out = Vec::with_capacity(self.units());59 let mut at = (a, b);60 for _ in 0..self.units() {61 out.push(at);62 at = self.turn(at.0, at.1);63 }64 out65 }66 /// Returns the conjugate: the mirror image in the real axis.67 pub fn conjugate(self, a: i64, b: i64) -> (i64, i64) {68 match self {69 Ring::Gaussian => (a, -b),70 Ring::Eisenstein => (a - b, -b),71 }72 }73 /// Returns the quotient and the remainder of a point by a nonzero point: `z = q w + r` with the norm of `r` below the norm of `w`.74 pub fn div_rem(self, z: (i64, i64), w: (i64, i64)) -> ((i64, i64), (i64, i64)) {75 let n = self.norm(w.0, w.1) as i64;76 let p = self.mul(z, self.conjugate(w.0, w.1));77 let q = (78 (2 * p.0 + n).div_euclid(2 * n),79 (2 * p.1 + n).div_euclid(2 * n),80 );81 let s = self.mul(q, w);82 (q, (z.0 - s.0, z.1 - s.1))83 }84 /// Returns the canonical associate of a point: the one with `a > 0` and `b >= 0` on the square lattice, the one with `a > 0` and `0 <= b < a` on the hexagonal, the origin for the origin.85 pub fn canon(self, a: i64, b: i64) -> (i64, i64) {86 self.associates(a, b)87 .into_iter()88 .find(|&(x, y)| match self {89 Ring::Gaussian => x > 0 && y >= 0,90 Ring::Eisenstein => x > 0 && y >= 0 && y < x,91 })92 .unwrap_or((0, 0))93 }94 /// Returns the greatest common divisor of two points as its canonical associate, by the nearest-point Euclidean algorithm, the origin for two origins.95 ///96 /// ```97 /// use mrlynum::gauss::Ring;98 /// assert_eq!(Ring::Gaussian.gcd((5, 0), (2, 1)), (2, 1));99 /// assert_eq!(Ring::Gaussian.gcd((3, 0), (0, 7)), (1, 0));100 /// ```101 pub fn gcd(self, z: (i64, i64), w: (i64, i64)) -> (i64, i64) {102 let (mut z, mut w) = (z, w);103 while w != (0, 0) {104 let (_, r) = self.div_rem(z, w);105 z = w;106 w = r;107 }108 self.canon(z.0, z.1)109 }110 /// Returns the whole number an associate of the point lies on, when one lies on the positive real axis.111 pub fn whole(self, a: i64, b: i64) -> Option<u64> {112 self.associates(a, b)113 .into_iter()114 .find(|&(x, y)| y == 0 && x > 0)115 .map(|(x, _)| x as u64)116 }117 /// Returns the one rational prime that ramifies: 2 or 3.118 pub fn ramified(self) -> u64 {119 match self {120 Ring::Gaussian => 2,121 Ring::Eisenstein => 3,122 }123 }124 /// Returns whether a rational prime stays prime in the ring: 3 mod 4, or 2 mod 3.125 pub fn inert(self, p: u64) -> bool {126 match self {127 Ring::Gaussian => p % 4 == 3,128 Ring::Eisenstein => p % 3 == 2,129 }130 }131 /// Returns the fate of a whole number as a prime of the ring: split, inert or ramified, unit for one, zero for zero, composite otherwise.132 pub fn fate(self, n: u64) -> Class {133 match n {134 0 => Class::Zero,135 1 => Class::Unit,136 _ if !is_prime(n as usize) => Class::Composite,137 _ if n == self.ramified() => Class::Ramified,138 _ if self.inert(n) => Class::Inert,139 _ => Class::Split,140 }141 }142 /// Returns the reach of a point: the ring of the window it sits on, the Chebyshev distance or the hex distance.143 pub fn reach(self, a: i64, b: i64) -> u64 {144 match self {145 Ring::Gaussian => a.abs().max(b.abs()) as u64,146 Ring::Eisenstein => a.abs().max(b.abs()).max((a - b).abs()) as u64,147 }148 }149 /// Returns the count of points within the reach: the square or the hexagon.150 pub fn count(self, radius: u64) -> usize {151 let r = radius as usize;152 match self {153 Ring::Gaussian => (2 * r + 1) * (2 * r + 1),154 Ring::Eisenstein => 3 * r * r + 3 * r + 1,155 }156 }157 /// Returns the largest norm within the reach: 2 r^2 at the square's corner, r^2 at the hexagon's.158 pub fn top(self, radius: u64) -> u64 {159 match self {160 Ring::Gaussian => 2 * radius * radius,161 Ring::Eisenstein => radius * radius,162 }163 }164 /// Returns the place of a point in the plane, x right and y up, one unit between neighbours.165 pub fn place(self, a: i64, b: i64) -> (f64, f64) {166 let (a, b) = (a as f64, b as f64);167 match self {168 Ring::Gaussian => (a, b),169 Ring::Eisenstein => (a - b / 2.0, b * ROOT3 / 2.0),170 }171 }172 /// Returns the point nearest a place in the plane.173 pub fn nearest(self, x: f64, y: f64) -> (i64, i64) {174 match self {175 Ring::Gaussian => (x.round() as i64, y.round() as i64),176 Ring::Eisenstein => {177 let v = -2.0 * y / ROOT3;178 let u = x - v / 2.0;179 let w = -u - v;180 let (mut ru, mut rv, rw) = (u.round(), v.round(), w.round());181 let (du, dv, dw) = ((ru - u).abs(), (rv - v).abs(), (rw - w).abs());182 if du > dv && du > dw {183 ru = -rv - rw;184 } else if dv > dw {185 rv = -ru - rw;186 }187 (ru as i64, -rv as i64)188 }189 }190 }191}192193/// What a point of the ring is.194#[derive(Clone, Copy, Debug, PartialEq, Eq)]195pub enum Class {196 /// The origin.197 Zero,198 /// A unit: norm one.199 Unit,200 /// A prime over the one rational prime that ramifies, 2 or 3.201 Ramified,202 /// A prime whose norm is a rational prime that splits into it and its conjugate.203 Split,204 /// A rational prime that stays prime in the ring, times a unit.205 Inert,206 /// A product of two points of norm above one.207 Composite,208}209210impl Class {211 /// Returns the class as a word.212 pub fn word(self) -> &'static str {213 match self {214 Class::Zero => "zero",215 Class::Unit => "unit",216 Class::Ramified => "ramified",217 Class::Split => "split",218 Class::Inert => "inert",219 Class::Composite => "composite",220 }221 }222 /// Returns whether the class is prime.223 pub fn prime(self) -> bool {224 matches!(self, Class::Ramified | Class::Split | Class::Inert)225 }226}227228/// The tallies of a window: every class counted and the share of primes.229#[derive(Clone, Debug, PartialEq)]230pub struct Census {231 /// The count of points.232 pub points: usize,233 /// The count of primes.234 pub primes: usize,235 /// The split primes.236 pub split: usize,237 /// The inert primes.238 pub inert: usize,239 /// The ramified primes.240 pub ramified: usize,241 /// The units.242 pub units: usize,243 /// The composites.244 pub composites: usize,245 /// The primes over the points.246 pub density: f64,247}248249/// The symmetric window of one ring: every point within a reach, with the norms sieved once.250#[derive(Clone, Debug)]251pub struct Window {252 ring: Ring,253 radius: u64,254 prime: Vec<bool>,255}256257impl Window {258 /// Opens the window of a ring out to a reach, sieving every norm inside it.259 pub fn new(ring: Ring, radius: u64) -> Window {260 Window {261 ring,262 radius,263 prime: flags(ring.top(radius) as usize),264 }265 }266 /// Returns the ring.267 pub fn ring(&self) -> Ring {268 self.ring269 }270 /// Returns the reach.271 pub fn radius(&self) -> u64 {272 self.radius273 }274 /// Returns whether a point lies inside.275 pub fn holds(&self, a: i64, b: i64) -> bool {276 self.ring.reach(a, b) <= self.radius277 }278 fn is_prime(&self, n: u64) -> bool {279 match self.prime.get(n as usize) {280 Some(&p) => p,281 None => is_prime(n as usize),282 }283 }284 /// Classifies a point: prime when its norm is a rational prime, or when it is a unit times a rational prime that stays prime.285 ///286 /// ```287 /// use mrlynum::gauss::{Class, Ring, Window};288 /// let window = Window::new(Ring::Gaussian, 3);289 /// assert_eq!(window.class(2, 1), Class::Split);290 /// assert_eq!(window.class(0, -3), Class::Inert);291 /// ```292 pub fn class(&self, a: i64, b: i64) -> Class {293 let n = self.ring.norm(a, b);294 if n == 0 {295 return Class::Zero;296 }297 if n == 1 {298 return Class::Unit;299 }300 if self.is_prime(n) {301 return if n == self.ring.ramified() {302 Class::Ramified303 } else {304 Class::Split305 };306 }307 match self.ring.whole(a, b) {308 Some(p) if self.ring.inert(p) && self.is_prime(p) => Class::Inert,309 _ => Class::Composite,310 }311 }312 /// Lists every point inside, row by row from the bottom left of the bounding square.313 pub fn points(&self) -> Vec<(i64, i64)> {314 let r = self.radius as i64;315 let mut out = Vec::with_capacity(self.ring.count(self.radius));316 for b in -r..=r {317 for a in -r..=r {318 if self.holds(a, b) {319 out.push((a, b));320 }321 }322 }323 out324 }325 /// Counts every class inside.326 pub fn census(&self) -> Census {327 let mut census = Census {328 points: 0,329 primes: 0,330 split: 0,331 inert: 0,332 ramified: 0,333 units: 0,334 composites: 0,335 density: 0.0,336 };337 for (a, b) in self.points() {338 census.points += 1;339 match self.class(a, b) {340 Class::Split => census.split += 1,341 Class::Inert => census.inert += 1,342 Class::Ramified => census.ramified += 1,343 Class::Unit => census.units += 1,344 Class::Composite => census.composites += 1,345 Class::Zero => {}346 }347 }348 census.primes = census.split + census.inert + census.ramified;349 census.density = census.primes as f64 / census.points as f64;350 census351 }352}353354/// Counts the points of every norm from zero through the limit, by enumeration: the ring weights of the lattice.355///356/// ```357/// assert_eq!(mrlynum::gauss::shells(mrlynum::gauss::Ring::Gaussian, 5), vec![1, 4, 4, 0, 4, 8]);358/// ```359pub fn shells(ring: Ring, limit: usize) -> Vec<u32> {360 let mut out = vec![0u32; limit + 1];361 let reach = (4 * limit / 3).isqrt() as i64 + 1;362 for a in -reach..=reach {363 for b in -reach..=reach {364 let n = ring.norm(a, b) as usize;365 if n <= limit {366 out[n] += 1;367 }368 }369 }370 out371}372373/// Returns the norm from one through the limit with the most points and that count, the earliest on a tie.374///375/// ```376/// assert_eq!(mrlynum::gauss::peak(mrlynum::gauss::Ring::Gaussian, 60), (25, 12));377/// ```378pub fn peak(ring: Ring, limit: usize) -> (usize, u32) {379 shells(ring, limit)380 .into_iter()381 .enumerate()382 .skip(1)383 .fold(384 (0, 0),385 |best, (n, count)| if count > best.1 { (n, count) } else { best },386 )387}388389/// Lists one point per associate class of the nonzero points of norm at most the bound: canonical associates, in order of norm and then of coordinates.390///391/// ```392/// use mrlynum::gauss::{classes, Ring};393/// assert_eq!(classes(Ring::Gaussian, 5), vec![(1, 0), (1, 1), (2, 0), (1, 2), (2, 1)]);394/// ```395pub fn classes(ring: Ring, bound: u64) -> Vec<(i64, i64)> {396 let reach = (4 * bound / 3).isqrt() as i64 + 1;397 let mut out = Vec::new();398 for a in 1..=reach {399 let wide = match ring {400 Ring::Gaussian => reach,401 Ring::Eisenstein => a - 1,402 };403 for b in 0..=wide {404 if ring.norm(a, b) <= bound {405 out.push((a, b));406 }407 }408 }409 out.sort_by_key(|&(a, b)| (ring.norm(a, b), a, b));410 out411}412413#[cfg(test)]414mod tests {415 use super::*;416 use crate::factor::divisors;417418 fn divides(ring: Ring, z: (i64, i64), w: (i64, i64)) -> bool {419 let n = ring.norm(w.0, w.1) as i64;420 let (x, y) = ring.mul(z, ring.conjugate(w.0, w.1));421 x % n == 0 && y % n == 0422 }423424 fn brute(ring: Ring, z: (i64, i64)) -> bool {425 let n = ring.norm(z.0, z.1);426 if n < 2 {427 return false;428 }429 let r = n.isqrt() as i64 + 1;430 for c in -r..=r {431 for d in -r..=r {432 let m = ring.norm(c, d);433 if m > 1 && m < n && divides(ring, z, (c, d)) {434 return false;435 }436 }437 }438 true439 }440441 #[test]442 fn the_gaussian_window_pins_its_primes() {443 let two = Window::new(Ring::Gaussian, 2).census();444 assert_eq!((two.points, two.primes), (25, 12));445 assert_eq!((two.ramified, two.split, two.inert), (4, 8, 0));446 let three = Window::new(Ring::Gaussian, 3).census();447 assert_eq!((three.points, three.primes, three.units), (49, 24, 4));448 assert_eq!((three.ramified, three.split, three.inert), (4, 16, 4));449 assert_eq!(three.composites, 49 - 24 - 4 - 1);450 assert!((three.density - 24.0 / 49.0).abs() < 1e-12);451 }452453 #[test]454 fn the_eisenstein_window_pins_its_primes() {455 let two = Window::new(Ring::Eisenstein, 2).census();456 assert_eq!((two.points, two.primes, two.units), (19, 12, 6));457 assert_eq!((two.ramified, two.split, two.inert), (6, 0, 6));458 let three = Window::new(Ring::Eisenstein, 3).census();459 assert_eq!((three.points, three.primes), (37, 24));460 assert_eq!((three.ramified, three.split, three.inert), (6, 12, 6));461 }462463 #[test]464 fn the_classes_follow_the_norm_rules() {465 let g = Window::new(Ring::Gaussian, 5);466 assert_eq!(g.class(0, 0), Class::Zero);467 assert_eq!(g.class(-1, 0), Class::Unit);468 assert_eq!(g.class(1, 1), Class::Ramified);469 assert_eq!(g.class(2, 1), Class::Split);470 assert_eq!(g.class(3, 0), Class::Inert);471 assert_eq!(g.class(0, -3), Class::Inert);472 assert_eq!(g.class(5, 0), Class::Composite);473 assert_eq!(g.class(1, 3), Class::Composite);474 let e = Window::new(Ring::Eisenstein, 5);475 assert_eq!(e.class(1, 1), Class::Unit);476 assert_eq!(e.class(1, -1), Class::Ramified);477 assert_eq!(e.class(2, 0), Class::Inert);478 assert_eq!(e.class(2, 2), Class::Inert);479 assert_eq!(e.class(3, 0), Class::Composite);480 assert_eq!(e.class(2, -1), Class::Split);481 assert_eq!(e.class(3, 1), Class::Split);482 assert_eq!(e.class(4, 1), Class::Split);483 for ring in [Ring::Gaussian, Ring::Eisenstein] {484 let window = Window::new(ring, 9);485 for (a, b) in window.points() {486 assert_eq!(487 window.class(a, b).prime(),488 brute(ring, (a, b)),489 "{ring:?} {a} {b}"490 );491 }492 for n in 0..40 {493 let axis = window.class(n, 0);494 let fate = ring.fate(n as u64);495 let expect = match fate {496 Class::Split | Class::Ramified => Class::Composite,497 other => other,498 };499 assert_eq!(axis, expect, "{ring:?} {n}");500 }501 }502 assert_eq!(Ring::Gaussian.fate(5), Class::Split);503 assert_eq!(Ring::Gaussian.fate(7), Class::Inert);504 assert_eq!(Ring::Gaussian.fate(2), Class::Ramified);505 assert_eq!(Ring::Eisenstein.fate(7), Class::Split);506 assert_eq!(Ring::Eisenstein.fate(5), Class::Inert);507 assert_eq!(Ring::Eisenstein.fate(3), Class::Ramified);508 }509510 #[test]511 fn the_units_and_the_mirror_keep_the_norm() {512 assert_eq!(513 Ring::Gaussian.associates(2, 1),514 vec![(2, 1), (-1, 2), (-2, -1), (1, -2)]515 );516 assert_eq!(517 Ring::Eisenstein.associates(2, 0),518 vec![(2, 0), (2, 2), (0, 2), (-2, 0), (-2, -2), (0, -2)]519 );520 assert_eq!(Ring::Gaussian.conjugate(2, 1), (2, -1));521 assert_eq!(Ring::Eisenstein.conjugate(2, -1), (3, 1));522 assert_eq!(Ring::Eisenstein.conjugate(0, 1), (-1, -1));523 for ring in [Ring::Gaussian, Ring::Eisenstein] {524 assert_eq!(ring.symmetry(), 2 * ring.units());525 for a in -4..=4 {526 for b in -4..=4 {527 let n = ring.norm(a, b);528 let (c, d) = ring.conjugate(a, b);529 assert_eq!(ring.norm(c, d), n);530 assert_eq!(ring.mul((a, b), (c, d)), (n as i64, 0));531 for (x, y) in ring.associates(a, b) {532 assert_eq!(ring.norm(x, y), n);533 assert_eq!(ring.reach(x, y), ring.reach(a, b));534 }535 let (x, y) = ring.place(a, b);536 assert!((x * x + y * y - n as f64).abs() < 1e-9);537 assert_eq!(ring.nearest(x + 0.45, y), (a, b));538 assert_eq!(ring.nearest(x, y - 0.45), (a, b));539 }540 }541 }542 assert_eq!(Ring::Eisenstein.reach(1, -1), 2);543 assert_eq!(Ring::Eisenstein.reach(2, -1), 3);544 for r in 0..6 {545 let window = Window::new(Ring::Eisenstein, r);546 assert_eq!(window.points().len(), Ring::Eisenstein.count(r));547 let peak = window548 .points()549 .iter()550 .map(|&(a, b)| Ring::Eisenstein.norm(a, b))551 .max()552 .unwrap();553 assert_eq!(peak, Ring::Eisenstein.top(r));554 }555 }556557 #[test]558 fn the_shells_are_the_ring_weights() {559 let r2 = shells(Ring::Gaussian, 25);560 assert_eq!(561 (r2[0], r2[1], r2[2], r2[3], r2[5], r2[25]),562 (1, 4, 4, 0, 8, 12)563 );564 for (n, &count) in shells(Ring::Gaussian, 2000).iter().enumerate().skip(1) {565 let (d1, d3) = divisors(n).iter().fold((0, 0), |(d1, d3), &d| {566 (d1 + (d % 4 == 1) as u32, d3 + (d % 4 == 3) as u32)567 });568 assert_eq!(count, 4 * (d1 - d3), "{n}");569 }570 let hex = shells(Ring::Eisenstein, 7);571 assert_eq!(hex, vec![1, 6, 0, 6, 6, 0, 0, 12]);572 for (n, &count) in shells(Ring::Eisenstein, 300).iter().enumerate().skip(1) {573 let chi: i32 = divisors(n)574 .iter()575 .map(|&d| match d % 3 {576 1 => 1,577 2 => -1,578 _ => 0,579 })580 .sum();581 assert_eq!(count as i32, 6 * chi, "{n}");582 }583 }584585 #[test]586 fn the_peak_is_the_busiest_norm() {587 assert_eq!(peak(Ring::Gaussian, 60), (25, 12));588 assert_eq!(peak(Ring::Eisenstein, 60), (49, 18));589 assert_eq!(peak(Ring::Gaussian, 4), (1, 4));590 }591592 #[test]593 fn the_division_leaves_a_remainder_under_the_divisor() {594 for ring in [Ring::Gaussian, Ring::Eisenstein] {595 for a in -6..=6 {596 for b in -6..=6 {597 for c in -4..=4 {598 for d in -4..=4 {599 if (c, d) == (0, 0) {600 continue;601 }602 let (q, r) = ring.div_rem((a, b), (c, d));603 let s = ring.mul(q, (c, d));604 assert_eq!((s.0 + r.0, s.1 + r.1), (a, b), "{ring:?} {a} {b} {c} {d}");605 assert!(606 ring.norm(r.0, r.1) < ring.norm(c, d),607 "{ring:?} {a} {b} {c} {d}"608 );609 }610 }611 }612 }613 }614 }615616 #[test]617 fn the_gcd_swallows_every_common_divisor() {618 for ring in [Ring::Gaussian, Ring::Eisenstein] {619 for a in -4..=4 {620 for b in -4..=4 {621 for c in -4..=4 {622 for d in -4..=4 {623 if (a, b) == (0, 0) || (c, d) == (0, 0) {624 continue;625 }626 let g = ring.gcd((a, b), (c, d));627 assert_eq!(ring.canon(g.0, g.1), g, "{ring:?} {a} {b} {c} {d}");628 assert!(divides(ring, (a, b), g), "{ring:?} {a} {b} {c} {d}");629 assert!(divides(ring, (c, d), g), "{ring:?} {a} {b} {c} {d}");630 for x in -8..=8 {631 for y in -8..=8 {632 if (x, y) == (0, 0) {633 continue;634 }635 if divides(ring, (a, b), (x, y))636 && divides(ring, (c, d), (x, y))637 {638 assert!(639 divides(ring, g, (x, y)),640 "{ring:?} {a} {b} {c} {d} {x} {y}"641 );642 }643 }644 }645 }646 }647 }648 }649 }650 }651652 #[test]653 fn the_canonical_associate_is_the_only_one_of_its_class() {654 for ring in [Ring::Gaussian, Ring::Eisenstein] {655 assert_eq!(ring.canon(0, 0), (0, 0));656 for a in -6..=6 {657 for b in -6..=6 {658 if (a, b) == (0, 0) {659 continue;660 }661 let c = ring.canon(a, b);662 assert_eq!(ring.canon(c.0, c.1), c, "{ring:?} {a} {b}");663 assert_eq!(ring.norm(c.0, c.1), ring.norm(a, b), "{ring:?} {a} {b}");664 let fixed = ring665 .associates(a, b)666 .into_iter()667 .filter(|&(x, y)| ring.canon(x, y) == (x, y))668 .count();669 assert_eq!(fixed, 1, "{ring:?} {a} {b}");670 }671 }672 }673 }674675 #[test]676 fn the_classes_cover_every_shell_once_per_unit() {677 for ring in [Ring::Gaussian, Ring::Eisenstein] {678 let list = classes(ring, 200);679 let total: u32 = shells(ring, 200).into_iter().skip(1).sum();680 assert_eq!(list.len() * ring.units(), total as usize, "{ring:?}");681 for &(a, b) in &list {682 assert_eq!(ring.canon(a, b), (a, b), "{ring:?} {a} {b}");683 assert!(ring.norm(a, b) <= 200, "{ring:?} {a} {b}");684 }685 for pair in list.windows(2) {686 assert!(ring.norm(pair[0].0, pair[0].1) <= ring.norm(pair[1].0, pair[1].1));687 }688 }689 assert_eq!(classes(Ring::Gaussian, 50).len(), 40);690 assert_eq!(classes(Ring::Eisenstein, 50).len(), 31);691 }692}