Store squared limit to reduce operations.

More readable.

.gitignore

@@ -1,2 +1,3 @@
 /target
 /indivisible.1.gz
+callgrind.out.*

Cargo.lock

@@ -43,6 +43,37 @@ dependencies = [
  "vec_map",
 ]
 
+[[package]]
+name = "crossbeam-deque"
+version = "0.8.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51"
+dependencies = [
+ "crossbeam-epoch",
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-epoch"
+version = "0.9.18"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e"
+dependencies = [
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-utils"
+version = "0.8.21"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28"
+
+[[package]]
+name = "either"
+version = "1.15.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719"
+
 [[package]]
 name = "heck"
 version = "0.3.3"
@@ -65,6 +96,7 @@ dependencies = [
 name = "indivisible"
 version = "2.0.0"
 dependencies = [
+ "rayon",
  "structopt",
 ]
 
@@ -76,9 +108,9 @@ checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe"
 
 [[package]]
 name = "libc"
-version = "0.2.178"
+version = "0.2.182"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "37c93d8daa9d8a012fd8ab92f088405fb202ea0b6ab73ee2482ae66af4f42091"
+checksum = "6800badb6cb2082ffd7b6a67e6125bb39f18782f793520caee8cb8846be06112"
 
 [[package]]
 name = "proc-macro-error"
@@ -106,22 +138,42 @@ dependencies = [
 
 [[package]]
 name = "proc-macro2"
-version = "1.0.103"
+version = "1.0.106"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5ee95bc4ef87b8d5ba32e8b7714ccc834865276eab0aed5c9958d00ec45f49e8"
+checksum = "8fd00f0bb2e90d81d1044c2b32617f68fcb9fa3bb7640c23e9c748e53fb30934"
 dependencies = [
  "unicode-ident",
 ]
 
 [[package]]
 name = "quote"
-version = "1.0.42"
+version = "1.0.45"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a338cc41d27e6cc6dce6cefc13a0729dfbb81c262b1f519331575dd80ef3067f"
+checksum = "41f2619966050689382d2b44f664f4bc593e129785a36d6ee376ddf37259b924"
 dependencies = [
  "proc-macro2",
 ]
 
+[[package]]
+name = "rayon"
+version = "1.11.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "368f01d005bf8fd9b1206fb6fa653e6c4a81ceb1466406b81792d87c5677a58f"
+dependencies = [
+ "either",
+ "rayon-core",
+]
+
+[[package]]
+name = "rayon-core"
+version = "1.13.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91"
+dependencies = [
+ "crossbeam-deque",
+ "crossbeam-utils",
+]
+
 [[package]]
 name = "strsim"
 version = "0.8.0"
@@ -174,9 +226,9 @@ dependencies = [
 
 [[package]]
 name = "unicode-ident"
-version = "1.0.22"
+version = "1.0.24"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9312f7c4f6ff9069b165498234ce8be658059c6728633667c526e27dc2cf1df5"
+checksum = "e6e4313cd5fcd3dad5cafa179702e2b244f760991f45397d14d4ebf38247da75"
 
 [[package]]
 name = "unicode-segmentation"

Cargo.toml

@@ -10,3 +10,4 @@ keywords = [ "primes", "math" ]
 
 [dependencies]
 structopt = "0.3"
+rayon = "1.11"

README.md

@@ -1,38 +1,29 @@
 # Indivisible
 
-Indivisible is an optimized prime number generator and tester written in Rust.
+Indivisible is an optimized prime number generator and tester.
 
 ## Build & Installation
 
-To build the project you will require the Rust compiler and build
-system, `cargo`. At which point you simply run `cargo build` in the root
-directory of the project. To create an optimized release build append
-the `--release` flag to the previous command.
+Build a binary release using the following command:
 
-Once a release build has been compiled, you may install Indivisible and
-the manpage documentation by running the `install` script. In a similar
-manner, you can run the `uninstall` script to remove the previously
-installed binary and documentation.
+```console
+$ cargo build -r
+```
 
-## Usage
+Once compiled, you can install it using the `install` script. By default this
+will install to `/usr/local/{bin,share}`, but you can optionally change the
+`PREFIX` used as follows:
 
-The purpose of Indivisible is to find the nth prime and all the primes
-before it. The basic usage is `indivisible <n>` where `n` is the ordinal
-of the prime you'd like to find. To display all primes before `n`, you
-can run verbose mode by using the `--verbose` or simply `-v` option.
+```console
+$ PREFIX="~/.local" ./install
+```
 
-Since Indivisible generates primes using previously computed primes, you
-can also import prime numbers previously computed with the `--import` or
-`-i` option. To store already computed primes you are expected to use
-piping like any UNIX user would expect. Here is an example:
+Additionally there is an `uninstall` script available. If you have used a custom
+`PREFIX` this must also be specified to uninstall in the same manner. To see
+where you've installed it you can always use the following command:
 
-```bash
-# store first 100 primes in ./primes
-indivisible -v 100 > ./primes
-# appends next 400 primes
-indivisible -i ./primes -v 500 >> ./primes
-# display the 600th prime
-indivisible -i ./primes 600
+```console
+$ which indivisible
 ```
 
 ## Legacy

benchmark.sh

@@ -0,0 +1,53 @@
+#!/bin/bash
+
+EXE="./target/release/indivisible"
+OPTIONS=()
+TRIALS=20
+
+if ! [ -f "$EXE" ]
+then
+	>&2 echo "Release build not available. Please run 'cargo build -r'."
+	exit 1
+fi
+
+if ! command -v calc &>/dev/null
+then
+	>&2 echo "Missing 'calc' program. Please install it for this script."
+	exit 1
+fi
+
+while getopts "t:s:j:" opt
+do
+	case "$opt" in
+		s)
+			OPTIONS=("${OPTIONS[@]}" -s "$OPTARG")
+			SIEVE=$OPTARG
+			;;
+		t)
+			TRIALS="$OPTARG"
+			;;
+		j)
+			OPTIONS=("${OPTIONS[@]}" -j "$OPTARG")
+			JOBS=$OPTARG
+			;;
+		*)
+			>&2 echo "Uknown option $opt"
+			exit 1
+			;;
+	esac
+done
+
+echo "Calculating primes up to 1,000,000,000"
+echo "Trials: $TRIALS"
+echo "Sieve segment size: ${SIEVE:-"default"}"
+echo "Number of jobs: ${JOBS:-"default"}"
+TOTAL="0"
+for _ in $(seq "$TRIALS")
+do
+	TIME=$(command time -f "%e" "$EXE" "${OPTIONS[@]}" 1000000000 2>&1 >/dev/null)
+	TOTAL=$(calc "$TOTAL + $TIME")
+done
+
+AVG=$(calc "$TOTAL / $TRIALS")
+
+echo "Average time: ${AVG}s"

indivisible.1

@@ -6,7 +6,9 @@ indivisible \- prime number generator and tester
 [\fIOPTIONS\fR] <\fIn\fR>
 
 .SH "DESCRIPTION"
-Indivisible is a program for working with prime numbers. Its default behaviour will attempt to find the \fInth\fR prime.
+Indivisible is a program for working with prime numbers. It can either generate
+primes up to a limit \fIn\fR or test to see if \fIn\fR is a prime. It makes use
+of a segmented sieve of Eratosthenes.
 
 .SH "OPTIONS"
 .TP
@@ -19,6 +21,9 @@ Import prime numbers from \fIFILE\fR
 \fB\-t\fR, \fB\-\-test\fR
 Test if n is a prime instead of generation
 .TP
+\fB\-s\fR, \fB\-\-sieve\fR <\fISIZE\fR>
+Set a custom sieve size
+.TP
 \fB\-v\fR, \fB\-\-verbose\fR
 Print all found primes
 .TP

src/main.rs

@@ -21,7 +21,13 @@ use std::fs::File;
 use std::io::{BufRead, BufReader};
 use std::path::PathBuf;
 use std::process;
+use std::sync::{Arc, RwLock};
+use std::sync::mpsc::{Sender, Receiver};
+use std::sync::mpsc;
 use structopt::StructOpt;
+use rayon::ThreadPoolBuilder;
+
+mod worker;
 
 #[derive(StructOpt)]
 #[structopt(about = "A prime number generator and tester.")]
@@ -32,26 +38,33 @@ struct Opt {
 	import:Option<PathBuf>,
 	#[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,
+	#[structopt(help = "Max of the prime to generate or number to test for primality")]
+	num:usize,
+	#[structopt(short, long, name = "SIZE", default_value = "10485760", help = "Set a custom sieve size")]
+	sieve:usize,
 	#[structopt(short, long, name = "n", default_value = "1", help = "Number of threads to spawn")]
-	jobs:u64,
+	jobs:usize,
+}
+
+struct BatchResult {
+	batch_id: u64,
+	primes: Vec<u64>,
 }
 
 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);
-	}
+	let mut prime_list_raw = Vec::new();
 
 	if opts.import.is_some() {
 		let in_file = File::open(opts.import.unwrap()).unwrap();
 		let reader = BufReader::new(in_file);
 		for p in reader.lines().into_iter() {
-			prime_list.push_back(p.unwrap().parse().unwrap());
+			let prime:u64 = p.unwrap().parse().unwrap();
+			if (prime as usize) > opts.num {
+				break;
+			}
+
+			prime_list_raw.push(prime);
 		}
 	}
 
@@ -59,57 +72,102 @@ fn main() {
 		eprintln!("Invalid value for num: {}", opts.num);
 		process::exit(1);
 	}
-	if opts.num > 2 {
-		arr[2] = true;
+
+	let pool = ThreadPoolBuilder::new().num_threads(opts.jobs).build().unwrap();
+	let mut pending_tasks = 0;
+	/* force u64 to avoid overflow if sieve size / opts.num is larger than i32::MAX.
+	 * This is unlikely, but possible.
+	 */
+	let mut batch_num: u64 = 0;
+	let mut last_batch_id: u64 = 0;
+	let mut remaining_batches = VecDeque::<BatchResult>::new();
+	let (tx, rx): (Sender<BatchResult>, Receiver<BatchResult>) = mpsc::channel();
+
+	let mut start:usize = if prime_list_raw.is_empty() {
+		2
+	} else {
+		(*prime_list_raw.last().unwrap() + 1) as usize
+	};
+	let mut end = if start + opts.sieve < opts.num {
+		start + opts.sieve
+	} else {
+		opts.num + 1
+	};
+	let mut squared_limit = 0;
+
+	let prime_list = Arc::new(RwLock::new(prime_list_raw));
+	loop {
+		if (start < opts.num) && (squared_limit >= (end as u64) ||
+			(prime_list.read().unwrap().is_empty() && pending_tasks == 0))
+		{
+			let prime_list_clone = Arc::clone(&prime_list);
+			let tx_clone = tx.clone();
+			batch_num += 1;
+			let batch_id = batch_num;
+			pool.spawn(move || {
+				let res = worker::work_segment(&prime_list_clone, start, end);
+				tx_clone.send(BatchResult{batch_id, primes: res}).unwrap();
+			});
+			pending_tasks += 1;
+			start += opts.sieve;
+			end = if end + opts.sieve < opts.num {
+				end + opts.sieve
+			} else {
+				opts.num + 1
+			};
+		} else if pending_tasks > 0 {
+			let res = rx.recv().unwrap();
+			pending_tasks -= 1;
+			if res.batch_id != last_batch_id + 1 {
+				if remaining_batches.is_empty() || (res.batch_id > remaining_batches.back().unwrap().batch_id) {
+					remaining_batches.push_back(res);
+				} else if res.batch_id < remaining_batches.front().unwrap().batch_id {
+					remaining_batches.push_front(res);
+				} else {
+					let mut i = 0;
+					while i < remaining_batches.len() {
+						if res.batch_id < remaining_batches[i].batch_id {
+							remaining_batches.insert(i, res);
+							break;
 						}
-	if opts.num > 3 {
-		arr[3] = true;
+						i += 1;
+					}
+				}
+			} else {
+				last_batch_id = res.batch_id;
+
+				if opts.verbose {
+					for p in &res.primes {
+						println!("{}", *p);
+					}
+				}
+				prime_list.write().unwrap().append(&mut res.primes.clone());
+
+				loop {
+					match remaining_batches.pop_front_if(|br| br.batch_id == last_batch_id + 1) {
+						Some(br) => {
+							last_batch_id = br.batch_id;
+							if opts.verbose {
+								for p in &br.primes {
+									println!("{}", *p);
+								}
+							}
+							prime_list.write().unwrap().append(&mut br.primes.clone());
+						},
+						None => break,
+					};
 				}
 
-	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 as usize) && (n1 % 12 == 1 || n1 % 12 == 5) {
-				arr[n1] = !arr[n1];
-			}
-
-			let n2 = ((3 * x * x) + (y * y)) as usize;
-			if n2 <= (opts.num as usize) && n2 % 12 == 7 {
-				arr[n2] = !arr[n2];
-			}
-
-			let n3 = ((3 * x * x) - (y * y)) as usize;
-			if x > y && n3 <= (opts.num as usize) && n3 % 12 == 11 {
-				arr[n3] = !arr[n3];
+				squared_limit = prime_list.read().unwrap().last().unwrap().pow(2);
 			}
+		} else {
+			// no more tasks to be added, and all tasks have been processed
+			break;
 		}
 	}
 
-	for i in 5..=(f64::sqrt(opts.num as f64) as u64 + 1) as usize {
-		if !arr[i as usize] {
-			continue;
-		}
-
-		let mut j = i * i;
-		while j <= (opts.num as usize) {
-			arr[j] = false;
-			j += i * i;
-		}
-	}
-
-	for i in 2..=(opts.num as usize) {
-		if arr[i] {
-			if opts.verbose && !opts.test {
-				println!("{}", i);
-			}
-			prime_list.push_back(i as u64);
-		}
-	}
-
-	if !opts.verbose && !opts.test {
-		println!("{}", prime_list.back().unwrap());
-	} else if opts.test {
-		if *prime_list.back().unwrap() == opts.num {
+	if opts.test {
+		if *prime_list.read().unwrap().last().unwrap() == (opts.num as u64) {
 			if opts.verbose {
 				println!("{} is prime", opts.num);
 			}
@@ -120,5 +178,7 @@ fn main() {
 			}
 			process::exit(1);
 		}
+	} else if !opts.verbose {
+		println!("{}", prime_list.read().unwrap().last().unwrap());
 	}
 }

src/worker.rs

@@ -0,0 +1,68 @@
+/*
+ * 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::sync::RwLock;
+
+/**
+ * @brief Work on a segment.
+ *
+ * @param known_primes:&Vec<u64> List of known primes at least until sqrt(end).
+ * @param start:usize Beginning of the segment (inclusive).
+ * @param end:usize End of the segment (exclusive).
+ *
+ * @return List of primes found in segment.
+ */
+pub fn work_segment(known_primes:&RwLock<Vec<u64>>, start:usize, end:usize) -> Vec<u64> {
+	let mut sieve = vec![true; end - start];
+	let mut found_primes = Vec::new();
+
+	let sqrt_end = (end as f64).sqrt() as usize;
+
+	for p in known_primes.read().unwrap().iter() {
+		let prime = *p as usize;
+		if prime > sqrt_end {
+			break;
+		}
+
+		let modu = start % prime;
+		let mut mult = if modu == 0 {
+			start
+		} else {
+			start + prime - modu
+		};
+		while mult < end {
+			sieve[mult - start] = false;
+			mult += prime;
+		}
+	}
+
+	for i in 0..(end - start) {
+		if sieve[i] {
+			let prime = i + start;
+			found_primes.push(prime as u64);
+
+			let mut mult = prime * prime;
+			while mult < end {
+				sieve[mult - start] = false;
+				mult += prime;
+			}
+		}
+	}
+
+	found_primes
+}

test.sh

@@ -17,9 +17,12 @@ fi
 tests=0
 passed=0
 
+## TEST
 ((tests++))
-echo -n "${tests}: Find 5th prime number..."
-if [[ $("$BINARY" 5) == 11 ]]
+echo -n "${tests}: Find all prime numbers before 70..."
+expect="2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 "
+
+if [ "$("$BINARY" -v 70 | tr '\n' ' ')" = "$expect" ]
 then
 	echo " pass"
 	((passed++))
@@ -27,6 +30,20 @@ else
 	echo " FAIL"
 fi
 
+## TEST
+((tests++))
+echo -n "${tests}: Find last prime before 1,000,000,000..."
+expect="999999937"
+
+if [ "$("$BINARY" 1000000000)" = "$expect" ]
+then
+	echo " pass"
+	((passed++))
+else
+	echo " FAIL"
+fi
+
+## TEST
 ((tests++))
 echo -n "${tests}: 11 is prime..."
 if "$BINARY" -t 11
@@ -37,6 +54,7 @@ else
 	echo " FAIL"
 fi
 
+## TEST
 ((tests++))
 echo -n "${tests}: 9 is not prime..."
 if ! "$BINARY" -t 9
@@ -47,4 +65,10 @@ else
 	echo " FAIL"
 fi
 
+## RESULTS
 echo "Results: $passed/$tests"
+
+if [ $passed -ne $tests ]
+then
+	exit 1
+fi