derive.rs

1use crate::root_floor;2use crate::transfer::{levels, masks_seen, rung, show};34// PATTERN LAW56fn mask_of(steps: [i64; 4], reach: usize) -> u16 {7    let mut mask = 0u16;8    for column in 0..reach.min(3) {9        let low = steps[column + 1].max(0);10        let high = steps[column].min(2);11        let mut row = low;12        while row <= high {13            mask |= 1 << (3 * column as u16 + row as u16);14            row += 1;15        }16    }17    mask18}1920fn window(steps: [i64; 4]) -> Option<(i64, i64)> {21    let mut low = i64::MIN;22    let mut high = i64::MAX;23    for near in 0..4usize {24        for far in 0..4usize {25            if far > near {26                let gap = (far - near) as i64;27                low = low.max(6 * (steps[near] - steps[far] - 1) / gap);28            }29            if far < near {30                let gap = (near - far) as i64;31                high = high.min(6 * (steps[far] - steps[near] + 1) / gap);32            }33        }34    }35    let floor = low.max(0);36    if high > floor {37        Some((floor, high))38    } else {39        None40    }41}4243// LINE ALPHABET4445fn line_alphabet(span: i64) -> Vec<(u16, i64, i64)> {46    let mut found: Vec<Option<(i64, i64)>> = vec![None; 512];47    for a in -span..=span {48        for b in -span..=span {49            for c in -span..=span {50                for d in -span..=span {51                    let steps = [a, b, c, d];52                    let Some((low, high)) = window(steps) else {53                        continue;54                    };55                    let mask = mask_of(steps, 3);56                    if mask == 0 {57                        continue;58                    }59                    let slot = &mut found[mask as usize];60                    match slot {61                        None => *slot = Some((low, high)),62                        Some(range) => {63                            range.0 = range.0.min(low);64                            range.1 = range.1.max(high);65                        }66                    }67                }68            }69        }70    }71    (0..512)72        .filter_map(|mask| found[mask].map(|range| (mask as u16, range.0, range.1)))73        .collect()74}7576// SHELL AGAINST THE LINE7778fn compare(label: &str, radius: u64) {79    let deep = levels(radius);80    let seen = masks_seen(radius, 1, deep - 3);81    let all = masks_seen(radius, 1, deep);82    let alphabet = line_alphabet(9);83    let line: Vec<u16> = alphabet.iter().map(|entry| entry.0).collect();84    let outside: Vec<u16> = seen.iter().cloned().filter(|m| !line.contains(m)).collect();85    let missing: Vec<u16> = line.iter().cloned().filter(|m| !seen.contains(m)).collect();86    let curved: Vec<u16> = all.iter().cloned().filter(|m| !line.contains(m)).collect();87    assert!(outside.is_empty() && missing.is_empty());88    assert_eq!(curved.len(), 1);89    println!(90        "alphabet {} r={} line={} shell={} shell_outside_line=[{}] line_missing_from_shell=[{}] every_level_outside_line=[{}]",91        label,92        radius,93        line.len(),94        seen.len(),95        show(&outside),96        show(&missing),97        show(&curved)98    );99}100101// FROZEN SLOPE102103fn frozen(label: &str, radius: u64) {104    let deep = levels(radius);105    for level in 0..deep.saturating_sub(1) {106        let below = rung(radius, level);107        let here = rung(radius, level + 1);108        let scale = 3u64.pow(level) as f64;109        let far = (radius * radius) as f64;110        let top_here = here.high.len() - 2;111        let top_below = below.high.len() - 2;112        let mut columns = 0u64;113        let mut column_hits = 0u64;114        let mut boxes = 0u64;115        let mut faults = 0u64;116        let mut chord_faults = 0u64;117        let mut shallow = 0u64;118        let mut shallow_faults = 0u64;119        let mut worst = 0usize;120        for i in 0..=top_here {121            columns += 1;122            let edge = (3 * i) as f64 * scale;123            let height = (far - edge * edge).max(0.0).sqrt() / scale;124            if height <= 0.0 {125                continue;126            }127            let slope = (3 * i) as f64 / height;128            let far_edge = (3 * (i + 1)) as f64 * scale;129            let end = (far - far_edge * far_edge).max(0.0).sqrt() / scale;130            let chord = (height - end) / 3.0;131            let mut guess = [0i64; 4];132            let mut rope = [0i64; 4];133            guess[0] = below.high[3 * i] as i64;134            rope[0] = guess[0];135            for k in 1..4usize {136                guess[k] = (height - k as f64 * slope).floor() as i64;137                rope[k] = (height - k as f64 * chord).floor() as i64;138            }139            let reach = (top_below + 1).saturating_sub(3 * i).min(3);140            let mut exact = true;141            for k in 1..4usize {142                if 3 * i + k <= top_below + 1 && guess[k] != below.high[3 * i + k] as i64 {143                    exact = false;144                }145            }146            if exact {147                column_hits += 1;148            }149            let low = here.high[i + 1];150            let high = here.high[i];151            for y in low..=high {152                boxes += 1;153                let mut truth = [0i64; 4];154                for k in 0..4usize {155                    let column = 3 * i + k;156                    truth[k] = if column <= top_below + 1 {157                        below.high[column] as i64 - 3 * y as i64158                    } else {159                        -1160                    };161                }162                let shifted = [163                    guess[0] - 3 * y as i64,164                    guess[1] - 3 * y as i64,165                    guess[2] - 3 * y as i64,166                    guess[3] - 3 * y as i64,167                ];168                let roped = [169                    rope[0] - 3 * y as i64,170                    rope[1] - 3 * y as i64,171                    rope[2] - 3 * y as i64,172                    rope[3] - 3 * y as i64,173                ];174                let seen = mask_of(truth, reach);175                if seen != mask_of(shifted, reach) {176                    faults += 1;177                    worst = worst.max(i);178                }179                let missed = seen != mask_of(roped, reach);180                if missed {181                    chord_faults += 1;182                }183                if slope <= 1.0 {184                    shallow += 1;185                    if missed {186                        shallow_faults += 1;187                    }188                }189            }190        }191        println!(192            "frozen {} r={} level={} columns={} column_exact={} boxes={} tangent_faults={} tangent_rate={:.6} chord_faults={} chord_rate={:.6} shallow_boxes={} shallow_chord_faults={} shallow_rate={:.6} last_fault_column={} of {}",193            label,194            radius,195            level,196            columns,197            column_hits,198            boxes,199            faults,200            faults as f64 / boxes.max(1) as f64,201            chord_faults,202            chord_faults as f64 / boxes.max(1) as f64,203            shallow,204            shallow_faults,205            shallow_faults as f64 / shallow.max(1) as f64,206            worst,207            top_here208        );209    }210}211212// CENTRE SEAT213214fn centre(label: &str, radius: u64) {215    let deep = levels(radius);216    for level in 0..deep.saturating_sub(1) {217        let below = rung(radius, level);218        let here = rung(radius, level + 1);219        let top_here = here.high.len() - 2;220        let top_below = below.high.len() - 2;221        let mut boxes = 0u64;222        let mut seated = 0u64;223        for i in 0..=top_here {224            let low = here.high[i + 1];225            let high = here.high[i];226            for y in low..=high {227                boxes += 1;228                let column = 3 * i + 1;229                if column <= top_below230                    && below.high[column + 1] <= 3 * y + 1231                    && 3 * y + 1 <= below.high[column]232                {233                    seated += 1;234                }235            }236        }237        println!(238            "centre {} r={} level={} boxes={} centre_crossed={} rate={:.6} derived={:.6} ratio={:.6}",239            label,240            radius,241            level,242            boxes,243            seated,244            seated as f64 / boxes as f64,245            1.0 / 3.0,246            3.0 * seated as f64 / boxes as f64247        );248    }249}250251// STRAIGHT LINE LADDER252253fn ones(value: u64, depth: u32) -> u32 {254    let mut mask = 0u32;255    let mut left = value;256    for slot in 0..depth {257        if left % 3 == 1 {258            mask |= 1 << slot;259        }260        left /= 3;261    }262    mask263}264265fn ladder(name: &str, num: i128, den: i128, off: i128, depth: u32) {266    let width = 3i128.pow(depth);267    let height = |x: i128| -> i128 { (num * (width - x) + off) / den };268    let mut total = 0u64;269    let mut alive = 0u64;270    let mut hits = vec![0u64; depth as usize];271    for x in 0..width {272        let column = ones(x as u64, depth);273        let low = height(x + 1);274        let high = height(x);275        for z in low..=high {276            total += 1;277            let both = column & ones(z as u64, depth);278            if both == 0 {279                alive += 1;280            }281            let mut left = both;282            while left != 0 {283                let slot = left.trailing_zeros() as usize;284                hits[slot] += 1;285                left &= left - 1;286            }287        }288    }289    let mut product = 1.0f64;290    for slot in 0..depth as usize {291        product *= 1.0 - hits[slot] as f64 / total as f64;292    }293    let survival = alive as f64 / total as f64;294    let psi = survival / product;295    let independent = (8.0f64 / 9.0).powi(depth as i32);296    println!(297        "line {} slope={:.12} shift={:.6} depth={} cells={} alive={} survival={:.9} marginal={:.9} psi={:.9} logpsi={:.6} rate={:.9} eight_ninths={:.9}",298        name,299        num as f64 / den as f64,300        off as f64 / den as f64,301        depth,302        total,303        alive,304        survival,305        product,306        psi,307        psi.ln(),308        survival.powf(1.0 / depth as f64),309        independent.powf(1.0 / depth as f64)310    );311}312313// RESONANCE TRACKING314315fn coprime(a: i128, b: i128) -> bool {316    let mut x = a;317    let mut y = b;318    while y != 0 {319        let t = x % y;320        x = y;321        y = t;322    }323    x == 1324}325326fn steeper(u: i128, r: i128, p: i128, q: i128) -> bool {327    if p < 0 {328        return true;329    }330    u * u * (q * q + p * p) > p * p * r * r331}332333fn shallower(u: i128, r: i128, p: i128, q: i128) -> bool {334    u * u * (q * q + p * p) < p * p * r * r335}336337fn cap_levels(r: i128, b: i128) -> u32 {338    let mut count = 0u32;339    let mut scale = 1i128;340    while scale * b * b < 2 * r {341        count += 1;342        scale *= 3;343    }344    count345}346347fn floor_levels(r: i128, a: i128, b: i128) -> u32 {348    if a < 1 || a * b + 1 > b * b {349        return 0;350    }351    let mut count = 0u32;352    let mut scale = 1i128;353    while 8 * scale * scale * b * b * b * b < r * r {354        count += 1;355        scale *= 3;356    }357    count358}359360fn track(radius: u64, a: i128, b: i128) -> (u32, u64, [u64; 3]) {361    let r = radius as i128;362    let low = a * b - 1;363    let high = a * b + 1;364    let den = b * b;365    let deep = levels(radius);366    let mut run = 0u32;367    let mut boxes = 0u64;368    let mut seats = [0u64; 3];369    for level in 0..deep {370        let scale = 3i128.pow(level);371        let top = r / scale;372        let mut here = 0u64;373        let mut mine = [0u64; 3];374        for i in 0..=top {375            let left = scale * i;376            let right = scale * (i + 1);377            if 2 * right * right > r * r {378                break;379            }380            if steeper(left, r, low, den) && shallower(right, r, high, den) {381                here += 1;382                mine[(i % 3) as usize] += 1;383            }384        }385        if here > 0 {386            assert_eq!(run, level);387            run = level + 1;388            boxes = here;389            seats = mine;390        }391    }392    (run, boxes, seats)393}394395fn resonance(label: &str, radius: u64) -> u64 {396    let r = radius as i128;397    let deep = levels(radius);398    let mut rows = 0u64;399    let mut live = 0u64;400    let mut at_cap = 0u64;401    let mut at_floor = 0u64;402    let mut budget = vec![0i128; deep as usize];403    for b in 1..=30i128 {404        for a in 0..=b {405            if !coprime(a, b) {406                continue;407            }408            let cap = cap_levels(r, b);409            let low = floor_levels(r, a, b);410            let (run, boxes, seats) = track(radius, a, b);411            assert!(run <= cap);412            assert!(run >= low);413            rows += 1;414            if cap < deep {415                live += 1;416            }417            if run == cap {418                at_cap += 1;419            }420            if low > 0 && run == low {421                at_floor += 1;422            }423            for level in 0..run as usize {424                budget[level] = budget[level].max(b);425            }426            println!(427                "track {label} r={radius} slope={a}/{b} eps=1/{den} run={run} cap={cap} floor={low} boxes={boxes} seats_mod3={s0},{s1},{s2}",428                den = b * b,429                s0 = seats[0],430                s1 = seats[1],431                s2 = seats[2]432            );433        }434    }435    println!(436        "track {label} r={radius} levels={deep} slopes={rows} live_caps={live} run_equals_cap={at_cap} run_equals_floor={at_floor}"437    );438    for level in 0..deep as usize {439        let scale = 3i128.pow(level as u32);440        let seen = budget[level];441        assert!(seen * seen * scale < 2 * r);442        println!(443            "budget {label} r={radius} level={level} rank={rank} biggest_denominator={seen} denominator_square_times_scale={weight} twice_r={twice}",444            rank = deep as usize - 1 - level,445            weight = seen * seen * scale,446            twice = 2 * r447        );448    }449    rows + deep as u64 + 1450}451452fn secant(label: &str, radius: u64) -> u64 {453    let r = radius as i128;454    let square = r * r;455    let deep = levels(radius);456    let mut best = 0u32;457    let mut witness = 0i128;458    for level in 1..deep {459        let width = 3i128.pow(level);460        let step = (width - 1) / 2;461        let top = r / width;462        for i in 0..=top {463            let u = width * i;464            if u + 2 * step >= r {465                break;466            }467            let v0 = root_floor((square - u * u) as u64) as i128;468            let v1 = root_floor((square - (u + step) * (u + step)) as u64) as i128;469            let v2 = root_floor((square - (u + 2 * step) * (u + 2 * step)) as u64) as i128;470            if (v0 - 2 * v1 + v2).abs() <= 1 {471                best = level;472                witness = u;473                break;474            }475        }476    }477    let step = (3i128.pow(best) - 1) / 2;478    let next = (3i128.pow(best + 1) - 1) / 2;479    assert!(step * step <= 3 * r);480    assert!(next * next > 3 * r);481    println!(482        "secant {label} r={radius} levels={deep} line_depth={best} column={witness} spacing={step} spacing_square={sq} three_r={three} next_spacing={next} next_square={nsq}",483        sq = step * step,484        three = 3 * r,485        nsq = next * next486    );487    1488}489490fn blind(sigma_num: i128, sigma_den: i128, back: u32) -> u64 {491    let scale = 3i128.pow(back);492    assert_eq!(scale % sigma_den, 0);493    let lift = scale / sigma_den;494    let slow = sigma_num * lift;495    let fast = slow + 1;496    let mut differ = 0i128;497    let samples = 3 * scale;498    for t in 0..samples {499        let mut a = [0i64; 4];500        let mut b = [0i64; 4];501        for k in 0..4i128 {502            let one = 2 * t + 1 - 2 * k * slow;503            let two = 2 * t + 1 - 2 * k * fast;504            a[k as usize] = one.div_euclid(2 * scale) as i64;505            b[k as usize] = two.div_euclid(2 * scale) as i64;506        }507        if mask_of(a, 3) != mask_of(b, 3) {508            differ += 1;509        }510    }511    assert!(differ * scale <= 6 * samples);512    println!(513        "blind slope={sigma_num}/{sigma_den} eps=1/{scale} samples={samples} differ={differ} fraction={frac:.9} bound={bound:.9}",514        frac = differ as f64 / samples as f64,515        bound = 6.0 / scale as f64516    );517    1518}519520// BOX AGAINST BLOCK521522fn heights(radius: u64) -> Vec<i128> {523    let r = radius as i128;524    let square = r * r;525    let mut out = Vec::with_capacity(radius as usize + 2);526    for x in 0..=r + 1 {527        if x > r {528            out.push(0);529        } else {530            out.push(root_floor((square - x * x) as u64) as i128);531        }532    }533    out534}535536fn line_bracket(v: &[i128]) -> Option<((i128, i128), (i128, i128))> {537    let span = v.len() as i128;538    let mut low: Option<(i128, i128)> = None;539    let mut high: Option<(i128, i128)> = None;540    for a in 0..span {541        for b in 0..span {542            if a == b {543                continue;544            }545            let gap = b - a;546            let rise = v[b as usize] - v[a as usize] + 1;547            if gap > 0 {548                high = Some(match high {549                    None => (rise, gap),550                    Some(seen) => {551                        if rise * seen.1 < seen.0 * gap {552                            (rise, gap)553                        } else {554                            seen555                        }556                    }557                });558            } else {559                low = Some(match low {560                    None => (-rise, -gap),561                    Some(seen) => {562                        if -rise * seen.1 > seen.0 * -gap {563                            (-rise, -gap)564                        } else {565                            seen566                        }567                    }568                });569            }570        }571    }572    let (a, b) = (low?, high?);573    if a.0 * b.1 < b.0 * a.1 {574        Some((a, b))575    } else {576        None577    }578}579580fn finest(radius: i128, left: i128, right: i128) -> Option<(i128, i128)> {581    let mut best: Option<(i128, i128)> = None;582    for b in 1..=30i128 {583        for a in 0..=b {584            if !coprime(a, b) {585                continue;586            }587            if steeper(left, radius, a * b - 1, b * b) && shallower(right, radius, a * b + 1, b * b)588            {589                if best.is_none_or(|seen| b > seen.1) {590                    best = Some((a, b));591                }592            }593        }594    }595    best596}597598fn blocks(label: &str, radius: u64, level: u32) -> u64 {599    let r = radius as i128;600    let scale = 3i128.pow(level);601    let v = heights(radius);602    let side = (r / scale + 2) as usize;603    let mut extent = vec![(i128::MAX, i128::MIN, i128::MAX, i128::MIN); side * side];604    for x in 0..=r {605        let top = v[x as usize];606        let bottom = v[(x + 1) as usize];607        let column = x / scale;608        let mut row = bottom;609        while row <= top {610            let band = row / scale;611            let slot = &mut extent[column as usize * side + band as usize];612            slot.0 = slot.0.min(x);613            slot.1 = slot.1.max(x);614            slot.2 = slot.2.min(row);615            slot.3 = slot.3.max(row);616            row = (band + 1) * scale;617        }618    }619    let mut boxes = 0u64;620    let mut across = 0u64;621    let mut short = 0u64;622    let mut lines = 0u64;623    let mut rows = 0u64;624    for column in 0..side as i128 {625        for band in 0..side as i128 {626            let slot = extent[column as usize * side + band as usize];627            if slot.1 < slot.0 {628                continue;629            }630            boxes += 1;631            let wide = slot.1 - slot.0 + 1;632            let tall = slot.3 - slot.2 + 1;633            let step = (wide - 1) / 2;634            let lift = (tall - 1) / 2;635            if wide == scale {636                across += 1;637            }638            if step * step <= 3 * r && lift * lift <= 3 * r {639                short += 1;640            }641            let stop = (slot.1 + 1).min(r);642            let cut: Vec<i128> = (slot.0..=stop).map(|x| v[x as usize]).collect();643            let Some((low, high)) = line_bracket(&cut) else {644                continue;645            };646            lines += 1;647            assert!(step * step <= 3 * r);648            let window = finest(r, slot.0, slot.1);649            let (a, b) = window.unwrap();650            let reach = t_of(r, slot.1) - t_of(r, slot.0);651            println!(652                "boxline {label} r={radius} level={level} box={column},{band} columns={from}..{to} wide={wide} tall={tall} t_in=({sl:.7},{sh:.7}) span={reach:.9} block_floor={floor:.9} share={share:.6} finest_window={a}/{b}",653                from = slot.0,654                to = slot.1,655                sl = -(high.0 as f64) / high.1 as f64,656                sh = -(low.0 as f64) / low.1 as f64,657                floor = scale as f64 / r as f64,658                share = reach / (scale as f64 / r as f64)659            );660            rows += 1;661        }662    }663    println!(664        "boxes {label} r={radius} level={level} boxes={boxes} side_to_side={across} both_extents_short={short} line_boxes={lines} block_cap_level={cap}",665        cap = 3i128.pow(level)666    );667    rows + 1668}669670fn t_of(r: i128, u: i128) -> f64 {671    u as f64 / ((r * r - u * u) as f64).sqrt()672}673674// PASS675676pub fn derive() {677    let alphabet = line_alphabet(9);678    let masks: Vec<u16> = alphabet.iter().map(|entry| entry.0).collect();679    for span in [12i64, 15] {680        let other: Vec<u16> = line_alphabet(span).iter().map(|entry| entry.0).collect();681        assert_eq!(other, masks);682    }683    println!(684        "line_alphabet states={} list=[{}]",685        masks.len(),686        show(&masks)687    );688    for radius in [6560u64, 19682, 12345] {689        compare("carpet", radius);690    }691    for radius in [6560u64, 19682] {692        frozen("carpet", radius);693        centre("carpet", radius);694    }695    let unit: i128 = 300000000000;696    let slopes: [(&str, i128, i128); 10] = [697        ("half", unit / 2, 0),698        ("half_shifted", unit / 2, 123456789012),699        ("third", unit / 3, 0),700        ("third_shifted", unit / 3, 123456789012),701        ("third_again", unit / 3, 271828182845),702        ("seventh", unit / 7, 0),703        ("seventh_shifted", unit / 7, 123456789012),704        ("root2m1", 124264068712, 123456789012),705        ("golden", 185410196625, 123456789012),706        ("quarterpi", 235619449020, 123456789012),707    ];708    for (name, num, off) in slopes {709        for depth in 4..=12u32 {710            ladder(name, num, unit, off, depth);711        }712    }713    for step in 2..=9u32 {714        let drift = unit / 3i128.pow(step);715        for depth in 4..=12u32 {716            ladder(717                &format!("near_third_3^-{step}"),718                unit / 3 + drift,719                unit,720                123456789012,721                depth,722            );723        }724    }725    let mut rows = 0u64;726    for radius in [728u64, 2186, 6560, 19682] {727        rows += resonance("carpet", radius);728        rows += secant("carpet", radius);729    }730    for back in 2..=7u32 {731        rows += blind(1, 3, back);732    }733    rows += blocks("carpet", 19682, 6);734    println!("track carpet rows={rows}");735}