coprime.rs

1use mrlynum::factor::coprime;2use std::thread;34pub const THREADS: u64 = 8;56fn row_slice(n: u32, rem: u64, step: u64) -> u64 {7    let mut count = 0u64;8    let mut m = 2 * rem + 1;9    while m < 1u64 << n {10        let rest = m - 1;11        let mut sub = rest;12        loop {13            if coprime(sub as usize + 1, m as usize) {14                count += 1;15            }16            if sub == 0 {17                break;18            }19            sub = (sub - 1) & rest;20        }21        m += 2 * step;22    }23    count24}2526fn complement_slice(n: u32, rem: u64, step: u64) -> u64 {27    let size = 1u64 << n;28    let mut count = 0u64;29    let mut i = rem;30    while i < size {31        let free = !i & (size - 1);32        let mut j = free;33        while j > i {34            if coprime(i as usize, j as usize) {35                count += 1;36            }37            j = (j - 1) & free;38        }39        i += step;40    }41    count42}4344fn spread(n: u32, slice: fn(u32, u64, u64) -> u64) -> u64 {45    let handles: Vec<thread::JoinHandle<u64>> = (0..THREADS)46        .map(|rem| thread::spawn(move || slice(n, rem, THREADS)))47        .collect();48    2 * handles.into_iter().map(|h| h.join().unwrap()).sum::<u64>()49}5051pub fn by_rows(n: u32) -> u64 {52    spread(n, row_slice)53}5455pub fn by_complement(n: u32) -> u64 {56    spread(n, complement_slice)57}5859pub fn by_pascal(n: u32) -> u64 {60    let mut row: Vec<u8> = vec![1];61    let mut total = 0u64;62    for m in 0..1u64 << n {63        for (k, entry) in row.iter().enumerate() {64            if *entry == 1 && coprime(k, m as usize - k) {65                total += 1;66            }67        }68        let mut next: Vec<u8> = Vec::with_capacity(row.len() + 1);69        next.push(1);70        for k in 1..row.len() {71            next.push((row[k - 1] + row[k]) & 1);72        }73        next.push(1);74        row = next;75    }76    total77}7879pub fn density(term: u64, n: u32) -> f64 {80    term as f64 / 3f64.powi(n as i32)81}8283pub fn limit() -> f64 {84    16.0 / (3.0 * std::f64::consts::PI * std::f64::consts::PI)85}