sumset.rs

1use super::design::elements;2use crate::core::error::{value_error, Result};34/// The deepest level a [`Sumset`] is built to: `S meet [0, 3^20]`, a bit array of `436` MB.5pub const DEEPEST: u32 = 20;67/// The widest level a [`Pair`] names in either base, so `3^k`, `4^m` and every difference string fit a signed 64-bit integer.8pub const WIDEST: u32 = 30;910const BLOCK: usize = 8;1112// THE SUMSET1314/// The sumset `S = A + B` of Erdos problem 125 up to `3^level`: `A` the integers whose base-3 digits are all `0` or `1`, `B` those whose base-4 digits are.15///16/// One bit per integer of `[0, 3^level]`: the members of `A` are set directly and each power `4^j <= 3^level` is folded in by one shift-or pass, `S |= S << 4^j`.17/// A running count every eight words makes `card(S meet [1, x])` one table read and at most eight popcounts.18pub struct Sumset {19    level: u32,20    top: u64,21    words: Vec<u64>,22    runs: Vec<u64>,23}2425impl Sumset {26    /// Builds `S meet [0, 3^level]`.27    ///28    /// ```29    /// let s = mrlyrs::num::sumset::Sumset::new(4).unwrap();30    /// assert_eq!(s.count(81), Some(79));31    /// ```32    ///33    /// # Errors34    ///35    /// Errs when the level is zero or past [`DEEPEST`].36    pub fn new(level: u32) -> Result<Sumset> {37        if !(1..=DEEPEST).contains(&level) {38            return value_error(format!("the level runs from 1 to {DEEPEST}, not {level}."));39        }40        let top = 3u64.pow(level);41        let mut words = vec![0u64; (top / 64 + 1) as usize];42        for a in elements(3, &[0, 1], level as usize)43            .into_iter()44            .chain([0, top])45        {46            words[(a / 64) as usize] |= 1 << (a % 64);47        }48        let mut power = 1u64;49        while power <= top {50            shift_or(&mut words, power);51            power *= 4;52        }53        let spare = top % 64;54        if spare != 63 {55            let last = words.len() - 1;56            words[last] &= (1u64 << (spare + 1)) - 1;57        }58        let mut runs = Vec::with_capacity(words.len() / BLOCK + 2);59        let mut total = 0u64;60        for chunk in words.chunks(BLOCK) {61            runs.push(total);62            total += chunk.iter().map(|w| u64::from(w.count_ones())).sum::<u64>();63        }64        runs.push(total);65        Ok(Sumset {66            level,67            top,68            words,69            runs,70        })71    }7273    /// The level the array was built to.74    pub fn level(&self) -> u32 {75        self.level76    }7778    /// The largest integer the array holds, `3^level`.79    pub fn top(&self) -> u64 {80        self.top81    }8283    /// Whether `x` is in `S`, or `None` past the top.84    pub fn contains(&self, x: u64) -> Option<bool> {85        (x <= self.top).then(|| self.words[(x / 64) as usize] >> (x % 64) & 1 == 1)86    }8788    /// Counts `card(S meet [1, x])`, or `None` past the top.89    pub fn count(&self, x: u64) -> Option<u64> {90        (x <= self.top).then(|| self.upto(x) - 1)91    }9293    /// Reads the density `D(x) = card(S meet [1, x])/x`, or `None` at zero and past the top.94    pub fn density(&self, x: u64) -> Option<f64> {95        (1..=self.top)96            .contains(&x)97            .then(|| (self.upto(x) - 1) as f64 / x as f64)98    }99100    /// Reads the share of members in each of `cells` equal runs of the integers `[low, high)`, the strip of `S` a page draws.101    ///102    /// Run `i` is `[low + floor(i n/cells), low + floor((i + 1) n/cells))` with `n = high - low`.103    ///104    /// # Errors105    ///106    /// Errs when the window is empty or runs past the top, or when `cells` is zero or more than the integers in the window.107    pub fn fills(&self, low: u64, high: u64, cells: usize) -> Result<Vec<f64>> {108        if low >= high || high > self.top + 1 {109            return value_error(format!(110                "the window [{low}, {high}) must be a nonempty part of [0, {}].",111                self.top112            ));113        }114        let span = high - low;115        if cells == 0 || cells as u64 > span {116            return value_error(format!(117                "a window of {span} integers splits into 1 to {span} cells, not {cells}."118            ));119        }120        let edge = |i: u64| low + (u128::from(i) * u128::from(span) / cells as u128) as u64;121        Ok((0..cells as u64)122            .map(|i| {123                let (a, b) = (edge(i), edge(i + 1));124                (self.below(b) - self.below(a)) as f64 / (b - a) as f64125            })126            .collect())127    }128129    /// Reads the least and the greatest `D(x)` over each window `[edges[i], edges[i + 1])`, `None` for an empty window.130    ///131    /// # Errors132    ///133    /// Errs when the edges decrease, start at zero or run past `top + 1`.134    pub fn extremes(&self, edges: &[u64]) -> Result<Vec<Option<(f64, f64)>>> {135        if edges.first() == Some(&0)136            || edges.last().is_some_and(|&e| e > self.top + 1)137            || edges.windows(2).any(|w| w[0] > w[1])138        {139            return value_error(format!(140                "the edges must climb inside [1, {}].",141                self.top + 1142            ));143        }144        Ok(edges145            .windows(2)146            .map(|w| {147                let (a, b) = (w[0], w[1]);148                if a == b {149                    return None;150                }151                let mut members = self.below(a);152                let (mut low, mut high) = (f64::INFINITY, f64::NEG_INFINITY);153                for x in a..b {154                    members += self.words[(x / 64) as usize] >> (x % 64) & 1;155                    let d = (members - 1) as f64 / x as f64;156                    low = low.min(d);157                    high = high.max(d);158                }159                Some((low, high))160            })161            .collect())162    }163164    fn upto(&self, x: u64) -> u64 {165        self.below(x + 1)166    }167168    fn below(&self, x: u64) -> u64 {169        let word = (x / 64) as usize;170        let block = word / BLOCK;171        let whole: u64 = self.words[block * BLOCK..word]172            .iter()173            .map(|w| u64::from(w.count_ones()))174            .sum();175        let part = match (x % 64, self.words.get(word)) {176            (0, _) | (_, None) => 0,177            (bits, Some(w)) => u64::from((w & ((1u64 << bits) - 1)).count_ones()),178        };179        self.runs[block] + whole + part180    }181}182183fn shift_or(words: &mut [u64], shift: u64) {184    let (jump, bits) = ((shift / 64) as usize, (shift % 64) as u32);185    for i in (jump..words.len()).rev() {186        let from = i - jump;187        let mut v = words[from] << bits;188        if bits > 0 && from > 0 {189            v |= words[from - 1] >> (64 - bits);190        }191        words[i] |= v;192    }193}194195// THE PAIRS196197/// A pair of levels: the base-3 level `A_k = A meet [0, 3^k)` against the base-4 level `B_m = B meet [0, 4^m)`, `k` the field `three` and `m` the field `four`.198#[derive(Clone, Copy, Debug, PartialEq, Eq)]199pub struct Pair {200    /// The base-3 level `k`.201    pub three: u32,202    /// The base-4 level `m`.203    pub four: u32,204}205206impl Pair {207    /// Names the pair `(k, m)`.208    ///209    /// # Errors210    ///211    /// Errs when either level is zero or past [`WIDEST`].212    pub fn new(three: u32, four: u32) -> Result<Pair> {213        if !(1..=WIDEST).contains(&three) || !(1..=WIDEST).contains(&four) {214            return value_error(format!(215                "both levels run from 1 to {WIDEST}, not ({three}, {four})."216            ));217        }218        Ok(Pair { three, four })219    }220221    fn powers(&self) -> (u64, u64) {222        (3u64.pow(self.three), 4u64.pow(self.four))223    }224225    /// The largest element `d(k, m) = (3^k - 1)/2 + (4^m - 1)/3` of `A_k + B_m`.226    pub fn largest(&self) -> u64 {227        let (p, q) = self.powers();228        (p - 1) / 2 + (q - 1) / 3229    }230231    /// The first and the last integer of the open interval `(d(k, m), min(3^k, 4^m))`, which `S` misses, or `None` when it holds none.232    ///233    /// A sum with `a < 3^k` and `b < 4^m` is at most `d(k, m)`, and a sum with `a >= 3^k` or `b >= 4^m` is at least `min(3^k, 4^m)`.234    pub fn gap(&self) -> Option<(u64, u64)> {235        let (p, q) = self.powers();236        let d = self.largest();237        (p.min(q) > d + 1).then(|| (d + 1, p.min(q) - 1))238    }239240    /// Whether the pair is clean, `3^k > d(k, m)` and `4^m > d(k, m)`, so that `S meet [0, d] = A_k + B_m`.241    pub fn clean(&self) -> bool {242        let (p, q) = self.powers();243        p.min(q) > self.largest()244    }245246    /// Whether the pair is a gap copy, `2 4^m < 3^k + 5`: `A_k + B_m` is then two disjoint translates of `A_(k-1) + B_m` and its energy is twice theirs.247    pub fn copy(&self) -> bool {248        let (p, q) = self.powers();249        2 * q < p + 5250    }251252    /// The scaling `tau = 4^m/3^k`.253    pub fn scale(&self) -> f64 {254        let (p, q) = self.powers();255        q as f64 / p as f64256    }257258    /// The additive energy `E(k, m) = sum_x r(x)^2`, `r(x)` the number of ways `x = a + b` with `a` in `A_k` and `b` in `B_m`.259    ///260    /// Summed as `sum_t 2^(z_3(t) + z_4(t))` over the `3^m` integers `t` with an `m`-digit base-4 string in `{-1, 0, 1}`, `z_4(t)` its zero digits and `z_3(t)` the zero digits of the `k`-digit balanced ternary expansion of `t`, a `t` past `(3^k - 1)/2` adding nothing: a difference of two digit strings from `{0, 1}` is a digit string from `{-1, 0, 1}`, which fixes `t` in either base, and each zero difference arises twice. The cost is `3^m` strings.261    pub fn energy(&self) -> u128 {262        let powers: Vec<i64> = (0..self.four).map(|i| 4i64.pow(i)).collect();263        strings(&powers, 0, 0, 0, self.three)264    }265266    /// The energy ratio `Q(k, m) = E(k, m) (d + 1)/4^(k+m)` of the energy [`Pair::energy`] returns, the energy against its flat value, at least `1`; `card(A_k + B_m) >= (d + 1)/Q` by Cauchy-Schwarz.267    pub fn ratio(&self, energy: u128) -> f64 {268        let flat = 4f64.powi((self.three + self.four) as i32);269        energy as f64 * (self.largest() + 1) as f64 / flat270    }271}272273fn strings(powers: &[i64], at: usize, t: i64, zeros: u32, three: u32) -> u128 {274    if at == powers.len() {275        return balanced_zeros(t, three).map_or(0, |z| 1u128 << (z + zeros));276    }277    (-1..=1)278        .map(|digit| {279            strings(280                powers,281                at + 1,282                t + digit * powers[at],283                zeros + u32::from(digit == 0),284                three,285            )286        })287        .sum()288}289290fn balanced_zeros(mut t: i64, digits: u32) -> Option<u32> {291    let mut zeros = 0;292    for _ in 0..digits {293        let r = t.rem_euclid(3);294        zeros += u32::from(r == 0);295        t = (t - r) / 3 + i64::from(r == 2);296    }297    (t == 0).then_some(zeros)298}299300/// Lists the pairs of the census: every `(k, m)` with `4^m` within a factor `3` of `3^k` and `d(k, m) <= 3^level`, by `d`.301///302/// # Errors303///304/// Errs when the level is zero or past [`WIDEST`].305pub fn pairs(level: u32) -> Result<Vec<Pair>> {306    if !(1..=WIDEST).contains(&level) {307        return value_error(format!("the level runs from 1 to {WIDEST}, not {level}."));308    }309    let top = 3u64.pow(level);310    let mut out = Vec::new();311    for three in 1..=level {312        for four in 1..=WIDEST {313            let pair = Pair { three, four };314            let (p, q) = pair.powers();315            if q >= 3 * p {316                break;317            }318            if p < 3 * q && pair.largest() <= top {319                out.push(pair);320            }321        }322    }323    out.sort_by_key(Pair::largest);324    Ok(out)325}326327#[cfg(test)]328mod tests {329    use super::*;330331    fn level(k: u32) -> Vec<u64> {332        let mut out = elements(3, &[0, 1], k as usize);333        out.push(0);334        out335    }336337    fn quarter(m: u32) -> Vec<u64> {338        let mut out = elements(4, &[0, 1], m as usize);339        out.push(0);340        out341    }342343    fn histogram(pair: Pair) -> u128 {344        let mut r = vec![0u128; pair.largest() as usize + 1];345        for a in level(pair.three) {346            for b in quarter(pair.four) {347                r[(a + b) as usize] += 1;348            }349        }350        r.iter().map(|v| v * v).sum()351    }352353    fn up(pair: Pair) -> u128 {354        let flat = 4u128.pow(pair.three + pair.four);355        (pair.energy() * u128::from(pair.largest() + 1) * 1_000_000).div_ceil(flat)356    }357358    #[test]359    fn the_array_is_the_double_loop() {360        let s = Sumset::new(8).unwrap();361        let mut sums = vec![false; s.top() as usize + 1];362        for a in level(9) {363            for b in quarter(7) {364                if a + b <= s.top() {365                    sums[(a + b) as usize] = true;366                }367            }368        }369        for (x, &member) in sums.iter().enumerate() {370            assert_eq!(s.contains(x as u64), Some(member), "x = {x}");371        }372        assert_eq!(373            s.count(s.top()),374            Some(sums.iter().filter(|&&v| v).count() as u64 - 1)375        );376        assert_eq!(s.contains(s.top() + 1), None);377    }378379    #[test]380    fn the_first_non_members_are_a367090() {381        let s = Sumset::new(6).unwrap();382        let missed: Vec<u64> = (0..=480)383            .filter(|&x| s.contains(x) == Some(false))384            .collect();385        let mut want = vec![62, 63, 143, 144];386        want.extend(207..=242);387        want.extend(463..=480);388        assert_eq!(missed, want);389    }390391    #[test]392    fn the_densities_at_the_powers_of_three_read_as_printed() {393        let s = Sumset::new(10).unwrap();394        let read: Vec<u64> = (4..=10)395            .map(|k| {396                let x = 3u64.pow(k);397                s.count(x).unwrap() * 1_000_000 / x398            })399            .collect();400        assert_eq!(401            read,402            [975308, 835390, 858710, 887517, 908855, 864959, 778472]403        );404    }405406    #[test]407    fn the_window_extremes_read_as_printed() {408        let s = Sumset::new(10).unwrap();409        let edges: Vec<u64> = (5..=10).map(|k| 3u64.pow(k)).collect();410        let six = |v: f64| (v * 1e6).floor() as u64;411        let read = s.extremes(&edges).unwrap();412        let lows: Vec<u64> = read.iter().map(|w| six(w.unwrap().0)).collect();413        let highs: Vec<u64> = read.iter().map(|w| six(w.unwrap().1)).collect();414        assert_eq!(lows, [835390, 852729, 887517, 858945, 778468]);415        assert_eq!(highs, [913419, 903768, 912038, 931596, 913781]);416    }417418    #[test]419    fn the_strip_is_the_share_of_members() {420        let s = Sumset::new(6).unwrap();421        let strip = s.fills(200, 250, 5).unwrap();422        assert_eq!(strip, [0.7, 0.0, 0.0, 0.0, 0.7]);423        assert!(s.fills(0, 731, 5).is_err() && s.fills(10, 12, 3).is_err());424    }425426    #[test]427    fn the_energy_is_the_representation_histogram() {428        for (k, m) in [(3, 2), (4, 3), (5, 4), (6, 5), (7, 5), (9, 7)] {429            let pair = Pair::new(k, m).unwrap();430            assert_eq!(pair.energy(), histogram(pair), "({k}, {m})");431        }432    }433434    #[test]435    fn the_energy_ratios_read_as_printed() {436        let read: Vec<u128> = pairs(22)437            .unwrap()438            .into_iter()439            .filter(|p| p.three >= 6)440            .take(4)441            .map(up)442            .collect();443        assert_eq!(read, [1467705, 1638125, 1664808, 1724517]);444    }445446    #[test]447    fn the_census_holds_twenty_seven_pairs() {448        let census: Vec<Pair> = pairs(22)449            .unwrap()450            .into_iter()451            .filter(|p| p.three >= 6)452            .collect();453        assert_eq!(census.len(), 27);454        assert!(!census.contains(&Pair::new(22, 18).unwrap()));455    }456457    #[test]458    fn a_gap_holds_no_member() {459        let s = Sumset::new(12).unwrap();460        let mut seen = 0;461        for pair in pairs(12).unwrap() {462            if let Some((a, b)) = pair.gap() {463                assert!(pair.clean());464                seen += 1;465                assert!((a..=b).all(|x| s.contains(x) == Some(false)));466                assert_eq!(s.contains(a - 1), Some(true));467            }468        }469        assert!(seen >= 4);470    }471472    #[test]473    fn a_gap_copy_doubles_the_energy() {474        for (k, m) in [(6, 4), (7, 5), (11, 8), (12, 9)] {475            let pair = Pair::new(k, m).unwrap();476            assert!(pair.copy());477            assert_eq!(pair.energy(), 2 * Pair::new(k - 1, m).unwrap().energy());478        }479    }480}