nekedome / project-euler

solutions/problem27.py

"""
Project Euler Problem 27: Quadratic Primes
==========================================

.. module:: solutions.problem27
   :synopsis: My solution to problem #27.

The source code for this problem can be
`found here <https://bitbucket.org/nekedome/project-euler/src/master/solutions/problem27.py>`_.

Problem Statement
#################

Euler discovered the remarkable quadratic formula:

.. math::

    n^2 + n + 41

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

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

Considering quadratics of the form:

.. math::

    &n^2 + an + b, \\mbox{ where } |a| \\lt 1000 \\mbox{ and } |b| \\le 1000 \\\\
    & \\\\
    &\\mbox{where } |n| \\mbox{ is the modulus/absolute value of } n \\\\
    &\\mbox{e.g. } |11| = 11 \\mbox{ and } |-4| = 4

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

Solution Discussion
###################

If :math:`n^2 + an + b` is prime for :math:`0 \\le n \\le n^{\\prime}` where :math:`n^{\\prime} \\gt 40`, then we can
assume that the solution maintains :math:`n^2 + an + b` as prime for the two cases :math:`n = 0, 1`. First, we'll
consider these cases.

Case (:math:`n=0`):

.. math::

    &n^2 + an + b = b \\\\
    \\Rightarrow &b \\mbox{ must be a prime}

Case (:math:`n=1`):

.. math::

    &n^2 + an + b = 1 + a + b \\\\
    \\Rightarrow &1 + a + b \\mbox{ must be a prime} \\\\
    \\Rightarrow &a = p - b - 1 \\mbox{ for some primes } b, p

Now, this has set up a search space on :math:`a,b`. For each such :math:`a,b` we must determine the :math:`n^{\\prime}`
s.t. :math:`n^2 + an + b` is prime for :math:`0 \\le n \\le n^{\\prime}`.

The answer is simply the maximal value of :math:`n^{\\prime}` where :math:`|a| \\lt 1000 \\mbox{ and } |b| \\le 1000`.

.. note:: the work in this algorithm is dominated by repeatably checking whether :math:`n^2 + an + b` is prime for
          various values of :math:`a,b,n`. Many tuples in this search will result in repeated values for
          :math:`n^2 + an + b` and so this algorithm benefits heavily from applying memoization to cache these results,
          avoiding redundant calculations.

Solution Implementation
#######################

.. literalinclude:: ../../solutions/problem27.py
   :language: python
   :lines: 80-
"""

from itertools import product
from math import ceil, log
from typing import Callable

from lib.numbertheory import is_probably_prime
from lib.sequence import Primes
from lib.util import memoize


def find_n_prime(is_prime: Callable[[int], bool], a: int, b: int) -> int:
    """ Find :math:`n^{\\prime}` s.t. :math:`n^2 + an + b` is prime for :math:`0 \\le n \\le n^{\\prime}`

    :param is_prime: a primality testing function
    :param a: the parameter :math:`a`
    :param b: the parameter :math:`b`
    :return: :math:`n^{\\prime}`
    """

    n = 1  # we have already tested n=0 since 0^2 + a*0 + b = b, and b is a prime
    while is_prime(n ** 2 + a * n + b):
        n += 1
    return n


def solve():
    """ Compute the answer to Project Euler's problem #27 """

    range_limit = 1000

    maxsize = 2 ** int(ceil(log(2 * range_limit, 2)))  # round 2 * range_limit up to a power of two
    is_prime = memoize(is_probably_prime, maxsize=maxsize)  # memoization to avoid redundant calculations

    primes = list(Primes(upper_bound=range_limit))  # build a list of primes in the search space

    # Perform the search for the maximal n^{\\prime} given by an a, b in the search space
    max_n_prime = 0
    for b, p in product(primes, repeat=2):
        a = p - b - 1  # a is determined by b, p
        new_n_prime = find_n_prime(is_prime, a, b)
        if new_n_prime > max_n_prime:
            max_n_prime = new_n_prime
            best = a, b

    # Calculate the answer (product of the coefficients)
    a, b = best

The variable best does not seem to be defined for all execution paths.

    answer = a * b

    return answer


expected_answer = -59231