gaussian.rs

1use std::collections::HashSet;23pub fn least_factors(limit: usize) -> Vec<u32> {4    let mut out = vec![0u32; limit + 1];5    for value in 2..=limit {6        if out[value] == 0 {7            let mut multiple = value;8            while multiple <= limit {9                if out[multiple] == 0 {10                    out[multiple] = value as u32;11                }12                multiple += value;13            }14        }15    }16    out17}1819pub fn factorise(mut value: usize, least: &[u32]) -> Vec<(usize, usize)> {20    let mut out: Vec<(usize, usize)> = Vec::new();21    while value > 1 {22        let prime = least[value] as usize;23        let mut power = 0;24        while value % prime == 0 {25            value /= prime;26            power += 1;27        }28        out.push((prime, power));29    }30    out31}3233pub fn two_square(value: usize, least: &[u32]) -> bool {34    factorise(value, least)35        .iter()36        .all(|(prime, power)| prime % 4 != 3 || power % 2 == 0)37}3839pub fn primitive_norm(value: usize, least: &[u32]) -> bool {40    if value == 1 || value == 2 {41        return true;42    }43    let parts = factorise(value, least);44    parts.iter().all(|(prime, power)| match prime % 4 {45        1 => true,46        2 => *power == 1,47        _ => false,48    })49}5051pub fn new_disc(scale: usize, least: &[u32]) -> usize {52    let cap = 2 * scale * scale;53    let primes: Vec<usize> = factorise(scale, least).iter().map(|(p, _)| *p).collect();54    (1..=cap)55        .filter(|norm| two_square(*norm, least))56        .filter(|norm| primes.iter().all(|prime| norm % (prime * prime) != 0))57        .count()58}5960pub fn two_square_prefix(cap: usize, least: &[u32]) -> Vec<usize> {61    let mut out = vec![0usize; cap + 1];62    for norm in 1..=cap {63        out[norm] = out[norm - 1] + usize::from(two_square(norm, least));64    }65    out66}6768pub fn mobius_count(scale: usize, prefix: &[usize], least: &[u32]) -> usize {69    let primes: Vec<usize> = factorise(scale, least).iter().map(|(p, _)| *p).collect();70    let cap = 2 * scale * scale;71    let mut total = 0i64;72    for mask in 0..1usize << primes.len() {73        let mut divisor = 1usize;74        let mut sign = 1i64;75        for (index, prime) in primes.iter().enumerate() {76            if mask >> index & 1 == 1 {77                divisor *= prime;78                sign = -sign;79            }80        }81        total += sign * prefix[cap / (divisor * divisor)] as i64;82    }83    total as usize84}8586pub fn jordan(scale: usize, least: &[u32]) -> f64 {87    factorise(scale, least)88        .iter()89        .map(|(prime, _)| 1.0 - 1.0 / (prime * prime) as f64)90        .product()91}9293fn reduce(mut top: u64, mut bottom: u64) -> (u64, u64) {94    let mut a = top;95    let mut b = bottom;96    while b != 0 {97        let next = a % b;98        a = b;99        b = next;100    }101    top /= a;102    bottom /= a;103    (top, bottom)104}105106pub fn disc_norms(scale: usize) -> Vec<u64> {107    let cap = 2 * scale * scale;108    let mut seen = vec![false; cap + 1];109    let reach = (cap as f64).sqrt() as usize + 1;110    for a in 0..=reach {111        for b in 0..=reach {112            let norm = a * a + b * b;113            if norm >= 1 && norm <= cap {114                seen[norm] = true;115            }116        }117    }118    (1..=cap)119        .filter(|norm| seen[*norm])120        .map(|norm| norm as u64)121        .collect()122}123124pub fn box_norms(scale: usize) -> Vec<u64> {125    let cap = 2 * scale * scale;126    let mut seen = vec![false; cap + 1];127    for a in 0..=scale {128        for b in 0..=scale {129            let norm = a * a + b * b;130            if norm >= 1 {131                seen[norm] = true;132            }133        }134    }135    (1..=cap)136        .filter(|norm| seen[*norm])137        .map(|norm| norm as u64)138        .collect()139}140141pub fn union_counts(top: usize, boxed: bool) -> Vec<usize> {142    let mut seen: HashSet<(u64, u64)> = HashSet::new();143    let mut out = Vec::new();144    for scale in 1..=top {145        let norms = if boxed {146            box_norms(scale)147        } else {148            disc_norms(scale)149        };150        let square = (scale * scale) as u64;151        let mut fresh = 0;152        for norm in norms {153            if seen.insert(reduce(norm, square)) {154                fresh += 1;155            }156        }157        out.push(fresh);158    }159    out160}