eulerbooks/notebooks/problem0027.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e4ea692f",
   "metadata": {},
   "source": [
    "# [Quadratic Primes](https://projecteuler.net/problem=27)\n",
    "\n",
    "We can easily write a function to test how many consecutive primes $n^2 + an + b$ generates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a9751bc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import count\n",
    "\n",
    "def consecutive_primes(a, b):\n",
    "    f = lambda x: x^2 + a*x + b\n",
    "    for n in count(0):\n",
    "        if not is_prime(f(n)):\n",
    "            return n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce6a48ba",
   "metadata": {},
   "source": [
    "Couple of observations to speed up iterating through values of $a$ and $b$:\n",
    "1. $0^2 + 0a + b = b$ must be prime. Additionally, since $|b| \\leq 1000$, the largest possible value for $b$ is 997.\n",
    "2. $1^2 + 1a + b = 1 + a + b$ must be prime. Therefore $a = p - b - 1$ for some prime $p$. Furthermore, since $|a| < 1000$, it follows that $-1000 < p - b - 1 < 1000 \\implies b - 999 < p < b + 1001$, so $p$ must be less than $997 + 1001 = 1998$ in the most extreme case.\n",
    "\n",
    "These facts allows us to iterate over just prime numbers instead of the larger range of integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "03a9f808",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-61, 971)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coeffs = dict()\n",
    "primes = prime_range(2000)\n",
    "for b in primes:\n",
    "    if b > 1000:\n",
    "        break\n",
    "        \n",
    "    for p in primes:\n",
    "        a = p - b - 1\n",
    "        if not (abs(a) < 1000):\n",
    "            break\n",
    "            \n",
    "        coeffs[(a, b)] = consecutive_primes(a, b)\n",
    "\n",
    "a, b = max(coeffs, key=coeffs.get)\n",
    "a, b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbd83361",
   "metadata": {},
   "source": [
    "Their product is -59231.\n",
    "\n",
    "If you iterate through $n=0,1,2,\\ldots$ for the resulting quadratic and the other formula given, $n^2 - 79n + 1601$, you'll see that these polynomials don't actually generate any new prime numbers compared to Euler's formula - they just repeat primes that Euler's formula already gives.\n",
    "\n",
    "Furthermore, both of these polynomials are actually just shifts of Euler's formula:\n",
    "$$n^2 - 61n + 971 = (n-31)^2 + (n-31) + 41$$\n",
    "$$n^2 - 79n + 1601 = (n-40)^2 + (n-40) + 41$$\n",
    "\n",
    "What's going on? Let $f(n) = n^2 + n + 41$. An interesting property of [prime-generating functions](https://mathworld.wolfram.com/Prime-GeneratingPolynomial.html) is if $p(n)$ generates primes for $0 \\leq n \\leq u$, then $p(u - n)$ also will. To understand why, try evaluating $p(u-n)$ at $n=0,1,2,\\ldots,u$; you get $p(u),p(u-1),p(u-2),\\ldots,p(0)$, which are the same values given by $p(n)$, just in reverse.\n",
    "\n",
    "However, that property alone doesn't explain why shifting $f$ right generates *more* (non-distinct) primes than $f(n)$. To understand that aspect, let's look at a plot of $f(n)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "93c85dea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABGDUlEQVR4nO3deZzN1ePH8fc1MtYZpRiKaZJJ1yjFRKlI1lSitKhvlihRSLKkbMlYCsUlVFKytKqUwrcvLZYkW8oWMlkzMdc6GPf3x/mZ7GbMvffc+7mv5+Pxeci9zH03Y2bec87nnOPy+Xw+AQAAIOzlsR0AAAAA/kGxAwAAcAiKHQAAgENQ7AAAAByCYgcAAOAQFDsAAACHoNgBAAA4BMUOAADAISh2ABzP5/PJ6/WK/dgBOB3FDoDj7dmzR7GxsdqzZ4/tKAAQUBQ7AAAAh6DYAQAAOATFDgAAwCEodgAAAA5BsQMAAHAIih0AAIBDUOwAAACC6PPPpTp1pL17/f+2KXYAHMvj8cjtdis5Odl2FADIMmqU5PVKhQv7/227fGzFDsDhvF6vYmNjlZ6erpiYGNtxAESwP/6QrrxSGj9eatHC/2+fETsAAIAgGTNGuvBC6YEHAvP2KXYAAABBcPCg9PbbZqSuQIHAvAbFDgAAIAg+/FBKS5Patg3ca1DsAAAAgmD0aOn226XExMC9Rt7AvWkAAABI0rJl0vz50kcfBfZ1GLEDAAAIsNGjpVKlpLvvDuzrUOwAAAACyOuVJk6U2rSRLrggsK9FsQNg1eWXXy6Xy3XK1b59e0mSz+dTnz59VKpUKRUoUEA1a9bUypUrLacGgOybONGsiG3TJvCvRbEDYNWiRYu0devWrGvWrFmSpKZNm0qSBg8erKFDh2rkyJFatGiR4uLiVKdOHe3Zs8dmbADIFp/PnDRx993SpZcG/vUodgCsuuSSSxQXF5d1TZ8+XWXLllWNGjXk8/k0fPhw9ezZU02aNFFSUpImTJig/fv3a9KkSbajA8A5/fCDtHKl9OSTwXk9ih2AkHHo0CFNnDhRrVq1ksvl0oYNG7Rt2zbVrVs3689ER0erRo0amjdv3hnfTkZGhrxe7wkXANgwerQ5Quz224PzehQ7ACFj2rRp2r17t1r8/wGK27ZtkySVKFHihD9XokSJrOdOJyUlRbGxsVlX6dKlA5YZAM5kxw6zvUnbtlKeIDUuih2AkPHWW2+pQYMGKlWq1AmPu1yuE37v8/lOeex4PXr0UHp6etaVmpoakLwAcDZvv20K3f//rBoUbFAMICT8+eefmj17tj755JOsx+Li4iSZkbuSJUtmPb5jx45TRvGOFx0drejo6MCFBYBzyMw007APPigVKxa812XEDkBIGD9+vIoXL66GDRtmPZaQkKC4uLislbKSuQ9v7ty5uummm2zEBIBsmT5d2rRJeuqp4L4uI3YArDt69KjGjx+v5s2bK2/ef78suVwuderUSQMGDFC5cuVUrlw5DRgwQAULFlSzZs0sJgaAs/N4pKpVpSpVgvu6FDsA1s2ePVubNm1Sq1atTnmua9euOnDggNq1a6ddu3apatWqmjlzpooUKWIhKQCc2+rV0qxZ0rvvBv+1XT6fzxf8lwWA4PF6vYqNjVV6erpiYmJsxwHgcB07SpMnm6nY/PmD+9rcYwcAAOAne/dK77wjtW4d/FInUewAAAD8ZuJEU+7atrXz+hQ7AAAAP/D5zKKJu++WypSxk4FiBwAA4AfffSf9+qvUvr29DBQ7AAAAP/B4pKuuCt65sKdDsQMAAMilzZulTz4xo3VnOfEw4Ch2ABzL4/HI7XYrOTnZdhQADjd2rFSggNS8ud0c7GMHwPHYxw5AIB06JMXHS40bS6NG2c3CiB0AAEAufPKJtG2b3UUTx1DsAAAAcsHjkWrWlCpUsJ2Es2IBAADO27Jl0g8/SB99ZDuJwYgdAADAefJ4pEsvlRo1sp3EoNgBAACch127zBFibdtKeUNkDpRiBwAAcB7Gj5eOHJFat7ad5F8UOwAAgBzKzJRGjJAefFCKi7Od5l8hMnAIAAAQPqZPlzZulD74wHaSEzFiBwAAkEOvvSbdeKMUagfbMGIHAACQAytWSP/7nzRliu0kp2LEDgAAIAdGjJBKlZKaNLGd5FQUOwCO5fF45Ha7lRxqcyUAwlZamvTee1K7dtIFF9hOcyqXz+fz2Q4BAIHk9XoVGxur9PR0xcTE2I4DIIwNGiT17i2lpkqXXGI7zakYsQMAAMiGI0fMSRPNmoVmqZModgAAANkybZoZqevQwXaSM2MqFoDjMRULwB9uvVVyuaS5c20nOTO2OwEAADiHJUuk77+XPvrIdpKzYyoWAADgHF5/XSpTRmrUyHaSs6PYAQAAnMWOHdKkSVL79lLeEJ/rpNgBAACcxdixUlSU1Lq17STnRrEDAAA4g8OHpVGjpEcekS66yHaac6PYAQAAnMFHH0lbt4b2FifHo9gBAACcweuvS7VqSUlJtpNkT4jfAggAAGDHTz9JCxaYjYnDBSN2ABzL4/HI7XYrOTnZdhQAYej116WEBOnOO20nyT5OngDgeJw8ASCntmyRLr9cGjhQ6tzZdprsY8QOAADgJB6PFB0tPfaY7SQ5Q7EDAAA4zv790htvmH3rYmNtp8kZih0AAMBx3n1X2r07fLY4OR7FDgAA4P8dPSoNGyY1bmwWToQbtjsBAAD4fzNmSGvWSOPH205yfhixAwAA+H9Dh0o33CDdeKPtJOeHETsAAABJS5dK334rTZkiuVy205wfRuwAAAAkDR8ulSkj3Xuv7STnj2IHAAAi3tat0qRJ0tNPS3nDeD6TYgcAACLeqFFmQ+LWrW0nyR2KHQAAiGj790ujR0utWklFi9pOkzsUOwCO5fF45Ha7lZycbDsKgBD23nvSP/9IHTvaTpJ7Lp/P57MdAgACyev1KjY2Vunp6YqJibEdB0AIOXpUqlBBcruljz+2nSb3wvj2QAAAgNz5+mtp1Spp3DjbSfyDqVgAABCxhg2TkpOl6tVtJ/EPRuwAAEBEWr5cmj1bmjw5fDckPhkjdgAAICINGyZddll4b0h8MoodAACIONu2mQ2JO3SQLrjAdhr/odgBAICIM2qUKXRt2thO4l8UOwAAEFH27zfFzgkbEp+MYgfAus2bN+uRRx5RsWLFVLBgQVWqVEmLFy/Oer5FixZyuVwnXNWqVbOYGEA4Gz9e2rVLeuYZ20n8j1WxAKzatWuXqlevrttuu00zZsxQ8eLF9ccff6joST9G169fX+PHj8/6fb58+YKcFIATZGZKQ4dKTZtKCQm20/gfxQ6AVYMGDVLp0qVPKG2XX375KX8uOjpacXFxQUwGwIk++URav1764APbSQKDqVgAVn3++eeqUqWKmjZtquLFi+u6667TuNNsAT9nzhwVL15ciYmJatOmjXbs2GEhLYBw5vNJQ4ZIt90mVa5sO01gcFYsAKvy588vSercubOaNm2qn376SZ06ddKYMWP06KOPSpKmTp2qwoULKz4+Xhs2bNCLL76oI0eOaPHixYqOjj7lbWZkZCgjIyPr916vV6VLl+asWCDCzZ0r1awpffWV1KCB7TSBQbEDYFW+fPlUpUoVzZs3L+uxDh06aNGiRZo/f/5p/87WrVsVHx+vKVOmqEmTJqc836dPH/Xt2/eUxyl2QGS7807pzz/NiRNOOWniZEzFArCqZMmScrvdJzx29dVXa9OmTWf9O/Hx8Vq7du1pn+/Ro4fS09OzrtTUVL9mBhB+fvtN+vJLqUsX55Y6icUTACyrXr26Vq9efcJja9asUXx8/Bn/TlpamlJTU1WyZMnTPh8dHX3aKVoAkevVV6VSpaSHHrKdJLAYsQNg1TPPPKMFCxZowIABWrdunSZNmqSxY8eqffv2kqS9e/eqS5cumj9/vjZu3Kg5c+borrvu0sUXX6zGjRtbTg8gHGzdKk2cKHXsKDl9pySKHQCrkpOT9emnn2ry5MlKSkrSSy+9pOHDh+vhhx+WJEVFRWnFihVq1KiREhMT1bx5cyUmJmr+/PkqUqSI5fQAwsHrr0vR0dI
      "text/plain": [
       "Graphics object consisting of 1 graphics primitive"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(n) = n^2 + n + 41\n",
    "plot(f, (-5, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af51cf45",
   "metadata": {},
   "source": [
    "This quadratic is symmetric across $x=-0.5$. Because of this, $f(-1) = f(0)$, $f(-2) = f(1)$, and so on; in general, $f(n) = f(-1-n)$. However, since the axis of symmetry is just to the left of the $y$-axis, when we iterate through $n=0,1,2,\\ldots$, we only output values on one side of the axis, so we don't output any duplicated values.\n",
    "\n",
    "But when we transform this function by shifting it $k$ units to the right, $f(n - k)$, we cause the part of the function that's left of the axis of symmetry to *also* be output when we iterate through $n=0,1,2,\\ldots$. Since these values are all equal to values we were already outputting, they'll also be prime, but they won't be *new* primes. Since $f$ outputs primes from $n=0$ to $39$ (40 distinct values), we can shift it up to 40 units to the right and have every value from $n=0$ to $79$ be prime.\n",
    "\n",
    "Because of the symmetry of $f$, the transformation $f(u-n)$ is actually the *same* as shifting $f$ to the right $u+1$ units (this can be shown algebraically); then, after outputting the primes in reverse, the vertex is reached and the function starts increasing again, repeating the primes that were already output.\n",
    "\n",
    "## Relevant sequences\n",
    "* Primes generated by Euler's formula: [A005846](https://oeis.org/A005846)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.5",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}