shape.rs
37.1 kB · rust · 1094 lines
1use mrlycore::errors::{value_error, Result};2use mrlycore::tensor::Tensor;3use mrlynum::classics::gcd;45// FRACTIONS67/// An exact rational number with a positive, reduced denominator.8#[derive(Clone, Copy, Debug, PartialEq, Eq)]9pub struct Frac {10 /// The numerator, carrying the sign.11 pub num: i64,12 /// The denominator, always positive.13 pub den: i64,14}1516fn reduce(num: i128, den: i128) -> Frac {17 assert!(den != 0, "Frac denominator must be nonzero");18 let sign = if den < 0 { -1 } else { 1 };19 let g = (gcd(num.unsigned_abs(), den.unsigned_abs()) as i128).max(1);20 Frac {21 num: i64::try_from(sign * num / g).expect("Frac numerator overflow"),22 den: i64::try_from(sign * den / g).expect("Frac denominator overflow"),23 }24}2526impl Frac {27 /// Builds the reduced fraction num over den, panicking on a zero denominator.28 pub fn new(num: i64, den: i64) -> Frac {29 reduce(num as i128, den as i128)30 }31 /// Wraps an integer as a fraction over one.32 pub fn whole(num: i64) -> Frac {33 Frac { num, den: 1 }34 }35}3637impl std::ops::Add for Frac {38 type Output = Frac;39 /// Returns the exact sum, panicking when the reduced result overflows i64.40 fn add(self, other: Frac) -> Frac {41 reduce(42 self.num as i128 * other.den as i128 + other.num as i128 * self.den as i128,43 self.den as i128 * other.den as i128,44 )45 }46}4748impl std::ops::Sub for Frac {49 type Output = Frac;50 /// Returns the exact difference, panicking when the reduced result overflows i64.51 fn sub(self, other: Frac) -> Frac {52 reduce(53 self.num as i128 * other.den as i128 - other.num as i128 * self.den as i128,54 self.den as i128 * other.den as i128,55 )56 }57}5859impl std::ops::Mul for Frac {60 type Output = Frac;61 /// Returns the exact product, panicking when the reduced result overflows i64.62 fn mul(self, other: Frac) -> Frac {63 reduce(64 self.num as i128 * other.num as i128,65 self.den as i128 * other.den as i128,66 )67 }68}6970fn lcm(a: i64, b: i64) -> i64 {71 let g = gcd(a.unsigned_abs() as u128, b.unsigned_abs() as u128) as i128;72 i64::try_from(a as i128 / g * b as i128).expect("lcm overflow")73}7475// SHAPES7677/// A closed half-space: the points x with normal dot x at most offset.78#[derive(Clone, Debug, PartialEq, Eq)]79pub struct Half {80 /// The integer outward normal.81 pub normal: Vec<i64>,82 /// The rational offset the linear form stays under.83 pub offset: Frac,84}8586/// An exact region of the unit box, scaled onto the lattice by the side.87#[derive(Clone, Debug, PartialEq, Eq)]88pub enum Shape {89 /// The closed ball of the given rational center and radius.90 Ball {91 /// The rational center, one coordinate per axis.92 center: Vec<Frac>,93 /// The rational radius.94 radius: Frac,95 },96 /// The intersection of closed half-spaces.97 Polytope {98 /// The bounding walls.99 walls: Vec<Half>,100 },101 /// The complement of the inner shape: In and Out swap, Cut stays.102 Anti(Box<Shape>),103}104105/// Where one lattice cell sits relative to a shape.106#[derive(Clone, Copy, Debug, PartialEq, Eq)]107#[repr(u8)]108pub enum Region {109 /// The cell lies fully outside the shape.110 Out = 0,111 /// The cell crosses the shape's boundary.112 Cut = 1,113 /// The cell lies fully inside the shape.114 In = 2,115}116117impl Region {118 /// Swaps In and Out, keeping Cut.119 pub fn flip(self) -> Region {120 match self {121 Region::Out => Region::In,122 Region::Cut => Region::Cut,123 Region::In => Region::Out,124 }125 }126}127128// CLASSIFICATION129130fn classify_half(wall: &Half, side: usize, index: &[usize]) -> Region {131 let bound = wall.offset.num as i128 * side as i128;132 let den = wall.offset.den as i128;133 let mut low: i128 = 0;134 let mut high: i128 = 0;135 for (axis, &n) in wall.normal.iter().enumerate() {136 let n = n as i128;137 let i = index[axis] as i128;138 if n >= 0 {139 low += n * i;140 high += n * (i + 1);141 } else {142 low += n * (i + 1);143 high += n * i;144 }145 }146 if high * den <= bound {147 Region::In148 } else if low * den > bound {149 Region::Out150 } else {151 Region::Cut152 }153}154155fn ball_scale(center: &[Frac], radius: Frac) -> i64 {156 center.iter().fold(radius.den, |l, c| lcm(l, c.den))157}158159fn classify_ball(center: &[Frac], radius: Frac, side: usize, index: &[usize]) -> Region {160 let l = ball_scale(center, radius) as i128;161 let s = side as i128;162 let r = radius.num as i128 * (l / radius.den as i128) * s;163 if r < 0 {164 return Region::Out;165 }166 let rr = r * r;167 let mut near: i128 = 0;168 let mut far: i128 = 0;169 for (axis, c) in center.iter().enumerate() {170 let cc = c.num as i128 * (l / c.den as i128) * s;171 let lo = index[axis] as i128 * l;172 let hi = lo + l;173 let dn = cc.clamp(lo, hi) - cc;174 near += dn * dn;175 let df = (lo - cc).abs().max((hi - cc).abs());176 far += df * df;177 }178 if far <= rr {179 Region::In180 } else if near > rr {181 Region::Out182 } else {183 Region::Cut184 }185}186187/// Places one lattice cell relative to the shape, exactly, with no floats.188///189/// The cell at the index occupies the closed box from the index to the index plus one on each axis, and the shape's unit-box coordinates are scaled by the side.190///191/// A polytope is judged wall by wall: Out means some wall excludes the cell, In means every wall contains it, and Cut means neither - so a cell that touches each wall's feasible side separately reads Cut even when the wall intersection misses it, a conservative call that never mislabels In or Out.192pub fn classify(shape: &Shape, side: usize, index: &[usize]) -> Region {193 match shape {194 Shape::Ball { center, radius } => classify_ball(center, *radius, side, index),195 Shape::Polytope { walls } => {196 let mut cut = false;197 for wall in walls {198 match classify_half(wall, side, index) {199 Region::Out => return Region::Out,200 Region::Cut => cut = true,201 Region::In => {}202 }203 }204 if cut {205 Region::Cut206 } else {207 Region::In208 }209 }210 Shape::Anti(inner) => classify(inner, side, index).flip(),211 }212}213214/// Classifies every cell of the grid, packing Out, Cut and In as 0, 1 and 2; the first extent sets the lattice side.215pub fn regions(shape: &Shape, dims: &[usize]) -> Tensor {216 let side = dims.first().copied().unwrap_or(0);217 let rank = dims.len();218 let mut out = Tensor::new(dims.to_vec());219 let mut index = vec![0usize; rank];220 for flat in 0..out.size() {221 let mut rem = flat;222 for axis in (0..rank).rev() {223 index[axis] = rem % dims[axis];224 rem /= dims[axis];225 }226 out.bytes_mut()[flat] = classify(shape, side, &index) as u8;227 }228 out229}230231// CATALOG232233/// Lists the named shapes of a dimension.234pub fn shapes(dimension: usize) -> Vec<&'static str> {235 let mut out = vec!["ball", "box", "diamond"];236 if dimension == 2 {237 out.extend(["triangle", "octagon"]);238 }239 if dimension == 3 {240 out.extend(["octahedron", "tetrahedron", "pyramid"]);241 }242 out243}244245fn half(normal: Vec<i64>, offset: Frac) -> Half {246 Half { normal, offset }247}248249fn axis_normal(dimension: usize, axis: usize, sign: i64) -> Vec<i64> {250 let mut n = vec![0i64; dimension];251 n[axis] = sign;252 n253}254255fn box_shape(dimension: usize, radius: Frac) -> Shape {256 let h = Frac::new(1, 2);257 let mut walls = Vec::with_capacity(2 * dimension);258 for axis in 0..dimension {259 walls.push(half(axis_normal(dimension, axis, 1), h + radius));260 walls.push(half(axis_normal(dimension, axis, -1), radius - h));261 }262 Shape::Polytope { walls }263}264265fn diamond_walls(dimension: usize, radius: Frac) -> Vec<Half> {266 (0..1usize << dimension)267 .map(|bits| {268 let normal: Vec<i64> = (0..dimension)269 .map(|axis| if (bits >> axis) & 1 == 1 { -1 } else { 1 })270 .collect();271 let toward = Frac::new(normal.iter().sum(), 2);272 half(normal, toward + radius)273 })274 .collect()275}276277fn diamond(dimension: usize, radius: Frac) -> Shape {278 Shape::Polytope {279 walls: diamond_walls(dimension, radius),280 }281}282283fn triangle(radius: Frac) -> Shape {284 let h = Frac::new(1, 2);285 let low = radius - Frac::new(3, 2);286 Shape::Polytope {287 walls: vec![288 half(vec![1, 0], h + radius),289 half(vec![-1, -2], low),290 half(vec![-1, 2], h + radius),291 ],292 }293}294295fn octagon(radius: Frac) -> Shape {296 let cut = radius * Frac::new(3, 2);297 let mut walls = match box_shape(2, radius) {298 Shape::Polytope { walls } => walls,299 _ => unreachable!(),300 };301 walls.extend(diamond_walls(2, cut));302 Shape::Polytope { walls }303}304305fn tetrahedron(radius: Frac) -> Shape {306 let normals = [[1, 1, 1], [-1, -1, 1], [-1, 1, -1], [1, -1, -1]];307 let walls = normals308 .iter()309 .map(|n| {310 let toward = Frac::new(n.iter().sum(), 2);311 half(n.to_vec(), toward + radius)312 })313 .collect();314 Shape::Polytope { walls }315}316317fn pyramid(radius: Frac) -> Shape {318 let h = Frac::new(1, 2);319 let low = radius - Frac::new(3, 2);320 Shape::Polytope {321 walls: vec![322 half(vec![1, 0, 0], h + radius),323 half(vec![-1, -2, 0], low),324 half(vec![-1, 2, 0], h + radius),325 half(vec![-1, 0, -2], low),326 half(vec![-1, 0, 2], h + radius),327 ],328 }329}330331/// Builds a named shape of the dimension, centered at one half on every axis, or an error for an unknown or irrational name.332///333/// Regular hexagons and equilateral triangles have irrational walls in the grid frame, so they are excluded on purpose rather than approximated.334pub fn named(name: &str, dimension: usize, radius: Frac) -> Result<Shape> {335 match (name, dimension) {336 ("ball", _) => Ok(Shape::Ball {337 center: vec![Frac::new(1, 2); dimension],338 radius,339 }),340 ("box", _) => Ok(box_shape(dimension, radius)),341 ("diamond", _) => Ok(diamond(dimension, radius)),342 ("triangle", 2) => Ok(triangle(radius)),343 ("octagon", 2) => Ok(octagon(radius)),344 ("octahedron", 3) => Ok(diamond(3, radius)),345 ("tetrahedron", 3) => Ok(tetrahedron(radius)),346 ("pyramid", 3) => Ok(pyramid(radius)),347 ("hexagon", 2) | ("equilateral", 2) => {348 value_error("a regular hexagon or equilateral triangle has irrational walls on the grid; no exact crop exists.")349 }350 _ => value_error(format!("unknown shape {name} in dimension {dimension}.")),351 }352}353354// CROPPING355356/// Zeroes every cell of the design outside the shape, keeping Cut cells on request; anti-crop is Shape::Anti.357pub fn crop(types: &Tensor, shape: &Shape, keep_cut: bool) -> Tensor {358 let map = regions(shape, &types.shape);359 let mut out = types.clone();360 for flat in 0..out.size() {361 let region = map.bytes()[flat];362 let keep = region == Region::In as u8 || (keep_cut && region == Region::Cut as u8);363 if !keep {364 out.put(flat, 0);365 }366 }367 out368}369370/// The refine output ceiling in cells.371pub const REFINE_LIMIT: usize = 20_000_000;372373/// Replicates each design cell base to the extra per axis and keeps a sub-cell only where its own region passes, or an error past the cell ceiling.374pub fn refine(375 types: &Tensor,376 shape: &Shape,377 base: usize,378 extra: usize,379 keep_cut: bool,380) -> Result<Tensor> {381 if base < 1 {382 return value_error("base must be at least 1.");383 }384 let exp = match u32::try_from(extra) {385 Ok(e) => e,386 Err(_) => return value_error(format!("refine output would exceed {REFINE_LIMIT} cells.")),387 };388 let factor = match base.checked_pow(exp) {389 Some(f) => f,390 None => return value_error(format!("refine output would exceed {REFINE_LIMIT} cells.")),391 };392 let mut cells: usize = 1;393 let mut dims = Vec::with_capacity(types.shape.len());394 for &n in &types.shape {395 let grown = match n.checked_mul(factor) {396 Some(g) => g,397 None => {398 return value_error(format!("refine output would exceed {REFINE_LIMIT} cells."))399 }400 };401 cells = match cells.checked_mul(grown) {402 Some(c) if c <= REFINE_LIMIT => c,403 _ => return value_error(format!("refine output would exceed {REFINE_LIMIT} cells.")),404 };405 dims.push(grown);406 }407 let side = dims.first().copied().unwrap_or(0);408 let rank = dims.len();409 let mut out = Tensor::typed(dims.clone(), types.dtype());410 let mut index = vec![0usize; rank];411 let mut parent = vec![0usize; rank];412 for flat in 0..out.size() {413 let mut rem = flat;414 for axis in (0..rank).rev() {415 index[axis] = rem % dims[axis];416 rem /= dims[axis];417 }418 for axis in 0..rank {419 parent[axis] = index[axis] / factor;420 }421 let value = types.at(types.index(&parent));422 if value == 0 {423 continue;424 }425 let region = classify(shape, side, &index);426 let keep = region == Region::In || (keep_cut && region == Region::Cut);427 if keep {428 out.put(flat, value);429 }430 }431 Ok(out)432}433434// CENSUS435436/// The per-region tallies of a shape over a design, indexed Out, Cut, In.437#[derive(Clone, Copy, Debug, PartialEq, Eq)]438pub struct ShapeCensus {439 /// The cell count of each region.440 pub cells: [usize; 3],441 /// The filled-cell count of each region.442 pub filled: [usize; 3],443}444445/// Tallies the design's cells and filled cells per region of the shape.446pub fn census(shape: &Shape, types: &Tensor) -> ShapeCensus {447 let map = regions(shape, &types.shape);448 let mut out = ShapeCensus {449 cells: [0; 3],450 filled: [0; 3],451 };452 for flat in 0..types.size() {453 let region = map.bytes()[flat] as usize;454 out.cells[region] += 1;455 if types.at(flat) != 0 {456 out.filled[region] += 1;457 }458 }459 out460}461462// RADIAL463464/// The tallies of one design against a single integer radius about a centre.465#[derive(Clone, Copy, Debug, PartialEq, Eq)]466pub struct RadialCounts {467 /// The filled cells whose own centre lies within the radius.468 pub seen: u64,469 /// The filled cells lying wholly within the radius.470 pub inside: u64,471 /// The filled cells the sphere of the radius crosses.472 pub cut: u64,473}474475fn root_up(square: u64) -> u64 {476 if square == 0 {477 return 0;478 }479 let mut root = (square as f64).sqrt() as u64;480 while root.saturating_mul(root) < square {481 root += 1;482 }483 while root > 0 && (root - 1) * (root - 1) >= square {484 root -= 1;485 }486 root487}488489fn shell(square: u64) -> u64 {490 root_up(square).div_ceil(2)491}492493/// Counts a design's filled cells against every integer radius about one centre, in exact integer arithmetic.494///495/// Cell centres sit at the index plus one half, so every length is doubled and the centre is given in those doubled units: the lattice corner is zero on each axis and the grid centre is the side. Entry r of the result holds the filled cells whose centre lies within radius r, the filled cells whose whole cell lies within it, and the filled cells the sphere of radius r crosses; squared distances are compared as integers and never as floats.496///497/// The tallies are prefix sums over one pass of the grid: each cell lands in the smallest radius that sees it, the smallest that swallows it and the smallest that clears it, and the columns are summed afterwards.498pub fn radial_census(types: &Tensor, centre: &[i64], r_max: u64) -> Vec<RadialCounts> {499 let rank = types.shape.len();500 let width = r_max as usize + 1;501 let mut seen = vec![0u64; width];502 let mut inside = vec![0u64; width];503 let mut touch = vec![0u64; width];504 let mut index = vec![0i64; rank];505 for flat in 0..types.size() {506 if types.at(flat) != 0 {507 let (mut near, mut mid, mut far) = (0u64, 0u64, 0u64);508 for (axis, coordinate) in index.iter().enumerate() {509 let gap = (2 * coordinate + 1 - centre[axis]).unsigned_abs();510 mid += gap * gap;511 far += (gap + 1) * (gap + 1);512 near += gap.saturating_sub(1) * gap.saturating_sub(1);513 }514 for (column, square) in [(&mut seen, mid), (&mut inside, far), (&mut touch, near)] {515 let at = shell(square);516 if at < width as u64 {517 column[at as usize] += 1;518 }519 }520 }521 for axis in (0..rank).rev() {522 index[axis] += 1;523 if (index[axis] as usize) < types.shape[axis] {524 break;525 }526 index[axis] = 0;527 }528 }529 for column in [&mut seen, &mut inside, &mut touch] {530 for r in 1..width {531 column[r] += column[r - 1];532 }533 }534 (0..width)535 .map(|r| RadialCounts {536 seen: seen[r],537 inside: inside[r],538 cut: touch[r] - inside[r],539 })540 .collect()541}542543// CROSSING SHELL544545/// One box of a crossing shell: where it sits, the seat it takes in its parent and whether the design keeps its path.546#[derive(Clone, Copy, Debug, PartialEq, Eq)]547pub struct ShellBox {548 /// The box's first coordinate in its own level's grid.549 pub x: u64,550 /// The box's second coordinate in its own level's grid.551 pub y: u64,552 /// The seat the box takes in its parent, row-major over the side, and the side squared at the root.553 pub seat: usize,554 /// The parent's place in the level above, and `usize::MAX` at the root or when the level above holds no such box.555 pub parent: usize,556 /// Whether every seat from the root down to this box is one the design keeps.557 pub live: bool,558}559560/// The rooted tree of the boxes one circle crosses, level by level.561#[derive(Clone, Debug)]562pub struct Shell {563 /// The radius in cells.564 pub radius: u64,565 /// The side of a box measured in the boxes one level below it.566 pub number: u64,567 /// The boxes of level `j`, the crossed cells at `0` and the single root last, each level in the arc's own order.568 pub levels: Vec<Vec<ShellBox>>,569 /// The boxes whose parent is missing from the level above, which the crossing identity forbids.570 pub orphans: usize,571}572573/// Lists the level-`level` boxes the circle of radius `radius` crosses, in the arc's own order.574///575/// A box `X` of side `number^level` is crossed when its near corner lies in the closed disc and its far corner outside it, `|X| <= radius / number^level < |X + 1|`, so the level is the whole grid's crossing shell read at that real radius. The columns telescope, the low row of one column being the high row of the next, so the walk is one pass over `floor(radius / number^level) + 1` columns and every comparison is on integers. The order runs along the arc, the first coordinate rising while the second falls, which is why a parent's children are contiguous among its own.576pub fn crossing_shell(radius: u64, number: u64, level: u32) -> Vec<(u64, u64)> {577 let scale = match number.checked_pow(level) {578 Some(scale) if scale <= radius => scale,579 _ => return vec![(0, 0)],580 };581 let square = radius * radius;582 let top = radius / scale;583 let high: Vec<u64> = (0..=top + 1)584 .map(|i| {585 let reach = scale * i;586 if reach > radius {587 0588 } else {589 (square - reach * reach).isqrt() / scale590 }591 })592 .collect();593 let mut out = Vec::with_capacity(2 * top as usize + 1);594 for i in 0..=top as usize {595 for y in (high[i + 1]..=high[i]).rev() {596 out.push((i as u64, y));597 }598 }599 out600}601602/// Builds the whole crossing tree of one radius, pruned by the seats the design keeps.603///604/// The depth is the least level whose grid holds the circle inside one box, so the tree is rooted there and its leaves are the `2 * radius + 1` crossed cells. `keep` reads the design's level-one tile row-major, one flag per seat, and a box is live when every seat from the root down to it is kept, so the live leaves are exactly the crossed cells the design fills.605pub fn crossing_tree(radius: u64, number: u64, keep: &[bool]) -> Shell {606 let mut depth = 0u32;607 while number608 .checked_pow(depth)609 .is_some_and(|scale| scale <= radius)610 {611 depth += 1;612 }613 let seats = (number * number) as usize;614 let mut levels: Vec<Vec<ShellBox>> = (0..=depth)615 .map(|level| {616 crossing_shell(radius, number, level)617 .into_iter()618 .map(|(x, y)| ShellBox {619 x,620 y,621 seat: seats,622 parent: usize::MAX,623 live: true,624 })625 .collect()626 })627 .collect();628 let mut orphans = 0;629 for level in (0..depth as usize).rev() {630 let mut start = 0usize;631 for k in 0..levels[level].len() {632 let here = levels[level][k];633 let (px, py) = (here.x / number, here.y / number);634 let mut at = start;635 while at < levels[level + 1].len()636 && (levels[level + 1][at].x, levels[level + 1][at].y) != (px, py)637 {638 at += 1;639 }640 if at == levels[level + 1].len() {641 orphans += 1;642 levels[level][k].live = false;643 continue;644 }645 start = at;646 let seat = ((here.x % number) * number + here.y % number) as usize;647 let live = levels[level + 1][at].live && keep.get(seat).copied().unwrap_or(false);648 let cell = &mut levels[level][k];649 cell.seat = seat;650 cell.parent = at;651 cell.live = live;652 }653 }654 Shell {655 radius,656 number,657 levels,658 orphans,659 }660}661662#[cfg(test)]663mod tests {664 use super::*;665 use crate::bang::factory::create;666 use mrlycore::atoms;667668 struct Lcg(u64);669670 impl Lcg {671 fn next(&mut self) -> u64 {672 self.0 = self673 .0674 .wrapping_mul(6364136223846793005)675 .wrapping_add(1442695040888963407);676 self.0 >> 33677 }678 fn pick(&mut self, lo: i64, hi: i64) -> i64 {679 lo + (self.next() % (hi - lo + 1) as u64) as i64680 }681 }682683 fn corners(index: &[usize]) -> Vec<Vec<i128>> {684 let d = index.len();685 (0..1usize << d)686 .map(|bits| {687 index688 .iter()689 .enumerate()690 .map(|(k, &i)| (i + ((bits >> k) & 1)) as i128)691 .collect()692 })693 .collect()694 }695696 fn oracle_half(wall: &Half, side: usize, index: &[usize]) -> Region {697 let bound = wall.offset.num as i128 * side as i128;698 let den = wall.offset.den as i128;699 let mut ins = 0;700 let all = corners(index);701 for corner in &all {702 let value: i128 = wall703 .normal704 .iter()705 .zip(corner)706 .map(|(&n, &c)| n as i128 * c)707 .sum();708 if value * den <= bound {709 ins += 1;710 }711 }712 if ins == all.len() {713 Region::In714 } else if ins == 0 {715 Region::Out716 } else {717 Region::Cut718 }719 }720721 fn oracle_ball(center: &[Frac], radius: Frac, side: usize, index: &[usize]) -> Region {722 let l = ball_scale(center, radius) as i128;723 let s = side as i128;724 let r = radius.num as i128 * (l / radius.den as i128) * s;725 if r < 0 {726 return Region::Out;727 }728 let scaled: Vec<i128> = center729 .iter()730 .map(|c| c.num as i128 * (l / c.den as i128) * s)731 .collect();732 let far = corners(index)733 .iter()734 .map(|corner| {735 corner736 .iter()737 .zip(&scaled)738 .map(|(&p, &c)| (p * l - c) * (p * l - c))739 .sum::<i128>()740 })741 .max()742 .unwrap();743 let near: i128 = scaled744 .iter()745 .enumerate()746 .map(|(axis, &c)| {747 let lo = index[axis] as i128 * l;748 let d = c.clamp(lo, lo + l) - c;749 d * d750 })751 .sum();752 if far <= r * r {753 Region::In754 } else if near > r * r {755 Region::Out756 } else {757 Region::Cut758 }759 }760761 fn oracle(shape: &Shape, side: usize, index: &[usize]) -> Region {762 match shape {763 Shape::Ball { center, radius } => oracle_ball(center, *radius, side, index),764 Shape::Polytope { walls } => {765 let mut cut = false;766 for wall in walls {767 match oracle_half(wall, side, index) {768 Region::Out => return Region::Out,769 Region::Cut => cut = true,770 Region::In => {}771 }772 }773 if cut {774 Region::Cut775 } else {776 Region::In777 }778 }779 Shape::Anti(inner) => oracle(inner, side, index).flip(),780 }781 }782783 fn random_shape(rng: &mut Lcg, dimension: usize) -> Shape {784 let core = if rng.pick(0, 1) == 0 {785 Shape::Ball {786 center: (0..dimension)787 .map(|_| Frac::new(rng.pick(-4, 8), rng.pick(1, 4)))788 .collect(),789 radius: Frac::new(rng.pick(0, 6), rng.pick(1, 3)),790 }791 } else {792 Shape::Polytope {793 walls: (0..rng.pick(1, 4))794 .map(|_| Half {795 normal: (0..dimension).map(|_| rng.pick(-3, 3)).collect(),796 offset: Frac::new(rng.pick(-6, 6), rng.pick(1, 4)),797 })798 .collect(),799 }800 };801 if rng.pick(0, 3) == 0 {802 Shape::Anti(Box::new(core))803 } else {804 core805 }806 }807808 fn every_index(dims: &[usize]) -> Vec<Vec<usize>> {809 let size: usize = dims.iter().product();810 (0..size)811 .map(|flat| {812 let mut rem = flat;813 let mut index = vec![0usize; dims.len()];814 for axis in (0..dims.len()).rev() {815 index[axis] = rem % dims[axis];816 rem /= dims[axis];817 }818 index819 })820 .collect()821 }822823 #[test]824 fn classify_matches_the_corner_oracle() {825 let mut rng = Lcg(9);826 for dimension in 1..=3 {827 for _ in 0..40 {828 let shape = random_shape(&mut rng, dimension);829 let side = rng.pick(1, 8) as usize;830 for index in every_index(&vec![side; dimension]) {831 assert_eq!(832 classify(&shape, side, &index),833 oracle(&shape, side, &index),834 "{shape:?} side {side} at {index:?}"835 );836 }837 }838 }839 }840841 #[test]842 fn census_and_crops_partition_the_grid() {843 let types = create(7, 3, 2, 2, 2).unwrap();844 let ball = named("ball", 2, Frac::new(1, 2)).unwrap();845 let tally = census(&ball, &types);846 assert_eq!(tally.cells.iter().sum::<usize>(), 81);847 assert_eq!(tally.filled.iter().sum::<usize>(), 64);848 let kept = crop(&types, &ball, true);849 let anti = crop(&types, &Shape::Anti(Box::new(ball.clone())), false);850 let filled = |t: &Tensor| t.bytes().iter().filter(|&&v| v != 0).count();851 assert_eq!(filled(&kept) + filled(&anti), filled(&types));852 assert_eq!(kept.get(&[0, 0]), 0);853 assert!(filled(&kept) > 0);854 }855856 #[test]857 fn diamond_crop_matches_the_closed_form() {858 let half = Frac::new(1, 2);859 for m in [1usize, 2, 3, 4, 6] {860 let ones = Tensor::full(vec![2 * m, 2 * m], 1);861 let d = named("diamond", 2, half).unwrap();862 let inside = crop(&ones, &d, false);863 assert_eq!(inside.sum() as usize, 2 * m * (m - 1), "side {}", 2 * m);864 }865 }866867 #[test]868 fn sponge_ball_crop_counts() {869 let sponge = create(23, 3, 3, 2, 1).unwrap();870 assert_eq!(sponge.sum(), 20);871 let ball = named("ball", 3, Frac::new(1, 2)).unwrap();872 assert_eq!(crop(&sponge, &ball, true).sum(), 20);873 assert_eq!(crop(&sponge, &ball, false).sum(), 0);874 let tally = census(&ball, &sponge);875 assert_eq!(tally.cells, [0, 26, 1]);876 }877878 #[test]879 fn carpet_ball_crop_trims_the_corners() {880 let carpet = create(7, 3, 2, 2, 2).unwrap();881 let ball = named("ball", 2, Frac::new(1, 2)).unwrap();882 let kept = crop(&carpet, &ball, true);883 assert_eq!(kept.get(&[0, 0]), 0);884 assert_eq!(kept.get(&[8, 8]), 0);885 assert!(kept.sum() > 0);886 assert!(kept.sum() < carpet.sum());887 }888889 #[test]890 fn refine_of_ones_is_a_crop_at_the_finer_side() {891 let ball = named("ball", 2, Frac::new(1, 2)).unwrap();892 for keep_cut in [false, true] {893 let refined = refine(&atoms::ones_2d(2), &ball, 2, 2, keep_cut).unwrap();894 let cropped = crop(&atoms::ones_2d(8), &ball, keep_cut);895 assert_eq!(refined, cropped);896 }897 }898899 #[test]900 fn refine_replicates_the_parent_under_a_covering_shape() {901 let carpet = atoms::carpet_2d(3);902 let everything = named("ball", 2, Frac::whole(2)).unwrap();903 let refined = refine(&carpet, &everything, 3, 1, false).unwrap();904 assert_eq!(refined, carpet.kron(&atoms::ones_2d(3)));905 }906907 #[test]908 fn refine_guards_the_cell_ceiling() {909 let ball = named("ball", 2, Frac::new(1, 2)).unwrap();910 assert!(refine(&atoms::ones_2d(10), &ball, 10, 4, false).is_err());911 assert!(refine(&atoms::ones_3d(3), &ball, 3, 9, false).is_err());912 assert!(refine(&atoms::ones_2d(3), &ball, 0, 1, false).is_err());913 }914915 #[test]916 fn named_catalog_covers_its_dimensions() {917 assert_eq!(918 shapes(2),919 vec!["ball", "box", "diamond", "triangle", "octagon"]920 );921 assert_eq!(922 shapes(3),923 vec![924 "ball",925 "box",926 "diamond",927 "octahedron",928 "tetrahedron",929 "pyramid"930 ]931 );932 let r = Frac::new(1, 2);933 for dimension in 2..=3 {934 for name in shapes(dimension) {935 assert!(named(name, dimension, r).is_ok(), "{name} d{dimension}");936 }937 }938 assert!(named("hexagon", 2, r).is_err());939 assert!(named("triangle", 3, r).is_err());940 assert!(named("blob", 2, r).is_err());941 }942943 #[test]944 fn named_solids_fill_a_sane_share() {945 let r = Frac::new(1, 2);946 for (dimension, side) in [(2usize, 12usize), (3, 8)] {947 let ones = Tensor::full(vec![side; dimension], 1);948 for name in shapes(dimension) {949 let shape = named(name, dimension, r).unwrap();950 let inside = crop(&ones, &shape, false).sum();951 let touched = crop(&ones, &shape, true).sum();952 assert!(inside > 0, "{name} d{dimension} inside");953 assert!(touched >= inside, "{name} d{dimension} touched");954 assert!(touched <= ones.sum(), "{name} d{dimension} bounded");955 }956 }957 }958959 #[test]960 fn frac_reduces_and_computes() {961 assert_eq!(Frac::new(2, 4), Frac::new(1, 2));962 assert_eq!(Frac::new(1, -2), Frac::new(-1, 2));963 assert_eq!(Frac::new(1, 2) + Frac::new(1, 3), Frac::new(5, 6));964 assert_eq!(Frac::new(1, 2) - Frac::new(1, 2), Frac::whole(0));965 assert_eq!(Frac::new(2, 3) * Frac::new(3, 4), Frac::new(1, 2));966 }967968 fn corner_seen(code: u128, dimension: usize, level: usize, r_max: u64) -> Vec<u64> {969 let types = create(code, 3, dimension, 2, level).unwrap();970 let table = radial_census(&types, &vec![0; dimension], r_max);971 (1..=r_max as usize).map(|r| table[r].seen).collect()972 }973974 #[test]975 fn radial_census_counts_the_corner_disc() {976 assert_eq!(977 corner_seen(7, 2, 6, 12),978 [1, 3, 7, 12, 16, 22, 30, 38, 48, 63, 77, 91]979 );980 assert_eq!(corner_seen(23, 3, 4, 8), [1, 4, 13, 28, 47, 65, 95, 137]);981 }982983 #[test]984 fn radial_census_does_not_depend_on_the_level() {985 let shallow = create(7, 3, 2, 2, 4).unwrap();986 let deep = create(7, 3, 2, 2, 5).unwrap();987 let reach = 80;988 let near = radial_census(&shallow, &[0, 0], reach);989 let far = radial_census(&deep, &[0, 0], reach);990 assert_eq!(near, far);991 assert!(near[reach as usize].seen > 0);992 }993994 #[test]995 fn radial_census_matches_the_ball_census() {996 for (code, dimension, level) in [(7u128, 2usize, 4usize), (23, 3, 3)] {997 let types = create(code, 3, dimension, 2, level).unwrap();998 let side = types.shape[0];999 let table = radial_census(&types, &vec![0; dimension], side as u64 - 1);1000 for r in [1usize, 5, 17, side - 1] {1001 let ball = Shape::Ball {1002 center: vec![Frac::whole(0); dimension],1003 radius: Frac::new(r as i64, side as i64),1004 };1005 let tally = census(&ball, &types);1006 assert_eq!(tally.filled[2] as u64, table[r].inside, "in r={r}");1007 assert_eq!(tally.filled[1] as u64, table[r].cut, "cut r={r}");1008 assert!(table[r].inside <= table[r].seen);1009 assert!(table[r].seen <= table[r].inside + table[r].cut);1010 }1011 }1012 }10131014 fn brute_shell(radius: u64, scale: u64) -> Vec<(u64, u64)> {1015 let mut out = Vec::new();1016 let reach = radius / scale + 1;1017 for x in 0..=reach {1018 for y in 0..=reach {1019 let near = (scale * x).pow(2) + (scale * y).pow(2);1020 let far = (scale * (x + 1)).pow(2) + (scale * (y + 1)).pow(2);1021 if near <= radius * radius && radius * radius < far {1022 out.push((x, y));1023 }1024 }1025 }1026 out.sort_by_key(|&(x, y)| (x, std::cmp::Reverse(y)));1027 out1028 }10291030 #[test]1031 fn crossing_shell_matches_the_brute_sweep() {1032 for radius in 1..=90u64 {1033 for level in 0..5u32 {1034 let scale = 3u64.pow(level);1035 assert_eq!(1036 crossing_shell(radius, 3, level),1037 brute_shell(radius, scale),1038 "r={radius} j={level}"1039 );1040 }1041 }1042 }10431044 #[test]1045 fn crossing_shell_is_two_floor_plus_one() {1046 for radius in 1..=400u64 {1047 for level in 0..8u32 {1048 let want = 2 * (radius / 3u64.pow(level)) + 1;1049 assert_eq!(crossing_shell(radius, 3, level).len() as u64, want);1050 let five = 2 * (radius / 5u64.pow(level.min(5))) + 1;1051 assert_eq!(crossing_shell(radius, 5, level.min(5)).len() as u64, five);1052 }1053 }1054 }10551056 #[test]1057 fn crossing_tree_hangs_every_box_on_a_crossed_parent() {1058 let tile = create(7, 3, 2, 2, 1).unwrap();1059 let keep: Vec<bool> = tile.bytes().iter().map(|&b| b != 0).collect();1060 for radius in 1..=120u64 {1061 let tree = crossing_tree(radius, 3, &keep);1062 assert_eq!(tree.orphans, 0, "r={radius}");1063 assert_eq!(tree.levels.last().unwrap().len(), 1);1064 assert_eq!(tree.levels[0].len() as u64, 2 * radius + 1);1065 for level in 0..tree.levels.len() - 1 {1066 for cell in &tree.levels[level] {1067 let parent = tree.levels[level + 1][cell.parent];1068 assert_eq!((parent.x, parent.y), (cell.x / 3, cell.y / 3));1069 assert_eq!(cell.seat, ((cell.x % 3) * 3 + cell.y % 3) as usize);1070 assert!(!cell.live || parent.live);1071 }1072 }1073 }1074 }10751076 #[test]1077 fn the_live_leaves_are_the_designs_crossed_cells() {1078 for code in [7u128, 5, 11, 15] {1079 let tile = create(code, 3, 2, 2, 1).unwrap();1080 let keep: Vec<bool> = tile.bytes().iter().map(|&b| b != 0).collect();1081 for depth in 1..=4usize {1082 let grid = create(code, 3, 2, 2, depth).unwrap();1083 let side = 3u64.pow(depth as u32);1084 let table = radial_census(&grid, &[0, 0], side - 1);1085 for radius in side / 3..side {1086 let tree = crossing_tree(radius.max(1), 3, &keep);1087 assert_eq!(tree.levels.len(), depth + 1, "code={code} r={radius}");1088 let live = tree.levels[0].iter().filter(|cell| cell.live).count() as u64;1089 assert_eq!(live, table[radius as usize].cut, "code={code} r={radius}");1090 }1091 }1092 }1093 }1094}