rays.rs

1use crate::sums::{mean, odds};2use mrlynum::series::chi4;34fn triangle(value: f64) -> f64 {5    let rest = value.rem_euclid(2.0);6    1.0 - 2.0 * rest.min(2.0 - rest)7}89fn average(point: f64, cap: usize, twisted: bool) -> f64 {10    let terms: Vec<f64> = odds(cap)11        .map(|n| {12            let sign = if twisted { f64::from(chi4(n)) } else { 1.0 };13            sign * triangle(n as f64 * point)14        })15        .collect();16    mean(&terms)17}1819pub fn twisted() {20    println!("(1/M) sum over odd n <= N of chi4(n) T(n x), the untwisted average beside it");21    let probes = [22        0.0,23        1.0 / 3.0,24        2.0 / 3.0,25        0.2,26        1.0 / 7.0,27        0.5,28        0.25,29        1.0 / 9.0,30        2f64.sqrt() - 1.0,31    ];32    for point in probes {33        println!(34            "  x = {point:.6}: N=55 {:+.5}  N=2001 {:+.6}  N=40001 {:+.7}  untwisted N=40001 {:+.6}",35            average(point, 55, true),36            average(point, 2001, true),37            average(point, 40001, true),38            average(point, 40001, false)39        );40    }41}4243pub const RESOLUTION: usize = 200003;4445fn square(value: f64) -> f64 {46    if (value.floor() as i64) % 2 == 0 {47        1.048    } else {49        -1.050    }51}5253fn field(cap: usize, first: &[f64], second: &[f64]) -> f64 {54    let mut running = vec![0.0f64; first.len()];55    let mut layers = 0;56    for n in odds(cap) {57        let twist = -f64::from(chi4(n));58        let scale = n as f64;59        for (slot, cell) in running.iter_mut().enumerate() {60            let left = square(scale * first[slot]);61            let right = square(scale * second[slot]);62            *cell += 0.5 + twist / 8.0 + (left + right) / 4.0 - twist * left * right / 8.0;63        }64        layers += 1;65    }66    mean(&running) / layers as f6467}6869pub fn crosshairs() {70    let points: Vec<f64> = (0..RESOLUTION)71        .map(|slot| (slot as f64 + 0.5) / RESOLUTION as f64)72        .collect();73    let generic: Vec<f64> = points74        .iter()75        .map(|point| (point * 2f64.sqrt()).rem_euclid(1.0))76        .collect();77    let shifted: Vec<f64> = points78        .iter()79        .map(|point| (point + 1.0 / 3.0).rem_euclid(1.0))80        .collect();81    println!("coarse cut model, {RESOLUTION} samples along a line, excess over the background");82    for cap in [55usize, 555, 5555] {83        let background = field(cap, &points, &generic);84        let diagonal = field(cap, &points, &points) - background;85        let offset = field(cap, &points, &shifted) - background;86        println!("  N = {cap:4}: background {background:.6}  A=C line {diagonal:+.6}  A-C=1/3 line {offset:+.6}");87    }88    let control: Vec<f64> = points89        .iter()90        .map(|point| (point * 3f64.sqrt()).rem_euclid(1.0))91        .collect();92    println!("coarse crosshair at A = a/q, predicted (-1)^a/(4q) for odd q and 0 for even q");93    for cap in [55usize, 5555] {94        let background = field(cap, &control, &points);95        for (a, q) in [(1, 3), (2, 3), (1, 5), (2, 5), (1, 7), (1, 2), (1, 4)] {96            let line = vec![a as f64 / q as f64 + 1e-12; RESOLUTION];97            let excess = field(cap, &line, &points) - background;98            let predicted = if q % 2 == 1 {99                (if a % 2 == 0 { 1.0 } else { -1.0 }) / (4 * q) as f64100            } else {101                0.0102            };103            println!("  N = {cap:4}  A = {a}/{q}: excess {excess:+.6}  predicted {predicted:+.6}");104        }105    }106}