main.rs

26.5 kB · rust · 969 lines

1use std::env;2use std::sync::atomic::{AtomicUsize, Ordering};3use std::sync::Mutex;4use std::time::Instant;56// SETS78fn gcd(a: u128, b: u128) -> u128 {9    if b == 0 {10        a11    } else {12        gcd(b, a % b)13    }14}1516fn sigma(bases: &[u64]) -> (u128, u128) {17    let den = bases.iter().fold(1u128, |l, &d| {18        let e = (d - 1) as u128;19        l / gcd(l, e) * e20    });21    let num = bases.iter().map(|&d| den / (d - 1) as u128).sum();22    (num, den)23}2425fn coprime(bases: &[u64]) -> bool {26    bases.iter().fold(0u128, |g, &d| gcd(g, d as u128)) == 127}2829fn family(bases: &[u64]) -> bool {30    let (num, den) = sigma(bases);31    !bases.is_empty() && num > den && coprime(bases)32}3334fn subsets(r: u64) -> Vec<Vec<u64>> {35    let pool: Vec<u64> = (3..=r).collect();36    (1u64..(1 << pool.len()))37        .map(|mask| {38            pool.iter()39                .enumerate()40                .filter(|(i, _)| mask >> i & 1 == 1)41                .map(|(_, &d)| d)42                .collect()43        })44        .collect()45}4647fn order(sets: &mut [Vec<u64>]) {48    sets.sort_by(|a, b| {49        a.iter()50            .max()51            .cmp(&b.iter().max())52            .then(a.len().cmp(&b.len()))53            .then(a.cmp(b))54    });55}5657fn minimal(r: u64) -> Vec<Vec<u64>> {58    let mut out: Vec<Vec<u64>> = subsets(r)59        .into_iter()60        .filter(|s| family(s))61        .filter(|s| {62            (0..s.len()).all(|i| {63                let mut t = s.clone();64                t.remove(i);65                !family(&t)66            })67        })68        .collect();69    order(&mut out);70    out71}7273fn unit(r: u64) -> Vec<Vec<u64>> {74    let mut out: Vec<Vec<u64>> = subsets(r)75        .into_iter()76        .filter(|s| {77            let (num, den) = sigma(s);78            num == den && coprime(s)79        })80        .collect();81    order(&mut out);82    out83}8485// ELEMENTS8687fn elements(bases: &[u64], k: u32, cap: u128) -> Vec<u128> {88    let mut v = vec![];89    for &d in bases {90        let mut x = (d as u128).pow(k);91        while x <= cap {92            v.push(x);93            x *= d as u128;94        }95    }96    v.sort();97    v98}99100// BITS101102struct Bits {103    w: Vec<u64>,104}105106impl Bits {107    fn new(bits: usize) -> Bits {108        let mut w = vec![0u64; bits / 64 + 2];109        w[0] = 1;110        Bits { w }111    }112113    fn fold(&mut self, a: usize, top: usize) {114        let q = a / 64;115        let r = (a % 64) as u32;116        for i in (q..=top / 64).rev() {117            let v = if r == 0 {118                self.w[i - q]119            } else {120                let lo = if i > q {121                    self.w[i - q - 1] >> (64 - r)122                } else {123                    0124                };125                (self.w[i - q] << r) | lo126            };127            self.w[i] |= v;128        }129    }130131    fn hole(&self, h: usize) -> Option<usize> {132        let mut i = h / 64;133        let keep = h % 64;134        let mut word = !self.w[i];135        if keep < 63 {136            word &= (1u64 << (keep + 1)) - 1;137        }138        loop {139            if word != 0 {140                return Some(i * 64 + 63 - word.leading_zeros() as usize);141            }142            if i == 0 {143                return None;144            }145            i -= 1;146            word = !self.w[i];147        }148    }149150    fn gaps(&self, below: usize) -> usize {151        (0..below)152            .filter(|&x| self.w[x / 64] >> (x % 64) & 1 == 0)153            .count()154    }155}156157// CERTIFY158159#[derive(Debug, Clone, PartialEq)]160enum Verdict {161    Found {162        f: u128,163        gaps: usize,164        n: usize,165        last: u128,166        sum: u128,167        reach: u128,168    },169    Cap,170}171172fn far() -> u128 {173    10u128.pow(30)174}175176fn grows(a: &[u128], n: usize, s: u128, t: u128, surplus: (u128, u128, u128)) -> Option<u128> {177    let (gain, den, c) = surplus;178    let mut sm = s;179    for &x in &a[n + 1..] {180        if gain * x >= c + den * (2 * t).saturating_sub(1) {181            return Some(x);182        }183        if x + 2 * t > sm + 1 {184            return None;185        }186        sm += x;187    }188    None189}190191fn certify(bases: &[u64], k: u32, cap: usize) -> Verdict {192    let (num, den) = sigma(bases);193    if num <= den {194        return Verdict::Cap;195    }196    let c: u128 = bases197        .iter()198        .map(|&d| (d as u128).pow(k) * (den / (d - 1) as u128))199        .sum();200    let a = elements(bases, k, far());201    let mut bits = Bits::new(cap + 64);202    let mut s: u128 = 0;203    for n in 0..a.len() - 1 {204        if s + a[n] > cap as u128 {205            return Verdict::Cap;206        }207        bits.fold(a[n] as usize, (s + a[n]) as usize);208        s += a[n];209        let t = bits.hole((s / 2) as usize).map_or(0, |x| x as u128 + 1);210        if t > a[n + 1] || t == 0 {211            continue;212        }213        if let Some(reach) = grows(&a, n, s, t, (num - den, den, c)) {214            return Verdict::Found {215                f: t - 1,216                gaps: bits.gaps(t as usize),217                n: n + 1,218                last: a[n],219                sum: s,220                reach,221            };222        }223    }224    Verdict::Cap225}226227// CONTROL228229fn knapsack(bases: &[u64], k: u32, b: usize) -> Vec<bool> {230    let mut r = vec![false; b + 1];231    r[0] = true;232    for a in elements(bases, k, b as u128) {233        let a = a as usize;234        for x in (a..=b).rev() {235            if r[x - a] {236                r[x] = true;237            }238        }239    }240    r241}242243fn control() {244    let cells: [(&[u64], u32); 9] = [245        (&[3, 4, 5], 1),246        (&[3, 4, 5], 2),247        (&[3, 4, 5], 3),248        (&[3, 4, 6], 2),249        (&[3, 5, 6, 7], 3),250        (&[4, 5, 6, 7, 8], 3),251        (&[3, 4, 7, 8], 3),252        (&[3, 5, 6, 9], 2),253        (&[3, 4, 9, 10], 2),254    ];255    for (bases, k) in cells {256        let clock = Instant::now();257        let v = certify(bases, k, 1 << 30);258        let Verdict::Found { f, gaps, .. } = v else {259            println!("{:?} k {} cap", bases, k);260            continue;261        };262        let b = (4 * f as usize).max(1000);263        let r = knapsack(bases, k, b);264        let last = (0..=b).rev().find(|&x| !r[x]).unwrap();265        let count = (0..=b).filter(|&x| !r[x]).count();266        println!(267            "{:?} k {} certificate F {} gaps {} | knapsack to {} F {} gaps {} | {} | {:.2} s",268            bases,269            k,270            f,271            gaps,272            b,273            last,274            count,275            if last as u128 == f && count == gaps {276                "agree"277            } else {278                "DISAGREE"279            },280            clock.elapsed().as_secs_f64()281        );282    }283    for bases in [&[3u64, 4][..], &[3, 6, 9, 12, 15, 21][..]] {284        println!("{:?} k 1 {:?} to 2^24", bases, certify(bases, 1, 1 << 24));285    }286}287288// CENSUS289290fn show(bases: &[u64]) -> String {291    format!(292        "{{{}}}",293        bases294            .iter()295            .map(|d| d.to_string())296            .collect::<Vec<_>>()297            .join(",")298    )299}300301fn census(r: u64, kmax: u32, bits: u32, threads: usize) {302    let clock = Instant::now();303    let sets = minimal(r);304    let next = AtomicUsize::new(0);305    let rows: Mutex<Vec<(usize, Vec<(u32, Verdict, f64)>)>> = Mutex::new(vec![]);306    std::thread::scope(|scope| {307        for _ in 0..threads {308            scope.spawn(|| loop {309                let i = next.fetch_add(1, Ordering::SeqCst);310                if i >= sets.len() {311                    break;312                }313                let mut row = vec![];314                for k in 1..=kmax {315                    let t = Instant::now();316                    let v = certify(&sets[i], k, 1usize << bits);317                    let stop = v == Verdict::Cap;318                    row.push((k, v, t.elapsed().as_secs_f64()));319                    if stop {320                        break;321                    }322                }323                rows.lock().unwrap().push((i, row));324            });325        }326    });327    let mut rows = rows.into_inner().unwrap();328    rows.sort_by_key(|x| x.0);329    let mut depth = vec![0usize; kmax as usize + 1];330    for (i, row) in &rows {331        let (num, den) = sigma(&sets[*i]);332        let mut line = format!(333            "{} sigma {}/{}",334            show(&sets[*i]),335            num / gcd(num, den),336            den / gcd(num, den)337        );338        let mut d = 0;339        for (k, v, t) in row {340            match v {341                Verdict::Found {342                    f,343                    gaps,344                    n,345                    last,346                    sum,347                    reach,348                } => {349                    d = *k;350                    line += &format!(351                        " | k {} F {} gaps {} n {} last {} sum {} reach {} {:.2}s",352                        k, f, gaps, n, last, sum, reach, t353                    );354                }355                Verdict::Cap => line += &format!(" | k {} cap 2^{}", k, bits),356            }357        }358        depth[d as usize] += 1;359        println!("{}", line);360    }361    println!("minimal sets {} below {}", sets.len(), r + 1);362    for (d, count) in depth.iter().enumerate() {363        if *count > 0 {364            println!("certified to k = {}: {} sets", d, count);365        }366    }367    let floor = rows368        .iter()369        .map(|(_, row)| row.iter().filter(|x| x.1 != Verdict::Cap).count())370        .min()371        .unwrap_or(0);372    println!(373        "every D in [3, {}] with sigma > 1 and gcd 1 is complete at k = 1..{}",374        r, floor375    );376    for s in unit(r) {377        println!("sigma = 1, gcd 1, outside the certificate: {}", show(&s));378    }379    println!("{:.1} s", clock.elapsed().as_secs_f64());380}381382fn cell(bases: &[u64], k: u32, bits: u32) {383    let clock = Instant::now();384    let v = certify(bases, k, 1usize << bits);385    println!(386        "{} k {} {:?} {:.1} s",387        show(bases),388        k,389        v,390        clock.elapsed().as_secs_f64()391    );392}393394fn sets(bases: &[u64], k: u32, b: usize) {395    let mut terms = elements(bases, k, b as u128);396    terms.dedup();397    let mut r = vec![false; b + 1];398    r[0] = true;399    for a in terms {400        let a = a as usize;401        for x in (a..=b).rev() {402            if r[x - a] {403                r[x] = true;404            }405        }406    }407    let last = (0..=b).rev().find(|&x| !r[x]).unwrap();408    println!(409        "{} k {} as a set: largest non-sum below {} is {}",410        show(bases),411        k,412        b,413        last414    );415}416417// WINDOWS418419fn windows(bases: &[u64], k: u32, lo: u128, hi: u128) {420    let a = elements(bases, k, far());421    let mut sums = vec![];422    let mut s: u128 = 0;423    for &x in &a {424        s += x;425        sums.push(s);426    }427    let open: Vec<usize> = (0..a.len() - 2)428        .filter(|&n| a[n + 1] >= lo && a[n + 1] <= hi && a[n + 2] > sums[n])429        .collect();430    let seen = (0..a.len() - 2)431        .filter(|&n| a[n + 1] >= lo && a[n + 1] <= hi)432        .count();433    let mut chains = 0;434    for &m in &open {435        for &n in &open {436            if n > m && sums[n].saturating_sub(a[n + 1]) < a[m + 2] && sums[m] < a[n + 2] - a[n + 1]437            {438                chains += 1;439            }440        }441    }442    println!(443        "{} k {}: {} of {} terms a_(n+1) in [{}, {}] open a window a_(n+2) > S_n; {} pairs m < n of open windows where (S_m, a_(m+2)) meets (S_n - a_(n+1), a_(n+2) - a_(n+1))",444        show(bases), k, open.len(), seen, lo, hi, chains445    );446}447448fn graham(t: u128, p: u128, q: u128, bits: u32) {449    let cap = 1usize << bits;450    let mut terms = vec![];451    let (mut num, mut den) = (t, 1u128);452    loop {453        num *= p;454        den *= q;455        let x = num / den;456        if x as usize > cap {457            break;458        }459        if x > 0 {460            terms.push(x as usize);461        }462    }463    terms.dedup();464    let mut bits_ = Bits::new(terms.iter().sum::<usize>() + 64);465    let mut s = 0usize;466    let mut sums = vec![];467    for &x in &terms {468        bits_.fold(x, s + x);469        s += x;470        sums.push(s);471    }472    let mut chained = 0;473    let mut windows = 0;474    for n in 0..terms.len().saturating_sub(2) {475        if terms[n + 2] > sums[n] && terms[n + 2] <= cap {476            windows += 1;477            if (sums[n] + 1..terms[n + 2]).any(|x| bits_.w[x / 64] >> (x % 64) & 1 == 0) {478                chained += 1;479            }480        }481    }482    let last = bits_483        .hole(cap / 2)484        .map_or("none".to_string(), |x| x.to_string());485    println!(486        "floor({} ({}/{})^n), n >= 1, {} terms to 2^{}: {} windows a_(n+2) > S_n, {} of them holding a non-sum; largest non-sum up to 2^{} is {}",487        t, p, q, terms.len(), bits, windows, chained, bits - 1, last488    );489}490491fn fibonacci(cap: usize) -> Vec<usize> {492    let mut v = vec![1usize, 3];493    while v[v.len() - 1] + v[v.len() - 2] < cap {494        let n = v.len();495        v.push(v[n - 1] + v[n - 2] + 1);496    }497    v498}499500fn missing(terms: &[usize]) -> Bits {501    let mut bits = Bits::new(terms.iter().sum::<usize>() + 64);502    let mut s = 0;503    for &x in terms {504        bits.fold(x, s + x);505        s += x;506    }507    bits508}509510fn refute(bits_: u32) {511    let terms = fibonacci(1 << bits_);512    let bits = missing(&terms);513    let mut s = 0i128;514    let mut least = i128::MAX;515    let (mut windows, mut held) = (0, 0);516    let mut sums = vec![];517    for &x in &terms {518        s += x as i128;519        sums.push(s as usize);520    }521    for n in 0..terms.len() - 1 {522        least = least.min(2 * (sums[n] as i128 - terms[n + 1] as i128) - terms[n + 1] as i128);523    }524    for n in 0..terms.len() - 2 {525        if terms[n + 2] > sums[n] {526            windows += 1;527            if (sums[n] + 1..terms[n + 2]).any(|x| bits.w[x / 64] >> (x % 64) & 1 == 0) {528                held += 1;529            }530        }531    }532    let low: Vec<usize> = (1..1_000_000)533        .filter(|&x| bits.w[x / 64] >> (x % 64) & 1 == 0)534        .collect();535    println!(536        "2 F_m - 1, m >= 2, {} terms to 2^{}: least 2 (S_n - a_(n+1)) - a_(n+1) = {}; {} windows, {} holding a non-sum; {} non-sums below 10^6: {:?}",537        terms.len(), bits_, least, windows, held, low.len(), low538    );539}540541// TERNARY542543const WORDS: usize = 128;544545#[derive(Clone, Copy)]546struct Sums {547    w: [u64; WORDS],548    top: usize,549}550551impl Sums {552    fn zero() -> Sums {553        let mut w = [0u64; WORDS];554        w[0] = 1;555        Sums { w, top: 0 }556    }557558    fn meets(&self, a: usize) -> bool {559        let q = a / 64;560        let r = (a % 64) as u32;561        for i in q..=(self.top + a) / 64 {562            let lo = if r > 0 && i > q {563                self.w[i - q - 1] >> (64 - r)564            } else {565                0566            };567            let v = if r == 0 {568                self.w[i - q]569            } else {570                (self.w[i - q] << r) | lo571            };572            if v & self.w[i] != 0 {573                return true;574            }575        }576        false577    }578579    fn with(&self, a: usize) -> Sums {580        let mut out = *self;581        for shift in [a, 2 * a] {582            let q = shift / 64;583            let r = (shift % 64) as u32;584            for i in q..=(self.top + shift) / 64 {585                let lo = if r > 0 && i > q {586                    self.w[i - q - 1] >> (64 - r)587                } else {588                    0589                };590                let v = if r == 0 {591                    self.w[i - q]592                } else {593                    (self.w[i - q] << r) | lo594                };595                out.w[i] |= v;596            }597        }598        out.top = self.top + 2 * a;599        out600    }601}602603struct Hunt<'a> {604    floor: &'a [usize],605    near: usize,606    need: u128,607    nodes: u64,608    found: Option<Vec<usize>>,609}610611impl Hunt<'_> {612    fn dfs(&mut self, sums: &Sums, chosen: &mut Vec<usize>, sq: u128, left: usize) {613        self.nodes += 1;614        if self.found.is_some() {615            return;616        }617        if left == 0 {618            self.found = Some(chosen.clone());619            return;620        }621        let top = *chosen.last().unwrap();622        let mut lo = self.floor[left - 1] + 1;623        if left >= 2 {624            lo = lo.max(self.near);625        }626        for a in (lo..top).rev() {627            let rest = (left - 1) as u128 * ((a - 1) as u128).pow(2);628            if sq + (a as u128).pow(2) + rest < self.need {629                break;630            }631            if sums.meets(a) || sums.meets(2 * a) {632                continue;633            }634            chosen.push(a);635            let next = sums.with(a);636            self.dfs(&next, chosen, sq + (a as u128).pow(2), left - 1);637            chosen.pop();638            if self.found.is_some() {639                return;640            }641        }642    }643}644645fn hunt(n: usize, top: usize, floor: &[usize], threads: usize) -> (Option<Vec<usize>>, u64) {646    near(n, top, floor, threads, 0)647}648649fn near(650    n: usize,651    top: usize,652    floor: &[usize],653    threads: usize,654    band: usize,655) -> (Option<Vec<usize>>, u64) {656    let need = (9u128.pow(n as u32) - 1).div_ceil(8);657    let root = Sums::zero().with(top);658    let lo = floor[n - 2] + 1;659    let seconds: Vec<usize> = (lo..top).rev().collect();660    let next = AtomicUsize::new(0);661    let out: Mutex<(Option<Vec<usize>>, u64)> = Mutex::new((None, 0));662    std::thread::scope(|scope| {663        for _ in 0..threads {664            scope.spawn(|| loop {665                let i = next.fetch_add(1, Ordering::SeqCst);666                if i >= seconds.len() || out.lock().unwrap().0.is_some() {667                    break;668                }669                let a = seconds[i];670                let sq = (top as u128).pow(2) + (a as u128).pow(2);671                if sq + (n as u128 - 2) * ((a - 1) as u128).pow(2) < need672                    || root.meets(a)673                    || root.meets(2 * a)674                {675                    continue;676                }677                if band > 0 && a + band < top {678                    continue;679                }680                let mut h = Hunt {681                    floor,682                    near: if band > 0 { top - band } else { 0 },683                    need,684                    nodes: 0,685                    found: None,686                };687                let mut chosen = vec![top, a];688                h.dfs(&root.with(a), &mut chosen, sq, n - 2);689                let mut o = out.lock().unwrap();690                o.1 += h.nodes;691                if h.found.is_some() && o.0.is_none() {692                    o.0 = h.found;693                }694            });695        }696    });697    out.into_inner().unwrap()698}699700fn admissible(a: &[usize]) -> bool {701    let n = a.len();702    let mut seen = std::collections::HashSet::new();703    for c in 0..3usize.pow(n as u32) {704        let (mut x, mut s) = (c, 0);705        for &v in a {706            s += (x % 3) * v;707            x /= 3;708        }709        if !seen.insert(s) {710            return false;711        }712    }713    true714}715716fn probe(n: usize, top: usize, threads: usize) {717    let floor = [0usize, 1, 3, 8, 22, 60, 168];718    let clock = Instant::now();719    let (found, nodes) = hunt(n, top, &floor[..n], threads);720    println!(721        "n {} top {}: {:?}, {} nodes, {:.1} s",722        n,723        top,724        found,725        nodes,726        clock.elapsed().as_secs_f64()727    );728}729730fn band(lo: usize, hi: usize, width: usize, threads: usize) {731    let floor = [0usize, 1, 3, 8, 22, 60, 168];732    let clock = Instant::now();733    for top in lo..=hi {734        let (found, nodes) = near(7, top, &floor, threads, width);735        if let Some(mut set) = found {736            set.reverse();737            assert!(admissible(&set));738            println!(739                "band {}: first top {} with six terms within {} of it: {:?}, {} nodes, {:.1} s",740                width,741                top,742                width,743                set,744                nodes,745                clock.elapsed().as_secs_f64()746            );747            return;748        }749    }750    println!(751        "band {}: none with top in [{}, {}], {:.1} s",752        width,753        lo,754        hi,755        clock.elapsed().as_secs_f64()756    );757}758759fn offsets(lo: usize, hi: usize) {760    let families: [[usize; 6]; 3] = [761        [0, 1, 3, 8, 22, 60],762        [0, 2, 6, 9, 23, 61],763        [0, 2, 5, 7, 21, 60],764    ];765    for fam in families {766        let mut best: Option<Vec<usize>> = None;767        'outer: for top in lo..=hi {768            for x in fam[5] + 1..top {769                let mut set: Vec<usize> = fam.iter().map(|o| top - o).collect();770                set.push(top - x);771                set.sort();772                if admissible(&set) {773                    best = Some(set);774                    break 'outer;775                }776            }777        }778        println!("offsets {:?} plus one: {:?}", fam, best);779    }780}781782fn ternary(nmax: usize, hi: usize, threads: usize) {783    let mut floor = vec![0usize, 1];784    println!("g_3(1) = 1");785    for n in 2..=nmax {786        let clock = Instant::now();787        let mut total = 0;788        let mut hit = None;789        for top in floor[n - 1] + 1..=hi {790            let (found, nodes) = hunt(n, top, &floor, threads);791            total += nodes;792            if let Some(set) = found {793                hit = Some((top, set));794                break;795            }796        }797        match hit {798            Some((top, mut set)) => {799                set.reverse();800                assert!(admissible(&set));801                println!(802                    "g_3({}) = {}, witness {:?}, {} nodes, {:.1} s",803                    n,804                    top,805                    set,806                    total,807                    clock.elapsed().as_secs_f64()808                );809                floor.push(top);810            }811            None => {812                println!(813                    "g_3({}) > {}, {} nodes, {:.1} s",814                    n,815                    hi,816                    total,817                    clock.elapsed().as_secs_f64()818                );819                return;820            }821        }822    }823}824825// ROUTE826827fn route(r: u64, lo: u128) {828    for bases in minimal(r) {829        let a = elements(&bases, 1, far());830        let mut below: u128 = 0;831        let mut best = (u128::MAX, 0u128);832        for &x in &a {833            let key = below * 1_000_000 / x;834            if x >= lo && key < best.0 {835                best = (key, x);836            }837            below += x;838        }839        println!(840            "{} least sum(M below N)/N over N in M, {} <= N <= 10^30: {}.{:06} at N = {}",841            show(&bases),842            lo,843            best.0 / 1_000_000,844            best.0 % 1_000_000,845            best.1846        );847    }848}849850fn main() {851    let args: Vec<String> = env::args().collect();852    let num = |i: usize, d: u64| args.get(i).map(|s| s.parse().unwrap()).unwrap_or(d);853    match args.get(1).map(|s| s.as_str()).unwrap_or("") {854        "census" => census(num(2, 10), num(3, 4) as u32, num(4, 31) as u32, num(5, 4) as usize),855        "control" => control(),856        "cell" => {857            let bases: Vec<u64> = args[2].split(',').map(|x| x.parse().unwrap()).collect();858            cell(&bases, num(3, 1) as u32, num(4, 31) as u32)859        }860        "set" => {861            let bases: Vec<u64> = args[2].split(',').map(|x| x.parse().unwrap()).collect();862            sets(&bases, num(3, 1) as u32, num(4, 20000) as usize)863        }864        "windows" => {865            let bases: Vec<u64> = args[2].split(',').map(|x| x.parse().unwrap()).collect();866            windows(&bases, num(3, 1) as u32, 10u128.pow(num(4, 6) as u32), 10u128.pow(num(5, 30) as u32))867        }868        "graham" => graham(num(2, 2) as u128, num(3, 5) as u128, num(4, 3) as u128, num(5, 22) as u32),869        "refute" => refute(num(2, 27) as u32),870        "ternary" => ternary(num(2, 6) as usize, num(3, 200) as usize, num(4, 8) as usize),871        "probe" => probe(num(2, 7) as usize, num(3, 420) as usize, num(4, 8) as usize),872        "offsets" => offsets(num(2, 420) as usize, num(3, 504) as usize),873        "band" => band(num(2, 419) as usize, num(3, 504) as usize, num(4, 70) as usize, num(5, 8) as usize),874        "route" => route(num(2, 10), 10u128.pow(num(3, 12) as u32)),875        _ => println!("verbs census R K BITS THREADS, cell D K BITS, set D K B, windows D K E F, graham T P Q BITS, refute BITS, ternary N HI THREADS, probe N TOP THREADS, band LO HI W THREADS, offsets LO HI, control, route R E"),876    }877}878879#[cfg(test)]880mod tests {881    use super::*;882883    #[test]884    fn the_unit_sets_below_eight_are_one() {885        assert_eq!(unit(8), vec![vec![3, 4, 7]]);886        assert_eq!(sigma(&[3, 4, 7]), (6, 6));887    }888889    #[test]890    fn the_minimal_sets_below_nine() {891        let want: Vec<Vec<u64>> = vec![892            vec![3, 4, 5],893            vec![3, 4, 6],894            vec![3, 4, 7, 8],895            vec![3, 5, 6, 7],896            vec![3, 5, 6, 8],897            vec![3, 5, 7, 8],898            vec![3, 6, 7, 8],899            vec![4, 5, 6, 7, 8],900        ];901        let mut got = minimal(8);902        got.sort();903        assert_eq!(got, want);904    }905906    #[test]907    fn the_bit_array_matches_a_plain_array() {908        let mut bits = Bits::new(4096);909        let mut plain = vec![false; 4096];910        plain[0] = true;911        let mut s = 0;912        for a in [3usize, 64, 5, 130, 77, 1000, 128] {913            bits.fold(a, s + a);914            for x in (a..=s + a).rev() {915                if plain[x - a] {916                    plain[x] = true;917                }918            }919            s += a;920            for h in [1usize, 63, 64, 65, s / 2, s] {921                let want = (0..=h).rev().find(|&x| !plain[x]);922                assert_eq!(bits.hole(h), want);923            }924        }925    }926927    #[test]928    fn three_four_five_at_k_one_misses_seventy_nine_last() {929        let Verdict::Found { f, .. } = certify(&[3, 4, 5], 1, 1 << 20) else {930            panic!()931        };932        assert_eq!(f, 79);933        let r = knapsack(&[3, 4, 5], 1, 5000);934        assert!(!r[79]);935        assert!((80..=5000).all(|x| r[x]));936    }937938    #[test]939    fn surplus_and_a_full_spectrum_leave_a_hole_chain() {940        let terms = fibonacci(1 << 20);941        assert_eq!(&terms[..8], &[1, 3, 5, 9, 15, 25, 41, 67]);942        let bits = missing(&terms);943        for x in [944            2usize, 7, 22, 63, 172, 459, 1212, 3185, 11, 36, 103, 280, 745, 1964, 5157,945        ] {946            assert_eq!(bits.w[x / 64] >> (x % 64) & 1, 0);947        }948        let mut s = 0i64;949        for n in 0..terms.len() - 1 {950            s += terms[n] as i64;951            assert!(2 * (s - terms[n + 1] as i64) - terms[n + 1] as i64 >= -9);952        }953    }954955    #[test]956    fn the_seven_set_under_four_hundred_seventy_five_is_admissible() {957        assert!(admissible(&[302, 409, 447, 459, 465, 466, 474]));958        assert!(!admissible(&[302, 409, 447, 459, 465, 467, 474]));959        let floor = [0usize, 1, 3, 8, 22];960        assert_eq!(hunt(5, 59, &floor, 2).0, None);961        assert!(hunt(5, 60, &floor, 2).0.is_some());962    }963964    #[test]965    fn a_set_below_the_threshold_never_certifies() {966        assert_eq!(certify(&[3, 4], 1, 1 << 22), Verdict::Cap);967        assert_eq!(certify(&[3, 6, 9, 12, 15, 21], 1, 1 << 22), Verdict::Cap);968    }969}