Performance testing.

Dummy commit which I will revert later to see which implementation of the sieve of Atkin is faster.

src/main.rs

@@ -168,6 +168,11 @@ fn add_skels(problem_number: i32) {
 }
 
 fn main() {
+    for _ in 0..10000 {
+        solve_and_time_one(1);
+        solve_and_time_one(2);
+    }
+    return;
     let args = std::env::args().collect::<Vec<String>>();
     if args.len() <= 1 {
         solve_and_time_all();

src/solutions/even_fibonacci_numbers.rs

@@ -1,10 +1,4 @@
 use crate::utils;
-
 pub fn solve() -> i64 {
-    // Since we need only even terms, start from 2 and take every third term.
-    let fibonacci = utils::Fibonacci::new(2, 3);
-    let sum: i64 = fibonacci.step_by(3).take_while(|&num| num < 4000000).sum();
-
-    assert_eq!(sum, 4613732);
-    sum
+    utils::SieveOfAtkin2::new(100000000).is_prime(99999989) as i64
 }

src/solutions/multiples_of_3_or_5.rs

@@ -1,6 +1,4 @@
+use crate::utils;
 pub fn solve() -> i64 {
-    let sum = (1..1000).filter(|num| num % 3 == 0 || num % 5 == 0).sum();
-
-    assert_eq!(sum, 233168);
-    sum
+    utils::SieveOfAtkin::new(100000000).is_prime(99999989) as i64
 }

src/utils.rs

@@ -620,49 +620,39 @@ impl SieveOfAtkin {
     }
     fn algorithm_4_1(&mut self, delta: i32, f: i32, g: i32, h: i32) {
         let (mut x, mut y0, mut k0) = (f as i64, g as i64, h as i64);
-        let sieve_len = self.sieve.len() as i64;
-        if k0 < sieve_len {
-            let radicand = (x.pow(2) - 15 * (k0 - sieve_len)) as f64;
-            let n = ((radicand.sqrt() - x as f64) / 15.0).ceil() as i64;
-            (k0, x) = (k0 + 2 * x * n + 15 * n.pow(2), x + 15 * n);
+        while k0 < self.sieve.len() as i64 {
+            (k0, x) = (k0 + 2 * x + 15, x + 15);
         }
         loop {
             (k0, x) = (k0 - 2 * x + 15, x - 15);
             if x <= 0 {
                 return;
             }
-            if k0 < 0 {
-                let radicand = (y0.pow(2) - 60 * k0) as f64;
-                let n = ((radicand.sqrt() - y0 as f64) / 30.0).ceil() as i64;
-                (k0, y0) = (k0 + y0 * n + 15 * n.pow(2), y0 + 30 * n);
+            while k0 < 0 {
+                (k0, y0) = (k0 + y0 + 15, y0 + 30);
             }
             let (mut k, mut y) = (k0, y0);
-            while k < sieve_len {
+            while k < self.sieve.len() as i64 {
                 self.sieve[k as usize] ^= 1u16 << SieveOfAtkin::SHIFTS[delta as usize];
                 (k, y) = (k + y + 15, y + 30);
             }
         }
     }
     fn algorithm_4_2(&mut self, delta: i32, f: i32, g: i32, h: i32) {
         let (mut x, mut y0, mut k0) = (f as i64, g as i64, h as i64);
-        let sieve_len = self.sieve.len() as i64;
-        if k0 < sieve_len {
-            let radicand = (x.pow(2) - 20 * (k0 - sieve_len)) as f64;
-            let n = ((radicand.sqrt() - x as f64) / 10.0).ceil() as i64;
-            (k0, x) = (k0 + x * n + 5 * n.pow(2), x + 10 * n);
+        while k0 < self.sieve.len() as i64 {
+            (k0, x) = (k0 + x + 5, x + 10);
         }
         loop {
             (k0, x) = (k0 - x + 5, x - 10);
             if x <= 0 {
                 return;
             }
-            if k0 < 0 {
-                let radicand = (y0.pow(2) - 60 * k0) as f64;
-                let n = ((radicand.sqrt() - y0 as f64) / 30.0).ceil() as i64;
-                (k0, y0) = (k0 + y0 * n + 15 * n.pow(2), y0 + 30 * n);
+            while k0 < 0 {
+                (k0, y0) = (k0 + y0 + 15, y0 + 30);
             }
             let (mut k, mut y) = (k0, y0);
-            while k < sieve_len {
+            while k < self.sieve.len() as i64 {
                 self.sieve[k as usize] ^= 1u16 << SieveOfAtkin::SHIFTS[delta as usize];
                 (k, y) = (k + y + 15, y + 30);
             }
@@ -671,11 +661,7 @@ impl SieveOfAtkin {
     fn algorithm_4_3(&mut self, delta: i32, f: i32, g: i32, h: i32) {
         let (mut x, mut y0, mut k0) = (f as i64, g as i64, h as i64);
         loop {
-            let sieve_len = self.sieve.len() as i64;
-            if k0 >= sieve_len {
-                let radicand = (y0.pow(2) + 60 * (k0 - sieve_len)) as f64;
-                let n = ((radicand.sqrt() - y0 as f64) / 30.0).floor() as i64;
-                (k0, y0) = (k0 - y0 * n - 15 * n.pow(2), y0 + 30 * n);
+            while k0 >= self.sieve.len() as i64 {
                 if x <= y0 {
                     return;
                 }
@@ -1394,3 +1380,222 @@ mod tests {
         }
     }
 }
+
+/// A. O. L. Atkin and D. J. Bernstein, "Prime Sieves Using Binary Quadratic
+/// Forms", in Mathematics of Computation, vol. 73, no. 246, pp. 1023–1030.
+/// Generate prime numbers using the sieve of Atkin. This sieves prime numbers
+/// up to 1000000000 nearly twice as fast as my wheel-factorised sieve of
+/// Eratosthenes implementation (which I have now removed). It only determines
+/// the primality of numbers coprime to 60, because other numbers are
+/// guaranteed to be composite. (Exceptions 2, 3 and 5 are handled separately.)
+pub struct SieveOfAtkin2 {
+    limit: usize,
+    limit_rounded: usize,
+    limit_rounded_isqrt: usize,
+    sieve: Vec<u16>,
+}
+impl SieveOfAtkin2 {
+    // Consecutive differences between coprime residues modulo 60: 1, 7, 11,
+    // 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53 and 59.
+    const OFFSETS: [usize; 16] = [6, 4, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 4, 6, 2];
+    // Position of the bit indicating the primality of a coprime residue modulo
+    // 60 in a 16-element bitfield. For non-coprime residues, the value is 16.
+    const SHIFTS: [u8; 60] = [
+        16, 0, 16, 16, 16, 16, 16, 1, 16, 16, 16, 2, 16, 3, 16, 16, 16, 4, 16, 5, 16, 16, 16, 6, 16, 16, 16, 16, 16,
+        7, 16, 8, 16, 16, 16, 16, 16, 9, 16, 16, 16, 10, 16, 11, 16, 16, 16, 12, 16, 13, 16, 16, 16, 14, 16, 16, 16,
+        16, 16, 15,
+    ];
+}
+impl SieveOfAtkin2 {
+    /// Construct the sieve of Atkin up to and including the given number.
+    ///
+    /// * `limit` - Non-strict upper bound.
+    ///
+    /// -> Sieve of Atkin.
+    pub fn new(limit: usize) -> SieveOfAtkin2 {
+        // Strict upper bound divisible by 60.
+        let limit_rounded = (limit - limit % 60)
+            .checked_add(60)
+            .expect("overflow detected; argument too large");
+        let mut sieve_of_atkin = SieveOfAtkin2 {
+            limit,
+            limit_rounded,
+            limit_rounded_isqrt: isqrt(limit_rounded as i64) as usize,
+            sieve: vec![0; limit / 60 + 1],
+        };
+        sieve_of_atkin.init();
+        sieve_of_atkin
+    }
+    fn init(&mut self) {
+        self.algorithm_3_1(1);
+        self.algorithm_3_1(13);
+        self.algorithm_3_1(17);
+        self.algorithm_3_1(29);
+        self.algorithm_3_1(37);
+        self.algorithm_3_1(41);
+        self.algorithm_3_1(49);
+        self.algorithm_3_1(53);
+        self.algorithm_3_2(7);
+        self.algorithm_3_2(19);
+        self.algorithm_3_2(31);
+        self.algorithm_3_2(43);
+        self.algorithm_3_3(11);
+        self.algorithm_3_3(23);
+        self.algorithm_3_3(47);
+        self.algorithm_3_3(59);
+
+        // Mark composite all numbers divisible by the squares of primes.
+        let mut num: usize = 1;
+        let mut offset = SieveOfAtkin2::OFFSETS.iter().cycle();
+        'sieve: for sieve_idx in 0..self.sieve.len() {
+            for shift in 0..16 {
+                if self.sieve[sieve_idx] >> shift & 1 == 1 {
+                    let num_sqr = num.pow(2);
+                    for multiple in (num_sqr..)
+                        .step_by(num_sqr)
+                        .take_while(|&multiple| multiple < self.limit_rounded)
+                    {
+                        self.sieve[multiple / 60] &= !(1u32 << SieveOfAtkin2::SHIFTS[multiple % 60]) as u16;
+                    }
+                }
+                num += offset.next().unwrap();
+                if num > self.limit_rounded_isqrt {
+                    break 'sieve;
+                }
+            }
+        }
+    }
+    fn algorithm_3_1(&mut self, delta: i32) {
+        for f in 1..=15 {
+            for g in (1..=30).step_by(2) {
+                let quadratic = 4 * f * f + g * g;
+                if delta == quadratic % 60 {
+                    self.algorithm_4_1(delta, f, g, quadratic / 60);
+                }
+            }
+        }
+    }
+    fn algorithm_3_2(&mut self, delta: i32) {
+        for f in (1..=10).step_by(2) {
+            for g in [2, 4, 8, 10, 14, 16, 20, 22, 26, 28] {
+                let quadratic = 3 * f * f + g * g;
+                if delta == quadratic % 60 {
+                    self.algorithm_4_2(delta, f, g, quadratic / 60);
+                }
+            }
+        }
+    }
+    fn algorithm_3_3(&mut self, delta: i32) {
+        for (f, gstart) in (1..=10).zip([2, 1].into_iter().cycle()) {
+            for g in (gstart..=30).step_by(2) {
+                let quadratic = 3i32 * f * f - g * g;
+                // Remainder can be negative, so perform modulo operation.
+                if delta == quadratic.rem_euclid(60) {
+                    self.algorithm_4_3(delta, f, g, quadratic.div_euclid(60));
+                }
+            }
+        }
+    }
+    fn algorithm_4_1(&mut self, delta: i32, f: i32, g: i32, h: i32) {
+        let (mut x, mut y0, mut k0) = (f as i64, g as i64, h as i64);
+        let sieve_len = self.sieve.len() as i64;
+        if k0 < sieve_len {
+            let radicand = (x.pow(2) - 15 * (k0 - sieve_len)) as f64;
+            let n = ((radicand.sqrt() - x as f64) / 15.0).ceil() as i64;
+            (k0, x) = (k0 + 2 * x * n + 15 * n.pow(2), x + 15 * n);
+        }
+        loop {
+            (k0, x) = (k0 - 2 * x + 15, x - 15);
+            if x <= 0 {
+                return;
+            }
+            if k0 < 0 {
+                let radicand = (y0.pow(2) - 60 * k0) as f64;
+                let n = ((radicand.sqrt() - y0 as f64) / 30.0).ceil() as i64;
+                (k0, y0) = (k0 + y0 * n + 15 * n.pow(2), y0 + 30 * n);
+            }
+            let (mut k, mut y) = (k0, y0);
+            while k < sieve_len {
+                self.sieve[k as usize] ^= 1u16 << SieveOfAtkin2::SHIFTS[delta as usize];
+                (k, y) = (k + y + 15, y + 30);
+            }
+        }
+    }
+    fn algorithm_4_2(&mut self, delta: i32, f: i32, g: i32, h: i32) {
+        let (mut x, mut y0, mut k0) = (f as i64, g as i64, h as i64);
+        let sieve_len = self.sieve.len() as i64;
+        if k0 < sieve_len {
+            let radicand = (x.pow(2) - 20 * (k0 - sieve_len)) as f64;
+            let n = ((radicand.sqrt() - x as f64) / 10.0).ceil() as i64;
+            (k0, x) = (k0 + x * n + 5 * n.pow(2), x + 10 * n);
+        }
+        loop {
+            (k0, x) = (k0 - x + 5, x - 10);
+            if x <= 0 {
+                return;
+            }
+            if k0 < 0 {
+                let radicand = (y0.pow(2) - 60 * k0) as f64;
+                let n = ((radicand.sqrt() - y0 as f64) / 30.0).ceil() as i64;
+                (k0, y0) = (k0 + y0 * n + 15 * n.pow(2), y0 + 30 * n);
+            }
+            let (mut k, mut y) = (k0, y0);
+            while k < sieve_len {
+                self.sieve[k as usize] ^= 1u16 << SieveOfAtkin2::SHIFTS[delta as usize];
+                (k, y) = (k + y + 15, y + 30);
+            }
+        }
+    }
+    fn algorithm_4_3(&mut self, delta: i32, f: i32, g: i32, h: i32) {
+        let (mut x, mut y0, mut k0) = (f as i64, g as i64, h as i64);
+        loop {
+            let sieve_len = self.sieve.len() as i64;
+            if k0 >= sieve_len {
+                let radicand = (y0.pow(2) + 60 * (k0 - sieve_len)) as f64;
+                let n = ((radicand.sqrt() - y0 as f64) / 30.0).floor() as i64;
+                (k0, y0) = (k0 - y0 * n - 15 * n.pow(2), y0 + 30 * n);
+                if x <= y0 {
+                    return;
+                }
+                (k0, y0) = (k0 - y0 - 15, y0 + 30);
+            }
+            let (mut k, mut y) = (k0, y0);
+            while k >= 0 && y < x {
+                self.sieve[k as usize] ^= 1u16 << SieveOfAtkin2::SHIFTS[delta as usize];
+                (k, y) = (k - y - 15, y + 30);
+            }
+            (k0, x) = (k0 + x + 5, x + 10);
+        }
+    }
+    pub fn is_prime(&self, num: usize) -> bool {
+        if num < 2 {
+            return false;
+        }
+        if num == 2 || num == 3 || num == 5 {
+            return true;
+        }
+        let (num_div_60, num_mod_60) = (num / 60, num % 60);
+        if num_div_60 >= self.sieve.len() {
+            panic!("queried number is out of range of the sieve")
+        }
+        self.sieve[num_div_60] & (1u32 << SieveOfAtkin2::SHIFTS[num_mod_60]) as u16 != 0
+    }
+    pub fn iter(&self) -> impl Iterator<Item = i64> + '_ {
+        let mut num: usize = 1;
+        let mut offset = SieveOfAtkin2::OFFSETS.iter().cycle();
+        [2, 3, 5]
+            .into_iter()
+            .chain(
+                self.sieve
+                    .iter()
+                    .flat_map(|bitfield| (0..16).map(move |shift| bitfield >> shift & 1 == 1))
+                    .filter_map(move |is_prime| {
+                        let filtered = if is_prime { Some(num) } else { None };
+                        num += offset.next().unwrap();
+                        filtered
+                    }),
+            )
+            .take_while(|&num| num <= self.limit)
+            .map(|num| num as i64)
+    }
+}