Implement sieve of Atkin.

This commit is contained in:
2025-12-05 13:16:37 +01:00
parent 65314d52ac
commit 05e18e5aef
3 changed files with 41 additions and 169 deletions

View File

@@ -1,82 +0,0 @@
/*
* Copyright (C) 2025 Nicolás Ortega Froysa <nicolas@ortegas.org>
* Author: Nicolás Ortega Froysa <nicolas@ortegas.org>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
pub struct CandidateGenerator {
base:u64,
// first use of the base
first_use:bool,
}
impl CandidateGenerator {
pub fn new() -> Self {
CandidateGenerator {
base: 0,
first_use: true,
}
}
pub fn calc_base(&mut self, last_prime:u64) {
if last_prime == 2 {
self.base = 0;
self.first_use = false;
} else if last_prime == 3 {
self.base = 6;
self.first_use = true;
} else {
let modulo = last_prime % 6;
if modulo == 1 {
self.base = last_prime + 5;
self.first_use = true;
} else if modulo == 5 {
self.base = last_prime + 1;
self.first_use = false;
} else {
panic!("Invalid last prime {}" , last_prime);
}
}
}
pub fn next(&mut self) -> u64 {
/*
* All primes, except 2 and 3, will be equal to (n * 6 ± 1). This avoids
* multiples of three, optimizing our counting.
*/
let val;
if self.base != 0 {
if self.first_use {
val = self.base - 1;
} else {
val = self.base + 1;
}
} else {
if self.first_use {
val = 2;
} else {
val = 3;
}
}
if !self.first_use {
self.base += 6;
}
self.first_use = !self.first_use;
val
}
}

View File

@@ -23,10 +23,6 @@ use std::path::PathBuf;
use std::process;
use structopt::StructOpt;
mod candidate;
use candidate::CandidateGenerator;
mod worker;
#[derive(StructOpt)]
#[structopt(about = "A prime number generator and tester.")]
struct Opt {
@@ -37,7 +33,7 @@ struct Opt {
#[structopt(short, long, help = "Test if num is prime instead of generation")]
test:bool,
#[structopt(help = "Ordinal of the prime to generate or number to test for primality")]
num:u64,
num:usize,
#[structopt(short, long, name = "n", default_value = "1", help = "Number of threads to spawn")]
jobs:u64,
}
@@ -46,6 +42,10 @@ fn main() {
let opts = Opt::from_args();
let mut prime_list = VecDeque::<u64>::new();
let mut arr = Vec::new();
for _ in 0..=opts.num {
arr.push(false);
}
if opts.import.is_some() {
let in_file = File::open(opts.import.unwrap()).unwrap();
@@ -55,62 +55,52 @@ fn main() {
}
}
if opts.num == 0 {
if opts.num < 2 {
eprintln!("Invalid value for num: {}", opts.num);
process::exit(1);
}
if opts.num > 2 {
arr[2] = true;
}
if opts.num > 3 {
arr[3] = true;
}
if opts.test && *prime_list.back().unwrap_or(&0) >= opts.num {
for i in prime_list.iter() {
if *i == opts.num {
process::exit(0)
for x in 1..=(f64::sqrt(opts.num as f64) as u64 + 1) {
for y in 1..=(f64::sqrt(opts.num as f64) as u64 + 1) {
let n1 = ((4 * x * x) + (y * y)) as usize;
if n1 <= opts.num && (n1 % 12 == 1 || n1 % 12 == 5) {
arr[n1] = !arr[n1];
}
let n2 = ((3 * x * x) + (y * y)) as usize;
if n2 <= opts.num && n2 % 12 == 7 {
arr[n2] = !arr[n2];
}
let n3 = ((3 * x * x) - (y * y)) as usize;
if x > y && n3 <= opts.num && n3 % 12 == 11 {
arr[n3] = !arr[n3];
}
}
process::exit(1)
} else if !opts.test && prime_list.len() >= opts.num as usize {
let res = *prime_list.get(opts.num as usize).unwrap();
println!("{}", res);
} else {
let mut cand_gen = CandidateGenerator::new();
if !prime_list.is_empty() {
cand_gen.calc_base(*prime_list.back().unwrap());
}
for i in 5..=(f64::sqrt(opts.num as f64) as u64 + 1) as usize {
if !arr[i as usize] {
continue;
}
loop {
let cand = cand_gen.next();
if opts.test && cand > opts.num {
break;
}
let mut is_prime = true;
for p in prime_list.iter() {
if cand % *p == 0 {
is_prime = false;
break;
}
}
if is_prime {
prime_list.push_back(cand);
if opts.verbose {
println!("{}", cand);
}
if !opts.test && prime_list.len() == opts.num as usize {
break;
}
}
let mut j = i * i;
while j <= opts.num {
arr[j] = false;
j += i * i;
}
}
if opts.test {
if *prime_list.back().unwrap() == opts.num {
process::exit(0)
} else {
process::exit(1)
}
} else if !opts.verbose {
let last_prime = *prime_list.back().unwrap();
println!("{}", last_prime);
for i in 2..=opts.num {
if arr[i] {
println!("{}", i);
prime_list.push_back(i as u64);
}
}
}

View File

@@ -1,36 +0,0 @@
/*
* Copyright (C) 2025 Nicolás Ortega Froysa <nicolas@ortegas.org>
* Author: Nicolás Ortega Froysa <nicolas@ortegas.org>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
use std::cell::RefCell;
use std::collections::VecDeque;
use std::rc::Rc;
pub struct Worker {
prime_list:Rc<RefCell<VecDeque<u64>>>,
}
impl Worker {
pub fn new(primes_list:Rc<RefCell<VecDeque<u64>>>) -> Worker {
Worker {
prime_list: primes_list,
}
}
pub fn run(&mut self) {
}
}