diff --git a/ipynb/TruncatablePrimes.ipynb b/ipynb/TruncatablePrimes.ipynb new file mode 100644 index 0000000..9133eab --- /dev/null +++ b/ipynb/TruncatablePrimes.ipynb @@ -0,0 +1,689 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8812d495-6014-4669-97ae-003be4f54605", + "metadata": {}, + "source": [ + "
Peter Norvig
Jan 2026
\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//c/portal/getImageAttachment?filename=IMG_20181212_120939.jpg&userId=143131)](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+0b1795g0aVLs2LGjyT5btmyJCRMmRNeuXfPzzJo1Kw4ePNgy7wgAaFsBZe3atTl8bNiwIVavXh0HDhyIsWPHxuuvv17dZ8aMGfHII4/EihUr8v7btm2LiRMnVrcfOnQoh5P9+/fHunXr4r777otly5bFvHnzWvadAQBtYxTPY4891mQ5BYtUA7Jx48a45JJLYvfu3XHvvffG8uXL47LLLsv7LF26NM4+++wcakaNGhWPP/54PP/88/HEE09Ev379YsSIEXHbbbfFjTfeGLfcckt06tSpZd8hANC2+qCkQJL07t07z1NQSbUqY8aMqe4zdOjQGDRoUKxfvz4vp/mwYcNyOKkYN25c7NmzJzZt2nTM19m3b1/efuQEANSvZgeUw4cPxw033BAf/ehH49xzz83rtm/fnmtAevXq1WTfFEbStso+R4aTyvbKtrfq+9KzZ8/qNHDgwOYWGwCo54CS+qL86U9/igcffDBOtjlz5uTamsq0devWk/6aAEAru5Ls9ddfHytXroynn346zjjjjOr6/v37586vu3btalKLkkbxpG2VfX772982eb7KKJ/KPkfr3LlzngCAtuG4alAaGxtzOHn44YfjySefjDPPPLPJ9pEjR0bHjh1jzZo11XVpGHIaVjx69Oi8nObPPfdc7Ny5s7pPGhHUo0ePOOecc078HQEAbasGJTXrpBE6P//5z/O1UCp9RlK/kFNOOSXPp0yZEjNnzswdZ1PomD59eg4laQRPkoYlpyBy7bXXxqJFi/JzzJ07Nz+3WhIA4LgDyt13353nl156aZP1aSjxF77whfzz4sWLo6GhIV+gLY2+SSN07rrrruq+7du3z81DU6dOzcGlW7duMXny5Jg/f75PpI3dgC1xEzYATjigpCaet9OlS5dYsmRJnt7K4MGD49FHHz2elwYA2hD34gEAiiOgAADFEVAAgOIIKABAcQQUAKA4AgoAUBwBBQAojoACABRHQAEAiiOgAADFEVAAgOIIKABAcQQUAKA4AgoAUBwBBQAojoACABRHQAEAiiOgAADFEVAAgOIIKABAcQQUAKA4AgoAUBwBBQAojoACABRHQAEAiiOgAADFEVAAgOIIKABAcQQUAKA4AgoAUBwBBQAojoACABRHQAEAiiOgAADFEVAAgOIIKABAcQQUAKA4AgoAUBwBBQAojoACABRHQAEAiiOgAADFEVAAgOIIKABAcQQUAKA4AgoAUBwBBQAoTodaFwCgtRkye1WzH/vywgktWhaoV2pQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgC0/oDy9NNPx5VXXhkDBgyIdu3axc9+9rMm2xsbG2PevHlx+umnxymnnBJjxoyJF154ock+//73v+Oaa66JHj16RK9evWLKlCnx2muvnfi7AQDaZkB5/fXXY/jw4bFkyZJjbl+0aFHccccdcc8998QzzzwT3bp1i3HjxsXevXur+6RwsmnTpli9enWsXLkyh54vfelLJ/ZOAIC2e6G28ePH5+lYUu3J7bffHnPnzo2rrroqr/vRj34U/fr1yzUtn/nMZ+LPf/5zPPbYY/G73/0uzj///LzPnXfeGVdccUX8z//8T66ZAQDathbtg/LSSy/F9u3bc7NORc+ePePCCy+M9evX5+U0T806lXCSpP0bGhpyjcux7Nu3L/bs2dNkAgDqV4te6j6FkyTVmBwpLVe2pXnfvn2bFqJDh+jdu3d1n6MtWLAgbr311pYsapt2IpfpTlyqG4CTrVWM4pkzZ07s3r27Om3durXWRQIAWktA6d+/f57v2LGjyfq0XNmW5jt37myy/eDBg3lkT2Wfo3Xu3DmP+DlyAgDqV4sGlDPPPDOHjDVr1lTXpf4iqW/J6NGj83Ka79q1KzZu3Fjd58knn4zDhw/nvioAAMfdByVdr+TFF19s0jH22WefzX1IBg0aFDfccEN885vfjPe///05sNx00015ZM7VV1+d9z/77LPjk5/8ZFx33XV5KPKBAwfi+uuvzyN8jOABAJoVUH7/+9/Hxz/+8eryzJkz83zy5MmxbNmy+PrXv56vlZKua5JqSi666KI8rLhLly7Vx9x///05lFx++eV59M6kSZPytVMAAJoVUC699NJ8vZO3kq4uO3/+/Dy9lVTbsnz5cp8AANB6R/EAAG2LgAIAFEdAAQCKI6AAAMURUACA4ggoAEBxBBQAoDgCCgBQHAEFACiOgAIAFEdAAQBa/714AGg5Q2avavZjX144wUdB3VKDAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiCCgAQHEEFACgOAIKAFAcAQUAKI6AAgAUR0ABAIrTodYFAKBlDJm9qtmPfXnhBB8DRVGDAgAURw1KK+K/IwDaCjUoAEBxBBQAoDgCCgBQHAEFACiOgAIAFEdAAQCKI6AAAMURUACA4ggoAEBxXEkWgDdx5WpqTQ0KAFAcAQUAKI6AAgAUR0ABAIojoAAAxRFQAIDiGGYMwEllyDLNIaCcZH4xAeD4aeIBAIojoAAAxalpE8+SJUviO9/5Tmzfvj2GDx8ed955Z1xwwQW1LBIABdNs3nbULKA89NBDMXPmzLjnnnviwgsvjNtvvz3GjRsXmzdvjr59+9aqWAC0EcJO2WoWUL773e/GddddF1/84hfzcgoqq1atih/+8Icxe/bsqCUnLQC1+LtxIs9z9HO1djUJKPv374+NGzfGnDlzqusaGhpizJgxsX79+jftv2/fvjxV7N69O8/37NlzUsp3eN8bzX7s0WUq8blO5HnawnO1hs/Qczlezon6/B1qye/BElXK19jY+PY7N9bA3//+91SyxnXr1jVZP2vWrMYLLrjgTfvffPPNeX+TY+AccA44B5wDzoFo9cdg69atb5sVWsV1UFJNS+qvUnH48OH497//HX369Il27dq9bVobOHBgbN26NXr06PEOlBbHvvac9457W+Ocbx3HPtWcvPrqqzFgwIC3fd6aBJR3vetd0b59+9ixY0eT9Wm5f//+b9q/c+fOeTpSr169jus100ETUGrDsa8dx95xb2uc8+Uf+549e5Z7HZROnTrFyJEjY82aNU1qRdLy6NGja1EkAKAgNWviSU02kydPjvPPPz9f+yQNM3799dero3oAgLarZgHl05/+dPzjH/+IefPm5Qu1jRgxIh577LHo169fi75Oahq6+eab39RExMnn2NeOY++4tzXO+fo79u1ST9kWfUYAgBPkXjwAQHEEFACgOAIKAFAcAQUAKE7dB5QlS5bEkCFDokuXLvmuyb/97W9rXaS6d8stt+Qr/B45DR06tNbFqjtPP/10XHnllfmKjOkY/+xnP2uyPfV/T6PkTj/99DjllFPyva5eeOGFmpW3LR37L3zhC2/6HfjkJz9Zs/LWiwULFsRHPvKROPXUU/Nd76+++urYvHlzk3327t0b06ZNy1ca7969e0yaNOlNFwXl5Bz7Sy+99E3n/Ze//OVorroOKA899FC+3koa/vSHP/whhg8fHuPGjYudO3fWumh170Mf+lC88sor1enXv/51rYtUd9J1g9I5nUL4sSxatCjuuOOOfKfwZ555Jrp165bP//QFzsk99kkKJEf+DjzwwAMO+wlau3ZtDh8bNmyI1atXx4EDB2Ls2LH586iYMWNGPPLII7FixYq8/7Zt22LixImO/Ttw7JPrrruuyXmfvoearbGOpRsPTps2rbp86NChxgEDBjQuWLCgpuWqd+nmjsOHD691MdqU9Kv88MMPV5cPHz7c2L9//8bvfOc71XW7du1q7Ny5c+MDDzxQo1K2jWOfTJ48ufGqq66qWZnaip07d+bjv3bt2uo53rFjx8YVK1ZU9/nzn/+c91m/fn0NS1r/xz752Mc+1vjVr361saXUbQ3K/v37Y+PGjblau6KhoSEvr1+/vqZlawtSU0Kq/n7ve98b11xzTWzZsqXWRWpTXnrppXwBxCPP/3T/i9TM6fx/Zzz11FO5KvyDH/xgTJ06Nf71r3+9Q6/cduzevTvPe/funefpOz/9Z3/keZ+alwcNGuS8P8nHvuL+++/P99s799xz841+33jjjWa/Rqu4m3Fz/POf/4xDhw696cq0afkvf/lLzcrVFqQ/gsuWLctfzKmK79Zbb42LL744/vSnP+X2S06+FE6SY53/lW2cPKl5JzUrnHnmmfHXv/41vvGNb8T48ePzH8l0o1ROXLp/2w033BAf/ehH8x/DJJ3b6V5vR99M1nl/8o998rnPfS4GDx6c/zn94x//GDfeeGPup/LTn/60Wa9TtwGF2klfxBXnnXdeDizppP3JT34SU6ZM8dFQ9z7zmc9Ufx42bFj+PXjf+96Xa1Uuv/zympatXqT+EOmfHv3byjn2X/rSl5qc96mDfjrfU0hP5//xqtsmnlTFlP5TObr3dlru379/zcrVFqX/Zj7wgQ/Eiy++WOuitBmVc9z5X4bU1Jm+k/wOtIzrr78+Vq5cGb/61a/ijDPOaHLep+b9Xbt2Ndnf9/7JP/bHkv45TZp73tdtQEnVfCNHjow1a9Y0qZZKy6NHj65p2dqa1157LSfolKZ5Z6SmhfRlfeT5v2fPnjyax/n/zvvb3/6W+6D4HTgxqU9y+gP58MMPx5NPPpnP8yOl7/yOHTs2Oe9TE0PqA+e8P7nH/lieffbZPG/ueV/XTTxpiPHkyZPj/PPPjwsuuCBuv/32PCTqi1/8Yq2LVte+9rWv5WtEpGadNMQvDfNOtVmf/exna120ugt+R/5nkjrGpi+E1GktdQpMbcTf/OY34/3vf3/+Mrnpppty23C6fgEn79inKfW7StffSCExhfOvf/3rcdZZZ+Vh3pxY08Ly5cvj5z//ee7PVulPlTqAp2v9pHlqRk7f/elz6NGjR0yfPj2Hk1GjRjn0J/HYp/M8bb/iiivyNWhSH5Q05PuSSy7JTZzN0ljn7rzzzsZBgwY1durUKQ873rBhQ62LVPc+/elPN55++un5mL/nPe/Jyy+++GKti1V3fvWrX+VhfkdPaYhrZajxTTfd1NivX788vPjyyy9v3Lx5c62LXffH/o033mgcO3Zs47vf/e485HXw4MGN1113XeP27dtrXexW71jHPE1Lly6t7vOf//yn8Stf+Urjaaed1ti1a9fGT33qU42vvPJKTcvdFo79li1bGi+55JLG3r175++bs846q3HWrFmNu3fvbvZrtvv/LwwAUIy67YMCALReAgoAUBwBBQAojoACABRHQAEAiiOgAADFEVAAgOIIKABAcQQUAKA4AgoAUBwBBQAojoACAERp/h8ihJnT0b/I0wAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "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+vWrFkTffr0KfemAIAKVPb4uOKKK2L58uX5YZd169bFggUL4tZbb42LL7643JsCACpQ2eNjyJAhsWjRoli4cGEMHDgwrr322pg9e3aMHz++3JsCACpQ2c/zkTn99NPzBQDg47y3CwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACTVLu3mAKA4+k57OCrRazNH7dLtG/kAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAihUfM2fOjKqqqrj88stbe1MAQFuPj+eeey5uueWWOOKII1pzMwBABWm1+NiyZUuMHz8+brvttth7771bazMAQIVptfi4+OKLY9SoUTF8+PD/uF5TU1M0Nja2WACA4mrXGt/07rvvjpUrV+aHXT5LbW1tzJgxI1LpO+3hqESvzRy1q3cBAL6YIx/19fVx2WWXxV133RXt27f/zPVramqioaGhecm+HgAorrKPfKxYsSI2b94cRx99dPN1H374YSxbtixuvvnm/DDL7rvv3vy56urqfAEA2oayx8epp54aL730UovrJk6cGP37948rr7yyRXgAAG1P2eOjU6dOMXDgwBbXdezYMbp16/aJ6wGAtscZTgGAyn+1y8ctXbo0xWYAgApg5AMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAKjs+KitrY0hQ4ZEp06donv37jF27NhYvXp1uTcDAFSossfHk08+GRdffHEsX748Hnvssfjggw/itNNOi61bt5Z7UwBABWpX7m+4ZMmSFpfnz5+fj4CsWLEiTj755HJvDgBo6/HxcQ0NDfnHrl277vDzTU1N+bJdY2Nja+8SAFDUCafbtm2Lyy+/PE488cQYOHDgp84R6dKlS/PSq1ev1twlAKDI8ZHN/Xj55Zfj7rvv/tR1ampq8tGR7Ut9fX1r7hIAUNTDLpdcckk89NBDsWzZsjjggAM+db3q6up8AQDahrLHR6lUiksvvTQWLVoUS5cujX79+pV7EwBABWvXGodaFixYEPfff39+ro9Nmzbl12fzOTp06FDuzQEAbX3Ox9y5c/O5G0OHDo0ePXo0L/fcc0+5NwUAVKBWOewCAPBpvLcLAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwBQjPiYM2dO9O3bN9q3bx/HHXdcPPvss621KQCgrcfHPffcE1OmTImrr746Vq5cGYMHD44RI0bE5s2bW2NzAEBbj48bb7wxLrjggpg4cWIMGDAg5s2bF1/60pfiN7/5TWtsDgCoIO3K/Q3ff//9WLFiRdTU1DRft9tuu8Xw4cPjmWee+cT6TU1N+bJdQ0ND/rGxsTFaw7amd6MS7ezPw+0szv1Zqfdlxu0szv3pd1Dbvj935nuWSqXPXrlUZq+//nq21dLTTz/d4vqpU6eWjj322E+sf/XVV+frW/wMPAY8BjwGPAY8BqLifwb19fWf2QplH/nYWdkISTY/ZLtt27bFv/71r+jWrVtUVVVFpciKr1evXlFfXx+dO3eOonI7i6Ut3J9t4TZm3M5iaazAx2024vHOO+9Ez549P3PdssfHPvvsE7vvvnu8+eabLa7PLu+///6fWL+6ujpfPmqvvfaKSpU9SCrlgfK/cDuLpS3cn23hNmbczmLpXGGP2y5duuyaCad77rlnfO1rX4vHH3+8xWhGdvn4448v9+YAgArTKoddssMoEyZMiGOOOSaOPfbYmD17dmzdujV/9QsA0La1Snycc8458c9//jOuuuqq2LRpUxx55JGxZMmS2G+//aKoskNH2XlNPn4IqWjczmJpC/dnW7iNGbezWKoL/ritymad7uqdAADaDu/tAgAkJT4AgKTEBwCQlPgAAJISH/+jZcuWxejRo/MzumVnZF28eHEUUW1tbQwZMiQ6deoU3bt3j7Fjx8bq1aujSObOnRtHHHFE80l9svPSPPLII1F0M2fOzB+7l19+eRTJNddck9+ujy79+/ePInr99dfjvPPOy88M3aFDhxg0aFDU1dVFkfTt2/cT92e2XHzxxVEUH374YUyfPj369euX348HHXRQXHvttZ/vvVIqzC4/vXqly85fMnjw4Jg0aVKceeaZUVRPPvlk/o88C5B///vf8eMf/zhOO+20+Nvf/hYdO3aMIjjggAPyP8QHH3xw/o/9jjvuiDFjxsTzzz8fhx9+eBTRc889F7fcckseXUWU3W9//OMfmy+3a1e8X3lvvfVWnHjiiTFs2LA8lvfdd99Yu3Zt7L333lG0x2r2x3m7l19+Ob7xjW/EWWedFUUxa9as/ElQ9rsne+xmAZmdHys7a+jkyZN39e6VVfH+JSY2cuTIfCm67DwtHzV//vx8BCR7B+OTTz45iiAbwfqo66+/Pv9FsHz58kLGx5YtW2L8+PFx2223xXXXXRdFlMXGjt7WoUiyP1jZe4Dcfvvtzddlz5yLJouqj8qeKGQjA6ecckoUxdNPP50/4Rk1alTzaM/ChQvj2WefjaJx2IX/SkNDQ/6xa9euhfwJZs+w7r777nxkq6hvC5CNZGW/5IYPHx5FlY0AZIdEDzzwwDy0Nm7cGEXzwAMP5GeTzkYAsicERx11VB6URfb+++/H7373u3zEuZLegPSznHDCCflbkaxZsya//MILL8RTTz1VyCe4Rj7Yadl79WTzA7Kh3oEDBxbqJ/jSSy/lsfHee+/Fl7/85Vi0aFEMGDAgiiYLq5UrV+ZD2UV13HHH5SN0hx56aPzjH/+IGTNmxP/93//lw/XZ3KWiePXVV/MRuuxtLbLDodl9mg3RZ++zlb3NRRFlc+vefvvt+M53vhNFMm3atPzdbLO5SdkbtGZPgrIR2CycCyc7wynlkf04Fy1aVPgf50UXXVTq06dPqb6+vlQ0TU1NpbVr15bq6upK06ZNK+2zzz6lv/71r6Ui2bhxY6l79+6lF154ofm6U045pXTZZZeViuytt94qde7cufTrX/+6VCR77LFH6fjjj29x3aWXXlr6+te/Xiqq0047rXT66aeXimbhwoWlAw44IP/44osvlu68885S165dS/Pnzy8VjZEPdsoll1wSDz30UP4qn2yCZtFkzxa/+tWv5v+fvTtz9izypptuyidlFkU2T2fz5s1x9NFHN1+XPcPK7tObb745mpqa8mddRbPXXnvFIYccEuvWrYsi6dGjxydG5w477LC49957o4g2bNiQTyK+7777omimTp2aj36MGzcuv5y9aim7vdmrDYs2iiU++FyygZ1LL700PwyxdOnSQk5o+7RDTNkf4yI59dRT88NLH5XNqM+Geq+88spChsf2CbavvPJKnH/++VEk2eHPj7/sPZsz0KdPnyiibGJtNrdl+6TMInn33Xdjt91aTsXM/j1mv4eKRnyU4RfaR59JrV+/PlatWpVPxOzdu3cUaXLiggUL4v7778+Pl2fvVpzJXgKWvR69CGpqavKJXdn99s477+S3NwutRx99NIoku/8+Plcne7l0do6IIs3h+eEPf5i/gin7I/zGG2/k7xCa/SI/99xzo0iuuOKKfKLiDTfcEGeffXb+yohbb701X4om+yOcxUc2ClDEl02PHj06n+OR/Q7KXmGXvcz/xhtvzCfWFs6uPu5T6Z544ol8rsfHlwkTJpSKZEe3MVtuv/32UlFMmjQpn8uy5557lvbdd9/SqaeeWvrDH/5QaguKOOfjnHPOKfXo0SO/P7/yla/kl9etW1cqogcffLA0cODAUnV1dal///6lW2+9tVREjz76aP57Z/Xq1aUiamxszP8d9u7du9S+ffvSgQceWPrJT36Sz0UrmqrsP7s6gACAtsN5PgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBApPT/ABTjgt70nMKKAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "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*(*n* - 1) ≡ 1 (mod *n*). \n", + "\n", + "That means we can check if *n* is prime by choosing a random *a* and testing if *a*(*n* - 1) ≡ 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 264 (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 +}