transfer.rs

1use crate::{list, ones_table, root_floor, trim};23// SHELL LADDER45pub(crate) struct Rung {6    pub(crate) high: Vec<u64>,7    start: Vec<usize>,8    size: usize,9}1011pub(crate) fn rung(radius: u64, level: u32) -> Rung {12    let scale = 3u64.pow(level);13    let square = radius * radius;14    let top = radius / scale;15    let mut high = Vec::with_capacity(top as usize + 2);16    for i in 0..=top + 1 {17        let reach = scale * i;18        if reach * reach > square {19            high.push(0);20        } else {21            high.push(root_floor(square - reach * reach) / scale);22        }23    }24    let mut start = Vec::with_capacity(top as usize + 1);25    let mut size = 0usize;26    for i in 0..=top as usize {27        start.push(size);28        size += (high[i] - high[i + 1] + 1) as usize;29    }30    assert_eq!(size as u64, 2 * top + 1);31    Rung { high, start, size }32}3334impl Rung {35    fn top(&self) -> usize {36        self.high.len() - 237    }3839    fn place(&self, i: usize, y: u64) -> usize {40        self.start[i] + (y - self.high[i + 1]) as usize41    }42}4344// LOCAL PATTERN4546fn pattern(here: &Rung, below: &Rung) -> Vec<u16> {47    let mut out = vec![0u16; here.size];48    let mut children = 0u64;49    for i in 0..=here.top() {50        let lo = here.high[i + 1];51        let hi = here.high[i];52        for y in lo..=hi {53            let mut mask = 0u16;54            for k in 0..3usize {55                let column = 3 * i + k;56                if column > below.top() {57                    continue;58                }59                let inner_lo = below.high[column + 1];60                let inner_hi = below.high[column];61                for t in 0..3u64 {62                    let row = 3 * y + t;63                    if row >= inner_lo && row <= inner_hi {64                        mask |= 1 << (3 * k as u16 + t as u16);65                    }66                }67            }68            assert!(mask != 0);69            children += u64::from(mask.count_ones());70            out[here.place(i, y)] = mask;71        }72    }73    assert_eq!(children, below.size as u64);74    out75}7677// EXACT RATIOS7879fn below(left: (u128, u128), right: (u128, u128)) -> bool {80    left.0 * right.1 < right.0 * left.181}8283fn ratio(pair: (u128, u128)) -> f64 {84    pair.0 as f64 / pair.1 as f6485}8687// PERRON BRACKET8889fn weights(matrix: &[Vec<u64>], rows: &[u64]) -> Vec<u128> {90    let size = matrix.len();91    let mut vector = vec![1.0f64; size];92    for _ in 0..400 {93        let mut next = vec![0.0f64; size];94        for source in 0..size {95            if rows[source] == 0 {96                continue;97            }98            let mut total = 0.0f64;99            for target in 0..size {100                total += matrix[source][target] as f64 * vector[target];101            }102            next[source] = total / rows[source] as f64;103        }104        let peak = next.iter().cloned().fold(0.0f64, f64::max);105        if peak <= 0.0 {106            break;107        }108        for slot in next.iter_mut() {109            *slot /= peak;110        }111        vector = next;112    }113    vector114        .iter()115        .map(|value| ((value * 1.0e9).round() as u128).max(1))116        .collect()117}118119fn bracket(matrix: &[Vec<u64>], rows: &[u64]) -> ((u128, u128), (u128, u128)) {120    let size = matrix.len();121    let vector = weights(matrix, rows);122    let mut low: Option<(u128, u128)> = None;123    let mut high: Option<(u128, u128)> = None;124    for source in 0..size {125        let mut total = 0u128;126        for target in 0..size {127            total += u128::from(matrix[source][target]) * vector[target];128        }129        let cell = (total, u128::from(rows[source]) * vector[source]);130        if low.is_none() || below(cell, low.unwrap()) {131            low = Some(cell);132        }133        if high.is_none() || below(high.unwrap(), cell) {134            high = Some(cell);135        }136    }137    (low.unwrap(), high.unwrap())138}139140fn living(matrix: &[Vec<u64>]) -> Vec<usize> {141    let size = matrix.len();142    let mut live: Vec<bool> = (0..size)143        .map(|source| matrix[source].iter().any(|&count| count > 0))144        .collect();145    loop {146        let mut moved = false;147        for source in 0..size {148            if !live[source] {149                continue;150            }151            if !(0..size).any(|target| live[target] && matrix[source][target] > 0) {152                live[source] = false;153                moved = true;154            }155        }156        if !moved {157            break;158        }159    }160    (0..size).filter(|&source| live[source]).collect()161}162163fn block(matrix: &[Vec<u64>], keep: &[usize]) -> Vec<Vec<u64>> {164    keep.iter()165        .map(|&source| keep.iter().map(|&target| matrix[source][target]).collect())166        .collect()167}168169// DOEBLIN170171fn doeblin(matrix: &[Vec<u64>], steps: usize) -> Vec<f64> {172    const UNIT: u128 = 1 << 40;173    let size = matrix.len();174    let rows: Vec<u128> = matrix175        .iter()176        .map(|row| row.iter().map(|&count| u128::from(count)).sum())177        .collect();178    let base: Vec<Vec<u128>> = (0..size)179        .map(|source| {180            (0..size)181                .map(|target| u128::from(matrix[source][target]) * UNIT / rows[source])182                .collect()183        })184        .collect();185    let mut power = base.clone();186    let mut out = Vec::new();187    for step in 1..=steps {188        if step > 1 {189            let mut next = vec![vec![0u128; size]; size];190            for source in 0..size {191                for middle in 0..size {192                    if power[source][middle] == 0 {193                        continue;194                    }195                    for target in 0..size {196                        next[source][target] += power[source][middle] * base[middle][target] / UNIT;197                    }198                }199            }200            power = next;201        }202        let mass: u128 = (0..size)203            .map(|target| (0..size).map(|source| power[source][target]).min().unwrap())204            .sum();205        out.push(mass as f64 / UNIT as f64);206    }207    out208}209210// DOBRUSHIN211212fn dobrushin(matrix: &[Vec<u64>]) -> (u128, u128) {213    let size = matrix.len();214    let rows: Vec<u64> = matrix.iter().map(|row| row.iter().sum()).collect();215    let mut worst = (0u128, 1u128);216    for first in 0..size {217        if rows[first] == 0 {218            continue;219        }220        for second in first + 1..size {221            if rows[second] == 0 {222                continue;223            }224            let (one, two) = (u128::from(rows[first]), u128::from(rows[second]));225            let mut gap = 0u128;226            for target in 0..size {227                let left = u128::from(matrix[first][target]) * two;228                let right = u128::from(matrix[second][target]) * one;229                gap += left.max(right) - left.min(right);230            }231            let cell = (gap, 2 * one * two);232            if below(worst, cell) {233                worst = cell;234            }235        }236    }237    worst238}239240// OPERATOR241242struct Operator {243    order: Vec<u16>,244    slot: Vec<usize>,245    all: Vec<Vec<u64>>,246    keep: Vec<Vec<u64>>,247    rows: Vec<u64>,248    wide_all: Vec<Vec<u64>>,249    wide_keep: Vec<Vec<u64>>,250    wide_rows: Vec<u64>,251    wide_live: usize,252    centres: u64,253    parents: u64,254}255256fn operator(radius: u64, deep: u32, shallow: u32) -> Operator {257    let mut rungs: Vec<Rung> = Vec::new();258    for level in 0..=shallow + 1 {259        rungs.push(rung(radius, level));260    }261    let mut masks: Vec<Vec<u16>> = vec![Vec::new()];262    for level in 1..=shallow + 1 {263        masks.push(pattern(&rungs[level as usize], &rungs[level as usize - 1]));264    }265    let mut seen = [false; 512];266    for level in deep..=shallow + 1 {267        for &mask in masks[level as usize].iter() {268            seen[mask as usize] = true;269        }270    }271    let order: Vec<u16> = (0..512u16).filter(|&mask| seen[mask as usize]).collect();272    let mut slot = vec![usize::MAX; 512];273    for (index, &mask) in order.iter().enumerate() {274        slot[mask as usize] = index;275    }276    let size = order.len();277    let mut all = vec![vec![0u64; size]; size];278    let mut keep = vec![vec![0u64; size]; size];279    let mut rows = vec![0u64; size];280    let broad = size * 9;281    let mut wide_all = vec![vec![0u64; broad]; broad];282    let mut wide_keep = vec![vec![0u64; broad]; broad];283    let mut wide_rows = vec![0u64; broad];284    let mut centres = 0u64;285    let mut parents = 0u64;286    for level in deep..=shallow {287        let here = &rungs[level as usize + 1];288        let inner = &rungs[level as usize];289        let outer_mask = &masks[level as usize + 1];290        let inner_mask = &masks[level as usize];291        for i in 0..=here.top() {292            let lo = here.high[i + 1];293            let hi = here.high[i];294            for y in lo..=hi {295                let mask = outer_mask[here.place(i, y)];296                let source = slot[mask as usize];297                let seat = 3 * (i % 3) + (y % 3) as usize;298                let wide_source = source * 9 + seat;299                rows[source] += 1;300                wide_rows[wide_source] += 1;301                parents += 1;302                if mask & (1 << 4) != 0 {303                    centres += 1;304                }305                for bit in 0..9u16 {306                    if mask & (1 << bit) == 0 {307                        continue;308                    }309                    let column = 3 * i + (bit / 3) as usize;310                    let row = 3 * y + u64::from(bit % 3);311                    let target = slot[inner_mask[inner.place(column, row)] as usize];312                    let wide_target = target * 9 + 3 * (bit / 3) as usize + (bit % 3) as usize;313                    all[source][target] += 1;314                    wide_all[wide_source][wide_target] += 1;315                    if bit != 4 {316                        keep[source][target] += 1;317                        wide_keep[wide_source][wide_target] += 1;318                    }319                }320            }321        }322    }323    let seen: Vec<usize> = (0..broad).filter(|&index| wide_rows[index] > 0).collect();324    let wide_live = seen.len();325    let wide_all = block(&wide_all, &seen);326    let wide_keep = block(&wide_keep, &seen);327    let wide_rows: Vec<u64> = seen.iter().map(|&index| wide_rows[index]).collect();328    Operator {329        order,330        slot,331        all,332        keep,333        rows,334        wide_all,335        wide_keep,336        wide_rows,337        wide_live,338        centres,339        parents,340    }341}342343// SURVIVAL LADDER344345struct Ladder {346    filled: u64,347    total: u64,348    margin: Vec<f64>,349    condition: Vec<f64>,350    gain: Vec<f64>,351}352353fn ladder(radius: u64, depth: usize, ones: &[u32]) -> Ladder {354    assert!(3u64.pow(depth as u32) > radius);355    let square = radius * radius;356    let mut filled = 0u64;357    let mut bad = vec![0u64; depth];358    let mut crest = vec![0u64; depth];359    let mut total = 0u64;360    for i in 0..=radius {361        let hi = root_floor(square - i * i);362        let step = (i + 1) * (i + 1);363        let lo = if step >= square {364            0365        } else {366            root_floor(square - step)367        };368        let first = ones[i as usize];369        for y in lo..=hi {370            let mask = first & ones[y as usize];371            total += 1;372            if mask == 0 {373                filled += 1;374                continue;375            }376            let mut rest = mask;377            while rest != 0 {378                bad[rest.trailing_zeros() as usize] += 1;379                rest &= rest - 1;380            }381            crest[(31 - mask.leading_zeros()) as usize] += 1;382        }383    }384    assert_eq!(total, 2 * radius + 1);385    let mut nested = vec![0u64; depth + 1];386    nested[0] = total;387    for cut in 1..=depth {388        let mut kept = filled;389        for position in 0..depth - cut {390            kept += crest[position];391        }392        nested[cut] = kept;393    }394    assert_eq!(nested[depth], filled);395    let mass = total as f64;396    let mut margin = Vec::with_capacity(depth);397    let mut condition = Vec::with_capacity(depth);398    let mut gain = Vec::with_capacity(depth);399    let mut product = 1.0f64;400    for position in 0..depth {401        let cut = depth - position;402        let free = (total - bad[position]) as f64 / mass;403        let tight = nested[cut] as f64 / nested[cut - 1] as f64;404        margin.push(free);405        condition.push(tight);406        gain.push(tight / free);407        product *= tight / free;408    }409    let direct = (filled as f64 / mass) / margin.iter().product::<f64>();410    assert!((product - direct).abs() <= 1.0e-12 * direct);411    assert_eq!(gain[depth - 1], 1.0);412    Ladder {413        filled,414        total,415        margin,416        condition,417        gain,418    }419}420421fn survival(label: &str, radius: u64, depth: usize) {422    let ones = ones_table(depth);423    let read = ladder(radius, depth, &ones);424    let (filled, total) = (read.filled, read.total);425    let product: f64 = read.gain.iter().product();426    let logs: Vec<f64> = read.gain.iter().map(|value| value.ln()).collect();427    let sum: f64 = logs.iter().map(|value| value.abs()).sum();428    let tail: f64 = logs[..depth.saturating_sub(3)]429        .iter()430        .map(|value| value.abs())431        .sum();432    println!(433        "circle-crop ladder {label} r={radius} L={depth} C={filled} Cfull={total} psi={:.6} logsum={:.6} logtail={:.6} share={:.6}",434        product,435        trim(sum, true),436        trim(tail, true),437        trim(tail / sum, true)438    );439    println!(440        "circle-crop ladder {label} r={radius} marginal={} conditional={} gain={}",441        list(&read.margin),442        list(&read.condition),443        list(&read.gain)444    );445    let scaled: Vec<f64> = logs.iter().map(|value| value * 1.0e3).collect();446    println!(447        "circle-crop ladder {label} r={radius} loggain_x1000={}",448        list(&scaled)449    );450    let top = read.margin[depth - 1];451    let next = read.margin[depth - 2];452    println!(453        "circle-crop ladder {label} r={radius} g0={:.12} g1={:.12} top_marginal={:.12} next_marginal={:.12} both_centres_crossed={}",454        read.gain[depth - 1],455        read.gain[depth - 2],456        top,457        next,458        top < 1.0 && next < 1.0459    );460}461462// PROFILE463464struct Gather {465    depth: usize,466    rows: u64,467    logs: Vec<f64>,468    sums: Vec<f64>,469    sizes: Vec<f64>,470    counts: Vec<u64>,471    flat: u64,472}473474fn gather(window: u32, step: u64) -> Gather {475    let depth = window as usize + 1;476    let ones = ones_table(depth);477    let low = 3u64.pow(window);478    let high = 3u64.pow(window + 1) - 1;479    let mut sums = vec![0.0f64; depth];480    let mut sizes = vec![0.0f64; depth];481    let mut counts = vec![0u64; depth];482    let mut logs = Vec::new();483    let mut rows = 0u64;484    let mut flat = 0u64;485    let mut radius = low;486    while radius <= high {487        let read = ladder(radius, depth, &ones);488        let mut total = 0.0f64;489        for position in 0..depth {490            let value = read.gain[position].ln();491            let rank = depth - 1 - position;492            sums[rank] += value;493            sizes[rank] += value.abs();494            counts[rank] += 1;495            total += value;496        }497        if read.gain[depth - 2] == 1.0 {498            flat += 1;499        }500        logs.push(total);501        rows += 1;502        radius += step;503    }504    Gather {505        depth,506        rows,507        logs,508        sums,509        sizes,510        counts,511        flat,512    }513}514515fn bulk(read: &Gather) -> f64 {516    (0..read.depth)517        .map(|rank| read.sizes[rank] / read.counts[rank] as f64 * 1.0e3)518        .sum()519}520521fn profile(label: &str, window: u32, samples: u64) -> (f64, f64, f64) {522    let low = 3u64.pow(window);523    let high = 3u64.pow(window + 1) - 1;524    let step = ((high - low + 1) / samples).max(1);525    let read = gather(window, step);526    let depth = read.depth;527    let rows = read.rows;528    let mean: f64 = read.logs.iter().sum::<f64>() / rows as f64;529    let spread = (read530        .logs531        .iter()532        .map(|value| (value - mean).powi(2))533        .sum::<f64>()534        / rows as f64)535        .sqrt();536    let worst = read.logs.iter().cloned().fold(f64::MIN, f64::max);537    let best = read.logs.iter().cloned().fold(f64::MAX, f64::min);538    let signed: Vec<f64> = (0..depth)539        .map(|rank| read.sums[rank] / read.counts[rank] as f64 * 1.0e3)540        .collect();541    let sized: Vec<f64> = (0..depth)542        .map(|rank| read.sizes[rank] / read.counts[rank] as f64 * 1.0e3)543        .collect();544    let deep: f64 = if depth > 6 {545        read.sums[6..].iter().sum::<f64>() / read.counts[6..].iter().sum::<u64>() as f64 * 1.0e3546    } else {547        0.0548    };549    let deep_size: f64 = if depth > 6 {550        read.sizes[6..].iter().sum::<f64>() / read.counts[6..].iter().sum::<u64>() as f64 * 1.0e3551    } else {552        0.0553    };554    println!(555        "circle-crop profile {label} window={low}..{high} L={depth} radii={rows} logpsi mean={:.6} sd={:.6} min={:.6} max={:.6} psi=[{:.6},{:.6}]",556        mean,557        trim(spread, true),558        trim(best, false),559        trim(worst, true),560        trim(best.exp(), false),561        trim(worst.exp(), true)562    );563    println!(564        "circle-crop profile {label} window={low}..{high} loggain_x1000={} absgain_x1000={}",565        list(&signed),566        list(&sized)567    );568    println!(569        "circle-crop profile {label} window={low}..{high} deep_signed_x1000={:.6} deep_abs_x1000={:.6} deep_ranks={}",570        deep,571        deep_size,572        depth.saturating_sub(6)573    );574    let total = bulk(&read);575    let steps: Vec<f64> = (2..depth)576        .map(|rank| sized[rank] / sized[rank - 1])577        .collect();578    println!(579        "circle-crop profile {label} window={low}..{high} decay={} tail_decay={:.6} rank1_flat={} of {rows}",580        list(&steps),581        trim(steps[1..].iter().cloned().fold(0.0f64, f64::max), true),582        read.flat583    );584    (total, trim(best.exp(), false), trim(worst.exp(), true))585}586587fn tail(totals: &[f64]) -> (Vec<f64>, Vec<f64>, f64, f64) {588    let steps: Vec<f64> = (1..totals.len())589        .map(|index| totals[index] - totals[index - 1])590        .collect();591    let decay: Vec<f64> = (1..steps.len())592        .map(|index| steps[index] / steps[index - 1])593        .collect();594    let worst = decay.iter().cloned().fold(0.0f64, f64::max);595    let last = *steps.last().unwrap();596    let cap = totals.last().unwrap() + last * worst / (1.0 - worst);597    (steps, decay, worst, cap)598}599600fn census(label: &str, totals: &[f64], floor: f64, roof: f64) {601    let (steps, decay, _, cap) = tail(totals);602    let band = (603        trim((-cap / 1000.0).exp(), false),604        trim((cap / 1000.0).exp(), true),605    );606    println!(607        "circle-crop profile {label} totals_x1000={} steps_x1000={} decay={} cap_x1000={:.6} psi=[{:.6},{:.6}] extremes=[{:.6},{:.6}] band_holds={}",608        list(totals),609        list(&steps),610        list(&decay),611        trim(cap, true),612        band.0,613        band.1,614        floor,615        roof,616        band.0 <= floor && roof <= band.1617    );618}619620fn resample(label: &str, windows: &[u32], strides: &[u64], head: &[f64]) {621    for &step in strides {622        let mut totals: Vec<f64> = head.to_vec();623        let mut radii = Vec::new();624        for &window in windows {625            let read = gather(window, step);626            radii.push(read.rows as f64);627            totals.push(bulk(&read));628        }629        let (steps, _, worst, cap) = tail(&totals);630        println!(631            "circle-crop profile {label} stride={step} radii={} totals_x1000={} last_step_x1000={:.6} worst_decay={:.6} cap_x1000={:.6} converges={}",632            list(&radii),633            list(&totals),634            *steps.last().unwrap(),635            trim(worst, true),636            trim(cap, true),637            worst < 1.0638        );639    }640}641642// REPORT643644fn spectrum(label: &str, tag: &str, matrix: &[Vec<u64>], rows: &[u64], target: f64) -> (f64, f64) {645    let (low, high) = bracket(matrix, rows);646    let (bottom, top) = (trim(ratio(low), false), trim(ratio(high), true));647    let holds = bottom <= target && target <= top;648    println!(649        "circle-crop operator {label} {tag} states={} perron=[{bottom:.6},{top:.6}] target={target:.6} covers={holds}",650        matrix.len()651    );652    (bottom, top)653}654655fn report(label: &str, radius: u64, deep: u32, shallow: u32) -> Vec<u16> {656    let built = operator(radius, deep, shallow);657    let size = built.order.len();658    let branches: Vec<u64> = built659        .order660        .iter()661        .map(|&mask| u64::from(mask.count_ones()))662        .collect();663    let mean: f64 = built664        .rows665        .iter()666        .zip(branches.iter())667        .map(|(&count, &arms)| (count * arms) as f64)668        .sum::<f64>()669        / built.parents as f64;670    let centre = built.centres as f64 / built.parents as f64;671    println!(672        "circle-crop operator {label} r={radius} levels={deep}..{shallow} states={size} parents={} meanb={:.6} centre={:.6} centre_x3={:.6}",673        built.parents,674        mean,675        centre,676        centre * 3.0677    );678    let widest = built679        .order680        .iter()681        .max_by_key(|&&mask| built.rows[built.slot[mask as usize]])682        .unwrap();683    println!(684        "circle-crop operator {label} r={radius} masks={} shares={}",685        built686            .order687            .iter()688            .map(|mask| format!("{mask:03o}"))689            .collect::<Vec<String>>()690            .join(","),691        list(692            &built693                .rows694                .iter()695                .map(|&count| count as f64 / built.parents as f64)696                .collect::<Vec<f64>>()697        )698    );699    println!(700        "circle-crop operator {label} r={radius} commonest={widest:03o} share={:.6}",701        built.rows[built.slot[*widest as usize]] as f64 / built.parents as f64702    );703    let whole = spectrum(label, "all", &built.all, &built.rows, 3.0);704    let alive = living(&built.keep);705    let trimmed = block(&built.keep, &alive);706    let kept_rows: Vec<u64> = alive.iter().map(|&index| built.rows[index]).collect();707    let pruned = spectrum(label, "pruned", &trimmed, &kept_rows, 8.0 / 3.0);708    let share = (709        trim(pruned.0 / whole.1, false),710        trim(pruned.1 / whole.0, true),711    );712    let drift = (713        trim((share.0 * 9.0 / 8.0).ln() / 3.0f64.ln(), false),714        trim((share.1 * 9.0 / 8.0).ln() / 3.0f64.ln(), true),715    );716    println!(717        "circle-crop operator {label} r={radius} share=[{:.6},{:.6}] target={:.6} drift=[{:.6},{:.6}]",718        share.0,719        share.1,720        8.0 / 9.0,721        drift.0,722        drift.1723    );724    let wide_whole = spectrum(label, "wide-all", &built.wide_all, &built.wide_rows, 3.0);725    let wide_alive = living(&built.wide_keep);726    let wide_trimmed = block(&built.wide_keep, &wide_alive);727    let wide_kept: Vec<u64> = wide_alive728        .iter()729        .map(|&index| built.wide_rows[index])730        .collect();731    let wide_pruned = spectrum(label, "wide-pruned", &wide_trimmed, &wide_kept, 8.0 / 3.0);732    let wide_share = (733        trim(wide_pruned.0 / wide_whole.1, false),734        trim(wide_pruned.1 / wide_whole.0, true),735    );736    println!(737        "circle-crop operator {label} r={radius} wide_states={} wide_share=[{:.6},{:.6}]",738        built.wide_live, wide_share.0, wide_share.1739    );740    let mixing = dobrushin(&built.all);741    let floors = doeblin(&built.all, 6);742    let least = floors743        .iter()744        .enumerate()745        .map(|(step, &mass)| (1.0 - mass).powf(1.0 / (step + 1) as f64))746        .fold(f64::MAX, f64::min);747    println!(748        "circle-crop operator {label} r={radius} dobrushin={:.6} doeblin={} rate={:.6}",749        trim(ratio(mixing), true),750        list(751            &floors752                .iter()753                .map(|mass| trim(*mass, false))754                .collect::<Vec<f64>>()755        ),756        trim(least, true)757    );758    built.order759}760761// ALPHABET762763pub(crate) fn levels(radius: u64) -> u32 {764    let mut level = 0u32;765    while 3u64.pow(level) <= radius {766        level += 1;767    }768    level769}770771pub(crate) fn masks_seen(radius: u64, deep: u32, shallow: u32) -> Vec<u16> {772    let mut seen = [false; 512];773    for level in deep..=shallow {774        let here = rung(radius, level);775        let under = rung(radius, level - 1);776        for &mask in pattern(&here, &under).iter() {777            seen[mask as usize] = true;778        }779    }780    (0..512u16).filter(|&mask| seen[mask as usize]).collect()781}782783pub(crate) fn show(masks: &[u16]) -> String {784    masks785        .iter()786        .map(|mask| format!("{mask:03o}"))787        .collect::<Vec<String>>()788        .join(",")789}790791fn extra(seen: &[u16], base: &[u16]) -> Vec<u16> {792    seen.iter()793        .filter(|mask| !base.contains(mask))794        .cloned()795        .collect()796}797798fn whole(label: &str, radius: u64, base: &[u16]) {799    let depth = levels(radius);800    let seen = masks_seen(radius, 1, depth);801    let cut = masks_seen(radius, 1, depth - 3);802    assert_eq!(cut, base);803    let more = extra(&seen, base);804    let root = pattern(&rung(radius, depth), &rung(radius, depth - 1))[0];805    println!(806        "circle-crop operator {label} r={radius} L={depth} every_level_states={} truncated_states={} extra={} root_state={root:03o}",807        seen.len(),808        cut.len(),809        show(&more)810    );811}812813fn pinned(label: &str, radius: u64, base: &[u16], want: &[f64]) {814    let depth = levels(radius);815    let mut widths = Vec::new();816    let mut seats = [false; 512];817    for back in 1..=4u32 {818        let seen = masks_seen(radius, 1, depth - back);819        widths.push(seen.len() as f64);820        for &mask in extra(&seen, base).iter() {821            seats[mask as usize] = true;822        }823    }824    assert_eq!(widths, want);825    let found: Vec<u16> = (0..512u16).filter(|&mask| seats[mask as usize]).collect();826    let cut = masks_seen(radius, 1, depth - 3);827    let more = extra(&cut, base);828    let mut first = 0u32;829    if !more.is_empty() {830        for level in 1..=depth - 3 {831            if pattern(&rung(radius, level), &rung(radius, level - 1)).contains(&more[0]) {832                first = level;833                break;834            }835        }836    }837    println!(838        "circle-crop operator {label} r={radius} L={depth} states_at_L_less_1_to_4={} truncated_states={} extra={} first_level={first} all_extras={}",839        list(&widths),840        cut.len(),841        show(&more),842        show(&found)843    );844}845846fn scan(label: &str, low: u64, high: u64, back: u32, base: &[u16]) {847    let mut rows = 0u64;848    let mut over = 0u64;849    let mut under = 0u64;850    let mut seats = [false; 512];851    let mut witness: Vec<u64> = Vec::new();852    for radius in low..=high {853        let depth = levels(radius);854        if depth < back + 2 {855            continue;856        }857        let seen = masks_seen(radius, 1, depth - back);858        rows += 1;859        let more = extra(&seen, base);860        if !more.is_empty() {861            over += 1;862            for &mask in more.iter() {863                seats[mask as usize] = true;864            }865            if witness.len() < 6 {866                witness.push(radius);867            }868        }869        if seen.len() - more.len() < base.len() {870            under += 1;871        }872    }873    let found: Vec<u16> = (0..512u16).filter(|&mask| seats[mask as usize]).collect();874    println!(875        "circle-crop operator {label} scan={low}..{high} levels=1..L-{back} radii={rows} over={over} under={under} extras={} first={}",876        show(&found),877        witness878            .iter()879            .map(|value| value.to_string())880            .collect::<Vec<String>>()881            .join(",")882    );883}884885// MEMORY ONE886887fn memory(label: &str, radius: u64, deep: u32, shallow: u32, start: u32) {888    let built = operator(radius, deep, shallow);889    let size = built.order.len();890    for source in 0..size {891        let arms = u64::from(built.order[source].count_ones());892        let row: u64 = built.all[source].iter().sum();893        assert_eq!(row, built.rows[source] * arms);894    }895    let masks = pattern(&rung(radius, start), &rung(radius, start - 1));896    let mut count = vec![0.0f64; size];897    for &mask in masks.iter() {898        count[built.slot[mask as usize]] += 1.0;899    }900    let mut model = Vec::new();901    let mut truth = Vec::new();902    let mut branch = Vec::new();903    let mut level = start;904    while level >= 1 {905        let mut next = vec![0.0f64; size];906        for source in 0..size {907            if count[source] == 0.0 || built.rows[source] == 0 {908                continue;909            }910            let share = count[source] / built.rows[source] as f64;911            for target in 0..size {912                next[target] += share * built.all[source][target] as f64;913            }914        }915        model.push(next.iter().sum::<f64>());916        let big = (2 * (radius / 3u64.pow(level - 1)) + 1) as f64;917        let small = (2 * (radius / 3u64.pow(level)) + 1) as f64;918        truth.push(big);919        branch.push(big / small);920        count = next;921        level -= 1;922    }923    let error = model924        .iter()925        .zip(truth.iter())926        .map(|(read, want)| (read - want).abs())927        .fold(0.0f64, f64::max);928    let mut step = Vec::new();929    let mut slip = 0.0f64;930    for back in (1..=start).rev() {931        let seen = pattern(&rung(radius, back), &rung(radius, back - 1));932        let mut heads = vec![0.0f64; size];933        for &mask in seen.iter() {934            heads[built.slot[mask as usize]] += 1.0;935        }936        let mass: f64 = (0..size)937            .map(|source| heads[source] * f64::from(built.order[source].count_ones()))938            .sum();939        step.push(mass);940        slip = slip.max((mass - (2 * (radius / 3u64.pow(back - 1)) + 1) as f64).abs());941    }942    println!(943        "circle-crop operator {label} r={radius} rowsum_is_popcount=true start={start} onestep={} truth={} onestep_err={:.9} iterated={} iterated_maxerr={:.9} branching={}",944        list(&step),945        list(&truth),946        slip,947        list(&model),948        error,949        list(&branch)950    );951}952953fn deepen(label: &str, radius: u64, top: u32) {954    let mut gaps = Vec::new();955    for back in (1..=top).rev() {956        let built = operator(radius, back, back);957        let (low, high) = bracket(&built.all, &built.rows);958        let (bottom, roof) = (trim(ratio(low), false), trim(ratio(high), true));959        let gap = if bottom <= 3.0 && 3.0 <= roof {960            0.0961        } else {962            (bottom - 3.0).abs().min((roof - 3.0).abs())963        };964        gaps.push(trim(gap, false));965        println!(966            "circle-crop operator {label} r={radius} pair={}..{back} parents={} states={} perron=[{bottom:.6},{roof:.6}] gap={:.6}",967            back + 1,968            built.parents,969            built.all.len(),970            trim(gap, false)971        );972    }973    println!(974        "circle-crop operator {label} r={radius} gap_by_pair={} monotone={}",975        list(&gaps),976        gaps.windows(2).all(|pair| pair[1] <= pair[0])977    );978}979980fn signs(label: &str, radii: &[u64]) {981    let mut above = 0u32;982    let mut misses = 0u32;983    let mut least = f64::MAX;984    let mut most = 0.0f64;985    let mut rows = Vec::new();986    for &radius in radii {987        let depth = levels(radius);988        let built = operator(radius, 1, depth - 4);989        let (low, high) = bracket(&built.all, &built.rows);990        let (bottom, roof) = (trim(ratio(low), false), trim(ratio(high), true));991        if !(bottom <= 3.0 && 3.0 <= roof) {992            misses += 1;993        }994        let middle = (bottom + roof) / 2.0 - 3.0;995        if middle > 0.0 {996            above += 1;997        }998        least = least.min(middle.abs());999        most = most.max(middle.abs());1000        rows.push(middle * 1.0e3);1001    }1002    println!(1003        "circle-crop operator {label} sweep_radii={} perron_less_3_x1000={} above={above} below={} excludes_3={misses} size_x1000=[{:.6},{:.6}]",1004        radii.len(),1005        list(&rows),1006        radii.len() as u32 - above,1007        trim(least * 1.0e3, false),1008        trim(most * 1.0e3, true)1009    );1010}10111012// INDEX BOUND10131014fn floor_shift(square: i64, back: i64, step: i64) -> i64 {1015    let root = if square <= 0 {1016        01017    } else {1018        root_floor(square as u64) as i641019    };1020    (root - back).div_euclid(step)1021}10221023fn index_cap(reach: f64) -> f64 {1024    13.60 * reach.powf(-1.0 / 3.0) + 305.08 / reach.sqrt() + 9.84 / reach1025}10261027struct Index {1028    rate: Vec<f64>,1029    cap: Vec<f64>,1030    slack: f64,1031    drift: f64,1032    logind: f64,1033    live: usize,1034}10351036fn index(radius: u64) -> Index {1037    let square = radius * radius;1038    let mut high = Vec::with_capacity(radius as usize + 2);1039    for x in 0..=radius + 1 {1040        high.push(if x * x > square {1041            01042        } else {1043            root_floor(square - x * x)1044        });1045    }1046    let depth = levels(radius);1047    let shell = 2 * radius + 1;1048    let mut rate = Vec::new();1049    let mut cap = Vec::new();1050    let mut slack = 0.0f64;1051    let mut drift = 0.0f64;1052    let mut logind = 0.0f64;1053    let mut live = 0usize;1054    for level in 0..depth {1055        let scale = 3u64.pow(level);1056        let step = 3 * scale;1057        let last = radius / scale;1058        let mut cells = 0u64;1059        let mut wide = 0u64;1060        let mut tall = 0u64;1061        let mut seats = 0u64;1062        let mut seat_cells = 0u64;1063        let mut seat_wide = 0u64;1064        let mut seat_tall = 0u64;1065        for c in 0..=last {1066            let head = high[(scale * c) as usize];1067            let foot = if scale * (c + 1) > radius {1068                01069            } else {1070                high[(scale * (c + 1)) as usize]1071            };1072            let low = foot / scale;1073            let top = head / scale;1074            let span = (top - low + 1) as usize;1075            let mut runs = vec![0u64; span];1076            let mut load = vec![0u64; span];1077            for x in scale * c..=(scale * (c + 1) - 1).min(radius) {1078                let up = high[x as usize] / scale;1079                let down = high[(x + 1) as usize] / scale;1080                for t in down..=up {1081                    let seat = (t - low) as usize;1082                    runs[seat] += 1;1083                    let roof = high[x as usize].min(scale * (t + 1) - 1);1084                    let base = high[(x + 1) as usize].max(scale * t);1085                    load[seat] += roof - base + 1;1086                }1087            }1088            for t in low..=top {1089                let seat = (t - low) as usize;1090                let rows = head.min(scale * (t + 1) - 1) - foot.max(scale * t) + 1;1091                assert_eq!(load[seat], runs[seat] + rows - 1);1092                cells += load[seat];1093                wide += runs[seat];1094                tall += rows;1095                if c % 3 == 1 && t % 3 == 1 {1096                    seats += 1;1097                    seat_cells += load[seat];1098                    seat_wide += runs[seat];1099                    seat_tall += rows;1100                }1101            }1102        }1103        assert_eq!(cells, shell);1104        assert_eq!(wide, last + radius + 1);1105        assert_eq!(tall, last + radius + 1);1106        assert_eq!(seat_wide, seat_tall);1107        assert_eq!(seat_cells, 2 * seat_wide - seats);1108        let mut form_wide = 0i64;1109        for x in 0..=radius {1110            if (x / scale) % 3 != 1 {1111                continue;1112            }1113            let near = square as i64 - (x * x) as i64;1114            let far = square as i64 - ((x + 1) * (x + 1)) as i64;1115            form_wide += floor_shift(near, scale as i64, step as i64)1116                - floor_shift(far, 2 * scale as i64, step as i64);1117        }1118        assert_eq!(form_wide, seat_wide as i64);1119        let mut form_seats = 0i64;1120        for c in 0..=last {1121            if c % 3 != 1 {1122                continue;1123            }1124            let near = square as i64 - (scale * c).pow(2) as i64;1125            let far = square as i64 - (scale * (c + 1)).pow(2) as i64;1126            form_seats += floor_shift(near, scale as i64, step as i64)1127                - floor_shift(far, 2 * scale as i64, step as i64);1128        }1129        assert_eq!(form_seats, seats as i64);1130        assert!(shell - seat_cells >= scale);1131        let share = seat_cells as f64 / shell as f64;1132        let reach = radius as f64 / scale as f64;1133        let roof = index_cap(reach);1134        let gap = (share - 1.0 / 9.0).abs();1135        assert!(gap <= roof);1136        if roof < 8.0 / 9.0 {1137            live += 1;1138        }1139        slack = slack.max(gap / roof);1140        drift += gap;1141        logind += (1.0 - share).ln() - (8.0f64 / 9.0).ln();1142        rate.push(share);1143        cap.push(roof);1144    }1145    assert!(drift <= 781.0);1146    assert!(logind.abs() <= 1191.0);1147    Index {1148        rate,1149        cap,1150        slack,1151        drift,1152        logind,1153        live,1154    }1155}11561157fn live(label: &str, radius: u64) {1158    let square = radius * radius;1159    let mut seats = 0i64;1160    let mut x = 1u64;1161    while x <= radius {1162        let near = square as i64 - (x * x) as i64;1163        let far = square as i64 - ((x + 1) * (x + 1)) as i64;1164        seats += floor_shift(near, 1, 3) - floor_shift(far, 2, 3);1165        x += 3;1166    }1167    assert!(seats >= 0);1168    let shell = 2 * radius + 1;1169    let share = seats as f64 / shell as f64;1170    let roof = index_cap(radius as f64);1171    let gap = (share - 1.0 / 9.0).abs();1172    assert!(roof < 8.0 / 9.0);1173    assert!(gap <= roof);1174    println!(1175        "circle-crop index {label} live r={radius} j=0 seats={seats} p0={share:.9} gap={} cap={:.9} ratio={} headroom={:.9}",1176        trim(gap, true),1177        (roof * 1e9).ceil() / 1e9,1178        trim(gap / roof, true),1179        ((8.0 / 9.0 - roof) * 1e9).floor() / 1e91180    );1181}11821183fn indexed(label: &str, radius: u64) -> usize {1184    let read = index(radius);1185    println!(1186        "circle-crop index {label} r={radius} L={} live_levels={} of={} drift={} cap=781.000000 logind={} cap=1191.000000 slack={} rate={} bound={}",1187        read.rate.len(),1188        read.live,1189        read.rate.len(),1190        trim(read.drift, true),1191        trim(read.logind.abs(), true),1192        trim(read.slack, true),1193        list(&read.rate),1194        list(&read.cap)1195    );1196    read.live1197}11981199// ENTRY12001201pub fn transfer() {1202    println!("circle-crop transfer convention: every line below is the carpet, the crossing shell at level j is the whole grid's shell at the real radius r/3^j, a box's state is the 9-bit pattern of its crossed children printed in octal with bit 4 the centre, every matrix entry is an exact count of parent-child pairs over the printed levels, a gain is indexed by the depth from the top k = L - 1 - j so g_k = (T_(k+1)/T_k)/(1 - p_(L-1-k)), and a profile line samples its window at an even stride");1203    for &radius in &[80u64, 242, 1000, 2186, 6560, 19682, 12345] {1204        survival("carpet", radius, levels(radius) as usize);1205    }1206    let mut totals = Vec::new();1207    let mut floor = f64::MAX;1208    let mut roof = 0.0f64;1209    for window in 3..=8u32 {1210        let (total, best, worst) = profile("carpet", window, 96);1211        totals.push(total);1212        floor = floor.min(best);1213        roof = roof.max(worst);1214    }1215    census("carpet", &totals, floor, roof);1216    resample("carpet", &[6, 7, 8], &[24, 48, 72], &totals[..3]);1217    let mut alphabet: Option<Vec<u16>> = None;1218    for &(radius, deep, shallow) in &[(6560u64, 1u32, 4u32), (19682, 1, 5), (12345, 1, 5)] {1219        let order = report("carpet", radius, deep, shallow);1220        match alphabet {1221            None => alphabet = Some(order),1222            Some(ref first) => assert_eq!(first, &order),1223        }1224    }1225    let base = alphabet.unwrap();1226    println!(1227        "circle-crop operator carpet alphabet={} states={} levels=1..L-3",1228        show(&base),1229        base.len()1230    );1231    memory("carpet", 6560, 1, 4, 5);1232    deepen("carpet", 6560, 4);1233    signs(1234        "carpet",1235        &[1236            3001, 3901, 4801, 5701, 6601, 7501, 8401, 9301, 10201, 11101, 12001, 12901, 13801,1237            14701, 15601, 16501, 17401, 18301,1238        ],1239    );1240    for &radius in &[6560u64, 19682, 12345] {1241        whole("carpet", radius, &base);1242    }1243    for &(radius, want) in &[1244        (1395u64, [31.0f64, 30.0, 30.0, 30.0]),1245        (1739, [31.0, 30.0, 30.0, 30.0]),1246        (6570, [31.0, 30.0, 30.0, 30.0]),1247        (15122, [31.0, 31.0, 31.0, 30.0]),1248        (3182, [31.0, 31.0, 31.0, 31.0]),1249    ] {1250        pinned("carpet", radius, &base, &want);1251    }1252    scan("carpet", 3000, 19682, 3, &base);1253    scan("carpet", 3000, 19682, 1, &base);1254    let mut live_levels = 0usize;1255    for &radius in &[80u64, 242, 1000, 2186, 6560, 12345, 19682] {1256        live_levels += indexed("carpet", radius) as usize;1257    }1258    println!(1259        "circle-crop index carpet totals radii=7 live_levels={live_levels} note=every level of these seven radii has cap >= 8/9, so the assert cannot fail there and the numbers are a check on the identities and not on the bound; the bound bites only at R >= 212957 and beats 1/9 only at R >= 23157375"1260    );1261    for &radius in &[212957u64, 531441, 2000000] {1262        live("carpet", radius);1263    }1264}