apollonian.rs
14.7 kB · rust · 446 lines
1use crate::factor::gcd;2use crate::lattice::farey;3use mrlycore::errors::{value_error, Result};45/// The largest curvature a packing is grown to.6pub const CURVATURE_CAP: i64 = 8192;7/// The most circles one growth makes before it gives up.8pub const CIRCLE_CAP: usize = 200_000;9/// The deepest the Farey stack is read against a packing.10pub const ORDER_CAP: usize = 64;11/// The root quadruples on offer: the strip first, then the bounded packings named by their curvatures.12pub const ROOTS: [&str; 4] = ["strip", "-1,2,2,3", "-2,3,6,7", "-3,4,12,13"];1314// THE OBJECTS1516/// A circle in the integer coordinates `(k, k x, k y)`: a line is `k = 0` with `(k x, k y)` its outward unit normal, and the curvature is negative on the circle that contains a bounded packing.17#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]18pub struct Circle {19 /// The curvature.20 pub k: i64,21 /// The curvature times the centre's abscissa.22 pub x: i64,23 /// The curvature times the centre's ordinate.24 pub y: i64,25}2627impl Circle {28 /// Whether the circle is a line.29 pub fn is_line(self) -> bool {30 self.k == 031 }3233 /// The radius, none on a line.34 pub fn radius(self) -> Option<f64> {35 (self.k != 0).then(|| 1.0 / self.k.abs() as f64)36 }3738 /// The centre, none on a line.39 pub fn centre(self) -> Option<(f64, f64)> {40 (self.k != 0).then(|| (self.x as f64 / self.k as f64, self.y as f64 / self.k as f64))41 }42}4344/// Four mutually tangent circles.45pub type Quad = [Circle; 4];4647fn column(q: &Quad, at: usize) -> [i64; 4] {48 let read = |c: &Circle| match at {49 0 => c.k,50 1 => c.x,51 _ => c.y,52 };53 [read(&q[0]), read(&q[1]), read(&q[2]), read(&q[3])]54}5556/// The bilinear form `B(u, v) = (sum u)(sum v) - 2 sum u v` that the reflection preserves.57///58/// ```59/// use mrlynum::apollonian::form;60/// assert_eq!(form([-1, 2, 2, 3], [-1, 2, 2, 3]), 0);61/// assert_eq!(form([0, 1, -1, 0], [0, 1, -1, 0]), -4);62/// ```63pub fn form(u: [i64; 4], v: [i64; 4]) -> i128 {64 let su: i128 = u.iter().map(|&a| i128::from(a)).sum();65 let sv: i128 = v.iter().map(|&a| i128::from(a)).sum();66 let dot: i128 = (0..4).map(|i| i128::from(u[i]) * i128::from(v[i])).sum();67 su * sv - 2 * dot68}6970/// Whether the quadruple carries all six exact invariants: Descartes `B(k, k) = 0`, the position half `B(k, kx) = B(k, ky) = B(kx, ky) = 0`, and the frame `B(kx, kx) = B(ky, ky) = -4`.71pub fn sound(q: &Quad) -> bool {72 let (k, x, y) = (column(q, 0), column(q, 1), column(q, 2));73 form(k, k) == 074 && form(k, x) == 075 && form(k, y) == 076 && form(x, y) == 077 && form(x, x) == -478 && form(y, y) == -479}8081/// Reflects the circle at the seat through the other three, `v' = 2(v_1 + v_2 + v_3) - v` on all three coordinates at once, which is the second root of the Descartes quadratic and needs no square root.82pub fn reflect(q: &Quad, at: usize) -> Circle {83 let (mut k, mut x, mut y) = (0i64, 0i64, 0i64);84 for (i, c) in q.iter().enumerate() {85 if i != at {86 k += c.k;87 x += c.x;88 y += c.y;89 }90 }91 Circle {92 k: 2 * k - q[at].k,93 x: 2 * x - q[at].x,94 y: 2 * y - q[at].y,95 }96}9798/// The quadruple with the circle at the seat replaced by its reflection.99pub fn swap(q: &Quad, at: usize) -> Quad {100 let mut out = *q;101 out[at] = reflect(q, at);102 out103}104105fn circle(k: i64, x: i64, y: i64) -> Circle {106 Circle { k, x, y }107}108109/// The named root quadruple: `strip` is the two lines a unit apart holding the circles at `0` and `1`, and the rest are bounded packings named by their four curvatures.110///111/// ```112/// use mrlynum::apollonian::{root, sound};113/// for name in mrlynum::apollonian::ROOTS {114/// assert!(sound(&root(name).unwrap()));115/// }116/// ```117pub fn root(name: &str) -> Result<Quad> {118 Ok(match name {119 "strip" => [120 circle(0, 0, -1),121 circle(0, 0, 1),122 circle(2, 0, 1),123 circle(2, 2, 1),124 ],125 "-1,2,2,3" => [126 circle(-1, 0, 0),127 circle(2, 1, 0),128 circle(2, -1, 0),129 circle(3, 0, 2),130 ],131 "-2,3,6,7" => [132 circle(-2, -1, 0),133 circle(3, 1, 0),134 circle(6, 5, 0),135 circle(7, 5, 2),136 ],137 "-3,4,12,13" => [138 circle(-3, -1, 0),139 circle(4, 1, 0),140 circle(12, 7, 0),141 circle(13, 7, 2),142 ],143 _ => return value_error(format!("the root must be one of {}.", ROOTS.join(", "))),144 })145}146147// THE GROWTH148149/// A packing grown from its root in exact integers.150pub struct Packing {151 /// The root quadruple the growth started from.152 pub root: Quad,153 /// Whether the root carries a line, so the packing is the strip and the growth keeps one period.154 pub strip: bool,155 /// The curvature the growth stopped at.156 pub cap: i64,157 /// The circles the growth made, the root excluded, in curvature order.158 pub circles: Vec<Circle>,159 /// The quadruples the growth made, the root counted.160 pub quads: u64,161 /// The quadruples that failed one of the six invariants.162 pub broken: u64,163 /// The circles centred outside the open period, which the strip must have none of.164 pub strayed: u64,165}166167/// Grows the named packing to the curvature cap, one circle per node of the reflection tree and the root quadruple excluded, so `circles.len()` is the census `N(T)`. On the strip only the two root swaps that replace a line are taken, which are exactly the two that stay inside one period.168pub fn grow(name: &str, cap: i64) -> Result<Packing> {169 let seed = root(name)?;170 if !(2..=CURVATURE_CAP).contains(&cap) {171 return value_error(format!(172 "the curvature cap must be between 2 and {CURVATURE_CAP}."173 ));174 }175 let strip = seed.iter().any(|c| c.is_line());176 let mut out = Packing {177 root: seed,178 strip,179 cap,180 circles: Vec::new(),181 quads: 1,182 broken: u64::from(!sound(&seed)),183 strayed: 0,184 };185 let mut stack: Vec<(Quad, usize)> = Vec::new();186 for i in 0..4 {187 if strip && !seed[i].is_line() {188 continue;189 }190 if reflect(&seed, i).k <= cap {191 stack.push((swap(&seed, i), i));192 }193 }194 while let Some((q, last)) = stack.pop() {195 out.quads += 1;196 if !sound(&q) {197 out.broken += 1;198 }199 let made = q[last];200 if out.circles.len() >= CIRCLE_CAP {201 return value_error(format!(202 "the packing passes {CIRCLE_CAP} circles below that curvature; lower the cap."203 ));204 }205 out.circles.push(made);206 if strip && (made.x <= 0 || made.x >= made.k) {207 out.strayed += 1;208 }209 for j in 0..4 {210 if j == last {211 continue;212 }213 let next = reflect(&q, j);214 if next.k > cap || next.k <= q[j].k {215 continue;216 }217 stack.push((swap(&q, j), j));218 }219 }220 out.circles.sort_unstable();221 Ok(out)222}223224/// The box the packing is drawn in: one period of the strip, or the box of the circle that contains a bounded packing.225pub fn frame(p: &Packing) -> [f64; 4] {226 if p.strip {227 return [0.0, 0.0, 1.0, 1.0];228 }229 let outer = p230 .root231 .iter()232 .find(|c| c.k < 0)233 .copied()234 .unwrap_or(p.root[0]);235 let (x, y) = outer.centre().unwrap_or((0.0, 0.0));236 let r = outer.radius().unwrap_or(1.0);237 [x - r, y - r, x + r, y + r]238}239240// THE FORD CIRCLES241242/// A tangency point on the line `y = 0`: the reduced fraction the circle rests at and the curvature it carries.243#[derive(Clone, Copy, Debug, PartialEq, Eq)]244pub struct Touch {245 /// The numerator of the reduced fraction.246 pub num: i64,247 /// The denominator of the reduced fraction.248 pub den: i64,249 /// The curvature of the circle resting there, which the Ford identification forces to be `2 den^2`.250 pub k: i64,251}252253fn reduce(c: Circle) -> (i64, i64) {254 let g = gcd(c.x.unsigned_abs() as usize, c.k.unsigned_abs() as usize) as i64;255 (c.x / g, c.k / g)256}257258/// Whether the circle has positive curvature and is tangent to the line `y = 0`, which in these coordinates reads `k > 0` and `k y = 1`: the curvature guard is what excludes the line `y = 1`, which is `(0, 0, 1)`.259pub fn on_line(c: Circle) -> bool {260 c.k > 0 && c.y == 1261}262263/// Whether the circle is the Ford circle over its own tangency point: curvature `2 b^2` and abscissa `2 a b` at the reduced `a/b`.264pub fn is_ford(c: Circle) -> bool {265 if !on_line(c) {266 return false;267 }268 let (a, b) = reduce(c);269 c.k == 2 * b * b && c.x == 2 * a * b270}271272/// The tangency points on the line `y = 0`, ascending: one per circle of the packing with `k y = 1`, the root excluded. Empty off the strip.273pub fn touches(p: &Packing) -> Vec<Touch> {274 let mut out: Vec<Touch> = p275 .circles276 .iter()277 .filter(|c| on_line(**c))278 .map(|&c| {279 let (num, den) = reduce(c);280 Touch { num, den, k: c.k }281 })282 .collect();283 out.sort_unstable_by(|a, b| {284 (a.num as i128 * b.den as i128).cmp(&(b.num as i128 * a.den as i128))285 });286 out287}288289// THE SHADOW290291/// The Farey stack read against the packing's tangency points.292pub struct Shadow {293 /// The depth the stack is read at.294 pub order: usize,295 /// The curvature a circle of denominator the order carries, twice the order squared.296 pub reach: i64,297 /// Whether the packing was grown far enough to carry every node of that depth.298 pub covered: bool,299 /// The stack's nodes inside the open period.300 pub nodes: usize,301 /// The packing's tangency points of denominator at most the order.302 pub touched: usize,303 /// The nodes no tangency point rests on, plus the tangency points no node lights.304 pub missed: usize,305 /// The circles below the reach tangent to the line that are not Ford circles.306 pub offford: usize,307 /// The brightness of the period summed node by node, the node `0/1` on the period's edge counted.308 pub bright: u128,309 /// The closed form that brightness lands on, `Q(Q + 1)/2`.310 pub want: u128,311}312313/// Reads the Farey stack of the order against the packing: the nodes lit inside the open period against the tangency points of the line-tangent circles of curvature at most `2 Q^2`, and the brightness `floor(Q/b)` summed on the nodes against `Q(Q + 1)/2`. Off the strip there is no line and every count is zero.314pub fn shadow(p: &Packing, order: usize) -> Result<Shadow> {315 if order == 0 || order > ORDER_CAP {316 return value_error(format!("the depth must be between 1 and {ORDER_CAP}."));317 }318 let reach = 2 * (order as i64) * (order as i64);319 let mut out = Shadow {320 order,321 reach,322 covered: p.strip && p.cap >= reach,323 nodes: 0,324 touched: 0,325 missed: 0,326 offford: 0,327 bright: 0,328 want: (order as u128) * (order as u128 + 1) / 2,329 };330 if !p.strip {331 return Ok(out);332 }333 out.offford = p334 .circles335 .iter()336 .filter(|c| on_line(**c) && c.k <= reach && !is_ford(**c))337 .count();338 let mut held: Vec<(i64, i64)> = touches(p)339 .into_iter()340 .filter(|t| t.den <= order as i64)341 .map(|t| (t.num, t.den))342 .collect();343 let mut lit: Vec<(i64, i64)> = farey(order)344 .iter()345 .filter(|n| n.num > 0 && n.num < n.den)346 .map(|n| (n.num as i64, n.den as i64))347 .collect();348 out.touched = held.len();349 out.nodes = lit.len();350 held.sort_unstable();351 lit.sort_unstable();352 out.missed = held353 .iter()354 .filter(|pair| lit.binary_search(pair).is_err())355 .count()356 + lit357 .iter()358 .filter(|pair| held.binary_search(pair).is_err())359 .count();360 out.bright = order as u128361 + lit362 .iter()363 .map(|&(_, den)| (order as u128) / den as u128)364 .sum::<u128>();365 Ok(out)366}367368#[cfg(test)]369mod tests {370 use super::*;371372 #[test]373 fn the_reflection_keeps_the_six_invariants_on_every_quadruple() {374 for name in ROOTS {375 let p = grow(name, 512).unwrap();376 assert_eq!(p.broken, 0);377 assert_eq!(p.strayed, 0);378 assert!(p.quads > p.circles.len() as u64);379 }380 let seed = root("strip").unwrap();381 assert_eq!(reflect(&seed, 0), Circle { k: 8, x: 4, y: 7 });382 assert_eq!(reflect(&seed, 1), Circle { k: 8, x: 4, y: 1 });383 assert_eq!(reflect(&seed, 2), Circle { k: 2, x: 4, y: 1 });384 }385386 #[test]387 fn the_census_lands_on_the_counts_the_generator_prints() {388 assert_eq!(grow("strip", 1000).unwrap().circles.len(), 950);389 assert_eq!(grow("-1,2,2,3", 1000).unwrap().circles.len(), 3325);390 assert_eq!(grow("strip", 2048).unwrap().circles.len(), 2448);391 }392393 #[test]394 fn every_line_tangent_circle_is_the_ford_circle_over_its_own_fraction() {395 let p = grow("strip", 2048).unwrap();396 let marks = touches(&p);397 assert_eq!(marks.len(), 323);398 assert!(p399 .circles400 .iter()401 .filter(|c| on_line(**c))402 .all(|&c| is_ford(c)));403 for mark in &marks {404 assert_eq!(mark.k, 2 * mark.den * mark.den);405 assert_eq!(gcd(mark.num as usize, mark.den as usize), 1);406 }407 assert!(marks408 .windows(2)409 .all(|pair| { pair[0].num * pair[1].den < pair[1].num * pair[0].den }));410 assert!(touches(&grow("-1,2,2,3", 2048).unwrap()).is_empty());411 }412413 #[test]414 fn the_stack_is_the_shadow_of_the_line_tangent_circles() {415 let p = grow("strip", 2048).unwrap();416 let read = shadow(&p, 32).unwrap();417 assert!(read.covered);418 assert_eq!(read.reach, 2048);419 assert_eq!((read.nodes, read.touched, read.missed), (323, 323, 0));420 assert_eq!(read.offford, 0);421 assert_eq!((read.bright, read.want), (528, 528));422 let shallow = shadow(&p, 16).unwrap();423 assert_eq!(424 (425 shallow.nodes,426 shallow.touched,427 shallow.missed,428 shallow.bright429 ),430 (79, 79, 0, 136)431 );432 assert!(!shadow(&grow("strip", 512).unwrap(), 32).unwrap().covered);433 assert!(shadow(&p, ORDER_CAP + 1).is_err());434 }435436 #[test]437 fn a_bounded_root_frames_its_own_outer_circle() {438 assert_eq!(frame(&grow("strip", 512).unwrap()), [0.0, 0.0, 1.0, 1.0]);439 assert_eq!(440 frame(&grow("-1,2,2,3", 512).unwrap()),441 [-1.0, -1.0, 1.0, 1.0]442 );443 assert!(grow("gasket", 512).is_err());444 assert!(grow("strip", CURVATURE_CAP + 1).is_err());445 }446}