79 lines
2.0 KiB
Rust
79 lines
2.0 KiB
Rust
fn moebius(mut x: u64) -> i8 {
|
|
let mut prime_count = 0;
|
|
|
|
// If x is divisible by the given factor this macro counts the factor and divides it out.
|
|
// It then returns zero if x is still divisible by the factor.
|
|
macro_rules! divide_x_by {
|
|
($factor:expr) => {
|
|
if x % $factor == 0 {
|
|
x /= $factor;
|
|
prime_count += 1;
|
|
if x % $factor == 0 {
|
|
return 0;
|
|
}
|
|
}
|
|
};
|
|
}
|
|
|
|
// Handle 2 and 3 separately,
|
|
divide_x_by!(2);
|
|
divide_x_by!(3);
|
|
|
|
// then use a wheel sieve to check the remaining factors <= √x.
|
|
for i in (5..=isqrt(x)).step_by(6) {
|
|
divide_x_by!(i);
|
|
divide_x_by!(i + 2);
|
|
}
|
|
|
|
// There can exist one prime factor larger than √x,
|
|
// in that case we can check if x is still larger than one, and then count it.
|
|
if x > 1 {
|
|
prime_count += 1;
|
|
}
|
|
|
|
if prime_count % 2 == 0 {
|
|
1
|
|
} else {
|
|
-1
|
|
}
|
|
}
|
|
|
|
/// Returns the largest integer smaller than or equal to `√n`
|
|
const fn isqrt(n: u64) -> u64 {
|
|
if n <= 1 {
|
|
n
|
|
} else {
|
|
let mut x0 = u64::pow(2, n.ilog2() / 2 + 1);
|
|
let mut x1 = (x0 + n / x0) / 2;
|
|
while x1 < x0 {
|
|
x0 = x1;
|
|
x1 = (x0 + n / x0) / 2;
|
|
}
|
|
x0
|
|
}
|
|
}
|
|
|
|
fn main() {
|
|
const ROWS: u64 = 10;
|
|
const COLS: u64 = 20;
|
|
println!(
|
|
"Values of the Möbius function, μ(x), for x between 0 and {}:",
|
|
COLS * ROWS
|
|
);
|
|
for i in 0..ROWS {
|
|
for j in 0..=COLS {
|
|
let x = COLS * i + j;
|
|
let μ = moebius(x);
|
|
if μ >= 0 {
|
|
// Print an extra space if there's no minus sign in front of the output
|
|
// in order to align the numbers in a nice grid.
|
|
print!(" ");
|
|
}
|
|
print!("{μ} ");
|
|
}
|
|
println!();
|
|
}
|
|
let x = u64::MAX;
|
|
println!("\nμ({x}) = {}", moebius(x));
|
|
}
|