main.rs
36.5 kB · rust · 1109 lines
1mod derive;2mod transfer;34use mrlycore::Tensor;5use mrlymath::bang::factory::create;6use mrlymath::shape::{census, named, Frac, Shape};78// ROOTS910fn ceil_sqrt(value: u64) -> u64 {11 if value == 0 {12 return 0;13 }14 let mut root = (value as f64).sqrt() as u64;15 while root.saturating_mul(root) < value {16 root += 1;17 }18 while root > 0 && (root - 1) * (root - 1) >= value {19 root -= 1;20 }21 root22}2324fn bucket(quad: u64) -> u64 {25 ceil_sqrt(quad).div_ceil(2)26}2728// DESIGNS2930fn design(dimension: usize, level: usize) -> Tensor {31 match dimension {32 2 => create(7, 3, 2, 2, level).unwrap(),33 _ => create(23, 3, 3, 2, level).unwrap(),34 }35}3637// SWEEP3839struct Sweep {40 r_max: u64,41 seen: Vec<u64>,42 inside: Vec<u64>,43 touch: Vec<u64>,44 all_seen: Vec<u64>,45 all_inside: Vec<u64>,46 all_touch: Vec<u64>,47}4849fn sweep(types: &Tensor, shift: i64, r_max: u64) -> Sweep {50 let dims = types.shape.clone();51 let rank = dims.len();52 let width = (r_max + 2) as usize;53 let mut out = Sweep {54 r_max,55 seen: vec![0; width],56 inside: vec![0; width],57 touch: vec![0; width],58 all_seen: vec![0; width],59 all_inside: vec![0; width],60 all_touch: vec![0; width],61 };62 let mut index = vec![0i64; rank];63 let bytes = types.bytes();64 for cell in bytes {65 let mut centre = 0u64;66 let mut far = 0u64;67 let mut near = 0u64;68 for coordinate in &index {69 let offset = (2 * coordinate + 1 - shift).unsigned_abs();70 centre += offset * offset;71 far += (offset + 1) * (offset + 1);72 let low = offset.saturating_sub(1);73 near += low * low;74 }75 let slots = [bucket(centre), bucket(far), bucket(near)];76 let filled = *cell != 0;77 for (which, slot) in slots.iter().enumerate() {78 if *slot > r_max {79 continue;80 }81 let at = *slot as usize;82 match which {83 0 => {84 out.all_seen[at] += 1;85 if filled {86 out.seen[at] += 1;87 }88 }89 1 => {90 out.all_inside[at] += 1;91 if filled {92 out.inside[at] += 1;93 }94 }95 _ => {96 out.all_touch[at] += 1;97 if filled {98 out.touch[at] += 1;99 }100 }101 }102 }103 for axis in (0..rank).rev() {104 index[axis] += 1;105 if (index[axis] as usize) < dims[axis] {106 break;107 }108 index[axis] = 0;109 }110 }111 for column in [112 &mut out.seen,113 &mut out.inside,114 &mut out.touch,115 &mut out.all_seen,116 &mut out.all_inside,117 &mut out.all_touch,118 ] {119 for r in 1..column.len() {120 column[r] += column[r - 1];121 }122 }123 out124}125126impl Sweep {127 fn cut(&self, r: u64) -> u64 {128 self.touch[r as usize] - self.inside[r as usize]129 }130 fn all_cut(&self, r: u64) -> u64 {131 self.all_touch[r as usize] - self.all_inside[r as usize]132 }133}134135// ORACLE136137fn ball(dimension: usize, shift: i64, radius: Frac) -> Shape {138 if shift == 0 {139 Shape::Ball {140 center: vec![Frac::whole(0); dimension],141 radius,142 }143 } else {144 named("ball", dimension, radius).unwrap()145 }146}147148fn oracle(label: &str, types: &Tensor, shift: i64, side: usize, table: &Sweep, radii: &[u64]) {149 for r in radii {150 let shape = ball(types.shape.len(), shift, Frac::new(*r as i64, side as i64));151 let tally = census(&shape, types);152 assert_eq!(tally.filled[2] as u64, table.inside[*r as usize]);153 assert_eq!(tally.filled[1] as u64, table.cut(*r));154 assert_eq!(tally.cells[2] as u64, table.all_inside[*r as usize]);155 assert_eq!(tally.cells[1] as u64, table.all_cut(*r));156 println!(157 "circle-crop oracle {label} r={r} census_in={} census_cut={} sweep_in={} sweep_cut={}",158 tally.filled[2],159 tally.filled[1],160 table.inside[*r as usize],161 table.cut(*r)162 );163 }164}165166// TREND167168fn trend(label: &str, name: &str, running: &[f64], from: u64, to: u64) {169 if to < from * 3 {170 return;171 }172 let mut least = f64::INFINITY;173 let mut most = f64::NEG_INFINITY;174 let mut total = 0.0f64;175 let mut count = 0u64;176 for r in from..=(to / 3) {177 let here = running[r as usize];178 let there = running[(r * 3) as usize];179 if here <= 0.0 || there <= 0.0 {180 continue;181 }182 let step = (there / here).ln() / 3.0f64.ln();183 least = least.min(step);184 most = most.max(step);185 total += step;186 count += 1;187 }188 if count == 0 {189 return;190 }191 let mut sum_x = 0.0f64;192 let mut sum_y = 0.0f64;193 let mut sum_xx = 0.0f64;194 let mut sum_xy = 0.0f64;195 let mut points = 0.0f64;196 for r in from..=to {197 if running[r as usize] <= 0.0 {198 continue;199 }200 let x = (r as f64).ln();201 let y = running[r as usize].ln();202 sum_x += x;203 sum_y += y;204 sum_xx += x * x;205 sum_xy += x * y;206 points += 1.0;207 }208 let fitted = (points * sum_xy - sum_x * sum_y) / (points * sum_xx - sum_x * sum_x);209 let ends =210 (running[to as usize] / running[from as usize]).ln() / (to as f64 / from as f64).ln();211 println!(212 "circle-crop trend {label} {name} r={from}..{} steps={count} min={:.6} mean={:.6} max={:.6} fitted={fitted:.6} ends={ends:.6}",213 to,214 least,215 total / count as f64,216 most217 );218}219220// WINDOWS221222fn windows(label: &str, series: &[(&str, Vec<f64>, u64)]) {223 for (name, values, limit) in series {224 let mut lows: Vec<(u64, f64, u64)> = Vec::new();225 let mut power = 1u64;226 while power <= *limit {227 let top = (power * 3 - 1).min(*limit);228 let mut best = 0.0f64;229 let mut spot = power;230 for r in power..=top {231 if values[r as usize] > best {232 best = values[r as usize];233 spot = r;234 }235 }236 lows.push((power, best, spot));237 power *= 3;238 }239 for step in 0..lows.len() {240 let (start, best, spot) = lows[step];241 let raw = if start > 1 && best > 0.0 {242 best.ln() / (start as f64).ln()243 } else {244 f64::NAN245 };246 let slope = if step + 1 < lows.len() && best > 0.0 && lows[step + 1].1 > 0.0 {247 (lows[step + 1].1 / best).ln() / 3.0f64.ln()248 } else {249 f64::NAN250 };251 println!(252 "circle-crop window {label} {name} k={step} r={start}..{} max={best:.3} at={spot} at_over_start={:.4} raw={raw:.6} slope={slope:.6}",253 (start * 3 - 1).min(*limit),254 spot as f64 / start as f64255 );256 }257 }258}259260// MEANS261262fn cap(dimension: usize, r: u64) -> f64 {263 match dimension {264 2 => (3 * r + 5) as f64,265 _ => std::f64::consts::PI * 3.0f64.sqrt() * (r * r + 1) as f64,266 }267}268269fn means(label: &str, table: &Sweep, fill: u64, dimension: usize) {270 let r_max = table.r_max;271 let mass = fill as f64;272 let share = 3.0f64.powi(dimension as i32) / mass;273 let edge = dimension as i32 - 1;274 let mut prior_least = f64::NAN;275 let mut prior_most = f64::NAN;276 let mut power = 1u64;277 while power * 3 <= r_max {278 let start = power;279 let stop = power * 3 - 1;280 let span = (stop - start + 1) as f64;281 let level = power.ilog(3) + 1;282 let weight = share.powi(level as i32);283 let mut total = 0u64;284 let mut least = f64::INFINITY;285 let mut most = f64::NEG_INFINITY;286 let mut full_least = f64::INFINITY;287 let mut full_most = f64::NEG_INFINITY;288 let mut carried = 0.0f64;289 for r in start..=stop {290 assert!(table.cut(r) <= table.all_cut(r));291 total += table.cut(r);292 let phi = table.cut(r) as f64 * weight / table.all_cut(r) as f64;293 least = least.min(phi);294 most = most.max(phi);295 carried += phi;296 let density = table.all_cut(r) as f64 / (r as f64).powi(edge);297 full_least = full_least.min(density);298 full_most = full_most.max(density);299 }300 let exact_low =301 table.inside[(power * 3) as usize] as i64 - table.touch[start as usize] as i64;302 let exact_high =303 2 * (table.touch[(power * 3) as usize] as i64 - table.inside[start as usize] as i64);304 assert!(exact_low <= total as i64);305 assert!(total as i64 <= exact_high);306 let mut depth = 0u32;307 let mut reach = power;308 while reach * 3 <= r_max {309 reach *= 3;310 depth += 1;311 }312 let scale = mass.powi(depth as i32);313 let main_low = table.inside[reach as usize] as f64 / scale;314 let main_high = table.touch[reach as usize] as f64 / scale;315 let slack = (table.cut(power * 3) + table.cut(start)) as f64;316 let form_low = (mass - 1.0) * main_low - slack;317 let form_high = 2.0 * ((mass - 1.0) * main_high + slack);318 assert!(form_low <= total as f64);319 assert!(total as f64 <= form_high);320 let kappa_low = total as f64 / ((mass - 1.0) * main_high);321 let kappa_high = total as f64 / ((mass - 1.0) * main_low);322 println!(323 "circle-crop mean {label} k={} r={start}..{stop} sum={total} mean={:.6} exact_low={exact_low} exact_high={exact_high} form_low={:.6} form_high={:.6} kappa_low={:.6} kappa_high={:.6}",324 level - 1,325 total as f64 / span,326 (form_low * 1e6).floor() / 1e6,327 (form_high * 1e6).ceil() / 1e6,328 (kappa_low * 1e6).floor() / 1e6,329 (kappa_high * 1e6).ceil() / 1e6330 );331 let mean_phi = carried / span;332 let ground = weight * form_low.max(0.0) / span / cap(dimension, stop);333 assert!(mean_phi >= ground);334 println!(335 "circle-crop factor {label} k={} r={start}..{stop} level={level} min={:.6} mean={:.6} max={:.6} ground={:.6} min_step={:.6} max_step={:.6} full_low={:.6} full_high={:.6}",336 level - 1,337 (least * 1e6).floor() / 1e6,338 mean_phi,339 (most * 1e6).ceil() / 1e6,340 (ground * 1e6).floor() / 1e6,341 least - prior_least,342 most - prior_most,343 (full_least * 1e6).floor() / 1e6,344 (full_most * 1e6).ceil() / 1e6345 );346 prior_least = least;347 prior_most = most;348 power *= 3;349 }350}351352// DIGITS353354fn root_floor(value: u64) -> u64 {355 if value == 0 {356 return 0;357 }358 let mut root = (value as f64).sqrt() as u64;359 while root * root > value {360 root -= 1;361 }362 while (root + 1) * (root + 1) <= value {363 root += 1;364 }365 root366}367368fn ones_table(level: usize) -> Vec<u32> {369 let size = 3usize.pow(level as u32);370 let mut out = vec![0u32; size];371 for value in 1..size {372 out[value] = (out[value / 3] << 1) | u32::from(value % 3 == 1);373 }374 out375}376377fn trim(value: f64, up: bool) -> f64 {378 if up {379 (value * 1e6).ceil() / 1e6380 } else {381 (value * 1e6).floor() / 1e6382 }383}384385fn list(values: &[f64]) -> String {386 values387 .iter()388 .map(|value| format!("{value:.6}"))389 .collect::<Vec<String>>()390 .join(",")391}392393struct Shell {394 total: u64,395 filled: u64,396 bad: Vec<u64>,397 pair: Vec<u64>,398}399400fn record(mask: u32, level: usize, out: &mut Shell) {401 out.total += 1;402 if mask == 0 {403 out.filled += 1;404 return;405 }406 let mut rest = mask;407 while rest != 0 {408 let low = rest.trailing_zeros() as usize;409 out.bad[low] += 1;410 let mut other = rest & (rest - 1);411 while other != 0 {412 out.pair[low * level + other.trailing_zeros() as usize] += 1;413 other &= other - 1;414 }415 rest &= rest - 1;416 }417}418419fn shell(ones: &[u32], dimension: usize, radius: u64, level: usize, out: &mut Shell) {420 out.total = 0;421 out.filled = 0;422 for slot in out.bad.iter_mut() {423 *slot = 0;424 }425 for slot in out.pair.iter_mut() {426 *slot = 0;427 }428 let square = radius * radius;429 if dimension == 2 {430 for i in 0..=radius {431 let high = root_floor(square - i * i);432 let step = (i + 1) * (i + 1);433 let low = if step >= square {434 0435 } else {436 root_floor(square - step)437 };438 let first = ones[i as usize];439 for y in low..=high {440 record(first & ones[y as usize], level, out);441 }442 }443 } else {444 for i in 0..=radius {445 let first = ones[i as usize];446 let span = root_floor(square - i * i);447 for j in 0..=span {448 let base = i * i + j * j;449 let high = root_floor(square - base);450 let step = (i + 1) * (i + 1) + (j + 1) * (j + 1);451 let low = if step >= square {452 0453 } else {454 root_floor(square - step)455 };456 let second = ones[j as usize];457 let both = first & second;458 let either = first | second;459 for z in low..=high {460 record(both | (either & ones[z as usize]), level, out);461 }462 }463 }464 }465}466467fn digit_census(label: &str, table: &Sweep, fill: u64, dimension: usize) {468 let r_max = table.r_max;469 let mass = fill as f64;470 let room = 3.0f64.powi(dimension as i32);471 let share = room / mass;472 let null = 1.0 - mass / room;473 let ones = ones_table((r_max + 1).ilog(3) as usize);474 let mut fine_low = f64::INFINITY;475 let mut fine_high = f64::NEG_INFINITY;476 let mut scale_low = f64::INFINITY;477 let mut scale_high = f64::NEG_INFINITY;478 let mut near_low = f64::INFINITY;479 let mut near_high = f64::NEG_INFINITY;480 let mut far_low = f64::INFINITY;481 let mut far_high = f64::NEG_INFINITY;482 let mut load_high = f64::NEG_INFINITY;483 let mut ind_least = f64::INFINITY;484 let mut ind_most = f64::NEG_INFINITY;485 let mut psi_least = f64::INFINITY;486 let mut psi_most = f64::NEG_INFINITY;487 let mut count = 0u64;488 let mut power = 1u64;489 while power * 3 <= r_max {490 let start = power;491 let stop = power * 3 - 1;492 let span = (stop - start + 1) as f64;493 let level = power.ilog(3) as usize + 1;494 let mut state = Shell {495 total: 0,496 filled: 0,497 bad: vec![0; level],498 pair: vec![0; level * level],499 };500 let mut sum_total = 0u64;501 let mut sum_bad = vec![0u64; level];502 let mut sum_pair = vec![0u64; level * level];503 let mut rate_low = vec![f64::INFINITY; level];504 let mut rate_high = vec![f64::NEG_INFINITY; level];505 let mut rate_sum = vec![0.0f64; level];506 let mut ind_low = f64::INFINITY;507 let mut ind_high = f64::NEG_INFINITY;508 let mut ind_sum = 0.0f64;509 let mut psi_low = f64::INFINITY;510 let mut psi_high = f64::NEG_INFINITY;511 let mut psi_sum = 0.0f64;512 let mut load_sum = 0.0f64;513 let mut load_most = f64::NEG_INFINITY;514 for r in start..=stop {515 shell(&ones, dimension, r, level, &mut state);516 assert_eq!(state.total, table.all_cut(r));517 assert_eq!(state.filled, table.cut(r));518 if dimension == 2 {519 assert_eq!(state.total, 2 * r + 1);520 }521 let mut union = 0u64;522 let mut product = 1.0f64;523 for j in 0..level {524 if dimension == 2 {525 let parents = 2 * (r / 3u64.pow(j as u32 + 1)) + 1;526 assert!(state.bad[j] <= 2 * 3u64.pow(j as u32) * parents);527 }528 union += state.bad[j];529 sum_bad[j] += state.bad[j];530 let rate = state.bad[j] as f64 / state.total as f64;531 rate_low[j] = rate_low[j].min(rate);532 rate_high[j] = rate_high[j].max(rate);533 rate_sum[j] += rate;534 product *= 1.0 - rate;535 }536 assert!(state.filled + union >= state.total);537 assert!(product > 0.0);538 sum_total += state.total;539 for slot in 0..level * level {540 sum_pair[slot] += state.pair[slot];541 }542 let load = union as f64 / state.total as f64;543 load_sum += load;544 load_most = load_most.max(load);545 let independent = product * share.powi(level as i32);546 let correction = state.filled as f64 / state.total as f64 / product;547 ind_low = ind_low.min(independent);548 ind_high = ind_high.max(independent);549 ind_sum += independent;550 psi_low = psi_low.min(correction);551 psi_high = psi_high.max(correction);552 psi_sum += correction;553 }554 let scale = sum_total as f64;555 let mean: Vec<f64> = rate_sum.iter().map(|value| value / span).collect();556 let scaled: Vec<f64> = mean557 .iter()558 .enumerate()559 .map(|(j, value)| (value - null) * start as f64 / 3.0f64.powi(j as i32))560 .collect();561 let least: Vec<f64> = rate_low.iter().map(|value| trim(*value, false)).collect();562 let most: Vec<f64> = rate_high.iter().map(|value| trim(*value, true)).collect();563 let mut near: Vec<f64> = Vec::new();564 let mut far: Vec<f64> = Vec::new();565 for j in 0..level {566 for gap in 1..=2usize {567 if j + gap >= level {568 continue;569 }570 let joint = sum_pair[j * level + j + gap] as f64 / scale;571 let apart = (sum_bad[j] as f64 / scale) * (sum_bad[j + gap] as f64 / scale);572 let value = if apart > 0.0 { joint / apart } else { f64::NAN };573 if gap == 1 {574 near_low = near_low.min(value);575 near_high = near_high.max(value);576 near.push(value);577 } else {578 far_low = far_low.min(value);579 far_high = far_high.max(value);580 far.push(value);581 }582 }583 }584 for j in 0..level {585 scale_low = scale_low.min(scaled[j]);586 scale_high = scale_high.max(scaled[j]);587 if j + 3 < level {588 fine_low = fine_low.min(mean[j]);589 fine_high = fine_high.max(mean[j]);590 }591 }592 load_high = load_high.max(load_most);593 ind_least = ind_least.min(ind_low);594 ind_most = ind_most.max(ind_high);595 psi_least = psi_least.min(psi_low);596 psi_most = psi_most.max(psi_high);597 count += 1;598 println!(599 "circle-crop digits {label} k={} r={start}..{stop} level={level} null={null:.6} ind_low={:.6} ind_mean={:.6} ind_high={:.6} psi_low={:.6} psi_mean={:.6} psi_high={:.6} union_mean={:.6} union_high={:.6}",600 level - 1,601 trim(ind_low, false),602 ind_sum / span,603 trim(ind_high, true),604 trim(psi_low, false),605 psi_sum / span,606 trim(psi_high, true),607 load_sum / span,608 trim(load_most, true)609 );610 println!(611 "circle-crop digitrate {label} k={} r={start}..{stop} p_low={} p_mean={} p_high={} p_scaled={}",612 level - 1,613 list(&least),614 list(&mean),615 list(&most),616 list(&scaled)617 );618 println!(619 "circle-crop digitpair {label} k={} r={start}..{stop} rho1={} rho2={}",620 level - 1,621 list(&near),622 list(&far)623 );624 power *= 3;625 }626 println!(627 "circle-crop digittotal {label} windows={count} fine_low={:.6} fine_high={:.6} scaled_low={:.6} scaled_high={:.6} rho1_low={:.6} rho1_high={:.6} rho2_low={:.6} rho2_high={:.6} union_high={:.6} ind_low={:.6} ind_high={:.6} psi_low={:.6} psi_high={:.6}",628 trim(fine_low, false),629 trim(fine_high, true),630 trim(scale_low, false),631 trim(scale_high, true),632 trim(near_low, false),633 trim(near_high, true),634 trim(far_low, false),635 trim(far_high, true),636 trim(load_high, true),637 trim(ind_least, false),638 trim(ind_most, true),639 trim(psi_least, false),640 trim(psi_most, true)641 );642}643644// CORNER645646fn corner(tag: &str, dimension: usize, levels: &[usize], fill: u64) -> (u64, u64) {647 let mut prior: Option<Sweep> = None;648 let mut tally = (0u64, 0u64);649 for (step, level) in levels.iter().enumerate() {650 let types = design(dimension, *level);651 let side = types.shape[0];652 let r_max = (side - 1) as u64;653 let table = sweep(&types, 0, r_max);654 let label = format!("{tag} corner L={level}");655 for r in 1..=r_max {656 let at = r as usize;657 assert!(table.inside[at] <= table.seen[at]);658 assert!(table.seen[at] <= table.touch[at]);659 assert!(table.all_inside[at] <= table.all_seen[at]);660 assert!(table.all_seen[at] <= table.all_touch[at]);661 assert!(table.all_cut(r) as f64 <= cap(dimension, r));662 if dimension == 2 {663 assert_eq!(table.all_cut(r), 2 * r + 1);664 }665 let column = (r as f64 / ((dimension - 1) as f64).sqrt()).powi(dimension as i32 - 1);666 assert!(table.all_cut(r) as f64 >= column);667 }668 if let Some(past) = &prior {669 for r in 1..=past.r_max {670 let at = r as usize;671 assert_eq!(table.seen[at], past.seen[at]);672 assert_eq!(table.inside[at], past.inside[at]);673 assert_eq!(table.touch[at], past.touch[at]);674 }675 println!(676 "circle-crop levelfree {label} matches L={} on r=1..{}",677 levels[step - 1],678 past.r_max679 );680 }681 let width = match (dimension, side) {682 (_, 0..=243) => 5,683 (2, 244..=729) => 4,684 (2, 730..=2187) => 3,685 (2, 2188..=6561) => 2,686 (2, _) => 0,687 (_, _) => 0,688 };689 let sample: Vec<u64> = [13u64, 7, 5, 3, 2][..width]690 .iter()691 .map(|d| (r_max / d).max(1))692 .collect();693 oracle(&label, &types, 0, side, &table, &sample);694 if step + 1 == levels.len() {695 tally = report(&label, &table, fill, dimension);696 digit_census(&label, &table, fill, dimension);697 }698 prior = Some(table);699 }700 tally701}702703fn report(label: &str, table: &Sweep, fill: u64, dimension: usize) -> (u64, u64) {704 let r_max = table.r_max;705 let reach_max = r_max / 3;706 let width = (r_max + 2) as usize;707 let mut low = vec![0.0f64; width];708 let mut high = vec![0.0f64; width];709 let mut swing = vec![0.0f64; width];710 let mut delta = vec![0.0f64; width];711 let mut rows: Vec<String> = Vec::new();712 let mut banded = 0u64;713 for r in 1..=r_max {714 let at = r as usize;715 let mut depth = 0u32;716 let mut reach = r;717 while reach * 3 <= r_max {718 reach *= 3;719 depth += 1;720 }721 if depth > 0 {722 banded += 1;723 }724 let scale = (fill as f64).powi(depth as i32);725 let centre = table.seen[at] as f64 - table.seen[reach as usize] as f64 / scale;726 let band = table.cut(reach) as f64 / scale;727 swing[at] = centre;728 low[at] = (centre.abs() - band).max(0.0);729 high[at] = centre.abs() + band;730 assert!(low[at] <= table.cut(r) as f64);731 let mut line = format!(732 "circle-crop row {label} r={r} res={} N={} in={} cut={} Nfull={} Nfull_in={} Nfull_cut={}",733 u64::from(r.is_power_of_three()),734 table.seen[at],735 table.inside[at],736 table.cut(r),737 table.all_seen[at],738 table.all_inside[at],739 table.all_cut(r)740 );741 if r <= reach_max {742 let jump = table.seen[(r * 3) as usize] as i64 - fill as i64 * table.seen[at] as i64;743 assert!(jump.unsigned_abs() <= table.cut(r * 3) + fill * table.cut(r));744 delta[at] = jump.abs() as f64;745 line.push_str(&format!(" delta={jump}"));746 } else {747 line.push_str(" delta=na");748 }749 line.push_str(&format!(750 " E={centre:.6} Elow={:.6} Ehigh={:.6} depth={depth}",751 (low[at] * 1e6).floor() / 1e6,752 (high[at] * 1e6).ceil() / 1e6753 ));754 rows.push(line);755 }756 let mut running_delta = vec![0.0f64; width];757 let mut running_cut = vec![0.0f64; width];758 for r in 1..=r_max {759 let at = r as usize;760 running_delta[at] = running_delta[at - 1].max(delta[at]);761 running_cut[at] = running_cut[at - 1].max(table.cut(r) as f64);762 println!(763 "{} run_delta={:.0} run_cut={:.0}",764 rows[at - 1],765 running_delta[at],766 running_cut[at]767 );768 }769 trend(label, "delta", &running_delta, 27, reach_max);770 trend(label, "cut", &running_cut, 27, r_max);771 let mut power = 1u64;772 while power <= r_max {773 let at = power as usize;774 let over = if power >= 3 {775 table.seen[at] as f64 / table.seen[(power / 3) as usize] as f64776 } else {777 f64::NAN778 };779 let mut best = 0.0f64;780 let mut under = 0u64;781 let mut span = 0u64;782 let top = (power * 3 - 1).min(reach_max);783 if power <= reach_max {784 for r in power..=top {785 best = best.max(delta[r as usize]);786 span += 1;787 if delta[r as usize] <= delta[at] {788 under += 1;789 }790 }791 }792 let rank = if span > 0 {793 under as f64 / span as f64794 } else {795 f64::NAN796 };797 let mut quiet = 0u64;798 let mut reach_span = 0u64;799 for r in power..=(power * 3 - 1).min(r_max) {800 reach_span += 1;801 if table.cut(r) <= table.cut(power) {802 quiet += 1;803 }804 }805 let cut_rank = quiet as f64 / reach_span as f64;806 let defect = if power <= reach_max {807 format!(808 "delta={:.0} window_max_delta={best:.0} delta_rank={rank:.4} E={:.6}",809 delta[at], swing[at]810 )811 } else {812 String::from(813 "delta=beyond_reach window_max_delta=beyond_reach delta_rank=beyond_reach E=beyond_reach",814 )815 };816 let mass = (fill as f64).powi(power.ilog(3) as i32);817 let mu_low = (table.inside[at] as f64 / mass * 1e9).floor() / 1e9;818 let mu_high = ((table.inside[at] + table.cut(power)) as f64 / mass * 1e9).ceil() / 1e9;819 println!(820 "circle-crop resonance {label} r={power} N={} ratio={over:.6} in={} cut={} {defect} cut_rank={cut_rank:.4} mu_low={mu_low:.9} mu_high={mu_high:.9}",821 table.seen[at],822 table.inside[at],823 table.cut(power)824 );825 power *= 3;826 }827 let cuts: Vec<f64> = (0..width)828 .map(|r| table.cut((r as u64).min(r_max)) as f64)829 .collect();830 windows(831 label,832 &[833 ("cut", cuts, r_max),834 ("delta", delta, reach_max),835 ("errorlow", low, reach_max),836 ("errorhigh", high, reach_max),837 ],838 );839 means(label, table, fill, dimension);840 println!("circle-crop bands {label} rows={r_max} banded={banded}");841 (r_max, banded)842}843844trait Triadic {845 fn is_power_of_three(&self) -> bool;846}847848impl Triadic for u64 {849 fn is_power_of_three(&self) -> bool {850 let mut value = *self;851 if value == 0 {852 return false;853 }854 while value.is_multiple_of(3) {855 value /= 3;856 }857 value == 1858 }859}860861// CENTRE862863fn centre(tag: &str, dimension: usize, level: usize) -> u64 {864 let types = design(dimension, level);865 let side = types.shape[0];866 let r_max = ((side - 1) / 2) as u64;867 let table = sweep(&types, side as i64, r_max);868 let label = format!("{tag} centre L={level}");869 let cells = (side as u64).pow(dimension as u32);870 let filled = types.sum();871 println!("circle-crop density {label} side={side} cells={cells} fill={filled}");872 let sample: Vec<u64> = [11u64, 5, 3, 2]873 .iter()874 .map(|d| (r_max / d).max(1))875 .collect();876 let sample = if level >= 7 || (dimension == 3 && level >= 5) {877 sample[..2].to_vec()878 } else {879 sample880 };881 oracle(&label, &types, side as i64, side, &table, &sample);882 let mut first = 0u64;883 let mut relative = vec![0.0f64; (r_max + 2) as usize];884 for r in 1..=r_max {885 let at = r as usize;886 assert!(table.inside[at] <= table.seen[at]);887 assert!(table.seen[at] <= table.touch[at]);888 if first == 0 && table.seen[at] > 0 {889 first = r;890 }891 let exact =892 cells as i64 * table.seen[at] as i64 - filled as i64 * table.all_seen[at] as i64;893 let error = exact as f64 / cells as f64;894 let main = filled as f64 / cells as f64 * table.all_seen[at] as f64;895 relative[at] = if main > 0.0 {896 (error / main).abs()897 } else {898 0.0899 };900 println!(901 "circle-crop row {label} r={r} res={} N={} in={} cut={} Nfull={} err_num={exact} err={error:.6} rel={:.6}",902 u64::from(r.is_power_of_three()),903 table.seen[at],904 table.inside[at],905 table.cut(r),906 table.all_seen[at],907 relative[at]908 );909 }910 let hole = (side / 3 - 1) as u64 / 2;911 assert!(first > hole);912 if dimension == 2 {913 assert_eq!(first, hole + 1);914 }915 println!("circle-crop hole {label} block_inradius={hole} first_hit={first}");916 let cuts: Vec<f64> = (0..=(r_max + 1))917 .map(|r| table.cut(r.min(r_max)) as f64)918 .collect();919 windows(920 &label,921 &[("cut", cuts, r_max), ("relative", relative, r_max)],922 );923 r_max924}925926// GAUSS927928fn gauss(dimension: usize, level: usize) {929 let types = design(dimension, level);930 let side = types.shape[0];931 let r_max = (side - 1) as u64;932 let table = sweep(&types, 0, r_max);933 let pi = std::f64::consts::PI;934 let mut worst = 0.0f64;935 for r in 1..=r_max {936 let volume = if dimension == 2 {937 pi * (r as f64) * (r as f64) / 4.0938 } else {939 pi * (r as f64).powi(3) / 6.0940 };941 let gap = (table.all_seen[r as usize] as f64 - volume).abs()942 / (r as f64).powi(dimension as i32 - 1);943 worst = worst.max(gap);944 }945 println!(946 "circle-crop gauss D={dimension} L={level} r_max={r_max} max|Nfull-vol|/r^(D-1)={worst:.6}"947 );948}949950// TRANSFORM951952fn digits(dimension: usize) -> Vec<Vec<f64>> {953 let types = design(dimension, 1);954 let mut out: Vec<Vec<f64>> = Vec::new();955 let mut index = vec![0usize; dimension];956 for cell in types.bytes() {957 if *cell != 0 {958 out.push(index.iter().map(|value| *value as f64).collect());959 }960 for axis in (0..dimension).rev() {961 index[axis] += 1;962 if index[axis] < 3 {963 break;964 }965 index[axis] = 0;966 }967 }968 out969}970971fn factor(set: &[Vec<f64>], point: &[f64]) -> (f64, f64) {972 let mut real = 0.0f64;973 let mut imaginary = 0.0f64;974 for digit in set {975 let dot: f64 = digit.iter().zip(point).map(|(a, b)| a * b).sum();976 let angle = -2.0 * std::f64::consts::PI * dot;977 real += angle.cos();978 imaginary += angle.sin();979 }980 let mass = set.len() as f64;981 (real / mass, imaginary / mass)982}983984fn hat(set: &[Vec<f64>], point: &[f64], terms: u32) -> (f64, f64) {985 let mut real = 1.0f64;986 let mut imaginary = 0.0f64;987 for term in 1..=terms {988 let scale = 3.0f64.powi(term as i32);989 let scaled: Vec<f64> = point.iter().map(|value| value / scale).collect();990 let (pr, pi) = factor(set, &scaled);991 let next = (real * pr - imaginary * pi, real * pi + imaginary * pr);992 real = next.0;993 imaginary = next.1;994 }995 (real, imaginary)996}997998fn step(set: &[Vec<f64>], point: &[f64]) -> f64 {999 let dimension = point.len();1000 let mut total = 0.0f64;1001 let count = 3usize.pow(dimension as u32);1002 for code in 0..count {1003 let mut shifted = point.to_vec();1004 let mut rest = code;1005 for value in shifted.iter_mut() {1006 *value += (rest % 3) as f64 / 3.0;1007 rest /= 3;1008 }1009 let (re, im) = factor(set, &shifted);1010 total += (re * re + im * im).sqrt();1011 }1012 total1013}10141015fn mass(tag: &str, dimension: usize, top: u32, sharp: f64) {1016 let set = digits(dimension);1017 let zero = vec![0.0f64; dimension];1018 let ground = step(&set, &zero);1019 let witness: Vec<f64> = vec![155.0 / 243.0; dimension];1020 let broken = step(&set, &witness);1021 assert!((ground - 2.0).abs() < 1e-9);1022 assert!(broken < 2.0);1023 println!(1024 "circle-crop transform {tag} step h_lattice={ground:.6} h_witness={:.6} at=155/243",1025 (broken * 1e6).ceil() / 1e61026 );1027 let mut prior = f64::NAN;1028 for level in 1..=top {1029 let side = 3usize.pow(level);1030 let count = side.pow(dimension as u32);1031 let mut total = 0.0f64;1032 let mut index = vec![0usize; dimension];1033 for _ in 0..count {1034 let point: Vec<f64> = index.iter().map(|value| *value as f64).collect();1035 let mut piece = 1.0f64;1036 for term in 1..=level {1037 let scale = 3.0f64.powi(term as i32);1038 let scaled: Vec<f64> = point.iter().map(|value| value / scale).collect();1039 let (re, im) = factor(&set, &scaled);1040 piece *= (re * re + im * im).sqrt();1041 }1042 total += piece;1043 for axis in (0..dimension).rev() {1044 index[axis] += 1;1045 if index[axis] < side {1046 break;1047 }1048 index[axis] = 0;1049 }1050 }1051 assert!(total >= 2.0f64.powi(level as i32) - 1e-9);1052 let bulk = total - 1.0;1053 let ratio = bulk / prior;1054 let stride = ratio.ln() / 3.0f64.ln();1055 println!(1056 "circle-crop transform {tag} mass L={level} lambda={bulk:.6} step={ratio:.6} log3={stride:.6} cost={:.6} error={:.6}",1057 stride + 0.5,1058 stride + 0.5 + sharp - 2.01059 );1060 prior = bulk;1061 }1062}10631064fn transform(tag: &str, dimension: usize, point: &[f64]) {1065 let set = digits(dimension);1066 let tripled: Vec<f64> = point.iter().map(|value| value * 3.0).collect();1067 let (ar, ai) = hat(&set, point, 60);1068 let (br, bi) = hat(&set, &tripled, 60);1069 let (pr, pi) = factor(&set, point);1070 assert!((br - (pr * ar - pi * ai)).abs() < 1e-12);1071 assert!((bi - (pr * ai + pi * ar)).abs() < 1e-12);1072 let integral = point1073 .iter()1074 .all(|value| (value - value.round()).abs() < 1e-12);1075 let size = (pr * pr + pi * pi).sqrt();1076 assert_eq!(integral, (size - 1.0).abs() < 1e-12);1077 let name: Vec<String> = point.iter().map(|value| format!("{value}")).collect();1078 println!(1079 "circle-crop transform {tag} t=({}) integral={integral} hat_mu={:.8} hat_mu_3t={:.8} P={size:.8}",1080 name.join(","),1081 (ar * ar + ai * ai).sqrt(),1082 (br * br + bi * bi).sqrt()1083 );1084}10851086// MAIN10871088fn main() {1089 println!("circle-crop generator: CARGO_BUILD_JOBS=4 cargo run --release -p circle-crop");1090 println!("circle-crop convention: cells are indexed x in [0,3^L)^D, a cell counts when its centre x+1/2 lies in the closed ball, corner balls sit at the lattice corner 0 and centre balls at the grid centre 3^L/2");1091 gauss(2, 5);1092 gauss(3, 3);1093 transform("carpet", 2, &[0.5, 0.0]);1094 transform("carpet", 2, &[1.0, 0.0]);1095 transform("sponge", 3, &[0.5, 0.0, 0.0]);1096 transform("sponge", 3, &[1.0, 0.0, 0.0]);1097 mass("carpet", 2, 6, 8.0f64.ln() / 3.0f64.ln());1098 let carpet = corner("carpet", 2, &[5, 6, 7, 8, 9], 8);1099 let sponge = corner("sponge", 3, &[3, 4, 5, 6], 20);1100 let mut rows = carpet.0 + sponge.0;1101 let banded = carpet.1 + sponge.1;1102 rows += centre("carpet", 2, 6);1103 rows += centre("carpet", 2, 7);1104 rows += centre("sponge", 3, 4);1105 rows += centre("sponge", 3, 5);1106 transfer::transfer();1107 derive::derive();1108 println!("circle-crop totals rows={rows} banded={banded}");1109}