main.rs
40.5 kB · rust · 1418 lines
1use std::f64::consts::PI;23// SETUP45const GAMMA: f64 = 0.577_215_664_901_532_9;6const GUARD: f64 = 1e-12;7const SCI_GUARD: f64 = 1e-10;8const BASES: [u64; 10] = [9 1_000,10 2_000,11 3_000,12 3_689,13 3_690,14 5_000,15 10_000,16 100_000,17 1_000_000,18 1_000_000_000,19];20const SWEEP_LO: u64 = 3_690;21const SWEEP_HI: u64 = 100_000;22const CORO: [u64; 3] = [10_000, 100_000, 1_000_000];23const COST: [u64; 7] = [24 3_689,25 3_690,26 5_000,27 10_000,28 100_000,29 1_000_000,30 1_000_000_000,31];32const WALL: [u64; 2] = [3_689, 3_690];33const SHARP_HI: u64 = 4_000_000;34const SHARP_COST_HI: u64 = 100_000;3536// KERNEL CONSTANTS3738fn harmonic_bound(n: f64) -> f64 {39 n.ln() + GAMMA + 0.5 / n40}4142fn phi(q: f64) -> f64 {43 let n = ((q - 2.0) / 2.0).ceil();44 (4.0 / PI) * q + (2.0 * q / PI) * harmonic_bound(n) + (1.0 - 2.0 / PI) * (q - 2.0) + 0.72745}4647fn pb(q: f64, m: f64) -> f64 {48 m.sqrt() + phi(q) / q49}5051fn alpha(q: f64, m: f64) -> f64 {52 (q - m).ln() / q.ln()53}5455fn c_exp(q: f64, m: f64) -> f64 {56 pb(q, m).ln() / q.ln()57}5859fn delta(q: f64, m: f64) -> f64 {60 let a = alpha(q, m);61 (a - 0.75 - c_exp(q, m)) / a62}6364fn delta_stable(q: f64, m: f64) -> f64 {65 (gap(q, m) / pb(q, m)).ln_1p() / (q.ln() * alpha(q, m))66}6768fn defect(q: f64, m: f64) -> f64 {69 m / (2.0 * (q - m) * q.ln())70}7172fn psi(q: f64) -> f64 {73 let n = ((q - 2.0) / 2.0).ceil();74 (2.0 * q / PI) * harmonic_bound(n) + (1.0 - 2.0 / PI) * q75}7677fn psi_chord(q: f64) -> f64 {78 let p = (q / 2.0).floor();79 let h = harmonic_bound(p - 1.0);80 if (q as u64) % 2 == 0 {81 (q / PI) * (2.0 * h - 1.0 + 1.0 / p) + (1.0 - 2.0 / PI) * q / 2.082 } else {83 (q / PI) * (2.0 * h - 1.0 + 2.0 / p) + (1.0 - 2.0 / PI) * (q / 2.0 + 1.0 / (2.0 * q))84 }85}8687fn sec_off(q: f64, extreme: bool) -> f64 {88 let c = if extreme {89 (q - 1.0) / 2.090 } else if (q as u64) % 2 == 1 {91 0.092 } else {93 0.594 };95 1.0 / (PI * c / q).cos()96}9798#[derive(Clone, Copy)]99enum Cst {100 Step3,101 Sharp(bool),102 Chord(bool),103}104105fn onestep(q: f64, c: Cst) -> f64 {106 match c {107 Cst::Step3 => 1.0 + phi(q) / q,108 Cst::Sharp(e) => 4.0 / PI + psi(q) / q + 0.5 - sec_off(q, e) / (2.0 * q),109 Cst::Chord(e) => 4.0 / PI + psi_chord(q) / q + 0.5 - sec_off(q, e) / (2.0 * q),110 }111}112113fn onestep_floor(c: Cst) -> u64 {114 match c {115 Cst::Step3 => 3,116 Cst::Sharp(_) => 17,117 Cst::Chord(_) => 36,118 }119}120121fn onestep_delta(q: f64, c: Cst) -> f64 {122 let a = alpha(q, 1.0);123 (a - 0.75 - onestep(q, c).ln() / q.ln()) / a124}125126fn onestep_wall(c: Cst, hi: u64) -> (u64, bool) {127 let (mut first, mut held) = (0u64, 0u64);128 for q in onestep_floor(c)..=hi {129 let qf = q as f64;130 if (qf - 1.0) * qf.powf(-0.75) > onestep(qf, c) {131 held += 1;132 if first == 0 {133 first = q;134 }135 }136 }137 (first, first > 0 && held == hi - first + 1)138}139140fn onestep_cross(c: Cst, hi: u64) -> (u64, bool, f64, f64) {141 let (mut first, mut held) = (0u64, 0u64);142 for q in onestep_floor(c)..=hi {143 let qf = q as f64;144 if onestep_delta(qf, c) > defect(qf, 1.0) {145 held += 1;146 if first == 0 {147 first = q;148 }149 }150 }151 let qf = first as f64;152 (153 first,154 first > 0 && held == hi - first + 1,155 onestep_delta(qf, c),156 defect(qf, 1.0),157 )158}159160fn closes(q: f64, m: f64) -> bool {161 0.75 + c_exp(q, m) < alpha(q, m)162}163164fn gap_b(q: f64, b: f64, m: f64) -> f64 {165 (q - m) * q.powf(-b) - pb(q, m)166}167168fn gap(q: f64, m: f64) -> f64 {169 gap_b(q, 0.75, m)170}171172// SAFE ROUNDING173174fn ceil_units(x: f64, digits: u32) -> i64 {175 let s = 10f64.powi(digits as i32);176 ((x + GUARD) * s).ceil() as i64177}178179fn floor_units(x: f64, digits: u32) -> i64 {180 let s = 10f64.powi(digits as i32);181 ((x - GUARD) * s).floor() as i64182}183184fn fixed(units: i64, digits: u32) -> String {185 let s = 10i64.pow(digits);186 let sign = if units < 0 { "-" } else { "" };187 let a = units.abs();188 format!(189 "{}{}.{:0width$}",190 sign,191 a / s,192 a % s,193 width = digits as usize194 )195}196197fn sci(x: f64, decimals: u32, up: bool) -> String {198 let neg = x < 0.0;199 let a = x.abs();200 assert!(a > 0.0 && a.is_finite(), "sci needs a nonzero finite value");201 let mag_up = up != neg;202 let b = if mag_up {203 a * (1.0 + SCI_GUARD)204 } else {205 a * (1.0 - SCI_GUARD)206 };207 let mut e = b.log10().floor() as i32;208 let mut u = {209 let s = 10f64.powi(decimals as i32 - e);210 if mag_up {211 (b * s).ceil() as i64212 } else {213 (b * s).floor() as i64214 }215 };216 let lo = 10i64.pow(decimals);217 if u >= 10 * lo {218 u /= 10;219 e += 1;220 }221 if u < lo {222 u *= 10;223 e -= 1;224 }225 format!("{}{}e{}", if neg { "-" } else { "" }, fixed(u, decimals), e)226}227228fn label(q: u64) -> String {229 let mut p = q;230 let mut e = 0u32;231 while p % 10 == 0 {232 p /= 10;233 e += 1;234 }235 if p == 1 && e >= 4 {236 format!("10^{e}")237 } else {238 q.to_string()239 }240}241242// TABLE ROWS243244struct Row {245 q: u64,246 m: u64,247 alpha: i64,248 c: i64,249 delta: i64,250 closes: bool,251}252253fn row(q: u64, m: u64) -> Row {254 let (qf, mf) = (q as f64, m as f64);255 Row {256 q,257 m,258 alpha: floor_units(alpha(qf, mf), 6),259 c: ceil_units(c_exp(qf, mf), 5),260 delta: floor_units(delta(qf, mf), 5),261 closes: closes(qf, mf),262 }263}264265fn render(r: &Row) -> String {266 format!(267 "| {} | {} | {} | {} | {} |",268 label(r.q),269 fixed(r.alpha, 6),270 fixed(r.c, 5),271 fixed(r.delta, 5),272 if r.closes { "yes" } else { "no" }273 )274}275276fn render_cost(q: u64, m: u64) -> String {277 let (qf, mf) = (q as f64, m as f64);278 let d = delta_stable(qf, mf);279 let f = defect(qf, mf);280 format!(281 "| {} | {} | {} | {} | {} |",282 label(q),283 m,284 sci(d, 5, false),285 sci(f, 5, true),286 if d > f { "yes" } else { "no" }287 )288}289290fn render_sharp(name: &str, digits: &str, c: Cst) -> String {291 let (w, w_up) = onestep_wall(c, SHARP_HI);292 let (x, x_up, d, f) = onestep_cross(c, SHARP_COST_HI);293 format!(294 "| {} | {} | {} | {} | {} | {} | {} |",295 name,296 digits,297 label(w),298 label(x),299 sci(d, 5, false),300 sci(f, 5, true),301 if w_up && x_up { "yes" } else { "no" }302 )303}304305fn render_coro(r: &Row) -> String {306 format!(307 "| {} | {} | {} | {} | {} |",308 label(r.q),309 r.m,310 fixed(r.alpha, 6),311 fixed(r.c, 5),312 fixed(r.delta, 5)313 )314}315316// CHECKS317318fn gap_sweep(lo: u64, hi: u64, m: u64) -> (u64, f64, f64) {319 let mf = m as f64;320 let mut prev = gap(lo as f64, mf);321 let mut worst = f64::INFINITY;322 let mut arg = lo;323 let mut min_gap = prev;324 for q in (lo + 1)..=hi {325 let g = gap(q as f64, mf);326 let step = g - prev;327 if step < worst {328 worst = step;329 arg = q - 1;330 }331 if g < min_gap {332 min_gap = g;333 }334 prev = g;335 }336 (arg, worst, min_gap)337}338339fn crossover(lo: u64, hi: u64) -> (u64, u64) {340 let mut first = 0u64;341 let mut down = 0u64;342 let mut prev = f64::NEG_INFINITY;343 for q in lo..=hi {344 let qf = q as f64;345 let d = delta_stable(qf, 1.0) - defect(qf, 1.0);346 if first == 0 && d > 0.0 {347 first = q;348 }349 if q > lo && d <= prev {350 down += 1;351 }352 prev = d;353 }354 (first, down)355}356357fn max_excluded(q: u64) -> u64 {358 let qf = q as f64;359 let mut m = 0u64;360 while gap(qf, (m + 1) as f64) > 0.0 {361 m += 1;362 }363 m364}365366// LADDER367368const LADDER: [(i64, i64); 10] = [369 (1, 2),370 (13, 25),371 (11, 20),372 (4, 7),373 (3, 5),374 (2, 3),375 (3, 4),376 (4, 5),377 (9, 10),378 (19, 20),379];380const MONO_C: f64 = 1.291;381const MONO_LO: f64 = 40.0;382const EXACT_ULPS: f64 = 1024.0;383const DYADIC_TOP: f64 = 9_007_199_254_740_992.0;384const CORO_Q: f64 = 1e7;385386fn gcd(mut a: i64, mut b: i64) -> i64 {387 while b != 0 {388 let t = a % b;389 a = b;390 b = t;391 }392 a.abs()393}394395#[derive(Clone, Copy)]396struct Rat {397 n: i64,398 d: i64,399}400401impl Rat {402 fn new(n: i64, d: i64) -> Rat {403 let g = gcd(n, d);404 Rat { n: n / g, d: d / g }405 }406407 fn val(self) -> f64 {408 self.n as f64 / self.d as f64409 }410411 fn lt(self, o: Rat) -> bool {412 self.n * o.d < o.n * self.d413 }414415 fn same(self, o: Rat) -> bool {416 self.n * o.d == o.n * self.d417 }418419 fn show(self) -> String {420 if self.d == 1 {421 self.n.to_string()422 } else {423 format!("{}/{}", self.n, self.d)424 }425 }426}427428fn zhang(a: Rat) -> Option<Rat> {429 if a.lt(Rat::new(1, 2)) || Rat::new(4, 7).lt(a) {430 return None;431 }432 Some(Rat::new(433 a.n * (8 * a.d - 7 * a.n),434 2 * a.d * (2 * a.d - a.n),435 ))436}437438fn baker_harman(a: Rat) -> Rat {439 if a.lt(Rat::new(11, 20)) {440 Rat::new(4 * a.n + a.d, 4 * a.d)441 } else if a.lt(Rat::new(3, 5)) {442 Rat::new(4, 5)443 } else {444 Rat::new(a.n + a.d, 2 * a.d)445 }446}447448fn ladder_b(a: Rat) -> (Rat, &'static str) {449 let h = baker_harman(a);450 match zhang(a) {451 None => (h, "BH"),452 Some(z) if z.lt(h) => (z, "Zhang"),453 Some(z) if h.lt(z) => (h, "BH"),454 Some(z) => (z, "both"),455 }456}457458fn mono_ok(q: f64, b: f64) -> bool {459 (1.0 - b) * (q - 2.0) * (q + 1.0).powf(-b) >= MONO_C460}461462fn mono_floor(b: f64) -> f64 {463 let mut hi = MONO_LO;464 while !mono_ok(hi, b) {465 hi *= 2.0;466 }467 let mut lo = MONO_LO;468 for _ in 0..400 {469 let mid = lo + (hi - lo) / 2.0;470 if mid <= lo || mid >= hi {471 break;472 }473 if mono_ok(mid, b) {474 hi = mid;475 } else {476 lo = mid;477 }478 }479 hi.ceil()480}481482fn pb_low(q: f64) -> f64 {483 1.0 + 4.0 / PI484 + (2.0 / PI) * (((q - 2.0) / 2.0).ln() + GAMMA)485 + (1.0 - 2.0 / PI) * (q - 2.0) / q486}487488fn u_bound(q: f64, b: f64) -> f64 {489 q.powf(1.0 - b) - pb_low(q)490}491492fn ladder_wall(b: f64) -> f64 {493 let floor = mono_floor(b).max(SWEEP_LO as f64);494 if gap_b(floor, b, 1.0) > 0.0 {495 return floor;496 }497 let mut hi = floor;498 while gap_b(hi, b, 1.0) <= 0.0 {499 hi *= 2.0;500 }501 let mut lo = floor;502 if hi <= DYADIC_TOP {503 while hi - lo > 1.0 {504 let mid = (lo + (hi - lo) / 2.0).floor();505 if gap_b(mid, b, 1.0) > 0.0 {506 hi = mid;507 } else {508 lo = mid;509 }510 }511 return hi;512 }513 for _ in 0..400 {514 let mid = lo + (hi - lo) / 2.0;515 if mid <= lo || mid >= hi {516 break;517 }518 if gap_b(mid, b, 1.0) > 0.0 {519 hi = mid;520 } else {521 lo = mid;522 }523 }524 hi525}526527fn ladder_exact(q: f64, b: f64) -> bool {528 if q >= DYADIC_TOP {529 return false;530 }531 let noise = EXACT_ULPS * f64::EPSILON * ((q - 1.0) * q.powf(-b)).max(pb(q, 1.0));532 gap_b(q, b, 1.0) > noise && gap_b(q - 1.0, b, 1.0) < -noise533}534535fn integer_or_bound(q: f64, exact: bool) -> String {536 if exact {537 format!("{}", q as u64)538 } else {539 format!("<= {}", sci(q, 5, true))540 }541}542543fn max_excluded_b(q: f64, b: f64) -> u64 {544 let mut m = 0u64;545 while gap_b(q, b, (m + 1) as f64) > 0.0 {546 m += 1;547 }548 m549}550551fn render_ladder(a: Rat) -> String {552 let (bq, src) = ladder_b(a);553 let b = bq.val();554 let floor = mono_floor(b);555 let q0 = ladder_wall(b);556 format!(557 "| {} | {} | {} | {} | {} |",558 a.show(),559 bq.show(),560 src,561 integer_or_bound(q0, ladder_exact(q0, b)),562 integer_or_bound(floor, floor < DYADIC_TOP)563 )564}565566fn ladder_trend() -> String {567 let mut parts: Vec<String> = Vec::new();568 for (n, d) in LADDER {569 let b = ladder_b(Rat::new(n, d)).0.val();570 parts.push(fixed(ceil_units(ladder_wall(b).log10(), 2), 2));571 }572 parts.join(" ")573}574575fn render_ladder_coro(a: Rat) -> String {576 let (bq, src) = ladder_b(a);577 format!(578 "| {} | {} | {} | {} |",579 a.show(),580 bq.show(),581 src,582 max_excluded_b(CORO_Q, bq.val())583 )584}585586// MAIN587588fn main() {589 println!("mertens-numerology");590 println!(591 "guard {GUARD:e}; alpha_q truncated down at 6 digits, c_q rounded up at 5, delta_q rounded down at 5"592 );593 println!();594 println!("TABLE, one excluded digit");595 println!("| `q` | `alpha_q` | `c_q` (proved, up) | `delta_q` (down) | closes |");596 println!("|---|---|---|---|---|");597 for q in BASES {598 let r = row(q, 1);599 assert!(c_exp(q as f64, 1.0) >= 0.0, "l^1 floor broken at q = {q}");600 assert_eq!(601 r.closes,602 gap(q as f64, 1.0) > 0.0,603 "tests disagree at q = {q}"604 );605 println!("{}", render(&r));606 }607 println!();608 println!("WALL, one excluded digit");609 assert!(!closes(3689.0, 1.0), "3689 must not close");610 assert!(closes(3690.0, 1.0), "3690 must close");611 println!(612 "- closes fails at `q = 3689`, holds at `q = 3690`; gap `<= {}` and `>= {}`",613 fixed(ceil_units(gap(3689.0, 1.0), 6), 6),614 fixed(floor_units(gap(3690.0, 1.0), 6), 6)615 );616 let (arg, worst, min_gap) = gap_sweep(SWEEP_LO, SWEEP_HI, 1);617 assert!(worst > 0.0, "gap not stepwise increasing");618 assert!(min_gap > 0.0, "gap not positive on the sweep");619 println!(620 "- gap `(q-1) q^(-3/4) - PB_q(1)` steps up at all {} steps of `{SWEEP_LO}..{SWEEP_HI}`; smallest step `>= {}` at `q = {arg}`",621 SWEEP_HI - SWEEP_LO,622 fixed(floor_units(worst, 8), 8)623 );624 println!();625 println!("MARGIN AT THE WALL, one excluded digit, 10 significant digits");626 println!(627 "- relative guard {SCI_GUARD:e}; `delta_q` from the cancellation-free form `ln(1 + gap_q/PB_q) / (alpha_q ln q)`"628 );629 for q in WALL {630 let qf = q as f64;631 let up = !closes(qf, 1.0);632 let rel = if up { "<=" } else { ">=" };633 println!(634 "- `q = {q}`: `delta_q {rel} {}`, `gap_q(1) {rel} {}`",635 sci(delta_stable(qf, 1.0), 9, up),636 sci(gap(qf, 1.0), 9, up)637 );638 }639 println!("- the fifth-digit table columns above cannot display either sign");640 println!();641 println!("COST-OUT, `delta_q` against the defect exponent `m/(2(q-m) ln q)` of the level-`x^(alpha/2)` bound");642 println!("| `q` | `m` | `delta_q` (down) | defect (up) | saving beats defect |");643 println!("|---|---|---|---|---|");644 for q in COST {645 println!("{}", render_cost(q, 1));646 }647 for q in CORO {648 println!("{}", render_cost(q, max_excluded(q)));649 }650 let (first, down) = crossover(SWEEP_LO, SWEEP_HI);651 assert!(first > 0, "no crossover in the scan");652 println!(653 "- least `q` in `{SWEEP_LO}..{SWEEP_HI}` with `delta_q > 1/(2(q-1) ln q)`: `q = {first}`"654 );655 println!(656 "- the difference `delta_q - 1/(2(q-1) ln q)` rises at every one of the {} steps of that scan ({down} exceptions); monotonicity beyond the scan is not proved",657 SWEEP_HI - SWEEP_LO658 );659 println!();660 println!(661 "SHARPENED COST-OUT, the GRH rung re-costed at every one-step constant this desk proves"662 );663 println!(664 "`1 + Phi_q/q` is step 3; `Psi_q` is the sharpening with the phase identity, `q >= 17`;"665 );666 println!("`Psi'_q` is the same with the chord kernel bound in place of `Phi_q`, `q >= 36`");667 println!("`H` is the harmonic upper bound `ln n + gamma + 1/(2n)` in both, and the `36` is that convention:");668 println!("with the harmonic number itself the monotone step's floor `Psi'_q >= (1 + pi) q/2` first holds at `q = 37`");669 println!("the wall is the least `q` with `(q-1) q^(-3/4)` above the constant, the crossing the least `q` where");670 println!("`delta_q` exceeds the defect `1/(2(q-1) ln q)`; no wall is inherited and no crossing is inherited");671 println!("| constant | excluded digit | wall `q_0(1/2)` | crossing `q` | `delta_q` there (down) | defect (up) | both up-sets |");672 println!("|---|---|---|---|---|---|---|");673 for (name, digits, c) in [674 ("`1 + Phi_q/q`", "any", Cst::Step3),675 ("`Psi_q`", "any", Cst::Sharp(false)),676 ("`Psi_q`", "`0` or `q-1`", Cst::Sharp(true)),677 ("`Psi'_q`", "any", Cst::Chord(false)),678 ("`Psi'_q`", "`0` or `q-1`", Cst::Chord(true)),679 ] {680 println!("{}", render_sharp(name, digits, c));681 }682 println!(683 "- each row is scanned exhaustively from its own floor, `q >= 3`, `17` and `36`, to `{SHARP_HI}` for the wall and `{SHARP_COST_HI}` for the crossing"684 );685 println!("- every crossing sits within five steps of its own wall, so a lower wall costs out at once and is never inherited");686 println!("- the rungs above `a = 1/2` are not costed here and stay owed");687 println!();688 println!("M-COROLLARY, largest `m` with `PB_q(m) < (q-m) q^(-3/4)`");689 println!("| `q` | max `m` | `alpha_q` | `c_q` (proved, up) | `delta_q` (down) |");690 println!("|---|---|---|---|---|");691 for q in CORO {692 let m = max_excluded(q);693 let r = row(q, m);694 assert!(r.closes, "corollary maximum must close at q = {q}");695 assert!(696 !closes(q as f64, (m + 1) as f64),697 "maximum not maximal at q = {q}"698 );699 println!("{}", render_coro(&r));700 }701 println!();702 println!("LADDER, one excluded digit, `q_0(a)` = least `q >= 3` with `(q-1) q^(-b(a)) - PB_q(1) > 0`");703 println!("- `b(a)` is the smaller of Baker-Harman and Zhang where both apply; `Q(b)` is the proved monotone floor");704 println!("| `a` | `b(a)` | source | `q_0(a)` | `Q(b)` |");705 println!("|---|---|---|---|---|");706 for (n, d) in LADDER {707 let a = Rat::new(n, d);708 let b = ladder_b(a).0.val();709 let floor = mono_floor(b);710 let q0 = ladder_wall(b);711 assert!(712 floor < q0,713 "monotone floor above the wall at a = {}",714 a.show()715 );716 for q in 3..SWEEP_LO {717 assert!(718 gap_b(q as f64, b, 1.0) < 0.0,719 "close below 3690 at a = {}",720 a.show()721 );722 }723 if floor > SWEEP_LO as f64 {724 assert!(725 u_bound(SWEEP_LO as f64, b) < 0.0,726 "floor range unclear at a = {}",727 a.show()728 );729 assert!(730 u_bound(floor, b) < 0.0,731 "floor range unclear at a = {}",732 a.show()733 );734 }735 println!("{}", render_ladder(a));736 }737 assert!(738 ladder_b(Rat::new(1, 2)).0.same(Rat::new(3, 4)),739 "the GRH rung must read b = 3/4"740 );741 assert_eq!(742 ladder_wall(0.75),743 SWEEP_LO as f64,744 "the GRH rung must be the wall"745 );746 println!(747 "- `log10 q_0(a)` across the rungs, rounded up: {}",748 ladder_trend()749 );750 println!("- `gap_q(a, 1)` steps up at every `q >= Q(b)`, and `Q(b) < q_0(a)` at every rung, so it steps up from `q_0(a)` on");751 println!(752 "- no `q < {SWEEP_LO}` closes at any printed rung; `<=` marks a rung whose wall is past `2^53` or whose neighbouring gap steps fall under the {EXACT_ULPS}-ulp noise floor"753 );754 println!();755 println!("LADDER M-COROLLARY, largest `m` with `PB_q(m) < (q-m) q^(-b(a))` at `q = 10^7`");756 println!("| `a` | `b(a)` | source | max `m` |");757 println!("|---|---|---|---|");758 for (n, d) in LADDER {759 let a = Rat::new(n, d);760 let b = ladder_b(a).0.val();761 if ladder_wall(b) >= CORO_Q {762 continue;763 }764 let m = max_excluded_b(CORO_Q, b);765 assert!(766 m >= 1,767 "rung past its wall must admit m = 1 at a = {}",768 a.show()769 );770 assert!(771 gap_b(CORO_Q, b, (m + 1) as f64) <= 0.0,772 "maximum not maximal at a = {}",773 a.show()774 );775 println!("{}", render_ladder_coro(a));776 }777}778779// TESTS780781#[cfg(test)]782mod tests {783 use super::*;784785 fn kernel_sum(q: u64, t: f64) -> f64 {786 let qf = q as f64;787 let num = (PI * t).sin().abs();788 let mut s = 0.0;789 for r in 0..q {790 let x = (t + r as f64) / qf;791 let d = if x <= 0.5 { x } else { 1.0 - x };792 let den = (PI * d).sin();793 s += if den < 1e-12 { qf } else { num / den };794 }795 s796 }797798 #[test]799 fn table_rows_pinned() {800 let got: Vec<String> = BASES.iter().map(|&q| render(&row(q, 1))).collect();801 let want = [802 "| 1000 | 0.999855 | 0.28087 | -0.03102 | no |",803 "| 2000 | 0.999934 | 0.26335 | -0.01342 | no |",804 "| 3000 | 0.999958 | 0.25430 | -0.00434 | no |",805 "| 3689 | 0.999966 | 0.24997 | -0.00001 | no |",806 "| 3690 | 0.999967 | 0.24997 | 0.00000 | yes |",807 "| 5000 | 0.999976 | 0.24393 | 0.00605 | yes |",808 "| 10^4 | 0.999989 | 0.23141 | 0.01858 | yes |",809 "| 10^5 | 0.999999 | 0.19906 | 0.05094 | yes |",810 "| 10^6 | 0.999999 | 0.17589 | 0.07411 | yes |",811 "| 10^9 | 0.999999 | 0.13305 | 0.11695 | yes |",812 ];813 assert_eq!(got, want);814 }815816 #[test]817 fn wall_sign_flip() {818 assert!(!closes(3689.0, 1.0));819 assert!(closes(3690.0, 1.0));820 assert!(gap(3689.0, 1.0) < 0.0);821 assert!(gap(3690.0, 1.0) > 0.0);822 for q in 3..3689u64 {823 assert!(824 !closes(q as f64, 1.0),825 "unexpected close below the wall at {q}"826 );827 }828 }829830 #[test]831 fn gap_steps_up_and_stays_positive() {832 let (arg, worst, min_gap) = gap_sweep(SWEEP_LO, SWEEP_HI, 1);833 assert!(worst > 0.0);834 assert!(min_gap > 0.0);835 assert_eq!(arg, SWEEP_HI - 2);836 assert_eq!(fixed(floor_units(worst, 8), 8), "0.00003172");837 for q in SWEEP_LO..5_000 {838 assert!(839 gap((q + 1) as f64, 1.0) > gap(q as f64, 1.0),840 "step down at {q}"841 );842 }843 }844845 #[test]846 fn closes_matches_gap_sign() {847 for q in 3..20_000u64 {848 assert_eq!(849 closes(q as f64, 1.0),850 gap(q as f64, 1.0) > 0.0,851 "at q = {q}"852 );853 }854 }855856 #[test]857 fn corollary_maxima() {858 assert_eq!(max_excluded(10_000), 6);859 assert_eq!(max_excluded(100_000), 78);860 assert_eq!(max_excluded(1_000_000), 451);861 let got: Vec<String> = CORO862 .iter()863 .map(|&q| render_coro(&row(q, max_excluded(q))))864 .collect();865 let want = [866 "| 10^4 | 6 | 0.999934 | 0.24865 | 0.00129 |",867 "| 10^5 | 78 | 0.999932 | 0.24972 | 0.00022 |",868 "| 10^6 | 451 | 0.999967 | 0.24994 | 0.00003 |",869 ];870 assert_eq!(got, want);871 }872873 #[test]874 fn margin_at_the_wall_pinned() {875 let got: Vec<String> = WALL876 .iter()877 .flat_map(|&q| {878 let qf = q as f64;879 let up = !closes(qf, 1.0);880 [sci(delta_stable(qf, 1.0), 9, up), sci(gap(qf, 1.0), 9, up)]881 })882 .collect();883 let want = [884 "-2.395807653e-6",885 "-1.533059397e-4",886 "5.863425182e-6",887 "3.752213034e-4",888 ];889 assert_eq!(got, want);890 }891892 #[test]893 fn margin_bounds_are_safe() {894 for q in WALL {895 let qf = q as f64;896 let up = !closes(qf, 1.0);897 for (v, s) in [898 (delta_stable(qf, 1.0), sci(delta_stable(qf, 1.0), 9, up)),899 (gap(qf, 1.0), sci(gap(qf, 1.0), 9, up)),900 ] {901 let p: f64 = s.parse().unwrap();902 if up {903 assert!(p >= v, "not an upper bound at q = {q}: {s}");904 } else {905 assert!(p <= v, "not a lower bound at q = {q}: {s}");906 }907 assert!(908 (p - v).abs() <= 1e-8 * v.abs(),909 "bound too loose at q = {q}"910 );911 }912 }913 }914915 #[test]916 fn delta_stable_matches_delta() {917 for q in BASES.iter().chain(COST.iter()) {918 let qf = *q as f64;919 let (a, b) = (delta(qf, 1.0), delta_stable(qf, 1.0));920 assert!(921 (a - b).abs() <= 1e-9 * b.abs(),922 "delta forms disagree at q = {q}"923 );924 }925 for q in CORO {926 let (qf, mf) = (q as f64, max_excluded(q) as f64);927 let (a, b) = (delta(qf, mf), delta_stable(qf, mf));928 assert!(929 (a - b).abs() <= 1e-9 * b.abs(),930 "delta forms disagree at q = {q}"931 );932 }933 }934935 #[test]936 fn cost_out_rows_pinned() {937 let mut got: Vec<String> = COST.iter().map(|&q| render_cost(q, 1)).collect();938 got.extend(CORO.iter().map(|&q| render_cost(q, max_excluded(q))));939 let want = [940 "| 3689 | 1 | -2.39581e-6 | 1.65072e-5 | no |",941 "| 3690 | 1 | 5.86342e-6 | 1.65022e-5 | no |",942 "| 5000 | 1 | 6.05211e-3 | 1.17434e-5 | yes |",943 "| 10^4 | 1 | 1.85809e-2 | 5.42923e-6 | yes |",944 "| 10^5 | 1 | 5.09409e-2 | 4.34299e-7 | yes |",945 "| 10^6 | 1 | 7.41160e-2 | 3.61913e-8 | yes |",946 "| 10^9 | 1 | 1.16951e-1 | 2.41275e-11 | yes |",947 "| 10^4 | 6 | 1.29265e-3 | 3.25917e-5 | yes |",948 "| 10^5 | 78 | 2.20253e-4 | 3.39015e-5 | yes |",949 "| 10^6 | 451 | 3.14081e-5 | 1.63296e-5 | yes |",950 ];951 assert_eq!(got, want);952 }953954 #[test]955 fn sharpened_cost_out_is_pinned() {956 let got: Vec<String> = [957 ("`1 + Phi_q/q`", "any", Cst::Step3),958 ("`Psi_q`", "any", Cst::Sharp(false)),959 ("`Psi_q`", "`0` or `q-1`", Cst::Sharp(true)),960 ("`Psi'_q`", "any", Cst::Chord(false)),961 ("`Psi'_q`", "`0` or `q-1`", Cst::Chord(true)),962 ]963 .iter()964 .map(|(n, d, c)| render_sharp(n, d, *c))965 .collect();966 let want = [967 "| `1 + Phi_q/q` | any | 3690 | 3692 | 1.69819e-5 | 1.64921e-5 | yes |",968 "| `Psi_q` | any | 2446 | 2450 | 3.90025e-5 | 2.61622e-5 | yes |",969 "| `Psi_q` | `0` or `q-1` | 1812 | 1815 | 3.85795e-5 | 3.67324e-5 | yes |",970 "| `Psi'_q` | any | 1499 | 1502 | 4.80703e-5 | 4.55409e-5 | yes |",971 "| `Psi'_q` | `0` or `q-1` | 1032 | 1036 | 8.24697e-5 | 6.95785e-5 | yes |",972 ];973 assert_eq!(got, want);974 }975976 #[test]977 fn step3_row_reproduces_the_ladder() {978 let (w, up) = onestep_wall(Cst::Step3, SHARP_HI);979 assert_eq!((w, up), (3_690, true));980 assert_eq!(onestep_cross(Cst::Step3, SHARP_COST_HI).0, 3_692);981 for q in [3_690u64, 10_000, 100_000] {982 let qf = q as f64;983 assert!((onestep(qf, Cst::Step3) - pb(qf, 1.0)).abs() < 1e-12 * pb(qf, 1.0));984 }985 }986987 #[test]988 fn the_chord_constant_saves_half_a_base() {989 for q in (36u64..4_000).step_by(1) {990 let qf = q as f64;991 let want = if q % 2 == 0 {992 psi(qf) - qf / 2.0 + 2.0 / PI993 } else {994 psi_chord(qf)995 };996 assert!(997 (psi_chord(qf) - want).abs() <= 1e-9 * psi(qf),998 "chord identity fails at q = {q}"999 );1000 assert!(psi_chord(qf) < psi(qf), "chord not sharper at q = {q}");1001 assert!(1002 psi_chord(qf) >= (1.0 + PI) * qf / 2.0,1003 "monotone floor fails at q = {q}"1004 );1005 }1006 assert!(psi_chord(35.0) < (1.0 + PI) * 35.0 / 2.0);1007 }10081009 #[test]1010 fn every_sharpened_wall_lies_under_its_predecessor() {1011 let walls: Vec<u64> = [1012 Cst::Step3,1013 Cst::Sharp(false),1014 Cst::Sharp(true),1015 Cst::Chord(false),1016 Cst::Chord(true),1017 ]1018 .iter()1019 .map(|&c| onestep_wall(c, SHARP_HI).0)1020 .collect();1021 for w in walls.windows(2) {1022 assert!(w[1] < w[0], "walls not decreasing: {walls:?}");1023 }1024 for (i, &c) in [1025 Cst::Step3,1026 Cst::Sharp(false),1027 Cst::Sharp(true),1028 Cst::Chord(false),1029 Cst::Chord(true),1030 ]1031 .iter()1032 .enumerate()1033 {1034 let (x, up, _, _) = onestep_cross(c, SHARP_COST_HI);1035 assert!(up, "crossing not an up-set at row {i}");1036 assert!(1037 x >= walls[i] && x - walls[i] <= 5,1038 "crossing far from wall at row {i}"1039 );1040 }1041 }10421043 #[test]1044 fn crossover_is_pinned() {1045 let (first, down) = crossover(SWEEP_LO, SWEEP_HI);1046 assert_eq!(first, 3_692);1047 assert_eq!(down, 0);1048 for q in [3_690u64, 3_691] {1049 let qf = q as f64;1050 assert!(1051 delta_stable(qf, 1.0) < defect(qf, 1.0),1052 "premature crossing at {q}"1053 );1054 }1055 assert!(delta_stable(3_692.0, 1.0) > defect(3_692.0, 1.0));1056 }10571058 #[test]1059 fn sci_is_directional() {1060 for x in [1.0e-6, -1.0e-6, 9.999999999e-3, 1.23456789e7, -7.5e-11] {1061 for d in [3u32, 5, 9] {1062 let up: f64 = sci(x, d, true).parse().unwrap();1063 let dn: f64 = sci(x, d, false).parse().unwrap();1064 assert!(up >= x, "sci up failed on {x} at {d}");1065 assert!(dn <= x, "sci down failed on {x} at {d}");1066 let w = 10f64.powi(-(d as i32)) * x.abs() * 20.0;1067 assert!(up - dn <= w, "sci band too wide on {x} at {d}");1068 }1069 }1070 }10711072 #[test]1073 fn rounding_is_directional() {1074 for q in BASES {1075 let (qf, mf) = (q as f64, 1.0);1076 let r = row(q, 1);1077 assert!(fixed(r.c, 5).parse::<f64>().unwrap() >= c_exp(qf, mf) - GUARD);1078 assert!(fixed(r.delta, 5).parse::<f64>().unwrap() <= delta(qf, mf) + GUARD);1079 assert!(fixed(r.alpha, 6).parse::<f64>().unwrap() <= alpha(qf, mf) + GUARD);1080 assert!(c_exp(qf, mf) - fixed(r.c, 5).parse::<f64>().unwrap() > -1e-5);1081 }1082 }10831084 #[test]1085 fn l1_floor_holds() {1086 for q in 3..5_000u64 {1087 assert!(c_exp(q as f64, 1.0) >= 0.0, "negative c_q at {q}");1088 assert!(pb(q as f64, 1.0) >= 1.0, "PB below the l^1 floor at {q}");1089 }1090 }10911092 #[test]1093 fn the_chord_floor_carries_its_harmonic_convention() {1094 let exact = |n: f64| (1..=(n as u64)).map(|j| 1.0 / j as f64).sum::<f64>();1095 let psi_exact = |q: f64| {1096 let p = (q / 2.0).floor();1097 let h = exact(p - 1.0);1098 if (q as u64) % 2 == 0 {1099 (q / PI) * (2.0 * h - 1.0 + 1.0 / p) + (1.0 - 2.0 / PI) * q / 2.01100 } else {1101 (q / PI) * (2.0 * h - 1.0 + 2.0 / p)1102 + (1.0 - 2.0 / PI) * (q / 2.0 + 1.0 / (2.0 * q))1103 }1104 };1105 assert!(psi_chord(36.0) >= (1.0 + PI) * 18.0);1106 assert!(psi_exact(36.0) < (1.0 + PI) * 18.0);1107 for q in 37u64..4_000 {1108 let qf = q as f64;1109 assert!(1110 psi_exact(qf) >= (1.0 + PI) * qf / 2.0,1111 "harmonic-number floor fails at q = {q}"1112 );1113 }1114 }11151116 #[test]1117 fn harmonic_bound_dominates() {1118 let mut h = 0.0;1119 for n in 1..=2_000u64 {1120 h += 1.0 / n as f64;1121 assert!(harmonic_bound(n as f64) >= h, "harmonic bound fails at {n}");1122 }1123 }11241125 #[test]1126 fn kernel_bound_dominates_sampled_sup() {1127 for q in [50u64, 101, 200] {1128 let bound = phi(q as f64);1129 let mut worst: f64 = 0.0;1130 for i in 0..4_001u64 {1131 let t = i as f64 / 4_001.0;1132 worst = worst.max(kernel_sum(q, t));1133 }1134 assert!(1135 worst <= bound,1136 "kernel bound fails at q = {q}: {worst} > {bound}"1137 );1138 assert!(1139 worst > 0.8 * bound,1140 "kernel bound absurdly loose at q = {q}"1141 );1142 }1143 }11441145 #[test]1146 fn ladder_rows_pinned() {1147 let got: Vec<String> = LADDER1148 .iter()1149 .map(|&(n, d)| render_ladder(Rat::new(n, d)))1150 .collect();1151 let want = [1152 "| 1/2 | 3/4 | both | 3690 | 723 |",1153 "| 13/25 | 1417/1850 | Zhang | 8578 | 1486 |",1154 "| 11/20 | 913/1160 | Zhang | 33547 | 4754 |",1155 "| 4/7 | 4/5 | both | 92317 | 11221 |",1156 "| 3/5 | 4/5 | BH | 92317 | 11221 |",1157 "| 2/3 | 5/6 | BH | 3107080 | 216023 |",1158 "| 3/4 | 7/8 | BH | 6939524168 | 129458304 |",1159 "| 4/5 | 9/10 | BH | <= 3.09358e13 | 128606353005 |",1160 "| 9/10 | 19/20 | BH | <= 3.23663e34 | <= 1.73431e28 |",1161 "| 19/20 | 39/40 | BH | <= 9.24614e83 | <= 3.30712e68 |",1162 ];1163 assert_eq!(got, want);1164 }11651166 #[test]1167 fn ladder_trend_pinned() {1168 assert_eq!(1169 ladder_trend(),1170 "3.57 3.94 4.53 4.97 4.97 6.50 9.85 13.50 34.52 83.97"1171 );1172 }11731174 #[test]1175 fn last_down_step_below_the_floor() {1176 assert!((4.0 / PI) / 662.0 > 0.25 * 664f64.powf(-0.75));1177 let mut last = 0u64;1178 for q in 3..SWEEP_LO {1179 if gap((q + 1) as f64, 1.0) <= gap(q as f64, 1.0) {1180 last = q;1181 }1182 }1183 assert_eq!(last, 662);1184 }11851186 #[test]1187 fn general_b_floor_bound() {1188 let k = 2.0 / PI;1189 let slack = 2.0 * 723f64.powf(-0.75);1190 assert!(slack < 0.015, "step-A slack above 0.015: {slack}");1191 let konst = slack - 1.0 - 4.0 / PI - k * GAMMA + k * (1448.0f64 / 721.0).ln()1192 - (1.0 - k) * (721.0 / 723.0);1193 assert!(1194 konst < -2.544,1195 "general-b constant weaker than -2.544: {konst}"1196 );1197 assert!(4.0 * k - 2.544 > 0.0, "bracket not increasing on (0, 1/4]");1198 let bracket = 1.291 - k * (1.291 * 4.0f64).ln() - 2.544 / 4.0;1199 assert!(bracket < -0.39014, "bracket above -0.39014: {bracket}");1200 assert!(1201 4.0 * bracket < -1.56,1202 "floor gap bound weaker than -1.56: {bracket}"1203 );1204 let over =1205 3691f64.powf(-0.75) + k * ((3690.0f64 / 3689.0).ln() + 1.0 / 3689.0) + 0.727 / 3691.0;1206 assert!(over < 0.004, "majorant overshoot above 0.004: {over}");1207 assert!(u_bound(3690.0, 1417.0 / 1850.0) < -0.95);1208 let mut b = 0.75;1209 while b < 0.98 {1210 let floor = mono_floor(b);1211 assert!(1212 gap_b(floor, b, 1.0) < -1.56,1213 "gap at the floor above the bound at b = {b}"1214 );1215 if floor > SWEEP_LO as f64 {1216 assert!(1217 u_bound(SWEEP_LO as f64, b) < 0.0,1218 "majorant positive at b = {b}"1219 );1220 assert!(1221 u_bound(floor, b) < 0.0,1222 "majorant positive at the floor, b = {b}"1223 );1224 }1225 b += 0.005;1226 }1227 }12281229 #[test]1230 fn ladder_grh_rung_is_the_wall() {1231 let (b, src) = ladder_b(Rat::new(1, 2));1232 assert!(b.same(Rat::new(3, 4)));1233 assert_eq!(src, "both");1234 assert_eq!(ladder_wall(0.75), 3_690.0);1235 assert!(gap(3_690.0, 1.0) > 0.0 && gap(3_689.0, 1.0) < 0.0);1236 }12371238 #[test]1239 fn zhang_beats_baker_harman_strictly_inside() {1240 assert!(zhang(Rat::new(1, 2))1241 .unwrap()1242 .same(baker_harman(Rat::new(1, 2))));1243 assert!(zhang(Rat::new(4, 7))1244 .unwrap()1245 .same(baker_harman(Rat::new(4, 7))));1246 let (lo, hi) = (Rat::new(1, 2), Rat::new(4, 7));1247 let mut seen = 0u64;1248 for d in 2..=200i64 {1249 for n in 1..d {1250 let a = Rat::new(n, d);1251 if !lo.lt(a) || !a.lt(hi) {1252 continue;1253 }1254 let z = zhang(a).unwrap();1255 assert!(1256 z.lt(baker_harman(a)),1257 "Zhang not smaller at a = {}",1258 a.show()1259 );1260 assert_eq!(ladder_b(a).1, "Zhang");1261 seen += 1;1262 }1263 }1264 assert!(seen > 400, "too few interior rationals tested");1265 assert!(zhang(Rat::new(3, 5)).is_none());1266 }12671268 #[test]1269 fn ladder_floor_forces_steps_up() {1270 for (n, d) in LADDER {1271 let a = Rat::new(n, d);1272 let b = ladder_b(a).0.val();1273 let floor = mono_floor(b);1274 assert!(1275 mono_ok(floor, b),1276 "floor fails its own test at a = {}",1277 a.show()1278 );1279 if floor < DYADIC_TOP {1280 assert!(1281 !mono_ok(floor - 1.0, b),1282 "floor not least at a = {}",1283 a.show()1284 );1285 }1286 assert!(1287 floor < ladder_wall(b),1288 "floor above the wall at a = {}",1289 a.show()1290 );1291 let lo = floor.max(3_690.0);1292 if lo < 4.0e6 {1293 let mut q = lo;1294 while q < (lo + 20_000.0).min(4.0e6) {1295 assert!(1296 gap_b(q + 1.0, b, 1.0) > gap_b(q, b, 1.0),1297 "step down at q = {q}"1298 );1299 q += 1.0;1300 }1301 }1302 }1303 }13041305 #[test]1306 fn pb_step_below_the_constant() {1307 let mut worst = 0.0f64;1308 for q in 40..400_000u64 {1309 let qf = q as f64;1310 worst = worst.max((pb(qf + 1.0, 1.0) - pb(qf, 1.0)) * (qf - 2.0));1311 }1312 assert!(1313 worst < MONO_C,1314 "PB step above the monotone constant: {worst}"1315 );1316 assert!(worst > 1.27, "PB step bound absurdly loose: {worst}");1317 }13181319 #[test]1320 fn pb_low_is_a_lower_bound() {1321 for q in 3..20_000u64 {1322 let qf = q as f64;1323 assert!(pb_low(qf) <= pb(qf, 1.0), "pb_low above pb at q = {q}");1324 assert!(1325 gap_b(qf, 0.9, 1.0) <= u_bound(qf, 0.9),1326 "u_bound below gap at q = {q}"1327 );1328 }1329 }13301331 #[test]1332 fn ladder_below_the_floor_is_clear() {1333 for (n, d) in LADDER {1334 let a = Rat::new(n, d);1335 let b = ladder_b(a).0.val();1336 for q in 3..3_690u64 {1337 assert!(1338 gap_b(q as f64, b, 1.0) < 0.0,1339 "close below the wall at a = {}",1340 a.show()1341 );1342 }1343 let floor = mono_floor(b);1344 if floor > 3_690.0 {1345 assert!(u_bound(3_690.0, b) < 0.0);1346 assert!(u_bound(floor, b) < 0.0);1347 assert!(u_bound(3_690.0, b).max(u_bound(floor, b)) < -1.5);1348 }1349 }1350 }13511352 #[test]1353 fn ladder_wall_matches_exhaustive_scan() {1354 for (n, d) in LADDER {1355 let a = Rat::new(n, d);1356 let b = ladder_b(a).0.val();1357 let q0 = ladder_wall(b);1358 if q0 > 4.0e6 {1359 continue;1360 }1361 let mut found = 0.0f64;1362 let mut q = 3.0;1363 while q <= q0 {1364 if gap_b(q, b, 1.0) > 0.0 {1365 found = q;1366 break;1367 }1368 q += 1.0;1369 }1370 assert_eq!(1371 found,1372 q0,1373 "scan disagrees with bisection at a = {}",1374 a.show()1375 );1376 }1377 }13781379 #[test]1380 fn ladder_m_corollary_pinned() {1381 let got: Vec<String> = LADDER1382 .iter()1383 .map(|&(n, d)| Rat::new(n, d))1384 .filter(|a| ladder_wall(ladder_b(*a).0.val()) < CORO_Q)1385 .map(render_ladder_coro)1386 .collect();1387 let want = [1388 "| 1/2 | 3/4 | both | 1971 |",1389 "| 13/25 | 1417/1850 | Zhang | 1002 |",1390 "| 11/20 | 913/1160 | Zhang | 365 |",1391 "| 4/7 | 4/5 | both | 176 |",1392 "| 3/5 | 4/5 | BH | 176 |",1393 "| 2/3 | 5/6 | BH | 8 |",1394 ];1395 assert_eq!(got, want);1396 }13971398 #[test]1399 fn ladder_wall_rises_with_a() {1400 let mut prev = 0.0f64;1401 for (n, d) in LADDER {1402 let b = ladder_b(Rat::new(n, d)).0.val();1403 let q0 = ladder_wall(b);1404 assert!(q0 >= prev, "wall not rising at {n}/{d}");1405 prev = q0;1406 }1407 assert!(prev > 1.0e83);1408 }14091410 #[test]1411 fn asymptotic_shape() {1412 let q = 1e12;1413 let c = c_exp(q, 1.0);1414 let model = (q.ln().ln() + (2.0 / PI).ln()) / q.ln();1415 assert!((c - model).abs() < 0.01);1416 assert!(delta(q, 1.0) > 0.13 && delta(q, 1.0) < 0.25);1417 }1418}