signed.rs

1use crate::menergy::{bits_of, box_row, boxes, column, exps, qpow, random_column, Bits, Rng};23// SIEVES45pub fn spf_table(lim: usize) -> Vec<u32> {6    let mut spf = vec![0u32; lim + 1];7    let mut i = 2usize;8    while i <= lim {9        if spf[i] == 0 {10            let mut j = i;11            while j <= lim {12                if spf[j] == 0 {13                    spf[j] = i as u32;14                }15                j += i;16            }17        }18        i += 1;19    }20    spf21}2223pub fn mobius(lim: usize, spf: &[u32]) -> Vec<i8> {24    let mut mu = vec![0i8; lim + 1];25    if lim >= 1 {26        mu[1] = 1;27    }28    for n in 2..=lim {29        let p = spf[n] as usize;30        let m = n / p;31        mu[n] = if m % p == 0 { 0 } else { -mu[m] };32    }33    mu34}3536fn liouville(lim: usize, spf: &[u32]) -> Vec<i8> {37    let mut lam = vec![0i8; lim + 1];38    if lim >= 1 {39        lam[1] = 1;40    }41    for n in 2..=lim {42        lam[n] = -lam[n / spf[n] as usize];43    }44    lam45}4647fn sign_draw(lim: usize, seed: u64) -> Vec<i8> {48    let mut rng = Rng::new(seed);49    (0..=lim)50        .map(|_| if rng.sign() > 0.0 { 1i8 } else { -1i8 })51        .collect()52}5354pub fn coefficients(lim: usize) -> (Vec<&'static str>, Vec<Vec<i8>>) {55    let spf = spf_table(lim);56    let mu = mobius(lim, &spf);57    let names = vec![58        "1", "mu", "lambda", "rand1", "rand2", "rand3", "sfr1", "sfr2", "sfr3",59    ];60    let masked = |seed: u64| -> Vec<i8> {61        let d = sign_draw(lim, seed);62        (0..=lim).map(|n| d[n] * mu[n] * mu[n]).collect()63    };64    let vals = vec![65        vec![1i8; lim + 1],66        mu.clone(),67        liouville(lim, &spf),68        sign_draw(lim, 0x51_6e_ed_01),69        sign_draw(lim, 0x51_6e_ed_02),70        sign_draw(lim, 0x51_6e_ed_03),71        masked(0x5f_6e_ed_01),72        masked(0x5f_6e_ed_02),73        masked(0x5f_6e_ed_03),74    ];75    (names, vals)76}7778// THE ROW7980pub struct SignedRow {81    pub m: u64,82    pub n: u64,83    pub r: u128,84    pub e: u128,85    pub num: Vec<i128>,86    pub diag: Vec<f64>,87    pub rms: Option<f64>,88    pub csloss: Vec<f64>,89}9091pub fn signed_row(bits: &Bits, mm: u64, nn: u64, g: u128, x: u128, coefs: &[Vec<i8>]) -> SignedRow {92    let nb = coefs.len();93    let mut cprime = vec![0i64; nn as usize];94    let mut qq = vec![0i128; nb];95    let mut cc = vec![0i128; nb];96    let mut u = vec![0i64; nb];97    let mut r: u128 = 0;98    let mut e: u128 = 0;99    let dl = g as f64 / x as f64;100    let bsum: Vec<f64> = coefs101        .iter()102        .map(|b| (nn..2 * nn).map(|l| b[l as usize] as f64).sum())103        .collect();104    let mut l1 = vec![0.0f64; nb];105    let mut sq = vec![0.0f64; nb];106    for m in mm..2 * mm {107        u.iter_mut().for_each(|v| *v = 0);108        let mut c = 0i64;109        for l in nn..2 * nn {110            if bits.get((m * l) as usize) {111                c += 1;112                cprime[(l - nn) as usize] += 1;113                for (j, b) in coefs.iter().enumerate() {114                    u[j] += b[l as usize] as i64;115                }116            }117        }118        r += c as u128;119        e += (c as u128) * (c as u128);120        for j in 0..nb {121            qq[j] += (u[j] as i128) * (u[j] as i128);122            cc[j] += u[j] as i128;123            let v = u[j] as f64 - dl * bsum[j];124            l1[j] += v.abs();125            sq[j] += v * v;126        }127    }128    let csloss: Vec<f64> = (0..nb)129        .map(|j| {130            let bd = (mm as f64 * sq[j]).sqrt();131            if bd > 0.0 {132                l1[j] / bd133            } else {134                0.0135            }136        })137        .collect();138    let xi = x as i128;139    let gi = g as i128;140    let mi = mm as i128;141    let delta = g as f64 / x as f64;142    let mut num = Vec::with_capacity(nb);143    let mut diag = Vec::with_capacity(nb);144    for (j, b) in coefs.iter().enumerate() {145        let mut bs = 0i128;146        let mut s2 = 0i128;147        let mut pp = 0i128;148        for l in nn..2 * nn {149            let v = b[l as usize] as i128;150            bs += v;151            s2 += v * v;152            pp += v * v * cprime[(l - nn) as usize] as i128;153        }154        num.push(155            xi * xi * (qq[j] - pp) - 2 * gi * xi * (bs * cc[j] - pp)156                + mi * gi * gi * (bs * bs - s2),157        );158        diag.push(159            pp as f64 * (1.0 - delta) * (1.0 - delta) + (mi * s2 - pp) as f64 * delta * delta,160        );161    }162    SignedRow {163        m: mm,164        n: nn,165        r,166        e,167        num,168        diag,169        rms: None,170        csloss,171    }172}173174const RMSCAP: u64 = 2048;175176pub fn pair_rms(bits: &Bits, mm: u64, nn: u64, g: u128, x: u128) -> Option<f64> {177    if mm > RMSCAP {178        return None;179    }180    let ms = mm as usize;181    let mut tt = vec![0u32; ms * ms];182    let mut cm = vec![0u32; ms];183    let mut cprime = vec![0u32; nn as usize];184    let mut hits: Vec<u32> = Vec::with_capacity(ms);185    for l in nn..2 * nn {186        hits.clear();187        for m in mm..2 * mm {188            if bits.get((m * l) as usize) {189                hits.push((m - mm) as u32);190            }191        }192        cprime[(l - nn) as usize] = hits.len() as u32;193        for &a in hits.iter() {194            cm[a as usize] += 1;195            for &b in hits.iter() {196                tt[a as usize * ms + b as usize] += 1;197            }198        }199    }200    let delta = g as f64 / x as f64;201    let nf = nn as f64;202    let mut all = 0.0f64;203    for a in 0..ms {204        for b in 0..ms {205            let gg =206                tt[a * ms + b] as f64 - delta * (cm[a] as f64 + cm[b] as f64) + nf * delta * delta;207            all += gg * gg;208        }209    }210    let mut dg = 0.0f64;211    for &c in cprime.iter() {212        let w = c as f64 * (1.0 - delta) * (1.0 - delta) + (mm as f64 - c as f64) * delta * delta;213        dg += w * w;214    }215    Some((2.0 * (all - dg)).max(0.0).sqrt())216}217218pub fn engineered(219    bits: &Bits,220    mm: u64,221    nn: u64,222    g: u128,223    x: u128,224    rounds: usize,225) -> (f64, f64, usize) {226    let delta = g as f64 / x as f64;227    let ms = mm as usize;228    let ns = nn as usize;229    let mut b = vec![1.0f64; ns];230    let mut v = vec![0.0f64; ms];231    let mut cprime = vec![0u32; ns];232    for (i, m) in (mm..2 * mm).enumerate() {233        let mut acc = 0.0f64;234        for (j, l) in (nn..2 * nn).enumerate() {235            if bits.get((m * l) as usize) {236                acc += 1.0 - delta;237                cprime[j] += 1;238            } else {239                acc -= delta;240            }241        }242        v[i] = acc;243    }244    let w: Vec<f64> = cprime245        .iter()246        .map(|&c| c as f64 * (1.0 - delta) * (1.0 - delta) + (mm as f64 - c as f64) * delta * delta)247        .collect();248    let mut flips = 0usize;249    let mut grad = vec![0.0f64; ns];250    for _ in 0..rounds {251        grad.iter_mut().for_each(|z| *z = 0.0);252        for (i, m) in (mm..2 * mm).enumerate() {253            for (j, l) in (nn..2 * nn).enumerate() {254                let ps = if bits.get((m * l) as usize) {255                    1.0 - delta256                } else {257                    -delta258                };259                grad[j] += v[i] * ps;260            }261        }262        let mut best = 0usize;263        let mut gain = 0.0f64;264        for j in 0..ns {265            let d = -4.0 * b[j] * grad[j] + 4.0 * w[j];266            if d < gain {267                gain = d;268                best = j;269            }270        }271        if gain >= 0.0 {272            break;273        }274        let l = nn + best as u64;275        for (i, m) in (mm..2 * mm).enumerate() {276            let ps = if bits.get((m * l) as usize) {277                1.0 - delta278            } else {279                -delta280            };281            v[i] -= 2.0 * b[best] * ps;282        }283        b[best] = -b[best];284        flips += 1;285    }286    let form: f64 = v.iter().map(|z| z * z).sum();287    let dg: f64 = w.iter().sum();288    (form, dg, flips)289}290291fn round_div(num: i128, den: i128) -> i128 {292    if num >= 0 {293        (num + den / 2) / den294    } else {295        -((-num + den / 2) / den)296    }297}298299// THE SWEEP300301pub struct SignedSweep {302    pub top: SignedRow,303    pub l5: Option<SignedRow>,304    pub seen: usize,305    pub live: usize,306}307308pub fn sweep(bits: &Bits, x: u128, g: u128, coefs: &[Vec<i8>]) -> Option<SignedSweep> {309    let all = boxes(x);310    let edge = (x as f64).powf(0.4);311    let mut top: Option<SignedRow> = None;312    let mut l5: Option<SignedRow> = None;313    let mut live = 0usize;314    for &(mm, nn) in all.iter() {315        let row = signed_row(bits, mm, nn, g, x, coefs);316        if row.r == 0 {317            continue;318        }319        live += 1;320        let score = row.num[0].abs();321        if top.as_ref().map(|b| score > b.num[0].abs()).unwrap_or(true) {322            top = Some(signed_row(bits, mm, nn, g, x, coefs));323        }324        if mm as f64 >= edge325            && nn as f64 >= edge326            && l5.as_ref().map(|b| score > b.num[0].abs()).unwrap_or(true)327        {328            l5 = Some(row);329        }330    }331    top.map(|t| SignedSweep {332        top: t,333        l5,334        seen: all.len(),335        live,336    })337}338339// STUDY340341pub struct Cell {342    pub q: u64,343    pub digits: Vec<u64>,344    pub label: &'static str,345    pub depths: [usize; 2],346}347348fn ex(q: u64, e: u64) -> Vec<u64> {349    (0..q).filter(|&f| f != e).collect()350}351352pub fn cells() -> Vec<Cell> {353    vec![354        Cell {355            q: 3,356            digits: vec![0, 1],357            label: "01",358            depths: [12, 14],359        },360        Cell {361            q: 4,362            digits: vec![0, 1, 2],363            label: "012",364            depths: [8, 9],365        },366        Cell {367            q: 5,368            digits: vec![0, 1, 2, 3],369            label: "0123",370            depths: [7, 8],371        },372        Cell {373            q: 10,374            digits: ex(10, 7),375            label: "ex7",376            depths: [5, 6],377        },378    ]379}380381const PAIRCAP: usize = 8_000_000;382383pub struct Tally {384    pub cells: usize,385    pub rand: [usize; 3],386    pub sf: [usize; 3],387    pub mu_under_lam: usize,388    pub musq: (f64, f64),389    pub lamsq: (f64, f64),390    pub under_r: usize,391    pub one_over_r: usize,392    pub fd: (f64, f64),393    pub murms: (f64, f64),394    pub sfrms: (f64, f64),395    pub cs: (f64, f64),396}397398fn place(v: f64, lo: f64, hi: f64) -> usize {399    if v < lo {400        0401    } else if v > hi {402        2403    } else {404        1405    }406}407408fn emit(409    q: u64,410    label: &str,411    l: usize,412    x: u128,413    kind: &str,414    tag: &str,415    sw: &SignedSweep,416    row: &SignedRow,417    tally: &mut Tally,418) {419    let den = (x as i128) * (x as i128);420    let sig: Vec<i128> = row.num.iter().map(|&n| round_div(n, den)).collect();421    let val = |n: i128| (n as f64) / (den as f64);422    let base = row.num[0].abs() as f64;423    let ratio = |n: i128| {424        if base == 0.0 {425            0.0426        } else {427            (n.abs() as f64) / base428        }429    };430    let band = |a: usize, b: usize, c: usize| -> (f64, f64) {431        let v = [ratio(row.num[a]), ratio(row.num[b]), ratio(row.num[c])];432        (433            v.iter().cloned().fold(f64::INFINITY, f64::min),434            v.iter().cloned().fold(f64::NEG_INFINITY, f64::max),435        )436    };437    let (rlo, rhi) = band(3, 4, 5);438    let (slo, shi) = band(6, 7, 8);439    let mut fdlo = f64::INFINITY;440    let mut fdhi = f64::NEG_INFINITY;441    for j in 1..row.num.len() {442        let f = (val(row.num[j]) + row.diag[j]) / row.diag[j];443        fdlo = fdlo.min(f);444        fdhi = fdhi.max(f);445    }446    let (rtxt, mrms, lrms) = match row.rms {447        Some(v) if v > 0.0 => (448            format!("{:.4e}", v),449            format!("{:.4}", val(row.num[1]).abs() / v),450            format!("{:.4}", val(row.num[2]).abs() / v),451        ),452        _ => ("-".to_string(), "-".to_string(), "-".to_string()),453    };454    let root = val(row.num[0]).abs().sqrt();455    let sq = |n: i128| {456        if root == 0.0 {457            0.0458        } else {459            val(n).abs() / root460        }461    };462    let sqband = [sq(row.num[6]), sq(row.num[7]), sq(row.num[8])];463    let sqlo = sqband.iter().cloned().fold(f64::INFINITY, f64::min);464    let sqhi = sqband.iter().cloned().fold(f64::NEG_INFINITY, f64::max);465    if kind == "digit" {466        tally.cells += 1;467        tally.rand[place(ratio(row.num[1]), rlo, rhi)] += 1;468        tally.sf[place(ratio(row.num[1]), slo, shi)] += 1;469        if row.num[1].abs() < row.num[2].abs() {470            tally.mu_under_lam += 1;471        }472        tally.musq.0 = tally.musq.0.min(sq(row.num[1]));473        tally.musq.1 = tally.musq.1.max(sq(row.num[1]));474        tally.lamsq.0 = tally.lamsq.0.min(sq(row.num[2]));475        tally.lamsq.1 = tally.lamsq.1.max(sq(row.num[2]));476        let rr = row.r as f64;477        if val(row.num[1]).abs() < rr && val(row.num[2]).abs() < rr {478            tally.under_r += 1;479        }480        if val(row.num[0]).abs() > rr {481            tally.one_over_r += 1;482        }483        for j in 1..row.num.len() {484            let f = (val(row.num[j]) + row.diag[j]) / row.diag[j];485            tally.fd.0 = tally.fd.0.min(f);486            tally.fd.1 = tally.fd.1.max(f);487        }488        tally.cs.0 = tally.cs.0.min(row.csloss[1]);489        tally.cs.1 = tally.cs.1.max(row.csloss[1]);490        if let Some(rv) = row.rms {491            if rv > 0.0 {492                tally.murms.0 = tally.murms.0.min(val(row.num[1]).abs() / rv);493                tally.murms.1 = tally.murms.1.max(val(row.num[1]).abs() / rv);494                for j in 6..9 {495                    tally.sfrms.0 = tally.sfrms.0.min(val(row.num[j]).abs() / rv);496                    tally.sfrms.1 = tally.sfrms.1.max(val(row.num[j]).abs() / rv);497                }498            }499        }500    }501    println!(502        "| {} | {} | {} | {} | {} | {} | {}/{} | {} | {} | {} | {} | {} | {} | {:.6} | {:.6} | {:.6} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {} | {} | {} | {:.4} | {:.4} |",503        q,504        label,505        l,506        x,507        kind,508        tag,509        sw.live,510        sw.seen,511        row.m,512        row.n,513        row.r,514        sig[0],515        sig[1],516        sig[2],517        exps(val(row.num[0]).abs(), x),518        exps(val(row.num[1]).abs(), x),519        exps(val(row.num[2]).abs(), x),520        ratio(row.num[1]),521        ratio(row.num[2]),522        rlo,523        rhi,524        slo,525        shi,526        sq(row.num[1]),527        sq(row.num[2]),528        sqlo,529        sqhi,530        rtxt,531        mrms,532        lrms,533        fdlo,534        fdhi535    );536}537538pub fn run() {539    let cs = cells();540    let mut lim = 16usize;541    for c in cs.iter() {542        for &l in c.depths.iter() {543            let x = qpow(c.q, l);544            for &(_, nn) in boxes(x).iter() {545                lim = lim.max(2 * nn as usize);546            }547        }548    }549    let (_, coefs) = coefficients(lim);550    println!("menergy signed");551    println!("| q | F | L | x | kind | box | live/seen | M | N | R | Sigma_1 | Sigma_mu | Sigma_lam | exp_1 | exp_mu | exp_lam | mu/1 | lam/1 | rand lo | rand hi | sf lo | sf hi | mu/sq | lam/sq | sf sq lo | sf sq hi | rms | mu/rms | lam/rms | fd lo | fd hi |");552    let mut tally = Tally {553        cells: 0,554        rand: [0; 3],555        sf: [0; 3],556        mu_under_lam: 0,557        musq: (f64::INFINITY, f64::NEG_INFINITY),558        lamsq: (f64::INFINITY, f64::NEG_INFINITY),559        under_r: 0,560        one_over_r: 0,561        fd: (f64::INFINITY, f64::NEG_INFINITY),562        murms: (f64::INFINITY, f64::NEG_INFINITY),563        sfrms: (f64::INFINITY, f64::NEG_INFINITY),564        cs: (f64::INFINITY, f64::NEG_INFINITY),565    };566    for c in cs.iter() {567        for (d, &l) in c.depths.iter().enumerate() {568            let x = qpow(c.q, l);569            let vals = column(c.q, &c.digits, l);570            let g = vals.len() as u128;571            let bits = bits_of(&vals, x);572            let mut sw = match sweep(&bits, x, g, &coefs) {573                Some(s) => s,574                None => continue,575            };576            sw.top.rms = pair_rms(&bits, sw.top.m, sw.top.n, g, x);577            if let Some(r) = sw.l5.as_mut() {578                r.rms = pair_rms(&bits, r.m, r.n, g, x);579            }580            let delta = g as f64 / x as f64;581            for (tag, row) in [Some(("top", &sw.top)), sw.l5.as_ref().map(|r| ("l5", r))]582                .into_iter()583                .flatten()584            {585                if let Some(chk) = box_row(&bits, row.m, row.n, delta, PAIRCAP) {586                    assert_eq!(chk.r, row.r, "R disagrees with the census routine");587                    assert_eq!(chk.e, row.e, "E_x(M,N) disagrees with the census routine");588                }589                emit(c.q, c.label, l, x, "digit", tag, &sw, row, &mut tally);590            }591            if d + 1 == c.depths.len() {592                let mut rng = Rng::new(0x51_6e_ed_00 + c.q * 131 + l as u64);593                let rv = random_column(x, vals.len(), &mut rng);594                let rbits = bits_of(&rv, x);595                if let Some(mut rs) = sweep(&rbits, x, g, &coefs) {596                    rs.top.rms = pair_rms(&rbits, rs.top.m, rs.top.n, g, x);597                    if let Some(r) = rs.l5.as_mut() {598                        r.rms = pair_rms(&rbits, r.m, r.n, g, x);599                    }600                    let tag = if rs.l5.is_some() { "l5" } else { "top" };601                    let row = rs.l5.as_ref().unwrap_or(&rs.top);602                    emit(c.q, c.label, l, x, "random", tag, &rs, row, &mut tally);603                }604            }605        }606    }607    println!("menergy signed summary");608    println!("| digit cells | mu below rand | mu inside rand | mu above rand | mu below sf | mu inside sf | mu above sf | mu under lam | mu/sq lo | mu/sq hi | lam/sq lo | lam/sq hi | signed under R | unsigned over R | form/diag lo | form/diag hi | mu/rms lo | mu/rms hi | sf/rms lo | sf/rms hi | mu cs lo | mu cs hi |");609    println!(610        "| {} | {} | {} | {} | {} | {} | {} | {} | {:.4} | {:.4} | {:.4} | {:.4} | {} | {} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} |",611        tally.cells,612        tally.rand[0],613        tally.rand[1],614        tally.rand[2],615        tally.sf[0],616        tally.sf[1],617        tally.sf[2],618        tally.mu_under_lam,619        tally.musq.0,620        tally.musq.1,621        tally.lamsq.0,622        tally.lamsq.1,623        tally.under_r,624        tally.one_over_r,625        tally.fd.0,626        tally.fd.1,627        tally.murms.0,628        tally.murms.1,629        tally.sfrms.0,630        tally.sfrms.1,631        tally.cs.0,632        tally.cs.1633    );634    println!("menergy signed engineered");635    println!("| q | F | L | box | M | N | R | flips | form | diag | form/diag | bound | floor | bound/floor |");636    for c in cs.iter() {637        let l = c.depths[0];638        let x = qpow(c.q, l);639        let vals = column(c.q, &c.digits, l);640        let g = vals.len() as u128;641        let delta = g as f64 / x as f64;642        let bits = bits_of(&vals, x);643        let sw = match sweep(&bits, x, g, &coefs) {644            Some(s) => s,645            None => continue,646        };647        for (tag, row) in [Some(("top", &sw.top)), sw.l5.as_ref().map(|r| ("l5", r))]648            .into_iter()649            .flatten()650        {651            let (form, dg, flips) = engineered(&bits, row.m, row.n, g, x, 4000);652            let bound = (row.m as f64 * form).max(0.0).sqrt();653            let floor = (1.0 - delta) * (row.m as f64 * row.r as f64).sqrt();654            println!(655                "| {} | {} | {} | {} | {} | {} | {} | {} | {:.4e} | {:.4e} | {:.4} | {:.4e} | {:.4e} | {:.4} |",656                c.q, c.label, l, tag, row.m, row.n, row.r, flips, form, dg, form / dg, bound, floor, bound / floor657            );658        }659    }660}661662// TESTS663664#[cfg(test)]665mod tests {666    use super::*;667668    #[test]669    fn the_sieves_are_right() {670        let spf = spf_table(20);671        let mu = mobius(20, &spf);672        let lam = liouville(20, &spf);673        assert_eq!(&mu[1..=12], &[1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0][..]);674        assert_eq!(675            &lam[1..=12],676            &[1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1][..]677        );678    }679680    fn brute_num(bits: &Bits, mm: u64, nn: u64, g: u128, x: u128, b: &[i8]) -> i128 {681        let xi = x as i128;682        let gi = g as i128;683        let mut total = 0i128;684        for l1 in nn..2 * nn {685            for l2 in nn..2 * nn {686                if l1 == l2 {687                    continue;688                }689                let mut w = 0i128;690                for m in mm..2 * mm {691                    let p1 = if bits.get((m * l1) as usize) { xi } else { 0 } - gi;692                    let p2 = if bits.get((m * l2) as usize) { xi } else { 0 } - gi;693                    w += p1 * p2;694                }695                total += b[l1 as usize] as i128 * b[l2 as usize] as i128 * w;696            }697        }698        total699    }700701    #[test]702    fn signed_row_matches_the_definition() {703        let (_, coefs) = coefficients(256);704        for (q, digits, l, mm, nn) in [705            (3u64, vec![0u64, 1], 8usize, 8u64, 16u64),706            (3, vec![1, 2], 8, 16, 16),707            (4, vec![0, 1, 2], 6, 8, 32),708            (5, vec![0, 1, 2, 3], 5, 16, 32),709        ] {710            let x = qpow(q, l);711            let vals = column(q, &digits, l);712            let g = vals.len() as u128;713            let bits = bits_of(&vals, x);714            let row = signed_row(&bits, mm, nn, g, x, &coefs);715            for (j, b) in coefs.iter().enumerate() {716                assert_eq!(717                    row.num[j],718                    brute_num(&bits, mm, nn, g, x, b),719                    "coefficient {j} at q={q} L={l}"720                );721            }722        }723    }724725    #[test]726    fn the_unsigned_row_matches_the_census() {727        let (_, coefs) = coefficients(1024);728        for (q, digits, l, mm, nn) in [729            (3u64, vec![0u64, 1], 12usize, 256u64, 512u64),730            (4, vec![0, 1, 2], 9, 256u64, 256u64),731            (5, vec![0, 1, 2, 3], 8, 256u64, 256u64),732        ] {733            let x = qpow(q, l);734            let vals = column(q, &digits, l);735            let g = vals.len() as u128;736            let delta = g as f64 / x as f64;737            let bits = bits_of(&vals, x);738            let row = signed_row(&bits, mm, nn, g, x, &coefs);739            let chk = box_row(&bits, mm, nn, delta, PAIRCAP).unwrap();740            assert_eq!(chk.r, row.r);741            assert_eq!(chk.e, row.e);742        }743    }744745    #[test]746    fn the_random_column_kills_the_unsigned_sum() {747        let (_, coefs) = coefficients(1024);748        let (q, digits, l, mm, nn) = (3u64, vec![0u64, 1], 12usize, 128u64, 128u64);749        let x = qpow(q, l);750        let vals = column(q, &digits, l);751        let g = vals.len() as u128;752        let bits = bits_of(&vals, x);753        let digit = signed_row(&bits, mm, nn, g, x, &coefs);754        let mut rng = Rng::new(0x7e57);755        let rv = random_column(x, vals.len(), &mut rng);756        let rbits = bits_of(&rv, x);757        let rand = signed_row(&rbits, mm, nn, g, x, &coefs);758        assert!(759            digit.num[0] > 20 * rand.num[0].abs(),760            "digit {} against random {}",761            digit.num[0],762            rand.num[0]763        );764    }765}