Add files via upload
This commit is contained in:
Peter Norvig authored and 
GitHub committed
1 parent
051b5ed609
commit
8a40e42170
1 file changed
+195
-180
+195
-180
@@ -23,9 +23,9 @@
|
||||
" 37 is prime\n",
|
||||
" 7 is prime\n",
|
||||
"\n",
|
||||
"Numbers like this are called [**truncatable primes**](https://en.wikipedia.org/wiki/Truncatable_prime). I thought I would write a program to find other truncatable primes.\n",
|
||||
"Numbers like this are called [**truncatable primes**](https://en.wikipedia.org/wiki/Truncatable_prime). \n",
|
||||
"\n",
|
||||
"My function `left_truncatable_primes` below starts with the list of one-digit primes: [2, 3, 5, 7]. Then it places each possible digit 1–9 to the left of each of the one-digit primes, giving a list of two-digit candidates. From the candidates it filters out just the primes. If there are any new primes, it recursively adds digits to them, giving us three-digit primes, then four-digit primes, and so on, stopping when there are no more new primes. In the end, the function gathers up all the truncatable primes and returns them in sorted order."
|
||||
"I thought I would write a program for this. My function `left_truncatable_primes` below starts with the list of one-digit primes: [2, 3, 5, 7]. Then it creates a list of one-digit-longer new primes by placing each possible digit 1–9 to the left of each of the one-digit primes, and checking if each resulting number is prime. If there are any new primes, it recursively adds digits to them, giving us three-digit primes, then four-digit primes, and so on, stopping when there are no more new primes. In the end, the function gathers up all the truncatable primes and returns them in sorted order."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -35,13 +35,12 @@
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from sympy import isprime # isprime checks if a number is a prime\n",
|
||||
"from sympy import isprime # isprime(n) returns True if n is a prime\n",
|
||||
"\n",
|
||||
"def left_truncatable_primes(starting_primes=[2, 3, 5, 7]) -> list[int]:\n",
|
||||
"def left_truncatable_primes(primes=[2, 3, 5, 7]) -> list[int]:\n",
|
||||
" \"\"\"All left-truncatable primes, in ascending order.\"\"\"\n",
|
||||
" candidates = [int(d + str(p)) for d in \"123456789\" for p in starting_primes]\n",
|
||||
" new_primes = [n for n in candidates if isprime(n)]\n",
|
||||
" return starting_primes + (left_truncatable_primes(new_primes) if new_primes else [])"
|
||||
" new_primes = [dp for d in \"123456789\" for p in primes if isprime(dp := int(d + str(p)))]\n",
|
||||
" return (primes + left_truncatable_primes(new_primes)) if new_primes else primes"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -229,11 +228,10 @@
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def right_truncatable_primes(starting_primes=[2, 3, 5, 7]) -> list[int]:\n",
|
||||
"def right_truncatable_primes(primes=[2, 3, 5, 7]) -> list[int]:\n",
|
||||
" \"\"\"All right-truncatable primes, in ascending order.\"\"\"\n",
|
||||
" candidates = [10 * p + d for p in starting_primes for d in (1, 3, 7, 9)]\n",
|
||||
" new_primes = list(filter(isprime, candidates))\n",
|
||||
" return starting_primes + (right_truncatable_primes(new_primes) if new_primes else [])"
|
||||
" new_primes = [pd for p in primes for d in (1, 3, 7, 9) if isprime(pd := 10 * p + d)]\n",
|
||||
" return (primes + right_truncatable_primes(new_primes)) if new_primes else primes"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -241,7 +239,7 @@
|
||||
"id": "976c7f1f-79c5-48c3-8909-62b809993057",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"(Note that I only try placing the digits (1, 3, 7, 9) to the right: placing an even digit or a 5 would always result in a composite number. Also, note that I don't have to do string concatenation to form the new candidates; it is simpler to do `10 * p + d`.)\n",
|
||||
"(Note that I only try placing the digits (1, 3, 7, 9) to the right: placing an even digit or a 5 would always result in a composite number. Since there are fewer digits to add, I expect there to be fewer right-truncatable primes in total. Also, note that I don't have to do string concatenation to form the new candidates; it is simpler to do `10 * p + d`.)\n",
|
||||
"\n",
|
||||
"Let's find and count the right-truncatable primes:"
|
||||
]
|
||||
@@ -354,63 +352,16 @@
|
||||
"\n",
|
||||
"# Note on Primality Checking\n",
|
||||
"\n",
|
||||
"I was very impressed by the speed of `sympy.isprime`, so I wanted to look more deeply into the topic of primality checking. I'm going to write some prime-checking functions, so the first thing I will do is define a function to test that a prime checker correctly returns True when passed a prime and False when passed a composite number, for some given test cases."
|
||||
"I was very impressed by the speed of `sympy.isprime`, and I wanted to examine the topic of primality checking. I'll start by defining `isprime_simple`, which follows the definition of a prime number almost verbatim:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"id": "4f377075-e352-46f9-a716-f1de673a7086",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from typing import Callable, Iterable\n",
|
||||
"\n",
|
||||
"def test_prime_checker(checker: Callable, \n",
|
||||
" primes=(2, 3, 5, 7, 11, 37, 73, 101, 11939, 65537, 117223, 7629137),\n",
|
||||
" composites=(0, 1, 4, 6, 8, 9, 10, 256, 11939*11939, 11939*117223)) -> Callable:\n",
|
||||
" \"\"\"Test that a primality checking function correctly handles some primes and composites.\"\"\"\n",
|
||||
" for n in primes:\n",
|
||||
" assert checker(n) == True, f'{n} should be prime'\n",
|
||||
" for n in composites:\n",
|
||||
" assert checker(n) == False, f'{n} should be composite'\n",
|
||||
" return checker # This allows us to use @test_prime_checker as a @decorator"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c08ed813-2942-42da-aed3-d5044b0b403a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Here's how we would use this to check the `isprime` function:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "008cd6fd-86eb-4b45-bea9-c030273c4eb1",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"assert test_prime_checker(isprime)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6b307b24-a93e-4353-9bee-fba739072bd1",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"But I can also use `@test_prime_checker` as a decorator to test the function right where it is defined; I'll do that for `isprime_simple`, which follows the definition of a prime number almost verbatim:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"id": "1d460931-3d70-485d-807f-8fb360ef4ef2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"@test_prime_checker\n",
|
||||
"def isprime_simple(n: int) -> bool:\n",
|
||||
" \"\"\"Simple primality checker. 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",
|
||||
@@ -425,21 +376,20 @@
|
||||
"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 small and one big. The small idea is to cut the run time in half by only checking the odd divisors. Once we determine that 2 is a prime and all other even numbers are composite, we only need to check odd divisors. The big idea is that we only need to check 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."
|
||||
"We can speed things up a bit with two ideas, one small and one big. The small idea is to cut the run time in half by only checking the odd divisors (after determining that 2 is a prime and all other even numbers are composite). The big idea is that we only need to check 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": 13,
|
||||
"execution_count": 11,
|
||||
"id": "647fb66b-32ae-494a-8945-cf7358c2f958",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from math import sqrt\n",
|
||||
"\n",
|
||||
"@test_prime_checker\n",
|
||||
"def isprime_faster(n: int) -> bool:\n",
|
||||
" \"\"\"More sophisticated primality checker: go to square root, checking odd divisors only.\"\"\"\n",
|
||||
" \"\"\"More sophisticated primality checker: check only odd divisors up to √n.\"\"\"\n",
|
||||
" if n <= 10: # Handle small numbers up to 10\n",
|
||||
" return n in (2, 3, 5, 7)\n",
|
||||
" elif n % 2 == 0: # Even numbers other than 2 are composite\n",
|
||||
@@ -458,37 +408,37 @@
|
||||
"\n",
|
||||
"Fortunately, [Pierre de Fermat](https://en.wikipedia.org/wiki/Pierre_de_Fermat) provided that breakthrough in 1640, [showing](https://en.wikipedia.org/wiki/Fermat_primality_test) that if *n* is prime and *a* is not divisible by *n*, then *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*). \n",
|
||||
"\n",
|
||||
"That means we can check if *n* is prime by choosing a random *a* and testing 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 be composite. 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, and can sometimes be wrong.\n",
|
||||
"Some examples:\n",
|
||||
"That means we can check if *n* is prime by choosing a random *a* and testing 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 be composite. 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. \n",
|
||||
"For example:\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|*221 is definitely composite*|\n",
|
||||
"|221|2|16|*try again with new **a** value: 221 is definitely composite*|\n",
|
||||
"|5|2|1|*5 could be prime or composite*|\n",
|
||||
"|5|3|1|*5 could be prime or composite*|\n",
|
||||
"|5|4|1|*5 could be prime or composite, but we have a lot of evidence that it is prime*|\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Here is an implementation:"
|
||||
"This is called a [Monte Carlo algorithm](https://en.wikipedia.org/wiki/Monte_Carlo_algorithm); an algorithm that uses randomization, and can sometimes be wrong. Here is an implementation:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"execution_count": 12,
|
||||
"id": "476282dd-a48a-4eb5-8fed-8495d4878ec2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import random\n",
|
||||
"from typing import Iterable\n",
|
||||
"\n",
|
||||
"def sample(lo: int, hi: int, k: int) -> Iterable[int]:\n",
|
||||
" \"\"\"Randomly sample k integers from range(lo, hi), one by one.\"\"\"\n",
|
||||
" return (random.randrange(lo, hi) for _ in range(k))\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))\n",
|
||||
" \n",
|
||||
"@test_prime_checker\n",
|
||||
"def isprime_fermat(n: int, k=20) -> bool:\n",
|
||||
" \"\"\"n is probably a prime if this returns True; definitely composite if it returns False.\"\"\"\n",
|
||||
" if n <= 10:\n",
|
||||
@@ -502,7 +452,7 @@
|
||||
"id": "1ec7bba6-ea89-4255-a73c-bc25a32da883",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Note that the built-in [modular exponentiation](https://en.wikipedia.org/wiki/Modular_exponentiation) function `pow(b, e, m)` computes *b*<sup>*e*</sup> (mod *m*), and does it in an efficient way that doesn't have to multiply *b* by itself *e* times, and doesn't have to deal with any numbers larger than *m*."
|
||||
"(Note that `pow(a, n - 1, n)` efficiently computes *a*<sup>(*n* - 1)</sup> (mod *n*); more on this at the bottom of this notebook.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -510,23 +460,22 @@
|
||||
"id": "15ef6884-253d-4a3a-b58f-397ac73f552d",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**The problem** is that there are some composite numbers *n* for which *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*) for some values of *a*, and thus `isprime_fermat` can lie: it can incorrectly report a **false prime** for a composite number *n*. \n",
|
||||
"**The problem** is that there are some composite numbers *n* for which *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*) for some values of *a*. In other words `isprime_fermat` can lie: it can incorrectly report a **false prime** for a composite number *n*. \n",
|
||||
"\n",
|
||||
"How common are these false primes? Here's a function to estimate the percentage of false primes in the range `lo` to `hi` when `isprime_fermat` is allowed `k` choices of the `a` parameter:"
|
||||
"How common are these false primes? Here's a function to test a sample of numbers and report the the percentage of false primes when `isprime_fermat` is allowed `k` choices of the `a` parameter:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"execution_count": 13,
|
||||
"id": "45a11a93-9ef3-48d1-80fa-e0113041346e",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def false_prime_percent(lo: int, hi: int, k: int, repeat=100_000) -> float:\n",
|
||||
" \"\"\"The estimated probability of a false prime from isprime_fermat(n, k) for n in range(lo, hi).\"\"\"\n",
|
||||
" composites = (n for n in sample(lo, hi, repeat) if not isprime(n))\n",
|
||||
" counts = Counter(isprime_fermat(n, k) for n in composites) \n",
|
||||
" return 100 * counts[True] / (counts[True] + counts[False])"
|
||||
"def false_prime_percent(samples: Iterable[int], k: int) -> float:\n",
|
||||
" \"\"\"The estimated percentage 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 100 * sum(isprime_fermat(n, k) for n in composites) / len(composites)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -534,37 +483,37 @@
|
||||
"id": "778c3a9b-dcc7-4ad2-9092-d66aecd1d93c",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We'll choose ranges that represent digit-lengths, and for now use 1 for the value of `k`:"
|
||||
"We'll choose sample ranges that represent digit-lengths, sample 100,000 times for each digit-length, and for now use 1 for the value of `k`:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"execution_count": 14,
|
||||
"id": "9eba96df-c298-4d71-b68a-3bf36bf15095",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{3: 1.5219868933421978,\n",
|
||||
" 4: 0.44898977301072585,\n",
|
||||
" 5: 0.15980955110048164,\n",
|
||||
" 6: 0.02713174088102189,\n",
|
||||
" 7: 0.007483109552723852,\n",
|
||||
" 8: 0.003178740582981023,\n",
|
||||
"{3: 1.4432622005587588,\n",
|
||||
" 4: 0.46280016787849226,\n",
|
||||
" 5: 0.12571126108244013,\n",
|
||||
" 6: 0.0335378057620114,\n",
|
||||
" 7: 0.007482149728504852,\n",
|
||||
" 8: 0.0021168501270110076,\n",
|
||||
" 9: 0.0,\n",
|
||||
" 10: 0.0,\n",
|
||||
" 11: 0.0,\n",
|
||||
" 12: 0.0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 16,
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{d: false_prime_percent(10**(d-1), 10**d, k=1) for d in range(3, 13)}"
|
||||
"{d: false_prime_percent(sample(10**(d-1), 10**d, 100_000), k=1) for d in range(3, 13)}"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -572,12 +521,14 @@
|
||||
"id": "f3fe25d6-5f0a-478b-90d0-6c89e092e92e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We see that 3-digit composites are falsely called prime about 1.5% of the time, but onnce we get to 10-digit composites, false primes are very rare indeed. Let's see how we reduce false primes by allowing up to 30 choices of `a`:"
|
||||
"We see that 3-digit composites are falsely called prime about 1.4% of the time, and 4-digit comnposites about 0.5%, but once we get to 10-digit composites, false primes are very rare indeed. \n",
|
||||
"\n",
|
||||
"Let's see how we reduce false primes by allowing `k=25` choices of `a`:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"execution_count": 15,
|
||||
"id": "8880c588-32f5-46e8-b59c-f9b49d492d3b",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
@@ -587,7 +538,7 @@
|
||||
"{3: 0.0,\n",
|
||||
" 4: 0.0,\n",
|
||||
" 5: 0.0,\n",
|
||||
" 6: 0.001082461951462406,\n",
|
||||
" 6: 0.0,\n",
|
||||
" 7: 0.0,\n",
|
||||
" 8: 0.0,\n",
|
||||
" 9: 0.0,\n",
|
||||
@@ -596,13 +547,13 @@
|
||||
" 12: 0.0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 17,
|
||||
"execution_count": 15,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{d: false_prime_percent(10**(d-1), 10**d, k=30) for d in range(3, 13)}"
|
||||
"{d: false_prime_percent(sample(10**(d-1), 10**d, 100_000), k=25) for d in range(3, 13)}"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -617,27 +568,27 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"execution_count": 16,
|
||||
"id": "7871aaa7-af1c-462b-85bb-d46818ccabea",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{561: 55.4,\n",
|
||||
" 1105: 68.8,\n",
|
||||
" 1729: 77.0,\n",
|
||||
" 2465: 76.2,\n",
|
||||
" 2821: 76.6,\n",
|
||||
" 6601: 81.7,\n",
|
||||
" 8911: 80.0,\n",
|
||||
" 41041: 68.2,\n",
|
||||
" 101101: 70.3,\n",
|
||||
" 825265: 59.4,\n",
|
||||
" 321197185: 68.6}"
|
||||
"{561: 54.7,\n",
|
||||
" 1105: 69.8,\n",
|
||||
" 1729: 74.0,\n",
|
||||
" 2465: 72.0,\n",
|
||||
" 2821: 75.2,\n",
|
||||
" 6601: 79.7,\n",
|
||||
" 8911: 78.4,\n",
|
||||
" 41041: 69.8,\n",
|
||||
" 101101: 74.0,\n",
|
||||
" 825265: 61.4,\n",
|
||||
" 321197185: 63.5}"
|
||||
]
|
||||
},
|
||||
"execution_count": 18,
|
||||
"execution_count": 16,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -645,7 +596,7 @@
|
||||
"source": [
|
||||
"carmichael_numbers = (561, 1105, 1729, 2465, 2821, 6601, 8911, 41041, 101101, 825265, 321197185)\n",
|
||||
"\n",
|
||||
"{n: false_prime_percent(n, n + 1, k=1, repeat=1000) for n in carmichael_numbers}"
|
||||
"{n: false_prime_percent([n] * 1000, k=1) for n in carmichael_numbers}"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -653,12 +604,12 @@
|
||||
"id": "4ecd2289-66d5-4b0f-ac00-5db0ad92dde9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"That seems seriously bad; they all are above 50% error rate. But if we set `k=50`, we eliminate most (and on some runs all) of the errors:"
|
||||
"That seems seriously bad; they all are above 50% error rate. But if we set `k=40`, we eliminate most (and on some runs all) of the errors:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"execution_count": 17,
|
||||
"id": "de492047-77b6-46a3-8a77-46308f219efd",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
@@ -669,7 +620,7 @@
|
||||
" 1105: 0.0,\n",
|
||||
" 1729: 0.0,\n",
|
||||
" 2465: 0.0,\n",
|
||||
" 2821: 0.0,\n",
|
||||
" 2821: 0.01,\n",
|
||||
" 6601: 0.0,\n",
|
||||
" 8911: 0.0,\n",
|
||||
" 41041: 0.0,\n",
|
||||
@@ -678,13 +629,13 @@
|
||||
" 321197185: 0.0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{n: false_prime_percent(n, n + 1, k=50, repeat=10_000) for n in carmichael_numbers}"
|
||||
"{n: false_prime_percent([n] * 10_000, k=40) for n in carmichael_numbers}"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -698,25 +649,41 @@
|
||||
"\n",
|
||||
"# Speed of Primality Checking\n",
|
||||
"\n",
|
||||
"**The great thing** about `isprime_fermat` and `sympy.isprime` is that they are **very fast**, handling 32-digit primes in well under a millisecond. Let's make a chart:"
|
||||
"**The great thing** about `isprime_fermat` and `sympy.isprime` is that they are **very fast**, handling 32-digit primes in well under a millisecond. Let's create a list of big primes in ascending order of digit-length from 5 to 32, and measure how fast the functions are (while also verifying that they work correctly on Carmichael numbers):"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"execution_count": 18,
|
||||
"id": "be5f2fe7-b1e0-4a02-b7f0-7903e865c4b8",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"big_primes = ([int('357686312646216567629137'[-i:]) for i in range(5, 25)] # left-truncated primes\n",
|
||||
" + # The following primes sourced from t5k.org/curios/\n",
|
||||
" [1000000000000000035000061, 59999999999899999999999999, 100000109999990000011000001,\n",
|
||||
" 2728487949505050529272727777, 24444666666888888889999999991, 100003100019100043100057100069,\n",
|
||||
" 9999999999999999777777775555331, 55555555555555555555555555555559])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6b2cbdc7-e8a9-4b83-aa0d-1e7cc83de071",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Now let's make a chart of all the prime functions running on the big primes:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"id": "8109d9fc-bbae-4ef0-a75c-908d9746a557",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import timeit\n",
|
||||
"\n",
|
||||
"big_primes = ([int('357686312646216567629137'[-i:]) for i in range(5, 25)] +\n",
|
||||
" # The following primes sourced from t5k.org/curios/\n",
|
||||
" [1000000000000000035000061, 59999999999899999999999999, 100000109999990000011000001,\n",
|
||||
" 2728487949505050529272727777, 24444666666888888889999999991, 100003100019100043100057100069,\n",
|
||||
" 9999999999999999777777775555331, 55555555555555555555555555555559])\n",
|
||||
"\n",
|
||||
"def plot_run_times(functions=(isprime_simple, isprime_faster, isprime_fermat, isprime), primes=big_primes):\n",
|
||||
"def plot_run_times(primes, functions=(isprime_simple, isprime_faster, isprime_fermat, isprime)):\n",
|
||||
" \"\"\"For each primality-checking function, plot its run time on primes of different digit-lengths.\"\"\"\n",
|
||||
" plt.figure(figsize=(9, 6))\n",
|
||||
" for function in functions:\n",
|
||||
@@ -732,13 +699,13 @@
|
||||
" \"\"\"Time the function on each of the primes, stopping when one exceeds a second.\n",
|
||||
" Return a list of the digit sizes and a list of the corresponding run times.\"\"\"\n",
|
||||
" D, T = [], [] # D is length of primes in digits; T is time in seconds\n",
|
||||
" repeat = 100 # Start with 100 repeats, but reduce when times get longer\n",
|
||||
" repeat = (1 if function == isprime_simple else 500)\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:\n",
|
||||
" break\n",
|
||||
" break # Don't make a function do bigger primes if it took over a second on this one\n",
|
||||
" if time > 1/1000:\n",
|
||||
" repeat = 1\n",
|
||||
" return D, T"
|
||||
@@ -746,13 +713,13 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"execution_count": 20,
|
||||
"id": "b1ab0f3d-6ce6-4b8c-9d4d-bb78bb1ac432",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAINCAYAAACTVPhqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAv7lJREFUeJzs3Qd4VMXaB/D/pveQAgkJvQjSey8iSFFBRMUuiBf7FUQU9SrYLkURseDF8mHBBiKgSBEEpHcEQaT3QArpvezu97yz2SUJKZtkN9v+v/ucu7NnT87Ono1h3jMz72j0er0eREREREREpXArbScREREREREDBiIiIiIiKhd7GIiIiIiIqEwMGIiIiIiIqEwMGIiIiIiIqEwMGIiIiIiIqEwMGIiIiIiIqEwMGIiIiIiIqEweZb9EQqfT4dKlSwgMDIRGo+FFISIiIiKHJ2s3p6enIyoqCm5u5fchMGCogAQL9evXt+T3Q0RERERkFy5cuIB69eqVewwDhgpIz4LxYgYFBVnu2yEiIiIispG0tDR1U9zY1i0PA4YKGIchSbDAgIGIiIiInIk5Q+456ZmIiIiIiMrEgIGIiIiIiMrEgIGIiIiIiMrEOQxERERENkprWVBQAK1Wy+tPVuHp6Ql3d/dqn4cBQxnmzZunNv5HTERERJaWl5eHy5cvIysrixeXrDqhWVKmBgQEVO88eglvqdyUU8HBwUhNTWWWJCIiIrLIorAnTpxQd35r164NLy8vLg5LFidN/ISEBBWUNm/e/Jqehsq0cdnDQERERFTDvQsSNEgOfD8/P157shoJSM+ePYv8/PxqDU3ipGciIiIiG3BzYzOMbL/Ggjn4m0pERERERGViwEBERERERGViwEBERETkoLQ6PXacSsTPB2LUozy3phtuuAETJ0602vllvL0Mozlw4ADs0dkarF+jRo0wd+5c2ANOeiYiIiJyQGsOX8brK47gcmqOaV/dYB9MG94KQ9vUtcp7Ll26VOX2txaZCC7pZsPDw2GP6tt5/ayFPQxEREREDhgsPPHN/mLBgohNzVH75XVrCA0NRWBgoNWyR0kmn8jISHh42Oc9bXc7r5+1MGAgh1PT3a9EREQ1kTM/K6/ArC09Jx/Tfvkbpf3rZ9z32i9H1HHmnK8yS3IVHZL08ccfq/z+Pj4+iIiIwJ133lnsuKefflptkutf7si/+uqrxd5Lhty8+eabeOihh9Q6AI8++ug1Q37++OMP9fy3335Dx44d4evrixtvvBHx8fFYvXo1rr/+evWz9913X7FF8CRt7YwZM9C4cWP1M+3bt8eSJUvM+ozJycm4//77VUpS+dnmzZvjiy++UK9Zqn7mXJ+SUlJS8K9//UvVS84p73Pw4EHUBNcKj8jh2aL7lYiIyNqy87VoNfU3i5xLmpyxaTlo+9pas44/8sYQ+HlVrkm4d+9ePPPMM1i4cCF69eqFpKQkbNmypdgxX331FR555BHs3r1bHS8BQYMGDTB+/HjTMbNnz8bUqVMxbdq0ct/vtddew0cffaTWrRg9erTavL298d133yEjIwO33347PvzwQ0yZMkUdL8HCN998g/nz56sG/+bNm/HAAw+oxnb//v3LfS9puB85ckQ1+KUhf/LkSWRnZ1u0fuZen6LuuusuFZBIvSTI+OSTTzBw4EAcP35c9fxYEwOGMsybN09tWq3Wql8AVb77tWTsbex+/d8DnRg0EBER1YDz58/D398ft956qxqi1LBhQ3WHveR4//fee0/dgW/RogUOHTqknhdtEMtd8ueee870XO7gl+att95C7969VVka2S+99BJOnTqFJk2aqH3Su7Fx40bVIM/NzcX06dPx+++/o2fPnup1OW7r1q2qkV1RwCCfTT5Lly5dTD0hFalM/SpzfYyk7hJYSM+FBCLGYGv58uWq50SCDWtiwFCGp556Sm3GZbPJtmTYkfQslNX9KsuSyOs3tYqEu5tlFikhIiKqKb6e7upOvzl2n0nC2C/2VHjclw93RbfGoWa9d2XddNNNKkiQBvHQoUPVJnfRi65c3aNHj2ILh0nj/d1331U3Y42rDhsb5RVp166dqSzDn+R9jI1x4z5pUAvpEZDhP1LHknMkSgY1pXniiSdwxx13YP/+/Rg8eDBGjhypelEsVb/KXB8jGXokPRVhYWHF9kvPhwQm1saAgRyC/HEsObGrZNAgr8txPZsW/4+JiMqg0wLntgMZcUBABNCwF+BW+YYDEVWfNBzNHRbUt3ltNRxXethLu5EmTdDIYB91nLVuokmvgjSoZQz/2rVr1bAiGZazZ88e1KpVy+zzSC+FOYpmZpJrVTJTk+yTeQtCGtZi5cqViI6OLnac8e58eYYNG4Zz585h1apVWLdunRr2IzeR5Y6+JepXFfKZ6tatq653SZW53lXFgIEcQnx6jkWPI3J5R34B1kwB0i5dvRRBUcDQWUCrES5/eYjsmQQBMndPhuNKOFA0aDCGB/K6tXvcJVPQoEGD1CZzEKThumHDBowaNUq9vmvXrmLH79y5U80nKHn33NJatWqlAgMZWlTR8KOyyFyHMWPGqK1v3754/vnnyw0YqqIy16dTp06IjY1V19ycIVKWxoCBHEKdQB+LHkcEVw8WFj9UopkBIO2yYf/orxk0ENk5SfQhc/dKJgKJrKFEIL/++itOnz6Nfv36ISQkRN2NlzvoMhbfSBrskyZNwmOPPaZ6I2TSrwy5sTbp/Zg8eTKeffZZVac+ffogNTUV27ZtU9mFJAgoj/SWdO7cGa1bt1bzIX799VeV7cjSKnN9JCiTIUsyPOrtt9/Gddddh0uXLqleFBkKZu7QrqpiwEAOQcZgSvdrWcOSjN2v5ozVJIKrD0OSnoXyZgSteRFoeQuHJxHZOQkKZO6eDMeVHna5aSb/DtbEXD7pTZBF3GQYUk5Ojroz/v3336tGtpGkS5Ux9t26dVN3zSdMmGD1yblGkq5VegkkW5IENlJfuUv/8ssvV/izXl5eatKyTMCWrER9+/bFDz/8YPE6Vub6yJAmCcr+85//4OGHH0ZCQoJaD0ICNpkfYW0afWWS77og46RniUwlKiXbZkl6/Jv91+w3/llkliQiM5zZAnx1a8XHjfkVaNyXl5TICqSBfebMGbVGgKxh4IxknYEOHTpg7ty5tq6KS1+fnHJ+1yrTxuXCbeQwbmhRBz6e1/7KSs8CgwUiM8kEZ0seR0RETo9DkshhrD0Sh5x8HaKCfTD7rvZIyMit0e5XIqcg2ZAseRwRkYN5/PHH1aJupZHF3WSxNyqOAQM5jCX7LqrHOzvXQ69m4bauDpFjKsgtHMhX1mhUjSFbkqRYJSKqotLSf9qLN954Q02KLk1NDT//w46vT2kYMJBDkFzTW08kqPIdnevZujpEjunY6hLZkcpIyDh0Jic8E5HTqlOnjtrIfJzDQA5h2Z8x0OmBbo1C0TDMvEVeiKiIwz8Bix4AtHlAy1uBOxcAQSXSLkrPAlOqEhFRCexhILsnibyW7LtgGo5ERJW0fyHwy78NvQnt7gZu+xhw9wBajeRKz0REVCEGDGT3DlxIwamETPh6uuPmdtZdiIbI6ez6FFj9vKHceSxwy3uAW2Hnsps7U6cSEVGFXGJIkqzQJysPyqIin3/+ua2rQ1Wc7DysTSQCvBnjEplty5yrwUKPp4Bb514NFoiIiMzk9K2vgoICtez2xo0b1eIUstS3LKEdFhZm66qRGXLytfjl4CVV5nAkIjPJepwb3gK2zDY87/cCMOBlWSqUl5CIiCrN6W817d69Wy1THh0djYCAAAwbNgxr1661dbXITOuOxCE9pwDRtXzRowmDPCKzgoXfXr4aLAx6HbjxPwwWiJyVTmtYwf3QEsOjPLfyCsUTJ0602vnPnj0LjUaDAwcOwJa2bduGtm3bwtPTEyNHjoSrs/uAYfPmzRg+fDiioqLUL9Dy5cuvOWbevHlo1KiRWvK6e/fuKkgwunTpkgoWjKQcExNTY/UnywxHGtUpGm5cnI2ofNJQWPEMsPNjw/ObZwN9rPcPOxHZ2JFfgLltgK9uBX56xPAoz2W/lSxduhRvvvmm1c5fv359XL58GW3atIEtyeiUDh064MyZM/jyyy+rfb7XXntNnc9R2X3AkJmZifbt26ugoDSLFi1SX+q0adOwf/9+deyQIUMQHx9f43Uly4pLy8EW49oLnZgdiahc2gJg2WPA/q8BjZshE1K38bxoRM5KggJZVyXNMGzXJO2yYb+VgobQ0FAEBgZa5dx5eXlwd3dHZGQkPDxsO2r+1KlTuPHGG1GvXj3UqlUL9kKukS3YfcAgQ4jeeustNe+gNHPmzMH48ePx8MMPo1WrVmo5bz8/PyxYsEC9Lj0TRXsUpCz7ypKbm4u0tLRim8jPzzfNiZDNuK+0snyZWq221LJOpzO9T1llSSMqW8mykOPKKht/ieT9SitL/Yp+Dnv/TEv3X1RrL3RtFIL6IT5O8Zmc8XviZ7KD7yk/B1pZY+HQj9C7eSD/tk+Ajvfze+J/T/wbYed/9+Q88n6qrNVCn5sB5GVCl5NefjknDfrVL0Bf6orthr361VOgy0pRx8smPyvnKLVc+LmN16W8snFIkpTlZq4klJERHhEREbjzzjtNn0mOe+qpp9Qmc0jDw8PxyiuvmM4njzI6RFZdfuihh9QKy9Kekzv6MqJEbgLLMbIisjxfs2YNOnbsCF9fX9WQj4uLw8qVK3H99dern7333nuRlZVlOr9s06dPR+PGjdXPyA3lxYsXV/j5Tp8+rd4vMTER48aNU+UFCxao7+yRRx4xnU+S6cydO/fq96fTqfmy3bp1g7+/vwoyevfurT7PF198gddffx0HDx5U5zOeU342OTlZnbd27drqc8hnk+OM9TL2THz66afqveVam/M9FS0b61jWv09OETCUR/6j3LdvHwYNGmTa5+bmpp7v2LFDPZcv7vDhwypQyMjIwOrVq1UPRFlmzJiLine truncated
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAINCAYAAACTVPhqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAv4BJREFUeJzs3Qd4FNXaB/B/ek9ID4ReJfTeQYoCUhRQbAg2bHhFEUW9Kl5BihcRURD9vFiwImJDQKr03hHpoZMGpPfd/Z73bDYkIWU37Gbb/+cz7uzsZPbMbLKcd857znHR6XQ6EBERERERlcK1tI1EREREREQMGIiIiIiIqFxsYSAiIiIiojIxYCAiIiIiojIxYCAiIiIiojIxYCAiIiIiojIxYCAiIiIiojIxYCAiIiIiojK5l/0SCa1Wi0uXLiEgIAAuLi68KERERERk92Tu5rS0NNSoUQOuruW3ITBgKMO8efPUkpubi1OnTlnicyIiIiIisqrz58+jZs2a5e7jopPwgsqUkpKCatWqqYsZGBjIK0VEREREdi81NRW1atVCcnIygoKCyt2XLQwVMKQhSbDAgIGIiIiIHIkxKffs9ExERERERGViwEBERERERGViwEBERERERGViHwYiIiIiK5BxZ/Lz86HRaHj9ySI8PDzg5uZ208dhwEBERERUxWTY9suXLyMzM5PXnizaoVmGTPX397+p4zBgqGAeBkb9REREZO5JYWNjY9WdX5k0y9PTk5PDkkVasBITE3HhwgU0atToploaOA+DEWPUyti0Mh8Dh1UlIiKim5Wdna0Chjp16sDX15cXlCwmKysLZ86cQb169eDt7V3pOi47PRMRERFZgasrq2Fk/TkWjMHfVCIiIiIiKhMDBiIiIiIiKhMDBiIiIiI7pdHqsO3UFfy6/6J6lOeWdOutt+L555+32PEl317SaPbv3w9bdKYKy1e3bl3MmTMHtoCjJBERERHZoZWHL+M/vx/B5ZTswm3Vg7wxeUgMBjSvbpH3XLp0qRrb31Jq1aqlhpsNCwuDLapl4+WzFLYwEBEREdlhsPD013uLBQsiLiVbbZfXLSEkJAQBAQEWm5tChv6MioqCu7tt3tN2s/HyWQoDBrI7Vd38SkREVBVj5mfm5hu1pGXnYfJvf6O0f/0M29767Yjaz5jjyXtXJiVp/vz5anx/Ga4zMjISd999d7H9nn32WbXI0J1yR/6NN94o9l6ScjNlyhSMHj1aDev5xBNP3JDy89dff6nnf/75J9q0aQMfHx/06dMHCQkJWLFiBZo2bap+9oEHHig2CZ7MdTF9+nQ1nKj8TKtWrbBkyRKjzvHatWt48MEHER4ern62UaNG+Pzzz9Vr5iqfMdenpOTkZDz++OOqXHJMeZ8DBw6gKjhXeER2zxrNr0RERJaWladBzJt/muVYUuWMS81Gi7dWGbX/kbf7w9fTtCrh7t278dxzz2HRokXo2rUrrl69ik2bNhXb58svv8Rjjz2GnTt3qv0lIKhduzbGjh1buM+sWbPw5ptvYvLkyeW+31tvvYWPPvpIzVsxcuRItXh5eeHbb79Feno6hg0bhg8//BCTJk1S+0uw8PXXX2PBggWqwr9x40aMGjVKVbZ79epV7ntJxf3IkSOqwi8V+ZMnT6r5DMxZPmOvT1H33HOPCkikXBJkfPLJJ+jbty+OHz+uWn4siQED2V3za8nY29D8+vGotgwaiIiIqsC5c+fg5+eHwYMHqxQlmYRO7rCXzPd///331R34Jk2a4NChQ+p50Qqx3CV/8cUXC5/LHfzSTJ06Fd26dVPrUsl+9dVXcerUKdSvX19tk9aN9evXqwp5Tk4Opk2bhjVr1qBLly7qddlv8+bNqpJdUcAg5ybn0r59+8KWkIqYUj5Tro+BlF0CC2m5kEDEEGz98ssvquVEgg1LYsBAdkHSjqRloazmV5mWRF6/LSYKbq7mmaSEyOFpNcDZrUB6POAfCdTpCri6WbtURE7Jx8NN3ek3xs7Yq3j4810V7vfFIx3QsV6IUe9tqttuu00FCVIhHjBggFrkLnrRmas7d+5cbOIwqby/99570Gg0qi+AMFTKK9KyZcvCdUl/kvcxVMYN26RCLaRFQNJ/pIwl+0iUDGpK8/TTT2PEiBHYu3cvbr/9dtx1112qFcVc5TPl+hhI6pG0VISGhhbbLi0fEphYGgOGMsybN08t8qGR9cmXY8mOXSWDBnld9uvSoPgfExGV4shvwMpJQOql69sCawADZgIxQ3nJiKqYVByNTQvq0ShcpeNKC3tpN9KkChoV5K32s9RNNGlVkAq15PCvWrVKpRVJWs6uXbtQrVo1o48jrRTGKDoyk1yrkiM1yTbptyCkYi3++OMPREdHF9vPcHe+PAMHDsTZs2exfPlyrF69WqX9jBs3Tt3RN0f5KkPOqXr16up6l2TK9a4sdnoug/xiSP6a/OKT9SWkZZt1PyKnJsHC4tHFgwWRelm/XV4nIpslQYD03RMlwwHDc3nd0i3uMlJQv3798O677+LgwYMqnWjdunWFr+/YsaPY/tu3b1f9CUrePTe3mJgYFRhIalHDhg2LLZIGZAzp6zBmzBjVD2LOnDn49NNPzV5OU65P27ZtERcXp655yXOqiiFe2cJAdiEiwNus+xE5dRqStCyUl+C38hXglkFMTyKyYTLQh/TdKzkQSFQVDQSybNkynD59Gj179kRwcLC6Gy930CUX30Aq7BMmTMCTTz6pWiOk06+k3FiatH5MnDgRL7zwgipT9+7dkZKSgi1btqjRhSQQKI+0lrRr1w7NmjVT/SGWLVumRjsyN1OujwRmkrIk6VESoDVu3BiXLl1SrSiSCmZsaldlMWAguyA5mNL8WlZakqH51ZhcTSKnJn0WSrYsFKMDUi/q96vXowoLRkSmkqBA+u5JOq60sMtNM/l3sCr68kkajEziJmlI2dnZ6s74d999pyrZBjJcquTYd+zYUd01Hz9+vMU75xrIcK3SSiCjJUlgI+WVu/SvvfZahT/r6empOi1Li4mMStSjRw98//33Zi+jKddHUpokKPv3v/+NRx55BImJiWo+CAnYpH+EpbnoTBl81wmlpqaqoaskMpWolKw7StJTX++9Ybvha5GjJBEZ4dAS4KfHKt5vxP+AFtfHVCci85EKdmxsrJojQOYwcEQyz0Dr1q1VOg9Z7/qU97tmSh2XfRjIbvRsHA4v9xt/ZaVlgcECkZFkNCRz7kdERA6PKUlkN1YcikNOvha1Q3wwc0RLJKTlVGnzK5FDkKFTvasB2cll7OCiHy1J9iMickBPPfWU6sxcGpncTSZ7o+IYMJDd+HHPefV4T7ta6NLA8iMCEDmkbR+VHyyIATPY4ZmIbkppw3/airffflt1ii5NVaWf/2XD16c0DBjILpy/montp69C5jcZ3q6mtYtDZH+ku9q6KcCmghE4bhkCXNpTyjwMMzgPAxE5tIiICLWQ8RgwkF1YsueCeuzWIAzR1XysXRwi+yKTBa14Gdj1f/rn/d4Cur/AmZ6JiMgoDBjI5mm1Ovy0Vx8w3M3WBSLTaPKBX8cBB2VIQBdg0HtAh4JRklzdOHQqERFViAED2bztsVdw4VoWArzc0b9ZlLWLQ2Q/8rL1Q6geXQa4uAHDPgFa3mPtUhERkZ1hwEA2b8lufevC4FY14ONp2enkiRxGTjrw/QNA7AbAzQsY+SXQZKC1S0VERHbIKeZhkCm9ZapymYXws88+s3ZxyARp2XlYfviyWmc6EpGRsq4Bi+7SBwue/sCoJQwWiIio0hw+YMjPz8eECROwbt067Nu3D//9739x5coVaxeLjLT80GVk52lRP9wPbWtX43Ujqkh6AvDFYODCLv18C6N/Ber15HUjclRaDRC7ST+LuzzKcwvPUPz8889b7PhnzpyBi4sL9u/fD2vasmULWrRoAQ8PD9x1111wdg6fkrRz5040a9YM0dHR6vnAgQOxatUq3H///dYuGhnhx4J0JJl7Qb5AiKgcyeeAr+4Crp7Sz9T80M9AZDNeMiJHdeQ3YOWkUoZHnmmx4ZGXLl2qKtGWUqtWLVy+fBlhYdadb0luNrdu3RorVqyAv7//TR/vrbfewi+//GL1QMhhWxg2btyIIUOGoEaNGqrCKBe7pHnz5qFu3brw9vZGp06dVJBgcOnSpcJgQcj6xYsXq6z8VHmxSRnYffYaZBLn4W2vf4ZEVIqkE8DCgfpgIag28MgKBgtEjh4sLB5dPFgQqZf12+V1CwgJCUFAQIBFjp2bmws3NzdERUXB3d2697RPnTqFPn36oGbNmqhWzXYyHOQaWYPNBwwZGRlo1aqVCgpK88MPP6gocPLkydi7d6/at3///khISKjyspJ5LSmY2bln43BEBnrz8hKV5fJBYOEAIPUCENYYeHQlENqA14vI3iZXzM0wbslO1c+tAl1pB9I/SMuD7GfM8eS9K5GSNH/+fNU/VG7YRkZG4u677y6237PPPquWoKAg1WLwxhtvQFfkveRm75QpUzB69Gg1w/ITTzxxQ0qSzIgsz//880+0adMGPj4+qiIv9Ty5+9+0aVP1sw888AAyMzMLj63VajF9+nTUq1dP/YzUD5csWVLh+RneX9LXH330UbX+xRdfQKPR4LHHHis8nvSN/eCDD4r9rJS1Y8eO8PPzU0FGt27dcPbsWfXz//nPf3DgwAF1PMMxRXJyMh5//HGEh4er85Bzk/2KtkxIS4f0wZX3lmttDTafkiQpRLKUZfbs2Rg7diweeeQR9XzBggX4448/sHDhQrzyyiuqZaJoi4Ksy4dZlpycHLUYpKamqse8vLzCPhFCIl/ZJh96yXVDhCxLyXV53dXVVb2HNOmVtu7p6aneQ/Yvuu7l5aX+AOS9SluXssn+8kstS8l1eV3+UOW9Sq7b2jm5e3hi6d6LhZ2dHeGcHPFz4jnZwOcUvw/45h645KQCUS2hffAn5HkGwavgH0x+Tvx74neLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 900x600 with 1 Axes>"
|
||||
]
|
||||
@@ -762,58 +729,12 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"plot_run_times()"
|
||||
"plot_run_times(big_primes)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "408ba0fc-9bd9-4929-8ff3-1f9582ffa483",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can test that our two Fermat functions work on the big primes and on the carmichael numbers, and get some more timing data:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"id": "7eacb754-e6ee-482e-82fd-8725845aad7e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 1.25 ms, sys: 0 ns, total: 1.25 ms\n",
|
||||
"Wall time: 1.25 ms\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%time test_prime_checker(isprime, primes=big_primes, composites=carmichael_numbers);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 23,
|
||||
"id": "389cc455-6b9c-40f6-9fc2-2ae129159761",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 3.47 ms, sys: 0 ns, total: 3.47 ms\n",
|
||||
"Wall time: 3.47 ms\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%time test_prime_checker(isprime_fermat, primes=big_primes, composites=carmichael_numbers);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1347399a-24f2-4d5f-aa7f-822951126bf5",
|
||||
"id": "159833b0-afaf-486f-9c8e-1c8d2a11d9bf",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"What have we learned about timing?\n",
|
||||
@@ -822,10 +743,104 @@
|
||||
"- `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_faster` is more than 10,000 times faster than `isprime_simple` on 8-digit primes\n",
|
||||
"- `isprime_fermat` is more than 10,000 times faster than `isprime_faster` on 17 digit primes\n",
|
||||
"- `isprime` is about twice as fast as `isprime_fermat`; but that's quite good for `isprime_fermat` considering how simple it is\n",
|
||||
"- `isprime` produces no known incorrect answers; `isprime_fermat` will on rare occasion produce a false prime\n",
|
||||
"- `isprime` is about twice as fast as `isprime_fermat`\n",
|
||||
" - (but that's quite good for `isprime_fermat` considering how much simpler it is: 10 lines versus 800)\n",
|
||||
"- `isprime` produces no known incorrect answers; `isprime_fermat` will on rare occasion produce a false prime"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4c76de9d-4803-4db5-b6c5-811a6a1c6da1",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"It is good to get the timing results, but I should also run some test cases to demonstrate correctness:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"id": "4f377075-e352-46f9-a716-f1de673a7086",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"for fn in (isprime, isprime_simple, isprime_faster, isprime_fermat):\n",
|
||||
" assert all(fn(p) for p in (2, 3, 5, 7, 11, 37, 73, 101, 11939, 65537, 117223, 7629137)) # primes\n",
|
||||
" assert all(fn(n) is False for n in (0, 1, 4, 6, 8, 9, 10, 256, 11939*11939, 11939*65537)) # composites"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1347399a-24f2-4d5f-aa7f-822951126bf5",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This was fun for me! I remember learning about Fermat's Little Theorem and its application to primality testing in [Michael Rosen's](https://mathematics.brown.edu/people/michael-rosen) number theory class in college, and was impressed with it then, but I never worked through the details of how often the test gives a false prime result until now.\n",
|
||||
"\n",
|
||||
"I hope you have enjoyed this excursion into the land of prime numbers.\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)`, we'd have two problems: we'd need about 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 (1) apply the modulus to intermediate results, so we don't need petabytes of storage, and (2) cut the exponent in half each iteration, so we need only do *O*(log *n*) multiplications, not *O*(*n*). That's a big difference because log<sub>2</sub>(10<sup>24</sup>) is only 80. The basic idea is that *b*<sup>2*e*</sup> is equal to (*b* × *b*)<sup>*e*</sup>. Here is an implementation:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"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 follwoing:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 23,
|
||||
"id": "01838139-52dd-48a0-8f62-bcdef351345e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 24 μs, sys: 0 ns, total: 24 μs\n",
|
||||
"Wall time: 24.8 μs\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"1"
|
||||
]
|
||||
},
|
||||
"execution_count": 23,
|
||||
"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!"
|
||||
]
|
||||
}
|
||||
],
|
||||
|
||||
Reference in new issue
Block a user