research-franel.rs

1use figures::{ink, save, Board, Grid};2use mrlyrs::core::error::Result;3use mrlyrs::core::Color;4use mrlyrs::num::design::elements;5use mrlyrs::num::factor::{gcd, mobius_sieve};6use std::f64::consts::TAU;78const BASE: u64 = 3;9const DIGITS: [u64; 2] = [0, 1];10const ORDER: usize = 81;11const LIT: [usize; 16] = [1, 2, 6, 8, 14, 15, 18, 20, 28, 30, 31, 36, 37, 39, 40, 81];1213// THE FORM1415fn dilated(holds: &[bool], mu: &[i8]) -> Vec<i64> {16    let mut out = vec![0i64; ORDER + 1];17    for d in 1..=ORDER {18        for c in 1..=ORDER / d {19            if holds[d * c] {20                out[d] += mu[c] as i64;21            }22        }23    }24    out25}2627fn kernel(d: usize, e: usize) -> f64 {28    let g = gcd(d as u128, e as u128) as f64;29    g * g / (d as f64 * e as f64)30}3132fn nodes(dens: &[usize]) -> Vec<(usize, usize)> {33    let mut out = Vec::new();34    for &b in dens {35        for a in 1..=b {36            if gcd(a as u128, b as u128) == 1 {37                out.push((a, b));38            }39        }40    }41    out42}4344// THE MARKS4546fn wash(k: f64) -> Color {47    ink::mix(ink::ground(), ink::dim(), 0.015 + 0.7 * k.sqrt())48}4950fn lit(v: f64, low: f64, high: f64) -> Color {51    let s = (v.abs().ln() - low.ln()) / (high.ln() - low.ln());52    let hue = if v > 0.0 { ink::blue() } else { ink::orange() };53    ink::mix(ink::ground(), hue, 0.55 + 0.45 * s)54}5556fn main() -> Result<()> {57    let design: Vec<usize> = elements(BASE, &DIGITS, 5)58        .into_iter()59        .filter(|&n| n as usize <= ORDER)60        .map(|n| n as usize)61        .collect();62    assert_eq!(design.len(), 16);63    let mut holds = vec![false; ORDER + 1];64    for &b in &design {65        holds[b] = true;66    }67    let mu = mobius_sieve(ORDER);68    let x = dilated(&holds, &mu);69    let support: Vec<usize> = (1..=ORDER).filter(|&d| x[d] != 0).collect();70    assert_eq!(support, LIT);71    assert_eq!(x[1], -2);72    assert_eq!(x.iter().map(|v| v * v).sum::<i64>(), 19);7374    let mut terms = Vec::new();75    for &d in &support {76        for &e in &support {77            terms.push((d, e, kernel(d, e) * (x[d] * x[e]) as f64));78        }79    }80    assert_eq!(terms.iter().filter(|t| t.2 != 0.0).count(), 256);81    let form: f64 = terms.iter().map(|t| t.2).sum();82    assert!((form - 12.574_087).abs() < 1e-5);8384    let farey = nodes(&design);85    assert_eq!(farey.len(), 241);86    for k in 1..=12usize {87        let (mut re, mut im) = (0.0, 0.0);88        for &(a, b) in &farey {89            let angle = TAU * (k * a) as f64 / b as f64;90            re += angle.cos();91            im += angle.sin();92        }93        let divisor: i64 = (1..=k)94            .filter(|d| k % d == 0 && *d <= ORDER)95            .map(|d| d as i64 * x[d])96            .sum();97        assert!((re - divisor as f64).abs() < 1e-9);98        assert!(im.abs() < 1e-9);99    }100101    let low = terms.iter().fold(f64::MAX, |a, t| a.min(t.2.abs()));102    let high = terms.iter().fold(0.0f64, |a, t| a.max(t.2.abs()));103    assert!(low > 3e-4 && high == 4.0);104105    let mut board = Board::square();106    let frame = board.frame(0.08);107    let fine = Grid::new(frame, ORDER, ORDER, 0.34);108    let bold = Grid::new(frame, ORDER, ORDER, 0.08);109    for d in 1..=ORDER {110        for e in 1..=ORDER {111            fine.fill(&mut board, e - 1, d - 1, wash(kernel(d, e)));112        }113    }114    for &(d, e, v) in &terms {115        bold.fill(&mut board, e - 1, d - 1, lit(v, low, high));116    }117    save("research-franel", &board)?;118    Ok(())119}