main.rs

1use num_bigint::BigUint;2use std::env;3use std::time::Instant;45// DIGITS67fn ok_small(mut n: u128, base: u128, top: u128) -> bool {8    while n > 0 {9        if n % base > top {10            return false;11        }12        n /= base;13    }14    true15}1617fn ok_big(n: &BigUint, base: u32, top: u8) -> bool {18    n.to_radix_le(base).iter().all(|&d| d <= top)19}2021fn pow_big(base: u32, e: usize) -> BigUint {22    BigUint::from(base).pow(e as u32)23}2425// SEARCH2627struct Thin {28    cap: BigUint,29    top5: u8,30    top7: u8,31    p3: Vec<BigUint>,32    p5: Vec<BigUint>,33    p7: Vec<BigUint>,34    cut5: Vec<usize>,35    cut7: Vec<usize>,36    hits: Vec<BigUint>,37    nodes: u64,38}3940impl Thin {41    fn new(cap: BigUint, top5: u8, top7: u8) -> Thin {42        let mut p3 = vec![BigUint::from(1u32)];43        while *p3.last().unwrap() < cap {44            let next = p3.last().unwrap() * 3u32;45            p3.push(next);46        }47        let l = p3.len() - 1;48        let one = BigUint::from(1u32);49        let rest: Vec<BigUint> = (0..=l).map(|k| (&p3[k] - &one) / 2u32).collect();50        let top = &rest[l];51        let mut p5 = vec![BigUint::from(1u32)];52        while p5.last().unwrap() <= top {53            let next = p5.last().unwrap() * 5u32;54            p5.push(next);55        }56        let mut p7 = vec![BigUint::from(1u32)];57        while p7.last().unwrap() <= top {58            let next = p7.last().unwrap() * 7u32;59            p7.push(next);60        }61        let cut5: Vec<usize> = rest62            .iter()63            .map(|r| p5.iter().position(|q| q > r).unwrap())64            .collect();65        let cut7: Vec<usize> = rest66            .iter()67            .map(|r| p7.iter().position(|q| q > r).unwrap())68            .collect();69        Thin {70            cap,71            top5,72            top7,73            p3,74            p5,75            p7,76            cut5,77            cut7,78            hits: Vec::new(),79            nodes: 0,80        }81    }8283    fn walk(&mut self, k: usize, v: BigUint) {84        self.nodes += 1;85        if k == 0 {86            if ok_big(&v, 5, self.top5) && ok_big(&v, 7, self.top7) {87                self.hits.push(v);88            }89            return;90        }91        let one = BigUint::from(1u32);92        let h5 = &v / &self.p5[self.cut5[k]];93        if !ok_big(&h5, 5, self.top5) && !ok_big(&(&h5 + &one), 5, self.top5) {94            return;95        }96        let h7 = &v / &self.p7[self.cut7[k]];97        if !ok_big(&h7, 7, self.top7) && !ok_big(&(&h7 + &one), 7, self.top7) {98            return;99        }100        let w = &v + &self.p3[k - 1];101        self.walk(k - 1, v);102        if w < self.cap {103            self.walk(k - 1, w);104        }105    }106107    fn run(mut self) -> (Vec<BigUint>, u64, usize) {108        let l = self.p3.len() - 1;109        self.walk(l, BigUint::from(0u32));110        (self.hits, self.nodes, l)111    }112}113114fn members(cap: BigUint, top5: u8, top7: u8) -> Vec<BigUint> {115    Thin::new(cap, top5, top7).run().0116}117118fn scan(cap: u128, top5: u128, top7: u128) -> Vec<BigUint> {119    let mut out = Vec::new();120    let mut n = 0u128;121    while n < cap {122        if ok_small(n, 3, 1) && ok_small(n, 5, top5) && ok_small(n, 7, top7) {123            out.push(BigUint::from(n));124        }125        n += 1;126    }127    out128}129130fn shown(v: &[BigUint]) -> Vec<String> {131    v.iter().map(|x| x.to_string()).collect()132}133134// OEIS CONTROL135136const A030979: [u128; 23] = [137    0,138    1,139    10,140    756,141    757,142    3160,143    3186,144    3187,145    3250,146    7560,147    7561,148    7651,149    20007,150    59548377,151    59548401,152    45773612811,153    45775397187,154    237617431723407,155    24991943420078301,156    24991943420078302,157    24991943420078307,158    24991943715007536,159    24991943715007537,160];161162fn a030979() -> Vec<BigUint> {163    A030979.iter().map(|&x| BigUint::from(x)).collect()164}165166// BUDGETS167168fn dim(base: f64, alphabet: f64) -> f64 {169    alphabet.ln() / base.ln()170}171172fn up(x: f64) -> f64 {173    (x * 1e6).ceil() / 1e6174}175176fn near(x: f64) -> f64 {177    (x * 1e6).round() / 1e6178}179180fn down(x: f64) -> f64 {181    (x * 1e6).floor() / 1e6182}183184fn budget() {185    let d3 = dim(3.0, 2.0);186    let d5 = dim(5.0, 3.0);187    let d7 = dim(7.0, 3.0);188    let d7w = dim(7.0, 4.0);189    println!("dim log_3 2 {:.6}", near(d3));190    println!("dim log_5 3 {:.6}", near(d5));191    println!("dim log_7 3 {:.6}", near(d7));192    println!("dim log_7 4 {:.6}", near(d7w));193    println!("pair 3 5 {:.6}", up(d3 + d5 - 1.0));194    println!("pair 3 7 {:.6}", up(d3 + d7 - 1.0));195    println!("pair 5 7 {:.6}", up(d5 + d7 - 1.0));196    println!("triple 3 5 7 {:.6}", up(d3 + d5 + d7 - 2.0));197    println!("triple 3 5 7 wide {:.6}", up(d3 + d5 + d7w - 2.0));198    println!("erdos 3 5 {:.6}", 1.0 / 2.0 + 2.0 / 4.0);199    println!("erdos 3 7 {:.6}", 1.0 / 2.0 + 2.0 / 6.0);200    println!("erdos 5 7 {:.6}", 2.0 / 4.0 + 2.0 / 6.0);201    println!("erdos 3 7 wide {:.6}", 1.0 / 2.0 + 3.0 / 6.0);202    println!("erdos 5 7 wide {:.6}", 2.0 / 4.0 + 3.0 / 6.0);203}204205// VERBS206207fn control() {208    let cap = BigUint::from(100_000_000u128);209    let thin = members(cap.clone(), 2, 2);210    let brute = scan(100_000_000u128, 2, 2);211    println!("thin below 10^8 pruned {:?}", shown(&thin));212    println!("thin agrees with the scan {}", thin == brute);213    let wide = members(cap, 2, 3);214    let wide_brute = scan(100_000_000u128, 2, 3);215    println!("wide below 10^8 pruned {:?}", shown(&wide));216    println!("wide agrees with the scan {}", wide == wide_brute);217    let far = BigUint::from(*A030979.last().unwrap() + 1);218    let wide_far = members(far, 2, 3);219    println!("wide terms below the last A030979 term {}", wide_far.len());220    println!("wide rebuilds A030979 {}", wide_far == a030979());221    println!(222        "thin sits inside wide {}",223        thin.iter().all(|x| wide.contains(x))224    );225}226227fn report(base: u32, exp: usize, top7: u8, name: &str) {228    let cap = pow_big(base, exp);229    let t0 = Instant::now();230    let shownice = cap.to_string();231    let (hits, nodes, l) = Thin::new(cap, 2, top7).run();232    if shownice.len() <= 40 {233        println!("{} height {}^{} = {}", name, base, exp, shownice);234    } else {235        println!(236            "{} height {}^{} of {} decimal digits",237            name,238            base,239            exp,240            shownice.len()241        );242    }243    println!("{} powers of three {}", name, l);244    println!("{} nodes {}", name, nodes);245    println!("{} count {}", name, hits.len());246    if hits.len() <= 16 {247        println!("{} members {:?}", name, shown(&hits));248    } else {249        println!("{} largest {}", name, hits.last().unwrap());250    }251    let e = (hits.len() as f64).ln() / (exp as f64 * (base as f64).ln());252    println!("{} effective exponent {:.6}", name, down(e));253    println!("{} seconds {:.2}", name, t0.elapsed().as_secs_f64());254}255256fn main() {257    let args: Vec<String> = env::args().collect();258    let verb = args.get(1).map(|s| s.as_str()).unwrap_or("budget");259    let n: usize = args.get(2).map(|s| s.parse().unwrap()).unwrap_or(0);260    let arg = move || n;261    match verb {262        "budget" => budget(),263        "control" => control(),264        "reach" => big(move || report(3, arg(), 2, "thin")),265        "seven" => big(move || report(7, arg(), 2, "thin")),266        "wide" => big(move || report(3, arg(), 3, "wide")),267        "ten" => big(move || report(10, arg(), 3, "wide")),268        _ => println!("verbs budget control reach seven wide ten"),269    }270}271272fn big<F: FnOnce() + Send + 'static>(f: F) {273    std::thread::Builder::new()274        .stack_size(1 << 30)275        .spawn(f)276        .unwrap()277        .join()278        .unwrap();279}280281#[cfg(test)]282mod tests {283    use super::*;284285    fn list(v: &[u128]) -> Vec<BigUint> {286        v.iter().map(|&x| BigUint::from(x)).collect()287    }288289    #[test]290    fn the_thin_set_below_seven_to_the_seventeen() {291        let cap = pow_big(7, 17);292        assert_eq!(cap.to_string(), "232630513987207");293        assert_eq!(members(cap, 2, 2), list(&[0, 1, 3186, 3187, 20007]));294    }295296    #[test]297    fn the_thin_set_gains_no_member_by_three_to_the_two_hundred() {298        assert_eq!(299            members(pow_big(3, 200), 2, 2),300            list(&[0, 1, 3186, 3187, 20007])301        );302    }303304    #[test]305    fn the_pruned_walk_agrees_with_a_direct_scan() {306        let cap = 3u128.pow(15);307        assert_eq!(members(BigUint::from(cap), 2, 2), scan(cap, 2, 2));308        assert_eq!(members(BigUint::from(cap), 2, 3), scan(cap, 2, 3));309    }310311    #[test]312    fn the_wide_variant_rebuilds_a030979() {313        let far = BigUint::from(*A030979.last().unwrap() + 1);314        assert_eq!(members(far, 2, 3), a030979());315    }316317    #[test]318    fn every_member_passes_all_three_digit_rules() {319        for n in members(pow_big(3, 120), 2, 2) {320            assert!(ok_big(&n, 3, 1));321            assert!(ok_big(&n, 5, 2));322            assert!(ok_big(&n, 7, 2));323        }324    }325326    #[test]327    fn the_pair_budgets_are_positive_and_the_triple_budget_is_not() {328        let d3 = dim(3.0, 2.0);329        let d5 = dim(5.0, 3.0);330        let d7 = dim(7.0, 3.0);331        assert_eq!(format!("{:.6}", up(d3 + d5 - 1.0)), "0.313536");332        assert_eq!(format!("{:.6}", up(d3 + d7 - 1.0)), "0.195505");333        assert_eq!(format!("{:.6}", up(d5 + d7 - 1.0)), "0.247182");334        assert_eq!(format!("{:.6}", up(d3 + d5 + d7 - 2.0)), "-0.121889");335        assert_eq!(336            format!("{:.6}", up(d3 + d5 + dim(7.0, 4.0) - 2.0)),337            "0.025951"338        );339    }340341    #[test]342    fn the_dimensions_print_to_six_places() {343        assert_eq!(format!("{:.6}", near(dim(3.0, 2.0))), "0.630930");344        assert_eq!(format!("{:.6}", near(dim(5.0, 3.0))), "0.682606");345        assert_eq!(format!("{:.6}", near(dim(7.0, 3.0))), "0.564575");346        assert_eq!(format!("{:.6}", near(dim(7.0, 4.0))), "0.712414");347    }348}