Add files via upload

2026-01-23 11:47:12 -08:00
1 parent d134b4220c
commit a4b4f0aa09
@@ -25,7 +25,7 @@
    "\n",
    "Numbers like this are called [**truncatable primes**](https://en.wikipedia.org/wiki/Truncatable_prime). \n",
    "\n",
-    "I thought I would write a program to find all the truncatable primes. My function `truncatable_primes` below starts with the list of one-digit primes: [2, 3, 5, 7]. Then it creates a list of one-digit-longer new primes by placing each possible digit 1–9 to the left of each of the one-digit primes, and testing if each resulting number is prime. 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 returns all the truncatable primes in sorted order."
+    "I thought I would write a program to find all the truncatable primes. (It is not immediately obvious that there is a finite number of them, but it turns out to be true.) My function `truncatable_primes` below starts with the list of one-digit primes, [2, 3, 5, 7], and on each iteration creates a list of new numbers that are one digit longer by placing each possible digit 1–9 to the left of each prime. Then out of those numbers we find just the new primes.  If there are any new primes, the function recursively adds digits to them, giving us three-digit primes, then four-digit primes, and so on, stopping when there are no more new primes. In the end, the function returns all the truncatable primes in sorted order."
   ]
  },
  {
@@ -38,8 +38,9 @@
    "from sympy import isprime # isprime(n) returns True if n is a prime\n",
    "\n",
    "def truncatable_primes(primes=[2, 3, 5, 7]) -> list[int]:\n",
-    "    \"\"\"All truncatable primes, in ascending order.\"\"\"\n",
-    "    new_primes = [dp for d in \"123456789\" for p in primes if isprime(dp := int(d + str(p)))]\n",
+    "    \"\"\"Truncatable primes, in ascending order.\"\"\"\n",
+    "    new_numbers = (int(d + str(p)) for d in \"123456789\" for p in primes)\n",
+    "    new_primes = [n for n in new_numbers if isprime(n)]\n",
    "    if new_primes:\n",
    "        return primes + truncatable_primes(new_primes)\n",
    "    else:\n",
@@ -147,6 +148,8 @@
   "id": "e82700e8-8fb2-45cf-bcf0-7805e46174bb",
   "metadata": {},
   "source": [
+    "## Right-Truncatable Primes\n",
+    "\n",
    "What if you sharpen the pencil from the other end? The primes we found so far are called **left-truncatable primes**; there are also **right-trunctable primes**:"
   ]
  },
@@ -160,7 +163,8 @@
    "def right_truncatable_primes(primes=[2, 3, 5, 7]) -> list[int]:\n",
    "    \"\"\"All right-truncatable primes, in ascending order.\"\"\"\n",
    "    # Only consider (1, 3, 7, 9) as the digit on the right; placing any other digit forms a composite number\n",
-    "    new_primes = [pd for p in primes for d in (1, 3, 7, 9) if isprime(pd := 10 * p + d)]\n",
+    "    new_numbers = (10 * p + d for p in primes for d in (1, 3, 7, 9))\n",
+    "    new_primes = [n for n in new_numbers if isprime(n)]\n",
    "    if new_primes:\n",
    "        return primes + right_truncatable_primes(new_primes)\n",
    "    else:\n",
@@ -287,7 +291,9 @@
   "source": [
    "There are only 9 of them, with 373 the largest.\n",
    "\n",
-    "What about primes that are both truncatable and [**palindromic primes**](https://en.wikipedia.org/wiki/Palindromic_prime)?"
+    "## Palindromic Primes\n",
+    "\n",
+    "What about primes that are both truncatable and [**palindromic**](https://en.wikipedia.org/wiki/Palindromic_prime) (the same when you reverse the digits)?"
   ]
  },
  {
@@ -316,7 +322,139 @@
   "id": "cd0da3ad-54c0-47dc-979c-4eecc935425c",
   "metadata": {},
   "source": [
-    "There are 13 **palindromic truncatable primes**."
+    "There are 13 **palindromic truncatable primes**, with only 4 of them longer than 3 digits.\n",
+    "\n",
+    "**Aside**: Now I'm curious about [**palindromic primes**]((https://en.wikipedia.org/wiki/Palindromic_prime) that are not truncatable. To find them, I'll enumerate the \"first half\" of numbers, and then concatenate the reverse of the first half, forming both odd- and even-length numbers, e.g. from 123, the function `palindromes` produces the palindromes 12321 and 123321. Then all we have to do is check each palindrome for primality."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "fc208a86-c3f4-4426-913d-62ec0a639576",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def palindromic_primes(digit_length: int) -> list[int]:\n",
+    "    \"\"\"Generate all palindromic primes of up to `digit_length` digits.\"\"\"\n",
+    "    end = 10 ** ((digit_length + 1) // 2)\n",
+    "    return [n for half in range(1, end)\n",
+    "              for n in palindromes(half)\n",
+    "              if isprime(n)]\n",
+    "\n",
+    "def palindromes(half: int) -> tuple[int]:\n",
+    "    \"\"\"Given the left half of a palindrome, make two full palindromes,\n",
+    "    one without the middle digit repeated, one with.\"\"\"\n",
+    "    left = str(half)\n",
+    "    right = left[::-1]\n",
+    "    return int(left + right[1:]), int(left + right)\n",
+    "\n",
+    "assert palindromes(123) == (12321, 123321)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ec5f9d5a-c89e-4dc9-812b-232437e4801e",
+   "metadata": {},
+   "source": [
+    "Let's count the palindromic primes of up to 11 digits, and peek at some of them:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "3f9c6d25-8d40-4267-a45f-c0284b910c86",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "47995"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "PP = palindromic_primes(11)\n",
+    "\n",
+    "len(PP)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "69149eb4-d79d-4cfa-8e48-40946e54c48f",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[11,\n",
+       " 2,\n",
+       " 3,\n",
+       " 5,\n",
+       " 7,\n",
+       " 101,\n",
+       " 131,\n",
+       " 151,\n",
+       " 181,\n",
+       " 191,\n",
+       " 313,\n",
+       " 353,\n",
+       " 373,\n",
+       " 383,\n",
+       " 727,\n",
+       " 757,\n",
+       " 787,\n",
+       " 797,\n",
+       " 919,\n",
+       " 929,\n",
+       " 10301,\n",
+       " 10501,\n",
+       " 10601,\n",
+       " 11311,\n",
+       " 11411]"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "PP[:25]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "3bebaad1-8fb2-4892-87f8-8846062d92d4",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[99985958999,\n",
+       " 99986068999,\n",
+       " 99988688999,\n",
+       " 99988988999,\n",
+       " 99991819999,\n",
+       " 99994349999,\n",
+       " 99995759999,\n",
+       " 99998189999,\n",
+       " 99998989999,\n",
+       " 99999199999]"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "PP[-10:]"
   ]
  },
  {
@@ -326,12 +464,12 @@
   "source": [
    "# Digit Lengths\n",
    "\n",
-    "How many digits are in these truncatable primes? The function `digit_lengths` returns a {number_of_digits: count_of_primes_with_that_many_digits} dict, and plots a bar chart for that data:"
+    "How many digits are in these primes? The function `digit_lengths` returns a {number_of_digits: count_of_primes_with_that_many_digits} dict, and plots a bar chart for that data:"
   ]
  },
  {
   "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 15,
   "id": "28b0bc97-1b22-449f-b33a-6f8c84fd890d",
   "metadata": {},
   "outputs": [],
@@ -343,7 +481,8 @@
    "    \"\"\"Plot a bar chart and return a dict of the number of digits in these primes.\"\"\"\n",
    "    digits = Counter(len(str(p)) for p in primes)\n",
    "    plt.bar(list(digits), list(digits.values()))\n",
-    "    plt.xlabel('Number of digits'); plt.ylabel('Count')\n",
+    "    plt.xlabel('Number of digits')\n",
+    "    plt.ylabel('Count')\n",
    "    return dict(digits)"
   ]
  },
@@ -357,7 +496,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 16,
   "id": "defd769e-9f92-4e4f-a12b-88a653b6efee",
   "metadata": {},
   "outputs": [
@@ -390,7 +529,7 @@
       " 24: 1}"
      ]
     },
-     "execution_count": 12,
+     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
@@ -419,7 +558,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 17,
   "id": "09f29951-d8f2-47a3-b666-bfa6fd5d1c87",
   "metadata": {},
   "outputs": [
@@ -429,7 +568,7 @@
       "{1: 4, 2: 9, 3: 14, 4: 16, 5: 15, 6: 12, 7: 8, 8: 5}"
      ]
     },
-     "execution_count": 13,
+     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
@@ -448,12 +587,52 @@
    "digit_lengths(TPr)"
   ]
  },
+  {
+   "cell_type": "markdown",
+   "id": "8618e5bc-ea92-4116-aa74-2380fbfd0dc3",
+   "metadata": {},
+   "source": [
+    "Finally, for the palindromic primes, we find that \"11\" is the only even-digit-length palindromic prime; it turns out that any even-digit-length palindromic number is [divisible by 11](https://www.geeksforgeeks.org/dsa/g-fact-all-palindromes-with-an-even-number-of-digits-are-divisible-by-11/). The counts of primes grow so quickly that it is better to use a log scale on the y axis: "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "id": "ed392e12-7480-4d96-994f-96cdb786f8f9",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{2: 1, 1: 4, 3: 15, 5: 93, 7: 668, 9: 5172, 11: 42042}"
+      ]
+     },
+     "execution_count": 18,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGwCAYAAACgi8/jAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJJdJREFUeJzt3QuUVdV9P/DfAAJqBEVSEAExVlF8gOEl5iWRSImSqm1iuhKD2Jh0FSIJ0RZihJhowTQSGkNKtUXa1VqJbUQjlahIglqUl2gMYqSFhoi8jPKy8pr5r326Zv6OPISZgXvZ9/NZ6yzuOffOnX23I/Nl79/ep6qmpqYmAAAy1KzUDQAAOFQEHQAgW4IOAJAtQQcAyJagAwBkS9ABALIl6AAA2WoRFa66ujrWrFkTxx13XFRVVZW6OQDAAUjbAG7ZsiU6deoUzZrte9ym4oNOCjldunQ5kD4FAMrM6tWro3Pnzvt8vmKDzpQpU4pj165ddR3Vpk2bUjcLADgAmzdvLgYq0ozM/lRV+i0gUke1bds2Nm3aJOgAQGa/vxUjAwDZEnQAgGwJOgBAtio26KRC5B49ekTfvn1L3RQA4BBRjKwYGQCOOIqRAYCKV7FTVwBA/gQdACBbFRt0FCMDQP4UIytGBoAjjmJkAKDiVezUFQCQP0EHAMiWoAMAZKtFqRsAABwa3cbMKnnXrpp4aUm/f8WO6FheDgD5q9igM2LEiFi2bFksXLiw1E0BAA6Rig06AED+BB0AIFuCDgCQLUEHAMiWoAMAZEvQAQCyJegAANkSdACAbFVs0LEzMgDkr2KDjp2RASB/FRt0AID8CToAQLYEHQAgW4IOAJAtQQcAyJagAwBkS9ABALIl6AAA2RJ0AIBsCToAQLYqNui41xUA5K9ig457XQFA/io26AAA+RN0AIBsCToAQLYEHQAgW4IOAJAtQQcAyJagAwBkS9ABALIl6AAA2RJ0AIBsCToAQLYEHQAgW4IOAJAtQQcAyJagAwBkK5ug89Zbb8Upp5wSN9xwQ6mbAgCUiWyCzm233RYXXHBBqZsBAJSRLILOK6+8EsuXL48hQ4aUuikAQBkpedCZN29eDB06NDp16hRVVVUxc+bMPV4zZcqU6NatW7Ru3Tr69+8fCxYsqPd8mq6aMGHCYWw1AHAkKHnQ2bZtW/Ts2bMIM3szY8aMGD16dIwfPz6WLFlSvHbw4MGxfv364vkHH3wwzjjjjOI4ENu3b4/NmzfXOwCAPLUodQPSdNP+ppwmTZoU1113XQwfPrw4nzp1asyaNSumTZsWY8aMiWeeeSbuu+++uP/++2Pr1q2xc+fOaNOmTYwbN26v75dGfm655ZZD9nkAgPJR8hGd/dmxY0csXrw4Bg0aVHetWbNmxfn8+fPrgsvq1atj1apV8b3vfa8IRfsKOcnYsWNj06ZNdUf6WgAgTyUf0dmfjRs3xu7du6NDhw71rqfzVHzcEK1atSoOACB/ZR10DtY111xT6iYAAGWkrKeu2rdvH82bN49169bVu57OO3bs2Kj3TsXPPXr0iL59+zaylQBAuSrroNOyZcvo3bt3zJkzp+5adXV1cT5gwIBGvfeIESNi2bJlsXDhwiZoKQBQjko+dZVWSq1YsaLufOXKlbF06dJo165ddO3atVhaPmzYsOjTp0/069cvJk+eXCxJr12FBQBQtkFn0aJFMXDgwLrzFGySFG6mT58eV111VWzYsKFYSbV27dro1atXzJ49e48C5YZMXaUjFTsDAHmqqqmpqYkKljYMbNu2bbHUPO2/AwC56DZmVqmbEKsmXlrS399lXaMDANAYgg4AkC1BBwDIVsUGHfvoAED+Kjbo2EcHAPJXsUEHAMifoAMAZKtig44aHQDIX8UGHTU6AJC/ig06AED+BB0AIFuCDgCQrZLfvRwAyukmlIfyRpQcfhU7omPVFQDkr2KDjlVXAJC/ig06AED+BB0AIFuCDgCQLUEHAMiWoAMAZKtig47l5QCQv4oNOpaXA0D+KjboAAD5E3QAgGwJOgBAtgQdACBbgg4AkC1BBwDIlqADAGSrYoOODQMBIH8VG3RsGAgA+avYoAMA5E/QAQCyJegAANkSdACAbAk6AEC2BB0AIFuCDgCQLUEHAMiWoAMAZEvQAQCyJegAANmq2KDjpp4AkL+KDTpu6gkA+avYoAMA5E/QAQCyJegAANkSdACAbAk6AEC2BB0AIFuCDgCQLUEHAMiWoAMAZEvQAQCyJegAANkSdACAbAk6AEC2BB0AIFuCDgCQrSM+6Lz55pvRp0+f6NWrV5xzzjlx9913l7pJAECZaBFHuOOOOy7mzZsXxxxzTGzbtq0IO1deeWWceOKJpW4aAFBiR/yITvPmzYuQk2zfvj1qamqKAwCg5EEnjcYMHTo0OnXqFFVVVTFz5sw9XjNlypTo1q1btG7dOvr37x8LFizYY/qqZ8+e0blz57jxxhujffv2h/ETAADlquRBJ003pZCSwszezJgxI0aPHh3jx4+PJUuWFK8dPHhwrF+/vu41xx9/fDz//POxcuXKuPfee2PdunX7/H5p1Gfz5s31DgAgTyUPOkOGDIlbb701rrjiir0+P2nSpLjuuuti+PDh0aNHj5g6dWoxVTVt2rQ9XtuhQ4ciCD355JP7/H4TJkyItm3b1h1dunRp0s8DAJSPkged/dmxY0csXrw4Bg0aVHetWbNmxfn8+fOL8zR6s2XLluLxpk2biqmw7t277/M9x44dW7yu9li9evVh+CQAQCmU9aqrjRs3xu7du4uRmndK58uXLy8e/8///E986UtfqitC/spXvhLnnnvuPt+zVatWxQEA5K+sg86B6NevXyxduvSgvy7VBKUjBSkAIE9lPXWVVk+l5ePvLi5O5x07dmzUe48YMSKWLVsWCxcubGQrAYByVdZBp2XLltG7d++YM2dO3bXq6urifMCAASVtGwBQ/ko+dbV169ZYsWJF3XlaIp6motq1axddu3YtlpYPGzasuM1DmqaaPHlysSQ9rcICACjroLNo0aIYOHBg3XkKNkkKN9OnT4+rrroqNmzYEOPGjYu1a9cW97SaPXv2HgXKB0uNDgDkr6qmwu+XkDYMTPvppKXmbdq0KXVzACpatzGzohysmnhp5KBbGfTnoerLA/39XdY1OgAAjSHoAADZqtigk2p00i0l+vbtW+qmAACHSMUGHfvoAED+Sr7qCuBIl3PBJxzpKnZEBwDIn6ADAGSrYoOOYmQAyF/FBh3FyACQv4oNOgBA/gQdACBbgg4AkC1BBwDIVsUGHauuACB/FRt0rLoCgPxVbNABAPIn6AAA2RJ0AIBsCToAQLYEHQAgWxUbdCwvB4D8VWzQsbwcAPJXsUEHAMhfg4LOBz7wgXj99df3uP7mm28WzwEAHLFBZ9WqVbF79+49rm/fvj1effXVpmgXAECjtTiYFz/00EN1j3/2s59F27Zt685T8JkzZ05069at8a0CADjcQefyyy8v/qyqqophw4bVe+6oo44qQs4dd9zRFO0CADi8Qae6urr489RTT42FCxdG+/btG98CAIByCDq1Vq5c2fQtAQAoh6CTpHqcdKxfv75upKfWtGnT4kjYMDAdeyuqBgAqeNXVLbfcEpdcckkRdDZu3BhvvPFGveNIYMNAAMhfg0Z0pk6dGtOnT4+rr7666VsEAFDKEZ0dO3bEhRde2FRtAAAon6DzxS9+Me69996mbw0AQKmnrt5+++2466674vHHH4/zzjuv2EPnnSZNmtRU7QMAOLxB54UXXohevXoVj1988cV6z6XNBAEAjtigM3fu3KZvCQBAOdToAABkO6IzcODA/U5RPfHEE41pEwBA6YJObX1OrZ07d8bSpUuLep133+wTAOCICjrf//7393r9W9/6VmzdurWxbQIAKL8anc9//vNHxH2uAIDK0KRBZ/78+dG6des4EqQbevbo0SP69u1b6qYAAOU0dXXllVfWO6+pqYnXXnstFi1aFDfffHMcKTf1TMfmzZujbdu2pW4OAFAuQefdwaBZs2bRvXv3+Pa3v13c1RwA4IgNOvfcc0/TtwQAoByCTq3FixfHSy+9VDw+++yz4/zzz2+qdgEAlCborF+/Pj772c/Gz3/+8zj++OOLa2+++WaxkeB9990X73//+xvfMgCAUqy6+spXvhJbtmyJX/3qV/G73/2uONJmgamw9/rrr29smwAASjeiM3v27Hj88cfjrLPOqruWlmqnJduKkQGAI3pEp7q6Oo466qg9rqdr6TkAgCM26Hz84x+PUaNGxZo1a+quvfrqq/G1r30tLr744qZsHwDA4Q06P/zhD4t6nG7dusVpp51WHKeeempx7c4772x4awAASl2j06VLl1iyZElRp7N8+fLiWqrXGTRoUFO2DQDg8I3oPPHEE0XRcRq5qaqqik984hPFCqx0pHtGpb10nnzyyca1CACgFEFn8uTJcd1110WbNm32eluIL3/5yzFp0qSmahsAwOELOs8//3z8wR/8wT6fT0vL027JAABHXNBZt27dXpeV12rRokVs2LChKdoFAHB4g87JJ59c7IC8Ly+88EKcdNJJcTitXr06LrrooqJ26Lzzzov777//sH5/ACCToPPJT34ybr755nj77bf3eO5///d/Y/z48XHZZZfF4ZRGkVLt0LJly+LRRx+Nr371q7Ft27bD2gYAIIPl5d/85jfjJz/5SZxxxhkxcuTI6N69e3E9LTFPt3/YvXt33HTTTXE4pRGk2lGkjh07Rvv27Yt7bx177LGHtR0AwBE+otOhQ4f4z//8zzjnnHNi7NixccUVVxTHN77xjeLaU089VbzmYMybNy+GDh0anTp1Kpasz5w5c4/XpBCVNids3bp19O/fPxYsWLDX90qF0ClspX1+AAAOemfkU045Jf7jP/4jNm7cGM8++2w888wzxeN0Le2OfLDSNFPPnj2LMLM3M2bMiNGjRxfTYmmLine truncated
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.yscale('log')\n",
+    "digit_lengths(PP)"
+   ]
+  },
  {
   "cell_type": "markdown",
   "id": "f2323985-2518-41d1-8a3b-6e21e40d27be",
   "metadata": {},
   "source": [
-    "# Summary of Truncatable Prime Facts\n",
+    "# Summary \n",
    "\n",
    "Here's what we learned:\n",
    "\n",
@@ -474,6 +653,9 @@
    "- Palindromic truncatable primes:\n",
    "  - Count: 13\n",
    "  - Largest: 799636997 (9 digits)\n",
+    "- Palindromic primes:\n",
+    "  - Count: infinite\n",
+    "  - All except \"11\" have an odd number of digits\n",
    "\n",
    "\n",
    "# Primality Testing\n",
@@ -483,7 +665,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 19,
   "id": "1d460931-3d70-485d-807f-8fb360ef4ef2",
   "metadata": {},
   "outputs": [],
@@ -507,7 +689,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 20,
   "id": "647fb66b-32ae-494a-8945-cf7358c2f958",
   "metadata": {},
   "outputs": [],
@@ -559,7 +741,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 21,
   "id": "c8a59c5e-a01c-49be-a931-c1547c4ac712",
   "metadata": {},
   "outputs": [],
@@ -582,7 +764,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 22,
   "id": "476282dd-a48a-4eb5-8fed-8495d4878ec2",
   "metadata": {},
   "outputs": [],
@@ -615,7 +797,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 23,
   "id": "45a11a93-9ef3-48d1-80fa-e0113041346e",
   "metadata": {},
   "outputs": [],
@@ -638,26 +820,26 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 24,
   "id": "9eba96df-c298-4d71-b68a-3bf36bf15095",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
-       "{3: 1.433229785995569,\n",
-       " 4: 0.47277977521728054,\n",
-       " 5: 0.12450624459417024,\n",
-       " 6: 0.04328208013677137,\n",
-       " 7: 0.010154128989503837,\n",
-       " 8: 0.0021231309812580613,\n",
-       " 9: 0.0005263573440008422,\n",
-       " 10: 0,\n",
+       "{3: 1.4119561859956755,\n",
+       " 4: 0.44668416368593794,\n",
+       " 5: 0.12196467991169978,\n",
+       " 6: 0.028730958963517104,\n",
+       " 7: 0.008030451471981755,\n",
+       " 8: 0.0026523932544334754,\n",
+       " 9: 0.0005267205326198026,\n",
+       " 10: 0.0005232807610595388,\n",
       " 11: 0,\n",
       " 12: 0}"
      ]
     },
-     "execution_count": 19,
+     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
@@ -678,7 +860,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 25,
   "id": "8880c588-32f5-46e8-b59c-f9b49d492d3b",
   "metadata": {},
   "outputs": [
@@ -688,8 +870,8 @@
       "{3: 0,\n",
       " 4: 0,\n",
       " 5: 0,\n",
-       " 6: 0.0003613604498214879,\n",
-       " 7: 0.000356548185704557,\n",
+       " 6: 0.0003605006633212205,\n",
+       " 7: 0,\n",
       " 8: 0,\n",
       " 9: 0,\n",
       " 10: 0,\n",
@@ -697,7 +879,7 @@
       " 12: 0}"
      ]
     },
-     "execution_count": 20,
+     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
@@ -720,35 +902,35 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 26,
   "id": "7871aaa7-af1c-462b-85bb-d46818ccabea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
-       "{561: 56.84,\n",
-       " 1105: 69.72,\n",
-       " 1729: 74.68,\n",
-       " 2465: 72.79,\n",
-       " 2821: 76.78,\n",
-       " 6601: 80.1,\n",
-       " 8911: 80.1,\n",
-       " 41041: 70.28,\n",
-       " 101101: 71.63,\n",
-       " 294409: 95.18,\n",
-       " 340561: 81.13,\n",
-       " 825265: 59.82,\n",
-       " 56052361: 99.19,\n",
-       " 118901521: 99.24,\n",
-       " 172947529: 99.45,\n",
+       "{561: 56.85,\n",
+       " 1105: 68.86,\n",
+       " 1729: 76.46,\n",
+       " 2465: 72.11,\n",
+       " 2821: 76.03,\n",
+       " 6601: 80.03,\n",
+       " 8911: 79.98,\n",
+       " 41041: 70.17,\n",
+       " 101101: 70.69,\n",
+       " 294409: 95.11,\n",
+       " 340561: 81.57,\n",
+       " 825265: 60.17,\n",
+       " 56052361: 99.21,\n",
+       " 118901521: 99.29,\n",
+       " 172947529: 99.43,\n",
       " 216821881: 99.4,\n",
-       " 1299963601: 99.65,\n",
-       " 2301745249: 99.76,\n",
-       " 9624742921: 99.9}"
+       " 1299963601: 99.64,\n",
+       " 2301745249: 99.72,\n",
+       " 9624742921: 99.84}"
      ]
     },
-     "execution_count": 21,
+     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
@@ -770,7 +952,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 27,
   "id": "de492047-77b6-46a3-8a77-46308f219efd",
   "metadata": {},
   "outputs": [
@@ -786,19 +968,19 @@
       " 8911: 0,\n",
       " 41041: 0,\n",
       " 101101: 0,\n",
-       " 294409: 4,\n",
+       " 294409: 6,\n",
       " 340561: 0,\n",
       " 825265: 0,\n",
-       " 56052361: 58,\n",
-       " 118901521: 66,\n",
-       " 172947529: 80,\n",
-       " 216821881: 86,\n",
-       " 1299963601: 90,\n",
-       " 2301745249: 87,\n",
-       " 9624742921: 94}"
+       " 56052361: 62,\n",
+       " 118901521: 81,\n",
+       " 172947529: 75,\n",
+       " 216821881: 71,\n",
+       " 1299963601: 91,\n",
+       " 2301745249: 94,\n",
+       " 9624742921: 92}"
      ]
     },
-     "execution_count": 22,
+     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
@@ -823,7 +1005,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 28,
   "id": "be5f2fe7-b1e0-4a02-b7f0-7903e865c4b8",
   "metadata": {},
   "outputs": [],
@@ -845,7 +1027,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": 29,
   "id": "8109d9fc-bbae-4ef0-a75c-908d9746a557",
   "metadata": {},
   "outputs": [],
@@ -884,13 +1066,13 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 30,
   "id": "b1ab0f3d-6ce6-4b8c-9d4d-bb78bb1ac432",
   "metadata": {},
   "outputs": [
    {
     "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAINCAYAAACTVPhqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAwGpJREFUeJzs3Qd8U1X7B/BfuncpnbTsIVL2liUyFJQhQ/EVF+6BrwiC4MTXAbgQB/5Bfd0TEX2VJYjI3hsqsjedQPdM8v88J01oS0fSJs36ff1cc3Jzc3PuTWnvc89zztHo9Xo9iIiIiIiIyuFR3koiIiIiIiIGDEREREREVCm2MBARERERUYUYMBARERERUYUYMBARERERUYUYMBARERERUYUYMBARERERUYUYMBARERERUYW8Kn6JhE6nw7lz5xAcHAyNRsOTQkREREROT+ZuzszMRGxsLDw8Km9DYMBQBQkWGjRoYM3vh4iIiIjIIZw+fRr169evdBsGDFWQlgXjyQwJCbHet0NEREREZCcZGRnqprjxWrcyDBiqYExDkmCBAQMRERERuRJzUu7Z6ZmIiIiIiCrEgIGIiIiIiCrEgIGIiIiIiCrEPgxEREREdhrWsqioCFqtluefbMLb2xuenp413g8DhgrMnTtXLfxHTERERNZWUFCA8+fPIycnhyeXbNqhWYZMDQoKqtl+9BLeUqVDToWGhiI9PZ2jJBEREZFVJoU9fPiwuvMbGRkJHx8fTg5LVieX+CkpKSoobdGixRUtDZZc47KFgYiIiKiWWxckaJAx8AMCAnjuyWYkID1x4gQKCwtrlJrETs9EREREduDhwcswsv8cC+bgTyoREREREVWIAQMREREREVWIAQMRERGRk9Lq9Nh0NA3/231WPcpzW7ruuuvw5JNP2mz/km8vaTS7d++GIzpRi/Vr3Lgx5syZA0fATs9ERERETmj5/vP4z28JOJ+eZ1pXL9QP04fFY3Cbejb5zEWLFqmx/W1FOoLLcLMRERFwRA0cvH62whYGIiIiIicMFh79emepYEEkpuep9fK6LdStWxfBwcE2Gz1KRvKJiYmBl5dj3tP2dPD62QoDBnI6td38SkREVBtj5ucUFJm1ZOYVYvqvB1DeXz/jupd+TVDbmbM/S6bkKpmS9OGHH6rx/f38/BAdHY1bbrml1HaPP/64WmSsf7kj/8ILL5T6LEm5eeWVV3D33XereQAeeuihK1J+/vrrL/X8999/R8eOHeHv74/+/fsjOTkZy5YtQ6tWrdR7x44dW2oSPBm2dubMmWjSpIl6T/v27bFw4UKzjvHixYu444471JCk8t4WLVrgs88+U69Zq37mnJ+yLl26hAceeEDVS/Ypn7Nnzx7UBvcKj8jp2aP5lYiIyNZyC7WIf/F3q+xLLjkTM/LQ9qUVZm2f8PIgBPhYdkm4fft2PPHEE/jqq6/Qs2dPXLhwAevWrSu1zRdffIH7778fW7duVdtLQNCwYUM8+OCDpm3eeustvPjii5g+fXqln/fSSy/hgw8+UPNWjBkzRi2+vr749ttvkZWVhZEjR+L999/H1KlT1fYSLHz99deYN2+euuBfu3Yt7rzzTnWx3bdv30o/Sy7cExIS1AW/XMgfOXIEubm5Vq2fueenpFtvvVUFJFIvCTLmz5+PAQMG4NChQ6rlx5YYMJDTNb+Wjb2Nza//d2cnBg1ERES14NSpUwgMDMTQoUNVilKjRo3UHfay+f7vvPOOugPfsmVL7Nu3Tz0veUEsd8mfeuop03O5g1+eV199Fb169VJluch+5plncPToUTRt2lStk9aN1atXqwvy/Px8zJgxA3/88Qd69OihXpft1q9fry6yqwoY5NjkWLp06WJqCamKJfWz5PwYSd0lsJCWCwlEjMHWL7/8olpOJNiwJQYMFZg7d65atFqtTb8AMo+kHUnLQkXNrzItibx+fXwMPD2sM0kJkcvTaYGTG4GsJCAoGmjUE/Co/kygRFR9/t6e6k6/ObYev4Bxn22rcrvP7+2Kbk3qmvXZlrr++utVkCAXxIMHD1aL3EUvOXP1NddcU2riMLl4f/vtt9W1lXHWYeNFeVXatWtnKkv6k3yO8WLcuE4uqIW0CEj6j9SxbB+JskFNeR599FGMHj0aO3fuxA033IARI0aoVhRr1c+S82MkqUfSUhEeHl5qvbR8SGBiawwYKjB+/Hi1ZGRkqGYfsi/55Vi2Y1fZoEFel+16NCv9j4mIypHwK7B8KpBx7vK6kFhg8OtA/HCeMqJaJheO5qYF9WkRqdJxpYW9vBtpcgkaE+qntrPVTTRpVZALasnhX7FihUorkrScbdu2oU6dOmbvR1opzFFyZCY5V2VHapJ10m9ByIW1WLJkCeLi4kptZ7w7X5kbb7wRJ0+exNKlS7Fy5UqV9iPXhHJH3xr1qw45pnr16qnzXZYl57u62OmZnEJyZp5VtyNyaxIsLLi7dLAgMs4b1svrROSwJAiQvnuibDhgfC6v27rFXUYKGjhwIN544w3s3btXpRP9+eefpte3bNlSavvNmzer/gRl755bW3x8vAoMJLWoefPmpRZJAzKH9HW45557VD+IOXPm4KOPPrJ6PS05P506dUJiYqI652WPqTaGeGULAzmFqGA/q25H5NZpSNKyUFmC3/JpwNVDmJ5E5MBkoA/pu1d2IJCYWhoIZPHixTh27BiuvfZahIWFqbvxcgddcvGN5IJ90qRJePjhh1VrhHT6lZQbW5PWj8mTJ2PixImqTr1790Z6ejo2bNigRheSQKAy0lrSuXNntG7dWvWHWLx4sRrtyNosOT8SmEnKkqRHSYB21VVX4dy5c6oVRVLBzE3tqi4GDOQUJAdTml8rSksyNr+ak6tJ5Nakz0LZloVS9EDGWcN2TfrUYsWIyFISFEjfPUnHlRZ2uWkmfwdroy+fpMHIJG6ShpSXl6fujH/33XfqIttIhkuVHPtu3bqpu+YTJkyweedcIxmuVVoJZLQkCWykvnKX/tlnn63yvT4+PqrTsrSYyKhEffr0wffff2/1OlpyfiSlSYKy5557Dvfeey9SUlLUfBASsEn/CFvT6C0ZfNcNGfswSGQqUSnZd5SkR77eecV6469FjpJEZIZ9C4Gf7q96u9H/BdpeHlOdiKxHLrCPHz+u5giQOQxckcwz0KFDB5XOQ/Y7P5X9rFlyjcs+DOQ0ereIhI/XlT+y0rLAYIHITDIakjW3IyIil8eUJHIaS/eeR0GRDk3CAzBjVFskZ+bXavMrkUtocA3gHQAUXp5xtDSNYbQkGWKViMgFPfLII6ozc3lkcjeZ7I1KY8BATuPHHafV461dG6BHM9uPCEDkcrRFwG//rjxYEINnscMzEdVIecN/OoqXX35ZdYouT22ln//lwOenPAwYyCkcS8nCthMXIQ0JozvVt3d1iJyPthBY9CBw4GdA4wl0fxhI+KWceRhmcR4GInJpUVFRaiHzMWAgp7Bwxxn1eO1VkYgOcc0OYkQ2U5gH/DgOOLQM8PAGbv0MaDUMuOFVzvRMRERVYsBADk+r02PRzrOqfGtn8yZcIaJiBTnA92OBY6sBLz/gtm+AFgMNr3l4cuhUIiKqEgMGcnjrDqcgMSMPdQK8MTCeTYhEZsvPBL69DTi5AfAOBMZ+DzS5lieQiIgswoCBHN6P2w3pSCM6xMHXy7bTyRO5jNyLwNe3AGe3A74hwB0LgYbd7V0rIiJyQm4xD4NM6S1TlcsshJ988om9q0MWuJRTgJUJSap8S2d2diYyS3Yq8MUwQ7DgHwbc8yuDBSIiqjaXDxiKioowadIk/Pnnn9i1axfefPNNpKWl2btaZKb/7T6HAq0O8fVC0CYulOeNqCqZicDnQ4DEfUBgFDBuCRDbkeeNyFXptMDxdYZZ3OVRntt4huInn3zSZvs/ceIENBoNdu/eDXvasGED2rZtC29vb4wYMQLuzuVTkrZu3YrWrVsjLi5OPb/xxhuxYsUK3H777fauGlky90IXti4QVenSaeDL4cCFY0BwrKFlIaIFTxyRq0r4FVg+tZzhkV+32fDIixYtUhfRttKgQQOcP38eERH2nW9JbjZ36NABy5YtQ1BQUI3399JLL+GXX36xeyDksi0Ma9euxbBhwxAbG6siTjnZZc2dOxeNGzeGn58funfvroIEo3PnzpmCBSHls2cNI+6QY0s4l4H9ZzPg7anBzR0uf4dEVA4JEj67yfBYpyFw3zIGC0SuHiwsuLt0sCAyzhvWy+s2ULduXQQHB9tk3wUFBfD09ERMTAy8vOx7T/vo0aPo378/6tevjzp16sBRyDmyB4cPGLKzs9G+fXsVFJTnhx9+UFHg9OnTsXPnTrXtoEGDkJycXOt1Jdu0LlwfH426gT48vUQVSTlkCBbSTwHhzYF7lwNhjXm+iJyJXg8UZJu35GUAy56WN5W3I8ODtDzIdubsTz67GilJH374oeofKjdso6Ojccstt5Ta7vHHH1dLaGioajF44YUXoC/xWXKz95VXXsHdd9+tZlh+6KGHrkhJkhmR5fnvv/+Ojh07wt/fX13Iy3We3P1v1aqVeu/YsWORk3N5FnudToeZM2eiSZMm6j1yfbhw4cIqj8/4+ZK+ft9996ny559/Dq1Wi/vvv9+0P+kb++6775Z6r9S1W7duCAwMVEFGr169cPLkSfX+//znP9izZ4/an3Gf4tKlS3jggQcQGRmpjkOOTbYr2TIhLR3SB1c+W861PTh8SpKkEMlSkdmzZ+PBBx/Evffeq57PmzcPS5Yswaeffopp06aplomSLQpSli+zIvn5+WoxysjIUI+FhYWmPhFCIl9ZJ1962bIxQpalbFle9/DwUJ8hTXrllX18DBfHsn3Jsq+vr/oHIJ9VXlnqJtvLD7UsZcvyuvxDlc8qW3a0Y9J4euOXXZfnXnCFY3LF74nH5ADf04V/gK9GQpOTCkTFQ3fnIhT61oVv8R9Mfk/898TfEY73e8/Line truncated
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAINCAYAAACTVPhqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAv9BJREFUeJzs3Qd8U1X7B/BfuncpdFH2EmRvBBRkyBIQVHALoqiv+EdBBcerOAEX4sAXFXHgAhEcLBkiQ/YUqMjeHbSleyf5f54TEtrSkbRJs35fP9ec3NzenHvThvvc85xzNHq9Xg8iIiIiIqJSeJS2koiIiIiIiAEDERERERGViy0MRERERERUJgYMRERERERUJgYMRERERERUJgYMRERERERUJgYMRERERERUJgYMRERERERUJq+yXyKh0+lw4cIFBAcHQ6PR8KQQERERkdOTuZszMjIQExMDD4/y2xAYMFRAgoV69epZ8/MhIiIiInIIZ8+eRd26dcvdhgFDBaRlwXgyQ0JCrPfpEBERERHZSXp6uropbrzWLQ8DhgoY05AkWGDAQERERESuxJyUe3Z6JiIiIiKiMjFgICIiIiKiMjFgICIiIiKiMrEPAxEREZGdhrUsLCyEVqvl+Seb8Pb2hqenZ5X3w4CBiIiIqJrl5+cjLi4O2dnZPPdk0w7NMmRqUFBQlfbDgKEMc+bMUQujfiIiIrL2pLAnT55Ud35l0iwfHx9ODks2acG6ePEizp07h2bNmlWppUGjl71RuWPUhoaGIi0tjcOqEhERUZXl5uaqgKFBgwYICAjgGSWbycnJwalTp9CoUSP4+flV+hqXnZ6JiIiI7MDDg5dhZP85FszB31QiIiIiIioTAwYiIiIiIioTAwYiIiIiJ6XV6bH1eDJ+2XdePcpzW7rxxhvx5JNP2mz/km8vaTT79u2DIzpVjfVr2LAhZs+eDUfAUZKIiIiInNCqg3F45bdYxKXlmtbVDvXDtGEtMah1bZu855IlS9TY/rZSr149NdxseHg4HFE9B6+frbCFgYiIiMgJg4X/fLOnWLAg4tNy1Xp53RZq1qyJ4OBgm81NIUN/RkdHw8vLMe9pezp4/WyFAQM5nepufiUiIrI1GeU+O7/QrCUjtwDTfj2E0v71M657+ddYtZ05+7NkhP2iKUkff/yxGt9fhuuMiorC7bffXmy7xx9/XC0ydKfckX/xxReLvZek3Lz22mu4//771bCeDz/88FUpP3/++ad6/vvvv6NDhw7w9/dH3759kZiYiJUrV+Laa69VP3v33XcXmwRP5rqYMWOGGk5UfqZdu3ZYvHixWcd46dIl3HPPPYiIiFA/26xZM3zxxRfqNWvVz5zzU1JqaioeeughVS/Zp7zP/v37UR3cKzwip2eP5lciIiJbyynQouVLv1tlX3LJGZ+eizYvrzZr+9hXByLAx7JLwl27dmHixIlYsGABevTogZSUFGzatKnYNl999RUefPBB7NixQ20vAUH9+vUxfvx40zbvvPMOXnrpJUybNq3c93v55Zfx0UcfqXkrRo8erRZfX1989913yMzMxMiRI/Hhhx9i6tSpansJFr755hvMnTtXXfBv3LgR9957r7rY7t27d7nvJRfusbGx6oJfLuSPHTum5jOwZv3MPT9FjRo1SgUkUi8JMj755BP069cPR44cUS0/tsSAgZyu+bVk7G1sfv3fvR0ZNBAREVWDM2fOIDAwEEOHDlUpSjIJndxhL5nv/95776k78M2bN8eBAwfU86IXxHKX/KmnnjI9lzv4pXn99dfRs2dPVZaL7Oeeew7Hjx9H48aN1Tpp3Vi/fr26IM/Ly8P06dOxdu1adO/eXb0u223evFldZFcUMMixybF07tzZ1BJSEUvqZ8n5MZK6S2AhLRcSiBiDrZ9//lm1nEiwYUsMGMowZ84ctWi1Wpt+AGQeSTuSloWyml9lWhJ5/aaW0fD0sM4kJUQuT6cFTm8BMhOAoCigQQ/Aw9PetSJyS/7enupOvzl2nEzB2C92Vrjdlw90QddGNc16b0vddNNNKkiQC+JBgwapRe6iF525+rrrris2cZhcvL/77rvq2kr6AgjjRXlF2rZtaypL+pO8j/Fi3LhOLqiFtAhI+o/UsWQfiZJBTWn+85//4LbbbsOePXswYMAAjBgxQrWiWKt+lpwfI0k9kpaKWrVqFVsvLR8SmNgaA4YyTJgwQS3GabPJvuTLsWTHrpJBg7wu23VvUvyPiYhKEfsrsGoqkH7hyrqQGGDQm0DL4TxlRNVMLhzNTQu6oVmESseVFvbSbqTJJWh0qJ/azlY30aRVQS6oJYd/9erVKq1I0nJ27tyJGjVqmL0faaUwR9GRmeRclRypSdZJvwUhF9Zi+fLlqFOnTrHtjHfnyzN48GCcPn0aK1aswJo1a1Taj1wTyh19a9SvMuSYateurc53SZac78pip2dyCokZuVbdjsitSbCw6P7iwYJIjzOsl9eJyGFJECB990TJcMD4XF63dYu7jBTUv39/vPXWW/j7779VOtEff/xhen379u3Ftt+2bZvqT1Dy7rm1tWzZUgUGklrUtGnTYoukAZlD+jqMGTNG9YOYPXs2Pv30U6vX05Lz07FjR8THx6tzXvKYqmOIV7YwkFOIDPaz6nZEbp2GJC0L5SX4rXoWaHEz05OIHJgM9CF990oOBBJdTQOBLFu2DCdOnECvXr0QFham7sbLHXTJxTeSC/bJkyfjkUceUa0R0ulXUm5sTVo/nn76aUyaNEnV6frrr0daWhr++usvNbqQBALlkdaSTp06oVWrVqo/xLJly9RoR9ZmyfmRwExSliQ9SgK0a665BhcuXFCtKJIKZm5qV2UxYCCnIDmY0vxaVlqSsfnVnFxNIrcmfRZKtiwUowfSzxu2a3RDNVaMiCwlQYH03ZN0XGlhl5tm8u9gdfTlkzQYmcRN0pByc3PVnfHvv/9eXWQbyXCpkmPftWtXddf8iSeesHnnXCMZrlVaCWS0JAlspL5yl/7555+v8Gd9fHxUp2VpMZFRiW644Qb88MMPVq+jJedHUpokKHvhhRfwwAMP4OLFi2o+CAnYpH+ErWn0lgy+64aMfRgkMpWolOw7StKj3+y5ar3xa5GjJBGZ4cBi4KcHK97uts+BNlfGVCci65EL7JMnT6o5AmQOA1ck8wy0b99epfOQ/c5Peb9rllzjsg8DOQ3pvOXjdfWvrLQsMFggMpOMhmTN7YiIyOUxJYmcxvIDccgv1KFRrQBMv7UNEjPyqrX5lcgl1O8O+AQC+VllbKAxjJYkQ6wSEbmgRx99VHVmLo1M7iaTvVFxDBjIafy466x6vL1zPXRvYvsRAYhcjgzpt+Kp8oMFMWgmOzwTUZWUNvyno3j11VdVp+jSVFf6+Z8OfH5Kw4CBnMLJpCzsPHUJ0pBwW8e69q4OkfPRFgK/TAD+/gHQeACdHwT+XV7KPAwzOQ8DEbm0yMhItZD5GDCQU1i829C60OuaCNVngYgsUJgPLHkIiP0F0HgCt35q6NA8+E3O9ExERBViwEAOT6vT46fd51V5VCfzJlwhossKcoEfxwBHVgGePsCoLw1zLAgPTw6dSkREFWLAQA5v09GLiE/PRY0Ab/RvySZEIrNJX4Uf7gZO/Al4+QF3fgs07c8TSEREFmHAQA7vx93n1OMt7WLg62Xb6eSJXEZuOvDdaODMVsA7ELh7IVsTiIioUtxiHgaZ0lumKpdZCOfNm2fv6pAFUrPzseZQgiqP6sx0JCKzZKcAX99iCBZ8Q4H7f2GwQERElebyAUNhYSEmT56MP/74A3v37sXbb7+N5ORke1eLzPTr/gvI1+pwbe0QtK4TyvNGVJHMi8BXw4ELewD/msCYX4F6XXjeiFyVTguc3GSYxV0e5bmNZyh+8sknbbb/U6dOQaPRYN++fbCnv/76C23atIG3tzdGjBgBd+fyKUk7duxAq1atUKdOHfV88ODBWL16Ne666y57V43MsOjy3AujOnEoVaIKpccBXw8Hko4AgZGGloWoljxxRK4q9ldg1dRShkd+02bDIy9ZskRdRNtKvXr1EBcXh/Bw+863JDeb27dvj5UrVyIoKKjK+3v55Zfx888/2z0QctkWho0bN2LYsGGIiYlREaec7JLmzJmDhg0bws/PD926dVNBgtGFCxdMwYKQ8vnzhhF3yLH9E5eOg+fT4e2pwYgOVz5DIipF6hngi8GGYCGkDvDASgYLRK4eLCy6v3iwYLxxIOvldRuoWbMmgoODbbLv/Px8eHp6Ijo6Gl5e9r2nffz4cfTt2xd169ZFjRo14CjkHNmDwwcMWVlZaNeunQoKSrNw4UIVBU6bNg179uxR2w4cOBCJiYnVXleyrh93GTo79782CjUDfXh6icqSfByYPxi4dBKo0cAQLIQ35fkiciZ6vWFkM3MWGdRg5RT5odJ2ZHiQlgfZzpz9yXtXIiXp448/Vv1D5YZtVFQUbr/99mLbPf7442oJDQ1VLQYvvvgi9EXeS272vvbaa7j//vvVDMsPP/zwVSlJMiOyPP/999/RoUMH+Pv7qwt5uc6Tu//XXnut+tm7774b2dnZpn3rdDrMmDEDjRo1Uj8j14eLFy+u8PiM7y/p6+PGjVPlL7/8ElqtFg8++KBpf9I39v333y/2s1LXrl27IjAwUAUZPXv2xOnTp9XPv/LKK9i/f7/an3GfIjU1FQ899BAiIiLUccixyXZFWyakpUP64Mp7y7m2B4dPSZIUIlnKMmvWLIwfPx4PPPCAej537lwsX74c8+fPx7PPPqtaJoq2KEhZPsyy5OXlqcUoPT1dPRYUFJj6RAiJfGWdfOgly8YIWZaSZXndw8NDvYc06ZVW9vExXBzL9kXLvr6+6g9A3qu0stRNtpdfallKluV1+UOV9ypZdrRj0nh64+e9hoBhVOe6LnFMrvg58Zgc4HNKOwF8PQKazHigVjPo7vsZBf4R8L38DyY/J/498TvLine truncated
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
@@ -929,7 +1111,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 26,
+   "execution_count": 31,
   "id": "4f377075-e352-46f9-a716-f1de673a7086",
   "metadata": {},
   "outputs": [],
@@ -952,12 +1134,18 @@
    "\n",
    "Just one more thing: none of this would work unless we can efficiently compute *a*<sup>(*n* - 1)</sup> (mod *n*).  How does the `pow` builtin function do it? When *a* and *n* are 24-digit numbers, if we naively tried to compute `a ** (n - 1)`, we'd have two problems: we'd need nearly a billion petabytes to store the result, and we'd need centuries to compute it. The way around these problems is to use [modular exponentiation](https://en.wikipedia.org/wiki/Modular_exponentiation) where we apply the modulus to each intermediate result, and cut the exponent in half each iteration, so we need only do *O*(log *n*) multiplications, not *O*(*n*). That's a big difference: 10<sup>24</sup> is a trillion trillion, and log<sub>2</sub>(10<sup>24</sup>) is only 80. \n",
    "\n",
-    "The key idea is that *b*<sup>2*e*</sup> is equal to (*b* × *b*)<sup>*e*</sup>. Here is an implementation:"
+    "There are two key ides that make this work:\n",
+    "1) *b*<sup>2*e*</sup> = (*b* × *b*)<sup>*e*</sup>\n",
+    "2) *b*<sup>2</sup> (mod *m*) = (*b* (mod *m*))<sup>2</sup> (mod *m*)\n",
+    "    - Proof: let *r* = *b* (mod *m*); then *b* = (*km* + *r*) for some integer value of *k*, and:\n",
+    "    - *b*<sup>2</sup> (mod *m*) = (*km* + *r*)<sup>2</sup> (mod *m*) = ((*km*)<sup>2</sup> + 2*kmr*  + *r*<sup>2</sup>) (mod *m*)  = *r*<sup>2</sup> (mod *m*)\n",
+    "\n",
+    "Here is an implementation:"
   ]
  },
  {
   "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 32,
   "id": "e71eb140-7e1a-4146-87a7-8c4eebf61486",
   "metadata": {},
   "outputs": [],
@@ -980,7 +1168,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 33,
   "id": "01838139-52dd-48a0-8f62-bcdef351345e",
   "metadata": {},
   "outputs": [
@@ -988,8 +1176,8 @@
     "name": "stdout",
     "output_type": "stream",
     "text": [
-      "CPU times: user 25 μs, sys: 0 ns, total: 25 μs\n",
-      "Wall time: 27.2 μs\n"
+      "CPU times: user 23 μs, sys: 1 μs, total: 24 μs\n",
+      "Wall time: 24.1 μs\n"
     ]
    },
    {
@@ -998,7 +1186,7 @@
       "1"
      ]
     },
-     "execution_count": 28,
+     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }