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/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABGQUlEQVR4nO3de5zM1ePH8fdY7LrtynUJSyIGRdmIipJb6ivkm6jIJUouSYQUXcglFENIvlKiG5Ui9Iv4olxSiFCE2ETsuO5i5/fH+doIWbszc2Y+83o+Hp+HdmbtvLN2vfeczznH5fP5fAIAAEDYy2Y7AAAAAPyDYgcAAOAQFDsAAACHoNgBAAA4BMUOAADAISh2AAAADkGxAwAAcAiKHQAAgENQ7AA4ns/nk9frFfuxA3A6ih0Axzt8+LDi4uJ0+PBh21EAIKAodgAAAA5BsQMAAHAIih0AAIBDUOwAAAAcgmIHAADgEBQ7AAAAh6DYAQAABNEnn0gNGkhHjvj/Y1PsADiWx+OR2+1WYmKi7SgAkG78eCk5Wcqb1/8f2+VjK3YADuf1ehUXF6fk5GTFxsbajgMggm3bJpUrJ02dKrVr5/+Pz4gdAABAkEycKF1xhXTffYH5+BQ7AACAIDh+XHrzTenhh6VcuQLzGhQ7AACAIHj/fenPP6UuXQL3GhQ7AACAIBg/Xqpf39xjFyjZA/ehAQAAIElr10rffCPNnh3Y12HEDgAAIMAmTJBKlJDuuiuwr0OxAwAACKDkZGnGDOmRR6TsAZ4rpdgBsKp06dJyuVznXV27dpUk+Xw+DRo0SMWLF1euXLlUt25dbdy40XJqAMi4t96SUlOljh0D/1oUOwBWrVq1Snv37k2/Fi5cKElq2bKlJGn48OEaNWqUxo0bp1WrVik+Pl7169fX4cOHbcYGgAzx+cw0bLNmUrFigX89ih0AqwoXLqz4+Pj0a+7cuSpbtqzq1Kkjn8+nMWPGaMCAAWrevLkqV66sadOm6dixY5oxY4bt6ABwSUuWSJs2SY8+GpzXo9gBCBmpqal6++231b59e7lcLm3fvl1JSUlq0KBB+vtER0erTp06Wr58+UU/TkpKirxe7zkXANgwfrxUoYJUt25wXo9iByBkzJkzR4cOHVK7/x2gmJSUJEkqWrToOe9XtGjR9OcuZOjQoYqLi0u/SpYsGbDMAHAxe/ea7U0efVRyuYLzmhQ7ACFjypQpaty4sYoXL37O466/fUf0+XznPXa2fv36KTk5Of3atWtXQPICwD954w0pZ07poYeC95psUAwgJPz6669atGiRPvroo/TH4uPjJZmRu2Jn3XW8b9++80bxzhYdHa3o6OjAhQWASzh1Spo0SWrdWsqfP3ivy4gdgJAwdepUFSlSRE2aNEl/rEyZMoqPj09fKSuZ+/CWLFmiWrVq2YgJABkyd660e3fwFk2cwYgdAOvS0tI0depUtW3bVtnP2r3T5XKpZ8+eGjJkiMqVK6dy5cppyJAhyp07t1q3bm0xMQD8swkTpBo1pOuvD+7rUuwAWLdo0SLt3LlT7du3P++5Pn366Pjx43rsscd08OBB1ahRQwsWLFC+fPksJAWAS9u2TVqwQPrPf4L/2i6fz+cL/ssCQPB4vV7FxcUpOTlZsbGxtuMAcLjevaWpU81UbK5cwX1t7rEDAADwk2PHpClTpA4dgl/qJIodAACA37zzjpScLD32mJ3Xp9gBAAD4gc8njRsn3X23VLq0nQwUOwAAAD9Ytkz64Qfp8cftZaDYAQAA+MG4cdI110j16tnLQLEDAADIot9+kz76SOraVcpmsV1R7AA4lsfjkdvtVmJiou0oABxu0iQpJkZq29ZuDvaxA+B47GMHIJBSU6VSpaQWLSSPx24WRuwAAACy4MMPpd9/N9OwtlHsAAAAsmDsWOn22yW323YSzooFAADItDVrpBUrpNmzbScxGLEDAADIJI/H3F931122kxgUOwAAgEw4cECaMUN69FEpe4jMgVLsAAAAMmHKFPNrhw52c5yNYgcAAHCZTp+Wxo+XWrWSChe2neYvFDsAAIDL9Nln0q+/2j0X9kIodgAAAJdp3DipRg2penXbSc4VIrf6AQAAhIfNm6WFC6W33rKd5HyM2AEAAFyG116TihaV/v1v20nOR7ED4Fgej0dut1uJiYm2owBwiIMHpWnTzBYn0dG205zP5fP5fLZDAEAgeb1excXFKTk5WbGxsbbjAAhjI0dKAwZIO3eaUbtQw4gdAABABpw6ZRZNtGoVmqVOotgBAABkyCefmC1OevSwneTimIoF4HhMxQLwhzp1pLQ0aelS20kuju1OAAAALmHdOunrr6X337ed5J8xFQsAAHAJr74qlSol3XOP7ST/jGIHAADwD/btk2bMkLp2lbKH+FwnxQ4AAOAfTJwoRUVJHTvaTnJpFDsAAICLSE2Vxo+XHnpIKlDAdppLo9gBAABcxPvvS0lJUvfutpNkDMUOAADgAnw+s2iifn3J7badJmNC/BZAAAAAO1aulFatkubOtZ0k4xixA+BYHo9HbrdbiYmJtqMACEOvviqVKyc1bmw7ScZx8gQAx+PkCQCXa/duqXRpafRoqVs322kyjhE7AACAvxk/XsqTR2rXznaSy0OxAwAAOMvx49KkSVL79lK+fLbTXB6KHQAAwFneeUf688/wmoI9g2IHAADwPz6fua/u7rulq66ynebysd0JAADA/yxYIP34ozRhgu0kmcOIHQAAwP+MGiVVry7dcovtJJnDiB0AAICkDRvMiN2MGZLLZTtN5jBiBwAAIDNaV6KEdO+9tpNkHsUOAABEvKQksxq2e3cpRw7baTKPYgcAACLe+PGm0HXqZDtJ1lDsAABARDt+3BS7Dh2k/Pltp8kaih0Ax/J4PHK73UpMTLQdBUAImz7dbEjco4ftJFnn8vl8PtshACCQvF6v4uLilJycrNjYWNtxAISQtDSpUiXJ7ZY+/NB2mqxjuxMAABCx5s+XNm+W3njDdhL/YCoWAABErFGjpBtvlGrVsp3EPxixAwAAEen776Uvv5RmzgzfDYn/jhE7AAAQkUaPlkqVklq0sJ3Efyh2AAAg4uzZY44O695dyu6g+UuKHQAAiDgejxQTI3XsaDuJf1HsAABARDl6VHr9dVPq4uJsp/Evih0A63777Tc98MADKliwoHLnzq2qVatqzZo16c+3a9dOLpfrnKtmzZoWEwMIZ2+9JR06ZKZhncZBs8oAwtHBgwdVu3Zt3XbbbZo3b56KFCmin3/+Wfn/dq5Po0aNNHXq1PS3c+bMGeSkAJzg9GmzaKJ5c6l0adtp/I9iB8CqYcOGqWTJkueUttIX+G4bHR2t+Pj4ICYD4EQffyxt3Sq9/bbtJIHBVCwAqz755BNVr15dLVu2VJEiRVStWjVNnjz5vPdbvHixihQpovLly6tTp07at2+fhbQAwpnPJ40YId16q9mU2Ik4KxaAVTExMZKkXr16qWXLlvr222/Vs2dPTZw4UQ899JAkadasWcqbN68SEhK0fft2DRw4UKdOndKaNWsUHR193sdMSUlRSkpK+tter1clS5bkrFggwi1bJt1yi/Tpp9Jdd9lOExgUOwBW5cyZU9WrV9fy5cvTH+vevbtWrVqlFStWXPD37N27VwkJCZo5c6aaN29+3vODBg3S4MGDz3ucYgdEtqZNzTTshg1SNofOWTr0fwtAuChWrJjcbvc5j1WsWFE7d+78x9+TkJCgrVu3XvD5fv36KTk5Of3atWuXXzMDCD+bN0uffCI99ZRzS53E4gkAltWuXVs//fTTOY9t2bJFCQkJF/09Bw4c0K5du1SsWLELPh8dHX3BKVoAkWvkSKlYMal1a9tJAsvBnRVAOHjiiSe0cuVKDRkyRNu2bdOMGTM0adIkde3aVZJ05MgR9e7dWytWrNCOHTu0ePFi3X333SpUqJCaNWtmOT2AcLB3rzR9utSjh+T0n/kodgCsSkxM1OzZs/Xuu++qcuXKeuGFFzRmzBi1adNGkhQVFaX169eradOmKl++vNq2bavy5ctrxYoVypcvn+X0AMLB2LFSzpxS5862kwQeiycAOJ7X61VcXByLJ4AIdPiwVKqU1L699MorttMEHiN2AADAsaZMkY4ckXr2tJ0kOCh2AADAkU6elEaNklq1kkqWtJ0mOCh2AADAkd57T9q1S+rd23aS4OEeOwCOxz12QOTx+aRq1aT4eGn+fNtpgod97AAAgOMsWiR9/31kLJg4GyN2AByPETsg8jRoIO3fL61ZI7lcttMEDyN2AADAUdaskRYulN59N7JKncTiCQAO5vF45Ha7lZiYaDsKgCB6+WWpbFnp3nttJwk+pmIBOB5TsUDk2LJFqlBBmjAhMk6a+DtG7AAAgGMMHy4VLSq1bWs7iR0UOwAA4Ai7d0tvvSX16iXFxNhOYwfFDgAAOMLo0VKePJE5BXsGxQ4AAIS9AwekiROlrl2lSL6VlmIHAADCnscjnT4tde9uO4ldFDsAABDWjh6VXntN6thRKlLEdhq7KHYAACCsvfGGdOiQ9OSTtpPYR7EDAABhKzVVGjlSat1aKl3adhr7KHYAACBszZhhtjnp29d2ktBAsQMAAGEpLU0aNkz617+kSpVspwkN2W0HAAAAyIyPP5Y2b5amTrWdJHQwYgfAsTwej9xutxITE21HAeBnPp80dKhUp45Us6btNKHD5fP5fLZDAEAgeb1excXFKTk5WbGRvHMp4CD/939SvXrS/PlSw4a204QOih0Ax6PYAc5Tv745bWLNGsnlsp0mdHCPHQAACCurV0uLFkmzZlHq/o577AAAQFh5+WXp6qulFi1sJwk9jNgBAICw8eOP0kcfSRMnSlFRttOEHkbsAABA2Bg6VLrySqltW9tJQhMjdgAAICz8/LM5aeLVV6WcOW2nCU2M2AEAgLDw8stS4cJShw62k4Quih0AAAh5O3dK06ZJvXtLuXLZThO6KHYAACDkDR8uxcZKXbrYThLaKHYAACCk7d0rvfGG9MQTUt68ttOENoodAAAIaa+8IsXESI8/bjtJ6KPYAXAsj8cjt9utxMRE21EAZNL+/dKECVK3blJcnO00oY+zYgE4HmfFAuHrmWekMWOkHTukQoVspwl9jNgBAICQdOiQNHas9OijlLqMotgBAICQNG6clJoqPfmk7SThg2IHAABCzuHD0ujRUqdOUny87TThg2IHAABCzuuvm3L31FO2k4QXih0AAAgpx46ZLU7atpVKlrSdJrxQ7AAAQEiZOFE6cEDq1892kvBDsQMAACHj2DFp2DAzWnfVVbbThB+KHQAACBkTJ5pNifv3t50kPFHsAABASDh+XBo+nNG6rKDYAQCAkDBxovTHH9KAAbaThC+/F7vDh6XWraXFi/39kQEAgFMdP27urXvoIUbrssLvxS5vXmnzZum55yROoQVgk8fjkdvtVmJiou0oAC5h0iRG6/zB5fP5v3598onUtKn05ZfS7bf7+6MDwOXxer2Ki4tTcnKyYmNjbccB8DfHj5tRukaNpKlTbacJbwG5x+7uu6UbbmDUDgAAXNrkyYzW+UtAip3LJQ0aJC1bZkbtAAAALuT4cenll6UHH5Suvtp2mvAXsFWxTZpIiYmM2gEAgIubPFnat4/ROn8JWLE7M2q3fLm0cGGgXgUAAISrEyfMaN0DDzBa5y8B3ceucWOpRg1G7QAAwPnOjNY984ztJM4R0GJ3ZtRu5Urpiy8C+UoAACCcnBmta9OG0Tp/CvjJEw0bSjVrMmoHAAD+8sYbUlISo3X+FvBi53JJgwdL334rzZsX6FcDEI5+++03PfDAAypYsKBy586tqlWras2aNenP+3w+DRo0SMWLF1euXLlUt25dbdy40WJiAFlx/Lg0ZIgZrStXznYaZwnKWbH160u1ajFqB+B8Bw8eVO3atZUjRw7NmzdPP/74o1555RXlz58//X2GDx+uUaNGady4cVq1apXi4+NVv359HT582F5wAJk2YYK5t+7ZZ20ncZ6AnDxxIYsWmYL36afSXXcF4xUBhIOnn35a//3vf7V06dILPu/z+VS8eHH17NlTffv2lSSlpKSoaNGiGjZsmDp37nzJ1+DkCSB0HDliTplo2tQsnoB/BWXETpLq1ZNuvtkspmDUDsAZn3zyiapXr66WLVuqSJEiqlatmiaf9d1++/btSkpKUoMGDdIfi46OVp06dbR8+fILfsyUlBR5vd5zLgChYexYKTlZGjjQdhJnClqxO3Ov3Zo15ixZAJCkX375RRMmTFC5cuX0xRdfqEuXLurevbveeustSVJSUpIkqWjRouf8vqJFi6Y/93dDhw5VXFxc+lWyZMnA/k8AyJDkZGnECKlTJ6lUKdtpnCloxU6SbrtNqlvXzKmnpQXzlQGEqrS0NF1//fUaMmSIqlWrps6dO6tTp06aMGHCOe/ncrnOedvn85332Bn9+vVTcnJy+rVr166A5QeQcaNHm4UT/fvbTuJcQS12Lpf04ovSDz9I770XzFcGEKqKFSsmt9t9zmMVK1bUzp07JUnx8fGSdN7o3L59+84bxTsjOjpasbGx51wA7DpwQBo1SuraVSpe3HYa5wpqsZOk2rWlO+80K2RPnQr2qwMINbVr19ZPP/10zmNbtmxRQkKCJKlMmTKKj4/XwrPOJkxNTdWSJUtUq1atoGYFkHkjR5rZuv+tgUKABL3YSWbUbssW6X+30ACIYE888YRWrlypIUOGaNu2bZoxY4YmTZqkrl27SjJTsD179tSQIUM0e/ZsbdiwQe3atVPu3LnVunVry+kBZMTvv0uvvSb16CEVLmw7jbMFbbuTv2vZ0mxavGWLFB1tIwGAUDF37lz169dPW7duVZkyZdSrVy916tQp/Xmfz6fBgwdr4sSJOnjwoGrUqCGPx6PKlStn6OOz3Qlg1xNPSFOnStu3S1dcYTuNs1krdps2SZUrS6++Kj3+uI0EACIFxQ6w57ffpLJlzYIJNiQOPCtTsZJUsaL0wAPSSy9Jx47ZSgEAAALppZekPHmknj1tJ4kM1oqdZBZQ7N8vjRtnMwUAAAiEHTukN94wCyYYLA8Oq8Xuqqukjh2lYcPMpoUAAMA5XnjB3FP3v7VQCAKrxU6SnnnGTMWOHm07CQAA8JetW6Vp08y9dXny2E4TOawXuyuvlB57zGxaeOCA7TQAAMAfBg+W4uOlzp1tJ4ks1oudJD39tOTzScOH204CAACyasMGacYMMysXE2M7TWQJiWJXuLBZLTN2rLR3r+00AJzC4/HI7XYrMTHRdhQgogwYYO6j79DBdpLIY20fu787dEgqU8ZsgTJ2rO00AJyEfeyA4Fm+3Bwf+s47EofDBF9IjNhJUv78Up8+0sSJ0q+/2k4DAAAul88n9esnXXut1KqV7TSRKWSKnSR1726WRT//vO0kAADgcn3xhfT119LQoVK2kGoYkSNkpmLPePVVqVcv6ccfpWuusZ0GgBMwFQsEXlqaVL262drk668ll8t2osgUcn26c2epRAlz4yUAAAgP778vffedGa2j1NkTcsUuJsbsVP3hh9I339hOAwAALuXkSbO1SZMm0s03204T2UKu2ElSmzZS5cp/7W8HAABC19Sp0s8/S0OG2E6CkCx2UVHSyy9LixebGzEBAEBoOn7cnDLRurVZDQu7QrLYSdKdd0q33GJG7dLSbKcBAAAXMm6ctG+fKXewL2SLncslDRsmff+99O67ttMAAIC/O3TILJZ45BGpbFnbaSCFcLGTpJtuku65x9yQmZJiOw0AADjb8OHSiRPm32mEhpAudpK5EXPnTnMiBQAACA27d0ujR0tPPikVK2Y7Dc4I+WJXsaL08MNmCxSv13YaAOHE4/HI7XYrMTHRdhTAcZ59VsqXT3rqKdtJcLaQO3niQnbvlsqVM2fJcnMmgMvFyROAf61fL113nTR2rNS1q+00OFtYFDtJ6ttX8njMPjlFi9pOAyCcUOwA/7rzTmnbNmnjRilHDttpcLaQn4o94+mnzV+eF16wnQQAgMj15ZfSvHlmNSylLvSEzYidZFbfDBggbdokXX217TQAwgUjdoB/pKVJiYlSzpzS8uWcCRuKwmbETpK6dTPTsCyrBgAg+GbOlNauNQMtlLrQFFbFLlcu6fnnpVmzpG++sZ0GAIDIkZJiZs2aNjUnQyE0hVWxk6S2baUqVaTevaXwmUQGACC8jR8v7dplznJH6Aq7YhcVJY0cKS1bJn38se00AAA438GDZvFix45ShQq20+CfhF2xk6QGDaSGDc2+didP2k4DAICzvfyylJoqDRpkOwkuJSyLnWRu3Ny2jaPGAAAIpJ07pVdfNbdAxcfbToNLCavtTv6uQwczHfvzz1JcnO00AEIV250Amde2rfTFF9LWreYIMYS2sB2xk8wK2ePHzSaJAADAv9atk6ZPN1OwlLrwENbF7sorzdDwmDHSr7/aTgMg1Hg8HrndbiUmJtqOAoSlvn3NWe0dOthOgowK66lYSTpyxJxCUb+++akCAP6OqVjg8i1YYBYqfvSR1KyZ7TTIqLAvdpI0aZLUubO0erV0ww220wAINRQ74PKcPi1VrWruX1+6lFMmwklYT8We0b695HazaTEAAP4wZYq0YYM0ejSlLtw4othlzy6NGCEtXizNnWs7DQAA4cvrlQYOlB54QOL21PDjiGInSY0bS/XqSU89xabFAABk1tCh0uHD0pAhtpMgMxxT7FwuM2q3ZQubFgMAkBk7dpjp1969pZIlbadBZjhi8cTZOnSQ5swxGykWKGA7DYBQwOIJIGNatZK+/toMkuTNazsNMsMxI3ZnvPSSOc/u+edtJwEAIHwsXy7NmmX+HaXUhS/HjdhJ5rDigQOl9eulChVspwFgGyN2wD9LS5Nuusnco756tZTNccM+kcORn7qePc29Ab162U4CAEDomzlT+vZbadQoSl24c+SnLyZGGjlSmjfPXAAA4MKOHZOeflq65x6pbl3baZBVjix2kjn+pG5dM2rH9icAAFzY6NFSUpI0fLjtJPAHxxY7l8v8Zf3pJ2nCBNtpAAAIPXv2mH3rHn9cKlfOdhr4g2OLnWTOuevYURo0SDpwwHYaAMHm8XjkdruVyPb5wAX17Svlzi09+6ztJPAXR66KPdu+fdLVV0tt20pjx9pOA8AGVsUC51uxQqpVS5o0SerUyXYa+Ivji51kTqTo10/6/nupUiXbaQAEG8UOOFdamnTjjZLPZ1bDRkXZTgR/cfRU7Bndu0ulS5uFFM6vsQAA/LOpU6U1a6TXXqPUOU1EFLvoaOmVV6QFC6TPP7edBgAAe5KTpf79pTZtpNq1baeBv0VEsZOkf/1LqlfPjNqlptpOAwCAHc8/Lx09Kg0bZjsJAiFiit2Z7U+2bZPGjbOdBgCA4Nu0yUy/9u8vXXml7TQIhIhYPHG2rl2lt982+9vFx9tOAyAYWDwBmHvMGzeWtm6VNm40pzTBeSJmxO6MF16QcuQwx6cAsG/QoEFyuVznXPFn/dTVrl27856vWbOmxcRAeJo7V/riC3MeLKXOubLbDhBsBQpIQ4ZInTtLjzxi9vABYFelSpW0aNGi9Lej/rZMr1GjRpo6dWr62zlz5gxaNsAJUlKkJ56QGjQw95zDuSKu2ElShw5mQ8Zu3di/BwgF2bNnP2eU7u+io6P/8XkA/2z0aOnXX6VPPzX3nMO5Im4qVjJFbtw4ae1a6Y03bKcBsHXrVhUvXlxlypRRq1at9Msvv5zz/OLFi1WkSBGVL19enTp10r59+ywlBcLPnj3Siy+a82ArVrSdBoEWcYsnzta+vfTxx9KWLVLBgrbTAJFp3rx5OnbsmMqXL6/ff/9dL774ojZv3qyNGzeqYMGCmjVrlvLmzauEhARt375dAwcO1KlTp7RmzRpFR0df8GOmpKQoJSUl/W2v16uSJUuyeAIR6aGHpPnzzb91+fPbToNAi+hi9/vvUvnyZpPG8eNtpwEgSUePHlXZsmXVp08f9erV67zn9+7dq4SEBM2cOVPNmze/4McYNGiQBg8efN7jFDtEmuXLzSbEkydLHTvaToNgiMip2DOKFjUbNb7+upmWBWBfnjx5VKVKFW3duvWCzxcrVkwJCQkXfV6S+vXrp+Tk5PRr165dgYoLhKxTp6THHpOqV5cefth2GgRLRBc7yexrV6mSufcgLc12GgApKSnatGmTihUrdsHnDxw4oF27dl30eckstoiNjT3nAiLN+PHSDz9IEyawSDCSRHyxy55dGjtWWrHCbFwMILh69+6tJUuWaPv27frmm2907733yuv1qm3btjpy5Ih69+6tFStWaMeOHVq8eLHuvvtuFSpUSM2aNbMdHQhZe/dKAwearb2qV7edBsEU8cVOkurWle67T+rTR/J6bacBIsvu3bt1//3365prrlHz5s2VM2dOrVy5UgkJCYqKitL69evVtGlTlS9fXm3btlX58uW1YsUK5cuXz3Z0IGQ99ZSUM6f00ku2kyDYInrxxNl275YqVJA6dTL7/QBwDo4UQyRZvFi67TbpzTe5ty4SUezOMmKEOWpszRqpalXbaQD4C8UOkeLkSfPvV/780tKlUjbm5SIOn/Kz9OxpNm987DEWUgAAws+YMdLmzWbhBKUuMvFpP0uOHGb10IoVZggbAIBwsXu3NHiwOS7zuutsp4EtTMVeQLt25jy9n36SChWynQZAVjEVi0jQsqW0bJkZsYuLs50GtjBidwEjRkg+n9S3r+0kAABc2oIF0gcfSK+8QqmLdIzYXcTEiVKXLuann9q1bacBkBWM2MHJUlKkKlWkEiWkL7+UXC7biWATI3YX0amTdOONptydPGk7DQAAFzZihLR9uzRuHKUOFLuLypbNLKT48UfptddspwGQGR6PR263W4mJibajAAGxfbvZhLhXL8nttp0GoYCp2Evo3t2skN282QxzAwg/TMXCiXw+6c47pY0bzSBE3ry2EyEUMGJ3CS+8IOXLZ/a4AwAgVMyaJc2fL3k8lDr8hWJ3CXFx0qhR0ocfSnPn2k4DAIB08KDUo4fUooV099220yCUUOwyoFUrqWFDcyLF4cO20wAAIl3fvtKJE9wDjvNR7DLA5TILKQ4ckAYOtJ0GABDJli6VJk+Whg6Vihe3nQahhsUTl+GVV6SnnpJWrjRboQAIDyyegFOkpEjVqpnbhP77X86Dxfn4K3EZevSQqlaVHnmEve0AAME3fLi0das0aRKlDhfGX4vLkD27Gf5ev14aPdp2GgBAJNmyxexZ17u3OWkCuBCmYjPhySfNPXfr10tly9pOA+BSmIpFuPP5pHr1pF9/lTZskHLlsp0IoYoRu0wYPFgqUkR69FHzxQYAQCC99Zb01VfS669T6vDPKHaZkDevNH68tHCh9M47ttMAAJxs/34zU9SmjVS/vu00CHUUu0y6806zv90TT5gvOgAAAuHJJ6W0NLNZPnApFLssGDNGOnXKfNEBCD0ej0dut1uJiYm2owCZsnChmYYdMcLcAgRcCosnsmjKFKljR/PFd8cdttMAuBAWTyAcHTkiVa5sFuktWmQ2ywcuhRG7LGrfXqpTR+rcWTp2zHYaAIBT9O8v/fGH2WaLUoeMothlkctlNor87Tfp2WdtpwEAOMF//yuNG2f2rbvqKttpEE6YivWT4cOlfv3MF2PNmrbTADgbU7EIJydOmFOOChQw58JGRdlOhHDCiJ2f9Ool3XCDmZo9ccJ2GgBAuBo8WNq+3dzDTanD5aLY+Un27NKbb0rbtknPP287DQAgHK1da1bAPvusVLGi7TQIR0zF+tmLL0qDBknffGNG8ADYx1QswsHJk1Jiorl3+9tvpRw5bCdCOGLEzs/69pWuvVZ6+GEpNdV2GgBAuBg2zJwD++ablDpkHsXOz3LkMF+UmzaZ1UwAAFzKjz9KL7wg9ekjVatmOw3CGVOxAfLss9LQodLq1dJ119lOA0Q2pmIRyk6flmrXlpKTpe++k2JibCdCOGPELkCeeUaqUMFMyZ48aTsNACBUvfqquaduyhRKHbKOYhcgOXNKU6dKP/xg9rgDAODvNm82J0z06CHVqmU7DZyAqdgA69dPGjXKLGGvVMl2GiCyeDweeTwenT59Wlu2bGEqFiHl1CkzBXvokLRunZQrl+1EcAKKXYCdOGFuhM2XT1q+3Ox3ByC4uMcOoWjoUHPbDicWwZ+Yig2wmBizSnb1aumVV2ynAQCEgvXrpeeek556ilIH/2LELkj69DE3yK5aZfa5AxA8jNghlJw8KdWoYfY6XbNGio62nQhOwohdkDz/vFS+vPTgg1JKiu00AABbXnrJLKybNo1SB/+j2AVJTIw0fbrZuHjwYNtpAAA2rF1rit2AARw7icBgKjbIhgyRBg6Uli5laTsQLEzFIhSkpJgylyOHOU88Z07bieBEjNgFWZ8+0o03Sm3bSkeP2k4DAAiW556TtmyR3nqLUofAodgFWfbs5ot6zx5T8gAAzrdypTRihDRokFSliu00cDKmYi0ZP17q2lWaP19q2NB2GsDZmIqFTUeOSFWrSoUKScuWsZ8pAosRO0sefVSqX19q3146eNB2GgBAoDzxhJSUZBbQUeoQaBQ7S1wus3Hx0aPS44/bTgMACIQ5c6Q33pDGjJHKlbOdBpGAYmdRiRKSxyPNmCG9/77tNAAAf9q7V+rYUWraVOrQwXYaRAqKnWWtW0v33it16WK+CQDwH4/HI7fbrcTERNtREGF8PnOrTfbs0uTJZpYGCAYWT4SA/fvNKqkqVcxiimzUbcCvWDyBYPN4zG02n38uNW5sOw0iCRUiBBQqZLZAWbjQ3IcBAAhfmzZJvXubnQ8odQg2RuxCyJNPSuPGmR3Jq1a1nQZwDkbsECypqVLNmtLx49KaNVLu3LYTIdIwYhdChgyR3G7p/vulY8dspwEAXK7nnpPWr5feeYdSBzsodiEkOtqskP31VzN6BwAIH0uXSsOGSS+8IF1/ve00iFRMxYag1183GxjPmWOWyQPIGqZiEWjJydJ110mlSklffSVFRdlOhEjFiF0I6tz5r32P9uyxnQYA8E98PvN9++BBsxCOUgebKHYhyOUyO5XnzCk99JCUlmY7EQDgYt58U5o1y+xXV7q07TSIdBS7EFWokDRtmvTll9KoUbbTAIEzaNAguVyuc674+Pj0530+nwYNGqTixYsrV65cqlu3rjZu3GgxMfCXTZukbt3MCRP//rftNADFLqTVr2/2QurfX1q71nYaIHAqVaqkvXv3pl/r169Pf2748OEaNWqUxo0bp1WrVik+Pl7169fX4cOHLSYGpBMnpFatzCjdq6/aTgMYFLsQ99JL5kSK++6TvF7baYDAyJ49u+Lj49OvwoULSzKjdWPGjNGAAQPUvHlzVa5cWdOmTdOxY8c0Y8YMy6kR6Z56SvrpJ2nmTLY2Qeig2IW4nDnNvRu//25uzmUNM5xo69atKl68uMqUKaNWrVrpl19+kSRt375dSUlJatCgQfr7RkdHq06dOlq+fPlFP15KSoq8Xu85F+BPH39sNpR/5RXp2mttpwH+QrELA1dfLU2aZH4qfOMN22kA/6pRo4beeustffHFF5o8ebKSkpJUq1YtHThwQElJSZKkokWLnvN7ihYtmv7chQwdOlRxcXHpV8mSJQP6/4DIsnu31L69dM890mOP2U4DnIt97MJI585mKf033/ATIpzr6NGjKlu2rPr06aOaNWuqdu3a2rNnj4oVK5b+Pp06ddKuXbs0f/78C36MlJQUpaSkpL/t9XpVsmRJ9rFDlp0+Ld1+u/TLL9L330sFCthOBJyLEbswMmaMVK6cud/uyBHbaYDAyJMnj6pUqaKtW7emr479++jcvn37zhvFO1t0dLRiY2PPuQB/GDxYWrbMHBlGqUMootiFkVy5pPfek3btkrp2tZ0GCIyUlBRt2rRJxYoVU5kyZRQfH6+FCxemP5+amqolS5aoVq1aFlMiEn3xhfTii+bIsFtvtZ0GuDCKXZipUEGaMMFMyU6bZjsNkHW9e/fWkiVLtH37dn3zzTe699575fV61bZtW7lcLvXs2VNDhgzR7NmztWHDBrVr1065c+dW69atbUdHBNm9W3rgAalhQ+npp22nAS4uu+0AuHwPPij93/+Zm3YTEyW323YiIPN2796t+++/X/v371fhwoVVs2ZNrVy5UgkJCZKkPn366Pjx43rsscd08OBB1ahRQwsWLFC+fPksJ0ekOHnS7FcXEyNNny5lY0gEIYzFE2Hq6FFT6qKizGIK9lACLs7r9SouLo7FE8iUPn2k0aOlJUsk7gBAqOPnjjCVJ4+53+7nn6Xu3W2nAQBn+uQTacQIadgwSh3CA8UujFWuLHk80pQp5hBqAID/7NghtW0rNW0qPfGE7TRAxjAV6wCdOklvvy0tXy5Vq2Y7DRB6mIrF5UpJkW65Rdq/X1qzRrriCtuJgIxhxM4Bxo41CyhatJAOHrSdBgDCX+/eZgPi996j1CG8UOwcICZG+uAD6dAh6aGHpLQ024kAIHzNmmXOgR01Sqpe3XYa4PJQ7ByiTBkzHTt3rjR0qO00ABCeNmyQOnQw25twDizCEcXOQe68Uxo40FxnbdQPRCyPxyO3263ExETbURAGDh2SmjWTrrpKeuMNyeWynQi4fCyecJjTp03BW7vWXCVL2k4E2MfiCVxKWppZ/bpsmbR6tVS2rO1EQOYwYucwUVHmcOpcuaSWLc3KLgDAP3vxRemzz8z3T0odwhnFzoEKFTKLKb77TurWTWJMFgAu7rPPpEGDzHXnnbbTAFlDsXOoG2+UXn9dmjxZmjjRdhoACE3btklt2kh33SU984ztNEDWZbcdAIHz8MPmPrtu3aRKlcxmmwAA4+hRs1iiSBFp+nQpG0MdcAD+GjvcqFHSzTdL994r7dplOw0AhAafz2xrsn27NHu2FBdnOxHgHxQ7h8uRw+ycHhNjfjI9ftx2IgCwb9QosxHx1KlmRgNwCopdBChcWJozR/rxR+mRR1hMASCyzZsn9eljrpYtbacB/ItiFyGqVZOmTDGnU4wZYzsNANixaZM5VaJxY2nIENtpAP9j8UQEuf9+ad06c7h1lSrSHXfYTgQAwfPnn9K//iWVKCHNmGH2/QSchhG7CDNkiFS/vnTffdIvv9hOAwDBcfKk9O9/m3L3yScSB5DAqSh2ESYqSnr3XalAAenuu6XkZNuJACDwevWSliwxm7dzsgScjGIXga64Qpo7V9qzx/wEe+qU7URAYHg8HrndbiUmJtqOAosmTZLGjZNee0267TbbaYDAcvl8rJGMVF9+KTVqZFbKjhsnuVy2EwGB4fV6FRcXp+TkZMUyBxdRliwx9xM/8ojk8dhOAwQeI3YRrF49afx4c40dazsNAPjX9u1Sixbm1B12A0CkYFVshOvUSfrpJ+mJJ6Srr+YAbADOcPCg1KSJlD+/9P77ZrN2IBIwYgcNG2YOwG7VSlq/3nYaAMia1FQzUpeUJH32mVSwoO1EQPBQ7KCoKOmdd8xKsbvukn7/3XYiAMgcn8/cT7dsmTkD9pprbCcCgotiB0lS3rzSp5+avZ6aNuVMWQDh6aWXpGnTpDfflOrUsZ0GCD6KHdKVKGE27vzhB+nhh6W0NNuJACDjZsyQBg6UBg+WHnjAdhrADoodzlG9ujR9uvTee1K/frbTAEDGLF1qfiBt29aUOyBSUexwnhYtpFGjpOHDzf52ABDKtmyR7rlHql3bbEbMnpyIZGx3ggvq2VPauVPq3l268kqpWTPbiQDgfH/8YbZpKlJE+vBDKWdO24kAuyh2uKiRI6XffpNatzanVNSqZTsRAPzlyBGzV93hw9LKlea4RCDSMRWLi8qWzawuu/FG6e67zUbGABAKUlOle++VNm+W5s+XypSxnQgIDRQ7/KOYGGnOHCk+3pwrm5RkOxGASJeWJnXoIP3f/5nvT9Wq2U4EhA6KHS7piiukefPMT8hnpj2AcODxeOR2u5WYmGg7CvyoTx+zqfrbb0u33247DRBaXD6fz2c7BMLD999Lt95qtkT57DMzmgeEA6/Xq7i4OCUnJys2NtZ2HGTByJHSU09Jr70mdetmOw0QehixQ4Zdd505nWL5crOg4tQp24kARJLp002p69+fUgdcDMUOl+XWW6X33zcnVDzyiDmXEQACbd48qX17c734ou00QOii2OGy3XWXWS07darUuzflDkBgffONWQHbuLE0cSIbEAP/hH3skClt2kh//mk2MC5Y0EyNAIC/rV9vCl3VqtLMmVJ2/tUC/hFfIsi0bt1MuRswQCpQQOrSxXYiAE6ydatUv76UkGAWbOXObTsREPoodsiSZ5+VDhyQHntMyp9fatXKdiIATrBrl3THHWa7pS++MN9fAFwaxQ5Z4nJJY8ZIBw9KDz4oxcaacxsBILN+/92UumzZpIULzTmwADKGxRPIsmzZpDffNIWueXPzjRgAMuPPP8306+HD0qJFUokSthMB4YViB7/IkUN67z2zC/y//iV99ZXtRADCjddrfkDcs8f8gFi2rO1EQPih2MFvoqOljz6SbrnFbImydKntRADCxeHDZvXr5s3mnrpKlWwnAsITxQ5+FRNjDuWuUcP85L1ihe1EAELd0aPmHOoNG6QFC6QbbrCdCAhfFDv4Xe7c5uixatWkRo2kb7+1nQiRyuPxyO12KzEx0XYUXMSxY2aE/7vvpPnzpRtvtJ0ICG8un49zAxAYhw9LDRtKmzZJX34pXX+97USIVF6vV3FxcUpOTlZsbKztOPif48elu++WVq40pe7mm20nAsIfI3YImHz5zPmO5cubVW7ff287EYBQceKEdM895naNzz+n1AH+QrFDQMXFmRuhS5c2+1JR7gAcP25K3dKl5kSJW2+1nQhwDoodAi5/frN1QUKCdNtt0urVthMBsOXMQomlS6W5c6W6dW0nApyFYoegKFDAbDZ6zTVSvXqslgUikddrFlStXm1G8m+/3XYiwHkodgia/PnNVgbXXWfuuVuyxHYiAMFy8KD5ul+/3ozgc08dEBgUOwTVmQUVNWqYzUgXLbKdCECgHThgRuq3bZP+7//M1z+AwKDYIejy5DH31tSpY/av+vxz24kQKoYOHSqXy6WePXumP9auXTu5XK5zrpo1a9oLicuyb5+5t3b3bmnxYrY9AgKNYgcrcuUyJ1Q0bGhWx33wge1EsG3VqlWaNGmSrr322vOea9Sokfbu3Zt+fc5PA2Fh506z4nX/fnPrRZUqthMBzkexgzXR0abQtWgh/fvf0qRJthPBliNHjqhNmzaaPHmyrrjiivOej46OVnx8fPpVoEABCylxOTZtkmrXllJSpK+/lipWtJ0IiAwUO1iVI4f0zjvSY49JnTtLQ4ZInIUSebp27aomTZrojjvuuODzixcvVpEiRVS+fHl16tRJ+/btC3JCXI5Vq6RbbjELpv77X+nqq20nAiJHdtsBgGzZpLFjpcKFpQEDzLTNyJHmcTjfzJkztXbtWq1ateqCzzdu3FgtW7ZUQkKCtm/froEDB+r222/XmjVrFB0dfcHfk5KSopSUlPS3vV5vQLLjfF9+aW6vqFLFbD58gQFYAAFEsUNIcLmk556TChWSunUz5W7KFDOiB+fatWuXevTooQULFigmJuaC73Pfffel/3flypVVvXp1JSQk6LPPPlPz5s0v+HuGDh2qwYMHByQzLu6jj6T77zf7033wgVkoBSC4XD4fE18ILTNnSg8+aDYyfe89s9ACzjRnzhw1a9ZMUVFR6Y+dPn1aLpdL2bJlU0pKyjnPnVGuXDl17NhRffv2veDHvdCIXcmSJZWcnKzY2Fj//49Ab7xhbqe47z7pP/+Rcua0nQiITIzYIeS0amWmb5o3lxo0kD7+2JxcAeepV6+e1q9ff85jDz/8sCpUqKC+fftesNQdOHBAu3btUrFixS76caOjoy86TQv/8vmkF14wI+5du0qvvcZtFIBNFDuEpIYNzb06d90l1apl9rq76irbqeBv+fLlU+XKlc95LE+ePCpYsKAqV66sI0eOaNCgQWrRooWKFSumHTt2qH///ipUqJCaNWtmKTXOOHlS6tJFevNN6aWXpH79zG0VAOzh5yqErJo1zZmyp05JN90kffut7UQItqioKK1fv15NmzZV+fLl1bZtW5UvX14rVqxQvnz5bMeLaF6v+cFr+nRz9e9PqQNCAffYIeT98YfUtKm0bp307rvmv4HL4fV6FRcXxz12fvLbb1KTJtKOHdLs2eZkCQChgRE7hLzChc207J13Ss2ama1RANixfr0ZTf/zT2nZMkodEGoodggLuXKZFbK9ekndu5tfT5+2nQqILPPnSzffbLYlWrlS+tvtkQBCAMUOYSNbNrNx8bhx0quvmilZ9p0FAs/nM19zTZqYEyW+/loqXtx2KgAXQrFD2Ona1exov3SpWVTxyy+2EwHOlZoqPfKI1LOnGSn/+GOJdStA6KLYISw1amSmglJTpRtvlBYvtp0IcJ79+81ektOmSVOnSiNGSBfYWhBACKHYIWxVrCh9841UtapUv740caLtRIBzbNwo1agh/fij9NVXUrt2thMByAiKHcJagQLSvHlmk9QuXaRHH5XOOkkKQCZ88IFZ+Zonj9k/snZt24kAZBTFDmEvRw6zBcqkSWYH/Dp1pN27badCKPB4PHK73UpMTLQdJSycOiX16SO1bGm2F1q+XCpd2nYqAJeDDYrhKKtWSS1aSCdOSDNnSrffbjsRQgEbFF/aH3+Yc5qXLJGGD5eeeIKTJIBwxIgdHCUxUVq7VrruOnPf3fDhZqsGABe3apV0ww3Shg3SokVm9SulDghPFDs4TqFCZiPVvn3Nde+97HcHXIjPJ02ebDYdLl5cWrNGqlvXdioAWUGxgyNFRUlDhkhz5pgRiMRE6YcfbKcCQsfhw9IDD5g96tq3N1OwJUrYTgUgqyh2cLSmTaXVq6WYGLPf3fjxTM0Ca9dK118vffqp9O670oQJUnS07VQA/IFiB8crV87sd9exozm1okULc4A5EGl8PnMk3003SbGxpuC1amU7FQB/otghIsTEmH/QPvrInFJRtaq0bJntVEDwHDxofqjp1s3s+bh8uXT11bZTAfA3ih0iSrNm0rp1UkKC2e/uxRel06dtpwIC66uvzErxxYul2bOlV19l6hVwKoodIk6pUuYfugEDpGeflW69Vdq2zXYqwP9SUqSnnpLq1ZPKljU/1Nxzj+1UAAKJYoeIlD279Pzz0tdfS0lJZjTj9ddZWAHnWL/eLBh67TWzn+OXX5ofagA4G8UOEe3mm6Xvv5cefNCcM9u4sfTbb7ZTAZl3+rQ0cqTZ4ictzWw+3Lu3lI3v9kBE4EsdES9vXjNa9/nnZq+7ypWlGTMYvUP42bhRqlXLnPf6+OOm1F17re1UAIKJYgf8T+PG5kilxo2lNm3MiRV79thOhazweDxyu91KTEy0HSWgTp6UXnhBqlbNbDy8fLkZtYuJsZ0MQLC5fD7GJYC/e/99M+KRkiKNGCF16MBUVjjzer2Ki4tTcnKyYmNjbcfxq+++kx5+2PxQ8vTT0sCBrHgFIhn/VAEX0LKltGmT1Ly5OXLpttukn36ynQr4y4kT0jPPmHvpJDPt+uKLlDog0lHsgIsoUEB6802zmvC338zK2ZdeklJTbSdDpJs/39wLOny4NGiQKXXVqtlOBSAUUOyAS7j9drN1RM+e0nPPmVMrFi60nQqRaNcuc3pE48Zmk+0ffjCjdjly2E4GIFRQ7IAMyJVLevllc7ZmoUJSgwZmmnbHDtvJEAlSU83oXMWK0ooV0rvvSosWSRUq2E4GINRQ7IDLcO210pIlZjuUb74x/9AOGiQdP247GZzqyy/NKHH//uZ+z82bpVatJJfLdjIAoYhiB1wml0u6/36zmKJnT2noUFPwPvqIve/gP5s2SXfdJd1xh7nfc+1aadQoyWGLegH4GcUOyKS8eU2p27DB3MjeooW5H2/VKtvJEM727ZMee0yqUkX68UfpvfekpUvZaBhAxlDsgCwqV06aO9dcf/xhzuf897+lrVttJ0M4OX7c3Md59dXmHrphw8yoXcuWTLsCyDiKHeAnTZqYc2enTjU3uLvdZuQlKcl2MoSykyelSZOka64xmwu3by9t2yY9+SR70gG4fBQ7wI+ioqR27aQtW6QhQ8zIy9VXS/36Sfv3206HUHL6tPTWW2Zla5cuUu3aZup1zBipYEHb6QCEK4odEAC5cklPPSX98os5mmzsWKl0aXM4+759ttPBprQ0adYsqVIlqW1bs/H199+bHwLKlbOdDkC4o9gBAXTFFea+qR07pB49pNdfNwXvySeZog0Gj8cjt9utxDPnbll08qT0n/+YKfpWraSyZaXVq81q6ipVbKcD4BQun48NGoBg+fNP6dVXzZWSInXoYLZMufpq28mczev1Ki4uTsnJyYoN8n4hx45JU6ZII0aYkyOaNjVT8zVqBDUGgAhBsQMsOHRIeu01M0V74ID5x75XL+nmm1kBGQg2il1SkjRhgrn+/NPsfdi3r9kaBwAChalYwIL8+aVnn5V27pQmTjSnCdx6q9kqZeZMM22H8LR6tfTgg1KpUtIrr5itb7ZskaZPp9QBCDyKHWBRrlxSp07Sxo3S559LcXFmZOeqq6TBg6Xdu20nREakpJhCXru2lJgoLVtm7q3cvVsaN858PgEgGCh2QAjIlk1q3Ngc7L5undSokbknKyHBTNN+9pnZHgOhZf16c4/klVeaQp4jhzR7ttmHrlcvMzILAMHEPXZAiPJ6pRkzzFTtunVSyZLSww9LbdpI5cvbThde/HmPnddrRuemTJG+/VYqXNhsW9Khg9mTDgBsotgBIc7nk9asMQXvvfdMsbjxRlPwWrWSihSxnTD0ZbXYnTghzZ9v9p/75BPzdqNGpszddZeUM2cAQgNAJlDsgDBy/Lg5k/add8w9eWlp0h13SM2bmynbokVtJwxNmSl2qalmanzWLGnOHFOor71Wuu8+6aGHpBIlApsZADKDYgeEqQMHpPffN8Xj66/NyF6tWlKzZubihv2/ZLTYJSdLX3xhyvNnn5ltSipUMCOj993HVCuA0EexAxxg/37p00/NjfsLFphVmpUrS/Xrm+vWW6U8eWyntOdixc7nkzZtMtOsc+dKS5dKp06Zkbm77zZblVSpwt6CAMIHxQ5wmCNHzKjTZ59JCxeaLTdy5DCjeXfcYa7q1aXs2W0nDZ4zxe7QoWTt3RurxYulr76SFi82Z/fGxEj16klNmpirVCnbiQEgcyh2gIP5fGZz3IULzfXVV9Lhw1Lu3GYBxk03/XUVKmQ7rf8dOSJ99520bJlX/fvHqUiRZO3bF6uoKLPf3G23SXXrmhM/cue2nRYAso5iB0SQU6fMFh3//a+0YoW0fLn0++/mudKlpapV/7quvdbso5ctTHa7/OMPM636/ffm9IfVq6VNm3zy+VyS9kkqqp49k9WwYaxuvlnKm9d2YgDwP4odEMF8PmnHDlPy1q41++WtW2cWZkhSdLRUtqxUrpy5ypc3v5YqJRUrZk7OCJa0NFNCd+yQtm//69effjKFbv9+8345c0rXXWemm89cJUp4VbBgcM+KBQAbKHYAzuHzSXv2SD/8IG3d+te1ZYv066+mYJ2RP79UvLgpefHx5u38+c3RaPnzm1Gx7NnNPX5n/5o9u9lO5MQJs9DjxAlzHT8uHTxoStrfr337zPueUbCgGVEsX16qWFFyu82v5cqdv6+cPzcoBoBQRrEDkGEpKWaUbPduU/727v3r16Qks13IoUPmV6/XlMTLER0tXXGFud+vcGHz65mrSBEzXVy6tCl0+fJl/ONS7ABEighaFwcgq6KjzV5uGdnPLS3NjMCdOiWdPGmuM/996pT5WDExf/2aM2f43M8HAKGKYgcgILJli+y98wDABn4+BgAAcAiKHQAAgENQ7AAAAByCYgcAAOAQFDsAAACHoNgBAAA4BBsUA3A8n8+nw4cPK1++fHK5XLbjAEDAUOwAAAAcgqlYAAAAh6DYAQAAOATFDgAAwCEodgAAAA5BsQMAAHAIih0AAIBDUOwAAAAcgmIHAADgEBQ7AAAAh6DYAQAAOET2jLzTmXMWAQAAYEdGzrvOULE7fPiw4uLi/BIKAAAAly85OVmxsbH/+D4un8/nu9QHysyIXWJiolatWnVZvyergv2aXq9XJUuW1K5duy75B+1vkfDna+s1bX1eI+XP18brRtLn1NbrRsrnVIqMP19brxlJX6uZeU2/jdi5XK7L/gOOiooK+hebjdeUpNjY2Ij4f42U1zwj2J/XSPrzjZSv1Uj6842Uz6kUOX++fP8Nz9cM2OKJrl27BupDh9Rr2hIpf758Tp33mjZfN9gi6c83Uj6nUuT8+fI5Dc/XzNBULC7M6/UqLi4uQ3PeCB98Xp2Hz6nz8Dl1Jj6vWcd2J1kQHR2t5557TtHR0bajwI/4vDoPn1Pn4XPqTHxes44ROwAAAIdgxA4AAMAhKHYAAAAOQbEDAABwCIodAACAQ1DsAiAlJUVVq1aVy+XSunXrbMdBJu3YsUMdOnRQmTJllCtXLpUtW1bPPfecUlNTbUfDZRg/frzKlCmjmJgY3XDDDVq6dKntSMiCoUOHKjExUfny5VORIkV0zz336KeffrIdC340dOhQuVwu9ezZ03aUsESxC4A+ffqoePHitmMgizZv3qy0tDRNnDhRGzdu1OjRo/X666+rf//+tqMhg2bNmqWePXtqwIAB+u6773TLLbeocePG2rlzp+1oyKQlS5aoa9euWrlypRYuXKhTp06pQYMGOnr0qO1o8INVq1Zp0qRJuvbaa21HCVtsd+Jn8+bNU69evfThhx+qUqVK+u6771S1alXbseAnI0aM0IQJE/TLL7/YjoIMqFGjhq6//npNmDAh/bGKFSvqnnvu0dChQy0mg7/88ccfKlKkiJYsWaJbb73VdhxkwZEjR3T99ddr/PjxevHFF1W1alWNGTPGdqyww4idH/3+++/q1KmTpk+frty5c9uOgwBITk5WgQIFbMdABqSmpmrNmjVq0KDBOY83aNBAy5cvt5QK/pacnCxJfF06QNeuXdWkSRPdcccdtqOEtey2AziFz+dTu3bt1KVLF1WvXl07duywHQl+9vPPP2vs2LF65ZVXbEdBBuzfv1+nT59W0aJFz3m8aNGiSkpKspQK/uTz+dSrVy/dfPPNqly5su04yIKZM2dq7dq1WrVqle0oYY8Ru0sYNGiQXC7XP16rV6/W2LFj5fV61a9fP9uRcQkZ/Zyebc+ePWrUqJFatmypjh07WkqOzHC5XOe87fP5znsM4enxxx/XDz/8oHfffdd2FGTBrl271KNHD7399tuKiYmxHSfscY/dJezfv1/79+//x/cpXbq0WrVqpU8//fScfzBOnz6tqKgotWnTRtOmTQt0VGRQRj+nZ77B7NmzR7fddptq1Kih//znP8qWjZ+HwkFqaqpy586t999/X82aNUt/vEePHlq3bp2WLFliMR2yqlu3bpozZ46+/vprlSlTxnYcZMGcOXPUrFkzRUVFpT92+vRpuVwuZcuWTSkpKec8h39GsfOTnTt3yuv1pr+9Z88eNWzYUB988IFq1KihEiVKWEyHzPrtt99022236YYbbtDbb7/NN5cwU6NGDd1www0aP358+mNut1tNmzZl8USY8vl86tatm2bPnq3FixerXLlytiMhiw4fPqxff/31nMcefvhhVahQQX379mWa/TJxj52flCpV6py38+bNK0kqW7YspS5M7dmzR3Xr1lWpUqU0cuRI/fHHH+nPxcfHW0yGjOrVq5cefPBBVa9eXTfddJMmTZqknTt3qkuXLrajIZO6du2qGTNm6OOPP1a+fPnS75eMi4tTrly5LKdDZuTLl++88pYnTx4VLFiQUpcJFDvgIhYsWKBt27Zp27Zt55VzBrrDw3333acDBw7o+eef1969e1W5cmV9/vnnSkhIsB0NmXRm65q6deue8/jUqVPVrl274AcCQgxTsQAAAA7BXeAAAAAOQbEDAABwCIodAACAQ1DsAAAAHIJiBwAA4BAUOwAAAIeg2AEAADgExQ4AAMAhKHYAAAAOQbEDAABwCIodAACAQ1DsAAAAHOL/AeT1fjyJEoHMAAAAAElFTkSuQmCC\n",
      "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)\n",
    "\n",
    "#### Copyright (C) 2025 filifa\n",
    "\n",
    "This work is licensed under the [Creative Commons Attribution-ShareAlike 4.0 International license](https://creativecommons.org/licenses/by-sa/4.0/) and the [BSD Zero Clause license](https://spdx.org/licenses/0BSD.html)."
   ]
  }
 ],
 "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
}