constants.rs
8.6 kB · rust · 280 lines
1use crate::lattice::{Family, FAMILIES};2use crate::star::{law, slice};3use crate::sums::odds;4use mrlynum::series::{beta, chi4, dirichlet, zeta, CATALAN};5use num_bigint::BigInt;6use num_rational::{BigRational, Ratio};7use std::f64::consts::PI;89const TERMS: usize = 10_000_000;10const LAYERS: usize = 55;1112type Q = Ratio<i64>;1314struct Shape {15 a: Q,16 b: Q,17 c: Q,18 d: Q,19 e: Q,20 f: Q,21}2223fn shape(family: Family) -> Shape {24 let r = |p: i64, q: i64| Ratio::new(p, q);25 match family {26 Family::Carpet => Shape {27 a: r(1, 2),28 b: r(1, 8),29 c: r(1, 2),30 d: r(0, 1),31 e: r(0, 1),32 f: r(-1, 8),33 },34 Family::Net => Shape {35 a: r(1, 2),36 b: r(-1, 8),37 c: r(-1, 2),38 d: r(0, 1),39 e: r(0, 1),40 f: r(1, 8),41 },42 Family::Tree => Shape {43 a: r(1, 4),44 b: r(0, 1),45 c: r(1, 3),46 d: r(-1, 12),47 e: r(1, 6),48 f: r(-1, 6),49 },50 Family::Void => Shape {51 a: r(1, 4),52 b: r(0, 1),53 c: r(0, 1),54 d: r(-1, 4),55 e: r(1, 2),56 f: r(0, 1),57 },58 }59}6061fn chi(number: i64) -> i64 {62 -i64::from(chi4(number as usize))63}6465fn quasi(shape: &Shape, n: i64) -> Q {66 let (t, square) = (Ratio::from_integer(n), Ratio::from_integer(n * n));67 let s = Ratio::from_integer(chi(n));68 shape.a + shape.b * s + (shape.c + shape.d * s) / t + (shape.e + shape.f * s) / square69}7071fn text(value: Q) -> String {72 if *value.denom() == 1 {73 return value.numer().to_string();74 }75 format!("{}/{}", value.numer(), value.denom())76}7778fn big(value: Q) -> BigRational {79 BigRational::new(BigInt::from(*value.numer()), BigInt::from(*value.denom()))80}8182fn real(value: Q) -> f64 {83 *value.numer() as f64 / *value.denom() as f6484}8586fn sums(limit: usize) -> [f64; 4] {87 let mut out = [0.0; 4];88 for n in odds(limit) {89 let x = n as f64;90 let sign = f64::from(chi4(n));91 out[0] += sign;92 out[1] += sign / x;93 out[2] += sign / (x * x);94 out[3] += 1.0 / (x * x);95 }96 out97}9899fn shapes_match() {100 for family in FAMILIES {101 let form = shape(family);102 let mut matched = 0;103 let mut layers = 0;104 for n in odds(LAYERS) {105 matched += usize::from(law(family, n) == quasi(&form, n as i64));106 layers += 1;107 }108 println!(109 " {}: A = {} B = {} c = {} d = {} e = {} f = {}, rebuild the closed form {matched}/{layers}",110 family.name(),111 text(form.a),112 text(form.b),113 text(form.c),114 text(form.d),115 text(form.e),116 text(form.f)117 );118 }119}120121fn summed_identity() {122 println!(123 " the summed identity in exact rationals, counted hexagons against the character sums"124 );125 for family in FAMILIES {126 let form = shape(family);127 let zero = BigRational::from_integer(BigInt::from(0));128 let mut left = zero.clone();129 let (mut s0, mut s1, mut s2, mut s3) = (zero.clone(), zero.clone(), zero.clone(), zero);130 let mut classes = [0usize; 4];131 let mut counts = [0usize; 4];132 for n in odds(LAYERS) {133 let (side, sign) = (n as i64, Ratio::from_integer(i64::from(chi4(n))));134 let ink = slice(family, n).ink();135 left += big(ink - form.a - form.c / Ratio::from_integer(side));136 s0 += big(sign);137 s1 += big(sign / Ratio::from_integer(side));138 s2 += big(sign / Ratio::from_integer(side * side));139 s3 += big(Ratio::new(1, side * side));140 let right =141 -big(form.b) * &s0 - big(form.d) * &s1 + big(form.e) * &s3 - big(form.f) * &s2;142 let slot = (n % 8) / 2;143 counts[slot] += 1;144 classes[slot] += usize::from(left == right);145 }146 println!(147 " {}: n = 1,3,5,7 mod 8 read {}/{} {}/{} {}/{} {}/{}",148 family.name(),149 classes[0],150 counts[0],151 classes[1],152 counts[1],153 classes[2],154 counts[2],155 classes[3],156 counts[3]157 );158 }159}160161fn limits(catalan: f64) -> Vec<(Family, f64, f64)> {162 FAMILIES163 .into_iter()164 .map(|family| {165 let form = shape(family);166 let even =167 -real(form.d) * PI / 4.0 + real(form.e) * PI * PI / 8.0 - real(form.f) * catalan;168 (family, even, even - real(form.b))169 })170 .collect()171}172173fn ladder(catalan: f64) {174 let stops = [400usize, 1600, 3200];175 let deepest = stops[stops.len() - 1] + 3;176 for family in FAMILIES {177 let form = shape(family);178 let (b, d, e, f) = (real(form.b), real(form.d), real(form.e), real(form.f));179 let flat = d == 0.0 && e == 0.0;180 let base = -d * PI / 4.0 + e * PI * PI / 8.0 - f * catalan;181 let mut running = 0.0;182 let mut readings: Vec<(usize, f64)> = Vec::new();183 for count in 1..=deepest {184 let n = 2 * count as i64 - 1;185 running += real(law(family, n as usize) - form.a - form.c / Ratio::from_integer(n));186 if !stops.iter().any(|stop| count >= *stop && count <= stop + 3) {187 continue;188 }189 let target = base - if count % 2 == 0 { 0.0 } else { b };190 let scale = count as f64;191 let gap = (running - target) * scale;192 readings.push((count, if flat { gap * scale } else { gap }));193 }194 let power = if flat { "M^2" } else { "M" };195 let even = if flat { f / 8.0 } else { (d - e) / 4.0 };196 let odd = if flat { -f / 8.0 } else { (-d - e) / 4.0 };197 println!(198 " {}: gap * {power} -> {even:+.8} at even M and {odd:+.8} at odd M",199 family.name()200 );201 for chunk in readings.chunks(4) {202 let cells: Vec<String> = chunk203 .iter()204 .map(|(count, read)| format!("M = {count} {read:+.8}"))205 .collect();206 println!(" {}", cells.join(" "));207 }208 }209}210211fn split() {212 let mut cells: Vec<String> = Vec::new();213 for count in [400usize, 1600, 3200] {214 let [_, s1, _, s3] = sums(2 * count - 1);215 let gap = s1 + s3 / 4.0 - (PI / 4.0 + PI * PI / 32.0);216 cells.push(format!("M = {count} {:+.8}", gap * count as f64));217 }218 println!(219 " carpet split at even M only: gap * M -> {:+.8}, reading {}",220 -5.0 / 16.0,221 cells.join(" ")222 );223}224225pub fn run() {226 let catalan = beta(2.0, TERMS);227 let eisenstein = dirichlet(2.0, &[0, 1, -1], TERMS);228 let zeta3 = zeta(3.0, 2_000_000);229 println!("constants from their own series");230 println!(" G = {catalan:.10} L(2, chi_-3) = {eisenstein:.10} zeta(3) = {zeta3:.10}");231 println!(232 " carpet split pi/4 + pi^2/32 = {:.10}",233 PI / 4.0 + PI * PI / 32.0234 );235 println!(236 " carpet mean G/8 = {:.10} G/8 - 1/8 = {:.10}",237 catalan / 8.0,238 catalan / 8.0 - 0.125239 );240 println!(241 " void mean (pi + pi^2)/16 = {:.10}",242 (PI + PI * PI) / 16.0243 );244 println!(245 " flat stack pi^2 ln2/(7 zeta3) = {:.10}",246 PI * PI * 2f64.ln() / (7.0 * zeta3)247 );248 println!(249 " tree mean pi/48+pi^2/48+G/6 = {:.10}",250 PI / 48.0 + PI * PI / 48.0 + catalan / 6.0251 );252 println!("partial character sums over odd n <= N, M layers, from the exact ink laws");253 println!(" the carpet split is an identity at even M only, so it prints at even M only");254 for limit in [53usize, 55] {255 let [s0, s1, s2, s3] = sums(limit);256 let carpet = if (limit + 1) / 2 % 2 == 0 {257 format!("M(I1 - I3 + 1/4) = {:.10} ", s1 + s3 / 4.0)258 } else {259 String::new()260 };261 println!(262 " N = {limit}: {carpet}M(mean ink - 1/2 - eps/2) = {:.10} void M(mean - 1/4) = {:.10}",263 -s0 / 8.0 + s2 / 8.0,264 s1 / 4.0 + s3 / 2.0265 );266 }267 println!("the one-layer law M(mean I - A - c eps) = -B S - d s1 + e s3 - f s2");268 shapes_match();269 summed_identity();270 println!(" the limit -B[M odd] - d pi/4 + e pi^2/8 - f G, both classes of the layer count");271 for (family, even, odd) in limits(CATALAN) {272 println!(273 " {}: {even:.10} at even M, N = 3 mod 4 and {odd:.10} at odd M, N = 1 mod 4",274 family.name()275 );276 }277 println!(" the approach, every class of M mod 4 so every class of N mod 8");278 ladder(CATALAN);279 split();280}