Add files via upload
This commit is contained in:
Peter Norvig authored and 
GitHub committed
1 parent
9cff05c81b
commit
fe022118e6
1 file changed
+689
@@ -0,0 +1,689 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8812d495-6014-4669-97ae-003be4f54605",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"<div style=\"text-align:right\"><i>Peter Norvig<br>Jan 2026</i></div>\n",
|
||||
"\n",
|
||||
"# Truncatable Primes\n",
|
||||
"\n",
|
||||
"What's special about [this pencil](https://mathsgear.co.uk/products/truncatable-prime-pencil)?\n",
|
||||
"\n",
|
||||
"[](https://community.wolfram.com/groups/-/m/t/1569707)\n",
|
||||
"\n",
|
||||
"The number printed on it is a prime, and as you sharpen the pencil and remove digits one at a time from the left, the resulting numbers are all primes:\n",
|
||||
"\n",
|
||||
" 357686312646216567629137 is prime\n",
|
||||
" 57686312646216567629137 is prime\n",
|
||||
" 7686312646216567629137 is prime\n",
|
||||
" ...\n",
|
||||
" 137 is prime\n",
|
||||
" 37 is prime\n",
|
||||
" 7 is prime\n",
|
||||
"\n",
|
||||
"Numbers like this are called [**truncatable primes**](https://en.wikipedia.org/wiki/Truncatable_prime). I thought I would write a program to find other truncatable primes.\n",
|
||||
"\n",
|
||||
"My function `left_truncatable_primes` below starts with the list of one-digit primes: [2, 3, 5, 7]. Then it tries placing each possible digit 1–9 to the left of each of the one-digit primes, giving a list of two-digit candidates. From the candidates it filters out just the primes and recursively tries to add digits to them, giving us three-digit primes, then four-digit primes, and so on, stopping when none of the candidates are primes. It gathers up all the truncatable primes and returns them in sorted order."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "8ab95f58-0d60-4516-b19c-5f55614abb54",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from sympy import isprime # isprime checks if a number is a prime\n",
|
||||
"\n",
|
||||
"def left_truncatable_primes(starting_primes=[2, 3, 5, 7]) -> list[int]:\n",
|
||||
" \"\"\"All left-truncatable primes, in ascending order.\"\"\"\n",
|
||||
" if not starting_primes:\n",
|
||||
" return []\n",
|
||||
" else:\n",
|
||||
" candidates = [int(d + str(p)) for d in \"123456789\" for p in starting_primes]\n",
|
||||
" return starting_primes + left_truncatable_primes(list(filter(isprime, candidates)))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "90d8f036-da26-4975-a0da-ab637694a679",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's see how many truncatable primes there are:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "4c9a569c-af7f-4335-b10b-785e5709862b",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"4260"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"P = left_truncatable_primes()\n",
|
||||
"\n",
|
||||
"len(P)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a32b56c4-587b-48d5-91ec-1abcef2e909f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are 4260 left-truncatable primes. Here are the smallest and largest ones:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "e937655d-16eb-44ce-8698-0a0c09c9f6a9",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[2, 3, 5, 7, 13, 17, 23, 37, 43, 47, 53, 67, 73, 83, 97, 113]"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"P[:16]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "502be81c-fb57-431a-855a-21184a1aeb0d",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[66276812967623946997,\n",
|
||||
" 67986315421273233617,\n",
|
||||
" 86312646216567629137,\n",
|
||||
" 315396334245663786197,\n",
|
||||
" 367986315421273233617,\n",
|
||||
" 666276812967623946997,\n",
|
||||
" 686312646216567629137,\n",
|
||||
" 918918997653319693967,\n",
|
||||
" 5918918997653319693967,\n",
|
||||
" 6686312646216567629137,\n",
|
||||
" 7686312646216567629137,\n",
|
||||
" 9918918997653319693967,\n",
|
||||
" 57686312646216567629137,\n",
|
||||
" 95918918997653319693967,\n",
|
||||
" 96686312646216567629137,\n",
|
||||
" 357686312646216567629137]"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"P[-16:]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "853d5a6c-936e-4ab4-afc3-2082641063fc",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Below is the count for each digit-size, with a bar chart. We see there is only one 24-digit left-truncatable prime, the one on the pencil, and that the most common number of digits is 9.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"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": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHaxJREFUeJzt3X2s1nX9P/AXh1sBAaEASW4sKzQRFqaw1EwJQuY0+KMbZ9SYLUKWsEhoiIoVjL4L06G2ZmBL1NiyJpiJmLgCrGguw2LadFAIdDPAm7g/v73f+13Xl4P4NQ4Hr/e5zuOxffY5n5vrut7X5/qc6zzP++bzadfY2NgYAAAFaah1AQAAjiagAADFEVAAgOIIKABAcQQUAKA4AgoAUBwBBQAojoACABSnQ7RChw8fjm3btsWpp54a7dq1q3VxAID/Qro27KuvvhoDBgyIhoaG+gsoKZwMHDiw1sUAAJph69atccYZZ9RfQEk1J5U32KNHj1oXBwD4L+zZsydXMFT+jtddQKk066RwIqAAQOvy33TP0EkWACiOgAIAFEdAAQCKI6AAAMURUACA4ggoAEBxBBQAoDgCCgBQHAEFACiOgAIAFEdAAQCKI6AAAMURUACA4ggoAEBxOtS6ANCWDJm9qtmPfXnhhBYtC0DJ1KAAAMURUACA4ggoAEBxBBQAoDgCCgBQHAEFACiOgAIAFEdAAQCKI6AAAMURUACA4ggoAEBxBBQAoDgCCgBQHAEFACiOgAIAFKdDrQsApRsye1WzH/vywgktWhaAtkINCgBQHDUo0Eqp2QHqmRoUAKA4AgoAUBwBBQAojoACABRHQAEAiiOgAACtO6Dccsst0a5duybT0KFDq9v37t0b06ZNiz59+kT37t1j0qRJsWPHjibPsWXLlpgwYUJ07do1+vbtG7NmzYqDBw+23DsCANredVA+9KEPxRNPPPG/T9Dhf59ixowZsWrVqlixYkX07Nkzrr/++pg4cWL85je/ydsPHTqUw0n//v1j3bp18corr8TnP//56NixY3z7299uqfcEALS1gJICSQoYR9u9e3fce++9sXz58rjsssvyuqVLl8bZZ58dGzZsiFGjRsXjjz8ezz//fA44/fr1ixEjRsRtt90WN954Y66d6dSpU8u8KwCgbfVBeeGFF2LAgAHx3ve+N6655prcZJNs3LgxDhw4EGPGjKnum5p/Bg0aFOvXr8/LaT5s2LAcTirGjRsXe/bsiU2bNr3la+7bty/vc+QEANSv4wooF154YSxbtiwee+yxuPvuu+Oll16Kiy++OF599dXYvn17rgHp1atXk8ekMJK2JWl+ZDipbK9seysLFizITUaVaeDAgcdTbACgnpt4xo8fX/35vPPOy4Fl8ODB8ZOf/CROOeWUOFnmzJkTM2fOrC6nGhQhBQDq1wkNM061JR/4wAfixRdfzP1S9u/fH7t27WqyTxrFU+mzkuZHj+qpLB+rX0tF586do0ePHk0mAKB+nVBAee211+Kvf/1rnH766TFy5Mg8GmfNmjXV7Zs3b859VEaPHp2X0/y5556LnTt3VvdZvXp1DhznnHPOiRQFAGirTTxf+9rX4sorr8zNOtu2bYubb7452rdvH5/97Gdz35ApU6bkppjevXvn0DF9+vQcStIInmTs2LE5iFx77bWxaNGi3O9k7ty5+dopqZYEAOC4A8rf/va3HEb+9a9/xbvf/e646KKL8hDi9HOyePHiaGhoyBdoSyNv0gidu+66q/r4FGZWrlwZU6dOzcGlW7duMXny5Jg/f75PAwBoXkB58MEH/8/tXbp0iSVLluTpraTal0cfffR4XhYAaGPciwcAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIrTodYFAGpvyOxVzX7sywsntGhZABI1KABAcQQUAKA4AgoAUBwBBQAojoACABRHQAEAiiOgAADFEVAAgOIIKABAcQQUAKA4AgoAUBz34qEuubcMQOumBgUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQH0FlIULF0a7du3ihhtuqK7bu3dvTJs2Lfr06RPdu3ePSZMmxY4dO5o8bsuWLTFhwoTo2rVr9O3bN2bNmhUHDx48kaIAAHWk2QHld7/7XXz/+9+P8847r8n6GTNmxCOPPBIrVqyItWvXxrZt22LixInV7YcOHcrhZP/+/bFu3bq47777YtmyZTFv3rwTeycAQNsOKK+99lpcc8018YMf/CBOO+206vrdu3fHvffeG9/97nfjsssui5EjR8bSpUtzENmwYUPe5/HHH4/nn38+fvzjH8eIESNi/Pjxcdttt8WSJUtyaAEAaFZASU04qRZkzJgxTdZv3LgxDhw40GT90KFDY9CgQbF+/fq8nObDhg2Lfv36VfcZN25c7NmzJzZt2nTM19u3b1/efuQEANSvDsf7gAcffDD+8Ic/5Caeo23fvj06deoUvXr1arI+hZG0rbLPkeGksr2y7VgWLFgQt9566/EWFQBoCzUoW7duja9+9atx//33R5cuXeKdMmfOnNx8VJlSOQCA+nVcASU14ezcuTM+/OEPR4cOHfKUOsLecccd+edUE5L6kezatavJ49Ionv79++ef0/zoUT2V5co+R+vcuXP06NGjyQQA1K/jCiiXX355PPfcc/Hss89Wp/PPPz93mK383LFjx1izZk31MZs3b87DikePHp2X0zw9Rwo6FatXr86h45xzzmnJ9wYAtIU+KKeeemqce+65TdZ169YtX/Oksn7KlCkxc+bM6N27dw4d06dPz6Fk1KhRefvYsWNzELn22mtj0aJFud/J3Llzc8fbVFMCAHDcnWTfzuLFi6OhoSFfoC2NvkkjdO66667q9vbt28fKlStj6tSpObikgDN58uSYP3++TwMAaJmA8tRTTzVZTp1n0zVN0vRWBg8eHI8++uiJvjQAUKfciwcAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOB1qXQCgvgyZvarZj3154YQWLQvQeqlBAQCKI6AAAMURUACA4ggoAEBxBBQAoDgCCgBQHAEFACiOgAIAFEdAAQCKI6AAAMURUACA4ggoAEBxBBQAoDgCCgBQHAEFACiOgAIAFEdAAQCK06HWBYCKIbNXndDBeHnhBAcToE6oQQEAiiOgAADFEVAAgOIIKABAcQQUAKA4AgoAUBwBBQAojoACABRHQAEAiiOgAADFEVAAgNYdUO6+++4477zzokePHnkaPXp0/OIXv6hu37t3b0ybNi369OkT3bt3j0mTJsWOHTuaPMeWLVtiwoQJ0bVr1+jbt2/MmjUrDh482HLvCABoWwHljDPOiIULF8bGjRvj97//fVx22WVx1VVXxaZNm/L2GTNmxCOPPBIrVqyItWvXxrZt22LixInVxx86dCiHk/3798e6devivvvui2XLlsW8efNa/p0BAG3jbsZXXnllk+VvfetbuVZlw4YNObzce++9sXz58hxckqVLl8bZZ5+dt48aNSoef/zxeP755+OJJ56Ifv36xYgRI+K2226LG2+8MW655Zbo1KlTy747AKBt9UFJtSEPPvhgvP7667mpJ9WqHDhwIMaMGVPdZ+jQoTFo0KBYv359Xk7zYcOG5XBSMW7cuNizZ0+1FgYA4LhqUJLnnnsuB5LU3yT1M3n44YfjnHPOiWeffTbXgPTq1avJ/imMbN++Pf+c5keGk8r2yra3sm/fvjxVpEADANSv465B+eAHP5jDyDPPPBNTp06NyZMn52abk2nBggXRs2fP6jRw4MCT+noAQCsLKKmW5KyzzoqRI0fm4DB8+PD43ve+F/3798+dX3ft2tVk/zSKJ21L0vzoUT2V5co+xzJnzpzYvXt3ddq6devxFhsAaEvXQTl8+HBufkmBpWPHjrFmzZrqts2bN+dhxalJKEnz1ES0c+fO6j6rV6/OQ5ZTM9Fb6dy5c3Voc2UCAOrXcfVBSTUZ48ePzx1fX3311Txi56mnnopf/vKXuellypQpMXPmzOjdu3cOEdOnT8+hJI3gScaOHZuDyLXXXhuLFi3K/U7mzp2br52SQggAwHEHlFTz8fnPfz5eeeWVHEjSRdtSOPnEJz6Rty9evDgaGhryBdpSrUoaoXPXXXdVH9++fftYuXJl7ruSgku3bt1yH5b58+f7NACA5gWUdJ2T/0uXLl1iyZIleXorgwcPjkcfffR4XhYAaGPciwcAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAaN1XkgV4Jw2ZvarZj3154YQWLQvwzlKDAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKE6HWhcA4J0wZPaqZj/25YUTWrQswNtTgwIAFEdAAQCKI6AAAK07oCxYsCA+8pGPxKmnnhp9+/aNq6++OjZv3txkn71798a0adOiT58+0b1795g0aVLs2LGjyT5Line truncated
|
||||
"text/plain": [
|
||||
"<Figure size 640x480 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from collections import Counter\n",
|
||||
"\n",
|
||||
"digits = Counter(len(str(p)) for p in P)\n",
|
||||
"\n",
|
||||
"plt.bar(list(digits), list(digits.values()))\n",
|
||||
"dict(digits)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "e82700e8-8fb2-45cf-bcf0-7805e46174bb",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Right-Truncatable Primes\n",
|
||||
"\n",
|
||||
"What if you sharpen the pencil from the other end? For our 357686312646216567629137 pencil it wouldn't work; removing the \"7\" from the right results in a composite number. But it is possible to build up the **right-truncatable** primes:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"id": "3870bc47-5833-4879-b622-2094b4247239",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def right_truncatable_primes(starting_primes=[2, 3, 5, 7]) -> list[int]:\n",
|
||||
" \"\"\"All right-truncatable primes, in ascending order.\"\"\"\n",
|
||||
" if not starting_primes:\n",
|
||||
" return []\n",
|
||||
" else:\n",
|
||||
" candidates = [10 * p + d for p in starting_primes for d in (1, 3, 7, 9)]\n",
|
||||
" return starting_primes + right_truncatable_primes(list(filter(isprime, candidates)))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "976c7f1f-79c5-48c3-8909-62b809993057",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"(Note that I only try placing the digits (1, 3, 7, 9) to the right: placing an even digit or a 5 would always result in a composite number. Also, note that I don't have to do string concatenation to form the new candidates; it is simpler to do `10 * p + d`.)\n",
|
||||
"\n",
|
||||
"Let's find and count the right-truncatable primes:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"id": "7b39ffb5-ba52-42e6-a270-79bc10b26c2c",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"83"
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"Q = right_truncatable_primes()\n",
|
||||
"\n",
|
||||
"len(Q)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5323ae60-285d-49d6-a133-5c476eb4c5e5",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are only 83 right-truncatable primes, so we might as well see them all:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"id": "e3780212-2b86-469f-9df0-81bd6fcc9c94",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"[2, 3, 5, 7, 23, 29, 31, 37, 53, 59, 71, 73, 79, 233, 239, 293, 311, 313, 317, 373, 379, 593, 599, 719, 733, 739, 797, 2333, 2339, 2393, 2399, 2939, 3119, 3137, 3733, 3739, 3793, 3797, 5939, 7193, 7331, 7333, 7393, 23333, 23339, 23399, 23993, 29399, 31193, 31379, 37337, 37339, 37397, 59393, 59399, 71933, 73331, 73939, 233993, 239933, 293999, 373379, 373393, 593933, 593993, 719333, 739391, 739393, 739397, 739399, 2339933, 2399333, 2939999, 3733799, 5939333, 7393913, 7393931, 7393933, 23399339, 29399999, 37337999, 59393339, 73939133]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(Q)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d4abc561-5a2d-4e93-a18e-b78b104f7b7c",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Here is the count of digit sizes:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"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": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAGpBJREFUeJzt3QuMVOX9+OHvKrpQCigoCpWbVkUEUStaLz+FSCUEEWqqYtBSSLS2KiotlW2LSr0sNI3BVgJqU9FW0CYK3iLWWpEapbogXtqGiyJstEiT6q5gXK3MP+f8w8ZVrNLOvjhnnyc5rjN7ds+ZnWH3M+9550xVqVQqBQBAIrul2hAAgPgAAJIz8gEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBIql18wWzbti3eeOON6NSpU1RVVe3q3QEAPofsnKXvvPNO9OzZM3bbbbfKio8sPHr16rWrdwMA+C/U19fHAQccUFnxkY14bN/5zp077+rdAQA+h8bGxnzwYPvf8YqKj+2HWrLwEB8AUFk+z5QJE04BgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcA8MWOj2XLlsXo0aPzd63LTqG6ePHiT6zz97//Pc4444zo0qVLdOzYMYYMGRIbN24s1z4DAG0pPrZu3RqDBw+OOXPm7PDzr7zySpx00knRv3//WLp0abz44osxffr0aN++fTn2FwCocFWlUqn0X39xVVUsWrQoxo4d23zduHHjYo899ojf/va3//W74mUjJg0NDd5YDgAqxM78/S7rnI9t27bFww8/HIccckiMGDEiunfvHscdd9wOD81s19TUlO/wRxcAoLjalfObbd68ObZs2RIzZ86M6667LmbNmhVLliyJM888M5544ok45ZRTPvE1tbW1MWPGjHLuBhRK32kPR6V6beaoXb0LwBdQ2Uc+MmPGjIkrrrgijjzyyJg2bVqcfvrpMW/evB1+TU1NTT5Es32pr68v5y4BAEUe+dhnn32iXbt2MWDAgBbXH3bYYfHUU0/t8Guqq6vzBQBoG8o68rHnnnvmL6tdvXp1i+vXrFkTffr0KeemAIC2MvKRzelYt25d8+X169fHqlWromvXrtG7d++YOnVqnHPOOXHyySfHsGHD8jkfDz74YP6yWwCAnY6Purq6PCq2mzJlSv5xwoQJMX/+/PjmN7+Zz+/IJpJOnjw5Dj300Lj33nvzc38AAOx0fAwdOjQ+69QgkyZNyhcAgI/z3i4AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQVLu0mwPYsb7THq7YH81rM0ft6l2AimLkAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAPhix8eyZcti9OjR0bNnz6iqqorFixd/6roXXXRRvs7s2bP/1/0EANpqfGzdujUGDx4cc+bM+Y/rLVq0KJYvX55HCgDAdu1iJ40cOTJf/pPXX389Lr300nj00Udj1KhRO7sJAKDAdjo+Psu2bdvi/PPPj6lTp8bhhx/+mes3NTXly3aNjY3l3iUAoMjxMWvWrGjXrl1Mnjz5c61fW1sbM2bMKPdu0Eb0nfZwVKLXZhoRBNqusr7aZcWKFXHTTTfF/Pnz84mmn0dNTU00NDQ0L/X19eXcJQCgyPHx5z//OTZv3hy9e/fORz+yZcOGDfGDH/wg+vbtu8Ovqa6ujs6dO7dYAIDiKuthl2yux/Dhw1tcN2LEiPz6iRMnlnNTAEBbiY8tW7bEunXrmi+vX78+Vq1aFV27ds1HPLp169Zi/T322CP233//OPTQQ8uzxwBA24qPurq6GDZsWPPlKVOm5B8nTJiQz/UAAChrfAwdOjRKpdLnXv+1117b2U0AAAXmvV0AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AIAvdnwsW7YsRo8eHT179oyqqqpYvHhx8+c++OCDuPLKK2PQoEHRsWPHfJ1vf/vb8cYbb5R7vwGAthIfW7dujcGDB8ecOXM+8bl33303Vq5cGdOnT88/3nfffbF69eo444wzyrW/AECFa7ezXzBy5Mh82ZEuXbrEY4891uK6m2++OY499tjYuHFj9O7d+7/fUwCgbcbHzmpoaMgPz+y11147/HxTU1O+bNfY2NjauwQAFDU+3nvvvXwOyLnnnhudO3fe4Tq1tbUxY8aM1twNgC+MvtMejkr02sxRu3oXKJBWe7VLNvn07LPPjlKpFHPnzv3U9WpqavLRke1LfX19a+0SAFDUkY/t4bFhw4b405/+9KmjHpnq6up8AQDahnatFR5r166NJ554Irp161buTQAAbSk+tmzZEuvWrWu+vH79+li1alV07do1evToEd/61rfyl9k+9NBD8eGHH8amTZvy9bLP77nnnuXdewCg+PFRV1cXw4YNa748ZcqU/OOECRPimmuuiQceeCC/fOSRR7b4umwUZOjQof/7HgMAbSs+soDIJpF+mv/0OQAA7+0CACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQB8seNj2bJlMXr06OjZs2dUVVXF4sWLW3y+VCrFVVddFT169IgOHTrE8OHDY+3ateXcZwCgLcXH1q1bY/DgwTFnzpwdfv7nP/95/PKXv4x58+bFX/7yl+jYsWOMGDEi3nvvvXLsLwBQ4drt7BeMHDkyX3YkG/WYPXt2/PSnP40xY8bk1915552x33775SMk48aN+9/3GACoaGWd87F+/frYtGlTfqhluy5dusRxxx0XzzzzzA6/pqmpKRobG1ssAEBxlTU+svDIZCMdH5Vd3v65j6utrc0DZfvSq1evcu4SAPAFs8tf7VJTUxMNDQ3NS319/a7eJQCgUuJj//33zz+++eabLa7PLm//3MdVV1dH586dWywAQHGVNT769euXR8bjjz/efF02hyN71cvxxx9fzk0BAG3l1S5btmyJdevWtZhkumrVqujatWv07t07Lr/88rjuuuvi4IMPzmNk+vTp+TlBxo4dW+59BwDaQnzU1dXFsGHDmi9PmTIl/zhhwoSYP39+/OhHP8rPBXLhhRfG22+/HSeddFIsWbIk2rdvX949BwDaRnwMHTo0P5/Hp8nOevqzn/0sXwAAvnCvdgEA2hbxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBIqp2fdzH1nfZwVKLXZo7a1bsAQCsz8gEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMAEB8AQHEZ+QAAkhIfAEBlx8eHH34Y06dPj379+kWHDh3ioIMOimuvvTZKpVK5NwUAVKB25f6Gs2bNirlz58Ydd9wRhx9+eNTV1cXEiROjS5cuMXny5HJvDgBo6/Hx9NNPx5gxY2LUqFH55b59+8bChQvj2WefLfemAIAKVPbDLieccEI8/vjjsWbNmvzyCy+8EE899VSMHDlyh+s3NTVFY2NjiwUAKK6yj3xMmzYtD4j+/fvH7rvvns8Buf7662P8+PE7XL+2tjZmzJhR7t0AYBfqO+3hivz5vzbz/4/aU2EjH7///e/jrrvuigULFsTKlSvzuR+/+MUv8o87UlNTEw0NDc1LfX19uXcJACjyyMfUqVPz0Y9x48bllwcNGhQbNmzIRzgmTJjwifWrq6vzBQBoG8o+8vHuu+/Gbru1/LbZ4Zdt27aVe1MAQAUq+8jH6NGj8zkevXv3zl9q+/zzz8eNN94YkyZNKvemAIAKVPb4+NWvfpWfZOz73/9+bN68OXr27Bnf/e5346qrrir3pgCAClT2+OjUqVPMnj07XwAAPs57uwAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAJUfH6+//nqcd9550a1bt+jQoUMMGjQo6urqWmNTAECFaVfub/jWW2/FiSeeGMOGDYtHHnkk9t1331i7dm3svffe5d4UAFCByh4fs2bNil69esXtt9/efF2/fv3KvRkAoEKV/bDLAw88EMccc0ycddZZ0b179zjqqKPitttu+9T1m5qaorGxscUCABRX2ePj1Vdfjblz58bBBx8cjz76aHzve9+LyZMnxx133LHD9Wtra6NLly7NSzZqAgAUV9njY9u2bXH00UfHDTfckI96XHjhhXHBBRfEvHnzdrh+TU1NNDQ0NC/19fXl3iUAoMjx0aNHjxgwYECL6w477LDYuHHjDtevrq6Ozp07t1gAgOIqe3xkr3RZvXp1i+vWrFkTffr0KfemAIAKVPb4uOKKK2L58uX5YZd169bFggUL4tZbb42LL7643JsCACpQ2eNLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 640x480 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"digits = dict(Counter(len(str(p)) for p in Q))\n",
|
||||
"\n",
|
||||
"plt.bar(list(digits), list(digits.values()));\n",
|
||||
"digits"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6b307b24-a93e-4353-9bee-fba739072bd1",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Summary of Truncatable Primes\n",
|
||||
"\n",
|
||||
"Here's what we learned:\n",
|
||||
"\n",
|
||||
"- Left-trunctable primes:\n",
|
||||
" - There are 4260 of them\n",
|
||||
" - The largest has 24 digits: 357686312646216567629137\n",
|
||||
" - The plurality have 9 digits; few have more than 20 digits\n",
|
||||
"- Right-truncatable primes:\n",
|
||||
" - There are only 83 of them\n",
|
||||
" - The largest has 8 digits: 73939133\n",
|
||||
" - The plurality have 4 digits\n",
|
||||
"\n",
|
||||
"# Note on Primality Checking\n",
|
||||
"\n",
|
||||
"A prime number can be defined as *\"an integer greater than 1 that cannot be evenly divided by any whole number other than 1 and itself.\"* Here's a naive function to test whether a number is prime:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"id": "1d460931-3d70-485d-807f-8fb360ef4ef2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def isprime_naive(n: int) -> bool:\n",
|
||||
" \"\"\"Simple primality checker.\"\"\"\n",
|
||||
" divisors = range(2, n)\n",
|
||||
" return n > 1 and not any(n % d == 0 for d in divisors)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "f30802c2-2e39-4cdd-8e33-fc74ad50977f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"And here's a function to test it. By default we test some easy cases, along with one 8-digit prime and a couple of 8-digit composites."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "4f377075-e352-46f9-a716-f1de673a7086",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"tests = {True: [2, 3, 5, 7, 37, 73, 101, 11939, 117223, 73939133],\n",
|
||||
" False: [0, 1, 4, 6, 8, 9, 10, 256, 123456, 7333*7393, 7393**2]}\n",
|
||||
"\n",
|
||||
"def prime_test(predicate, tests=tests) -> None:\n",
|
||||
" \"\"\"Test that a primality checking predicate gets the right answers on all the test cases.\n",
|
||||
" `tests` is a dict of {True: [primes,...], False: [composites,...]}.\"\"\"\n",
|
||||
" for result in tests:\n",
|
||||
" for n in tests[result]:\n",
|
||||
" assert predicate(n) == result, f'isprime({n}) should be {result}'"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"id": "328285d6-9d5c-44ef-8f84-7dbf2b6fe6bc",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 2.1 s, sys: 18.7 ms, total: 2.12 s\n",
|
||||
"Wall time: 2.12 s\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%time prime_test(isprime_naive)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1834f8f0-3712-4302-ac56-79972f85d0db",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We see that `isprime_naive` passes the tests, but it takes 2 seconds, and we're only up to 8-digit numbers. That's too slow!\n",
|
||||
"\n",
|
||||
"We can make it faster with two ideas. First, we don't have to check divisors all the way up to *n*. If *p* is composite, then *p* is the product of two numbers, and one of them must be less than or equal to the square root of *n*, so that's as far as we need to go. Second, once we check to see if *n* is divisible by 2, we don't have to check whether it is divisible by any other even number. "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"id": "647fb66b-32ae-494a-8945-cf7358c2f958",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from math import sqrt\n",
|
||||
"\n",
|
||||
"def isprime_faster(n: int) -> bool:\n",
|
||||
" \"\"\"More sophisticated primality checker: go to square root, checking odd numbers only.\"\"\"\n",
|
||||
" if n <= 10:\n",
|
||||
" return n in (2, 3, 5, 7)\n",
|
||||
" elif n % 2 == 0:\n",
|
||||
" return False\n",
|
||||
" else:\n",
|
||||
" divisors = range(3, int(sqrt(n)) + 1, 2) # Check odd divisors up to sqrt(n)\n",
|
||||
" return not any(n % d == 0 for d in divisors)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"id": "3bad2a74-49a2-43ce-aaf5-96112ae3c8d9",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 405 μs, sys: 1 μs, total: 406 μs\n",
|
||||
"Wall time: 407 μs\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%time prime_test(isprime_faster)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "462fbfea-19ca-487e-b104-96eea821c001",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We see that this function is roughly a 5,000 times faster than the naive function on our specific test cases. But it would be wrong to conclude that it is always 5,000 times faster; more importantly it is *O*(√*n*) whereas the naive version is *O*(*n*). \n",
|
||||
"\n",
|
||||
"Let's try some tougher test cases, a pair of 16-digit numbers:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"id": "b8664d7a-2033-47a1-b94a-d1e1ad77e31f",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 2.51 s, sys: 21.2 ms, total: 2.53 s\n",
|
||||
"Wall time: 2.53 s\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"test16 = {True: [8637267627626947], False: [73939133**2]}\n",
|
||||
"\n",
|
||||
"%time prime_test(isprime_faster, test16)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "251e2b35-91b1-40cf-a57e-a9b86947feef",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"As expected, `isprime_faster` can handle 16-digit numbers in about the same time `isprime_naive` handles 8-digit numbers. \n",
|
||||
"\n",
|
||||
"What about `sympy.isprime`?"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"id": "e71701b0-6ad3-4a5d-b871-c433712ee99c",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 56 μs, sys: 0 ns, total: 56 μs\n",
|
||||
"Wall time: 57.2 μs\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%time prime_test(isprime, test16)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "05b8cbf7-ea10-4047-bcdb-a5ac7981384b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"That's about 40,000 times faster than `isprime_faster`! What makes `sympy.isprime` so fast? The answer is: [Number theory](https://en.wikipedia.org/wiki/Number_theory).\n",
|
||||
"\n",
|
||||
"The big breakthrough came in 1640, when Pierre de Fermat [showed](https://en.wikipedia.org/wiki/Fermat_primality_test) that if *n* is prime and *a* is not divisible by *n*, then *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*). \n",
|
||||
"\n",
|
||||
"That means we can check if *n* is prime by choosing a random *a* and testing if *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*). If the test is false then *n* is definitely composite. If the test is true, then we're not sure: *n* might be prime, might be composite. But if we repeat the test multiple times with different values of *a*, and they are all true, then that is stronger evidence that *n* is prime. This is called a [Monte Carlo algorithm](https://en.wikipedia.org/wiki/Monte_Carlo_algorithm); an algorithm that uses randomization, and can sometimes be wrong.\n",
|
||||
"\n",
|
||||
"Here is an implementation. (The optional third argument to the built-in `pow` function is the modulus.) "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"id": "476282dd-a48a-4eb5-8fed-8495d4878ec2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import random\n",
|
||||
"\n",
|
||||
"def isprime_fermat(n: int, k=20) -> bool:\n",
|
||||
" \"\"\"n is probably a prime if this returns True; definitely composite if it returns False.\"\"\"\n",
|
||||
" if n <= 10:\n",
|
||||
" return n in (2, 3, 5, 7)\n",
|
||||
" else:\n",
|
||||
" a_values = (random.randint(2, n - 1) for _ in range(k))\n",
|
||||
" return any(pow(a, n - 1, n) == 1 for a in a_values)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1ec7bba6-ea89-4255-a73c-bc25a32da883",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's test it out:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"id": "ef715472-cacf-4af0-a8e9-c54039d9013c",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 99 μs, sys: 2 μs, total: 101 μs\n",
|
||||
"Wall time: 102 μs\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"prime_test(isprime_fermat)\n",
|
||||
"\n",
|
||||
"%time prime_test(isprime_fermat, test16)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "b0b64505-706b-4b09-8f41-c8bfcda9b0c9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"That's pretty good! Only a few lines of code, it passes the tests, and it is almost as fast as the highly-tuned (and highly complex) `sympy.isprime`. \n",
|
||||
"\n",
|
||||
"The problem is that `isprime_fermat` can lie:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"id": "a6834bd8-b904-47ed-84d1-dae2c7227cd2",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"True"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"assert 561 == 3 * 11 * 17\n",
|
||||
"\n",
|
||||
"isprime_fermat(561)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "53980370-73e9-49c4-8d14-3d1e38e21e33",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We see that 561 is composite, but `isprime_fermat` claims it is a prime. We can increase the number of samples; that doesn't help:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"id": "cf7a60d0-0b03-4b89-bd53-fdd182743e96",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"True"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"isprime_fermat(561, k=500)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "83f69e44-4f3b-43b8-985c-0db82e983156",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"561 is a [Carmichael number](https://en.wikipedia.org/wiki/Carmichael_number), a number that is resistant to the Fermat test.\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.\n",
|
||||
"\n",
|
||||
"For numbers greater than 2<sup>64</sup> (20 digits), `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",
|
||||
"\n",
|
||||
"I hope you have enjoyed this excursion into the land of prime numbers!"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3 (ipykernel)",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.13.3"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
Reference in new issue
Block a user