Add files via upload
This commit is contained in:
Peter Norvig authored and 
GitHub committed
1 parent
c19ae6e8d5
commit
42960b72dd
1 file changed
+132
-573
+132
-573
@@ -25,7 +25,12 @@
|
||||
"\n",
|
||||
"Numbers like this are called [**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 starts with the list of one-digit primes, [2, 3, 5, 7], and 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. Then out of those numbers we find just the new primes. If there are any new primes, the function 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 returns all the truncatable primes in sorted order."
|
||||
"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",
|
||||
"- 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",
|
||||
"- 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."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -37,10 +42,10 @@
|
||||
"source": [
|
||||
"from sympy import isprime # isprime(n) returns True if n is a prime\n",
|
||||
"\n",
|
||||
"def truncatable_primes(primes=[2, 3, 5, 7]) -> list[int]:\n",
|
||||
"def truncatable_primes(primes=(2, 3, 5, 7)) -> tuple[int]:\n",
|
||||
" \"\"\"Truncatable primes, in ascending order.\"\"\"\n",
|
||||
" new_numbers = (int(d + str(p)) for d in \"123456789\" for p in primes)\n",
|
||||
" new_primes = [n for n in new_numbers if isprime(n)]\n",
|
||||
" new_primes = tuple(filter(isprime, new_numbers))\n",
|
||||
" if new_primes:\n",
|
||||
" return primes + truncatable_primes(new_primes)\n",
|
||||
" else:\n",
|
||||
@@ -83,7 +88,7 @@
|
||||
"id": "a32b56c4-587b-48d5-91ec-1abcef2e909f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are **4260** of these truncatable primes. Here are the smallest and largest of them:"
|
||||
"There are **4260** of these truncatable primes. Here are the smallest and largest 16 of them:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -95,7 +100,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[2, 3, 5, 7, 13, 17, 23, 37, 43, 47, 53, 67, 73, 83, 97, 113]"
|
||||
"(2, 3, 5, 7, 13, 17, 23, 37, 43, 47, 53, 67, 73, 83, 97, 113)"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
@@ -116,7 +121,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[66276812967623946997,\n",
|
||||
"(66276812967623946997,\n",
|
||||
" 67986315421273233617,\n",
|
||||
" 86312646216567629137,\n",
|
||||
" 315396334245663786197,\n",
|
||||
@@ -131,7 +136,7 @@
|
||||
" 57686312646216567629137,\n",
|
||||
" 95918918997653319693967,\n",
|
||||
" 96686312646216567629137,\n",
|
||||
" 357686312646216567629137]"
|
||||
" 357686312646216567629137)"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
@@ -150,7 +155,7 @@
|
||||
"source": [
|
||||
"## Right-Truncatable Primes\n",
|
||||
"\n",
|
||||
"What if you sharpen the pencil from the other end? The primes we found so far are called **left-truncatable primes**; there are also **right-trunctable primes**:"
|
||||
"What if you sharpen the pencil from the other end? That gives you **right-trunctable primes**:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -160,11 +165,11 @@
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def right_truncatable_primes(primes=[2, 3, 5, 7]) -> list[int]:\n",
|
||||
"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 = [n for n in new_numbers if isprime(n)]\n",
|
||||
" new_primes = tuple(filter(isprime, new_numbers))\n",
|
||||
" if new_primes:\n",
|
||||
" return primes + right_truncatable_primes(new_primes)\n",
|
||||
" else:\n",
|
||||
@@ -212,7 +217,7 @@
|
||||
"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"
|
||||
"(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"
|
||||
]
|
||||
}
|
||||
],
|
||||
@@ -281,350 +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 TP if is_omni_truncatable_prime(p)]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "7b316527-f437-439f-ac64-ef7980ecf6ff",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are only 9 of them, with 373 the largest.\n",
|
||||
"\n",
|
||||
"## Palindromic Primes\n",
|
||||
"\n",
|
||||
"What about primes that are both truncatable and [**palindromic**](https://en.wikipedia.org/wiki/Palindromic_prime) (the same when you reverse the digits)?"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"id": "5a705ef4-b2b7-4006-b518-6e030de8fd30",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{2, 3, 5, 7, 313, 353, 373, 383, 797, 76367, 79397, 7693967, 799636997}"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{p for p in TP + TPr if str(p) == str(p)[::-1]}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "cd0da3ad-54c0-47dc-979c-4eecc935425c",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are 13 **palindromic truncatable primes**, with only 4 of them longer than 3 digits.\n",
|
||||
"\n",
|
||||
"**Aside**: Now I'm curious about [**palindromic primes**](https://en.wikipedia.org/wiki/Palindromic_prime) in general. There are probably an infinite number of them, so I'll define `palindromic_primes(d)` to list all the palindromic primes up to `d` digits. To find them, I'll enumerate the \"first half\" of numbers, and then concatenate the reverse of the first half, forming both odd- and even-length numbers, e.g. when the first half is 123, `palindromes(123)` produces (12321, 123321). Then all we have to do is check each palindrome for primality."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "fc208a86-c3f4-4426-913d-62ec0a639576",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def palindromic_primes(digit_length: int) -> list[int]:\n",
|
||||
" \"\"\"Generate all palindromic primes of up to `digit_length` digits.\"\"\"\n",
|
||||
" end = 10 ** ((digit_length + 1) // 2)\n",
|
||||
" return [n for half in range(1, end)\n",
|
||||
" for n in palindromes(half)\n",
|
||||
" if isprime(n)]\n",
|
||||
"\n",
|
||||
"def palindromes(half: int) -> tuple[int]:\n",
|
||||
" \"\"\"Given the left half of a palindrome, make two full palindromes,\n",
|
||||
" one without the middle digit repeated, one with.\"\"\"\n",
|
||||
" left = str(half)\n",
|
||||
" right = left[::-1]\n",
|
||||
" return int(left + right[1:]), int(left + right)\n",
|
||||
"\n",
|
||||
"assert palindromes(123) == (12321, 123321)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "ec5f9d5a-c89e-4dc9-812b-232437e4801e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's count the palindromic primes of up to 11 digits, and peek at some of them:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"id": "3f9c6d25-8d40-4267-a45f-c0284b910c86",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"47995"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"PP = palindromic_primes(11)\n",
|
||||
"\n",
|
||||
"len(PP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"id": "69149eb4-d79d-4cfa-8e48-40946e54c48f",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[11,\n",
|
||||
" 2,\n",
|
||||
" 3,\n",
|
||||
" 5,\n",
|
||||
" 7,\n",
|
||||
" 101,\n",
|
||||
" 131,\n",
|
||||
" 151,\n",
|
||||
" 181,\n",
|
||||
" 191,\n",
|
||||
" 313,\n",
|
||||
" 353,\n",
|
||||
" 373,\n",
|
||||
" 383,\n",
|
||||
" 727,\n",
|
||||
" 757,\n",
|
||||
" 787,\n",
|
||||
" 797,\n",
|
||||
" 919,\n",
|
||||
" 929,\n",
|
||||
" 10301,\n",
|
||||
" 10501,\n",
|
||||
" 10601,\n",
|
||||
" 11311,\n",
|
||||
" 11411]"
|
||||
]
|
||||
},
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"PP[:25]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"id": "3bebaad1-8fb2-4892-87f8-8846062d92d4",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[99985958999,\n",
|
||||
" 99986068999,\n",
|
||||
" 99988688999,\n",
|
||||
" 99988988999,\n",
|
||||
" 99991819999,\n",
|
||||
" 99994349999,\n",
|
||||
" 99995759999,\n",
|
||||
" 99998189999,\n",
|
||||
" 99998989999,\n",
|
||||
" 99999199999]"
|
||||
]
|
||||
},
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"PP[-10:]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "853d5a6c-936e-4ab4-afc3-2082641063fc",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Digit Lengths\n",
|
||||
"\n",
|
||||
"How many digits are in the primes in all these lists? The function `digit_lengths` returns a {number_of_digits: count_of_primes_with_that_many_digits} dict, and plots a bar chart for that data:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"id": "28b0bc97-1b22-449f-b33a-6f8c84fd890d",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from collections import Counter\n",
|
||||
"\n",
|
||||
"def digit_lengths(primes) -> dict[int, int]:\n",
|
||||
" \"\"\"Plot a bar chart and return a dict of the number of digits in these primes.\"\"\"\n",
|
||||
" digits = Counter(len(str(p)) for p in primes)\n",
|
||||
" plt.bar(list(digits), list(digits.values()))\n",
|
||||
" plt.xlabel('Number of digits')\n",
|
||||
" plt.ylabel('Count')\n",
|
||||
" return dict(digits)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "09bc98cf-2f43-4f58-b414-f3eff2d834b3",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"First for the **left-truncatable primes**:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"id": "defd769e-9f92-4e4f-a12b-88a653b6efee",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{1: 4,\n",
|
||||
" 2: 11,\n",
|
||||
" 3: 39,\n",
|
||||
" 4: 99,\n",
|
||||
" 5: 192,\n",
|
||||
" 6: 326,\n",
|
||||
" 7: 429,\n",
|
||||
" 8: 521,\n",
|
||||
" 9: 545,\n",
|
||||
" 10: 517,\n",
|
||||
" 11: 448,\n",
|
||||
" 12: 354,\n",
|
||||
" 13: 276,\n",
|
||||
" 14: 212,\n",
|
||||
" 15: 117,\n",
|
||||
" 16: 72,\n",
|
||||
" 17: 42,\n",
|
||||
" 18: 24,\n",
|
||||
" 19: 13,\n",
|
||||
" 20: 6,\n",
|
||||
" 21: 5,\n",
|
||||
" 22: 4,\n",
|
||||
" 23: 3,\n",
|
||||
" 24: 1}"
|
||||
]
|
||||
},
|
||||
"execution_count": 16,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKBhJREFUeJzt3QlU1WX+x/EvyKKiuKWgiWi57+ZKpZkSpObo6JmpxjEqszI1l8mUyS2tMCs1HcypY+qcySxnsnLJxH1KXMLMNUY7eqQUKE1RS0W4//N9/v97/1wFQ7xwLw/v1zm/c/kt997n93CBD8/y+/k5HA6HAAAAWMrf2wUAAAAoToQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrBXi7AL4gNzdXTpw4IZUrVxY/Pz9vFwcAABSCXirw3LlzUqdOHfH3L7j9hrAjYoJOREREYeoVAAD4mLS0NKlbt26B+wk7IqZFx1lZoaGhJffdAQAARZaVlWUaK5x/xwtC2BFxdV1p0CHsAABQuvzWEBQGKAMAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsFuDtAgAouvoTVhf5ucdm9KHqAZQJtOwAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWC3A2wUAypr6E1YX+bnHZvTxaFkAoCygZQcAAFiNsAMAAKxGNxYAg+41ALaiZQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGpeDTtTp04VPz8/t6Vp06au/RcvXpThw4dLjRo1pFKlSjJw4EDJyMhwe43jx49Lnz59pGLFilKrVi0ZN26cXLlyxQtnAwAAfJHXr7PTokULWb9+vWs9IOD/izRmzBhZvXq1LF++XKpUqSIjRoyQAQMGyJdffmn25+TkmKATHh4u27Ztk5MnT8ojjzwigYGB8sorr3jlfAAAgG/xetjRcKNh5Wpnz56VhQsXytKlS6VHjx5m26JFi6RZs2ayfft26dKli6xbt04OHjxowlJYWJi0bdtWpk+fLuPHjzetRkFBQV44IwAA4Eu8Pmbn8OHDUqdOHbnttttk0KBBpltKpaSkSHZ2tkRHR7uO1S6uevXqSXJyslnXx1atWpmg4xQbGytZWVly4MCBAt/z0qVL5pi8CwAAsJNXw07nzp1l8eLFsnbtWnnrrbfk6NGj0rVrVzl37pykp6eblpmqVau6PUeDje5T+pg36Dj3O/cVJCEhwXSLOZeIiIhiOT8AAFDGu7F69erl+rp169Ym/ERGRsqHH34oFSpUKLb3jY+Pl7Fjx7rWtWWHwAMAgJ283o2Vl7biNG7cWI4cOWLG8Vy+fFnOnDnjdozOxnKO8dHHq2dnOdfzGwfkFBwcLKGhoW4LAACwk0+FnfPnz8t3330ntWvXlvbt25tZVRs2bHDtT01NNWN6oqKizLo+7tu3TzIzM13HJCUlmfDSvHlzr5wDAADwLV7txnruueekb9++puvqxIkTMmXKFClXrpw8/PDDZizNkCFDTHdT9erVTYAZOXKkCTg6E0vFxMSYUDN48GCZOXOmGaczceJEc20ebb0BAADwatj5/vvvTbA5deqU1KxZU+6++24zrVy/VrNnzxZ/f39zMUGdQaUzrebPn+96vgajVatWybBhw0wICgkJkbi4OJk2bZoXzwoAAPgSr4adZcuWXXd/+fLlJTEx0SwF0VahNWvWFEPpAACADXxqzA4AAICnEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVArxdAAD2qT9hdZGfe2xGH4+WBQBo2QEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAatwIFCgEbmwJAKUXLTsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAVvOZsDNjxgzx8/OT0aNHu7ZdvHhRhg8fLjVq1JBKlSrJwIEDJSMjw+15x48flz59+kjFihWlVq1aMm7cOLly5YoXzgAAAPginwg7u3btkr///e/SunVrt+1jxoyRlStXyvLly2XLli1y4sQJGTBggGt/Tk6OCTqXL1+Wbdu2yZIlS2Tx4sUyefJkL5wFAADwRV4PO+fPn5dBgwbJO++8I9WqVXNtP3v2rCxcuFBmzZolPXr0kPbt28uiRYtMqNm+fbs5Zt26dXLw4EH55z//KW3btpVevXrJ9OnTJTEx0QSggly6dEmysrLcFgAAYCevhx3tptLWmejoaLftKSkpkp2d7ba9adOmUq9ePUlOTjbr+tiqVSsJCwtzHRMbG2vCy4EDBwp8z4SEBKlSpYpriYiIKJZzAwAAZTzsLFu2THbv3m3Cx9XS09MlKChIqlat6rZdg43ucx6TN+g49zv3FSQ+Pt60HDmXtLQ0D50RAADwNQHeemMNGKNGjZKkpCQpX758ib53cHCwWQAAgP281rKj3VSZmZlyxx13SEBAgFl0EPLcuXPN19pCo+Nuzpw54/Y8nY0VHh5uvtbHq2dnOdedxwAAgLLNa2GnZ8+esm/fPtmzZ49r6dChgxms7Pw6MDBQNmzY4HpOamqqmWoeFRVl1vVRX0NDk5O2FIWGhkrz5s29cl4AAMC3eK0bq3LlytKyZUu3bSEhIeaaOs7tQ4YMkbFjx0r16tVNgBk5cqQJOF26dDH7Y2JiTKgZPHiwzJw504zTmThxohn0TDcVAADwatgpjNmzZ4u/v7+5mKBOF9eZVvPnz3ftL1eunKxatUqGDRtmQpCGpbi4OJk2bZpXyw0AAHyHT4WdzZs3u63rwGW9Zo4uBYmMjJQ1a9aUQOkAAEBp5PXr7AAAABQnwg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrBXi7AABwPfUnrC5yBR2b0YfKBUDLDgAAsBvdWAAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVArxdAKC41J+wusjPPTajj0fLAgDwHlp2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWK1LYue222+TUqVPXbD9z5ozZBwAAUKrDzrFjxyQnJ+ea7ZcuXZIffvjBE+UCAAAo+Ssof/rpp66vP//8c6lSpYprXcPPhg0bpH79+p4pGQAAQEmHnf79+5tHPz8/iYuLc9sXGBhogs4bb7zhiXIBAACUfNjJzc01jw0aNJBdu3bJLbfc4plSAAAA+NKYnaNHj3ok6Lz11lvSunVrCQ0NNUtUVJR89tlnrv0XL16U4cOHS40aNaRSpUoycOBAycjIcHuN48ePS58+faRixYpSq1YtGTdunFy5cuWmywYAAMr4Xc91fI4umZmZrhYfp3fffbdQr1G3bl2ZMWOGNGrUSBwOhyxZskT69esnX3/9tbRo0ULGjBkjq1evluXLl5vxQSNGjJABAwbIl19+6RonpEEnPDxctm3bJidPnpRHHnnEdKm98sorRT01AABQ1sPOiy++KNOmTZMOHTpI7dq1zRieoujbt6/b+ssvv2xae7Zv326C0MKFC2Xp0qXSo0cPs3/RokXSrFkzs79Lly6ybt06OXjwoKxfv17CwsKkbdu2Mn36dBk/frxMnTpVgoKCilQuAABQxsPOggULZPHixTJ48GCPFURbabQF58KFC6Y7KyUlRbKzsyU6Otp1TNOmTaVevXqSnJxswo4+tmrVygQdp9jYWBk2bJgcOHBA2rVrl+976RR5XZyysrI8dh4AAMCCMTuXL1+WO++80yMF2LdvnxmPExwcLE8//bSsWLFCmjdvLunp6aZlpmrVqm7Ha7DRfUof8wYd537nvoIkJCSYbjHnEhER4ZFzAQAAloSdJ554wnQveUKTJk1kz549smPHDtMio1PatWuqOMXHx8vZs2ddS1paWrG+HwAAKGXdWDpL6u233zZjZXQ2lQ4IzmvWrFmFfi1tvWnYsKH5un379mZK+5tvvikPPvigaUHSW1Dkbd3R2Vg6IFnp486dO91ezzlby3lMfrQVSRcAAGC/IoWdvXv3msHAav/+/W77ijpY2Ulndul4Gg0+GqJ0xpdOOVepqalmqrmO6VH6qIOadUaYTjtXSUlJZhq7doUBAAAUKexs2rTJY91JvXr1MoOOz507Z7rGNm/e7LoVxZAhQ2Ts2LFSvXp1E2BGjhxpAo4OTlYxMTEm1OhA6ZkzZ5pxOhMnTjTX5qHlBgAA3NR1djxBW2T0ujh6fRwNN9olpkHnvvvuM/tnz54t/v7+pmVHW3t0ptX8+fNdzy9XrpysWrXKjPXREBQSEmLG/Oi0eAAAgCKHnXvvvfe63VUbN24s1OvodXSup3z58pKYmGiWgkRGRsqaNWsK9X4AAKDsKVLLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 640x480 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"digit_lengths(TP)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "ec4c6f1d-d298-48de-884a-efee0b81fd19",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"And for the **right-truncatable** primes:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"id": "09f29951-d8f2-47a3-b666-bfa6fd5d1c87",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{1: 4, 2: 9, 3: 14, 4: 16, 5: 15, 6: 12, 7: 8, 8: 5}"
|
||||
]
|
||||
},
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJSxJREFUeJzt3Q1UVVX+//EvhoCZ4LNCglgpPoaaaGU+kKRjijquHmyZkfacz5ajTJmSFtCUWUmSrkmdKTNXiZmOmo8xTpqBmtkUSqEyPtGUgtCIBve/9v4v7o+roEjIOfvwfq21hXPOPZd97r3e+7l773O2l8vlcgkAAIChalldAQAAgN+DMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDRvcbji4mI5duyY1KtXT7y8vKyuDgAAqAB1GbwzZ85IUFCQ1KpVq2aHGRVkgoODra4GAACohOzsbGnRokXNDjOqRabkwfD397e6OgAAoALy8vJ0Y0TJ53iNDjMlXUsqyBBmAAAwS0WGiDAAGAAAGI0wAwAAjEaYAQAARiPMAAAAoxFmAACA0QgzAADAaIQZAABgNMIMAAAwGmEGAAAYjTADAACMRpgBAABGszTMpKamSnR0tJ7eW829sGrVqotu891338mQIUMkICBA6tatKxEREXLkyBFL6gsAAOzH0jBTUFAg4eHhkpSUVOb2H374Qe644w5p27atbNu2Tfbt2yczZswQPz+/aq8rAACwJy+Xy+USG1AtMykpKTJs2DD3uhEjRkjt2rXl73//+++aQly16uTm5jJrNgAAhriSz2/bjpkpLi6WtWvXSps2bWTAgAHStGlT6dGjR5ldUaUVFhbqB6B0AQAAzuUtNpWTkyP5+fmSkJAgc+bMkcTERFm/fr0MHz5ctm7dKn369Clzv/j4eImLi6v2+gKmCJ2+Vkx0KGGQ1VUAYFO2bplRhg4dKpMnT5bOnTvL9OnTZfDgwZKcnFzufrGxsbpJqqRkZ2dXY60BAEB1s23LTOPGjcXb21vat2/vsb5du3ayffv2cvfz9fXVBQAA1Ay2bZnx8fHRp2FnZGR4rD9w4IC0bNnSsnoBAAB7sbRlRo2JyczMdC9nZWXJ3r17pWHDhhISEiJTp06V+++/X3r37i2RkZF6zMynn36qT9MGAACwPMykpaXpkFJiypQp+mdMTIwsWbJE/vjHP+rxMWpQ74QJEyQsLEw+/vhjfe0ZAAAAy8NM37595XKXuRkzZowuAAAARo2ZAQAAqAjCDAAAMBphBgAAGI0wAwAAjEaYAQAARiPMAAAAoxFmAACA0QgzAADAaIQZAABgNMIMAAAwGmEGAAAYjTADAACMRpgBAABGI8wAAACjEWYAAIDRvK2uAABcDaHT1xr5wB5KGGR1FQDj0DIDAACMRpgBAABGI8wAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAAGI0wAwAAjGZpmElNTZXo6GgJCgoSLy8vWbVqVbm3ffLJJ/Vt5s2bV611BAAA9mZpmCkoKJDw8HBJSkq65O1SUlJk586dOvQAAACU5i0WGjhwoC6XcvToURk/frxs2LBBBg0aVG11AwAAZrA0zFxOcXGxjBo1SqZOnSodOnSo0D6FhYW6lMjLy7uKNQQAAFazdZhJTEwUb29vmTBhQoX3iY+Pl7i4uKtaLzhT6PS1YqJDCbRYAqjZbHs2U3p6urzxxhuyZMkSPfC3omJjYyU3N9ddsrOzr2o9AQCAtWwbZv75z39KTk6OhISE6NYZVQ4fPizPPPOMhIaGlrufr6+v+Pv7exQAAOBctu1mUmNloqKiPNYNGDBArx89erRl9QIAAPZiaZjJz8+XzMxM93JWVpbs3btXGjZsqFtkGjVq5HH72rVrS/PmzSUsLMyC2gIAADuyNMykpaVJZGSke3nKlCn6Z0xMjB4rAwAAYOsw07dvX3G5XBW+/aFDh65qfQAAgHlsOwAYAACgIggzAADAaIQZAABgNMIMAAAwGmEGAAAYjTADAACMRpgBAABGI8wAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAAGI0wAwAAjEaYAQAARiPMAAAAoxFmAACA0QgzAADAaIQZAABgNMIMAAAwGmEGAAAYjTADAACMRpgBAABGI8wAAACjEWYAAIDRLA0zqampEh0dLUFBQeLl5SWrVq1ybzt//rxMmzZNOnXqJHXr1tW3eeihh+TYsWNWVhkAANiMpWGmoKBAwsPDJSkp6aJtv/76q+zevVtmzJihf65cuVIyMjJkyJAhltQVAADYk7eVf3zgwIG6lCUgIEA2btzosW7+/PnSvXt3OXLkiISEhFRTLQEAgJ1ZGmauVG5uru6Oql+/frm3KSws1KVEXl5eNdUOAABYwZgwc/bsWT2G5oEHHhB/f/9ybxcfHy9xcXHVWjcAsEro9LVGPviHEgZZXQU4iBFnM6nBwPfdd5+4XC5ZsGDBJW8bGxurW3BKSnZ2drXVEwAAVD9vU4LM4cOHZcuWLZdslVF8fX11AQAANYO3CUHm4MGDsnXrVmnUqJHVVQIAADZjaZjJz8+XzMxM93JWVpbs3btXGjZsKIGBgXLPPffo07LXrFkjRUVFcuLECX07td3Hx8fCmgMAALuwNMykpaVJZGSke3nKlCn6Z0xMjMyaNUtWr16tlzt37uyxn2ql6du3bzXXFgAA2JGlYUYFEjWotzyX2gYAAGDM2UwAAADlIcwAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAAGI0wAwAAjEaYAQAARiPMAAAAoxFmAACA0QgzAADAaIQZAABgNMIMAAAwGmEGAAAYjTADAACMRpgBAABGI8wAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEazNMykpqZKdHS0BAUFiZeXl6xatcpju8vlkhdeeEECAwOlTp06EhUVJQcPHrSsvgAAwH4sDTMFBQUSHh4uSUlJZW5/5ZVX5M0335Tk5GT58ssvpW7dujJgwAA5e/ZstdcVAADYk7eVf3zgwIG6lEW1ysybN0+ef/55GTp0qF73t7/9TZo1a6ZbcEaMGFHNtQUAAHZk2zEzWVlZcuLECd21VCIgIEB69OghO3bsKHe/wsJCycvL8ygAAMC5bBtmVJBRVEtMaWq5ZFtZ4uPjdegpKcHBwVe9rgAAwDq2DTOVFRsbK7m5ue6SnZ1tdZUAAEBNDDPNmzfXP0+ePOmxXi2XbCuLr6+v+Pv7exQAAOBctg0zrVq10qFl8+bN7nVq/Is6q+m2226ztG4AAMA+LD2bKT8/XzIzMz0G/e7du1caNmwoISEhMmnSJJkzZ460bt1ah5sZM2boa9IMGzbMymoDAAAbsTTMpKWlSWRkpHt5ypQp+mdMTIwsWbJE/vSnP+lr0Tz++ONy+vRpueOOO2T9+vXi5+dnYa0BAICdWBpm+vbtq68nUx51VeAXX3xRFwAAAKPGzAAAAFQEYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAA1Lwwc8MNN8jPP/980Xo1f5LaBgAAYOswc+jQISkqKrpofWFhoRw9erQq6gUAAFD1E02uXr3a/fuGDRskICDAvazCzebNmyU0NPRK7hIAAKD6wsywYcPcs1nHxMR4bKtdu7YOMq+99trvqxEAAMDVCjPFxcX6Z6tWreSrr76Sxo0bX8nuAAAA1oaZEllZWVVfE9hW6PS1YqpDCYOsrgIAwI5hRlHjY1TJyclxt9iUePfdd6uibgAAAFcnzMTFxcmLL74o3bp1k8DAQD2GBgAAwJgwk5ycLEuWLJFRo0ZVfY0AAACu9nVmzp07J7fffntldgUAALA+zDz66KOybNmyqq0JAABAdXUznT17VhYuXCibNm2Sm2++WV9jprS5c+dW5m4BAACqJ8zs27dPOnfurH/fv3+/xzYGAwMAgOpUqTCzdevWqq8JAABAdY2ZAQAAMLplJjIy8pLdSVu2bPk9dQIAALi6YaZkvEyJ8+fPy969e/X4mQsnoAQAALBdmHn99dfLXD9r1izJz8//vXUCAACwZszMgw8+yLxMAADA3DCzY8cO8fPzq8q7BAAAqPpupuHDh3ssu1wuOX78uKSlpcmMGTMqc5cAAADVF2YCAgI8lmvVqiVhYWF6Ju3+/ftXriYAAADVFWYWL15cmd0AAADsNWYmPT1d3nvvPV327NkjVa2oqEh3W7Vq1Urq1KkjN954o8yePVt3awEAAFS6ZSYnJ0dGjBgh27Ztk/r16+t1p0+f1hfTW758uTRp0qRKHt3ExERZsGCBLF26VDp06KDH5IwePVp3c02YMIFnEAAAVK5lZvz48XLmzBn59ttv5ZdfftFFXTAvLy+vSkPGF198IUOHDpVBgwZJaGio3HPPPXpMzq5du3jqAABA5cPM+vXr5e2335Z27dq517Vv316SkpJk3bp1UlVuv/122bx5sxw4cEAvf/3117J9+3YZOHBgufsUFhbqUFW6AAAA56pUN1NxcbHUrl37ovVqndpWVaZPn67DSNu2beWaa67RY2heeuklGTlyZLn7xMfHS1xcXJXVAQBgrdDpa419Cg4lDLK6CjVCpVpm7rzzTpk4caIcO3bMve7o0aMyefJk6devX5VVbsWKFfL+++/LsmXLZPfu3XrszKuvvqp/lic2NlZyc3PdJTs7u8rqAwAAHNIyM3/+fBkyZIgexxILine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 640x480 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"digit_lengths(TPr)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8618e5bc-ea92-4116-aa74-2380fbfd0dc3",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Finally, for the **palindromic primes**, we find that \"11\" is the only even-digit-length palindromic prime. (It turns out that any even-digit-length palindromic number is [divisible by 11](https://www.geeksforgeeks.org/dsa/g-fact-all-palindromes-with-an-even-number-of-digits-are-divisible-by-11/).) The counts of primes grow so quickly that it is better to use a log scale on the y axis: "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"id": "ed392e12-7480-4d96-994f-96cdb786f8f9",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{2: 1, 1: 4, 3: 15, 5: 93, 7: 668, 9: 5172, 11: 42042}"
|
||||
]
|
||||
},
|
||||
"execution_count": 18,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGwCAYAAACgi8/jAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJJdJREFUeJzt3QuUVdV9P/DfAAJqBEVSEAExVlF8gOEl5iWRSImSqm1iuhKD2Jh0FSIJ0RZihJhowTQSGkNKtUXa1VqJbUQjlahIglqUl2gMYqSFhoi8jPKy8pr5r326Zv6OPISZgXvZ9/NZ6yzuOffOnX23I/Nl79/ep6qmpqYmAAAy1KzUDQAAOFQEHQAgW4IOAJAtQQcAyJagAwBkS9ABALIl6AAA2WoRFa66ujrWrFkTxx13XFRVVZW6OQDAAUjbAG7ZsiU6deoUzZrte9ym4oNOCjldunQ5kD4FAMrM6tWro3Pnzvt8vmKDzpQpU4pj165ddR3Vpk2bUjcLADgAmzdvLgYq0ozM/lRV+i0gUke1bds2Nm3aJOgAQGa/vxUjAwDZEnQAgGwJOgBAtio26KRC5B49ekTfvn1L3RQA4BBRjKwYGQCOOIqRAYCKV7FTVwBA/gQdACBbFRt0FCMDQP4UIytGBoAjjmJkAKDiVezUFQCQP0EHAMiWoAMAZKtFqRsAABwa3cbMKnnXrpp4aUm/f8WO6FheDgD5q9igM2LEiFi2bFksXLiw1E0BAA6Rig06AED+BB0AIFuCDgCQLUEHAMiWoAMAZEvQAQCyJegAANkSdACAbFVs0LEzMgDkr2KDjp2RASB/FRt0AID8CToAQLYEHQAgW4IOAJAtQQcAyJagAwBkS9ABALIl6AAA2RJ0AIBsCToAQLYqNui41xUA5K9ig457XQFA/io26AAA+RN0AIBsCToAQLYEHQAgW4IOAJAtQQcAyJagAwBkS9ABALIl6AAA2RJ0AIBsCToAQLYEHQAgW4IOAJAtQQcAyJagAwBkK5ug89Zbb8Upp5wSN9xwQ6mbAgCUiWyCzm233RYXXHBBqZsBAJSRLILOK6+8EsuXL48hQ4aUuikAQBkpedCZN29eDB06NDp16hRVVVUxc+bMPV4zZcqU6NatW7Ru3Tr69+8fCxYsqPd8mq6aMGHCYWw1AHAkKHnQ2bZtW/Ts2bMIM3szY8aMGD16dIwfPz6WLFlSvHbw4MGxfv364vkHH3wwzjjjjOI4ENu3b4/NmzfXOwCAPLUodQPSdNP+ppwmTZoU1113XQwfPrw4nzp1asyaNSumTZsWY8aMiWeeeSbuu+++uP/++2Pr1q2xc+fOaNOmTYwbN26v75dGfm655ZZD9nkAgPJR8hGd/dmxY0csXrw4Bg0aVHetWbNmxfn8+fPrgsvq1atj1apV8b3vfa8IRfsKOcnYsWNj06ZNdUf6WgAgTyUf0dmfjRs3xu7du6NDhw71rqfzVHzcEK1atSoOACB/ZR10DtY111xT6iYAAGWkrKeu2rdvH82bN49169bVu57OO3bs2Kj3TsXPPXr0iL59+zaylQBAuSrroNOyZcvo3bt3zJkzp+5adXV1cT5gwIBGvfeIESNi2bJlsXDhwiZoKQBQjko+dZVWSq1YsaLufOXKlbF06dJo165ddO3atVhaPmzYsOjTp0/069cvJk+eXCxJr12FBQBQtkFn0aJFMXDgwLrzFGySFG6mT58eV111VWzYsKFYSbV27dro1atXzJ49e48C5YZMXaUjFTsDAHmqqqmpqYkKljYMbNu2bbHUPO2/AwC56DZmVqmbEKsmXlrS399lXaMDANAYgg4AkC1BBwDIVsUGHfvoAED+Kjbo2EcHAPJXsUEHAMifoAMAZKtig44aHQDIX8UGHTU6AJC/ig06AED+BB0AIFuCDgCQrZLfvRwAyukmlIfyRpQcfhU7omPVFQDkr2KDjlVXAJC/ig06AED+BB0AIFuCDgCQLUEHAMiWoAMAZKtig47l5QCQv4oNOpaXA0D+KjboAAD5E3QAgGwJOgBAtgQdACBbgg4AkC1BBwDIlqADAGSrYoOODQMBIH8VG3RsGAgA+avYoAMA5E/QAQCyJegAANkSdACAbAk6AEC2BB0AIFuCDgCQLUEHAMiWoAMAZEvQAQCyJegAANmq2KDjpp4AkL+KDTpu6gkA+avYoAMA5E/QAQCyJegAANkSdACAbAk6AEC2BB0AIFuCDgCQLUEHAMiWoAMAZEvQAQCyJegAANkSdACAbAk6AEC2BB0AIFuCDgCQrSM+6Lz55pvRp0+f6NWrV5xzzjlx9913l7pJAECZaBFHuOOOOy7mzZsXxxxzTGzbtq0IO1deeWWceOKJpW4aAFBiR/yITvPmzYuQk2zfvj1qamqKAwCg5EEnjcYMHTo0OnXqFFVVVTFz5sw9XjNlypTo1q1btG7dOvr37x8LFizYY/qqZ8+e0blz57jxxhujffv2h/ETAADlquRBJ003pZCSwszezJgxI0aPHh3jx4+PJUuWFK8dPHhwrF+/vu41xx9/fDz//POxcuXKuPfee2PdunX7/H5p1Gfz5s31DgAgTyUPOkOGDIlbb701rrjiir0+P2nSpLjuuuti+PDh0aNHj5g6dWoxVTVt2rQ9XtuhQ4ciCD355JP7/H4TJkyItm3b1h1dunRp0s8DAJSPkged/dmxY0csXrw4Bg0aVHetWbNmxfn8+fOL8zR6s2XLluLxpk2biqmw7t277/M9x44dW7yu9li9evVh+CQAQCmU9aqrjRs3xu7du4uRmndK58uXLy8e/8///E986UtfqitC/spXvhLnnnvuPt+zVatWxQEA5K+sg86B6NevXyxduvSgvy7VBKUjBSkAIE9lPXWVVk+l5ePvLi5O5x07dmzUe48YMSKWLVsWCxcubGQrAYByVdZBp2XLltG7d++YM2dO3bXq6urifMCAASVtGwBQ/ko+dbV169ZYsWJF3XlaIp6motq1axddu3YtlpYPGzasuM1DmqaaPHlysSQ9rcICACjroLNo0aIYOHBg3XkKNkkKN9OnT4+rrroqNmzYEOPGjYu1a9cW97SaPXv2HgXKB0uNDgDkr6qmwu+XkDYMTPvppKXmbdq0KXVzACpatzGzohysmnhp5KBbGfTnoerLA/39XdY1OgAAjSHoAADZqtigk2p00i0l+vbtW+qmAACHSMUGHfvoAED+Sr7qCuBIl3PBJxzpKnZEBwDIn6ADAGSrYoOOYmQAyF/FBh3FyACQv4oNOgBA/gQdACBbgg4AkC1BBwDIVsUGHauuACB/FRt0rLoCgPxVbNABAPIn6AAA2RJ0AIBsCToAQLYEHQAgWxUbdCwvB4D8VWzQsbwcAPJXsUEHAMhfg4LOBz7wgXj99df3uP7mm28WzwEAHLFBZ9WqVbF79+49rm/fvj1effXVpmgXAECjtTiYFz/00EN1j3/2s59F27Zt685T8JkzZ05069at8a0CADjcQefyyy8v/qyqqophw4bVe+6oo44qQs4dd9zRFO0CADi8Qae6urr489RTT42FCxdG+/btG98CAIByCDq1Vq5c2fQtAQAoh6CTpHqcdKxfv75upKfWtGnT4kjYMDAdeyuqBgAqeNXVLbfcEpdcckkRdDZu3BhvvPFGveNIYMNAAMhfg0Z0pk6dGtOnT4+rr7666VsEAFDKEZ0dO3bEhRde2FRtAAAon6DzxS9+Me69996mbw0AQKmnrt5+++2466674vHHH4/zzjuv2EPnnSZNmtRU7QMAOLxB54UXXohevXoVj1988cV6z6XNBAEAjtigM3fu3KZvCQBAOdToAABkO6IzcODA/U5RPfHEE41pEwBA6YJObX1OrZ07d8bSpUuLep133+wTAOCICjrf//7393r9W9/6VmzdurWxbQIAKL8anc9//vNHxH2uAIDK0KRBZ/78+dG6des4EqQbevbo0SP69u1b6qYAAOU0dXXllVfWO6+pqYnXXnstFi1aFDfffHMcKTf1TMfmzZujbdu2pW4OAFAuQefdwaBZs2bRvXv3+Pa3v13c1RwA4IgNOvfcc0/TtwQAoByCTq3FixfHSy+9VDw+++yz4/zzz2+qdgEAlCborF+/Pj772c/Gz3/+8zj++OOLa2+++WaxkeB9990X73//+xvfMgCAUqy6+spXvhJbtmyJX/3qV/G73/2uONJmgamw9/rrr29smwAASjeiM3v27Hj88cfjrLPOqruWlmqnJduKkQGAI3pEp7q6Oo466qg9rqdr6TkAgCM26Hz84x+PUaNGxZo1a+quvfrqq/G1r30tLr744qZsHwDA4Q06P/zhD4t6nG7dusVpp51WHKeeempx7c4772x4awAASl2j06VLl1iyZElRp7N8+fLiWqrXGTRoUFO2DQDg8I3oPPHEE0XRcRq5qaqqik984hPFCqx0pHtGpb10nnzyyca1CACgFEFn8uTJcd1110WbNm32eluIL3/5yzFp0qSmahsAwOELOs8//3z8wR/8wT6fT0vL027JAABHXNBZt27dXpeV12rRokVs2LChKdoFAHB4g87JJ59c7IC8Ly+88EKcdNJJcTitXr06LrrooqJ26Lzzzov777//sH5/ACCToPPJT34ybr755nj77bf3eO5///d/Y/z48XHZZZfF4ZRGkVLt0LJly+LRRx+Nr371q7Ft27bD2gYAIIPl5d/85jfjJz/5SZxxxhkxcuTI6N69e3E9LTFPt3/YvXt33HTTTXE4pRGk2lGkjh07Rvv27Yt7bx177LGHtR0AwBE+otOhQ4f4z//8zzjnnHNi7NixccUVVxTHN77xjeLaU089VbzmYMybNy+GDh0anTp1Kpasz5w5c4/XpBCVNids3bp19O/fPxYsWLDX90qF0ClspX1+AAAOemfkU045Jf7jP/4jNm7cGM8++2w888wzxeN0Le2OfLDSNFPPnj2LMLM3M2bMiNGjRxfTYmmLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 640x480 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"plt.yscale('log')\n",
|
||||
"digit_lengths(PP)"
|
||||
"[p for p in TPr if is_omni_truncatable_prime(p)]"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -632,31 +294,16 @@
|
||||
"id": "f2323985-2518-41d1-8a3b-6e21e40d27be",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Summary \n",
|
||||
"There are only 9 of them: 4 one-digit, 4 two-digit, and 1 three-digit.\n",
|
||||
"\n",
|
||||
"Here's what we learned:\n",
|
||||
"\n",
|
||||
"- Left-trunctable primes:\n",
|
||||
" - Count: 4260\n",
|
||||
" - Largest: 357686312646216567629137 (24 digits)\n",
|
||||
" - The plurality have 9 digits; only 13 have more than 20 digits\n",
|
||||
"- Right-truncatable primes:\n",
|
||||
" - Count: 83\n",
|
||||
" - Largest: 73939133 (8 digits)\n",
|
||||
" - The plurality have 4 digits\n",
|
||||
"- Two-sided truncatable primes:\n",
|
||||
" - Count: 15\n",
|
||||
" - Largest: 739397 (6 digits)\n",
|
||||
"- Omni-truncatable primes:\n",
|
||||
" - Count: 9\n",
|
||||
" - Largest: 373\n",
|
||||
"- Palindromic truncatable primes:\n",
|
||||
" - Count: 13\n",
|
||||
" - Largest: 799636997 (9 digits)\n",
|
||||
"- Palindromic primes:\n",
|
||||
" - Count: probably infinite; 47995 of them up to 11 digits\n",
|
||||
" - All except \"11\" have an odd number of digits\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",
|
||||
@@ -665,7 +312,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"execution_count": 10,
|
||||
"id": "1d460931-3d70-485d-807f-8fb360ef4ef2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
@@ -684,12 +331,12 @@
|
||||
"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 bigger. The small 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 bigger 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."
|
||||
"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": 20,
|
||||
"execution_count": 11,
|
||||
"id": "647fb66b-32ae-494a-8945-cf7358c2f958",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
@@ -702,7 +349,7 @@
|
||||
" return n in (2, 3, 5, 7)\n",
|
||||
" elif n % 2 == 0: # Even numbers other than 2 are composite\n",
|
||||
" return False\n",
|
||||
" else: # test odd divisors up to sqrt(n)\n",
|
||||
" else: # test odd divisors up to √n\n",
|
||||
" divisors = range(11, int(sqrt(n)) + 1, 2) \n",
|
||||
" return not any(n % d == 0 for d in divisors)"
|
||||
]
|
||||
@@ -712,11 +359,11 @@
|
||||
"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. 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",
|
||||
"\n",
|
||||
"**Fortunately**, [Pierre de Fermat](https://en.wikipedia.org/wiki/Pierre_de_Fermat) provided the big breakthrough in 1640! \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",
|
||||
@@ -731,28 +378,29 @@
|
||||
"|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",
|
||||
"|5|2|1|*5 could be prime or composite*|\n",
|
||||
"|5|3|1|*try again: 5 could still be prime or composite*|\n",
|
||||
"|5|4|1|*try again: 5 could still be prime or composite, but a lot of evidence that it is prime*|\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",
|
||||
"\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. Here is a simple implementation:"
|
||||
"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:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"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 False then n is definitely composite; if True then n is probably prime.\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 for a in sample(2, n, k))"
|
||||
" return all(pow(a, n - 1, n) == 1 \n",
|
||||
" for a in sample(2, n, k))"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -767,14 +415,14 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"execution_count": 13,
|
||||
"id": "476282dd-a48a-4eb5-8fed-8495d4878ec2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import random\n",
|
||||
"from typing import Iterable\n",
|
||||
" \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))"
|
||||
@@ -785,62 +433,51 @@
|
||||
"id": "15ef6884-253d-4a3a-b58f-397ac73f552d",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**Unfortunately**, `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 test a sample of numbers and report the the percentage of false primes when `isprime_fermat` is allowed `k` choices of the `a` parameter:"
|
||||
"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*:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 23,
|
||||
"execution_count": 22,
|
||||
"id": "45a11a93-9ef3-48d1-80fa-e0113041346e",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from statistics import mean\n",
|
||||
"\n",
|
||||
"def false_prime_percent(samples: Iterable[int], k: int) -> float:\n",
|
||||
" \"\"\"The 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 mean(100 * isprime_fermat(n, k) for n in composites)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "778c3a9b-dcc7-4ad2-9092-d66aecd1d93c",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We'll choose sample ranges that represent digit-lengths, sample 200,000 times for each digit-length, and for now look at `k=1` values of *a*:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 24,
|
||||
"id": "9eba96df-c298-4d71-b68a-3bf36bf15095",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{3: 1.408066267526657,\n",
|
||||
" 4: 0.4848670443525749,\n",
|
||||
" 5: 0.1327436065501534,\n",
|
||||
" 6: 0.03568377685744872,\n",
|
||||
" 7: 0.011224724326115657,\n",
|
||||
" 8: 0.0015893702914905115,\n",
|
||||
" 9: 0,\n",
|
||||
" 10: 0,\n",
|
||||
" 11: 0,\n",
|
||||
" 12: 0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 24,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
"name": "stdout",
|
||||
"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",
|
||||
"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": [
|
||||
"{d: false_prime_percent(sample(10**(d-1), 10**d, 200_000), k=1) for d in range(3, 13)}"
|
||||
"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=200_000):\n",
|
||||
" \"\"\"Report false prime probability for different digit lengths, and k values.\"\"\"\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)])"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -848,93 +485,52 @@
|
||||
"id": "f3fe25d6-5f0a-478b-90d0-6c89e092e92e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"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 9- or 10-digit composites, false primes are very rare indeed. Remember these are percents, not probabilities, so 0.001% is one in a hundred thousand.\n",
|
||||
"\n",
|
||||
"Let's see how we reduce false primes by allowing `k=25` choices of `a`:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 25,
|
||||
"id": "8880c588-32f5-46e8-b59c-f9b49d492d3b",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{3: 0,\n",
|
||||
" 4: 0,\n",
|
||||
" 5: 0.00036752274965820385,\n",
|
||||
" 6: 0.0007220216606498195,\n",
|
||||
" 7: 0,\n",
|
||||
" 8: 0,\n",
|
||||
" 9: 0,\n",
|
||||
" 10: 0,\n",
|
||||
" 11: 0,\n",
|
||||
" 12: 0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 25,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{d: false_prime_percent(sample(10**(d-1), 10**d, 300_000), k=25) for d in range(3, 13)}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "61ee44c3-f963-4fb1-ac20-4a27adc90f83",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The results are randomized, and in some runs this produces no false primes; in some runs a few get through. \n",
|
||||
"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",
|
||||
"\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",
|
||||
"\n",
|
||||
"We can test some of them:"
|
||||
"We can test some of them:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 26,
|
||||
"execution_count": 23,
|
||||
"id": "7871aaa7-af1c-462b-85bb-d46818ccabea",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{561: 57.14,\n",
|
||||
" 1105: 69.55,\n",
|
||||
" 1729: 74.87,\n",
|
||||
" 2465: 72.6,\n",
|
||||
" 2821: 76.78,\n",
|
||||
" 6601: 80.01,\n",
|
||||
" 8911: 79.49,\n",
|
||||
" 41041: 70.78,\n",
|
||||
" 101101: 71.57,\n",
|
||||
" 294409: 95.12,\n",
|
||||
" 340561: 82.24,\n",
|
||||
" 825265: 60.79,\n",
|
||||
" 56052361: 99.04,\n",
|
||||
" 118901521: 99.29,\n",
|
||||
" 172947529: 99.32,\n",
|
||||
" 216821881: 99.31,\n",
|
||||
" 1299963601: 99.58,\n",
|
||||
" 2301745249: 99.88,\n",
|
||||
" 9624742921: 99.8}"
|
||||
]
|
||||
},
|
||||
"execution_count": 26,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
"name": "stdout",
|
||||
"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"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"carmichael_numbers = (561, 1105, 1729, 2465, 2821, 6601, 8911, 41041, 101101, 294409, 340561, 825265, \n",
|
||||
"carmichael_numbers = (1105, 1729, 2465, 2821, 6601, 8911, 41041, 101101, 294409, 340561, 825265, \n",
|
||||
" 56052361, 118901521, 172947529, 216821881, 1299963601, 2301745249, 9624742921)\n",
|
||||
"\n",
|
||||
"{n: false_prime_percent([n] * 10_000, k=1) for n in carmichael_numbers}"
|
||||
"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)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -942,46 +538,7 @@
|
||||
"id": "4ecd2289-66d5-4b0f-ac00-5db0ad92dde9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This is bad; some error rates are above 99%. Increasing `k` can fix the numbers with a 60% or 70% error rate, but not the ones with a 99% error rate. Below we increase `k` to 50, and some Carmichael numbers still claim to be prime most of the time:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 27,
|
||||
"id": "de492047-77b6-46a3-8a77-46308f219efd",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{561: 0,\n",
|
||||
" 1105: 0,\n",
|
||||
" 1729: 0,\n",
|
||||
" 2465: 0,\n",
|
||||
" 2821: 0,\n",
|
||||
" 6601: 0,\n",
|
||||
" 8911: 0,\n",
|
||||
" 41041: 0,\n",
|
||||
" 101101: 0,\n",
|
||||
" 294409: 9,\n",
|
||||
" 340561: 0,\n",
|
||||
" 825265: 0,\n",
|
||||
" 56052361: 61,\n",
|
||||
" 118901521: 64,\n",
|
||||
" 172947529: 80,\n",
|
||||
" 216821881: 77,\n",
|
||||
" 1299963601: 85,\n",
|
||||
" 2301745249: 85,\n",
|
||||
" 9624742921: 96}"
|
||||
]
|
||||
},
|
||||
"execution_count": 27,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{n: false_prime_percent([n] * 100, k=50) for n in carmichael_numbers}"
|
||||
"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."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -991,25 +548,26 @@
|
||||
"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 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` 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 counterexamples found so far (no equivalent of the Carmichael numbers for the Baillie-PSW test), but no proof that such numbers don't exist.\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",
|
||||
"\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 in ascending order of digit-length from 5 to 32:"
|
||||
"**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": 28,
|
||||
"id": "be5f2fe7-b1e0-4a02-b7f0-7903e865c4b8",
|
||||
"execution_count": 16,
|
||||
"id": "2a10367d-45c1-4425-a630-96892555c1b7",
|
||||
"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]"
|
||||
"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]"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -1017,16 +575,17 @@
|
||||
"id": "6b2cbdc7-e8a9-4b83-aa0d-1e7cc83de071",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Now let's make a chart of all the primality checking functions running on the big primes:"
|
||||
"Next make a chart of the speeds of all the primality checking functions running on the big primes:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 29,
|
||||
"execution_count": 17,
|
||||
"id": "8109d9fc-bbae-4ef0-a75c-908d9746a557",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import timeit\n",
|
||||
"\n",
|
||||
"prime_tests = (isprime_simple, isprime_faster, isprime_fermat, isprime)\n",
|
||||
@@ -1061,13 +620,13 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 30,
|
||||
"execution_count": 18,
|
||||
"id": "b1ab0f3d-6ce6-4b8c-9d4d-bb78bb1ac432",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAINCAYAAACTVPhqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAvs5JREFUeJzs3Qd8U1X7B/BfunfphpayUdkbBERkCcgQRMHJcL/iK0MF9VXcDEUEFf+4cOBEFBVkgzJkT0Fk70JbCnTvJP/Pc9KUtnQkbdKs39fPNTfJzc3JTVvuc89znqPR6/V6EBERERERlcKttAeJiIiIiIgYMBARERERUbnYw0BERERERGViwEBERERERGViwEBERERERGViwEBERERERGViwEBERERERGViwEBERERERGXyKPspEjqdDufPn0dgYCA0Gg0PChERERE5PJm7OS0tDdHR0XBzK78PgQFDGebOnauW3NxcHD9+3BrfExERERGRTZ09exa1a9cudxuNXsILKlNKSgpq1KihDmZQUBCPFBERERE5vNTUVMTGxiI5ORnBwcHlbssehgoY05AkWGDAQERERETOxJSUew56JiIiIiKiMjFgICIiIiKiMjFgICIiIiKiMnEMAxEREZENSN2Z/Px8aLVaHn+yCk9PT7i7u1d5PwwYiIiIiKqZlG2/cOECMjMzeezJqgOapWRqQEBAlfbDgKGCeRgY9RMREZGlJ4U9efKkuvIrk2Z5eXlxcliySg/WxYsXce7cOTRu3LhKPQ2ch8GEGrVSm1bmY2BZVSIiIqqq7OxsFTDUrVsXfn5+PKBkNVlZWTh16hTq168PHx+fSp/jctAzERERkQ24ufE0jGw/x4Ip+JNKRERERERlYsBARERERERlYsBARERE5KC0Oj22HL+EX/fGqVu5b0233HILxo8fb7X9S769pNHs3bsX9uhUNbavXr16mD17NuwBqyQREREROaAVBy7g1SUHcSElu/CxWsE+eHlQU/RrXssq7/nzzz+r2v7WEhsbq8rNhoeHwx7F2nn7rIU9DEREREQOGCz85+vdxYIFEZ+SrR6X560hNDQUgYGBVpubQkp/1qxZEx4e9nlN293O22ctDBjI4VR39ysREVF11MzPzM03aUnLzsPLv/2D0v71Mz72ym8H1Xam7E/euzIpSR9++KGq7y/lOqOionDnnXcW2+7JJ59Ui5TulCvyL730UrH3kpSb119/HSNHjlRlPR999NFrUn7+/PNPdX/lypVo06YNfH190bNnTyQmJmL58uVo0qSJeu29995bbBI8meti2rRpqpyovKZVq1ZYtGiRSZ/xypUruO+++xAREaFe27hxY3z++efqOUu1z5TjU1JycjIefvhh1S7Zp7zPvn37UB1cKzwih2eL7lciIiJry8rToumUlRbZl5xyxqdmo8Urq0za/uBrfeHnZd4p4c6dO/HUU09hwYIF6NKlCy5fvoyNGzcW2+bLL7/EQw89hO3bt6vtJSCoU6cOHnnkkcJtZs6ciSlTpuDll18u9/1eeeUVfPDBB2reiuHDh6vF29sb3377LdLT0zF06FC8//77mDx5stpegoWvv/4a8+bNUyf8GzZswP33369Otrt3717ue8mJ+8GDB9UJv5zIHzt2TM1nYMn2mXp8irrrrrtUQCLtkiDjo48+Qq9evXDkyBHV82NNDBjI4bpfS8bexu7X/7u/LYMGIiKianDmzBn4+/tj4MCBKkVJJqGTK+wl8/3fffdddQX++uuvx/79+9X9oifEcpX86aefLrwvV/BL88Ybb6Br165qXU6yn3/+eRw/fhwNGjRQj0nvxh9//KFOyHNycjB16lSsWbMGnTt3Vs/Ldps2bVIn2RUFDPLZ5LO0b9++sCekIua0z5zjYyRtl8BCei4kEDEGW7/88ovqOZFgw5oYMJRh7ty5atFqtVb9Asg0knYkPQtldb/KtCTyfJ+mNeHuZplJSoicnk4LnN4MpCcAAVFA3S6Am7utW0Xkknw93dWVflNsP3kZoz/fUeF2X4zpgI71Q016b3P16dNHBQlyQtyvXz+1yFX0ojNX33jjjcUmDpOT93feeUedW8lYAGE8Ka9Iy5YtC9cl/Unex3gybnxMTqiF9AhI+o+0seQYiZJBTWn+85//YNiwYdi9ezduvfVWDBkyRPWiWKp95hwfI0k9kp6KsLCwYo9Lz4cEJtbGgKEMY8eOVYtx2myyLfnjWHJgV8mgQZ6X7To3LP7LRESlOPgbsGIykHr+6mNB0UC/GUDTwTxkRNVMThxNTQvq1jhCpeNKD3tpF9LkFLRmsI/azloX0aRXQU6oJYd/1apVKq1I0nJ27NiBGjVqmLwf6aUwRdHKTHKsSlZqksdk3IKQE2vx+++/IyYmpth2xqvz5enfvz9Onz6NZcuWYfXq1SrtR84J5Yq+JdpXGfKZatWqpY53SeYc78rioGdyCIlp2RbdjsilSbCwcGTxYEGkXjA8Ls8Tkd2SIEDG7omS4YDxvjxv7R53qRTUu3dvvPXWW/j7779VOtG6desKn9+2bVux7bdu3arGE5S8em5pTZs2VYGBpBY1atSo2CJpQKaQsQ6jRo1S4yBmz56Njz/+2OLtNOf4tG3bFvHx8eqYl/xM1VHilT0M5BAiA30suh2RS6chSc9CeQl+K54DbhjA9CQiOyaFPmTsXslCIDWrqRDI0qVLceLECdx8880ICQlRV+PlCrrk4hvJCfvEiRPx2GOPqd4IGfQrKTfWJr0fzzzzDCZMmKDadNNNNyElJQV//fWXqi4kgUB5pLekXbt2aNasmRoPsXTpUlXtyNLMOT4SmEnKkqRHSYB23XXX4fz586oXRVLBTE3tqiwGDOQQJAdTul/LSksydr+akqtJ5NJkzELJnoVi9EBqnGG7+t2qsWFEZC4JCmTsnqTjSg+7XDSTfwerYyyfpMHIJG6ShpSdna2ujH/33XfqJNtIyqVKjn3Hjh3VVfNx48ZZfXCukZRrlV4CqZYkgY20V67Sv/DCCxW+1svLSw1alh4TqUrUrVs3fP/99xZvoznHR1KaJCj73//+hzFjxuDixYtqPggJ2GR8hLVp9OYU33VBxjEMEplKVEq2rZL0+Ne7r3nc+GeRVZKITLB/EfDTQxVvN+wzoMXVmupEZDlygn3y5Ek1R4DMYeCMZJ6B1q1bq3Qest3xKe9nzZxzXI5hIIdxU+MIeLlf+yMrPQsMFohMJNWQLLkdERE5PaYkkcNY9vcF5Gp1qB/uh6lDWyAxLadau1+JnIKUTvWpAWQnl7GBxlAtSbYjInJCjz/+uBrMXBqZ3E0me6PiGDCQw1i486y6Hd6+Djo3tH5FACKntPvL8oMF0W86BzwTUZWUVv7TXrz22mtqUHRpqiv9/E87Pj6lYcBADuHExXTsPH0F0pFwR9viNZWJyERb5gIrCwb8NeoDJP5TyjwM0zkPAxE5tcjISLWQ6RgwkEP4cdc5dXvL9ZGICnLOAWJEViO1LTbMBP54w3C/6zig96uAXseZnomIqEIMGMju5Wt1+Hm3IWC4q11tWzeHyPGChbWvApveNdzv8T/g5melRh+gcWfpVCIiqhADBrJ7G48mISE1B6H+XujVhJVbiEym0wErnwe2FQzgu/UNoMt/eQCJiMgsDBjI7v24yzDY+fbW0fDyYCVgIpNndF4yDtizwHB/wDtAh4d58IiIyGwucfYlU3rLVOUyC+Gnn35q6+aQGS5n5GL1wQS1fle7WB47IlNo84DFjxmCBY0bMOT/GCwQEVGlOX3AkJ+fj4kTJ2LdunXYs2cP3n77bVy6dMnWzSIT/bo3DnlaPZrHBKFpNGfaJqpQfg7w42hg/4+Amwdw53yg9b08cETO3Jt4cqNhFne5lftWnqF4/PjxVtv/qVOnoNFosHfvXtjSX3/9hRYtWsDT0xNDhgyBq3P6lKTt27ejWbNmiIkxlOLs378/Vq1ahXvuucfWTSMTLNxpGOw8vD17F4gqlJsJLHwAOLYGcPcGhn8FXN+PB47IWR38DVgxuZTyyDOsVh75559/VifR1hIbG4sLFy4gPNy28y3JxebWrVtj+fLlCAgIqPL+XnnlFfzyyy82D4Sctodhw4YNGDRoEKKjo1XEKQe7pLlz56JevXrw8fFBp06dVJBgdP78+cJgQch6XFxctbWfKu9AXAr+vZAKL3c3DG4VzUNJVJ6cNOCbuwzBgqcfcO8PDBaInD1YWDiyeLAgUi8YHpfnrSA0NBSBgYFW2Xdubi7c3d1Rs2ZNeHjY9pr28ePH0bNnT9SuXRs1atSAvZBjZAt2HzBkZGSgVatWKigozQ8//KCiwJdffhm7d+9W2/bt2xeJiYnV3layrB8LZna+tVkUavh58fASlSXrCrBgKHB6E+AVCNz/M9CwB48XkaOVQM7NMG3JTgWWT5IXlbYjw430PMh2puxP3rsSKUkffvihGh8qF2yjoqJw5513FtvuySefVEtwcLDqMXjppZegL/JecrH39ddfx8iRI9UMy48++ug1KUkyI7LcX7lyJdq0aQNfX191Ii/neXL1v0mTJuq19957LzIzMwv3rdPpMG3aNNSvX1+9Rs4PFy1aVOHnM76/pK8/+OCDav2LL76AVqvFQw89VLg/GRs7Z86cYq+Vtnbs2BH+/v4qyOjatStOnz6tXv/qq69i3759an/GfYrk5GQ8/PDDiIiIUJ9DPptsV7RnQno6ZAyuvLcca1uw+5QkSSGSpSyzZs3CI488gjFjxqj78+bNw++//4758+fjueeeUz0TRXsUZF2+zLLk5OSoxSg1NVXd5uXlFY6JEBL5ymPypZdcN0bIspRcl+fd3NzUe0iXXmnrXl6Gk2PZvui6t7e3+gWQ9yptXdom28sPtSwl1+V5+UWV9yq5bm+fSe/mgV/2Gq6a3NU+1ik+kzN+T/xMdvA95aUCXw+FJn4/4FMDuvt+Ql5kC3gX/IPJ74m/T/wLine truncated
|
||||
"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
|
||||
"text/plain": [
|
||||
"<Figure size 900x600 with 1 Axes>"
|
||||
]
|
||||
@@ -1091,7 +650,7 @@
|
||||
"- `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`\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"
|
||||
]
|
||||
@@ -1106,7 +665,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 31,
|
||||
"execution_count": 19,
|
||||
"id": "4f377075-e352-46f9-a716-f1de673a7086",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
@@ -1123,11 +682,11 @@
|
||||
"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. 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",
|
||||
"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 nearly 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",
|
||||
"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",
|
||||
@@ -1140,13 +699,13 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 32,
|
||||
"execution_count": 20,
|
||||
"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",
|
||||
" \"\"\"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",
|
||||
@@ -1163,7 +722,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 33,
|
||||
"execution_count": 21,
|
||||
"id": "01838139-52dd-48a0-8f62-bcdef351345e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
@@ -1171,8 +730,8 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 27 μs, sys: 0 ns, total: 27 μs\n",
|
||||
"Wall time: 28.8 μs\n"
|
||||
"CPU times: user 21 μs, sys: 1e+03 ns, total: 22 μs\n",
|
||||
"Wall time: 22.2 μs\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -1181,7 +740,7 @@
|
||||
"1"
|
||||
]
|
||||
},
|
||||
"execution_count": 33,
|
||||
"execution_count": 21,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -1218,7 +777,7 @@
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.13.3"
|
||||
"version": "3.13.9"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
||||
Reference in new issue
Block a user