diff options
Diffstat (limited to 'src/Lecture7/slides/.ipynb_checkpoints')
4 files changed, 1269 insertions, 0 deletions
diff --git a/src/Lecture7/slides/.ipynb_checkpoints/X1-ComputationalComplexity-checkpoint.ipynb b/src/Lecture7/slides/.ipynb_checkpoints/X1-ComputationalComplexity-checkpoint.ipynb new file mode 100644 index 0000000..ab7dd49 --- /dev/null +++ b/src/Lecture7/slides/.ipynb_checkpoints/X1-ComputationalComplexity-checkpoint.ipynb | |||
| @@ -0,0 +1,394 @@ | |||
| 1 | { | ||
| 2 | "cells": [ | ||
| 3 | { | ||
| 4 | "cell_type": "markdown", | ||
| 5 | "metadata": {}, | ||
| 6 | "source": [ | ||
| 7 | "# Nested loops\n", | ||
| 8 | "\n", | ||
| 9 | "The following two functions compute sum and product of matrices, respectively.\n", | ||
| 10 | "\n", | ||
| 11 | "By counting the nested loops it is easy to see that `add()` is $O(n^2)$ while `prod()` is $O(n^3)$." | ||
| 12 | ] | ||
| 13 | }, | ||
| 14 | { | ||
| 15 | "cell_type": "code", | ||
| 16 | "execution_count": 37, | ||
| 17 | "metadata": {}, | ||
| 18 | "outputs": [ | ||
| 19 | { | ||
| 20 | "name": "stdout", | ||
| 21 | "output_type": "stream", | ||
| 22 | "text": [ | ||
| 23 | "Time for add: 0.00012074300000008975\n", | ||
| 24 | "Time for prod: 0.00036587199999971176\n" | ||
| 25 | ] | ||
| 26 | } | ||
| 27 | ], | ||
| 28 | "source": [ | ||
| 29 | "from random import randint\n", | ||
| 30 | "import time\n", | ||
| 31 | "\n", | ||
| 32 | "def add(A, B):\n", | ||
| 33 | " S = [[0] * len(A) for i in range(len(A))]\n", | ||
| 34 | " for i in range(len(A)):\n", | ||
| 35 | " for j in range(len(A)):\n", | ||
| 36 | " S[i][j] = A[i][j] + B[i][j]\n", | ||
| 37 | " return S\n", | ||
| 38 | "\n", | ||
| 39 | "def prod(A, B):\n", | ||
| 40 | " S = [[0] * len(A) for i in range(len(A))]\n", | ||
| 41 | " for i in range(len(A)):\n", | ||
| 42 | " for j in range(len(A)):\n", | ||
| 43 | " for k in range(len(A)):\n", | ||
| 44 | " S[i][j] = S[i][j] + A[i][k] * B[k][j]\n", | ||
| 45 | " return S\n", | ||
| 46 | "\n", | ||
| 47 | "N = 10\n", | ||
| 48 | "A = [ [randint(0,100) for i in range(N)] for j in range(N) ]\n", | ||
| 49 | "B = [ [randint(0,100) for i in range(N)] for j in range(N) ]\n", | ||
| 50 | "\n", | ||
| 51 | "t0 = time.process_time()\n", | ||
| 52 | "add(A,B)\n", | ||
| 53 | "t1 = time.process_time()\n", | ||
| 54 | "prod(A,B)\n", | ||
| 55 | "t2 = time.process_time()\n", | ||
| 56 | "\n", | ||
| 57 | "print(\"Time for add: \", t1-t0)\n", | ||
| 58 | "print(\"Time for prod:\", t2-t1)" | ||
| 59 | ] | ||
| 60 | }, | ||
| 61 | { | ||
| 62 | "cell_type": "markdown", | ||
| 63 | "metadata": {}, | ||
| 64 | "source": [ | ||
| 65 | "# Sorting a list, slow version\n", | ||
| 66 | "\n", | ||
| 67 | "The following code implements a slow version of the so-called *insertion sort* alogithm\n", | ||
| 68 | "\n", | ||
| 69 | "Complexity: $O(n^2)$." | ||
| 70 | ] | ||
| 71 | }, | ||
| 72 | { | ||
| 73 | "cell_type": "code", | ||
| 74 | "execution_count": 61, | ||
| 75 | "metadata": {}, | ||
| 76 | "outputs": [ | ||
| 77 | { | ||
| 78 | "name": "stdout", | ||
| 79 | "output_type": "stream", | ||
| 80 | "text": [ | ||
| 81 | "Running time: 1.1012288430000012\n" | ||
| 82 | ] | ||
| 83 | } | ||
| 84 | ], | ||
| 85 | "source": [ | ||
| 86 | "from random import randint\n", | ||
| 87 | "import time\n", | ||
| 88 | "\n", | ||
| 89 | "def correct_position(e, S):\n", | ||
| 90 | " for i in range(len(S)):\n", | ||
| 91 | " if S[i] > e:\n", | ||
| 92 | " return i\n", | ||
| 93 | " return len(S)\n", | ||
| 94 | "\n", | ||
| 95 | "def sort_list(L):\n", | ||
| 96 | " S = []\n", | ||
| 97 | " for e in L:\n", | ||
| 98 | " cp = correct_position(e, S)\n", | ||
| 99 | " S.insert(cp, e)\n", | ||
| 100 | " return S\n", | ||
| 101 | "\n", | ||
| 102 | "N = 10000\n", | ||
| 103 | "L = [randint(0,10**9) for i in range(N)]\n", | ||
| 104 | "\n", | ||
| 105 | "t0 = time.process_time()\n", | ||
| 106 | "sort_list(L)\n", | ||
| 107 | "t1 = time.process_time()\n", | ||
| 108 | "\n", | ||
| 109 | "print(\"Running time:\", t1-t0)" | ||
| 110 | ] | ||
| 111 | }, | ||
| 112 | { | ||
| 113 | "cell_type": "markdown", | ||
| 114 | "metadata": {}, | ||
| 115 | "source": [ | ||
| 116 | "# Binary search\n", | ||
| 117 | "\n", | ||
| 118 | "The following code implements a binary search.\n", | ||
| 119 | "\n", | ||
| 120 | "Complexity: $O(\\log_2(n))$" | ||
| 121 | ] | ||
| 122 | }, | ||
| 123 | { | ||
| 124 | "cell_type": "code", | ||
| 125 | "execution_count": 53, | ||
| 126 | "metadata": {}, | ||
| 127 | "outputs": [ | ||
| 128 | { | ||
| 129 | "name": "stdout", | ||
| 130 | "output_type": "stream", | ||
| 131 | "text": [ | ||
| 132 | "The correct position of e = 658230309 in L is:\n", | ||
| 133 | "... 658211821 658224379 e 658234625 658246765 ...\n", | ||
| 134 | "\n", | ||
| 135 | "Time for sorting: 0.021211020999999164\n", | ||
| 136 | "Time for searching: 7.820199999741817e-05\n" | ||
| 137 | ] | ||
| 138 | } | ||
| 139 | ], | ||
| 140 | "source": [ | ||
| 141 | "from random import randint\n", | ||
| 142 | "import time\n", | ||
| 143 | "\n", | ||
| 144 | "def binary_search(e, S, start, end):\n", | ||
| 145 | " if start == end:\n", | ||
| 146 | " return start\n", | ||
| 147 | " midpoint = (start+end) // 2\n", | ||
| 148 | " if e < S[midpoint]:\n", | ||
| 149 | " return binary_search(e, S, start, midpoint)\n", | ||
| 150 | " else:\n", | ||
| 151 | " return binary_search(e, S, midpoint+1, end)\n", | ||
| 152 | " \n", | ||
| 153 | "N = 100000\n", | ||
| 154 | "L = [randint(0,10**9) for i in range(N)]\n", | ||
| 155 | "e = randint(0,10**9)\n", | ||
| 156 | "\n", | ||
| 157 | "t0 = time.process_time()\n", | ||
| 158 | "L.sort() # Using Python's sort()\n", | ||
| 159 | "t1 = time.process_time()\n", | ||
| 160 | "i = binary_search(e, L, 0, len(L))\n", | ||
| 161 | "t2 = time.process_time()\n", | ||
| 162 | "print(\"The correct position of e =\", e, \"in L is:\")\n", | ||
| 163 | "print(\"...\", L[i-2], L[i-1], \"e\", L[i], L[i+1], \"...\")\n", | ||
| 164 | "print(\"\")\n", | ||
| 165 | "print(\"Time for sorting: \", t1-t0)\n", | ||
| 166 | "print(\"Time for searching:\", t2-t1)\n" | ||
| 167 | ] | ||
| 168 | }, | ||
| 169 | { | ||
| 170 | "cell_type": "markdown", | ||
| 171 | "metadata": {}, | ||
| 172 | "source": [ | ||
| 173 | "# Sorting a list, fast version (with binary_search)\n", | ||
| 174 | "\n", | ||
| 175 | "The following code uses the function `binary_search()` above instead of `correct_position()` in our insertion sort algorithm.\n", | ||
| 176 | "\n", | ||
| 177 | "Complexity: $O(n\\log_2(n))$" | ||
| 178 | ] | ||
| 179 | }, | ||
| 180 | { | ||
| 181 | "cell_type": "code", | ||
| 182 | "execution_count": 69, | ||
| 183 | "metadata": {}, | ||
| 184 | "outputs": [ | ||
| 185 | { | ||
| 186 | "name": "stdout", | ||
| 187 | "output_type": "stream", | ||
| 188 | "text": [ | ||
| 189 | "Running time: 0.03710268399998995\n" | ||
| 190 | ] | ||
| 191 | } | ||
| 192 | ], | ||
| 193 | "source": [ | ||
| 194 | "from random import randint\n", | ||
| 195 | "import time\n", | ||
| 196 | "\n", | ||
| 197 | "def binary_search(e, S, start, end):\n", | ||
| 198 | " if start == end:\n", | ||
| 199 | " return start\n", | ||
| 200 | " midpoint = (start+end) // 2\n", | ||
| 201 | " if e < S[midpoint]:\n", | ||
| 202 | " return binary_search(e, S, start, midpoint)\n", | ||
| 203 | " else:\n", | ||
| 204 | " return binary_search(e, S, midpoint+1, end)\n", | ||
| 205 | " \n", | ||
| 206 | "def sort_list(L):\n", | ||
| 207 | " S = []\n", | ||
| 208 | " for e in L:\n", | ||
| 209 | " cp = binary_search(e, S, 0, len(S)) # Changed here\n", | ||
| 210 | " S.insert(cp, e)\n", | ||
| 211 | " return S\n", | ||
| 212 | " \n", | ||
| 213 | "N = 10000\n", | ||
| 214 | "L = [randint(0,10**9) for i in range(N)]\n", | ||
| 215 | "\n", | ||
| 216 | "t0 = time.process_time()\n", | ||
| 217 | "sort_list(L)\n", | ||
| 218 | "t1 = time.process_time()\n", | ||
| 219 | "\n", | ||
| 220 | "print(\"Running time:\", t1-t0)" | ||
| 221 | ] | ||
| 222 | }, | ||
| 223 | { | ||
| 224 | "cell_type": "markdown", | ||
| 225 | "metadata": {}, | ||
| 226 | "source": [ | ||
| 227 | "# Fast exponentiation\n", | ||
| 228 | "\n", | ||
| 229 | "The following cell contains two functions for computing $a^n$ ($n$ non-negative integer): a slow one that runs in $O(n)$ and a fast one that runs in $O(\\log_2(n))$. We compare these two also with Python's built-in operator `**`.\n", | ||
| 230 | "\n", | ||
| 231 | "Complexity: $O(n)$ for the slow algorithm, $O(\\log_2(n))$ for the other two." | ||
| 232 | ] | ||
| 233 | }, | ||
| 234 | { | ||
| 235 | "cell_type": "code", | ||
| 236 | "execution_count": 30, | ||
| 237 | "metadata": {}, | ||
| 238 | "outputs": [ | ||
| 239 | { | ||
| 240 | "name": "stdout", | ||
| 241 | "output_type": "stream", | ||
| 242 | "text": [ | ||
| 243 | "2.71828179834636\n", | ||
| 244 | "2.7182817863957984\n", | ||
| 245 | "2.7182817983473577\n", | ||
| 246 | "Time for slow_power(): 3.234879998000004\n", | ||
| 247 | "Time for fast_power(): 9.059099999575437e-05\n", | ||
| 248 | "Time for Python's **: 0.00010159500000384014\n" | ||
| 249 | ] | ||
| 250 | } | ||
| 251 | ], | ||
| 252 | "source": [ | ||
| 253 | "import time\n", | ||
| 254 | "\n", | ||
| 255 | "def slow_power(a, n):\n", | ||
| 256 | " r = 1\n", | ||
| 257 | " for i in range(n):\n", | ||
| 258 | " r = r * a\n", | ||
| 259 | " return r\n", | ||
| 260 | "\n", | ||
| 261 | "def fast_power(a, n):\n", | ||
| 262 | " if n == 0:\n", | ||
| 263 | " return 1\n", | ||
| 264 | " if n%2 == 0:\n", | ||
| 265 | " return fast_power(a*a, n//2)\n", | ||
| 266 | " else:\n", | ||
| 267 | " return a * fast_power(a, n-1)\n", | ||
| 268 | "\n", | ||
| 269 | "a = 1.00000001\n", | ||
| 270 | "n = 100000000\n", | ||
| 271 | "\n", | ||
| 272 | "t0 = time.process_time()\n", | ||
| 273 | "print(slow_power(a, n))\n", | ||
| 274 | "t1 = time.process_time()\n", | ||
| 275 | "print(fast_power(a, n))\n", | ||
| 276 | "t2 = time.process_time()\n", | ||
| 277 | "print(a**n)\n", | ||
| 278 | "t3 = time.process_time()\n", | ||
| 279 | "\n", | ||
| 280 | "print(\"Time for slow_power():\", t1-t0)\n", | ||
| 281 | "print(\"Time for fast_power():\", t2-t1)\n", | ||
| 282 | "print(\"Time for Python's **: \", t3-t2)" | ||
| 283 | ] | ||
| 284 | }, | ||
| 285 | { | ||
| 286 | "cell_type": "markdown", | ||
| 287 | "metadata": {}, | ||
| 288 | "source": [ | ||
| 289 | "# Fast gcd\n", | ||
| 290 | "\n", | ||
| 291 | "Complexity: $O(\\log_2(n))$" | ||
| 292 | ] | ||
| 293 | }, | ||
| 294 | { | ||
| 295 | "cell_type": "code", | ||
| 296 | "execution_count": 31, | ||
| 297 | "metadata": {}, | ||
| 298 | "outputs": [ | ||
| 299 | { | ||
| 300 | "name": "stdout", | ||
| 301 | "output_type": "stream", | ||
| 302 | "text": [ | ||
| 303 | "126\n", | ||
| 304 | "Running time: 0.00017707599999994272\n" | ||
| 305 | ] | ||
| 306 | } | ||
| 307 | ], | ||
| 308 | "source": [ | ||
| 309 | "import time\n", | ||
| 310 | "\n", | ||
| 311 | "def gcd(a, b):\n", | ||
| 312 | " if b == 0:\n", | ||
| 313 | " return a\n", | ||
| 314 | " else:\n", | ||
| 315 | " return gcd(b, a%b)\n", | ||
| 316 | "\n", | ||
| 317 | "t0 = time.process_time()\n", | ||
| 318 | "print(gcd(155275387236018, 572335397352432))\n", | ||
| 319 | "t1 = time.process_time()\n", | ||
| 320 | "\n", | ||
| 321 | "print(\"Running time:\", t1-t0)" | ||
| 322 | ] | ||
| 323 | }, | ||
| 324 | { | ||
| 325 | "cell_type": "markdown", | ||
| 326 | "metadata": {}, | ||
| 327 | "source": [ | ||
| 328 | "# Fibonacci numbers\n", | ||
| 329 | "\n", | ||
| 330 | "In the following cell there are two functions that compute the $n$-th Fibonacci number. They are almost the same, but the second one memorizes the results in a list to avoid computing them multiple times, and it is much much faster.\n", | ||
| 331 | "\n", | ||
| 332 | "Complexity: $O\\left(\\left(\\frac{1+\\sqrt 5}{2}\\right)^n\\right)\\sim O(1.6^n)$ for the slow version, $O(n)$ for the fast version." | ||
| 333 | ] | ||
| 334 | }, | ||
| 335 | { | ||
| 336 | "cell_type": "code", | ||
| 337 | "execution_count": null, | ||
| 338 | "metadata": {}, | ||
| 339 | "outputs": [], | ||
| 340 | "source": [ | ||
| 341 | "import time\n", | ||
| 342 | "\n", | ||
| 343 | "F_memorized = [-1] * (10**6)\n", | ||
| 344 | "\n", | ||
| 345 | "def F_slow(n):\n", | ||
| 346 | " if n <= 1:\n", | ||
| 347 | " return n\n", | ||
| 348 | " else:\n", | ||
| 349 | " return F_slow(n-1) + F_slow(n-2)\n", | ||
| 350 | " \n", | ||
| 351 | "def F_fast(n):\n", | ||
| 352 | " if F_memorized[n] == -1:\n", | ||
| 353 | " if n <= 1:\n", | ||
| 354 | " F_memorized[n] = n\n", | ||
| 355 | " else:\n", | ||
| 356 | " F_memorized[n] = F_fast(n-1) + F_fast(n-2)\n", | ||
| 357 | " \n", | ||
| 358 | " return F_memorized[n]\n", | ||
| 359 | "\n", | ||
| 360 | "n = 40\n", | ||
| 361 | "\n", | ||
| 362 | "t0 = time.process_time()\n", | ||
| 363 | "print(F_slow(n))\n", | ||
| 364 | "t1 = time.process_time()\n", | ||
| 365 | "print(F_fast(n))\n", | ||
| 366 | "t2 = time.process_time()\n", | ||
| 367 | "\n", | ||
| 368 | "print(\"Time for F_slow:\", t1-t0)\n", | ||
| 369 | "print(\"Time for F_fast:\", t2-t1)" | ||
| 370 | ] | ||
| 371 | } | ||
| 372 | ], | ||
| 373 | "metadata": { | ||
| 374 | "kernelspec": { | ||
| 375 | "display_name": "Python 3", | ||
| 376 | "language": "python", | ||
| 377 | "name": "python3" | ||
| 378 | }, | ||
| 379 | "language_info": { | ||
| 380 | "codemirror_mode": { | ||
| 381 | "name": "ipython", | ||
| 382 | "version": 3 | ||
| 383 | }, | ||
| 384 | "file_extension": ".py", | ||
| 385 | "mimetype": "text/x-python", | ||
| 386 | "name": "python", | ||
| 387 | "nbconvert_exporter": "python", | ||
| 388 | "pygments_lexer": "ipython3", | ||
| 389 | "version": "3.8.5" | ||
| 390 | } | ||
| 391 | }, | ||
| 392 | "nbformat": 4, | ||
| 393 | "nbformat_minor": 4 | ||
| 394 | } | ||
diff --git a/src/Lecture7/slides/.ipynb_checkpoints/X1-ComputationalComplexity-notebook-checkpoint.ipynb b/src/Lecture7/slides/.ipynb_checkpoints/X1-ComputationalComplexity-notebook-checkpoint.ipynb new file mode 100644 index 0000000..16a6d40 --- /dev/null +++ b/src/Lecture7/slides/.ipynb_checkpoints/X1-ComputationalComplexity-notebook-checkpoint.ipynb | |||
| @@ -0,0 +1,412 @@ | |||
| 1 | { | ||
| 2 | "cells": [ | ||
| 3 | { | ||
| 4 | "cell_type": "markdown", | ||
| 5 | "metadata": {}, | ||
| 6 | "source": [ | ||
| 7 | "# Nested loops\n", | ||
| 8 | "\n", | ||
| 9 | "The following two functions compute sum and product of matrices, respectively.\n", | ||
| 10 | "\n", | ||
| 11 | "By counting the nested loops it is easy to see that `add()` is $O(n^2)$ while `prod()` is $O(n^3)$." | ||
| 12 | ] | ||
| 13 | }, | ||
| 14 | { | ||
| 15 | "cell_type": "code", | ||
| 16 | "execution_count": 2, | ||
| 17 | "metadata": {}, | ||
| 18 | "outputs": [ | ||
| 19 | { | ||
| 20 | "name": "stdout", | ||
| 21 | "output_type": "stream", | ||
| 22 | "text": [ | ||
| 23 | "Time for add: 0.005766554000000035\n", | ||
| 24 | "Time for prod: 1.3871021639999999\n" | ||
| 25 | ] | ||
| 26 | } | ||
| 27 | ], | ||
| 28 | "source": [ | ||
| 29 | "from random import randint\n", | ||
| 30 | "import time\n", | ||
| 31 | "\n", | ||
| 32 | "def add(A, B):\n", | ||
| 33 | " S = [[0] * len(A) for i in range(len(A))]\n", | ||
| 34 | " for i in range(len(A)):\n", | ||
| 35 | " for j in range(len(A)):\n", | ||
| 36 | " S[i][j] = A[i][j] + B[i][j]\n", | ||
| 37 | " return S\n", | ||
| 38 | "\n", | ||
| 39 | "def prod(A, B):\n", | ||
| 40 | " S = [[0] * len(A) for i in range(len(A))]\n", | ||
| 41 | " for i in range(len(A)):\n", | ||
| 42 | " for j in range(len(A)):\n", | ||
| 43 | " for k in range(len(A)):\n", | ||
| 44 | " S[i][j] = S[i][j] + A[i][k] * B[k][j]\n", | ||
| 45 | " return S\n", | ||
| 46 | "\n", | ||
| 47 | "N = 200\n", | ||
| 48 | "A = [ [randint(0,100) for i in range(N)] for j in range(N) ]\n", | ||
| 49 | "B = [ [randint(0,100) for i in range(N)] for j in range(N) ]\n", | ||
| 50 | "\n", | ||
| 51 | "t0 = time.process_time()\n", | ||
| 52 | "add(A,B)\n", | ||
| 53 | "t1 = time.process_time()\n", | ||
| 54 | "prod(A,B)\n", | ||
| 55 | "t2 = time.process_time()\n", | ||
| 56 | "\n", | ||
| 57 | "print(\"Time for add: \", t1-t0)\n", | ||
| 58 | "print(\"Time for prod:\", t2-t1)" | ||
| 59 | ] | ||
| 60 | }, | ||
| 61 | { | ||
| 62 | "cell_type": "markdown", | ||
| 63 | "metadata": {}, | ||
| 64 | "source": [ | ||
| 65 | "# Sorting a list, slow version\n", | ||
| 66 | "\n", | ||
| 67 | "The following code implements a slow version of the so-called *insertion sort* alogithm\n", | ||
| 68 | "\n", | ||
| 69 | "Complexity: $O(n^2)$." | ||
| 70 | ] | ||
| 71 | }, | ||
| 72 | { | ||
| 73 | "cell_type": "code", | ||
| 74 | "execution_count": 1, | ||
| 75 | "metadata": {}, | ||
| 76 | "outputs": [ | ||
| 77 | { | ||
| 78 | "name": "stdout", | ||
| 79 | "output_type": "stream", | ||
| 80 | "text": [ | ||
| 81 | "Running time: 1.1191449070000001\n" | ||
| 82 | ] | ||
| 83 | } | ||
| 84 | ], | ||
| 85 | "source": [ | ||
| 86 | "from random import randint\n", | ||
| 87 | "import time\n", | ||
| 88 | "\n", | ||
| 89 | "def correct_position(e, S):\n", | ||
| 90 | " for i in range(len(S)):\n", | ||
| 91 | " if S[i] > e:\n", | ||
| 92 | " return i\n", | ||
| 93 | " return len(S)\n", | ||
| 94 | "\n", | ||
| 95 | "def sort_list(L):\n", | ||
| 96 | " S = []\n", | ||
| 97 | " for e in L:\n", | ||
| 98 | " cp = correct_position(e, S)\n", | ||
| 99 | " S.insert(cp, e)\n", | ||
| 100 | " return S\n", | ||
| 101 | "\n", | ||
| 102 | "N = 10000\n", | ||
| 103 | "L = [randint(0,10**9) for i in range(N)]\n", | ||
| 104 | "\n", | ||
| 105 | "t0 = time.process_time()\n", | ||
| 106 | "sort_list(L)\n", | ||
| 107 | "t1 = time.process_time()\n", | ||
| 108 | "\n", | ||
| 109 | "print(\"Running time:\", t1-t0)" | ||
| 110 | ] | ||
| 111 | }, | ||
| 112 | { | ||
| 113 | "cell_type": "markdown", | ||
| 114 | "metadata": {}, | ||
| 115 | "source": [ | ||
| 116 | "# Binary search\n", | ||
| 117 | "\n", | ||
| 118 | "The following code implements a binary search.\n", | ||
| 119 | "\n", | ||
| 120 | "Complexity: $O(\\log_2(n))$" | ||
| 121 | ] | ||
| 122 | }, | ||
| 123 | { | ||
| 124 | "cell_type": "code", | ||
| 125 | "execution_count": 3, | ||
| 126 | "metadata": {}, | ||
| 127 | "outputs": [ | ||
| 128 | { | ||
| 129 | "name": "stdout", | ||
| 130 | "output_type": "stream", | ||
| 131 | "text": [ | ||
| 132 | "The correct position of e = 216197744 in L is:\n", | ||
| 133 | "... 216196218 216197540 e 216198673 216198962 ...\n", | ||
| 134 | "\n", | ||
| 135 | "Time for sorting: 0.26413054400000036\n", | ||
| 136 | "Time for searching: 9.616099999965044e-05\n" | ||
| 137 | ] | ||
| 138 | } | ||
| 139 | ], | ||
| 140 | "source": [ | ||
| 141 | "from random import randint\n", | ||
| 142 | "import time\n", | ||
| 143 | "\n", | ||
| 144 | "def binary_search(e, S, start, end):\n", | ||
| 145 | " if start == end:\n", | ||
| 146 | " return start\n", | ||
| 147 | " midpoint = (start+end) // 2\n", | ||
| 148 | " if e < S[midpoint]:\n", | ||
| 149 | " return binary_search(e, S, start, midpoint)\n", | ||
| 150 | " else:\n", | ||
| 151 | " return binary_search(e, S, midpoint+1, end)\n", | ||
| 152 | " \n", | ||
| 153 | "N = 1000000\n", | ||
| 154 | "L = [randint(0,10**9) for i in range(N)]\n", | ||
| 155 | "e = randint(0,10**9)\n", | ||
| 156 | "\n", | ||
| 157 | "t0 = time.process_time()\n", | ||
| 158 | "L.sort() # Using Python's sort()\n", | ||
| 159 | "t1 = time.process_time()\n", | ||
| 160 | "i = binary_search(e, L, 0, len(L))\n", | ||
| 161 | "t2 = time.process_time()\n", | ||
| 162 | "print(\"The correct position of e =\", e, \"in L is:\")\n", | ||
| 163 | "print(\"...\", L[i-2], L[i-1], \"e\", L[i], L[i+1], \"...\")\n", | ||
| 164 | "print(\"\")\n", | ||
| 165 | "print(\"Time for sorting: \", t1-t0)\n", | ||
| 166 | "print(\"Time for searching:\", t2-t1)\n" | ||
| 167 | ] | ||
| 168 | }, | ||
| 169 | { | ||
| 170 | "cell_type": "markdown", | ||
| 171 | "metadata": {}, | ||
| 172 | "source": [ | ||
| 173 | "# Sorting a list, fast version (with binary_search)\n", | ||
| 174 | "\n", | ||
| 175 | "The following code uses the function `binary_search()` above instead of `correct_position()` in our insertion sort algorithm.\n", | ||
| 176 | "\n", | ||
| 177 | "Complexity: $O(n\\log_2(n))$" | ||
| 178 | ] | ||
| 179 | }, | ||
| 180 | { | ||
| 181 | "cell_type": "code", | ||
| 182 | "execution_count": 69, | ||
| 183 | "metadata": {}, | ||
| 184 | "outputs": [ | ||
| 185 | { | ||
| 186 | "name": "stdout", | ||
| 187 | "output_type": "stream", | ||
| 188 | "text": [ | ||
| 189 | "Running time: 0.03710268399998995\n" | ||
| 190 | ] | ||
| 191 | } | ||
| 192 | ], | ||
| 193 | "source": [ | ||
| 194 | "from random import randint\n", | ||
| 195 | "import time\n", | ||
| 196 | "\n", | ||
| 197 | "def binary_search(e, S, start, end):\n", | ||
| 198 | " if start == end:\n", | ||
| 199 | " return start\n", | ||
| 200 | " midpoint = (start+end) // 2\n", | ||
| 201 | " if e < S[midpoint]:\n", | ||
| 202 | " return binary_search(e, S, start, midpoint)\n", | ||
| 203 | " else:\n", | ||
| 204 | " return binary_search(e, S, midpoint+1, end)\n", | ||
| 205 | " \n", | ||
| 206 | "def sort_list(L):\n", | ||
| 207 | " S = []\n", | ||
| 208 | " for e in L:\n", | ||
| 209 | " cp = binary_search(e, S, 0, len(S)) # Changed here\n", | ||
| 210 | " S.insert(cp, e)\n", | ||
| 211 | " return S\n", | ||
| 212 | " \n", | ||
| 213 | "N = 10000\n", | ||
| 214 | "L = [randint(0,10**9) for i in range(N)]\n", | ||
| 215 | "\n", | ||
| 216 | "t0 = time.process_time()\n", | ||
| 217 | "sort_list(L)\n", | ||
| 218 | "t1 = time.process_time()\n", | ||
| 219 | "\n", | ||
| 220 | "print(\"Running time:\", t1-t0)" | ||
| 221 | ] | ||
| 222 | }, | ||
| 223 | { | ||
| 224 | "cell_type": "markdown", | ||
| 225 | "metadata": {}, | ||
| 226 | "source": [ | ||
| 227 | "# Fast exponentiation\n", | ||
| 228 | "\n", | ||
| 229 | "The following cell contains two functions for computing $a^n$ ($n$ non-negative integer): a slow one that runs in $O(n)$ and a fast one that runs in $O(\\log_2(n))$. We compare these two also with Python's built-in operator `**`.\n", | ||
| 230 | "\n", | ||
| 231 | "Complexity: $O(n)$ for the slow algorithm, $O(\\log_2(n))$ for the other two." | ||
| 232 | ] | ||
| 233 | }, | ||
| 234 | { | ||
| 235 | "cell_type": "code", | ||
| 236 | "execution_count": 30, | ||
| 237 | "metadata": {}, | ||
| 238 | "outputs": [ | ||
| 239 | { | ||
| 240 | "name": "stdout", | ||
| 241 | "output_type": "stream", | ||
| 242 | "text": [ | ||
| 243 | "2.71828179834636\n", | ||
| 244 | "2.7182817863957984\n", | ||
| 245 | "2.7182817983473577\n", | ||
| 246 | "Time for slow_power(): 3.234879998000004\n", | ||
| 247 | "Time for fast_power(): 9.059099999575437e-05\n", | ||
| 248 | "Time for Python's **: 0.00010159500000384014\n" | ||
| 249 | ] | ||
| 250 | } | ||
| 251 | ], | ||
| 252 | "source": [ | ||
| 253 | "import time\n", | ||
| 254 | "\n", | ||
| 255 | "def slow_power(a, n):\n", | ||
| 256 | " r = 1\n", | ||
| 257 | " for i in range(n):\n", | ||
| 258 | " r = r * a\n", | ||
| 259 | " return r\n", | ||
| 260 | "\n", | ||
| 261 | "def fast_power(a, n):\n", | ||
| 262 | " if n == 0:\n", | ||
| 263 | " return 1\n", | ||
| 264 | " if n%2 == 0:\n", | ||
| 265 | " return fast_power(a*a, n//2)\n", | ||
| 266 | " else:\n", | ||
| 267 | " return a * fast_power(a, n-1)\n", | ||
| 268 | "\n", | ||
| 269 | "a = 1.00000001\n", | ||
| 270 | "n = 100000000\n", | ||
| 271 | "\n", | ||
| 272 | "t0 = time.process_time()\n", | ||
| 273 | "print(slow_power(a, n))\n", | ||
| 274 | "t1 = time.process_time()\n", | ||
| 275 | "print(fast_power(a, n))\n", | ||
| 276 | "t2 = time.process_time()\n", | ||
| 277 | "print(a**n)\n", | ||
| 278 | "t3 = time.process_time()\n", | ||
| 279 | "\n", | ||
| 280 | "print(\"Time for slow_power():\", t1-t0)\n", | ||
| 281 | "print(\"Time for fast_power():\", t2-t1)\n", | ||
| 282 | "print(\"Time for Python's **: \", t3-t2)" | ||
| 283 | ] | ||
| 284 | }, | ||
| 285 | { | ||
| 286 | "cell_type": "markdown", | ||
| 287 | "metadata": {}, | ||
| 288 | "source": [ | ||
| 289 | "# Fast gcd\n", | ||
| 290 | "\n", | ||
| 291 | "Complexity: $O(\\log_2(n))$" | ||
| 292 | ] | ||
| 293 | }, | ||
| 294 | { | ||
| 295 | "cell_type": "code", | ||
| 296 | "execution_count": 31, | ||
| 297 | "metadata": {}, | ||
| 298 | "outputs": [ | ||
| 299 | { | ||
| 300 | "name": "stdout", | ||
| 301 | "output_type": "stream", | ||
| 302 | "text": [ | ||
| 303 | "126\n", | ||
| 304 | "Running time: 0.00017707599999994272\n" | ||
| 305 | ] | ||
| 306 | } | ||
| 307 | ], | ||
| 308 | "source": [ | ||
| 309 | "import time\n", | ||
| 310 | "\n", | ||
| 311 | "def gcd(a, b):\n", | ||
| 312 | " if b == 0:\n", | ||
| 313 | " return a\n", | ||
| 314 | " else:\n", | ||
| 315 | " return gcd(b, a%b)\n", | ||
| 316 | "\n", | ||
| 317 | "t0 = time.process_time()\n", | ||
| 318 | "print(gcd(155275387236018, 572335397352432))\n", | ||
| 319 | "t1 = time.process_time()\n", | ||
| 320 | "\n", | ||
| 321 | "print(\"Running time:\", t1-t0)" | ||
| 322 | ] | ||
| 323 | }, | ||
| 324 | { | ||
| 325 | "cell_type": "markdown", | ||
| 326 | "metadata": {}, | ||
| 327 | "source": [ | ||
| 328 | "# Fibonacci numbers\n", | ||
| 329 | "\n", | ||
| 330 | "In the following cell there are two functions that compute the $n$-th Fibonacci number. They are almost the same, but the second one memorizes the results in a list to avoid computing them multiple times, and it is much much faster.\n", | ||
| 331 | "\n", | ||
| 332 | "Complexity: $O\\left(\\left(\\frac{1+\\sqrt 5}{2}\\right)^n\\right)\\sim O(1.6^n)$ for the slow version, $O(n)$ for the fast version." | ||
| 333 | ] | ||
| 334 | }, | ||
| 335 | { | ||
| 336 | "cell_type": "code", | ||
| 337 | "execution_count": 37, | ||
| 338 | "metadata": {}, | ||
| 339 | "outputs": [ | ||
| 340 | { | ||
| 341 | "name": "stdout", | ||
| 342 | "output_type": "stream", | ||
| 343 | "text": [ | ||
| 344 | "9227465\n", | ||
| 345 | "9227465\n", | ||
| 346 | "Time for F_slow: 2.3301570169999906\n", | ||
| 347 | "Time for F_fast: 8.848800000293977e-05\n" | ||
| 348 | ] | ||
| 349 | } | ||
| 350 | ], | ||
| 351 | "source": [ | ||
| 352 | "import time\n", | ||
| 353 | "\n", | ||
| 354 | "F_memorized = [-1] * (10**6)\n", | ||
| 355 | "\n", | ||
| 356 | "def F_slow(n):\n", | ||
| 357 | " if n <= 1:\n", | ||
| 358 | " return n\n", | ||
| 359 | " else:\n", | ||
| 360 | " return F_slow(n-1) + F_slow(n-2)\n", | ||
| 361 | " \n", | ||
| 362 | "def F_fast(n):\n", | ||
| 363 | " if F_memorized[n] == -1:\n", | ||
| 364 | " if n <= 1:\n", | ||
| 365 | " F_memorized[n] = n\n", | ||
| 366 | " else:\n", | ||
| 367 | " F_memorized[n] = F_fast(n-1) + F_fast(n-2)\n", | ||
| 368 | " \n", | ||
| 369 | " return F_memorized[n]\n", | ||
| 370 | "\n", | ||
| 371 | "n = 35\n", | ||
| 372 | "\n", | ||
| 373 | "t0 = time.process_time()\n", | ||
| 374 | "print(F_slow(n))\n", | ||
| 375 | "t1 = time.process_time()\n", | ||
| 376 | "print(F_fast(n))\n", | ||
| 377 | "t2 = time.process_time()\n", | ||
| 378 | "\n", | ||
| 379 | "print(\"Time for F_slow:\", t1-t0)\n", | ||
| 380 | "print(\"Time for F_fast:\", t2-t1)" | ||
| 381 | ] | ||
| 382 | }, | ||
| 383 | { | ||
| 384 | "cell_type": "code", | ||
| 385 | "execution_count": null, | ||
| 386 | "metadata": {}, | ||
| 387 | "outputs": [], | ||
| 388 | "source": [] | ||
| 389 | } | ||
| 390 | ], | ||
| 391 | "metadata": { | ||
| 392 | "kernelspec": { | ||
| 393 | "display_name": "Python 3", | ||
| 394 | "language": "python", | ||
| 395 | "name": "python3" | ||
| 396 | }, | ||
| 397 | "language_info": { | ||
| 398 | "codemirror_mode": { | ||
| 399 | "name": "ipython", | ||
| 400 | "version": 3 | ||
| 401 | }, | ||
| 402 | "file_extension": ".py", | ||
| 403 | "mimetype": "text/x-python", | ||
| 404 | "name": "python", | ||
| 405 | "nbconvert_exporter": "python", | ||
| 406 | "pygments_lexer": "ipython3", | ||
| 407 | "version": "3.8.5" | ||
| 408 | } | ||
| 409 | }, | ||
| 410 | "nbformat": 4, | ||
| 411 | "nbformat_minor": 4 | ||
| 412 | } | ||
diff --git a/src/Lecture7/slides/.ipynb_checkpoints/X2-StudentsRequests-checkpoint.ipynb b/src/Lecture7/slides/.ipynb_checkpoints/X2-StudentsRequests-checkpoint.ipynb new file mode 100644 index 0000000..a35bb8d --- /dev/null +++ b/src/Lecture7/slides/.ipynb_checkpoints/X2-StudentsRequests-checkpoint.ipynb | |||
| @@ -0,0 +1,165 @@ | |||
| 1 | { | ||
| 2 | "cells": [ | ||
| 3 | { | ||
| 4 | "cell_type": "markdown", | ||
| 5 | "metadata": {}, | ||
| 6 | "source": [ | ||
| 7 | "# Diffie-Hellman key exchange\n", | ||
| 8 | "\n", | ||
| 9 | "The following is a simple implementation of the classic [Diffie-Hellman key exchange](https://en.wikipedia.org/wiki/Diffie%E2%80%93Hellman_key_exchange) cryptographic protocol." | ||
| 10 | ] | ||
| 11 | }, | ||
| 12 | { | ||
| 13 | "cell_type": "code", | ||
| 14 | "execution_count": 9, | ||
| 15 | "metadata": {}, | ||
| 16 | "outputs": [ | ||
| 17 | { | ||
| 18 | "name": "stdout", | ||
| 19 | "output_type": "stream", | ||
| 20 | "text": [ | ||
| 21 | "Public key: p = 20747 and g = 13428 \n", | ||
| 22 | "\n", | ||
| 23 | "[[ Alice's secret key: a = 12403 ]]\n", | ||
| 24 | "[[ Bob's secret key: b = 17642 ]] \n", | ||
| 25 | "\n", | ||
| 26 | "Alice sends h1 = 14710 to Bob\n", | ||
| 27 | "Bob sends h2 = 10680 to Alice \n", | ||
| 28 | "\n", | ||
| 29 | "Alice computed 10455 using h2 and her secret a\n", | ||
| 30 | "Bob computed 10455 using h1 and his secret b\n" | ||
| 31 | ] | ||
| 32 | } | ||
| 33 | ], | ||
| 34 | "source": [ | ||
| 35 | "# Public information:\n", | ||
| 36 | "p = Primes()[10^3 + randint(1,10000)] # random prime\n", | ||
| 37 | "g = randint(2, p-1) # random integer\n", | ||
| 38 | "\n", | ||
| 39 | "print(\"Public key: p =\", p, \"and g =\", g, \"\\n\")\n", | ||
| 40 | "\n", | ||
| 41 | "a = randint(2, p-1) # Only Alice knows this\n", | ||
| 42 | "b = randint(2, p-1) # Only Bob knows this\n", | ||
| 43 | "\n", | ||
| 44 | "print(\"[[ Alice's secret key: a =\", a, \"]]\")\n", | ||
| 45 | "print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n", | ||
| 46 | "\n", | ||
| 47 | "h1 = (g^a) % p # Alice sends this to Bob\n", | ||
| 48 | "h2 = (g^b) % p # Bob sends this to Alice\n", | ||
| 49 | "\n", | ||
| 50 | "print(\"Alice sends h1 =\", h1, \"to Bob\")\n", | ||
| 51 | "print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n", | ||
| 52 | "\n", | ||
| 53 | "secret_a = (h2^a) % p # Alice can compute this because she knows a\n", | ||
| 54 | "secret_b = (h1^b) % p # Bob can compute this because he knows b\n", | ||
| 55 | "\n", | ||
| 56 | "print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n", | ||
| 57 | "print(\"Bob computed\", secret_b, \"using h1 and his secret b\")" | ||
| 58 | ] | ||
| 59 | }, | ||
| 60 | { | ||
| 61 | "cell_type": "markdown", | ||
| 62 | "metadata": {}, | ||
| 63 | "source": [ | ||
| 64 | "## General Diffie-Hellman\n", | ||
| 65 | "\n", | ||
| 66 | "The following code is an implementation of a generic Diffie-Hellman key exchange protocol that uses a group $G$ instead of $(\\mathbb Z/p \\mathbb Z)^\\times$." | ||
| 67 | ] | ||
| 68 | }, | ||
| 69 | { | ||
| 70 | "cell_type": "code", | ||
| 71 | "execution_count": 24, | ||
| 72 | "metadata": {}, | ||
| 73 | "outputs": [ | ||
| 74 | { | ||
| 75 | "name": "stdout", | ||
| 76 | "output_type": "stream", | ||
| 77 | "text": [ | ||
| 78 | "Public key:\n", | ||
| 79 | "G = Additive abelian group isomorphic to Z/171 embedded in Abelian group of points on Elliptic Curve defined by y^2 = x^3 + x + 156 over Finite Field of size 157 \n", | ||
| 80 | "g = (155 : 60 : 1) \n", | ||
| 81 | "\n", | ||
| 82 | "[[ Alice's secret key: a = 141 ]]\n", | ||
| 83 | "[[ Bob's secret key: b = 158 ]] \n", | ||
| 84 | "\n", | ||
| 85 | "Alice sends h1 = (29 : 125 : 1) to Bob\n", | ||
| 86 | "Bob sends h2 = (60 : 59 : 1) to Alice \n", | ||
| 87 | "\n", | ||
| 88 | "Alice computed (109 : 94 : 1) using h2 and her secret a\n", | ||
| 89 | "Bob computed (109 : 94 : 1) using h1 and his secret b\n" | ||
| 90 | ] | ||
| 91 | } | ||
| 92 | ], | ||
| 93 | "source": [ | ||
| 94 | "def genericDH(G):\n", | ||
| 95 | " if G.cardinality() == 1:\n", | ||
| 96 | " print(\"Group is trivial, can't do anything\")\n", | ||
| 97 | " return\n", | ||
| 98 | " g = G.random_element()\n", | ||
| 99 | " while g == G.identity(): # Make sure g is not trivial\n", | ||
| 100 | " g = G.random_element()\n", | ||
| 101 | " \n", | ||
| 102 | " print(\"Public key:\\nG =\", G, \"\\ng =\", g, \"\\n\")\n", | ||
| 103 | " \n", | ||
| 104 | " a = randint(2, G.exponent()-1) # Only Alice knows this\n", | ||
| 105 | " b = randint(2, G.exponent()-1) # Only Bob knows this\n", | ||
| 106 | "\n", | ||
| 107 | " print(\"[[ Alice's secret key: a =\", a, \"]]\")\n", | ||
| 108 | " print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n", | ||
| 109 | " \n", | ||
| 110 | " # \"Ternary operator\", I did not explain this\n", | ||
| 111 | " # https://docs.python.org/3/reference/expressions.html#conditional-expressions\n", | ||
| 112 | " h1 = g^a if G.is_multiplicative() else a*g # Alice sends this to Bob\n", | ||
| 113 | " h2 = g^b if G.is_multiplicative() else b*g # Bob sends this to Alice\n", | ||
| 114 | "\n", | ||
| 115 | " print(\"Alice sends h1 =\", h1, \"to Bob\")\n", | ||
| 116 | " print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n", | ||
| 117 | " \n", | ||
| 118 | " secret_a = h2^a if G.is_multiplicative() else a*h2 # Alice can compute this because she knows a\n", | ||
| 119 | " secret_b = h1^b if G.is_multiplicative() else b*h1 # Bob can compute this because he knows b\n", | ||
| 120 | "\n", | ||
| 121 | " print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n", | ||
| 122 | " print(\"Bob computed\", secret_b, \"using h1 and his secret b\")\n", | ||
| 123 | " \n", | ||
| 124 | "E = EllipticCurve(GF(157), [1,-1])\n", | ||
| 125 | "G = E.abelian_group()\n", | ||
| 126 | "genericDH(G)" | ||
| 127 | ] | ||
| 128 | }, | ||
| 129 | { | ||
| 130 | "cell_type": "markdown", | ||
| 131 | "metadata": {}, | ||
| 132 | "source": [ | ||
| 133 | "# Numerical methods for PDEs" | ||
| 134 | ] | ||
| 135 | }, | ||
| 136 | { | ||
| 137 | "cell_type": "code", | ||
| 138 | "execution_count": null, | ||
| 139 | "metadata": {}, | ||
| 140 | "outputs": [], | ||
| 141 | "source": [] | ||
| 142 | } | ||
| 143 | ], | ||
| 144 | "metadata": { | ||
| 145 | "kernelspec": { | ||
| 146 | "display_name": "SageMath 9.0", | ||
| 147 | "language": "sage", | ||
| 148 | "name": "sagemath" | ||
| 149 | }, | ||
| 150 | "language_info": { | ||
| 151 | "codemirror_mode": { | ||
| 152 | "name": "ipython", | ||
| 153 | "version": 3 | ||
| 154 | }, | ||
| 155 | "file_extension": ".py", | ||
| 156 | "mimetype": "text/x-python", | ||
| 157 | "name": "python", | ||
| 158 | "nbconvert_exporter": "python", | ||
| 159 | "pygments_lexer": "ipython3", | ||
| 160 | "version": "3.8.5" | ||
| 161 | } | ||
| 162 | }, | ||
| 163 | "nbformat": 4, | ||
| 164 | "nbformat_minor": 4 | ||
| 165 | } | ||
diff --git a/src/Lecture7/slides/.ipynb_checkpoints/X2-StudentsRequests-notebook-checkpoint.ipynb b/src/Lecture7/slides/.ipynb_checkpoints/X2-StudentsRequests-notebook-checkpoint.ipynb new file mode 100644 index 0000000..ccd6a21 --- /dev/null +++ b/src/Lecture7/slides/.ipynb_checkpoints/X2-StudentsRequests-notebook-checkpoint.ipynb | |||
| @@ -0,0 +1,298 @@ | |||
| 1 | { | ||
| 2 | "cells": [ | ||
| 3 | { | ||
| 4 | "cell_type": "markdown", | ||
| 5 | "metadata": {}, | ||
| 6 | "source": [ | ||
| 7 | "# Diffie-Hellman key exchange\n", | ||
| 8 | "\n", | ||
| 9 | "The following is a simple implementation of the classic [Diffie-Hellman key exchange](https://en.wikipedia.org/wiki/Diffie%E2%80%93Hellman_key_exchange) cryptographic protocol." | ||
| 10 | ] | ||
| 11 | }, | ||
| 12 | { | ||
| 13 | "cell_type": "code", | ||
| 14 | "execution_count": 1, | ||
| 15 | "metadata": {}, | ||
| 16 | "outputs": [ | ||
| 17 | { | ||
| 18 | "name": "stdout", | ||
| 19 | "output_type": "stream", | ||
| 20 | "text": [ | ||
| 21 | "Public key: p = 75521 and g = 58258 \n", | ||
| 22 | "\n", | ||
| 23 | "[[ Alice's secret key: a = 22794 ]]\n", | ||
| 24 | "[[ Bob's secret key: b = 69773 ]] \n", | ||
| 25 | "\n", | ||
| 26 | "Alice sends h1 = 31067 to Bob\n", | ||
| 27 | "Bob sends h2 = 54398 to Alice \n", | ||
| 28 | "\n", | ||
| 29 | "Alice computed 30031 using h2 and her secret a\n", | ||
| 30 | "Bob computed 30031 using h1 and his secret b\n" | ||
| 31 | ] | ||
| 32 | } | ||
| 33 | ], | ||
| 34 | "source": [ | ||
| 35 | "# Public information:\n", | ||
| 36 | "p = Primes()[10^3 + randint(1,10000)] # random prime\n", | ||
| 37 | "g = randint(2, p-1) # random integer\n", | ||
| 38 | "\n", | ||
| 39 | "print(\"Public key: p =\", p, \"and g =\", g, \"\\n\")\n", | ||
| 40 | "\n", | ||
| 41 | "a = randint(2, p-1) # Only Alice knows this\n", | ||
| 42 | "b = randint(2, p-1) # Only Bob knows this\n", | ||
| 43 | "\n", | ||
| 44 | "print(\"[[ Alice's secret key: a =\", a, \"]]\")\n", | ||
| 45 | "print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n", | ||
| 46 | "\n", | ||
| 47 | "h1 = (g^a) % p # Alice sends this to Bob\n", | ||
| 48 | "h2 = (g^b) % p # Bob sends this to Alice\n", | ||
| 49 | "\n", | ||
| 50 | "print(\"Alice sends h1 =\", h1, \"to Bob\")\n", | ||
| 51 | "print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n", | ||
| 52 | "\n", | ||
| 53 | "secret_a = (h2^a) % p # Alice can compute this because she knows a\n", | ||
| 54 | "secret_b = (h1^b) % p # Bob can compute this because he knows b\n", | ||
| 55 | "\n", | ||
| 56 | "print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n", | ||
| 57 | "print(\"Bob computed\", secret_b, \"using h1 and his secret b\")" | ||
| 58 | ] | ||
| 59 | }, | ||
| 60 | { | ||
| 61 | "cell_type": "markdown", | ||
| 62 | "metadata": {}, | ||
| 63 | "source": [ | ||
| 64 | "## General Diffie-Hellman\n", | ||
| 65 | "\n", | ||
| 66 | "The following code is an implementation of a generic Diffie-Hellman key exchange protocol that uses a group $G$ instead of $(\\mathbb Z/p \\mathbb Z)^\\times$." | ||
| 67 | ] | ||
| 68 | }, | ||
| 69 | { | ||
| 70 | "cell_type": "code", | ||
| 71 | "execution_count": 2, | ||
| 72 | "metadata": {}, | ||
| 73 | "outputs": [ | ||
| 74 | { | ||
| 75 | "name": "stdout", | ||
| 76 | "output_type": "stream", | ||
| 77 | "text": [ | ||
| 78 | "Public key:\n", | ||
| 79 | "G = Additive abelian group isomorphic to Z/171 embedded in Abelian group of points on Elliptic Curve defined by y^2 = x^3 + x + 156 over Finite Field of size 157 \n", | ||
| 80 | "g = (53 : 90 : 1) \n", | ||
| 81 | "\n", | ||
| 82 | "[[ Alice's secret key: a = 145 ]]\n", | ||
| 83 | "[[ Bob's secret key: b = 65 ]] \n", | ||
| 84 | "\n", | ||
| 85 | "Alice sends h1 = (150 : 80 : 1) to Bob\n", | ||
| 86 | "Bob sends h2 = (4 : 58 : 1) to Alice \n", | ||
| 87 | "\n", | ||
| 88 | "Alice computed (28 : 28 : 1) using h2 and her secret a\n", | ||
| 89 | "Bob computed (28 : 28 : 1) using h1 and his secret b\n" | ||
| 90 | ] | ||
| 91 | } | ||
| 92 | ], | ||
| 93 | "source": [ | ||
| 94 | "def genericDH(G):\n", | ||
| 95 | " if G.cardinality() == 1:\n", | ||
| 96 | " print(\"Group is trivial, can't do anything\")\n", | ||
| 97 | " return\n", | ||
| 98 | " g = G.random_element()\n", | ||
| 99 | " while g == G.identity(): # Make sure g is not trivial\n", | ||
| 100 | " g = G.random_element()\n", | ||
| 101 | " \n", | ||
| 102 | " print(\"Public key:\\nG =\", G, \"\\ng =\", g, \"\\n\")\n", | ||
| 103 | " \n", | ||
| 104 | " a = randint(2, G.exponent()-1) # Only Alice knows this\n", | ||
| 105 | " b = randint(2, G.exponent()-1) # Only Bob knows this\n", | ||
| 106 | "\n", | ||
| 107 | " print(\"[[ Alice's secret key: a =\", a, \"]]\")\n", | ||
| 108 | " print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n", | ||
| 109 | " \n", | ||
| 110 | " # \"Ternary operator\", I did not explain this\n", | ||
| 111 | " # https://docs.python.org/3/reference/expressions.html#conditional-expressions\n", | ||
| 112 | " h1 = g^a if G.is_multiplicative() else a*g # Alice sends this to Bob\n", | ||
| 113 | " h2 = g^b if G.is_multiplicative() else b*g # Bob sends this to Alice\n", | ||
| 114 | "\n", | ||
| 115 | " print(\"Alice sends h1 =\", h1, \"to Bob\")\n", | ||
| 116 | " print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n", | ||
| 117 | " \n", | ||
| 118 | " secret_a = h2^a if G.is_multiplicative() else a*h2 # Alice can compute this because she knows a\n", | ||
| 119 | " secret_b = h1^b if G.is_multiplicative() else b*h1 # Bob can compute this because he knows b\n", | ||
| 120 | "\n", | ||
| 121 | " print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n", | ||
| 122 | " print(\"Bob computed\", secret_b, \"using h1 and his secret b\")\n", | ||
| 123 | " \n", | ||
| 124 | "E = EllipticCurve(GF(157), [1,-1])\n", | ||
| 125 | "G = E.abelian_group()\n", | ||
| 126 | "genericDH(G)" | ||
| 127 | ] | ||
| 128 | }, | ||
| 129 | { | ||
| 130 | "cell_type": "markdown", | ||
| 131 | "metadata": {}, | ||
| 132 | "source": [ | ||
| 133 | "# Numerical methods for differential equations\n", | ||
| 134 | "\n", | ||
| 135 | "## Euler's method (ODE)\n", | ||
| 136 | "\n", | ||
| 137 | "In sage you can use [`ode_solver()`](https://doc.sagemath.org/html/en/reference/calculus/sage/calculus/ode.html) to solve any ordinary differential equation by hand, but Euler's method is very simple to implement by hand:" | ||
| 138 | ] | ||
| 139 | }, | ||
| 140 | { | ||
| 141 | "cell_type": "code", | ||
| 142 | "execution_count": 4, | ||
| 143 | "metadata": {}, | ||
| 144 | "outputs": [ | ||
| 145 | { | ||
| 146 | "data": { | ||
| 147 | "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5yOdf7H8fcYzFCMs3EYzCLHQowIlYTVQee2ItRurUUl20ltRduGTmtpapeKWp1+5dhBpeQUyjGHEFKImRznhBlzz/X747NjHIa5cc993fd9vZ6Px/24m9M1n7lNeff9fq7PN8pxHEcAAAAoUgm3CwAAAAgXBCcAAAA/EZwAAAD8RHACAADwE8EJAADATwQnAAAAPxGcAAAA/ERwAgAA8BPBCfCT4zhKT08XM2MBwLsIToCfMjIyFBcXp4yMDLdLAQC4hOAEAADgJ4ITAACAnwhOAAAAfiI4AQAA+IngBAAAIs6XX0qzZwf+ugQnoAjJyclq2rSpkpKS3C4FAOCHnBypf3/phRcCf+0oh6E0gF/S09MVFxentLQ0lS9f3u1yAAAnMXasNHiw9P33UvPmgb02K04AACBipKVJTz8t9esX+NAkEZwAAEAEGTVKysqy8FQcCE4AACAibN8u/fOf0pAhUq1axfM9CE4AACAiPPmkVK6c9PDDxfc9ShbfpQEAAIJj1Spp4kRrDC/O+3e4qw7wE3fVAUBochypa1dp2zZpzRqpVKni+16sOAEAgLD26afSV19J06cXb2iS6HECisQATAAIXYcPSw8+KF1+uXTNNcX//diqA/zEVh0AhJ6XX5buu09avlxq2bL4vx8rTgAAICzt2ycNGybdeWdwQpNEcAIAAGHqH/+QDh2SnnkmeN+T4AQAAMLO5s3SmDHSo49KNWoE7/sSnOAJ9erVU1RU1AmPgQMHul0aAOAMPPKIVL26TQkPJsYRwBOWLFkin8935O01a9aoa9euuvnmm12sCgBwJubNkyZPliZNksqWDe735q46eNLgwYP18ccfa+PGjYqKivLra7irDgDcl5cntW0rlSghLV5sz8HEihM8JycnR5MmTdKQIUNOGZqys7OVnZ195O309PRglAcAOIW335aWLZPmzw9+aJLocYIHTZs2Tfv371e/fv1O+XkjRoxQXFzckUdCQkJwCgQAFCojw3qbbr5Z6tjRnRrYqoPndO/eXaVLl9ZHH310ys8rbMUpISGBrToAcMnQodLo0dL69VLduu7UwFYdPOWXX37Rl19+qSlTphT5uTExMYqJiQlCVQCAomzcKL30kvTYY+6FJomtOnjMhAkTVK1aNV111VVulwIAOA1Dhti8pocfdrcOVpzgGXl5eZowYYL69u2rkiX51QeAcDFzpvTxx9KHH0plyrhbCytO8Iwvv/xSW7du1V133eV2KQAAP+XkSIMHS5dfLt1wg9vVsOIED+nWrZu4FwIAwsuYMXa8yuTJkp9j94oVK04AACAk7dwpDR8uDRggNW/udjWG4AQAAELS0KFSTIyFp1BBcAKKkJycrKZNmyopKcntUgDAM779VnrzTenZZ6WKFd2upgADMAE/cVYdAARHXp7Urp2UmystWSJFR7tdUQGawwEAQEh5800LTPPnh1ZoktiqAwAAISQtTXr0Uen22907j+5UCE4AACBkPPGElJUljRrldiWFY6sOAACEhBUrpORk6bnnpNq13a6mcDSHA36iORwAik9ennTxxVJmpgWoUqXcrqhwrDgBAADXvfGGjSCYNy90Q5NEjxMAAHDZ7t3SI49IffpInTq5Xc2pEZyAIjAAEwCK19Chks9nvU2hjh4nwE/0OAFA4C1aZL1Nycl2Jl2oIzgBfiI4AUBg5eZKSUk25PLbb0Nv2GVhaA4HAACueOUV6fvvpcWLwyM0SfQ4AQAAF+zYYcMu77lHatvW7Wr8R3ACAABBd999Umys9OyzbldyetiqAwAAQfXRR9LkydK770qVKrldzemhORzwE83hAHD2MjKkZs3s8emnUlSU2xWdHrbqAABA0DzxhLRnjzWGh1tokghOQJEYgAkAgbFkiTR2rDR8uJSY6HY1Z4atOsBPbNUBwJnLn9kkWYAqGaZd1mFaNgAACCejR0urVtnMpnANTRJbdQAAoJht2SI99ZR0770Fq07hiuAEAACKjePYGXSVK0t//7vb1Zy9MF4sAwAAoe7996XPPpNmzJDKlXO7mrPHihMAACgW+/ZJ998v3XSTdM01blcTGAQnAABQLB5+WDp0SPrXv9yuJHDYqgMAAAH31VfSa69Jr74q1azpdjWBw4oTUAQGYALA6cnMlO6+W7rsMumee9yuJrAYgAn4iQGYAOCf+++Xxo+XVq+W6td3u5rAYqsOAAAEzIIFdqzKiy9GXmiSWHEC/MaKEwCc2sGDUsuWUqVKFqCio92uKPBYcQIAAAExfLj088/StGmRGZokmsPhIb/++qt69+6typUrq2zZsmrZsqWWLVvmdlkAEBGWLJGef14aNkxq0sTtaooPK07whH379qlDhw7q3LmzZs6cqWrVqmnz5s2qUKGC26UBQNjLyZHuuktq0UJ68EG3qyleBCd4wqhRo5SQkKAJEyYceV+9evXcKwgAIsizz0rr19uqU6lSbldTvNiqgyfMmDFDbdq00c0336xq1aqpVatWGj9+vNtlAUDYW7VK+sc/pKFDrTE80nFXHTwhNjZWkjRkyBDdfPPN+u677zR48GD95z//UZ8+fQr9muzsbGVnZx95Oz09XQkJCdxVBwD/k5MjtW0r+XzS0qVSTIzbFRU/ghM8oXTp0mrTpo0WLlx45H333XeflixZokWLFhX6NcOGDdPw4cNPeD/BCQDME09II0dK330ntWrldjXBwVYdPKFGjRpq2rTpMe9r0qSJtm7detKvGTp0qNLS0o48tm3bVtxlAkDY+O47acQI6cknvROaJJrD4REdOnTQhg0bjnnfjz/+qLp16570a2JiYhTjhXVnADhNBw9KfftaYHr0UberCS6CEzzhgQce0MUXX6xnn31Wt9xyi7777juNGzdO48aNc7s0AAg7f/ubtGWLtHx55N9Fdzx6nOAZH3/8sYYOHaqNGzcqMTFRQ4YM0d133+3313PkCgBI8+dLl14qPfdc5M9sKgzBCfATwQmA12Vm2pDLGjWkuXMj91iVU2GrDgAA+OXhh6WUFOnzz70ZmiSCEwAA8MOsWdKrr0rJyVKDBm5X4x626gA/sVUHwKv27ZMuuEBq3NhWm0p4eJiRh390wD/Jyclq2rSpkpKS3C4FAILOcaT+/a2/6Y03vB2aJFacAL+x4gTAi956y2Y2vf++dMstblfjPo/nRgAAcDJbtkiDBkl9+hCa8hGcAADACXJzpTvukCpXlsaOdbua0MFddQAA4AQjR0qLFknz5kl0JxRgxQkAABzju++kYcOkxx6TOnRwu5rQQnM44CeawwF4QWamHd5bsaL0zTfeO4uuKGzVAQCAIx54QNqxQ/rkE0JTYQhOAABAkjR1qvTaa9K4cdJ557ldTWiixwkoAgMwAXjBzp3S3XdL114r/elPblcTuuhxAvxEjxOASJWXJ3XrJq1dK61aJVWt6nZFoYutOgAAPG7UKGn2bOmLLwhNRWGrDgAAD1u4UHriCWnoUOmKK9yuJvSxVQf4ia06AJFm3z6pZUupdm1pzhzuovMHK04AAHiQ41gTeHq69M47hCZ/0eMEAIAH/fvf0pQp0uTJUt26blcTPlhxAgDAY1atskGXAwZIN9zgdjXhhR4nwE/0OAGIBFlZUps2UunS0rffSrGxblcUXtiqA4qQnJys5ORk+Xw+t0sBgLN2333S1q3S0qWEpjPBihPgJ1acAIS7t9+WeveW3nhDuvNOt6sJT/Q4AQDgAT/8IN1zjwWnfv3criZ8seIE+IkVJwDhKjNTattWKlHC+prOOcftisIXPU4AAEQwx5H697e+piVLCE1ni+AEAEAEGzfOepvefVdq0sTtasIfPU4AAESoZcvsLroBA6Rbb3W7mshAjxPgJ3qcAISTffuk1q2lypWlBQukmBi3K4oMbNUBABBhHMfGDezbJ331FaEpkNiqA4qQnJyspk2bKikpye1SAMAvL74oTZ8uvfWWlJjodjWRha06wE9s1QEIB3PnSl26SA8+KI0c6XY1kYfgBPiJ4AQg1G3fbn1NzZtLn38ulaQhJ+DYqgMAIAJkZ0s33mjnz733HqGpuPCyAgAQAe69V/r+e7uDrmpVt6uJXKw4wROGDRumqKioYx7x8fFulwUAATF+vD1efVVq08btaiIbK07wjGbNmunLL7888nZ0dLSL1QBAYHz7rTRokPSXv9gIAhQvghM8o2TJkqwyAYgoqanW19S6tTR6tNvVeANbdfCMjRs3qmbNmkpMTNStt96qn3766ZSfn52drfT09GMeABAqDh+WbrlF8vmkDz+USpd2uyJvIDjBEy666CK99dZb+vzzzzV+/HilpKTo4osv1p49e076NSNGjFBcXNyRR0JCQhArBoBTe+ghaeFC6YMPpJo13a7GO5jjBE/KyspS/fr19fDDD2vIkCGFfk52drays7OPvJ2enq6EhATmOAFw3ZtvSv36SWPHWn8TgoceJ3jSOeeco/PPP18bN2486efExMQohgOeAISYRYuke+6R/vhHaeBAt6vxHrbq4EnZ2dlat26datSo4XYpAOC37dul66+X2raVXnlFiopyuyLvITjBEx588EHNnTtXW7Zs0bfffqubbrpJ6enp6tu3r9ulAYBfDhyQrr1WiomRJk+mGdwtbNXBE7Zv367bbrtNu3fvVtWqVdWuXTstXrxYdevWdbs0ACiS40h33SWtXy99841UrZrbFXkXwQme8N5777ldAgCcsWefld5/38YOtGzpdjXexlYdAAAhbPp06W9/k556yoZdwl0EJwAAQtTq1VKvXhaYnnzS7WogEZyAIiUnJ6tp06ZKSkpyuxQAHpKaKvXsKTVoYHObSvA3dkhgACbgp/T0dMXFxTEAE0CxO3hQ6txZ+uUXO8S3Th23K0I+msMBAAgheXlS377SqlXSvHmEplBDcAIAIIQ88YTdPTd5stSmjdvV4HgEJwAAQsTEiTZ64LnnbEI4Qg+tZgAAhIA5c+wMurvvlh580O1qcDIEJwAAXLZhg3TDDdKll0rJyZxBF8oITgAAuGj3bumqq6QaNaQPPpBKlXK7IpwKPU4AALjk0CHrZUpPt7EDFSq4XRGKwooTUAQGYAIoDnl5Up8+0tKl0owZUmKi2xXBHwzABPzEAEwAgfTAA9KYMTZ24Lrr3K4G/mKrDgCAIHvpJWn0aGsEJzSFF7bqAAAIovffl/76V+nRR6UBA9yuBqeL4AQAQJDMmWN9Tb1726BLhB+CEwAAQbBmjW3Ldeokvf46s5rCFcEJAIBitn271KOHVLeuNYOXLu12RThTBCcAAIrRnj1St25SiRLSp59KcXFuV4SzwV11AAAUk8xMmwq+a5e0YIFUq5bbFeFsseIEFIEBmADORE6OdOON0tq10mefSY0auV0RAoEBmICfGIAJwF8+n9SrlzR1qjRzpnT55W5XhEBhqw4AgAByHOm+++zA3g8+IDRFGoITAAABNGyY9Mor0vjx0g03uF0NAo0eJwAAAmTsWOnpp6URI6Q//cntalAcCE4AAATApEm2RffXv0qPPOJ2NSguBCcAAM7SlClSv37SnXdKzz3HVPBIRnACAOAsfPqpdOut0s03W19TCf5mjWj88QIAcIZmz7YG8CuvlN56S4qOdrsiFDeCE1AEBmACKMw330g9e0qXXiq9/75UqpTbFSEYGIAJ+IkBmADyLV0qdekitWxpAy7LlnW7IgQLK04AAJyG1aul7t2lJk2kjz8mNHkNwQkAAD9t2CB17SrVqWPnz5Ur53ZFCDaCEwAAftiyxbbnKleWvvhCqlDB7YrgBoITPGnEiBGKiorS4MGD3S4FQBj4+Wepc2epTBnpyy+lqlXdrghuITjBc5YsWaJx48bpggsucLsUAGHg55+lyy6TSpa08QM1arhdEdxEcIKnZGZmqlevXho/frwqVqzodjkAQtyWLTZuoGRJac4cKSHB7YrgNoITPGXgwIG66qqrdMUVV7hdCoAQt2WLrTSVLm2hqXZttytCKCjpdgFAsLz33ntavny5lixZ4tfnZ2dnKzs7+8jb6enpxVUagBDz00/W01S6tPT114QmFGDFCZ6wbds23X///Zo0aZJiY2P9+poRI0YoLi7uyCOBNXrAE376iZUmnByTw+EJ06ZN0/XXX6/oow6S8vl8ioqKUokSJZSdnX3Mx6TCV5wSEhKYHA5EsM2bbaUpNtZWmmrVcrsihBq26uAJXbp00erVq49535133qnGjRvrkUceOSE0SVJMTIxiYmKCVSIAlxGa4A+CEzyhXLlyat68+THvO+ecc1S5cuUT3g/AezZskK64wuY0zZkj1azpdkUIVfQ4AQA8bdUq6ZJLpPLlCU0RZcYM6YEH7DmA6HEC/JSenq64uDh6nIAI8t130u9/L9WrZ8eoVKnidkUIiGnTpOuvl0qUkPLypOnTpZ49A3JpVpwAAJ40b56dPdekiU0EJzRFiCVLpP797Z/z8qToaFtKDBCCEwDAcz7/3Faa2ra1f+bA3giwZ48Fposusg5/yUKTz2fzJQKE4AQA8JSpU6VrrrHVpk8+kc491+2KcFby8qTXXpMaNZLefVcaPVratMm25+67L6DbdBI9TkCRkpOTlZycLJ/Ppx9//JEeJyCMTZok9esn3Xij/XOpUm5XhLOybJk0YIA1q/XpI40aJcXHF+u3JDgBfqI5HAhv48bZTk6/ftL48baLgzC1d6/0+OPSf/4jNW8uJSdLnToF5VuzVQcAiHgvvCD9+c/SoEG2q0NoClN5edLrr0vnnSe98470z39Ky5cHLTRJBCcAQATLy5MefFB66CFboPjXv+wOdYSh5culiy+W/vQn6corbWrp/fdLJYM7y5tfHwBARDp82LblXnxRGjNGeuYZKSrK7apw2vbutT6mNm2krCxp7lzprbeKvZfpZDhyBQAQcbKypFtukWbNshutbr3V7Ypw2vLypIkTpUcekbKzpZdekgYOdL2jn+AEAIgoe/ZIV18trV5t4wa6dnW7Ipy25cstJC1eLPXqJT3/vFSjhttVSWKrDgAQQbZulTp2tDE+c+YQmsLOvn3WwZ+UJGVk2B/ipEkhE5okVpwAABHihx+k7t2tV/ibb+zGK4SJvDzpzTdtW+7QIbsNctAg17flCsOKE1CE5ORkNW3aVElJSW6XAuAkFi2ylaZKlQhNYWfFCvvDu+suWyJcv1564IGQDE0SAzABvzEAEwhN06ZJt99uuzvTp3PuXNjYv1964gnplVekxo1tiGUAz5QrLqw4AQDC1ujR0g03WDM4h/WGify75c47z56ff15auTIsQpNEcAIAhCGfz85vfeAB6eGHpffek2Jj3a4KRVq50qZ833mndMUVNsRyyJCQ3ZYrDM3hAICwkpUl3Xab9OmndlTZPfe4XRGKtH+/9OSTth3XqJE0e7bUubPbVZ0RghMAIGzs3Cldc40tVHz8sfT737tdEU7JcaT//tfOvDlwQBo1yo5JCaMVpuMRnAAAYWHtWjuizOeTFiyQWrRwuyKc0vff20iBBQtsdPsLL0i1arld1VmjxwkAEPK+/NLOd61Y0YZJE5pCWFqarSpdeKGNcf/qKzv3JgJCk0RwAgCEuAkTpB49LDjNny/Vru12RShU/rZco0bS669LI0daM/jll7tdWUARnIAiMAATcIfPZ4Ok77pL+uMfpY8+ksqVc7sqFGrVKumSS6Q+fWyswPr11tdUurTblQUcAzABPzEAEwie9HQbajlzpvTii7bzExXldlU4QVqa9NRT0ssvSw0b2nOXLm5XVaxoDgcAhJSffrI753791UYOdO/udkU4geNIb78tPfiglJkpjRhh6TYCV5iOx1YdACBkzJkjtW0r5eRYEzihKQStXi1deql0xx32HMHbcoUhOAEAQsK4cXbGa8uW0rff2vFlCCHp6Tblu1Ur6bffpFmzpPff91y3PsEJAOCq3Fzp3nulP/9Z6t/f+poqVXK7KhyRvy3XqJGNav/HP6wZ/Ior3K7MFfQ4AQBcs3u3zUacO1f6978tPCGErFkjDRwozZsn3XyzdeonJLhdlasITgAAVyxfLt1wg53EMWuW3cWOEJGeLg0fLv3rX1L9+tIXX9g+KtiqAwAE33//K3XoIFWpIi1dSmgKGY4jvfOONZj9+9/SM8/Ythyh6QiCE1AEBmACgXP4sN213qePbdHNny/VqeN2VZBkhwF27iz16mWpdt066dFHpZgYtysLKQzABPzEAEzg7KSmSrfcIi1caDtAf/kLQy1DQkZGwbbc734njR0rdevmdlUhix4nAECxW7zYeosPH5a+/lrq2NHtiiDHsXECf/2rtG+f9PTTNm6AFaZTYqsOAFBsHMcWMjp1spuxli0jNIWEH36wo1Fuu01q186GWA4dSmjyA8EJAFAs0tJslWnwYJvTNHeuVKuW21V5XEaGTflu0ULavl367DNp8mQazU4DwQme8Oqrr+qCCy5Q+fLlVb58ebVv314zZ850uywgYn3/vdSmjY0ZmDxZeuklqVQpt6vysPxtucaNpeRk62lavZozbc4AwQmeULt2bY0cOVJLly7V0qVLdfnll+vaa6/V2rVr3S4NiCiOI73+uu3+nHuubc3dcIPbVXncunU25fvWW6WLLrK3H3uMbbkzxF118KxKlSrp+eef1x//+Ee/Pp+76oBTO3BAGjBAevNN6e67rbepTBm3q/KwzExr+P7nP6V69aQxY6QePdyuKuxxVx08x+fz6YMPPlBWVpbat29/0s/Lzs5Wdnb2kbfT09ODUR4QljZskG66Sdq82YJTnz5uV+RhjiN98IHdIbd3rzRsmN05FxvrdmURga06eMbq1at17rnnKiYmRv3799fUqVPVtGnTk37+iBEjFBcXd+SR4PHzmYCTefdd62c6fFj67jtCk6vWr7cp33/4g5SUZHfPPf44oSmA2KqDZ+Tk5Gjr1q3av3+/Jk+erNdee01z5849aXgqbMUpISGBrTrgfzIz7W65iROl22+3EzrKlXO7Ko/KzLTjUV56ye6QGzNGuvJKt6uKSAQneNYVV1yh+vXr6z//+Y9fn0+PE1Bg+XLrNd6xw27S6tOHKeCucBzpww9tW273bmv6fughVpiKEVt18CzHcY5ZUQJQtLw8W9Ro185Wl5Yvl/r2JTS5Yv16Oxrllluk1q1tW+6JJwhNxYzmcHjCY489ph49eighIUEZGRl67733NGfOHH322WdulwaEjdRUqV8/m5n4179K//gHd7S7IivLtuVefNHGsX/8sXTVVW5X5RkEJ3hCamqq7rjjDu3cuVNxcXG64IIL9Nlnn6lr165ulwaEhS++sO04x5FmzpR+/3u3K/Igx7FpokOGSLt2SX/7m/Tww6wwBRk9ToCf6HGCF2Vn2+7P88/brtCbb0rx8W5X5UEbNlgn/qxZUs+e0ujRUmKi21V5Ej1OAIBCrVljg6ZHj5aee85WmghNQZaVZQ3f558vbdokffSRNH06oclFBCcAwDHyG8Bbt5Zyc20200MPSSX4GyN4HEeaMkVq0sT+MB5/XFq7Vrr6arcr8zz+NQCKkJycrKZNmyopKcntUoBit3Wr1KWLNX8PGiQtXSq1bOl2VR7z4492NMqNN0otWtjdck89xfk1IYIeJ8BP9DghkjmO9Pbb0sCBUlycDbW8/HK3q/KYrCzp2WelF16QatWyw/6uucbtqnAcVpwAwON277YTOu64w/qOV60iNAWV40hTp0pNm9qIgaFDbVuO0BSSGEcAAB42darUv7+dM/d//yfdfLPbFXnMxo3SfffZcKwrr5Rmz5bq13e7KpwCK04A4EF79tj5cjfcILVvb200hKYgOnDA5jA1b24TwKdPt0GWhKaQx4oTAHjMtGm2ypSTI02aZAGKI1OCxHEsJA0eLKWkSI8+ag8av8MGK04A4BF790q9e0vXX2/zmdaulXr1IjQFzaZNdjTK9ddbP9OaNdLw4YSmMENwAgAPmDpVatZM+vRT6b//tVWnGjXcrsojDhyQnnzS/gB++MFe/E8+kRo0cLsynAGCEwBEsF9/tQWOG26Q2ra1RY7evVllCor8bbmmTaVRo+xcuR9+kK69lj+AMEZwAorAAEyEo7w86dVX7e/sxYulDz+0hY6aNd2uzCM2b7Yp39ddZ9O/16yR/v53qWxZtyvDWWIAJuAnBmAiXPzwg3TPPdI339jzqFFShQpuV+URBw9KI0fai169ug2xZIUporDiBAARIjtbGjbMjkjZtUuaO1f6z38ITUHz0Ue2xDdypPTgg9K6dbbiRGiKKIwjAIAIsGCBrS5t3Gh3tz/+uBQb63ZVHrF5s3T//dbw3b279MUXUsOGbleFYsKKEwCEsbQ06S9/kTp1sjPmVqywVhpCUxAcPGhLfM2aSatXS1OmSDNnEpoiHCtOABCGHMeOSHngASkjQxo71gJUdLTblXnExx/bUSnbt0sPPSQ99ph0zjluV4UgYMUJAMLM+vVS167SrbdK7dpZM/igQYSmoPjpJzsJ+ZprbGVpzRrpH/8gNHkIK04AECaysqRnnpFefFGqU8eGWfbo4XZVHjBjhjRrlrR7t00SrVZNmjzZBmTR+O05BCcACHGOYzOYBg+WUlPtbNiHH6aPKSj++1+pT5+Ct2+6SZo4kRUmDyM4AUVITk5WcnKyfD6f26XAgzZvlu6913qOr7pK+vpr6Xe/c7uqCHf4sC3nTZxok7/zRUdLCQmEJo9jACbgJwZgIpgyMqQRI6SXXpLi422OYs+e7AwVq1WrLCxNmmSDsC68UGrdWho/3kKTz2dBqmdPtyuFi1hxAoAQkpcnvfWWNHSotH+/9Mgjti3HIkcx2b1beucdC0wrVlj/0h13SH37ShdcYJ9z9dXSnDnSZZcRmsCKE+AvVpxQ3L75xvqYli61O+ZGjbImcATY4cPSZ59JEybYWPCxO0kAAB77SURBVAHHsbvk+vWzbvtSpdyuECGMFScAcNnWrbay9N57tjM0f77UsaPbVUWg1asLtuJ++01q1Up64QXpttukqlXdrg5hguAEAC45cEB67jl7xMVJb7xhO0QlmLAXOLt3S+++a4Fp+XILSL172wvdooXb1SEMEZwAIMgcx/4uf+QRW/gYMsQGT5cr53ZlEeLwYenzzy0szZhhL/jVV0tPPmlbcaVLu10hwhjBCQCCaP58a/ZevNjmJ77wAuMFAmbNmoKtuNRUqWVL6fnnpdtvZysOAUNwAoAgWLfO7pSbPt3ucp89W+rc2e2qIsCePQVbccuWSVWqSL16WaN3y5ZuV4cIRHACgGK0c6c0bJj02mt2h9w770h/+AN9TGclN/fYrbi8PJsO+re/SVdeyVYcihXBCSgCk8NxJjIybBvuhRfsaJQXXpAGDJBiYtyuLIytXWth6b//ta24Cy6wmQ23327zl4AgYI4T4CfmOMEfOTm2ujR8uJSWZnOZHn1UqlDB7crC1N69Nqdh4kRpyRKpcmXbirvzTrbi4ApWnAAgAHJzpbfftm25X36xc2GffpoBlmckN1f64ouCs+J8PtuKmzLFntmKg4sITgBwFvLypA8/lJ56Slq/XrrhBhtG3ayZ25WFoR9+KNiKS0mRzj9fGjnStuKqV3e7OkCSRHsiPGHEiBFKSkpSuXLlVK1aNV133XXasGGD22UhjDmOBaQLL7Rm78REOypl8mRC02nZt0969VXpoovshXv9denmm+0Oue+/lx54gNCEkEJwgifMnTtXAwcO1OLFizVr1izl5uaqW7duysrKcrs0hKGvvpIuvtiON6tQwWYzffqpHZcCP+TmSjNnWuKMj5fuvdfC0eTJ0o4d0pgxlkijotyuFDgBzeHwpF27dqlatWqaO3euLrnkEr++huZwLFokPf649PXXUtu20jPPSFdcwd/vflu3rmArbudOW2G68047AoVVJYQJepzgSWlpaZKkSpUquVwJwsHChdbo/fnndgf89Om22kRg8sO+fdL771tg+vZbqVIl61nq149VJYQlghM8x3EcDRkyRB07dlTz5s1P+nnZ2dnKzs4+8nZ6enowykMImTvXAtPs2bY48t571n7D8Moi+HzSrFkWlqZNs625Hj2si/7qqxlmhbBGcILnDBo0SKtWrdKCBQtO+XkjRozQ8OHDg1QVQoXjWA/T3/8uzZsntWhhrTfXXUdgKtL69QVbcTt2SE2b2n5mr15SjRpuVwcEBD1O8JR7771X06ZN07x585SYmHjKzy1sxSkhIYEepwjlONJnn9kK0+LFUps20pNP2gIJu0mnsH9/wVbc4sVSxYoFW3GtW/PiIeKw4gRPcBxH9957r6ZOnao5c+YUGZokKSYmRjFsKUS8vDwbK/D3v9s4gXbt7Iav7t35O/+kfD7pyy8tLE2dKh0+bFtxH3xgzV/8e4MIRnCCJwwcOFDvvPOOpk+frnLlyiklJUWSFBcXpzJlyrhcHdxw+LD1LI0aZUegdepkbTlduhCYTmrDBunNN6W33pJ+/VVq0sQSZ+/ebMXBM9iqgydEneRvwgkTJqhfv35+XYNxBJEhK8vOknvpJWnrVunKK6VHHpH8nErhPWlpBVtxixbZ4Kr8rbg2bUiZ8BxWnOAJ/P8Bdu+WXn5ZGjvWssBtt0kPP2yneuA4Pp91yOdvxeXk2N7l++9LPXtKsbFuVwi4huAEIKL98ov04ot2kofjSH/6kzRkiFSvntuVhaAffyzYitu+XWrcWBo+3LbiatZ0uzogJBCcAESk776TRo+W/u//pLg46aGHpEGDpCpV3K4sxKSl2Ys0caJN+qxQwZbj+vWTkpLYigOOQ3ACEDFyc21nafRoywC/+531Mv3xj9I557hdXQjx+ezcmAkTpClTbCuuWzfrlr/2WrbigFMgOAFFSE5OVnJysnw+n9ul4CT277eG77FjreH70kttYPXVV0vR0W5XF0I2bizYitu2TWrUSBo2zLbiatVyuzogLHBXHeAn7qoLPZs2SWPGSG+8YYsmt90mDR4stWrldmUhJD29YCvum29s3/LWW+1w3bZt2YoDThMrTgDCiuNIc+ZI//ynDa6sXNmavQcMkOLj3a4uROTl2VbcxIl2XsyhQ7YV9+67thXH7DLgjBGcAISFAwfs7/2xY6Xvv7dDd8ePt5FC5ID/2bTJtuLefNO24s47z86N6d1bql3b7eqAiEBwAhDSNm2SXn3V+pj377eBlS+8wITvIzIyCrbiFiyQype3rbh+/ez8GF4kIKAITgBCjs8nffqp9MordvBupUo2f6l/f7tTzvPy8my/Mn8r7uBBqWtX6Z13pOuuYwkOKEYEJwAhY+dOW1kaP176+Wc70WPCBOkPfyALSJI2by7Yitu6VWrYUHr8cemOO6SEBLerAzyB4ATAVT6f9PnnFpY++kgqXVq65RY73aNtW7erCwEZGdKHH9rq0rx5thX3hz/YVlz79mzFAUFGcALgim3bbIzA66/bP7doYaMFbr/dhld7Wl6eNHeuhaUPP7StuCuukN5+27biypZ1u0LAswhOQBEYgBk4ubnSJ5/Y6tLMmbb9dvvt0t1327ac5xdPfvqpYCvul1+kBg2kxx6zrbg6ddyuDoAYgAn4jQGYZ27LFltZeuMN62NKSrKwdOutUrlyblfnsszMgq24uXPtBcnfirv4YtIkEGJYcQJQLA4csHPjJk6UvvrK8kDv3haYWrZ0uzoXzZghzZ4tVa1qR6B8+KG9WJdfLk2aJF1/PVtxQAgjOAEIGMexUz0mTrTRQhkZ0iWX2ErTzTd7/KDd9HQ7cXj48IL3xcdLjz5qW3F167pXGwC/EZwAnLWtW+3c2DfftIGV9erZMSh9+nh47tKvv9pAyvzHqlXW9J2vRAnbq/zb39yrEcBpIzgBOCNpadKUKba79PXXtrt0003W+H3JJZYLPCMvT1q/3gLS/Pn2/PPP9rEGDaSOHaVBg+wk4gEDpOhom8PQubOrZQM4fQQnAH7Lzra74d5+22Yu5eRIl11mW3E33SSde67bFQZJdra0bFnBatI330h791ogatXKRgZ07Ch16HDiycO1atnU78suk3r2dKN6AGeBu+oAP3n1rrq8PJu7+Pbb1se8f781d/fqZTtNnjg7dv9+aeHCgqD03XcWns45x4ZQduxoj4su8lB6BLyJFScAJ3Aca8l5+23p3Xel7dutb2nAAAtMTZu6XWEx27bt2P6k1avtRale3QLSyJH23KKFVKqU29UCCCKCE1AELw3AXLtW+uADuyNu3TqpcmU7/qRXrwgeKZSXZz/40UFp61b7WKNGFpAeeMCe69eP0BcBgL/YqgP8FKlbdevWWVD6v/+TfvjBjkK79loLTN262dlxEeXQIWnp0oJG7oULbSuuZEnpwgsLtt06dJCqVXO7WgAhhhUnwIM2bCgIS2vW2HDKa6+1Hahu3aSYGLcrDKC9e4/tT1qyxLrazz3XltGGDLGg1LatxwdNAfAHK06An8J5xclxbGVpyhTbilu1ynJDz562stS9uxQb63aVAeA4dsbb0dtua9fax2rUsIDUqZM9n3++rTIBwGkgOAF+CrfglJdniytTp9rjxx9tQeXosFSmjNtVniWfz5bMjg5K27fbx5o0Kdh269hRSkykPwnAWeN/t4AIcviwnRM7dao0bZq0Y4dUpYqFpRdflK64IsxXlg4etFEA+SFp4UI7yqRUKalNG+m22ywkXXyx/eAAEGAEJyDMHTggff65haWPP5b27ZPq1LGz4a6/3nqcw3ZHavduC0f507iXLbN0WL68haOHH7aglJTEwbgAgiJc/3MKeNq2bdInn9jjq69sIaZZM2ngQAtLrVqF4a6U40hbthy77bZunX2sVi3rTerd24JS8+Y2pRsAgoweJ8BPbvY4+Xy2Q/XxxxaWvv/eckPHjtJVV9kJHw0bBrWks5eba13qRwelnTvtY82aFfQmdepkS2hhlwQBRCJWnIAiuDUAMy3NtuA++UT69FPbtapcWerRQxo61Jq7K1QIaklnJyvr2P6kRYukjAwbFJWUJPXpU9CfVKmS29UCQKFYcQL8VNwrTo5jd77lryrNn2+LMhdcYKtKV19tR6GFzQ7Vb7/Z4bf5QWn5cvuBKlSwxqv8FaU2bcK8Yx2Al7DiBLgoLU2aPdtWlr74wlp8YmOlyy+XxoyxwFSnjttV+sFxpM2bC6ZxL1hgKVCyH6BjR6lfP3tu1kwqUcLVcgHgTBGcgCDy+ey0jy++sLC0eLG9r2FDC0ndu1toCvkbxHJzpZUrj+1PSk21PqTzz7e5B8OG2cpSWCQ/APAPwQkoZtu2FQSlL7+0cQHly0tdukjJyXbESWKi21UWITNT+vbbgtWkxYutZykmxo4quesua+Ju3z7MGq8A4PTQ4wTPmDdvnp5//nktW7ZMO3fu1NSpU3Xdddf5/fX+9jjt3y/Nm2dbcF98YXfUlyhh/c/du1tQuuiiEJ+tlJJybH/SihW2NFax4rHTuFu3jrCD7QDg1EL5P91AQGVlZalFixa68847deONNwbwupYtZs+2x/LldtxJnTpS167S8OG2uhSyN4rld6Ufve22aZN9LDHRAtLdd9tz48b0JwHwNIITPKNHjx7q0aPHWV8nO9uONckPSt9+a8Os4+OtP6l/f6lz5xA+Gu3wYVtBOjoo7dplxbZoIf3+9xaSOnSQatd2u1oACCkEJ+AksrOzlZ2drZwcaeXKaM2adUiSlJBg4alSJQtI//ynBabGjUM0KKWnW09SfkhavNhGjcfGSu3aSX/+swWldu2kuDi3qwWAkEZwAo6TmWmzGZ95ZrHmzcuTdJGkspL2SZKefFK68kqbrxSSu1Y7dlh/Un4j9/ff295h5coWkJ5+2hq5W7Wy4ZMAAL8RnOB5u3YVjB+aP7+gD7pKlUt09dW5at/ep/btM5WYaNtvgwbZXXGumzHD9gobNrQAlL+i9NNP9vH69S0oDRhgz40aheiSGACED+6qg6c4jvTLL1Ji4h3q2vXv2ratntavt4/VrWsLMfmP47fe3DyrTjk5Nh3zxx/tMWuWzTfIFxVlK0idOhX0J9WoEdwaAcADCE6IaNafZG09ixbZgsz27faxOnXSdOWVcUeCUkLCqa9V7MHJ55O2bpU2brRwdPTzli223SbZdMyyZaU9eywJligh/eUv0ssvB74mAMAx2KpDRNmxwwJSflBatkw6dMhGDbVo4VOXLrvVqlWWBg9uo8GDn1Dnzp1VqVIlJSQEabq140g7dxYejjZtsqQnSaVK2VbbeedJ111nzw0b2nPNmtJHH0nXXmsH1/l8NhwKAFDsWHFC2MrOttWkRYsKwtLWrfaxOnVsiHW7dvbcsqW0aNEcde7c+YTr9O3bVxMnTizy+53WitOePSeGo/x/zsqyzylRwvYHjw5F+c916hQ9IXPGDGnOHOmyy6SePYusHwBw9ghOCAu5udIPP0hLlthZb0uWSKtW2UiimBipTZtjg1LNmoGv4YTglJl5YijKf967t+ALa9YsPBz97ndM3QaAMMNWHUJOXp5lj/yAtGSJ3el28KD1QDdtakHpzjvtGJOWLYvxrvpDh6TNm62gVavsfT16WM/Rzp0Fn1e5soWhRo2ka64pCEgNGkjnnltMxQEAgo3gBFc5jm2v5QekpUutLyktzT5ev76FoxtvtLB04YXFkENyc6Wffy687+iXX6xIqeAb16hhZ6jkh6OGDUP4PBUAQCCxVYeg8fmkDRts9WjlSntesaJgV6t2bQtHSUn2aN06gHkkL0/69dfC+45++snCk2RbZw0aHLOlNmXNGo3+5BOlSvpx40Z3xhEAAEICwQnF4uBBafXqY0PSqlX2fkmqV8+22Fq1skdSkp31dlYcx6ZZFtZ3tGlTwTePjrZJluedd2LvUULCSceBuzrHCQAQEghOOCv5d9evXm2P77+3kLR+va0wRUfbIMn8gNSqlQWmihXP4pumpZ24pZb/nL/HJ1kIKiwcJSba7f6nieAEACA4wW/p6dKaNfbID0qrVxdstZUtK51//rEhqXlzqUyZM/hmBw7YKlFh4ei33wo+r3r1E+9Wa9jQmqPKlg3Iz52P4AQAIDjhBIcPWy/S0eFo9Wrrk5ZsFalhQwtJRz8SE0/z0Nv8Y0QK6zvKH+8tSXFxBStHxwekIAYYghMAgODkYYcPW1/0unX2yF9JWr/ePiZJtWqdGJAaN5ZiY/38Jj6ftG1b4eHo55/t45ItSxU266hhQ6lKlZA4nJbgBAAgOHlAVpatIOUHpPzHpk0FASkuzrbVzj//2Ge/7mpzHCklpfCm7M2bbcS3VHCMSGHhqGbN01yuCj6CEwCA4BRBdu2y1aLjA1L+MSSSrSA1aWKrRk2aFDyqV/djUWfv3pNPys7MtM+JirJb5o4PR/4eIxLCCE4AAIJTmDl4sKBn+ugMs26dHY8mWQ9S/fonBqTGjU/REjRjhvT113ZmyXnnFd6Unf8NJE8eI0JwAgCE7//+R7CT9Uxv3GjtQvnKly/ILN26FQSkBg0KyS45OXY32sZUKTXVttbyn5cvlxYsOLGQ/GNEzjtPuvrqgnDksWNEkpOTlZycLF9+PxYAwLNYcXLJwYMWjjZvtgbt/OPQNm48sWc6/1SP4xd4qlbMVdTuXceGoJM9H33obL7KlW3qZHq63cXmONZndMstUnIyx4gchxUnAAArTsXEcWyBJz8Y5Yej/Oejz4eNibFb+Rs2lK7v6dP5NXarccUUJZZJVeXDKSqx63+rRFtTpCVHBaLduwvOUctXsaI1LMXH2/MFFxz7dv5ztWoFQyBnzJCuvdb2+Hw+6bbbCE0AABSCFaezkJ5us41++cVWiY4PSAcO2OdFKU+Nq+xRq5qpalYlReeVT1XdmBTFR1kwKpOeqqjU/4WhXbvsXLWjxcUVHn6Of1+1amfeXzRjhjRnjnTZZVLPnmfxqkQuVpwAAASnk3AcW9DJD0YnPH52FLV/r6orVfFKUa3oVDWplKIG56aqTukUVVeqKuak6JzMVJXc+5uiju+POffcE4NQYc/Vqp3h6G0EGsEJAODZrTqfz7bLjg9Duzen6cBPKfLtSFWF7JRjglHr2BTViE5VVV+Kyh/6TdE6fNQFJWWVlcrFS5Xyw0+7wleJqleXzjnHtZ8dAACcmYgMTrm51ga0fbv063ZHv23OUPqPKTr4c6p8O1JU4rdUxexLUVXHQlEDpapTiVRVd1JU2sk55lp5MbFyqlVXiZrxiqpeXYpvffLVIQ/daQYAgBeFXXDKzpZ27JB2/JipvT+kKGNTqg79kirfrymK2pWq2H0pKncgVdVlq0UtlKIyOnTMNXKjS+tgheryVYlXiZrVFVuvhUrXLryHqET58iFx3AcAAHBfyASn7Gwp5acD2vNDqtJ/TNGBLak6vM0apkvuTlFseqrKZ6WoUq6tEiXqwDFfnxtVUullqutQ+erK/V28omo0VUydy5VXv7qcxHhFxRcEopIVKqgcYQgAAJymwAen/AnUnTtLPXsqY9ch7VqTqv3rU5S5OVXZW1OVt8NusY/ZZ83TFbJTVDUvVXWVobpHXSpX0dpXqpoyylTXgbh4+RLPU1aNS7S9TnWd2yBeFRpV1zn1/xeGKlZUpRA/6wwAAIS3wN5V9795QI6kKEmZKqtzj1sZylOU9kRXU1pMdWWVq67sivFyqloPUcxRgahC43iVqFo55A9+ReQ7enL4jz/+yF11AOBhgQ1ODzygvNH/Ugk5ylOUdiR21Lbuf1LZ38Ur7rzqqtwsXufWq6KoktEB+5ZAsDCOAAAQ2K26zp1VYvRoKTpaJXw+1R79oGozTBEAAESIwAannj2l6dOZQA0AACISk8MBP7FVBwCg8xoAAMBPBCcAAAA/EZwAAAD8RI8T4CfHcZSRkaFy5copisnzAOBJBCcAAAA/sVUHAADgJ4ITAACAnwhOAAAAfiI4AQAA+IngBAAA4CeCEwAAgJ8ITgAAAH4iOAEAAPiJ4AQAAOAnghMAAICfSvrzSflndAEAAEQqf84i9Ss4ZWRkKC4uLiBFAQAAhKK0tDSVL1/+lJ/j1yG/p7PilJ6eroSEBG3btq3Ib16UpKQkLVmy5KyuEehrhdJ1AvlaB6qmSL5OKP5uh9p1AnWtUHytA3mtULoO/x0J7nVC8Xc71K4TqGudyWsdsBWnqKio0/4DLl++/Fn/UkRHRwfkX+RAXivUriMF5rWWQu9nC7Xr5Aul3+1Qu06grxVKr3UgrxVq15H470iwrpMvlH63Q+06gb5WoH6384V0c/jAgQND7lqhdp1ACrWfLdSuE0ih9rOF4r9rgRKKP1uoXSeQQu1nC7XrBFKo/Wyh+O9acfBrq+50pKenKy4uzq99QpwdXuvg4vUOHl7r4OG1Di5e7+Aprtc6etiwYcMCdrX8i0ZH67LLLlPJkn7tBOIs8FoHF6938PBaBw+vdXDxegdPcbzWAV9xAgAAiFQh3eMEAAAQSghOAAAAfiI4AQAA+IngBAAA4KeAB6cpU6aoe/fuqlKliqKiorRy5cpAfwtPcRxHw4YNU82aNVWmTBlddtllWrt27Sm/ZuLEiYqKijrhcejQoSBVHXleeeUVJSYmKjY2Vq1bt9b8+fPdLinsnc5rOmfOnEJ/p9evXx/EiiPLvHnzdM0116hmzZqKiorStGnT3C4p7J3ua8rvdfEYMWKEkpKSVK5cOVWrVk3XXXedNmzYELDrBzw4ZWVlqUOHDho5cmSgL+1Jzz33nF566SW9/PLLWrJkieLj49W1a9cij8ApX768du7cecwjNjY2SFVHlvfff1+DBw/W448/rhUrVqhTp07q0aOHtm7d6nZpYetMX9MNGzYc8zvdsGHDIFUcebKystSiRQu9/PLLbpcSMc70NeX3OrDmzp2rgQMHavHixZo1a5Zyc3PVrVs3ZWVlBeYbOMVky5YtjiRnxYoVxfUtIl5eXp4THx/vjBw58sj7Dh065MTFxTn//ve/T/p1EyZMcOLi4oJRoie0bdvW6d+//zHva9y4sfPoo4+6VFH4O93X9Ouvv3YkOfv27QtGeZ4jyZk6darbZUQUf15Tfq+D47fffnMkOXPnzg3I9ehxCmFbtmxRSkqKunXrduR9MTExuvTSS7Vw4cJTfm1mZqbq1q2r2rVr6+qrr9aKFSuKu9yIlJOTo2XLlh3zZyBJ3bp1K/LPAIU7m9e0VatWqlGjhrp06aKvv/66OMsEgobf6+KVlpYmSapUqVJArkdwCmEpKSmSpOrVqx/z/urVqx/5WGEaN26siRMnasaMGXr33XcVGxurDh06aOPGjcVabyTavXu3fD7faf8Z4OTO5DWtUaOGxo0bp8mTJ2vKlClq1KiRunTponnz5gWjZKBY8Htd/BzH0ZAhQ9SxY0c1b948INc8qxnkb7/9tv785z8feXvmzJnq1KnTWRflVce/np988okkKSoq6pjPcxznhPcdrV27dmrXrt2Rtzt06KALL7xQY8eO1ZgxYwJctTec7p8BinY6r2mjRo3UqFGjI2+3b99e27Zt0wsvvKBLLrmkWOsEigu/18Vv0KBBWrVqlRYsWBCwa55VcOrZs6cuuuiiI2/XqlXrrAvysuNfz+zsbEm28lSjRo0j7//tt99O+L/1UylRooSSkpJYcToDVapUUXR09AkrIaf7Z4ACgXpN27Vrp0mTJgW6PMBV/F4Hzr333qsZM2Zo3rx5ql27dsCue1ZbdeXKlVODBg2OPMqUKROoujzp+NezadOmio+P16xZs458Tk5OjubOnauLL77Y7+s6jqOVK1ceE77gn9KlS6t169bH/BlI0qxZs07rzwAFAvWarlixgt9pRBx+r8+e4zgaNGiQpkyZotmzZysxMTGg1w/40cx79+7V1q1btWPHDkk6MjshPj5e8fHxgf52ES0qKkqDBw/Ws88+q4YNG6phw4Z69tlnVbZsWd1+++1HPq9Pnz6qVauWRowYIUkaPny42rVrp4YNGyo9PV1jxozRypUrlZyc7NaPEtaGDBmiO+64Q23atFH79u01btw4bd26Vf3793e7tLBV1Gs6dOhQ/frrr3rrrbckSaNHj1a9evXUrFkz5eTkaNKkSZo8ebImT57s5o8R1jIzM7Vp06Yjb2/ZskUrV65UpUqVVKdOHRcrC19Fvab8XgfHwIED9c4772j69OkqV67ckdXtuLi4wCzwBOTevKNMmDDBkXTC46mnngr0t/KEvLw856mnnnLi4+OdmJgY55JLLnFWr159zOdceumlTt++fY+8PXjwYKdOnTpO6dKlnapVqzrdunVzFi5cGOTKI0tycrJTt25dp3Tp0s6FF14YsNtavexUr2nfvn2dSy+99Mjbo0aNcurXr+/ExsY6FStWdDp27Oh88sknLlQdOfJvhT/+cfR/S3B6inpN+b0OjsL+DCQ5EyZMCMj1o/73TQAAAFAExhEAAAD4ieAEAADgJ4ITAACAnwhOAAAAfiI4AQAA+IngBAAA4CeCEwAAgJ8ITgAAAH4iOAEAAPiJ4AQAAOAnghMAAICfCE4AAAB++n9rbrq56HnpOgAAAABJRU5ErkJggg==\n", | ||
| 148 | "text/plain": [ | ||
| 149 | "Graphics object consisting of 2 graphics primitives" | ||
| 150 | ] | ||
| 151 | }, | ||
| 152 | "execution_count": 4, | ||
| 153 | "metadata": {}, | ||
| 154 | "output_type": "execute_result" | ||
| 155 | } | ||
| 156 | ], | ||
| 157 | "source": [ | ||
| 158 | "var('y')\n", | ||
| 159 | "\n", | ||
| 160 | "def euler_desolve(f, x0, y0, x1):\n", | ||
| 161 | " n = 5\n", | ||
| 162 | " h = (x1-x0)/n\n", | ||
| 163 | " S = []\n", | ||
| 164 | " Y = [y0]\n", | ||
| 165 | " for i in range(n+1):\n", | ||
| 166 | " S.append(x0 + i*h)\n", | ||
| 167 | " Y.append(N( Y[i] + h*f(S[i], Y[i]) ))\n", | ||
| 168 | " return S, Y\n", | ||
| 169 | "\n", | ||
| 170 | "f(x,y) = y\n", | ||
| 171 | "x0 = -1\n", | ||
| 172 | "x1 = 2\n", | ||
| 173 | "y0 = e^(-1)\n", | ||
| 174 | "\n", | ||
| 175 | "S, Y = euler_desolve(f, x0, y0, x1)\n", | ||
| 176 | "plot(e^x, -1, 2) + line([(S[i], Y[i]) for i in range(len(S))], color='red', marker='o', markersize=2)" | ||
| 177 | ] | ||
| 178 | }, | ||
| 179 | { | ||
| 180 | "cell_type": "markdown", | ||
| 181 | "metadata": {}, | ||
| 182 | "source": [ | ||
| 183 | "Sage also has an `eulers_method()` function \"for pedagogical purposes only\":" | ||
| 184 | ] | ||
| 185 | }, | ||
| 186 | { | ||
| 187 | "cell_type": "code", | ||
| 188 | "execution_count": 5, | ||
| 189 | "metadata": {}, | ||
| 190 | "outputs": [ | ||
| 191 | { | ||
| 192 | "name": "stdout", | ||
| 193 | "output_type": "stream", | ||
| 194 | "text": [ | ||
| 195 | " x y h*f(x,y)\n", | ||
| 196 | " -1 0.367879441171442 0.0367879441171442\n", | ||
| 197 | "-0.900000000000000 0.404667385288587 0.0404667385288587\n", | ||
| 198 | "-0.800000000000000 0.445134123817445 0.0445134123817445\n", | ||
| 199 | "-0.700000000000000 0.489647536199190 0.0489647536199190\n", | ||
| 200 | "-0.600000000000000 0.538612289819109 0.0538612289819109\n", | ||
| 201 | "-0.500000000000000 0.592473518801020 0.0592473518801020\n", | ||
| 202 | "-0.400000000000000 0.651720870681122 0.0651720870681122\n", | ||
| 203 | "-0.300000000000000 0.716892957749234 0.0716892957749234\n", | ||
| 204 | "-0.200000000000000 0.788582253524157 0.0788582253524157\n", | ||
| 205 | "-0.100000000000000 0.867440478876573 0.0867440478876573\n", | ||
| 206 | "-1.38777878078145e-16 0.954184526764230 0.0954184526764230\n", | ||
| 207 | "0.0999999999999999 1.04960297944065 0.104960297944065\n", | ||
| 208 | "0.200000000000000 1.15456327738472 0.115456327738472\n", | ||
| 209 | "0.300000000000000 1.27001960512319 0.127001960512319\n", | ||
| 210 | "0.400000000000000 1.39702156563551 0.139702156563551\n", | ||
| 211 | "0.500000000000000 1.53672372219906 0.153672372219906\n", | ||
| 212 | "0.600000000000000 1.69039609441897 0.169039609441897\n", | ||
| 213 | "0.700000000000000 1.85943570386086 0.185943570386086\n", | ||
| 214 | "0.800000000000000 2.04537927424695 0.204537927424695\n", | ||
| 215 | "0.900000000000000 2.24991720167165 0.224991720167165\n", | ||
| 216 | "1.00000000000000 2.47490892183881 0.247490892183881\n", | ||
| 217 | "1.10000000000000 2.72239981402269 0.272239981402269\n", | ||
| 218 | "1.20000000000000 2.99463979542496 0.299463979542496\n", | ||
| 219 | "1.30000000000000 3.29410377496746 0.329410377496746\n", | ||
| 220 | "1.40000000000000 3.62351415246420 0.362351415246420\n", | ||
| 221 | "1.50000000000000 3.98586556771062 0.398586556771062\n", | ||
| 222 | "1.60000000000000 4.38445212448168 0.438445212448168\n", | ||
| 223 | "1.70000000000000 4.82289733692985 0.482289733692985\n", | ||
| 224 | "1.80000000000000 5.30518707062284 0.530518707062284\n", | ||
| 225 | "1.90000000000000 5.83570577768512 0.583570577768512\n", | ||
| 226 | "2.00000000000000 6.41927635545363 0.641927635545363\n" | ||
| 227 | ] | ||
| 228 | } | ||
| 229 | ], | ||
| 230 | "source": [ | ||
| 231 | "# Usage: eulers_method(f, x0, y0, h, x1)\n", | ||
| 232 | "eulers_method(f, -1, N(e^(-1)), 0.1, 2)" | ||
| 233 | ] | ||
| 234 | }, | ||
| 235 | { | ||
| 236 | "cell_type": "markdown", | ||
| 237 | "metadata": {}, | ||
| 238 | "source": [ | ||
| 239 | "## Solving the heat equation with a finite difference method" | ||
| 240 | ] | ||
| 241 | }, | ||
| 242 | { | ||
| 243 | "cell_type": "code", | ||
| 244 | "execution_count": null, | ||
| 245 | "metadata": {}, | ||
| 246 | "outputs": [], | ||
| 247 | "source": [ | ||
| 248 | "def heat_fdm(u0j, u1j, ui0):\n", | ||
| 249 | " m, n = len(u0j)-1, len(ui0)-1\n", | ||
| 250 | " k, h = 1/m, 1/n\n", | ||
| 251 | " \n", | ||
| 252 | " u = [[0] * (m+1) for i in range(n+1)]\n", | ||
| 253 | " for j in range(m+1):\n", | ||
| 254 | " u[0][j] = u0j[j]\n", | ||
| 255 | " for j in range(m+1):\n", | ||
| 256 | " u[n][j] = u1j[j]\n", | ||
| 257 | " for i in range(n+1):\n", | ||
| 258 | " u[i][0] = ui0[i]\n", | ||
| 259 | " \n", | ||
| 260 | " for j in range(0,m):\n", | ||
| 261 | " for i in range(1,n):\n", | ||
| 262 | " u[i][j+1] = (k/(h*h)) * (u[i+1][j] - 2*u[i][j] + u[i-1][j]) + u[i][j]\n", | ||
| 263 | " \n", | ||
| 264 | " return u\n", | ||
| 265 | "\n", | ||
| 266 | "n, m = 20, 20\n", | ||
| 267 | "u0j = [10 - (j/m)*10 for j in range(m+1)] # One extreme goes from hot to cold\n", | ||
| 268 | "u1j = [(j/m)*10 for j in range(m+1)] # The other does the opposite\n", | ||
| 269 | "ui0 = [10 - (i/m)*10 for i in range(0,n+1)]\n", | ||
| 270 | "\n", | ||
| 271 | "u = heat_fdm(u0j, u1j, ui0)\n", | ||
| 272 | "for t in range(m+1):\n", | ||
| 273 | " show(line([(i/n, u[i][t]) for i in range(n+1)], ymin=-1, ymax =12))" | ||
| 274 | ] | ||
| 275 | } | ||
| 276 | ], | ||
| 277 | "metadata": { | ||
| 278 | "kernelspec": { | ||
| 279 | "display_name": "SageMath 9.0", | ||
| 280 | "language": "sage", | ||
| 281 | "name": "sagemath" | ||
| 282 | }, | ||
| 283 | "language_info": { | ||
| 284 | "codemirror_mode": { | ||
| 285 | "name": "ipython", | ||
| 286 | "version": 3 | ||
| 287 | }, | ||
| 288 | "file_extension": ".py", | ||
| 289 | "mimetype": "text/x-python", | ||
| 290 | "name": "python", | ||
| 291 | "nbconvert_exporter": "python", | ||
| 292 | "pygments_lexer": "ipython3", | ||
| 293 | "version": "3.8.5" | ||
| 294 | } | ||
| 295 | }, | ||
| 296 | "nbformat": 4, | ||
| 297 | "nbformat_minor": 4 | ||
| 298 | } | ||
