Problem 27 "Quadratic primes"

Euler discovered the remarkable quadratic formula:

\[ n^2 + n + 41 \]

It turns out that the formula will produce 40 primes for the consecutive integer values \( 0 \le n \le 39 \). However, when \( n = 40, 40^2 + 40 + 41 = 40(40 + 1) + 41 \) is divisible by 41, and certainly when \(n = 41, 41^2 + 41 + 41\) is clearly divisible by 41.

The incredible formula \( n^2 - 79n + 1601 \) was discovered, which produces 80 primes for the consecutive values \( 0 \le n \le 79\). The product of the coefficients, −79 and 1601, is −126479.

Considering quadratics of the form:

\( n^2 + an + b \), where \( |a| < 1000 \) and \( |b| \le 1000 \)

where \( |n| \) is the modulus/absolute value of \( n \)
e.g. \( |11| = 11 \) and \( |-4| = 4 \)

Find the product of the coefficients, \( a \) and \( b \), for the quadratic expression that produces the maximum number of primes for consecutive values of \( n \), starting with \( n = 0 \).

問 27 「二次式素数」

オイラーは以下の二次式を考案している:

\[ n^2 + n + 41 \]

この式は, n を0から39までの連続する整数としたときに40個の素数を生成する. しかし, n = 40 のとき \( 40^2 + 40 + 41 = 40(40 + 1) + 41 \) となり41で割り切れる. また, n = 41 のときは \( 41^2 + 41 + 41 \) であり明らかに41で割り切れる.

計算機を用いて, 二次式 \( n^2 - 79n + 1601 \) という式が発見できた. これは n = 0 から 79 の連続する整数で80個の素数を生成する. 係数の積は, -79 × 1601 で -126479である.

さて, |a| < 1000, |b| ≤ 1000 として以下の二次式を考える (ここで |a| は絶対値): 例えば |11| = 11 |-4| = 4である.

\[ n^2 + an + b \]

n = 0 から始めて連続する整数で素数を生成したときに最長の長さとなる上の二次式の, 係数 a, b の積を答えよ.

struct Sieve {
    _sieve: Vec<bool>,
}

impl Sieve {
    fn rule_out(&mut self, prime: usize) {
        for i in (prime * prime..self._sieve.len()).step_by(prime) {
            self._sieve[i] = false;
        }
    }
    fn init(&mut self) {
        let sqrt = (self._sieve.len() as f64).sqrt() as usize;
        let mut index = 5usize;
        for &i in [2usize, 4].iter().cycle() {
            if index > sqrt {
                break;
            }
            if self._sieve[index] {
                self.rule_out(index);
            }
            index += i;
        }
    }
    fn new(below: u32) -> Self {
        assert!(below > 4);
        let sieve = vec![true; below as usize];
        let mut s = Self { _sieve: sieve };
        s.init();
        s
    }
    fn is_prime(&self, n: u32) -> bool {
        assert!(n < self._sieve.len() as u32);
        if n == 2 || n == 3 {
            return true;
        }
        if n == 0 || n == 1 || n % 2 == 0 || n % 3 == 0 {
            return false;
        }
        self._sieve[n as usize]
    }
    fn primes(&self, below: u32) -> Vec<u32> {
        let mut primes: Vec<u32> = vec![2u32, 3u32];
        let mut index = 5usize;
        for &i in [2usize, 4].iter().cycle() {
            if index >= below as usize {
                break;
            }
            if self._sieve[index] {
                primes.push(index as u32);
            }
            index += i;
        }
        primes
    }
}

fn quadratic_formula(n: u32, a: i32, b: u32) -> i32 {
    (n * n) as i32 + a * n as i32 + b as i32
}

fn main() {
    let (mut nmax, mut amax, mut bmax) = (1u32, 0i32, 0u32);
    let sieve = Sieve::new(2_000_000);
    for b in sieve.primes(1001) {
        for a in (-(b as i32) + 1)..=999 {
            let mut n = 1;
            let mut v = quadratic_formula(1, a, b);
            if !(v > 1 && sieve.is_prime(v as u32)) {
                continue;
            }
            while {
                n += 1;
                v = quadratic_formula(n, a, b);
                v > 1 && sieve.is_prime(v as u32)
            } {}
            if n > nmax {
                nmax = n;
                amax = a;
                bmax = b;
            }
        }
    }
    let product = amax * bmax as i32;

    println!("{}", product);
    assert_eq!(product, -59231);
}