diff --git a/notebooks/problem0007.ipynb b/notebooks/problem0007.ipynb index 0087ddf..bb86af5 100644 --- a/notebooks/problem0007.ipynb +++ b/notebooks/problem0007.ipynb @@ -30,8 +30,7 @@ } ], "source": [ - "# Primes.unrank is 0-indexed, so the 10001st prime is at index 10000\n", - "Primes().unrank(10000)" + "nth_prime(10001)" ] }, { @@ -97,7 +96,7 @@ "\n", "The test proceeds by repeatedly squaring $a^u$ modulo $n$ until we reach $a^{n-1}$. If at any point in this squaring, we find that $a^{2^{i-1} u}$ does not equal 1 or -1, but the following squaring gives $a^{2^i u} = 1$, then we have found a non-trivial square root of 1 modulo $n$. Therefore, $n$ must be composite.\n", "\n", - "After $t$ squarings, we will reach $a^{n-1}$. If this equals anything other than 1, then by the contrapositive of Fermat's last theorem, this also tells us that $n$ is composite.\n", + "After $t$ squarings, we will reach $a^{n-1}$. If this equals anything other than 1, then by the contrapositive of Fermat's little theorem, this also tells us that $n$ is composite.\n", "\n", "If, after $t$ squarings, we have not found a non-trivial square root of 1 modulo $n$ and $a^{n-1} \\equiv 1 \\pmod{n}$, then either $n$ is prime, or is a base $a$ strong pseudoprime." ] @@ -153,6 +152,149 @@ " return True" ] }, + { + "cell_type": "markdown", + "id": "45c9ca5f-97f7-4091-8a6b-6851a84b39a2", + "metadata": {}, + "source": [ + "With this implementation, we can generate the primes by iterating over the positive integers and testing if each is prime." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "9ea628f0-79f0-447d-a444-6e66565a3620", + "metadata": {}, + "outputs": [], + "source": [ + "from itertools import count\n", + "\n", + "def primes(start=1):\n", + " for n in count(start):\n", + " if is_prime(n):\n", + " yield n" + ] + }, + { + "cell_type": "markdown", + "id": "f8e9fad3-8ef0-469e-b381-43315b915559", + "metadata": {}, + "source": [ + "Then to solve the problem, we can just count the primes until we reach 10001.\n", + "\n", + "## Rosser's theorem\n", + "This strategy works... but we can go even deeper.\n", + "\n", + "Running Miller-Rabin on every integer until we reach the 10001st prime isn't super efficient. To avoid this, we can take advantage of [Rosser's theorem](https://en.wikipedia.org/wiki/Rosser%27s_theorem), which states that $p_n > n \\log n$. In other words, we can start our primality testing at $10001 \\log 10001$ instead of 1." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "d7f5255b-d21e-4cfa-837d-9289c44c156a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "92114" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lower = ceil(10001 * log(10001))\n", + "lower" + ] + }, + { + "cell_type": "markdown", + "id": "48cbf855-fe03-4f7f-9c07-45a1c1c0295f", + "metadata": {}, + "source": [ + "## Prime-counting function\n", + "There's one issue preventing us from jumping ahead to 92114: we don't know how many primes are less than 92114, so when we find the next prime (spoiler: it's 92119), we won't know what number prime that is (another spoiler: it's the 8897th prime).\n", + "\n", + "To answer this question, we can use the [prime-counting function](https://en.wikipedia.org/wiki/Prime-counting_function). SageMath has the prime-counting function available directly as `prime_pi`, but let's implement it ourselves. As a first pass, we can use a sieve. Check out [problem 10](https://projecteuler.net/problem=10) for a sieve implementation - to avoid duplicating work, we'll use SageMath's sieve here." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "43b0b9df-d745-4890-a835-aa0c3e3bced8", + "metadata": {}, + "outputs": [], + "source": [ + "def prime_counting_function(n):\n", + " return len(prime_range(n+1))" + ] + }, + { + "cell_type": "markdown", + "id": "0c1e85af-74ba-48ee-9f5e-9aab35cc2a55", + "metadata": {}, + "source": [ + "This will tell us where to start our tally if we start checking for primes at 92114." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "1636b854-6d24-4373-987a-8a66e3b91fe4", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "8896" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "tally = prime_counting_function(lower)\n", + "tally" + ] + }, + { + "cell_type": "markdown", + "id": "13c91512-6a23-4569-bdea-2e88f4414922", + "metadata": {}, + "source": [ + "This runs considerably quicker than starting from 1." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "1590869a-5d9d-4662-b07b-18511eb25b79", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "104743" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "for (i, p) in enumerate(primes(start=lower), start=tally+1):\n", + " if i == 10001:\n", + " break\n", + "\n", + "p" + ] + }, { "cell_type": "markdown", "id": "ec384441", @@ -161,6 +303,7 @@ "## Relevant sequences\n", "* Prime numbers: [A000040](https://oeis.org/A000040)\n", "* Carmichael numbers: [A002997](https://oeis.org/A002997)\n", + "* Prime-counting function: [A000720](https://oeis.org/A000720)\n", "\n", "#### Copyright (C) 2025 filifa\n", "\n", @@ -170,7 +313,7 @@ ], "metadata": { "kernelspec": { - "display_name": "SageMath 9.5", + "display_name": "SageMath 10.7", "language": "sage", "name": "sagemath" }, @@ -184,7 +327,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.2" + "version": "3.12.5" } }, "nbformat": 4,