Delete py/TruncatablePrimes.ipynb
This commit is contained in:
Peter Norvig authored and 
GitHub committed
1 parent
4578110de5
commit
276bb7246f
1 file changed
-780
@@ -1,780 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8812d495-6014-4669-97ae-003be4f54605",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"<div style=\"text-align:right\"><i>Peter Norvig<br>Jan 2026</i></div>\n",
|
||||
"\n",
|
||||
"# Truncatable Primes\n",
|
||||
"\n",
|
||||
"What's special about [this pencil](https://mathsgear.co.uk/products/truncatable-prime-pencil)?\n",
|
||||
"\n",
|
||||
"[](https://community.wolfram.com/groups/-/m/t/1569707)\n",
|
||||
"\n",
|
||||
"The 24-digit number printed on this pencil is a prime, and as you sharpen the pencil and remove digits one at a time from the left, the resulting numbers are all primes:\n",
|
||||
"\n",
|
||||
" 357686312646216567629137 is prime\n",
|
||||
" 57686312646216567629137 is prime\n",
|
||||
" 7686312646216567629137 is prime\n",
|
||||
" ...\n",
|
||||
" 137 is prime\n",
|
||||
" 37 is prime\n",
|
||||
" 7 is prime\n",
|
||||
"\n",
|
||||
"Numbers like this are called [**left-truncatable primes**](https://en.wikipedia.org/wiki/Truncatable_prime). \n",
|
||||
"\n",
|
||||
"My function `left_truncatable_primes` works as follows:\n",
|
||||
"- Start with the one-digit primes, (2, 3, 5, 7).\n",
|
||||
"- On each iteration create a list of new numbers that are one digit longer by placing each possible digit 1–9 to the left of each prime.\n",
|
||||
"- Out of those new numbers, find just the new primes.\n",
|
||||
"- If there are any new primes, the function recursively adds digits to them, giving us even longer primes.\n",
|
||||
"- In the end, the function returns all the left-truncatable primes."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "8ab95f58-0d60-4516-b19c-5f55614abb54",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from sympy import isprime # isprime(n) returns True if n is a prime\n",
|
||||
"\n",
|
||||
"def left_truncatable_primes(primes=(2, 3, 5, 7)) -> tuple[int]:\n",
|
||||
" \"\"\"All left-truncatable primes, in ascending order.\"\"\"\n",
|
||||
" new_numbers = (int(d + str(p)) for d in \"123456789\" for p in primes)\n",
|
||||
" new_primes = tuple(filter(isprime, new_numbers))\n",
|
||||
" if new_primes:\n",
|
||||
" return primes + left_truncatable_primes(new_primes)\n",
|
||||
" else:\n",
|
||||
" return primes"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "90d8f036-da26-4975-a0da-ab637694a679",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's see how many truncatable primes there are:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "4c9a569c-af7f-4335-b10b-785e5709862b",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"4260"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"LTP = left_truncatable_primes()\n",
|
||||
"\n",
|
||||
"len(LTP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a32b56c4-587b-48d5-91ec-1abcef2e909f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are **4260**; here are the smallest and largest 16 of them:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "e937655d-16eb-44ce-8698-0a0c09c9f6a9",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(2, 3, 5, 7, 13, 17, 23, 37, 43, 47, 53, 67, 73, 83, 97, 113)"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"LTP[:16]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "502be81c-fb57-431a-855a-21184a1aeb0d",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(66276812967623946997,\n",
|
||||
" 67986315421273233617,\n",
|
||||
" 86312646216567629137,\n",
|
||||
" 315396334245663786197,\n",
|
||||
" 367986315421273233617,\n",
|
||||
" 666276812967623946997,\n",
|
||||
" 686312646216567629137,\n",
|
||||
" 918918997653319693967,\n",
|
||||
" 5918918997653319693967,\n",
|
||||
" 6686312646216567629137,\n",
|
||||
" 7686312646216567629137,\n",
|
||||
" 9918918997653319693967,\n",
|
||||
" 57686312646216567629137,\n",
|
||||
" 95918918997653319693967,\n",
|
||||
" 96686312646216567629137,\n",
|
||||
" 357686312646216567629137)"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"LTP[-16:]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "e82700e8-8fb2-45cf-bcf0-7805e46174bb",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Right-Truncatable Primes\n",
|
||||
"\n",
|
||||
"What if you sharpen the pencil from the other end? You get `right-trunctable primes`:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"id": "3870bc47-5833-4879-b622-2094b4247239",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def right_truncatable_primes(primes=(2, 3, 5, 7)) -> tuple[int]:\n",
|
||||
" \"\"\"All right-truncatable primes, in ascending order.\"\"\"\n",
|
||||
" # Only consider (1, 3, 7, 9) as the digit on the right; placing any other digit forms a composite number\n",
|
||||
" new_numbers = (10 * p + d for p in primes for d in (1, 3, 7, 9))\n",
|
||||
" new_primes = tuple(filter(isprime, new_numbers))\n",
|
||||
" if new_primes:\n",
|
||||
" return primes + right_truncatable_primes(new_primes)\n",
|
||||
" else:\n",
|
||||
" return primes"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"id": "7b39ffb5-ba52-42e6-a270-79bc10b26c2c",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"83"
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"RTP = right_truncatable_primes()\n",
|
||||
"\n",
|
||||
"len(RTP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5323ae60-285d-49d6-a133-5c476eb4c5e5",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are only 83 right-truncatable primes. Here they are:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"id": "e3780212-2b86-469f-9df0-81bd6fcc9c94",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"(2, 3, 5, 7, 23, 29, 31, 37, 53, 59, 71, 73, 79, 233, 239, 293, 311, 313, 317, 373, 379, 593, 599, 719, 733, 739, 797, 2333, 2339, 2393, 2399, 2939, 3119, 3137, 3733, 3739, 3793, 3797, 5939, 7193, 7331, 7333, 7393, 23333, 23339, 23399, 23993, 29399, 31193, 31379, 37337, 37339, 37397, 59393, 59399, 71933, 73331, 73939, 233993, 239933, 293999, 373379, 373393, 593933, 593993, 719333, 739391, 739393, 739397, 739399, 2339933, 2399333, 2939999, 3733799, 5939333, 7393913, 7393931, 7393933, 23399339, 29399999, 37337999, 59393339, 73939133)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(RTP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "976c7f1f-79c5-48c3-8909-62b809993057",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can also find the **two-sided truncatable primes**, which are both right- and left-truncatable. It turns out there are only 15 of them:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"id": "b87ae26e-0dee-4c95-914f-c897fda17508",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{2, 3, 5, 7, 23, 37, 53, 73, 313, 317, 373, 797, 3137, 3797, 739397}"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"set(LTP) & set(RTP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "82fc3133-4a17-47ec-b742-7a37f78421a8",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There isn't an official name for a truncatable prime that remains prime for all possible choices of truncating either the rightmost or leftmost digit at each step. I'll call them **omni-truncatable primes**. I can compute them like this:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"id": "28f12184-b247-4cea-86d4-6b0e07e2a63f",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[2, 3, 5, 7, 23, 37, 53, 73, 373]"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def is_omni_truncatable_prime(n: int) -> bool:\n",
|
||||
" \"\"\"Is n a truncatable prime for all possible choices of truncating digits from left or right?\"\"\"\n",
|
||||
" return (isprime(n) and\n",
|
||||
" (n < 10 or (is_omni_truncatable_prime(int(str(n)[1:])) and\n",
|
||||
" is_omni_truncatable_prime(int(str(n)[:-1])))))\n",
|
||||
"\n",
|
||||
"[p for p in RTP if is_omni_truncatable_prime(p)]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "f2323985-2518-41d1-8a3b-6e21e40d27be",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are only 9 of them: 4 one-digit, 4 two-digit, and 1 three-digit.\n",
|
||||
"\n",
|
||||
"## Summary \n",
|
||||
"\n",
|
||||
"|Kind|Count|Largest|Digits|\n",
|
||||
"|---|---:|---:|---:|\n",
|
||||
"|Left-truncatable primes|4260|357686312646216567629137|24|\n",
|
||||
"|Right-truncatable primes|83|73939133 |8|\n",
|
||||
"|Two-sided primes|15|739397 |6|\n",
|
||||
"|Omni-truncatable primes|9|373 |3|\n",
|
||||
"\n",
|
||||
"# Primality Testing\n",
|
||||
"\n",
|
||||
"I was very impressed with the speed of the function [`sympy.isprime`](https://github.com/sympy/sympy/blob/master/sympy/ntheory/primetest.py), and I wanted to investigate the topic of primality testing. I'll start by defining `isprime_simple`, which follows the definition of a prime number almost verbatim:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"id": "1d460931-3d70-485d-807f-8fb360ef4ef2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def isprime_simple(n: int) -> bool:\n",
|
||||
" \"\"\"Is n prime? A prime number is defined as an integer greater than 1 \n",
|
||||
" that cannot be evenly divided by any whole number other than 1 and itself.\"\"\"\n",
|
||||
" divisors = range(2, n)\n",
|
||||
" return n > 1 and not any(n % d == 0 for d in divisors)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1834f8f0-3712-4302-ac56-79972f85d0db",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"To test a *d*-digit prime we have to iterate through nearly 10<sup>d</sup> divisors, which will be noticeably slow for anything larger than about 7-digit numbers.\n",
|
||||
"\n",
|
||||
"We can speed things up a bit with two ideas, one minor and one major. The minor idea is to cut the run time in half by only testing the odd divisors (after determining that 2 is a prime and all other even numbers are composite). The major idea is that we only need to test divisors up to √*n*, not up to *n*. If *n* is composite, then *n* is the product of two numbers, and one of them must be less than or equal to √*n*, so we can stop there."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "647fb66b-32ae-494a-8945-cf7358c2f958",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from math import sqrt\n",
|
||||
"\n",
|
||||
"def isprime_faster(n: int) -> bool:\n",
|
||||
" \"\"\"More sophisticated primality tester: test only odd divisors up to √n.\"\"\"\n",
|
||||
" if n <= 4: # Handle small numbers up to 4\n",
|
||||
" return n in (2, 3)\n",
|
||||
" elif n % 2 == 0: # Even numbers other than 2 are composite\n",
|
||||
" return False\n",
|
||||
" else: # Test odd divisors up to √n\n",
|
||||
" divisors = range(3, int(sqrt(n)) + 1, 2) \n",
|
||||
" return not any(n % d == 0 for d in divisors)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "05b8cbf7-ea10-4047-bcdb-a5ac7981384b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"For a 24-digit prime, `isprime_simple` would need to test septillions of divisors, while `isprime_faster` would \"only\" need to test trillions. That's the difference between never completing versus taking a few hours. But we would like to get the time down from hours to milliseconds. The major idea wasn't enough; we need a *massive* breakthrough. \n",
|
||||
"\n",
|
||||
"## Probabilistic Primality Testing\n",
|
||||
"\n",
|
||||
"**Fortunately**, [Pierre de Fermat](https://en.wikipedia.org/wiki/Pierre_de_Fermat) provided the massive breakthrough in 1640! \n",
|
||||
"\n",
|
||||
"[Fermat's little theorem](https://en.wikipedia.org/wiki/Fermat%27s_little_theorem) states that if *n* is prime and *a* is not divisible by *n*, then:\n",
|
||||
"\n",
|
||||
" *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*). \n",
|
||||
"\n",
|
||||
"We can use the theorem to create a [Fermat primality test](https://en.wikipedia.org/wiki/Fermat_primality_test): given an integer *n*, choose a random *a* and test if *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*). If the test is false then *n* is definitely composite. If the test is true, then we're not sure: *n* might be prime, might not. But if we choose multiple values of *a* and they all give a remainder of 1 (mod *n*), then that is stronger evidence that *n* is probably prime. This is called a [Monte Carlo algorithm](https://en.wikipedia.org/wiki/Monte_Carlo_algorithm); an algorithm that uses [randomization](https://en.wikipedia.org/wiki/Randomized_algorithm), and can sometimes be wrong: it can incorrectly report a **false prime** for a composite number *n*. \n",
|
||||
"\n",
|
||||
"Some examples of testing *n* by choosing random values of *a*:\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"|*n*|*a*|*a*<sup>(*n* - 1)</sup> (mod *n*)|Conclusion|\n",
|
||||
"|--:|--:|:--:|:--|\n",
|
||||
"|12|5|5|*12 is definitely composite*|\n",
|
||||
"|221|18|1|*221 could be prime or composite*|\n",
|
||||
"|221|2|16|*try again with new value for **a**; 221 is definitely composite*|\n",
|
||||
"|2137|2|1|*2137 could be prime or composite*|\n",
|
||||
"|2137|3|1|*try again: 2137 could still be prime or composite*|\n",
|
||||
"|2137|7|1|*try again: 2137 could still be prime or composite*|\n",
|
||||
"|2137|99|1|*try again: 2137 could still be prime or composite, but a lot of evidence that it is prime*|\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Here is a simple implementation:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"id": "c8a59c5e-a01c-49be-a931-c1547c4ac712",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def isprime_fermat(n: int, k=20) -> bool:\n",
|
||||
" \"\"\"If isprime_fermat(n) is True then n is probably prime; if False then n is definitely composite.\n",
|
||||
" Increase `k` to decrease the probability of false primes.\"\"\"\n",
|
||||
" if n <= 10:\n",
|
||||
" return n in (2, 3, 5, 7)\n",
|
||||
" else:\n",
|
||||
" return all(pow(a, n - 1, n) == 1 \n",
|
||||
" for a in sample(2, n, k))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6746894c-00e7-41df-b4a3-d5c85862cc94",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"(*Note* that `pow(a, n - 1, n)` computes *a*<sup>(*n* - 1)</sup> (mod *n*) very efficiently; see the last section of this notebook.)\n",
|
||||
"\n",
|
||||
"We need a way to sample *k* random integers from a range, one at a time:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"id": "476282dd-a48a-4eb5-8fed-8495d4878ec2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import random\n",
|
||||
"from typing import Iterable\n",
|
||||
"\n",
|
||||
"def sample(start: int, stop: int, k: int) -> Iterable[int]:\n",
|
||||
" \"\"\"Randomly sample `k` integers from range(start, stop), one by one.\"\"\"\n",
|
||||
" return (random.randrange(start, stop) for _ in range(k))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "15ef6884-253d-4a3a-b58f-397ac73f552d",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"How common are false primes? I define one function to estimate the probability of a false prime by testing a sample of integers, and another function to print a table of false prime probabilities, over various values of *k*, and various ranges of numbers *n*:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"id": "45a11a93-9ef3-48d1-80fa-e0113041346e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
" label k = 1 k = 5 k = 10 k = 20 \n",
|
||||
" 3 digits 1.431% 0.007% 0.000% 0.000%\n",
|
||||
" 4 digits 0.471% 0.027% 0.001% 0.000%\n",
|
||||
" 5 digits 0.120% 0.003% 0.003% 0.001%\n",
|
||||
" 6 digits 0.031% 0.002% 0.002% 0.000%\n",
|
||||
" 7 digits 0.011% 0.000% 0.000% 0.000%\n",
|
||||
" 8 digits 0.003% 0.000% 0.000% 0.000%\n",
|
||||
" 9 digits 0.001% 0.000% 0.000% 0.000%\n",
|
||||
"10 digits 0.000% 0.000% 0.000% 0.000%\n",
|
||||
"11 digits 0.000% 0.000% 0.000% 0.000%\n",
|
||||
"12 digits 0.000% 0.000% 0.000% 0.000%\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"from statistics import mean\n",
|
||||
"\n",
|
||||
"def false_prime_probability(samples: Iterable[int], k: int) -> float:\n",
|
||||
" \"\"\"The probability of false primes from isprime_fermat(n, k) for n across the samples.\"\"\"\n",
|
||||
" composites = (n for n in samples if not isprime(n))\n",
|
||||
" return mean(isprime_fermat(n, k) for n in composites)\n",
|
||||
" \n",
|
||||
"\n",
|
||||
"def report_false_primes(label_start_stop_triples: Iterable, k_values=(1, 5, 10, 20), num_samples=100_000):\n",
|
||||
" \"\"\"Report false prime probability for different digit ranges, and k values.\n",
|
||||
" First arg is a list of triples such as (\"3 digits\", 100, 1000).\"\"\"\n",
|
||||
" print(' label ', *[f'k = {k:<4}' for k in k_values])\n",
|
||||
" for (label, start, stop) in label_start_stop_triples:\n",
|
||||
" samples = list(sample(start, stop, num_samples))\n",
|
||||
" print(f'{label:10} ', *[f'{false_prime_probability(samples, k):8.3%}' for k in k_values])\n",
|
||||
"\n",
|
||||
"report_false_primes([(f'{d:2} digits', 10 ** (d - 1), 10 ** d) \n",
|
||||
" for d in range(3, 13)])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "f3fe25d6-5f0a-478b-90d0-6c89e092e92e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We see that with *k* = 1, 3-digit composites (such as 123 or 747) are falsely called prime about 1.4% of the time, and 4-digit comnposites about 0.5%, but once we get to 10-digit numbers (such as 1234567891), or increase *k* to 10 or more repetitions, false primes are very rare indeed.\n",
|
||||
"\n",
|
||||
"The so-called [**Carmichael numbers**](https://en.wikipedia.org/wiki/Carmichael_number) are composite numbers with a very high false prime probability. [Some Carmichael numbers](https://oeis.org/A033502) are worse than others. \n",
|
||||
"\n",
|
||||
"We can test some of them:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"id": "7871aaa7-af1c-462b-85bb-d46818ccabea",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
" label k = 1 k = 10 k = 30 k = 100 k = 1000\n",
|
||||
" 1105 69.560% 2.560% 0.000% 0.000% 0.000%\n",
|
||||
" 1729 75.350% 5.450% 0.010% 0.000% 0.000%\n",
|
||||
" 2465 72.460% 3.900% 0.020% 0.000% 0.000%\n",
|
||||
" 2821 76.310% 6.910% 0.090% 0.000% 0.000%\n",
|
||||
" 6601 80.270% 10.530% 0.130% 0.000% 0.000%\n",
|
||||
" 8911 79.820% 10.910% 0.150% 0.000% 0.000%\n",
|
||||
" 41041 70.160% 2.470% 0.000% 0.000% 0.000%\n",
|
||||
" 101101 71.590% 3.590% 0.010% 0.000% 0.000%\n",
|
||||
" 294409 95.080% 60.040% 22.570% 0.550% 0.000%\n",
|
||||
" 340561 81.990% 14.420% 0.160% 0.000% 0.000%\n",
|
||||
" 825265 60.780% 0.520% 0.000% 0.000% 0.000%\n",
|
||||
" 56052361 99.110% 91.720% 77.240% 41.920% 0.010%\n",
|
||||
" 118901521 99.350% 93.090% 81.470% 50.970% 0.150%\n",
|
||||
" 172947529 99.350% 94.220% 83.660% 53.970% 0.200%\n",
|
||||
" 216821881 99.470% 94.170% 84.640% 56.180% 0.430%\n",
|
||||
"1299963601 99.700% 97.030% 91.020% 73.300% 4.720%\n",
|
||||
"2301745249 99.790% 97.720% 92.210% 78.380% 8.150%\n",
|
||||
"9624742921 99.870% 98.340% 95.090% 85.480% 20.660%\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"carmichael_numbers = (1105, 1729, 2465, 2821, 6601, 8911, 41041, 101101, 294409, 340561, 825265, \n",
|
||||
" 56052361, 118901521, 172947529, 216821881, 1299963601, 2301745249, 9624742921)\n",
|
||||
"\n",
|
||||
"report_false_primes([(c, c, c + 1) \n",
|
||||
" for c in carmichael_numbers], \n",
|
||||
" k_values=(1, 10, 30, 100, 1000), num_samples=10_000)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4ecd2289-66d5-4b0f-ac00-5db0ad92dde9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This is bad; some Carmichael numbers have false prime probability above 99% with *k* = 1, and still over 50% with *k* = 100."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4c76de9d-4803-4db5-b6c5-811a6a1c6da1",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"So our Fermat test is mostly reliable, but a few composite numbers will consistently be identified as false primes.\n",
|
||||
"\n",
|
||||
"**Fortunately**, there are variations on Fermat's idea that always give the right answer. `sympy.isprime` ([source code here](https://github.com/sympy/sympy/blob/master/sympy/ntheory/primetest.py)) uses the [Miller-Rabin test](https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test) (sometimes called Rabin-Miller, but I give Gary Miller precedence because he is my former colleague). The algorithm is guaranteed to give the right answer, but it does get slow for very large numbers. `sympy.isprime` gains efficiency by breaking the range of *n* values into subranges and cleverly precomputing a list of *a* values that work for every *n* in a subrange. For numbers greater than 2<sup>64</sup>, `sympy.isprime` gives up on guaranteed correct answers and falls back on the [Baillie–PSW test](https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test), which is faster than Miller-Rabin, but is probabilistic. There are no known false positives for the Baillie-PSW test, but no proof that they don't exist.\n",
|
||||
"\n",
|
||||
"Below I test that all the primality test functions agree with `sympy.isprime` on the first 1,000 integers and a selection of larger primes and composites. (I don't use an `assert` because it is possible that `isprime_fermat` could fail, but a score of 100 indicates 100% success.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"id": "4f377075-e352-46f9-a716-f1de673a7086",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{'isprime_simple': 100, 'isprime_faster': 100, 'isprime_fermat': 100}"
|
||||
]
|
||||
},
|
||||
"execution_count": 16,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"integers = [*range(1000), 11939, 65537, 117223, 7629137, 11939*11939, 11939*65537]\n",
|
||||
"prime_fns = (isprime_simple, isprime_faster, isprime_fermat, isprime)\n",
|
||||
"\n",
|
||||
"{fn.__name__: 100 * mean(fn(n) == isprime(n) for n in integers) \n",
|
||||
" for fn in prime_fns if fn != isprime}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "bfd71dc7-7f6d-4c98-9886-5c242d9b297a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Speed of Primality testing\n",
|
||||
"\n",
|
||||
"**The great thing** about `isprime_fermat` and `sympy.isprime` is that they are **very fast**, handling even 32-digit primes in well under a millisecond. \n",
|
||||
"\n",
|
||||
"I'll make a chart of the speeds of all the primality checking functions running on primes with digit-lengths ranging from 5 to 32. For small *n* we repeat the test multiple time to get less noisy timing results; when a test takes over a second we skip the bigger primes."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"id": "8109d9fc-bbae-4ef0-a75c-908d9746a557",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import timeit\n",
|
||||
"import sympy\n",
|
||||
"\n",
|
||||
"big_primes = [sympy.nextprime(3 * 10 ** d) for d in range(4, 32)]\n",
|
||||
"\n",
|
||||
"def plot_run_times(primes, functions=prime_fns):\n",
|
||||
" \"\"\"For each primality-testing function, plot its run time on primes of different digit-lengths.\"\"\"\n",
|
||||
" plt.figure(figsize=(9, 6))\n",
|
||||
" for function in functions:\n",
|
||||
" plt.plot(*time_test(function, primes), 'o-', label=function.__name__)\n",
|
||||
" plt.grid(True, axis='y', which=\"major\", ls=':', color='gray')\n",
|
||||
" plt.xticks(range(len(str(min(primes))), len(str(max(primes))) + 1))\n",
|
||||
" plt.yscale('log')\n",
|
||||
" plt.xlabel('Number of digits in prime')\n",
|
||||
" plt.ylabel('Run time in seconds (log scale)')\n",
|
||||
" plt.legend()\n",
|
||||
"\n",
|
||||
"def time_test(function, primes, repeat=40) -> tuple[list[int], list[float]]:\n",
|
||||
" \"\"\"Time the function on each of the primes, repeated `repeat` times, stopping \n",
|
||||
" when one exceeds a second. Return two lists: digit sizes and corresponding run times.\"\"\"\n",
|
||||
" D, T = [], [] # D is length of primes in digits; T is time in seconds\n",
|
||||
" for p in primes:\n",
|
||||
" time = timeit.timeit(lambda: function(p), number=repeat) / repeat\n",
|
||||
" D.append(len(str(p)))\n",
|
||||
" T.append(time)\n",
|
||||
" if time > 1.0:\n",
|
||||
" break # Stop testing a function once it takes over a second of run time\n",
|
||||
" elif time > .01:\n",
|
||||
" repeat = 1 # Reset `repeat` to 1 if it took .01 second\n",
|
||||
" return D, T"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"id": "b1ab0f3d-6ce6-4b8c-9d4d-bb78bb1ac432",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAINCAYAAACTVPhqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAA0BJJREFUeJzs3XlclOX6+PHPMOyyiQqIIqLigvtaWq65Vpqay6lfLml1TDplfu3YrpkmbWanqDQrzU7l7jG3MpdcSEVQ00hNQVEEEVAWWWZ7fn+MTBIgDM4ww3C9v995cc8zN89z3TPYea65N5WiKApCCCGEEEIIUQYnWwcghBBCCCGEsF+SMAghhBBCCCHKJQmDEEIIIYQQolySMAghhBBCCCHKJQmDEEIIIYQQolySMAghhBBCCCHKJQmDEEIIIYQQolySMAghhBBCCCHK5WzrAOydwWDg8uXLeHt7o1KpbB2OEEIIIYQQd0xRFHJzcwkODsbJ6fZ9CJIwVODy5cuEhITYOgwhhBBCCCEs7uLFizRu3Pi2dSRhqIC3tzdgfDN9fHxsHI0QQgghhBB3Licnh5CQENO97u1IwlCB4mFIPj4+kjAIIYQQQgiHUpkh9zLpWQghhBBCCFEuSRiEEEIIIYQQ5ZKEQQghhBBCCFEumcMghBBCCGEDiqKg0+nQ6/W2DkU4KBcXF9Rq9R2fRxIGIYQQQohqptFoSE1NJT8/39ahCAemUqlo3LgxXl5ed3QeSRiEEEIIIaqRwWAgKSkJtVpNcHAwrq6usjmssDhFUbh69SqXLl0iPDz8jnoaakXCMGrUKPbs2cN9993H2rVrbR2OEEIIIWoxjUaDwWAgJCQET09PW4cjHFiDBg04f/48Wq32jhKGWjHp+dlnn+Xrr7+2dRhCCCGEECZOTrXiNkzYkKV6rmrFX2r//v0rtYudEEIIIYQQoiS7Txj27t3L8OHDCQ4ORqVSsXHjxlJ1PvnkE8LCwnB3d6dr167s27ev+gMVQgghhBDCAdl9wnDjxg06duzIxx9/XObrq1atYsaMGbzyyiscPXqU3r17M2zYMJKTk6s5UiGEEEKI6qU3KPx6LpP/HUvh13OZ6A2KVa/Xr18/ZsyYYbXznz9/HpVKxbFjx6x2jTtRnfE1bdqUxYsXW/06lWH3CcOwYcOYP38+o0ePLvP1RYsWMXXqVJ544gnatGnD4sWLCQkJ4dNPP63S9YqKisjJySnxANBqtQDodDp0Op3pWFlljUZjWlP572WDwWC6TnllRVFQFKVUGYwrK5RX1mg0AOj1+jLLOp2uRDukTdImaZO0SdokbZI22a5NBoMBRVGqXN52IpV7397FI58f5Lnvj/HI5we5J2oX206koiiKqX5xGbjj8tq1a3nzzTdLHb+Tdtxabty4MampqURERJSqY602lVcuK8ZGjRpx+fJl2rVrd8dtrahNxSzVjvL+PVWG3ScMt6PRaIiLi2Pw4MEljg8ePJiYmJgqnXPhwoX4+vqaHiEhIQDs2LEDgJ07d7Jz504Atm3bxv79+wHYuHEjsbGxAKxevZrjx48DsHLlSk6dOgXAsmXLSExMBCA6OpqUlBTAmPRkZGQAEBUVRW5uLhqNhqioKDQaDbm5uURFRQGQkZHBokWLAEhJSSE6OhqAxMREli1bBsCpU6dYuXIlAMePH2f16tUAxMbGmoZ07d+/n23btkmbpE3SJmmTtEnaJG2q5jYdPHgQgLy8PPLy8gC4fv06N27cAODatWum/RkyMzMpLCw0xVOctKyKOc30/8aTmm18rVhaTiHT/xvPthOppKWlmW4e09LSAONN45UrVwDjTWN6ejpgTJ6K35fCwkIyMzMByM/P59q1a4Bx1IeTkxPe3t7k5eWZvlS99QvWnJycKrfp0qVLGAwGgoKCyMzMNN3QpqWlmW6qrdGm69evmz6PitqUm5uLt7c3zs7OlWpTenq6KVm8cuWKWW0qTn7vpE35+fkUFBQApf89/frrr1SaUoMAyoYNG0zPU1JSFEA5cOBAiXoLFixQWrZsaXo+ePBgpX79+oqHh4fSqFEj5fDhw+Veo7CwUMnOzjY9Ll68qABKRkaGoiiKotVqFa1WqyiKomg0mjLLRUVFik6nK7Os1+tN1ymvbDAYFIPBUKqsKIqi1+vLLRcVFSmKoig6na7MslarVTQaTZllaZOxrNMblL2nUpUN8ReVmLMZyo38ghrfJkf8nKRN0iZpk7SpJrcpNzdXSUhIUG7cuKEYDAbT7+YVapQbRVolt6DotuWcAo3SY8EOJXT25jIfTWdvVu5asEO5fqNQySvUKHmFGiW3oEi5UaQts1zc7uL35Xblvn37Ks8995xiMBiUjz/+WGnRooXi5uamBAQEKA8//LCi1+tN9aZPn65Mnz5d8fX1Vfz9/ZWXX37ZdD6DwaCEhoYq8+bNUyZNmqT4+PgoEyZMUBITExVAiYuLUwwGg7J7924FULZt26Z06tRJcXd3V/r376+kpaUpmzdvVlq3bq14e3sr48ePN72fxeePiopSwsLCFHd3d6VDhw7KqlWrKmyfXq9XsrKylEcffVSpX7++4u7urrRo0UL54osvFEVRlHPnzimAcvToUUWv1yu7du1SAGXr1q2l4tu6daspvn/84x9Kbm6u6fMu7/0p/vstfn8++OADU1zXr19XnnjiCaVBgwaKt7e30r9/fyU+Pv62bcrPz1d+//13paCgoNS/p8zMTAVQsrOzlYo4xD4Mf18ySlGUEsd+/PHHSp/Lzc0NNze3UsddXFwAcHZ2LnXs72VXV9cKy7dew5yyk5NTueXi86vVatNau7eWb429vHJtbtP2k6m88UNCiW9rGvq6M2d4BEPbNayRbSqrXNM/J2mTtEnaJG2q6W0qfu7k5GS6XynSK7Sd8xOWoABpOUV0nPdzpeonzBuCp6uzKRaVSlVuuVhcXBzPPfccK1eupFevXmRlZbFv374SS8V+/fXXTJ06lUOHDnHkyBGeeuopmjZtypNPPmmq89577/Haa6/x6quvlrjGre8NwBtvvMHHH3+Mp6cn48aNY/z48bi5ufHtt9+Sl5fHqFGj+Oijj5g9ezYqlYpXXnmF9evX8+mnnxIeHs7evXuZOHEigYGB9O3b97Ztfe2110hISGDbtm3Ur1+fs2fPmr6lv7V9t8Y4b968CuOLjo5m9uzZlXp/bm17cfmBBx7A39+frVu34uvry5IlSxg0aBBnzpzB39+/ws/v7/+ebn1ekRqdMNSvXx+1Wm3qvimWnp5OYGDgHZ07Ojqa6OhoU3eQcGzbT6by9Dfx/H2qWFp2IU9/E8+nj3VhaLuGNolNCCGEsDfJycnUqVOHBx98EG9vb0JDQ+ncuXOJOiEhIXzwwQeoVCpatWrFiRMn+OCDD0okDAMGDGDWrFmm5+fPny/zevPnz+eee+4BYOrUqbz00kucO3eOZs2aATBmzBh2797N7NmzuXHjBosWLWLXrl307NkTgGbNmrF//36WLFlC3759K2xb586d6datG2CcfFwRc+Iz5/0ptnv3bk6cOEF6eropcX3vvffYuHEja9eu5amnnqowxjtRoxMGV1dXunbtyo4dOxg1apTp+I4dO3jooYfu6NyRkZFERkaSk5ODr6/vnYYq7JjeoPDGDwmlkgUwfkujAt74IYFBEUGonSyzAYoQQghxKw8XNQnzhlSq7uGkLCZ/FVthveWPd6dHmH+lrm2uQYMGERoaSrNmzRg6dChDhw5l1KhRJXauvvvuu0t8U96zZ0/ef/999Hq9qSem+Ka8Ih06dDCVAwMD8fT0NN2MFx87fPgwAAkJCRQWFjJo0KAS59BoNKWSmrI8/fTTPPzww8THxzN48GBGjhxJr169LBZfscq8P8Xi4uLIy8ujXr16JY4XFBRw7ty5Ctt0p+w+YcjLy+Ps2bOm50lJSRw7dgx/f3+aNGnCzJkzmTBhAt26daNnz54sXbqU5ORkpk2bZsOoRU1yOCmr1KSxWylAanYhh5Oy6Nm8Xrn1hBBCiKpSqVR4ulbutqx3eAMa+rqTll1Y5pddKiDI153e4Q2s9kWXt7c38fHx7Nmzh59++onXX3+duXPnEhsbi5+fX6XPU6dOnUrVu3UYmkqlKvG8+FjxykDFP7ds2UKjRo1K1Ctr2PnfDRs2jAsXLrBlyxZ+/vln7rvvPiIjI3nvvfcsEl9VGAwGGjZsyJ49e0q9Zs77XVV2nzAcOXKE/v37m57PnDkTgEmTJrF8+XLGjx9PZmYm8+bNIzU1lXbt2rF161ZCQ0Pv6LoyJKn2SM8tP1moSj0hhBDCmtROKuYMj+Dpb+JRQYmkoTg9mDM8wuq94s7OzgwcOJCBAwcyZ84c/Pz82LVrl2kp/OLVoIodPHiQ8PDwUt+eW1pERARubm4kJydXOPyoPA0aNGDy5MlMnjyZ3r1788ILL9w2YagKc96fLl26kJaWhrOzc6WGSFma3ScM/fr1K7EWbVmmT5/O9OnTLXpdGZJUewR4u1u0nhBCCGFtQ9s15NPHupRarCPolsU6rGnz5s0kJibSp08f6taty9atWzEYDLRq1cpU5+LFi8ycOZN//vOfxMfH89FHH/H+++9bNS4w9n7MmjWL559/HoPBwL333ktOTg4xMTF4eXkxadKk2/7+66+/TteuXWnbti1FRUVs3ryZNm3aWDxOc96fgQMH0rNnT0aOHMnbb79Nq1atuHz5Mlu3bmXkyJGVHtpVVXafMAhhbT3C/CvVtVuZcaBCCCFEdRnariGDIoI4nJRFem4hAd7G/62qjvl2fn5+rF+/nrlz51JYWEh4eDjfffcdbdu2NdWZOHEiBQUF9OjRA7Vazb/+9S+rT84t9uabbxIQEMDChQtJTEzEz8+PLl268PLLL1f4u66urrz00kucP38eDw8Pevfuzffff2/xGM15f1QqFVu3buWVV15hypQpXL16laCgIPr06XPHC/1Uhkqp6OvLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 900x600 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"plot_run_times(big_primes)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "159833b0-afaf-486f-9c8e-1c8d2a11d9bf",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"What have we learned about timing?\n",
|
||||
"- `isprime_simple` has a run time of *O*(*n*)\n",
|
||||
"- `isprime_faster` has a run time of *O*(√*n*); we see its slope is half as steep on the log-log graph\n",
|
||||
"- `isprime_fermat` and `isprime` ([according to Wikipedia](https://en.wikipedia.org/wiki/Fermat_primality_test)) have a run time of *O*(*k* log<sup>2</sup>*n*), where *k* is the number of repeats \n",
|
||||
"- `isprime_simple` takes more than a second for an 8-digit prime; `isprime_faster is 10,000 times faster for that prime\n",
|
||||
"- `isprime_faster` takes more than a second for a 17-digit prime; `isprime_fermat` is over 10,000 times faster for that prime\n",
|
||||
"- `isprime_fermat` and `isprime` are well under a millisecond run time all the way out to a 32-digit prime\n",
|
||||
"- `isprime_fermat` is about ten times slower than `isprime`\n",
|
||||
" - (But that's pretty good for `isprime_fermat` considering how much simpler it is: 8 lines versus about 800.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1347399a-24f2-4d5f-aa7f-822951126bf5",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This was fun for me! I remember being amazed by Fermat's Little Theorem when I learned about it in [Prof. Michael Rosen's](https://mathematics.brown.edu/people/michael-rosen) number theory class in college, but I never worked through the details of how often the test gives a false prime result until now. \n",
|
||||
"\n",
|
||||
"That number theory class was important for me in several ways. First, it was [hard fun](https://dailypapert.com/hard-fun/). Second, it was the first class that felt like actually *doing* math, rather than being taught some techniques to prepare you to do math some time in the future. And third, I learned a valuable lesson: there was a take-home final exam, handed out Friday and due Monday. I worked on it all day Saturday, and couldn't answer a single one of the ten problems. At that point, I thought to myself \"This never happens on my computer science assignments; I guess I'm more cut out for CS than math.\" On Sunday I came back to it and eventually got answers to most of the problems. (I remember that on one problem, Prof. Rosen said \"You had a very elegant solution for that one, better than the answer I expected. But then I realized it was so good, it could be extended to be a trivial solution to this known difficult problem, and so I looked more closely and noticed a flaw in your argument. Too bad.\")\n",
|
||||
"\n",
|
||||
"## Note on Modular Exponentiation\n",
|
||||
"\n",
|
||||
"Just one more thing: none of this would work unless we can efficiently compute *a*<sup>(*n* - 1)</sup> (mod *n*). How does the `pow` builtin function do it? When *a* and *n* are 24-digit numbers, if we naively tried to compute `a ** (n - 1)` and then apply (mod *n*), we'd have two problems: we'd need a billion petabytes to store the result, and we'd need centuries to compute it. The way around these problems is to use [modular exponentiation](https://en.wikipedia.org/wiki/Modular_exponentiation) where we apply the modulus to each intermediate result, and cut the exponent in half each iteration, so we need only do *O*(log *n*) multiplications, not *O*(*n*). That's a big difference: 10<sup>24</sup> is a trillion trillion, and log<sub>2</sub>(10<sup>24</sup>) is only 80. \n",
|
||||
"\n",
|
||||
"There are two key ides that make this work:\n",
|
||||
"1) *b*<sup>2*e*</sup> = (*b* × *b*)<sup>*e*</sup>\n",
|
||||
"2) *b*<sup>2</sup> (mod *m*) = (*b* (mod *m*))<sup>2</sup> (mod *m*)\n",
|
||||
" - Proof: let *r* = *b* (mod *m*); then *b* = (*km* + *r*) for some integer value of *k*, and:\n",
|
||||
" - *b*<sup>2</sup> (mod *m*) = (*km* + *r*)<sup>2</sup> (mod *m*) = ((*km*)<sup>2</sup> + 2*kmr* + *r*<sup>2</sup>) (mod *m*) = *r*<sup>2</sup> (mod *m*)\n",
|
||||
"\n",
|
||||
"Here is an implementation:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"id": "e71eb140-7e1a-4146-87a7-8c4eebf61486",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def modpow(b: int, e: int, m: int) -> int:\n",
|
||||
" \"\"\"Compute b ** e mod m, efficiently.\"\"\"\n",
|
||||
" return (1 if e == 0 else\n",
|
||||
" b if e == 1 else\n",
|
||||
" modpow((b * b) % m, e // 2, m) if e % 2 == 0 else\n",
|
||||
" modpow((b * b) % m, e // 2, m) * b % m)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6baf40c3-db36-40ba-b2d7-45afa21f811f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can see this function is fast and that it correctly computes the remainder of 1 for the following:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"id": "01838139-52dd-48a0-8f62-bcdef351345e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 25 μs, sys: 1 μs, total: 26 μs\n",
|
||||
"Wall time: 26.9 μs\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"1"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"a = 123456789012345678901234\n",
|
||||
"n = 357686312646216567629137\n",
|
||||
"\n",
|
||||
"%time modpow(a, n - 1, n)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4be033d1-93a0-4ef5-853e-81da0782ed20",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"I hope you have enjoyed this excursion into the world of primes!"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python [conda env:base] *",
|
||||
"language": "python",
|
||||
"name": "conda-base-py"
|
||||
},
|
||||
"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.13.9"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
Reference in new issue
Block a user