Add files via upload
This commit is contained in:
Peter Norvig authored and 
GitHub committed
1 parent
e9e1c4c43b
commit
4578110de5
1 file changed
+132
-137
+132
-137
@@ -23,14 +23,14 @@
|
||||
" 37 is prime\n",
|
||||
" 7 is prime\n",
|
||||
"\n",
|
||||
"Numbers like this are called [**truncatable primes**](https://en.wikipedia.org/wiki/Truncatable_prime). \n",
|
||||
"Numbers like this are called [**left-truncatable primes**](https://en.wikipedia.org/wiki/Truncatable_prime). \n",
|
||||
"\n",
|
||||
"I thought I would write a program to find all the truncatable primes. (It is not immediately obvious that there is a finite number of them, but it turns out to be true.) My function `truncatable_primes` below works as follows:\n",
|
||||
"My function `left_truncatable_primes` works as follows:\n",
|
||||
"- Start with the one-digit primes, (2, 3, 5, 7).\n",
|
||||
"- On each iteration creates a list of new numbers that are one digit longer by placing each possible digit 1–9 to the left of each prime.\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 longer primes.\n",
|
||||
"- In the end, the function returns all the truncatable primes."
|
||||
"- 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."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -42,12 +42,12 @@
|
||||
"source": [
|
||||
"from sympy import isprime # isprime(n) returns True if n is a prime\n",
|
||||
"\n",
|
||||
"def truncatable_primes(primes=(2, 3, 5, 7)) -> tuple[int]:\n",
|
||||
" \"\"\"Truncatable primes, in ascending order.\"\"\"\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 + truncatable_primes(new_primes)\n",
|
||||
" return primes + left_truncatable_primes(new_primes)\n",
|
||||
" else:\n",
|
||||
" return primes"
|
||||
]
|
||||
@@ -78,9 +78,9 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"TP = truncatable_primes()\n",
|
||||
"LTP = left_truncatable_primes()\n",
|
||||
"\n",
|
||||
"len(TP)"
|
||||
"len(LTP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -88,7 +88,7 @@
|
||||
"id": "a32b56c4-587b-48d5-91ec-1abcef2e909f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are **4260** of these truncatable primes. Here are the smallest and largest 16 of them:"
|
||||
"There are **4260**; here are the smallest and largest 16 of them:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -109,7 +109,7 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"TP[:16]"
|
||||
"LTP[:16]"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -145,7 +145,7 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"TP[-16:]"
|
||||
"LTP[-16:]"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -155,7 +155,7 @@
|
||||
"source": [
|
||||
"## Right-Truncatable Primes\n",
|
||||
"\n",
|
||||
"What if you sharpen the pencil from the other end? That gives you **right-trunctable primes**:"
|
||||
"What if you sharpen the pencil from the other end? You get `right-trunctable primes`:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -194,9 +194,9 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"TPr = right_truncatable_primes()\n",
|
||||
"RTP = right_truncatable_primes()\n",
|
||||
"\n",
|
||||
"len(TPr)"
|
||||
"len(RTP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -222,7 +222,7 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(TPr)"
|
||||
"print(RTP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -251,7 +251,7 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"set(TP) & set(TPr)"
|
||||
"set(LTP) & set(RTP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -286,7 +286,7 @@
|
||||
" (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 TPr if is_omni_truncatable_prime(p)]"
|
||||
"[p for p in RTP if is_omni_truncatable_prime(p)]"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -318,7 +318,7 @@
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def isprime_simple(n: int) -> bool:\n",
|
||||
" \"\"\"Simple primality tester. A prime number is defined as an integer greater than 1 \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)"
|
||||
@@ -345,12 +345,12 @@
|
||||
"\n",
|
||||
"def isprime_faster(n: int) -> bool:\n",
|
||||
" \"\"\"More sophisticated primality tester: test 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",
|
||||
" 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(11, int(sqrt(n)) + 1, 2) \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)"
|
||||
]
|
||||
},
|
||||
@@ -359,9 +359,9 @@
|
||||
"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",
|
||||
"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",
|
||||
"## 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",
|
||||
@@ -369,21 +369,23 @@
|
||||
"\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. \n",
|
||||
"For example:\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 **a** value: 221 is definitely 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, but a lot of evidence that it is prime*|\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",
|
||||
"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*. Here is a simple implementation:"
|
||||
"Here is a simple implementation:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -433,12 +435,12 @@
|
||||
"id": "15ef6884-253d-4a3a-b58f-397ac73f552d",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"How common are false primes? Here's a function to report a table of false prime probabilities, over various values of *k*, and various ranges of numbers *n*:"
|
||||
"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": 22,
|
||||
"execution_count": 14,
|
||||
"id": "45a11a93-9ef3-48d1-80fa-e0113041346e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
@@ -447,13 +449,13 @@
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
" label k = 1 k = 5 k = 10 k = 20 \n",
|
||||
" 3 digits 1.469% 0.010% 0.001% 0.000%\n",
|
||||
" 4 digits 0.453% 0.019% 0.003% 0.001%\n",
|
||||
" 5 digits 0.125% 0.006% 0.001% 0.001%\n",
|
||||
" 6 digits 0.033% 0.001% 0.001% 0.000%\n",
|
||||
" 7 digits 0.006% 0.000% 0.001% 0.000%\n",
|
||||
" 8 digits 0.001% 0.000% 0.000% 0.000%\n",
|
||||
" 9 digits 0.000% 0.000% 0.000% 0.000%\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"
|
||||
@@ -469,8 +471,9 @@
|
||||
" 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=200_000):\n",
|
||||
" \"\"\"Report false prime probability for different digit lengths, and k values.\"\"\"\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",
|
||||
@@ -485,16 +488,16 @@
|
||||
"id": "f3fe25d6-5f0a-478b-90d0-6c89e092e92e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We see that with *k* = 1, 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 numbers, and *k* at 10 or more, false primes are very rare indeed. \n",
|
||||
"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 percentage. [Some Carmichael numbers](https://oeis.org/A033502) are worse than others. \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": 23,
|
||||
"execution_count": 15,
|
||||
"id": "7871aaa7-af1c-462b-85bb-d46818ccabea",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
@@ -503,24 +506,24 @@
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
" label k = 1 k = 10 k = 30 k = 100 k = 1000\n",
|
||||
" 1105 69.790% 2.780% 0.000% 0.000% 0.000%\n",
|
||||
" 1729 75.240% 5.300% 0.020% 0.000% 0.000%\n",
|
||||
" 2465 72.190% 3.970% 0.010% 0.000% 0.000%\n",
|
||||
" 2821 76.750% 6.990% 0.020% 0.000% 0.000%\n",
|
||||
" 6601 79.830% 11.100% 0.130% 0.000% 0.000%\n",
|
||||
" 8911 80.370% 10.280% 0.130% 0.000% 0.000%\n",
|
||||
" 41041 69.760% 2.770% 0.010% 0.000% 0.000%\n",
|
||||
" 101101 70.840% 3.330% 0.000% 0.000% 0.000%\n",
|
||||
" 294409 95.410% 60.260% 21.680% 0.580% 0.000%\n",
|
||||
" 340561 81.760% 13.160% 0.250% 0.000% 0.000%\n",
|
||||
" 825265 60.220% 0.750% 0.000% 0.000% 0.000%\n",
|
||||
" 56052361 99.180% 91.410% 77.270% 41.040% 0.010%\n",
|
||||
" 118901521 99.330% 93.910% 80.950% 50.660% 0.180%\n",
|
||||
" 172947529 99.350% 94.020% 83.620% 54.680% 0.260%\n",
|
||||
" 216821881 99.570% 94.530% 84.460% 58.240% 0.420%\n",
|
||||
"1299963601 99.770% 97.190% 91.170% 73.480% 4.660%\n",
|
||||
"2301745249 99.720% 97.530% 92.990% 78.900% 8.150%\n",
|
||||
"9624742921 99.840% 98.420% 95.070% 85.410% 20.370%\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"
|
||||
]
|
||||
}
|
||||
],
|
||||
@@ -538,7 +541,44 @@
|
||||
"id": "4ecd2289-66d5-4b0f-ac00-5db0ad92dde9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This is bad; some Carmichael numbers have false prime probability above 99% with *k* = 1. Increasing `k` helps, but even with *k* = 1000, some Carmichael numbers still have a high false prime probability."
|
||||
"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}"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -546,36 +586,11 @@
|
||||
"id": "bfd71dc7-7f6d-4c98-9886-5c242d9b297a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"So our Fermat test is mostly reliable, but a few composite numbers will consistently be identified as false primes.\n",
|
||||
"## Speed of Primality testing\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 start to 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",
|
||||
"**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",
|
||||
"# 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. Let's create a list of big primes with different digit-lengths:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"id": "2a10367d-45c1-4425-a630-96892555c1b7",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"big_primes = [ # One prime of each digit length from 5 to 32\n",
|
||||
" 29137, 629137, 7629137, 67629137, 567629137, 6567629137, 16567629137, 216567629137, 6216567629137, 46216567629137, \n",
|
||||
" 646216567629137, 2646216567629137, 12646216567629137, 312646216567629137, 6312646216567629137, 86312646216567629137, \n",
|
||||
" 686312646216567629137, 7686312646216567629137, 57686312646216567629137, 357686312646216567629137, 1000000000000000035000061, \n",
|
||||
" 59999999999899999999999999, 100000109999990000011000001, 2728487949505050529272727777, 24444666666888888889999999991, \n",
|
||||
" 100003100019100043100057100069, 9999999999999999777777775555331, 55555555555555555555555555555559]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6b2cbdc7-e8a9-4b83-aa0d-1e7cc83de071",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Next make a chart of the speeds of all the primality checking functions running on the big primes:"
|
||||
"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."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -587,10 +602,11 @@
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import timeit\n",
|
||||
"import sympy\n",
|
||||
"\n",
|
||||
"prime_tests = (isprime_simple, isprime_faster, isprime_fermat, isprime)\n",
|
||||
"big_primes = [sympy.nextprime(3 * 10 ** d) for d in range(4, 32)]\n",
|
||||
"\n",
|
||||
"def plot_run_times(primes, functions=prime_tests):\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",
|
||||
@@ -602,19 +618,18 @@
|
||||
" plt.ylabel('Run time in seconds (log scale)')\n",
|
||||
" plt.legend()\n",
|
||||
"\n",
|
||||
"def time_test(function, primes) -> tuple[list[int], list[float]]:\n",
|
||||
" \"\"\"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",
|
||||
"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",
|
||||
" repeat = (1 if function == isprime_simple else 200)\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 # Stop tsting a function once it takes over a second of run time\n",
|
||||
" elif time > 1e-3:\n",
|
||||
" repeat = 1\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"
|
||||
]
|
||||
},
|
||||
@@ -626,7 +641,7 @@
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAINCAYAAACTVPhqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAA2XFJREFUeJzs3XlcVPX+x/HXMICACu64C2plrrmWS2ppLpnappVlWVlatpj3drWy9dbV+pVZN9ptubaYS5mVS+aaaIlrLqnlghuCCLLDMDPn98fIKAIyDDMM4Pv5eMyDL+ccvufzHbDOZ76byTAMAxERERERkUL4+ToAEREREREpv5QwiIiIiIhIkZQwiIiIiIhIkZQwiIiIiIhIkZQwiIiIiIhIkZQwiIiIiIhIkZQwiIiIiIhIkfx9HUB5Z7fbOX78ONWrV8dkMvk6HBERERGRUjMMg7S0NBo2bIif34X7EJQwFOP48eM0adLE12GIiIiIiHjckSNHaNy48QWvUcJQjOrVqwOONzM0NNTH0YiIiIiIlF5qaipNmjRxPuteiBKGYuQNQwoNDVXCICIiIiKViitD7jXpWUREREREiqSEQUREREREiqSEQUREREREiqQ5DEWIiooiKioKm83m61BERESkErHZbOTm5vo6DKnkAgICMJvNHqnLZBiG4ZGaKqnU1FTCwsJISUnRpGcRERFxm2EYnDhxgtOnT/s6FLlI1KhRg/r16xc6sbkkz7jqYRAREREpA3nJQr169QgJCdGGsOI1hmGQmZlJQkICAA0aNChVfUoYRERERLzMZrM5k4XatWv7Ohy5CAQHBwOQkJBAvXr1SjU8SZOeRURERLwsb85CSEiIjyORi0ne31tp58woYRAREREpIxqGJGXJU39vShhERERERKRIShhEREREpEht27Zl5syZleY+pVGWMZan90OTnkVEREQqEJvdYOPBJBLSsqlXPYhukbUw+2mok3iPEgYRERGRCmLpzjhe/GE3cSnZzmMNwoJ4fmhrBrUt3dKZRdm5c6dX6vXVfUqjIsToDRqSJBWOzW6wYf8pvt92jA37T2Gza+9BERGp/JbujOOhL7bkSxYATqRk89AXW1i6M84r9z1/aMznn3/O5ZdfTkhICO3ateN///uf81xERARPPfUUAwcOpGbNmjRs2LDAsJqIiAiefvppBgwYQEhICGPHji30PhERETzxxBMMGDCAWrVqUb9+fd5880127txJ3759CQ4Opnnz5ixcuDBf/RaLhaeeeorGjRsTEhJCp06d+PHHH11u74Xa56kYXXmfzlfadpWGehiKEBUVRVRUFDabzdehyDl88cmKiIiINxiGQVaua88ZNrvB84t2UdhHZAZgAl5YtJueLeu4NDwpOMDs1go6e/bs4f777+e7777juuuuIzY2lunTpzN8+HDCwsIAeOONN/jkk0+YN28eP//8M6NHj6Z58+YMGzbMWc9bb73FF198waJFiwgKCiryfp9++ilz5szh6quv5vPPP2fChAm0aNGC999/n+7du/Puu+9y9913c/jwYWrUqAHA2LFjOXz4MD///DMREREsXbqUkSNHsnr1arp161bq9nkiRlffp3OVpl2lZTIMQx/PXkBJts0W78r7ZOX8P9i8/9y9d1cnJQ0iIlIuZWdnc/DgQSIjI50PyJkWK62fW+aTeHa/NJCQQNc+N27bti1jx45l4sSJ/Pzzz9x6662cPHmSKlWqFLg2IiKCHj168NVXXzmPjR8/nj179rB69WrnNUOGDCEqKqrI++Rdd9ttt/Hqq68CjgSratWqPPnkk7z44ouA41P3oKAgVq1aRZ8+fThw4AAtW7bkwIEDREREOOu+++67qVq1Ku+9994F21pc+zwRo6vv07n3crddhf3d5SnJM66GJEmFYLMbvPjD7iI/WQF48YfdGp4kIiLiRb169SIyMpL27dszdepUVq5cidVqzXdNly5d8n3ftWtXdu/ene9YmzZtXLpfy5YtnWWTyURoaGi+Y4GBgYSEhJCcnAzA5s2bMQyDli1b4u/vj9lsxs/Pj9mzZ3PgwAGPtK+0MeZx5X3KU9p2lZaGJEmFsPFgUoExm+cygLiUbDYeTKJ7i9plF5hIRWa3Qex6SI+HauHQrAf4mX0dlchFIzjAzO6XBrp07caDSYz5NKbY6z67tyvdImu5dG93hISEsGnTJpYvX86KFSt44IEHqFKlCmvWrKFu3bqAa5uFBQYGunS/wuq6UP12ux0/Pz/S0tIIDg526R7ncqV9pY2xJNfkKW27Sks9DFIhJKQVnSy4c53IRW/3IpjZFj6/ARbc7/g6s63juIiUCZPJREigv0uvqy+pS4OwIIp6xDThmNN39SV1XaqvNDsABwQEcP311/PGG2/w559/curUKRYsWOA8HxOTP7GJiYnh8ssvd/t+JdGxY0fsdjs///yz23UU1z5PKcn75Il2lYYSBqkQ6lUvekKUO9eJXNR2L4K5d0Pq8fzHU+Mcx5U0iJQ7Zj8Tzw9tDVAgacj7/vmhrb2+H8O8efN45pln2Lt3Lzk5OURHR5OSkkKLFi2c18yfP58vvviC1NRUFixYwOeff84TTzzh1bjyXHrppdx555088sgjLFu2jPT0dP766y9eeeUVPvnkk2J/3pX2eUpJ3qfStqu0lDBIhdAtshYNwopOBvI+WXGlG1bkoma3wdLJcKEZQUunOK4TkXJlUNsGvHdXJ+qf9//D+mFBZbbwxw033EC1atUYPnw4NWrU4KGHHmLGjBlcd911zmsmTZrE7Nmzadq0KY888givvPIKN954o9djy/Ppp59y//3389BDD1G7dm1uuOEGLBYLI0aMKPZnXWmfp5T0fSpNu0pLqyQVQ6sklR9Ld8Yx/ostBY5rlSSREjj4q2P4UXHu+REir/Z+PCIXiQutVlNS5Xmn54iICKZOnercW0EKV1bvk6dWSdKkZ6kwBrapT3hoFeJTc/Idr699GERclx7v2etEpMyZ/Uxa4EPKlBIGqTC2HD5NfGoOQf5+vHtXZ9Kyc8vdJysi5V61cM9eJyJSgSxdupTBgwcXeV4DbwqnhEEqjPmbjwJwffsGXNuqno+jEamgTscWc4EJQhs6llgVESmhQ4cO+TqECxo0aFC5SArK+/t0Pk16lgohO9fGj9sdK7qM6NzEx9GIVFBb/gffP3LOgSLWWhk0XfsxiIiIkxIGqRCW7TpBWo6VxjWDuVIrIYmUXMwsWPQoYEDXB2Dk/yD0vHk/oQ0dx1sP80mIIiJSPl00Q5Jyc3Px9y/dRiXiO3nDkW7p1Bg/zVcQKZnfP4QlTzrKVz0MA/8DJhO0ukE7PYuISLEqfQ+DYRhMmDCBWrVqUatWLT777DNfhyQlFJeSxbq/EwFHwiAiJbAh6myy0OOxs8kCOJKDyKuh3a2Or0oWRESkEJW+h2Hx4sWsX7+e48ePc+LECXr16sWQIUOoW7eur0MTF3275RiGAVdG1qJp7RBfhyNScUS/Bcufc5Sv/gdc++zZZEFERMRFFaKHYevWrcyZM4eEhIRCz1ssFlatWsWiRYuIi4vLd2716tXcddddVK9enUsuuYRevXqxYcOGsghbPMAwDOdwpFs7q3dBxGVrXz+bLPSZomRBRETcVq4Thl9++YWrrrqK2267jTvuuIPdu3cXuObAgQNcfvnljBs3jjfffJOWLVvy0UcfOc8nJSVRq9bZSbK1a9cmKSmpTOKX0ttyOJmDiRmEBJq5vp02ZhNxyepXYeW/HeVrpsI1TylZEBG3jR8/noULF1aa+xTHMAzee+89RowYwaBBg8jIyPB1SD5XrhOG9PR03nzzTVauXFnkNQ8++CCRkZHs3r2bVatW8dZbbzFhwgQOHjwIQHh4OMePH3def+zYMcLDtSFRReHce6FdA6pWqfQj6ERKxzBg5cuw+j+O7/u/AH2e9GlIIuIFdhsc/BV2zHd8tdu8ert169aVyb4BZXWf4nz99dc8//zz3HLLLUycOJEqVaqUus4HHniAH374wQPR+Ua5ThhuvPFGunfvXuT5uLg4VqxYwWOPPYa/v+NhcsyYMVSvXp25c+cCMGTIEGbNmsW2bdv46aef2LRpE1dffXWRdebk5JCamprvBY5VlgCsVitWq9V5rLCyxWLBZrMVWrbb7c77FFU2DAPDMAqUAex2e5Fli8UCgM1mK7RstVrztaO8tykjO5cftjuGmN3csWGlaFNl/D2pTeWkTdnZ8MsLsPb/AGDAK9h7PF6x21QZf09q00Xfprx75NVT4vLu72FmW/j8BlhwP3x+A8bMto7j51xvGEaBct79S1p+7733uOmmmy54TanadKb83nvvceONNxZ5jSfbdKHyb7/9Ru/evbntttsYNGgQfn5+pW7fr7/+yqFDh3zSprx7FvXvyRXlOmEozs6dOwFo27at85i/vz+XX345O3bsAKBnz548/vjj3HHHHbzwwgt89dVXVKtWrcg6p02bRlhYmPPVpIljk7Dly5cDsGLFClasWAHAkiVLWLduHQALFy4kJiYGgLlz57J9+3YAZs+ezZ49ewD4+OOPOXDgAABRUVEcO3YMgBkzZpCY6FgFaPr06aSlpWGxWJg+fToWi4W0tDSmT58OQGJiIjNmzAAcvSVRUVGAY2jWxx9/DMCePXuYPXs2ANu3b3cmTzExMc6uvnXr1rFkyZJy3aYvVu8gPcdKk1rBVM+OrxRtqoy/J7WpHLQpNZUt06+H6JkALDf3hx6PVOw2Vcbfk9p0Ubfp119/BRyjJ/I+jExNTSU9PR2A06dPO4e+JCcnk5mZCcCpU6fIzs4GIO33L2HuPZB6duSEo6I4x/Hdizhx4oTzAfTEiRPOh8cTJ04AjofG+Ph4wPHQmDc/NCcnx/m+ZGdnc+rUKQAyMzOZNWsWW7duJSMjg9OnT7N3714eeeQRhg8fzpNPPskff/zhbNOdd97JvHnzePPNNxkxYgR3330Line truncated
|
||||
"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>"
|
||||
]
|
||||
@@ -648,31 +663,11 @@
|
||||
"- `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_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 ten times faster than `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": 19,
|
||||
"id": "4f377075-e352-46f9-a716-f1de673a7086",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"for fn in prime_tests:\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"
|
||||
"- `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.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -684,7 +679,7 @@
|
||||
"\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",
|
||||
"## 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",
|
||||
@@ -699,7 +694,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"execution_count": 19,
|
||||
"id": "e71eb140-7e1a-4146-87a7-8c4eebf61486",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
@@ -722,7 +717,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"execution_count": 20,
|
||||
"id": "01838139-52dd-48a0-8f62-bcdef351345e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
@@ -730,8 +725,8 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 21 μs, sys: 1e+03 ns, total: 22 μs\n",
|
||||
"Wall time: 22.2 μs\n"
|
||||
"CPU times: user 25 μs, sys: 1 μs, total: 26 μs\n",
|
||||
"Wall time: 26.9 μs\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -740,7 +735,7 @@
|
||||
"1"
|
||||
]
|
||||
},
|
||||
"execution_count": 21,
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -763,9 +758,9 @@
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3 (ipykernel)",
|
||||
"display_name": "Python [conda env:base] *",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
"name": "conda-base-py"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
|
||||
Reference in new issue
Block a user