135 lines
29 KiB
Plaintext
135 lines
29 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "e4ea692f",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"# [Quadratic Primes](https://projecteuler.net/problem=27)\n",
|
||
|
|
"\n",
|
||
|
|
"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": 1,
|
||
|
|
"id": "a9751bc6",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"-59231\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"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\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"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",
|
||
|
|
"print(a * b)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "cbd83361",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"It's worth noting that the coefficients end up being $a=-61$ and $b=971$. 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": 2,
|
||
|
|
"id": "93c85dea",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABGUElEQVR4nO3de5zVc+LH8depNF00g+i2klCYcm8UuWRXYW1L0bK5lEtCKEnEtopVbiU0RS0Rm6x1v+S2q1gVyT3XyK+hEmGmm0l1fn98tpSKLjPzOZfX8/E4j9E505y3OU29z+fz/Xw+iWQymUSSJElpr1LsAJIkSSobFjtJkqQMYbGTJEnKEBY7SZKkDGGxkyRJyhAWO0mSpAxhsZMkScoQFjtJkqQMYbGTlPGSySQlJSW4H7ukTGexk5TxFi5cSF5eHgsXLowdRZLKlcVOkiQpQ1jsJEmSMoTFTpIkKUNY7CRJkjKExU6SJClDWOwkSZIyhMVOkiSpAj3+OLRrB4sXl/3XtthJyliFhYXk5+dTUFAQO4okrVZYCMXFULNm2X/tRNKt2CVluJKSEvLy8iguLiY3Nzd2HElZbOZMaNIE7r4bunQp+6/viJ0kSVIFuf122HZb+NOfyufrW+wkSZIqwNKlMGYMnHkmVK9ePs9hsZMkSaoA//wnfPstnHtu+T2HxU6SJKkCjBgRVsPutlv5PUeV8vvSkiRJApg+HV57DR59tHyfxxE7SZKkcjZyJDRsCMceW77PY7GTJEkqR99/D+PGwTnnQJVyniu12EmKaueddyaRSKxz69GjBwDJZJIBAwbQoEEDqlevTps2bZgxY0bk1JK08e65B378Ec4+u/yfy2InKapp06Yxd+7c1bfnn38egE6dOgFwww03MHToUIYPH860adOoV68ebdu2ZeHChTFjS9JGSSbDNGzHjlCvXvk/n8VOUlQ77LAD9erVW3178skn2XXXXTn88MNJJpMMGzaMK6+8ko4dO9K8eXPuuecelixZwrhx42JHl6Rf9eKL8NFHcP75FfN8FjtJKWPZsmXcd999nHnmmSQSCWbNmsW8efNo167d6s/Jycnh8MMPZ/LkyRv8OqWlpZSUlKx1k6QYRoyA/Hw47LCKeT6LnaSU8eijj/L999/TtWtXAObNmwdA3bp11/q8unXrrn5sfQYPHkxeXt7qW8OGDcstsyRtyJw5YXuT886DRKJintNiJyll3HnnnRxzzDE0aNBgrfsTP/sbMZlMrnPfmvr160dxcfHqW1FRUbnklaRfMno0VKsGp51Wcc/pBsWSUsL//d//8cILL/Dwww+vvq/e/640njdvHvXr1199//z589cZxVtTTk4OOTk55RdWkn7Fjz/CqFFw6qmQl1dxz+uInaSUMGbMGOrUqcOxa+ze2bhxY+rVq7d6pSyE6/AmTZrEwQcfHCOmJG2UJ54IU7HnnVexz+uInaToVq5cyZgxY+jSpQtV1ti9M5FI0KtXLwYNGkSTJk1o0qQJgwYNokaNGnTu3DliYkn6ZSNGwMEHwz77VOzzWuwkRffCCy8we/ZszjzzzHUe69u3L0uXLuX888/nu+++o2XLljz33HPUqlUrQlJJ+nUffQT//jfce2/FP3cimUwmK/5pJanilJSUkJeXR3FxMbm5ubHjSMpwF18M990HRUVh8URF8ho7SZKkMrJkCdx9N5x1VsWXOrDYSZIklZn774fiYujePc7zW+wkSZLKQDIJw4fD738PjRvHyWCxkyRJKgOTJ8Nbb8GFF8bLYLGTJEkqA7fdBk2aQNu28TJY7CRJkrbQnDnw0ENwwQVQKWK7sthJyliFhYXk5+dTUFAQO4qkDHfHHZCTA126xM3hPnaSMp772EkqT8uWwU47wQknQGFh3CyO2EmSJG2Bf/0LvvoKevSIncRiJ0mStEWGD4ff/Q7y82Mn8axYSZKkzTZ9OkyZAo8+GjtJ4IidJEnSZho+HBo1gj/8IXaSwGInSZK0Gb7+Ohwhdv75ULly7DSBxU6SJGkz3HknJBJw1lmxk/zEYidJkrSJli+HESOgc2eoXTt2mp9Y7CRJkjbRE09AUVE4aSKVWOwkSZI20W23QevWsN9+sZOsze1OJEmSNsGMGfDii2HhRKpxxE6SJGkTDB8O9etDx46xk6zLYicpYxUWFpKfn09BQUHsKJIyxPffw9ix0L07VK0aO826EslkMhk7hCSVp5KSEvLy8iguLiY3Nzd2HElpbNgw6NsXZs+GevVip1mXI3aSJEkbYcWKMA3bqVNqljpw8YQkSdJGeeop+PRTGDcudpINc8ROkiRpI9xyCxx0EBx4YOwkG+aInSRJ0q945x34z39g/PjYSX6ZI3aSJEm/4tZb4Te/Sc0tTtZksZMkSfoFX38N990Xjg/baqvYaX6ZxU6SJOkXjBoFlSpBt26xk/w6i50kSdIGLFsGhYVw6qlQu3bsNL/OYidJkrQB//oXzJ0LPXvGTrJxLHaSJEnrkUyGkyaOPBKaNYudZuO43YkkSdJ6TJ0K06bBk0/GTrLxHLGTlLEKCwvJz8+noKAgdhRJaeiWW6BJEzjmmNhJNl4imUwmY4eQpPJUUlJCXl4excXF5Obmxo4jKQ0UFUHjxnDzzXDhhbHTbDxH7CRJkn5mxAioWRO6do2dZNNY7CRJktawZAnccQecdRbUqhU7zaax2EmSJK3hvvvg++/DSRPpxmInSZL0P8lkWDRx3HGwyy6x02w6tzuRJEn6nxdegPffD6dNpCNH7CRJkv5n6FDYd184/PDYSTaPI3aSJEnAjBnwzDMwdiwkErHTbB5H7CRJkgh71jVoACedFDvJ5rPYSZKkrPfVV3DvvWEz4qpVY6fZfBY7SZKU9UaMgCpVoHv32Em2jMVOkiRltaVLQ7E780zYdtvYabaMxU5SxiosLCQ/P5+CgoLYUSSlsHvvhQULoFev2Em2XCKZTCZjh5Ck8lRSUkJeXh7FxcXk5ubGjiMphaxcCc2awZ57wsMPx06z5dzuRJIkZa0JE+DDD+Hvf4+dpGw4Yicp4zliJ2lDfvc7WLQIpk5N373r1uSInSRJykpvvQX/+Q888EBmlDpw8YQkScpSQ4fCTjtBx46xk5Qdi50kSco6X34J998PPXuG/esyhcVOkiRlneHDoXp1OPvs2EnKlsVOkiRllUWL4PbboVs3yLT1VBY7SdF9+eWXnHrqqdSuXZsaNWqw7777Mn369NWPd+3alUQisdatVatWERNLSmd33w0LF8JFF8VOUvYyaFZZUjr67rvvaN26NUcccQQTJkygTp06fPrpp2yzzTZrfd7RRx/NmDFjVv+6ajqf0i0pmhUrYNgwOPFEaNQodpqyZ7GTFNX1119Pw4YN1yptO++88zqfl5OTQ7169SowmaRM9Pjj8OmnMG5c7CTlw6lYSVE9/vjjtGjRgk6dOlGnTh32228/Ro8evc7nTZw4kTp16tC0aVO6devG/PnzI6SVlO6GDIFDDoEDD4ydpHx48oSkqKpVqwZA79696dSpE6+99hq9evXijjvu4PTTTwfggQceYOutt6ZRo0bMmjWL/v37s3z5cqZPn05OTs46X7O0tJTS0tLVvy4pKaFhw4aePCFluVdeCaXuscfgj3+MnaZ8WOwkRVW1alVatGjB5MmTV9930UUXMW3aNKZMmbLe3zN37lwaNWrE+PHj6bienUUHDBjAwIED17nfYidlt+OPh48+ghkzoFKGzllm6P+WpHRRv3598vPz17pvzz33ZPbs2b/4exo1asQnn3yy3sf79etHcXHx6ltRUVGZZpaUfj78MFxfd+mlmVvqwMUTkiJr3bo1H3300Vr3ffzxxzT6heVqCxYsoKioiPr166/38ZycnPVO0UrKXkOGQL16cMopsZOUrwzurJLSwcUXX8zUqVMZNGgQM2fOZNy4cYwaNYoePXoAsGjRIvr06cOUKVP4/PPPmThxIu3bt2f77benQ4cOkdNLSgdz58LYseH4sEx/z2exkxRVQUEBjzzyCPfffz/NmzfnmmuuYdiwYZzyv7fVlStX5t133+W4446jadOmdOnShaZ
|
||
|
|
"text/plain": [
|
||
|
|
"Graphics object consisting of 1 graphics primitive"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 2,
|
||
|
|
"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
|
||
|
|
}
|