Add files via upload
This commit is contained in:
Peter Norvig authored and 
GitHub committed
1 parent
a69a980e7c
commit
23eefc32ff
1 file changed
+358
-194
+358
-194
@@ -13,7 +13,7 @@
|
||||
"\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",
|
||||
"The 24-digit 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",
|
||||
@@ -25,7 +25,7 @@
|
||||
"\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."
|
||||
"My function `left_truncatable_primes` below starts with the list of one-digit primes: [2, 3, 5, 7]. Then it places each possible digit 1–9 to the left of each of the one-digit primes, giving a list of two-digit candidates. From the candidates it filters out just the primes. If there are any new primes, it recursively adds digits to them, giving us three-digit primes, then four-digit primes, and so on, stopping when there are no more new primes. In the end, the function gathers up all the truncatable primes and returns them in sorted order."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -39,11 +39,9 @@
|
||||
"\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)))"
|
||||
" candidates = [int(d + str(p)) for d in \"123456789\" for p in starting_primes]\n",
|
||||
" new_primes = list(filter(isprime, candidates))\n",
|
||||
" return starting_primes + (left_truncatable_primes(new_primes) if new_primes else [])"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -82,7 +80,7 @@
|
||||
"id": "a32b56c4-587b-48d5-91ec-1abcef2e909f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There are 4260 left-truncatable primes. Here are the smallest and largest ones:"
|
||||
"There are 4260 left-truncatable primes. Here are the smallest and largest of them:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -147,7 +145,7 @@
|
||||
"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"
|
||||
"Below is the count for each digit-length, 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"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -191,7 +189,7 @@
|
||||
},
|
||||
{
|
||||
"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
|
||||
"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>"
|
||||
]
|
||||
@@ -204,10 +202,14 @@
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from collections import Counter\n",
|
||||
"\n",
|
||||
"digits = Counter(len(str(p)) for p in P)\n",
|
||||
"def digit_lengths(primes) -> Counter:\n",
|
||||
" \"\"\"Plot a bar chart and return a Counter 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'); plt.ylabel('Count')\n",
|
||||
" return digits\n",
|
||||
"\n",
|
||||
"plt.bar(list(digits), list(digits.values()))\n",
|
||||
"dict(digits)"
|
||||
"dict(digit_lengths(P))"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -229,11 +231,9 @@
|
||||
"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)))"
|
||||
" candidates = [10 * p + d for p in starting_primes for d in (1, 3, 7, 9)]\n",
|
||||
" new_primes = list(filter(isprime, candidates))\n",
|
||||
" return starting_primes + (right_truncatable_primes(new_primes) if new_primes else [])"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -300,7 +300,7 @@
|
||||
"id": "d4abc561-5a2d-4e93-a18e-b78b104f7b7c",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Here is the count of digit sizes:"
|
||||
"Here is the count of digit lengths:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -312,7 +312,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{1: 4, 2: 9, 3: 14, 4: 16, 5: 15, 6: 12, 7: 8, 8: 5}"
|
||||
"Counter({4: 16, 5: 15, 3: 14, 6: 12, 2: 9, 7: 8, 8: 5, 1: 4})"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
@@ -321,7 +321,7 @@
|
||||
},
|
||||
{
|
||||
"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
|
||||
"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>"
|
||||
]
|
||||
@@ -331,15 +331,12 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"digits = dict(Counter(len(str(p)) for p in Q))\n",
|
||||
"\n",
|
||||
"plt.bar(list(digits), list(digits.values()));\n",
|
||||
"digits"
|
||||
"digit_lengths(Q)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6b307b24-a93e-4353-9bee-fba739072bd1",
|
||||
"id": "f2323985-2518-41d1-8a3b-6e21e40d27be",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Summary of Truncatable Primes\n",
|
||||
@@ -357,65 +354,68 @@
|
||||
"\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:"
|
||||
"I was very impressed by the speed of `sympy.isprime`, so I wanted to look more deeply into the topic of primality checking. I'm going to write some prime-checking functions, so the first thing I will do is define a function to test that a prime checker correctly returns True when passed a prime and False when passed a composite number, for some given test cases."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"id": "1d460931-3d70-485d-807f-8fb360ef4ef2",
|
||||
"id": "4f377075-e352-46f9-a716-f1de673a7086",
|
||||
"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)"
|
||||
"from typing import Callable, Iterable\n",
|
||||
"\n",
|
||||
"def test_prime_checker(checker: Callable, \n",
|
||||
" primes=(2, 3, 5, 7, 11, 37, 73, 101, 11939, 65537, 117223, 7629137),\n",
|
||||
" composites=(0, 1, 4, 6, 8, 9, 10, 256, 11939*11939, 11939*117223)) -> Callable:\n",
|
||||
" \"\"\"Test that a primality checking function correctly handles some primes and composites.\"\"\"\n",
|
||||
" for n in primes:\n",
|
||||
" assert checker(n) == True, f'{n} should be prime'\n",
|
||||
" for n in composites:\n",
|
||||
" assert checker(n) == False, f'{n} should be composite'\n",
|
||||
" return checker # This allows us to use @test_prime_checker as a @decorator"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "f30802c2-2e39-4cdd-8e33-fc74ad50977f",
|
||||
"id": "c08ed813-2942-42da-aed3-d5044b0b403a",
|
||||
"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."
|
||||
"Here's how we would use this to check the `isprime` function:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "4f377075-e352-46f9-a716-f1de673a7086",
|
||||
"id": "008cd6fd-86eb-4b45-bea9-c030273c4eb1",
|
||||
"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}'"
|
||||
"assert test_prime_checker(isprime)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6b307b24-a93e-4353-9bee-fba739072bd1",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"But I can also use `@test_prime_checker` as a decorator to test the function right where it is defined; I'll do that for `isprime_simple`, which follows the definition of a prime number almost verbatim:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"id": "328285d6-9d5c-44ef-8f84-7dbf2b6fe6bc",
|
||||
"id": "1d460931-3d70-485d-807f-8fb360ef4ef2",
|
||||
"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"
|
||||
]
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%time prime_test(isprime_naive)"
|
||||
"@test_prime_checker\n",
|
||||
"def isprime_simple(n: int) -> bool:\n",
|
||||
" \"\"\"Simple primality checker. A prime number is defined as an integer greater than 1 \n",
|
||||
" that cannot be evenly divided by any whole number other than 1 and itself.\"\"\"\n",
|
||||
" divisors = range(2, n)\n",
|
||||
" return n > 1 and not any(n % d == 0 for d in divisors)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -423,9 +423,9 @@
|
||||
"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",
|
||||
"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 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. "
|
||||
"We can speed things up a bit with two ideas, one small and one big. The small idea is to cut the run time in half by only checking the odd divisors. Once we determine that 2 is a prime and all other even numbers are composite, we only need to check odd divisors. The big idea is that we only need to check divisors up to √*n*, not up to *n*. If *n* is composite, then *n* is the product of two numbers, and one of them must be less than or equal to √*n*, so we can stop there."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -437,126 +437,64 @@
|
||||
"source": [
|
||||
"from math import sqrt\n",
|
||||
"\n",
|
||||
"@test_prime_checker\n",
|
||||
"def isprime_faster(n: int) -> bool:\n",
|
||||
" \"\"\"More sophisticated primality checker: go to square root, checking odd numbers only.\"\"\"\n",
|
||||
" if n <= 10:\n",
|
||||
" \"\"\"More sophisticated primality checker: go to square root, checking odd divisors only.\"\"\"\n",
|
||||
" if n <= 10: # Handle small numbers up to 10\n",
|
||||
" return n in (2, 3, 5, 7)\n",
|
||||
" elif n % 2 == 0:\n",
|
||||
" elif n % 2 == 0: # Even numbers other than 2 are composite\n",
|
||||
" return False\n",
|
||||
" else:\n",
|
||||
" divisors = range(3, int(sqrt(n)) + 1, 2) # Check odd divisors up to sqrt(n)\n",
|
||||
" else: # Check odd divisors up to sqrt(n)\n",
|
||||
" divisors = range(11, int(sqrt(n)) + 1, 2) \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",
|
||||
"That's a noticeable improvement, but handling 24-digit numbers is still completely infeasible. We need a big breakthrough. \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",
|
||||
"Fortunately, [Pierre de Fermat](https://en.wikipedia.org/wiki/Pierre_de_Fermat) provided that breakthrough in 1640, [showing](https://en.wikipedia.org/wiki/Fermat_primality_test) that if *n* is prime and *a* is not divisible by *n*, then *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*). \n",
|
||||
"\n",
|
||||
"That means we can check if *n* is prime by choosing a random *a* and testing if *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*). If the test is false then *n* is definitely composite. If the test is true, then we're not sure: *n* might be prime, might be composite. But if we 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",
|
||||
"That means we can check if *n* is prime by choosing a random *a* and testing if *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*). If the test is false then *n* is definitely composite. If the test is true, then we're not sure: *n* might be prime, might be composite. But if we choose multiple values of *a* and they all give a remainder of 1 (mod *n*), then that is stronger evidence that *n* is probably prime. This is called a [Monte Carlo algorithm](https://en.wikipedia.org/wiki/Monte_Carlo_algorithm); an algorithm that uses randomization, and can sometimes be wrong.\n",
|
||||
"Some examples:\n",
|
||||
"\n",
|
||||
"Here is an implementation. (The optional third argument to the built-in `pow` function is the modulus.) "
|
||||
"\n",
|
||||
"|*n*|*a*|*a*<sup>(*n* - 1)</sup> (mod *n*)|Conclusion|\n",
|
||||
"|--:|--:|:--:|:--|\n",
|
||||
"|12|5|5|*12 is definitely composite*|\n",
|
||||
"|221|18|1|*221 could be prime or composite*|\n",
|
||||
"|221|2|16|*221 is definitely composite*|\n",
|
||||
"|5|2|1|*5 could be prime or composite*|\n",
|
||||
"|5|3|1|*5 could be prime or composite*|\n",
|
||||
"|5|4|1|*5 could be prime or composite, but we have a lot of evidence that it is prime*|\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Here is an implementation:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"execution_count": 14,
|
||||
"id": "476282dd-a48a-4eb5-8fed-8495d4878ec2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import random\n",
|
||||
"\n",
|
||||
"def sample(lo: int, hi: int, k: int) -> Iterable[int]:\n",
|
||||
" \"\"\"Randomly sample k integers from range(lo, hi), one by one.\"\"\"\n",
|
||||
" return (random.randrange(lo, hi) for _ in range(k))\n",
|
||||
" \n",
|
||||
"@test_prime_checker\n",
|
||||
"def isprime_fermat(n: int, k=20) -> bool:\n",
|
||||
" \"\"\"n is probably a prime if this returns True; definitely composite if it returns False.\"\"\"\n",
|
||||
" if n <= 10:\n",
|
||||
" 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)"
|
||||
" return all(pow(a, n - 1, n) == 1 for a in sample(2, n, k))"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -564,50 +502,180 @@
|
||||
"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)"
|
||||
"Note that the built-in [modular exponentiation](https://en.wikipedia.org/wiki/Modular_exponentiation) function `pow(b, e, m)` computes *b*<sup>*e*</sup> (mod *m*), and does it in an efficient way that doesn't have to multiply *b* by itself *e* times, and doesn't have to deal with any numbers larger than *m*."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "b0b64505-706b-4b09-8f41-c8bfcda9b0c9",
|
||||
"id": "15ef6884-253d-4a3a-b58f-397ac73f552d",
|
||||
"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",
|
||||
"**The problem** is that there are some composite numbers *n* for which *a*<sup>(*n* - 1)</sup> ≡ 1 (mod *n*) for some values of *a*, and thus `isprime_fermat` can lie: it can incorrectly report a **false prime** for a composite number *n*. \n",
|
||||
"\n",
|
||||
"The problem is that `isprime_fermat` can lie:"
|
||||
"How common are these false primes? Here's a function to estimate the percentage of false primes in the range `lo` to `hi` when `isprime_fermat` is allowed `k` choices of the `a` parameter:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"id": "a6834bd8-b904-47ed-84d1-dae2c7227cd2",
|
||||
"execution_count": 15,
|
||||
"id": "45a11a93-9ef3-48d1-80fa-e0113041346e",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def false_prime_percent(lo: int, hi: int, k: int, repeat=100_000) -> float:\n",
|
||||
" \"\"\"The estimated probability of a false prime from isprime_fermat(n, k) for n in range(lo, hi).\"\"\"\n",
|
||||
" composites = (n for n in sample(lo, hi, repeat) if not isprime(n))\n",
|
||||
" counts = Counter(isprime_fermat(n, k) for n in composites) \n",
|
||||
" return 100 * counts[True] / (counts[True] + counts[False])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "778c3a9b-dcc7-4ad2-9092-d66aecd1d93c",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We'll choose ranges that represent digit-lengths, and for now use 1 for the value of `k`:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"id": "9eba96df-c298-4d71-b68a-3bf36bf15095",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"True"
|
||||
"{3: 1.4553237025147137,\n",
|
||||
" 4: 0.5076430372271561,\n",
|
||||
" 5: 0.1345167870334638,\n",
|
||||
" 6: 0.04226359478965734,\n",
|
||||
" 7: 0.00642047704144418,\n",
|
||||
" 8: 0.003178908998431738,\n",
|
||||
" 9: 0.0,\n",
|
||||
" 10: 0.0,\n",
|
||||
" 11: 0.0,\n",
|
||||
" 12: 0.0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 16,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{d: false_prime_percent(10**(d-1), 10**d, k=1) for d in range(3, 13)}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "f3fe25d6-5f0a-478b-90d0-6c89e092e92e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We see that 3-digit composites are falsely called prime about 1.5% of the time, but onnce we get to 10-digit composites, false primes are very rare indeed. Let's see how we reduce false primes by allowing up to 30 choices of `a`:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"id": "8880c588-32f5-46e8-b59c-f9b49d492d3b",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{3: 0.0,\n",
|
||||
" 4: 0.0,\n",
|
||||
" 5: 0.0,\n",
|
||||
" 6: 0.0,\n",
|
||||
" 7: 0.0010682505261133842,\n",
|
||||
" 8: 0.0,\n",
|
||||
" 9: 0.0,\n",
|
||||
" 10: 0.0,\n",
|
||||
" 11: 0.0,\n",
|
||||
" 12: 0.0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{d: false_prime_percent(10**(d-1), 10**d, k=30) for d in range(3, 13)}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "61ee44c3-f963-4fb1-ac20-4a27adc90f83",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"In some runs this produces no false primes; in some runs a few get through. \n",
|
||||
"\n",
|
||||
"There are a few numbers, called [**Carmichael numbers**](https://en.wikipedia.org/wiki/Carmichael_number), that have a particularly high false prime percentage. We can test some of them:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"id": "7871aaa7-af1c-462b-85bb-d46818ccabea",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{561: 58.6,\n",
|
||||
" 1105: 68.8,\n",
|
||||
" 1729: 74.4,\n",
|
||||
" 2465: 72.7,\n",
|
||||
" 2821: 77.9,\n",
|
||||
" 6601: 81.2,\n",
|
||||
" 8911: 80.6,\n",
|
||||
" 41041: 71.3,\n",
|
||||
" 101101: 71.3,\n",
|
||||
" 825265: 61.4,\n",
|
||||
" 321197185: 66.5}"
|
||||
]
|
||||
},
|
||||
"execution_count": 18,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"carmichael_numbers = (561, 1105, 1729, 2465, 2821, 6601, 8911, 41041, 101101, 825265, 321197185)\n",
|
||||
"\n",
|
||||
"{n: false_prime_percent(n, n + 1, k=1, repeat=1000) for n in carmichael_numbers}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4ecd2289-66d5-4b0f-ac00-5db0ad92dde9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"That seems seriously bad; they all are above 50% error rate. But if we set `k=100`, we eliminate most (and on some runs all) of the errors:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"id": "de492047-77b6-46a3-8a77-46308f219efd",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{561: 0.0,\n",
|
||||
" 1105: 0.0,\n",
|
||||
" 1729: 0.0,\n",
|
||||
" 2465: 0.0,\n",
|
||||
" 2821: 0.0,\n",
|
||||
" 6601: 0.0,\n",
|
||||
" 8911: 0.0,\n",
|
||||
" 41041: 0.0,\n",
|
||||
" 101101: 0.0,\n",
|
||||
" 825265: 0.0,\n",
|
||||
" 321197185: 0.0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
@@ -616,52 +684,148 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"assert 561 == 3 * 11 * 17\n",
|
||||
"\n",
|
||||
"isprime_fermat(561)"
|
||||
"{n: false_prime_percent(n, n + 1, k=100, repeat=10_000) for n in carmichael_numbers}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "53980370-73e9-49c4-8d14-3d1e38e21e33",
|
||||
"id": "bfd71dc7-7f6d-4c98-9886-5c242d9b297a",
|
||||
"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:"
|
||||
"So our Fermat test is mostly reliable, but it will on rare occasion let a false prime slip through.\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",
|
||||
"\n",
|
||||
"# Speed of Primality Checking\n",
|
||||
"\n",
|
||||
"**The great thing** about `isprime_fermat` and `sympy.isprime` is that they are **very fast**, handling 32-digit primes in well under a millisecond. Let's make a chart:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"id": "cf7a60d0-0b03-4b89-bd53-fdd182743e96",
|
||||
"id": "8109d9fc-bbae-4ef0-a75c-908d9746a557",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import timeit\n",
|
||||
"\n",
|
||||
"big_primes = ([int('357686312646216567629137'[-i:]) for i in range(5, 25)] +\n",
|
||||
" # The following primes sourced from t5k.org/curios/\n",
|
||||
" [1000000000000000035000061, 59999999999899999999999999, 100000109999990000011000001,\n",
|
||||
" 2728487949505050529272727777, 24444666666888888889999999991, 100003100019100043100057100069,\n",
|
||||
" 9999999999999999777777775555331, 55555555555555555555555555555559])\n",
|
||||
"\n",
|
||||
"def plot_run_times(functions=(isprime_simple, isprime_faster, isprime_fermat, isprime), primes=big_primes):\n",
|
||||
" \"\"\"For each primality-checking function, plot its run time on primes of different digit-lengths.\"\"\"\n",
|
||||
" plt.figure(figsize=(9, 6))\n",
|
||||
" for function in functions:\n",
|
||||
" plt.plot(*time_test(function, primes), 'o-', label=function.__name__)\n",
|
||||
" plt.grid(True, axis='y', which=\"major\", ls=':', color='gray')\n",
|
||||
" plt.xticks(range(len(str(min(primes))), len(str(max(primes))) + 1))\n",
|
||||
" plt.yscale('log')\n",
|
||||
" plt.xlabel('Number of digits in prime')\n",
|
||||
" plt.ylabel('Run time in seconds (log scale)')\n",
|
||||
" plt.legend()\n",
|
||||
"\n",
|
||||
"def time_test(function, primes) -> tuple[list[int], list[float]]:\n",
|
||||
" \"\"\"Time the function on each of the primes, stopping when one exceeds a second.\n",
|
||||
" Return a list of the digit sizes and a list of the corresponding run times.\"\"\"\n",
|
||||
" D, T = [], [] # D is length of primes in digits; T is time in seconds\n",
|
||||
" repeat = 100 # Start with 100 repeats, but reduce when times get longer\n",
|
||||
" for p in primes:\n",
|
||||
" time = timeit.timeit(lambda: function(p), number=repeat) / repeat\n",
|
||||
" D.append(len(str(p)))\n",
|
||||
" T.append(time)\n",
|
||||
" if time > 1:\n",
|
||||
" break\n",
|
||||
" if time > 1/1000:\n",
|
||||
" repeat = 1\n",
|
||||
" return D, T"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"id": "b1ab0f3d-6ce6-4b8c-9d4d-bb78bb1ac432",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAINCAYAAACTVPhqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAv41JREFUeJzs3Qd8U1X7B/Bf926hLZ3sIVL2liUgKENAQAU34F5/B6Kor4qDF1AREMUXFyq4QASVKYiIIHuLiOzdXeieSf6f56SJbUnbtE2a9fu+//vPSXJ7c3JT6n1ynvMcN51OpwMREREREZEJ7qYeJCIiIiIiYsBAREREREQV4ggDERERERGViwEDERERERGViwEDERERERGViwEDERERERGViwEDERERERGViwEDERERERGVy7P8p0hotVpcvHgRQUFBcHNz40khIiIiIocnazdnZmYiJiYG7u4VjyEwYCjHvHnz1FZQUIATJ05Y43MiIiIiIrKpc+fOoX79+hXu46aT8ILKlZ6ejjp16qiTGRwczDNFRERERA4vIyMDDRo0wOXLlxESElLhvhxhqIQhDUmCBQYMRERERORMzEm556RnIiIiIiIqFwMGIiIiIiIqFwMGIiIiIiIqF+cwEBEREdmA1J0pKiqCRqPh+Ser8PLygoeHR42Pw4CBiIiIqJZJ2fb4+Hjk5OTw3JNVJzRLydTAwMAaHYcBQyXrMDDqJyIiIksvCnvq1Cn1za8smuXt7c3FYckqI1jJyck4f/48WrRoUaORBq7DYEaNWqlNK+sxsKwqERER1VReXp4KGBo1agR/f3+eULKa3NxcnD59Gk2aNIGvr2+1r3E56ZmIiIjIBtzdeRlGtl9jwRz8TSUiIiIionIxYCAiIiIionIxYCAiIiJyUBqtDttOpOLH/RfUrdy3pn79+uGpp56y2vEl317SaPbv3w97dLoW+9e4cWPMmTMH9oBVkoiIiIgc0NpD8XhtxWHEp+cZH4sO8cWU4XEY3CbaKq+5bNkyVdvfWho0aKDKzYaHh8MeNbDz/lkLRxiIiIiIHDBYeOTLvaWCBZGQnqcel+etITQ0FEFBQVZbm0JKf0ZFRcHT0z6/0/aw8/5ZCwMGcji1PfxKRERUGzXzcwqKzNoy8wox5ae/YOq/fobHXv3psNrPnOPJa1cnJemDDz5Q9f2lXGdkZCRuueWWUvs9/vjjapPSnfKN/Msvv1zqtSTl5o033sA999yjyno++OCDV6T8/Pbbb+r+zz//jI4dO8LPzw/XXXcdkpKSsGbNGrRq1Ur97B133FFqETxZ62L69OmqnKj8TPv27bF06VKz3uOlS5dw5513ol69eupnW7Rogc8++0w9Z6n+mXN+yrp8+TLuv/9+1S85przOgQMHUBtcKzwih2eL4VciIiJryy3UIO6Vny1yLLnkTMjIQ9tX15m1/+HXB8Hfu2qXhLt378YTTzyBRYsWoWfPnkhLS8PmzZtL7fPFF1/gvvvuw86dO9X+EhA0bNgQDzzwgHGfmTNn4pVXXsGUKVMqfL1XX30V77//vlq3YsyYMWrz8fHB119/jaysLIwaNQrvvfceJk+erPaXYOHLL7/E/Pnz1QX/77//jrvuuktdbPft27fC15IL98OHD6sLfrmQP378uFrPwJL9M/f8lHTrrbeqgET6JUHGhx9+iAEDBuDo0aNq5MeaGDCQww2/lo29DcOv/7urE4MGIiKiWnD27FkEBARg2LBhKkVJFqGTb9jL5vvPnj1bfQPfsmVL/Pnnn+p+yQti+Zb8mWeeMd6Xb/BNmTp1Knr16qXacpH9wgsv4MSJE2jatKl6TEY3Nm7cqC7I8/PzMW3aNPzyyy/o0aOHel7227Jli7rIrixgkPcm76VLly7GkZDKVKV/VTk/BtJ3CSxk5EICEUOw9cMPP6iREwk2rIkBAzkESTuSkYXyhl9lWRJ5/vq4KHi4W2aREiKnp9UAZ7YCWYlAYCTQqCfg7mHrXhG5JD8vD/VNvzl2nkrD+M92Vbrf5xO6oluTULNeu6quv/56FSTIBfHgwYPVJt+il1y5+pprrim1cJhcvL/zzjvQaDRqLoAwXJRXpl27dsa2pD/J6xguxg2PyQW1kBEBSf+RPpadI1E2qDHlkUcewc0334y9e/fihhtuwMiRI9UoiqX6V5XzYyCpRzJSERYWVupxGfmQwMTaGDCUY968eWqTD41sT/44lp3YVTZokOdlvx7NSv9jIiITDv8ErJ0MZFz897HgGGDwm0DcCJ4yolomF47mpgX1aVFPpePKCLupL9LkEjQqxFftZ60v0WRUQS6oJYd/3bp1Kq1I0nJ27dqFOnXqmH0cGaUwR8nKTHKuylZqksdk3oKQC2uxatUqxMbGltrP8O18RYYMGYIzZ85g9erVWL9+vUr7eeyxx9Q3+pboX3XIe4qOjlbnu6yqnO/q4qTncsgvhuSvyS8+2V5SZp5F9yNyaRIsLLmndLAgMuL1j8vzRGS3JAiQuXuibDhguC/PW3vEXSoFDRw4EG+99RYOHjyo0ol+/fVX4/M7duwotf/27dvVfIKy355bWlxcnAoMJLWoefPmpTZJAzKHzHUYN26cmgcxZ84cfPTRRxbvZ1XOT6dOnZCQkKDOedn3VBslXjnCQA4hIsjXovsRuXQakowsVJTgt/Z54OobmZ5EZMek0IfM3StbCCSqlgqBrFy5EidPnsS1116LunXrqm/j5Rt0ycU3kAv2iRMn4qGHHlKjETLpV1JurE1GPyZNmoSnn35a9al3795IT0/HH3/8oaoLSSBQERkt6dy5M1q3bq3mQ6xcuVJVO7K0qpwfCcwkZUnSoyRAu+qqq3Dx4kU1iiKpYOamdlUXAwZyCJKDKcOv5aUlGYZfzcnVJHJpMmeh7MhCKTog44J+vyZ9arFjRFRVEhTI3D1Jx5URdvnSTP47WBtz+SQNRhZxkzSkvLw89c34N998oy6yDaRcquTYd+vWTX1r/uSTT1p9cq6BlGuVUQKpliSBjfRXvqV/8cUXK/1Zb29vNWlZRkykKlGfPn3w7bffWryPVTk/ktIkQdl//vMfTJgwAcnJyWo9CAnYZH6EtbnpqlJ81wVlZGSo0lUSmUpUSratkvTwl3uveNzwZ5FVkojM8OdS4Pv7Kt/v5k+Btv/WVCciy5EL7FOnTqk1AmQNA2ck6wx06NBBpfOQ7c5PRb9rVbnG5RwGchj9WkbAz+vKX1kZWWCwQGQmqYZkyf2IiMjpMSWJHMYvfycit1CL6GAfvDOmA5Kz8mt1+JXIKUS1BTy8AU1BOTu46aslSYlVIiIn9PDDD6vJzKbI4m6y2BuVxoCBHMbSPefV7c2dG6Bnc+tXBCByOjlpwFe3VBwsiMEzOOGZiGrEVPlPe/H666+rSdGm1Fb6+W92fH5MYcBADiEpIw+/H01W7dGdStdUJiIzZCYCi0YCSYcBv7pAr6eAnR+aWIdhBtdhICKnFhERoTYyHwMGcgg/7L8ArQ7o3KgumtYLtHV3iBzL5bPAwpuAtJNAYBRwzw9ARCug5/9xpWciIqoUAwaye1LIy5iO1Km+rbtD5FhSjumDBSmVWqchcM+PQGhT/XPuHiydSkRElWLAQHbv0IUMHE3MgrenO25sZ92FaIicSvxBYNEoICcFCL9KHyxI2hEREVEVMGAgu7d0zzl1O6h1FEL8vGzdHSLHcHYH8NWtQH46ENUOuHs5EMBiAUREVHUusQ6DLOktS5XLKoSffPKJrbtDVZBfpMGPB/STMm/mZGci85zYqJ/gLMFCwx7A+JUMFoiIqNqcPmAoKirCxIkT8euvv2Lfvn14++23kZqaautukZk2HknC5ZxCRAb7oE+LejxvRJX5eyXw9RigMAdodh1w1zLAN4TnjchZaTXAqc36VdzlVu5beYXip556ymrHP336NNzc3LB//37Y0h9//IG2bdvCy8sLI0eOhKtz+pSknTt3onXr1oiN1ZfiHDJkCNatW4fbb7/d1l0jMyzdc0HdjuwYy8XZiCpzYDHwwyOATgO0Gg7c/Cng6cPzRuSsDv8ErJ1sojzym1Yrj7xs2TJ1EW0tDRo0QHx8PMLDbZtCKV82d+jQAWvWrEFgYM2rM7766qv44YcfbB4IOe0Iw++//47hw4cjJiZGRZxyssuaN28eGjduDF9fX3Tv3l0FCQYXL140BgtC2hcu6C9Cyb6lZOXjt3+SVPsWVkciqtiuT4DlD+qDhfZ3ALd8zmCByNmDhSX3lA4WREa8/nF53gpCQ0MRFBRklWMXFBTAw8MDUVFR8PS07XfaJ06cwHXXXYf69eujTp06sBdyjmzB7gOG7OxstG/fXgUFpixevFhFgVOmTMHevXvVvoMGDUJSkv5CkxzXj/svokirQ/v6IWgRaZ0/TkROYfMsYNUz+na3h4Cb5gEeTj+ATORcdDqgINu8LS8DWPOc/JCpA+lvZORB9jPnePLa1UhJ+uCDD9T8UPnCNjIyErfcckup/R5//HG1hYSEqBGDl19+WZVKN5Ave9944w3cc889aoXlBx988IqUJFkRWe7//PPP6NixI/z8/NSFvFznybf/rVq1Uj97xx13ICcnx3hsrVaL6dOno0mTJupn5Ppw6dKllb4/w+tL+vq9996r2p9//jk0Gg3uu+8+4/Fkbuy7775b6melr926dUNAQIAKMnr16oUzZ86on3/ttddw4MABdTzDMcXly5dx//33o169eup9yHuT/UqOTMhIh8zBldeWc20Ldv9fFEkhkq08s2bNwgMPPIAJEyao+/Pnz8eqVauwYMECPP/882pkouSIgrTlwyxPfn6+2gwyMjLUbWFhoXFOhJDIVx6TD71s2xAhy1a2Lc+7u7ur15AhPVNtb29v9Rqyf8m2j4+P+gcgr2WqLX2T/eWXWraybXle/qHKa5Vt2+N7MlRHuqVzfad5T874OfE92fBzkvbmGXD7Y7b+D1afSdD2e5GfE/8Line truncated
|
||||
"text/plain": [
|
||||
"True"
|
||||
"<Figure size 900x600 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"isprime_fermat(561, k=500)"
|
||||
"plot_run_times()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "83f69e44-4f3b-43b8-985c-0db82e983156",
|
||||
"id": "408ba0fc-9bd9-4929-8ff3-1f9582ffa483",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"561 is a [Carmichael number](https://en.wikipedia.org/wiki/Carmichael_number), a number that is resistant to the Fermat test.\n",
|
||||
"We can test that our two Fermat functions work on the big primes and on the carmichael numbers, and get some more timing data:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"id": "7eacb754-e6ee-482e-82fd-8725845aad7e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 1.2 ms, sys: 1 μs, total: 1.2 ms\n",
|
||||
"Wall time: 1.2 ms\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%time test_prime_checker(isprime, primes=big_primes, composites=carmichael_numbers);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 23,
|
||||
"id": "389cc455-6b9c-40f6-9fc2-2ae129159761",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 3.43 ms, sys: 0 ns, total: 3.43 ms\n",
|
||||
"Wall time: 3.43 ms\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%time test_prime_checker(isprime_fermat, primes=big_primes, composites=carmichael_numbers);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1347399a-24f2-4d5f-aa7f-822951126bf5",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"What have we learned about timing?\n",
|
||||
"- `isprime_simple` has a run time of *O*(*n*)\n",
|
||||
"- `isprime_faster` has a run time of *O*(√*n*); we see its slope is half as steep on the log-log graph\n",
|
||||
"- `isprime_fermat` and `isprime` ([according to Wikipedia](https://en.wikipedia.org/wiki/Fermat_primality_test)) have a run time of *O*(*k* log<sup>2</sup>*n*), where *k* is the number of repeats \n",
|
||||
"- `isprime_faster` is more than 10,000 times faster than `isprime_simple` on 8-digit primes\n",
|
||||
"- `isprime_fermat` is more than 10,000 times faster than `isprime_faster` on 17 digit primes\n",
|
||||
"- `isprime` is about twice as fast as `isprime_fermat`; but that's quite good for `isprime_fermat` considering how simple it is\n",
|
||||
"- `isprime` produces no known incorrect answers; `isprime_fermat` will on rare occasion produce a false prime\n",
|
||||
"\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!"
|
||||
"I hope you have enjoyed this excursion into the land of prime numbers.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
|
||||
Reference in new issue
Block a user