census.rs

1use mrlynum::factor::gcd;23pub struct Ray {4    pub a: u64,5    pub b: u64,6    pub weight: u64,7}89fn octave(a: u64, b: u64) -> usize {10    let top = a.max(b);11    if top == 1 {12        return 0;13    }14    let mut j = 1;15    while 3u64.pow(j as u32) <= top {16        j += 1;17    }18    j19}2021pub fn rays(level: u32) -> (u64, Vec<Ray>) {22    let size = 1usize << level;23    let value: Vec<u64> = (0..size)24        .map(|mask| {25            (0..level)26                .filter(|j| mask >> j & 1 == 1)27                .map(|j| 3u64.pow(j))28                .sum()29        })30        .collect();31    let mut keys: Vec<u64> = Vec::with_capacity(3usize.pow(level));32    let mut coprime = 0u64;33    for x in 0..size {34        let free = !x & (size - 1);35        let mut y = free;36        loop {37            let (a, b) = (value[x], value[y]);38            let g = gcd(a as usize, b as usize) as u64;39            if g == 1 {40                coprime += 1;41            }42            if g > 0 {43                keys.push((a / g) << 32 | (b / g));44            }45            if y == 0 {46                break;47            }48            y = (y - 1) & free;49        }50    }51    keys.sort_unstable();52    let mut out = Vec::new();53    let mut index = 0;54    while index < keys.len() {55        let key = keys[index];56        let mut end = index;57        while end < keys.len() && keys[end] == key {58            end += 1;59        }60        out.push(Ray {61            a: key >> 32,62            b: key & 0xffff_ffff,63            weight: (end - index) as u64,64        });65        index = end;66    }67    (coprime, out)68}6970pub fn report(level: u32, detail: bool) {71    let (coprime, rays) = rays(level);72    let points = 3f64.powi(level as i32);73    let mut bins = vec![0u64; level as usize + 1];74    let mut z = 0u128;75    let mut mass = 0u64;76    let mut peak = 0u64;77    let mut singles = 0u64;78    let mut light = 0u64;79    let mut light_z = 0u128;80    for ray in rays.iter() {81        bins[octave(ray.a, ray.b)] += 1;82        if ray.a == 0 || ray.b == 0 {83            continue;84        }85        z += (ray.weight as u128).pow(2);86        mass += ray.weight;87        peak = peak.max(ray.weight);88        if ray.weight == 1 {89            singles += 1;90        }91        if (2..=5).contains(&ray.weight) {92            light += 1;93            light_z += (ray.weight as u128).pow(2);94        }95    }96    println!(97        "census n {} gasket points {} primitive {}",98        level,99        3u64.pow(level),100        coprime101    );102    println!(103        "occupied rays {} with fibres {} without",104        rays.len(),105        rays.len() - 2106    );107    println!(108        "Z {} Z/3^n {:.6} sum M {} 3^n-2^(n+1)+1 {} max M {}",109        z,110        z as f64 / points,111        mass,112        3u64.pow(level) - 2u64.pow(level + 1) + 1,113        peak114    );115    let octaves: Vec<String> = bins116        .iter()117        .enumerate()118        .filter(|(_, c)| **c > 0)119        .map(|(j, c)| format!("{}:{}", j, c))120        .collect();121    println!("octaves {}", octaves.join(" "));122    println!("occ(8,n)/3^8 {:.3}", bins[8] as f64 / 6561.0);123    let cut = |j: usize| (bins[..=j].iter().sum::<u64>() as f64).ln() / 3f64.ln() / level as f64;124    if level % 2 == 0 {125        println!("band exponent j <= n/2 {:.4}", cut(level as usize / 2));126    } else {127        println!(128            "band exponent j <= floor(n/2) {:.4} or ceil {:.4}",129            cut(level as usize / 2),130            cut(level as usize / 2 + 1)131        );132    }133    if !detail {134        return;135    }136    let share = |part: u128| 100.0 * part as f64 / z as f64;137    println!(138        "M = 1 rays {} ({:.0}% of Z), M in [2,5] rays {} ({:.0}% of Z)",139        singles,140        share(singles as u128),141        light,142        share(light_z)143    );144    let mut heavy: Vec<&Ray> = rays.iter().filter(|r| r.a > 0 && r.b > 0).collect();145    heavy.sort_by(|p, q| q.weight.cmp(&p.weight).then(p.a.cmp(&q.a)));146    let top: Vec<String> = heavy[..10]147        .iter()148        .map(|r| format!("({},{}):{}", r.a, r.b, r.weight))149        .collect();150    let top_z: u128 = heavy[..10].iter().map(|r| (r.weight as u128).pow(2)).sum();151    println!(152        "ten heaviest {} carry {:.0}% of Z",153        top.join(" "),154        share(top_z)155    );156}