Add files via upload
This commit is contained in:
Peter Norvig authored and 
GitHub committed
1 parent
9333518be1
commit
6986eda3cb
1 file changed
+87
-71
@@ -3,9 +3,7 @@
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8812d495-6014-4669-97ae-003be4f54605",
|
||||
"metadata": {
|
||||
"jp-MarkdownHeadingCollapsed": true
|
||||
},
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"<div style=\"text-align:right\"><i>Peter Norvig<br>Jan 2026</i></div>\n",
|
||||
"\n",
|
||||
@@ -295,7 +293,7 @@
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"id": "403a6451-c486-4868-acdc-7dd80a4c8f9d",
|
||||
"id": "5a705ef4-b2b7-4006-b518-6e030de8fd30",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -310,7 +308,7 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{p for p in TP if str(p) == str(p)[::-1]}"
|
||||
"{p for p in TP + TPr if str(p) == str(p)[::-1]}"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -480,7 +478,7 @@
|
||||
"\n",
|
||||
"# Primality Testing\n",
|
||||
"\n",
|
||||
"I was very impressed with the function [`sympy.isprime`](https://github.com/sympy/sympy/blob/master/sympy/ntheory/primetest.py), and I wanted to investigate the topic of primality testing. I'll start by defining `isprime_simple`, which follows the definition of a prime number almost verbatim:"
|
||||
"I was very impressed with the speed of the function [`sympy.isprime`](https://github.com/sympy/sympy/blob/master/sympy/ntheory/primetest.py), and I wanted to investigate the topic of primality testing. I'll start by defining `isprime_simple`, which follows the definition of a prime number almost verbatim:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -532,7 +530,7 @@
|
||||
"id": "05b8cbf7-ea10-4047-bcdb-a5ac7981384b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"That's a noticeable improvement, but testing a 24-digit number would still take hours. We need a massive breakthrough. \n",
|
||||
"For a 24-digit prime, `isprime_simple` would need to test trillions of trillions of divisors, while `isprime_faster` would \"only\" need to test trillions. That's the difference between never completing and taking a few hours. But we would like the test to take milliseconds. We need a massive breakthrough. \n",
|
||||
"\n",
|
||||
"Fortunately, [Pierre de Fermat](https://en.wikipedia.org/wiki/Pierre_de_Fermat) provided that breakthrough in 1640! \n",
|
||||
"\n",
|
||||
@@ -562,7 +560,7 @@
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"id": "476282dd-a48a-4eb5-8fed-8495d4878ec2",
|
||||
"id": "c8a59c5e-a01c-49be-a931-c1547c4ac712",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
@@ -571,8 +569,24 @@
|
||||
" if n <= 10:\n",
|
||||
" return n in (2, 3, 5, 7)\n",
|
||||
" else:\n",
|
||||
" return all(pow(a, n - 1, n) == 1 for a in sample(2, n, k))\n",
|
||||
"\n",
|
||||
" return all(pow(a, n - 1, n) == 1 for a in sample(2, n, k))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6746894c-00e7-41df-b4a3-d5c85862cc94",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We need a way to take random samples of integers from a range:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"id": "476282dd-a48a-4eb5-8fed-8495d4878ec2",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import random\n",
|
||||
"from typing import Iterable\n",
|
||||
" \n",
|
||||
@@ -601,7 +615,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"execution_count": 18,
|
||||
"id": "45a11a93-9ef3-48d1-80fa-e0113041346e",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
@@ -624,26 +638,26 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"execution_count": 19,
|
||||
"id": "9eba96df-c298-4d71-b68a-3bf36bf15095",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{3: 1.420427496905646,\n",
|
||||
" 4: 0.4795347492645433,\n",
|
||||
" 5: 0.13674836646356595,\n",
|
||||
" 6: 0.03350572569618954,\n",
|
||||
" 7: 0.010703085699607196,\n",
|
||||
" 8: 0.0010595970352474955,\n",
|
||||
" 9: 0.0010529807251878254,\n",
|
||||
"{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",
|
||||
" 11: 0,\n",
|
||||
" 12: 0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 18,
|
||||
"execution_count": 19,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -664,7 +678,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"execution_count": 20,
|
||||
"id": "8880c588-32f5-46e8-b59c-f9b49d492d3b",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
@@ -674,8 +688,8 @@
|
||||
"{3: 0,\n",
|
||||
" 4: 0,\n",
|
||||
" 5: 0,\n",
|
||||
" 6: 0.0010824385174922064,\n",
|
||||
" 7: 0,\n",
|
||||
" 6: 0.0003613604498214879,\n",
|
||||
" 7: 0.000356548185704557,\n",
|
||||
" 8: 0,\n",
|
||||
" 9: 0,\n",
|
||||
" 10: 0,\n",
|
||||
@@ -683,13 +697,13 @@
|
||||
" 12: 0}"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{d: false_prime_percent(sample(10**(d-1), 10**d, 200_000), k=25) for d in range(3, 13)}"
|
||||
"{d: false_prime_percent(sample(10**(d-1), 10**d, 300_000), k=25) for d in range(3, 13)}"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -697,44 +711,44 @@
|
||||
"id": "61ee44c3-f963-4fb1-ac20-4a27adc90f83",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"In most runs this produces no false primes; in some runs one or two get through. \n",
|
||||
"The results are randomized, and 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 high false prime percentage. [Some Carmichael numbers](https://oeis.org/A033502) are worse than others. \n",
|
||||
"There are a few numbers, called [**Carmichael numbers**](https://en.wikipedia.org/wiki/Carmichael_number), that have a very high false prime percentage. [Some Carmichael numbers](https://oeis.org/A033502) are worse than others. \n",
|
||||
"\n",
|
||||
"We can test some of them:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"execution_count": 21,
|
||||
"id": "7871aaa7-af1c-462b-85bb-d46818ccabea",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{561: 56.32,\n",
|
||||
" 1105: 68.95,\n",
|
||||
" 1729: 74.45,\n",
|
||||
" 2465: 72.27,\n",
|
||||
" 2821: 77.05,\n",
|
||||
" 6601: 80.08,\n",
|
||||
" 8911: 80.47,\n",
|
||||
" 41041: 71.11,\n",
|
||||
" 101101: 71.25,\n",
|
||||
" 294409: 95.19,\n",
|
||||
" 340561: 82.67,\n",
|
||||
" 825265: 59.95,\n",
|
||||
" 56052361: 98.96,\n",
|
||||
" 118901521: 99.4,\n",
|
||||
" 172947529: 99.43,\n",
|
||||
" 216821881: 99.55,\n",
|
||||
"{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",
|
||||
" 216821881: 99.4,\n",
|
||||
" 1299963601: 99.65,\n",
|
||||
" 2301745249: 99.82,\n",
|
||||
" 9624742921: 99.86}"
|
||||
" 2301745249: 99.76,\n",
|
||||
" 9624742921: 99.9}"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"execution_count": 21,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -756,7 +770,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"execution_count": 22,
|
||||
"id": "de492047-77b6-46a3-8a77-46308f219efd",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
@@ -772,25 +786,25 @@
|
||||
" 8911: 0,\n",
|
||||
" 41041: 0,\n",
|
||||
" 101101: 0,\n",
|
||||
" 294409: 8.05,\n",
|
||||
" 294409: 4,\n",
|
||||
" 340561: 0,\n",
|
||||
" 825265: 0,\n",
|
||||
" 56052361: 64.72,\n",
|
||||
" 118901521: 71.33,\n",
|
||||
" 172947529: 74.18,\n",
|
||||
" 216821881: 75.86,\n",
|
||||
" 1299963601: 85.84,\n",
|
||||
" 2301745249: 87.53,\n",
|
||||
" 9624742921: 92.3}"
|
||||
" 56052361: 58,\n",
|
||||
" 118901521: 66,\n",
|
||||
" 172947529: 80,\n",
|
||||
" 216821881: 86,\n",
|
||||
" 1299963601: 90,\n",
|
||||
" 2301745249: 87,\n",
|
||||
" 9624742921: 94}"
|
||||
]
|
||||
},
|
||||
"execution_count": 21,
|
||||
"execution_count": 22,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"{n: false_prime_percent([n] * 10_000, k=50) for n in carmichael_numbers}"
|
||||
"{n: false_prime_percent([n] * 100, k=50) for n in carmichael_numbers}"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -809,7 +823,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"execution_count": 23,
|
||||
"id": "be5f2fe7-b1e0-4a02-b7f0-7903e865c4b8",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
@@ -831,7 +845,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 23,
|
||||
"execution_count": 24,
|
||||
"id": "8109d9fc-bbae-4ef0-a75c-908d9746a557",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
@@ -870,13 +884,13 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 24,
|
||||
"execution_count": 25,
|
||||
"id": "b1ab0f3d-6ce6-4b8c-9d4d-bb78bb1ac432",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAINCAYAAACTVPhqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAwCFJREFUeJzs3Qd8U1X7B/Bf925pKbS07L33RkAEBRQQFHG+Aiou/IsiS30FJ8OBOFBQX9wDREBlD5EhG2RZ2bvQlrZ07yT/z3NCQnfTNmnW7+vnmpObm5tzb9pyn3uec46LTqfTgYiIiIiIqBiuxa0kIiIiIiJiwEBERERERKViCwMREREREZWIAQMREREREZWIAQMREREREZWIAQMREREREZWIAQMREREREZWIAQMREREREZXIveSXSGi1Wly+fBkBAQFwcXHhSSEiIiIiuydzN6empiIiIgKurqW3ITBgKIMEC3Xq1DHn90NEREREZBMuXryI2rVrl7oNA4YSzJ8/Xy15eXnGkxkYGGj+b4mIiIiIqIqlpKSom+KSRVMWF520R1CpJzMoKAjJyckMGIiIiIjI6a5x2emZiIiIiIhKxICBiIiIiIhKxICBiIiIiIhKxE7PRERERFYg3UhlcBWNRsPzTxbh4eEBNze3Su+HAQMRERFRFcvJycGVK1eQkZHBc08WI3OIyZCp/v7+ldoPAwYiIiKiKp4U9uzZs+rOr0ya5enpyclhySItWFevXsWlS5fQpEmTSrU0MGAoYx4GNhMSERGRuVsXJGiQMfB9fX15csliatSogXPnziE3N7dSAQM7PZdg/PjxiIqKwt69eyt8comIiIhKvAhz5WUYWT4lyRz4k0pERERERCViwEBERERERCViwEBERERkpzRaHXaeTsCvB6PVozy3pJtvvhnPPfecxfYv+faSRnPw4EHYonNVWL/69etj3rx5sAXs9ExERERkh9YevYLXfo/CleQs47paQd6YMbQlBrWuZZHPXLZsmRrb31KkI7gMNxsaGgpbVMfG62cpbGEgIiIissNg4anvDhQIFkRMcpZaL69bQkhICAICAiw2epSM5BMeHg53d9u8p+1m4/WzFAYMZHequvmViIioKsbMz8jJM2lJzcrFjN/+QXH/+hnWvfpblNrOlP3JZ1ckJemTTz5R4/t7e3sjLCwMI0eOLLDdM888o5agoCB1R/6VV14p8FmScvPGG2/g4YcfRmBgIB5//PEiKT9//vmner5u3Tp06NABPj4+uOWWWxAXF4c1a9agRYsW6r0PPPBAgUnwZNjaWbNmoUGDBuo97dq1w9KlS006xmvXruHBBx9UQ5LKe5s0aYIvv/xSvWau+plyfgpLSkrCY489puol+5TPOXToEKqCc4VHZPes0fxKRERkaZm5GrScvs4s+5JLzpiULLR5db1J20e9PhC+nuW7JNy3bx+effZZfPvtt+jZsycSExOxbdu2Att8/fXXePTRR7Fnzx61vQQEdevWxbhx44zbvPvuu5g+fTpmzJhR6ue9+uqr+Pjjj9W8FaNGjVKLl5cXfvjhB6SlpWHEiBH46KOPMHXqVLW9BAvfffcdFixYoC74t27dioceekhdbPft27fUz5ILdxlaXy745UL+1KlTyMzMNGv9TD0/+d1zzz0qIJF6SZCxcOFC9O/fHydOnFAtP5bEgIHsrvm1cOxtaH799KGODBqIiIiqwIULF+Dn54chQ4aoFKV69eqpO+yF8/3ff/99dQe+WbNmOHLkiHqe/4JY7pK/8MILxudyB784b775Jnr16qXKcpH94osv4vTp02jYsKFaJ60bmzdvVhfk2dnZmDlzJjZu3IgePXqo12W77du3q4vssgIGOTY5ls6dOxtbQspSnvqV5/wYSN0lsJCWCwlEDMHWihUrVMuJBBuWxICB7IKkHUnLQknNrzItibx+a8twuLmaZ5ISIoen1QDndwBpsYB/GFCvJ+Ba8ZlAiajifDzc1J1+U+w5m4gxX5Y9sexXY7uga4MQkz67vG699VYVJMgF8aBBg9Qid9Hzz1zdvXv3AhOHycX7e++9B41GY5x12HBRXpa2bdsay5L+JJ9juBg3rJMLaiEtApL+I3Us3EeicFBTnKeeegp33303Dhw4gNtuuw3Dhw9XrSjmql95zo+BpB5JS0X16tULrJeWDwlMLI0BQwnmz5+vFvnSyPrkj2Phjl2FgwZ5Xbbr0ajgLxMRFSPqN2DtVCDl8o11gRHAoDlAy2E8ZURVTC4cTU0L6t2khkrHlRb24m6kySVoeJC32s5SN9GkVUEuqCWHf/369SqtSNJy9u7di2rVqpm8H2mlMEX+kZnkXBUeqUnWSb8FIRfWYtWqVYiMjCywneHufGkGDx6M8+fPY/Xq1diwYYNK+xk/fry6o2+O+lWEHFOtWrXU+S6sPOe7otjpuQTygyH5a/KDT9YXl5pl1u2InJoEC0seLhgsiJQr+vXyOhHZLAkCpO+eKBwOGJ7L65ZucZeRggYMGIC3334bhw8fVulEf/zxh/H13bt3F9h+165dqj9B4bvn5tayZUsVGEhqUePGjQsskgZkCunrMHr0aNUPYt68efjss8/MXs/ynJ+OHTsiJiZGnfPCx1QVQ7yyhYHsQs0Ab7NuR+TUaUjSslBagt/aaUDzO5ieRGTDZKAP6btXeCCQ8CoaCGTlypU4c+YM+vTpg+DgYHU3Xu6gSy6+gVywT5w4EU888YRqjZBOv5JyY2nS+jFp0iQ8//zzqk433XQTkpOT8ddff6nRhSQQKI20lnTq1AmtWrVS/SFWrlypRjsyt/KcHwnMJGVJ0qMkQGvatCkuX76sWlEkFczU1K6KYsBAdkFyMKX5taS0JEPzqym5mkROTfosFG5ZKEAHpETrt2vQuworRkTlJUGB9N2TdFxpYZebZvLvYFX05ZM0GJnETdKQsrKy1J3xH3/8UV1kG8hwqZJj37VrV3XXfMKECRbvnGsgw7VKK4GMliSBjdRX7tK/9NJLZb7X09NTdVqWFhMZlah379746aefzF7H8pwfSWmSoOzll1/G2LFjcfXqVTUfhARs0j/C0lx05Rl81wmlpKSooaskMpWolKw7StKT3x0ost7wZ5GjJBGZ4MhS4JdHy97u7v8BbW6MqU5E5iMX2GfPnlVzBMgcBo5I5hlo3769Such652f0n7WynONyz4MZDdualIDnm5Ff2SlZYHBApGJZDQkc25HREQOjylJZDdWH76CHI0WDUJ9MXNEG8SlZldp8yuRQ6jbA/D0B3L0o4gU5aIfLUmGWCUickBPPvmk6sxcHJncTSZ7o4IYMJDdWLLvonoc1bkuejSy/IgARA5HhvRbM6X0YEEMms0Oz0RUKcUN/2krXn/9ddUpujhVlX7+pw2fn+IwYCC7cPpqGvadv6ZaEu7uWHBMZSIygSYP+O0Z4NCP+sCg0xjg5Lpi5mGYzXkYiMih1axZUy1kOgYMZBeW7r+kHvs2rYGagY7ZQYzIYvJy9B2d//0NcHEDRiwE2t7DmZ6JiMgkDBjI5uVptPjlesAwqnNta1eHyL7kZgKL/wOc2gC4eQIjvwRaDNG/5urGoVOJiKhMDBjI5m07Ga86OIf4eeKW5hy5hchk2anAD/cB57cD7j7A/T8AjW7hCSQionJhwEB209l5ePtIeLpzJGAik2QkAt+PBKL3A54BwINLOPIRERFVCAMGsmmJ6TnY+G+sKt/DdCQi06TFAd+OAGKPAj7BwEPLgMiOPHtERFQhTnG7duXKlWjWrJmatvyLL76wdnWoHFb8HY1cjQ5tIoPQohZn2iYqU3I08OVgfbAgk6+NWc1ggciRaTXA2W36WdzlUZ5beIbi5557zmL7P3fuHFxcXHDw4EFY019//YU2bdrAw8MDw4cPh7Nz+BaGvLw8TJw4EZs3b1bTX3fq1AkjRoxA9erVrV01KoNOp8s39wI7OxOVKfEM8M2dQNIFILA2MPo3oHojnjgiRxX1G7B2ajHDI8+x2PDIy5YtUxfRllKnTh1cuXIFoaHWnW9Jrh3bt2+PNWvWwN/fv9L7e/XVV7FixQqrB0IV5fAtDHv27EGrVq0QGRmpvvDBgwdj/fr11q4WmeCfyyk4FpOq+i0Ma8e5F4hKFXcMWDRYHyyENAQeWctggcjRg4UlDxcMFkTKFf16ed0CQkJCEBAQYJF95+TkwM3NDeHh4XB3t+497dOnT+OWW25B7dq1Ua1aNdgKOUfWYPMBw9atWzF06FBERESoJiqJzgqbP38+6tevD29vb3Tr1k0FCQaXL19WwYKBlKOjo6us/lRxhtaFga3CEeRrubsZRHbv8kHgq9uBtBigZktg7FqgWh1r14qIykOnA3LSTVuyUvSztkNX3I70D9LyINuZsj/57AqkJH3yyScq3Vuuv8LCwjBy5MgC2z3zzDNqkQwPaTF45ZVXVPaAgVy7vfHGG3j44YfVDMuPP/54kZQkmRFZnq9btw4dOnSAj4+PupCPi4tTd/9btGih3vvAAw8gIyPDuG+tVotZs2ahQYMG6j3t2rXD0qVLyzw+w+cnJCTgkUceUeWvvvoKGo0Gjz76qHF/kur+wQcfFHiv1LVr167w8/NTQUavXr1w/vx59f7XXnsNhw4dUvsz7FMkJSXhscceQ40aNdRxyLHJdvlbJqSlQ1Lq5bPlXFuDzackpaenqy9ZvrS77rqryOuLFy9WzUYLFixQwcK8efMwcOBAHD9+nLP42bGsXA1+Pai/a3JPJ6YjEZXowm7g+3uA7GQgooO+g7NvCE8Ykb3JzQBmRphpZzp9y8NsE28cvHQZ8PQr1yfs27cPzz77LL799lv07NkTiYmJ2LZtW4Ftvv76a3WRLTdyZXsJCOrWrYtx48YZt3n33Xcxffp0zJgxo9TPkwvnjz/+GL6+vhg1apRavLy88MMPPyAtLU2lm3/00UeYOnWq2l6Che+++05dH0pQIzegH3roIXVLine truncated
|
||||
"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
|
||||
"text/plain": [
|
||||
"<Figure size 900x600 with 1 Axes>"
|
||||
]
|
||||
@@ -915,7 +929,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 25,
|
||||
"execution_count": 26,
|
||||
"id": "4f377075-e352-46f9-a716-f1de673a7086",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
@@ -930,18 +944,20 @@
|
||||
"id": "1347399a-24f2-4d5f-aa7f-822951126bf5",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This was fun for me! I remember learning being excited about Fermat's Little Theorem in [Prof. Michael Rosen's](https://mathematics.brown.edu/people/michael-rosen) number theory class in college, but I never worked through the details of how often the test gives a false prime result until now.\n",
|
||||
"This was fun for me! I remember being amazed by Fermat's Little Theorem when I learned about it in [Prof. Michael Rosen's](https://mathematics.brown.edu/people/michael-rosen) number theory class in college, but I never worked through the details of how often the test gives a false prime result until now. \n",
|
||||
"\n",
|
||||
"That number theory class was important for me in several ways. First, it was hard fun. Second, it was the first class that felt like actually *doing* math, rather than being taught some techniques to prepare you to do math some time in the future. And third, I learned a valuable lesson: there was a take-home final exam, handed out Friday and due Monday. I worked on it all day Saturday, and couldn't answer a single one of the ten problems. At that point, I thought to myself \"This never happens on my computer science assignments; I guess I'm more cut out for CS than math.\" On Sunday I came back to it and eventually got answers to most of the problems. (I remember that on one problem, Prof. Rosen said \"You had a very elegant solution for that one, better than the answer I expected. But then I realized it was so good, it could be extended to be a trivial solution to this known difficult problem, and so I looked more closely and noticed a flaw in your argument. Too bad.\")\n",
|
||||
"\n",
|
||||
"# Note on Modular Exponentiation\n",
|
||||
"\n",
|
||||
"Just one more thing: none of this would work unless we can efficiently compute *a*<sup>(*n* - 1)</sup> (mod *n*). How does the `pow` builtin function do it? When *a* and *n* are 24-digit numbers, if we naively tried to compute `a ** (n - 1)`, we'd have two problems: we'd need about 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 because log<sub>2</sub>(10<sup>24</sup>) is only 80. \n",
|
||||
"Just one more thing: none of this would work unless we can efficiently compute *a*<sup>(*n* - 1)</sup> (mod *n*). How does the `pow` builtin function do it? When *a* and *n* are 24-digit numbers, if we naively tried to compute `a ** (n - 1)`, 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:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 26,
|
||||
"execution_count": 27,
|
||||
"id": "e71eb140-7e1a-4146-87a7-8c4eebf61486",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
@@ -959,12 +975,12 @@
|
||||
"id": "6baf40c3-db36-40ba-b2d7-45afa21f811f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can see this function is fast and that it correctly computes the remainder of 1 for the follwoing:"
|
||||
"We can see this function is fast and that it correctly computes the remainder of 1 for the following:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 27,
|
||||
"execution_count": 28,
|
||||
"id": "01838139-52dd-48a0-8f62-bcdef351345e",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
@@ -972,8 +988,8 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"CPU times: user 25 μs, sys: 1 μs, total: 26 μs\n",
|
||||
"Wall time: 28.1 μs\n"
|
||||
"CPU times: user 25 μs, sys: 0 ns, total: 25 μs\n",
|
||||
"Wall time: 27.2 μs\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -982,7 +998,7 @@
|
||||
"1"
|
||||
]
|
||||
},
|
||||
"execution_count": 27,
|
||||
"execution_count": 28,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
|
||||
Reference in new issue
Block a user