From d6c61d988bfa4255baf9cdae42db59ebee38363f Mon Sep 17 00:00:00 2001 From: Sebastiano Tronto Date: Tue, 25 May 2021 17:10:49 +0200 Subject: Added files --- .../X1-ComputationalComplexity-checkpoint.ipynb | 394 ++++++++++++++++++++ ...mputationalComplexity-notebook-checkpoint.ipynb | 412 +++++++++++++++++++++ .../X2-StudentsRequests-checkpoint.ipynb | 165 +++++++++ .../X2-StudentsRequests-notebook-checkpoint.ipynb | 298 +++++++++++++++ 4 files changed, 1269 insertions(+) create mode 100644 src/Lecture7/slides/.ipynb_checkpoints/X1-ComputationalComplexity-checkpoint.ipynb create mode 100644 src/Lecture7/slides/.ipynb_checkpoints/X1-ComputationalComplexity-notebook-checkpoint.ipynb create mode 100644 src/Lecture7/slides/.ipynb_checkpoints/X2-StudentsRequests-checkpoint.ipynb create mode 100644 src/Lecture7/slides/.ipynb_checkpoints/X2-StudentsRequests-notebook-checkpoint.ipynb (limited to 'src/Lecture7/slides/.ipynb_checkpoints') 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 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Nested loops\n", + "\n", + "The following two functions compute sum and product of matrices, respectively.\n", + "\n", + "By counting the nested loops it is easy to see that `add()` is $O(n^2)$ while `prod()` is $O(n^3)$." + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time for add: 0.00012074300000008975\n", + "Time for prod: 0.00036587199999971176\n" + ] + } + ], + "source": [ + "from random import randint\n", + "import time\n", + "\n", + "def add(A, B):\n", + " S = [[0] * len(A) for i in range(len(A))]\n", + " for i in range(len(A)):\n", + " for j in range(len(A)):\n", + " S[i][j] = A[i][j] + B[i][j]\n", + " return S\n", + "\n", + "def prod(A, B):\n", + " S = [[0] * len(A) for i in range(len(A))]\n", + " for i in range(len(A)):\n", + " for j in range(len(A)):\n", + " for k in range(len(A)):\n", + " S[i][j] = S[i][j] + A[i][k] * B[k][j]\n", + " return S\n", + "\n", + "N = 10\n", + "A = [ [randint(0,100) for i in range(N)] for j in range(N) ]\n", + "B = [ [randint(0,100) for i in range(N)] for j in range(N) ]\n", + "\n", + "t0 = time.process_time()\n", + "add(A,B)\n", + "t1 = time.process_time()\n", + "prod(A,B)\n", + "t2 = time.process_time()\n", + "\n", + "print(\"Time for add: \", t1-t0)\n", + "print(\"Time for prod:\", t2-t1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sorting a list, slow version\n", + "\n", + "The following code implements a slow version of the so-called *insertion sort* alogithm\n", + "\n", + "Complexity: $O(n^2)$." + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running time: 1.1012288430000012\n" + ] + } + ], + "source": [ + "from random import randint\n", + "import time\n", + "\n", + "def correct_position(e, S):\n", + " for i in range(len(S)):\n", + " if S[i] > e:\n", + " return i\n", + " return len(S)\n", + "\n", + "def sort_list(L):\n", + " S = []\n", + " for e in L:\n", + " cp = correct_position(e, S)\n", + " S.insert(cp, e)\n", + " return S\n", + "\n", + "N = 10000\n", + "L = [randint(0,10**9) for i in range(N)]\n", + "\n", + "t0 = time.process_time()\n", + "sort_list(L)\n", + "t1 = time.process_time()\n", + "\n", + "print(\"Running time:\", t1-t0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Binary search\n", + "\n", + "The following code implements a binary search.\n", + "\n", + "Complexity: $O(\\log_2(n))$" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The correct position of e = 658230309 in L is:\n", + "... 658211821 658224379 e 658234625 658246765 ...\n", + "\n", + "Time for sorting: 0.021211020999999164\n", + "Time for searching: 7.820199999741817e-05\n" + ] + } + ], + "source": [ + "from random import randint\n", + "import time\n", + "\n", + "def binary_search(e, S, start, end):\n", + " if start == end:\n", + " return start\n", + " midpoint = (start+end) // 2\n", + " if e < S[midpoint]:\n", + " return binary_search(e, S, start, midpoint)\n", + " else:\n", + " return binary_search(e, S, midpoint+1, end)\n", + " \n", + "N = 100000\n", + "L = [randint(0,10**9) for i in range(N)]\n", + "e = randint(0,10**9)\n", + "\n", + "t0 = time.process_time()\n", + "L.sort() # Using Python's sort()\n", + "t1 = time.process_time()\n", + "i = binary_search(e, L, 0, len(L))\n", + "t2 = time.process_time()\n", + "print(\"The correct position of e =\", e, \"in L is:\")\n", + "print(\"...\", L[i-2], L[i-1], \"e\", L[i], L[i+1], \"...\")\n", + "print(\"\")\n", + "print(\"Time for sorting: \", t1-t0)\n", + "print(\"Time for searching:\", t2-t1)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sorting a list, fast version (with binary_search)\n", + "\n", + "The following code uses the function `binary_search()` above instead of `correct_position()` in our insertion sort algorithm.\n", + "\n", + "Complexity: $O(n\\log_2(n))$" + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running time: 0.03710268399998995\n" + ] + } + ], + "source": [ + "from random import randint\n", + "import time\n", + "\n", + "def binary_search(e, S, start, end):\n", + " if start == end:\n", + " return start\n", + " midpoint = (start+end) // 2\n", + " if e < S[midpoint]:\n", + " return binary_search(e, S, start, midpoint)\n", + " else:\n", + " return binary_search(e, S, midpoint+1, end)\n", + " \n", + "def sort_list(L):\n", + " S = []\n", + " for e in L:\n", + " cp = binary_search(e, S, 0, len(S)) # Changed here\n", + " S.insert(cp, e)\n", + " return S\n", + " \n", + "N = 10000\n", + "L = [randint(0,10**9) for i in range(N)]\n", + "\n", + "t0 = time.process_time()\n", + "sort_list(L)\n", + "t1 = time.process_time()\n", + "\n", + "print(\"Running time:\", t1-t0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fast exponentiation\n", + "\n", + "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", + "\n", + "Complexity: $O(n)$ for the slow algorithm, $O(\\log_2(n))$ for the other two." + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.71828179834636\n", + "2.7182817863957984\n", + "2.7182817983473577\n", + "Time for slow_power(): 3.234879998000004\n", + "Time for fast_power(): 9.059099999575437e-05\n", + "Time for Python's **: 0.00010159500000384014\n" + ] + } + ], + "source": [ + "import time\n", + "\n", + "def slow_power(a, n):\n", + " r = 1\n", + " for i in range(n):\n", + " r = r * a\n", + " return r\n", + "\n", + "def fast_power(a, n):\n", + " if n == 0:\n", + " return 1\n", + " if n%2 == 0:\n", + " return fast_power(a*a, n//2)\n", + " else:\n", + " return a * fast_power(a, n-1)\n", + "\n", + "a = 1.00000001\n", + "n = 100000000\n", + "\n", + "t0 = time.process_time()\n", + "print(slow_power(a, n))\n", + "t1 = time.process_time()\n", + "print(fast_power(a, n))\n", + "t2 = time.process_time()\n", + "print(a**n)\n", + "t3 = time.process_time()\n", + "\n", + "print(\"Time for slow_power():\", t1-t0)\n", + "print(\"Time for fast_power():\", t2-t1)\n", + "print(\"Time for Python's **: \", t3-t2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fast gcd\n", + "\n", + "Complexity: $O(\\log_2(n))$" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "126\n", + "Running time: 0.00017707599999994272\n" + ] + } + ], + "source": [ + "import time\n", + "\n", + "def gcd(a, b):\n", + " if b == 0:\n", + " return a\n", + " else:\n", + " return gcd(b, a%b)\n", + "\n", + "t0 = time.process_time()\n", + "print(gcd(155275387236018, 572335397352432))\n", + "t1 = time.process_time()\n", + "\n", + "print(\"Running time:\", t1-t0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fibonacci numbers\n", + "\n", + "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", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "F_memorized = [-1] * (10**6)\n", + "\n", + "def F_slow(n):\n", + " if n <= 1:\n", + " return n\n", + " else:\n", + " return F_slow(n-1) + F_slow(n-2)\n", + " \n", + "def F_fast(n):\n", + " if F_memorized[n] == -1:\n", + " if n <= 1:\n", + " F_memorized[n] = n\n", + " else:\n", + " F_memorized[n] = F_fast(n-1) + F_fast(n-2)\n", + " \n", + " return F_memorized[n]\n", + "\n", + "n = 40\n", + "\n", + "t0 = time.process_time()\n", + "print(F_slow(n))\n", + "t1 = time.process_time()\n", + "print(F_fast(n))\n", + "t2 = time.process_time()\n", + "\n", + "print(\"Time for F_slow:\", t1-t0)\n", + "print(\"Time for F_fast:\", t2-t1)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} 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 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Nested loops\n", + "\n", + "The following two functions compute sum and product of matrices, respectively.\n", + "\n", + "By counting the nested loops it is easy to see that `add()` is $O(n^2)$ while `prod()` is $O(n^3)$." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time for add: 0.005766554000000035\n", + "Time for prod: 1.3871021639999999\n" + ] + } + ], + "source": [ + "from random import randint\n", + "import time\n", + "\n", + "def add(A, B):\n", + " S = [[0] * len(A) for i in range(len(A))]\n", + " for i in range(len(A)):\n", + " for j in range(len(A)):\n", + " S[i][j] = A[i][j] + B[i][j]\n", + " return S\n", + "\n", + "def prod(A, B):\n", + " S = [[0] * len(A) for i in range(len(A))]\n", + " for i in range(len(A)):\n", + " for j in range(len(A)):\n", + " for k in range(len(A)):\n", + " S[i][j] = S[i][j] + A[i][k] * B[k][j]\n", + " return S\n", + "\n", + "N = 200\n", + "A = [ [randint(0,100) for i in range(N)] for j in range(N) ]\n", + "B = [ [randint(0,100) for i in range(N)] for j in range(N) ]\n", + "\n", + "t0 = time.process_time()\n", + "add(A,B)\n", + "t1 = time.process_time()\n", + "prod(A,B)\n", + "t2 = time.process_time()\n", + "\n", + "print(\"Time for add: \", t1-t0)\n", + "print(\"Time for prod:\", t2-t1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sorting a list, slow version\n", + "\n", + "The following code implements a slow version of the so-called *insertion sort* alogithm\n", + "\n", + "Complexity: $O(n^2)$." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running time: 1.1191449070000001\n" + ] + } + ], + "source": [ + "from random import randint\n", + "import time\n", + "\n", + "def correct_position(e, S):\n", + " for i in range(len(S)):\n", + " if S[i] > e:\n", + " return i\n", + " return len(S)\n", + "\n", + "def sort_list(L):\n", + " S = []\n", + " for e in L:\n", + " cp = correct_position(e, S)\n", + " S.insert(cp, e)\n", + " return S\n", + "\n", + "N = 10000\n", + "L = [randint(0,10**9) for i in range(N)]\n", + "\n", + "t0 = time.process_time()\n", + "sort_list(L)\n", + "t1 = time.process_time()\n", + "\n", + "print(\"Running time:\", t1-t0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Binary search\n", + "\n", + "The following code implements a binary search.\n", + "\n", + "Complexity: $O(\\log_2(n))$" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The correct position of e = 216197744 in L is:\n", + "... 216196218 216197540 e 216198673 216198962 ...\n", + "\n", + "Time for sorting: 0.26413054400000036\n", + "Time for searching: 9.616099999965044e-05\n" + ] + } + ], + "source": [ + "from random import randint\n", + "import time\n", + "\n", + "def binary_search(e, S, start, end):\n", + " if start == end:\n", + " return start\n", + " midpoint = (start+end) // 2\n", + " if e < S[midpoint]:\n", + " return binary_search(e, S, start, midpoint)\n", + " else:\n", + " return binary_search(e, S, midpoint+1, end)\n", + " \n", + "N = 1000000\n", + "L = [randint(0,10**9) for i in range(N)]\n", + "e = randint(0,10**9)\n", + "\n", + "t0 = time.process_time()\n", + "L.sort() # Using Python's sort()\n", + "t1 = time.process_time()\n", + "i = binary_search(e, L, 0, len(L))\n", + "t2 = time.process_time()\n", + "print(\"The correct position of e =\", e, \"in L is:\")\n", + "print(\"...\", L[i-2], L[i-1], \"e\", L[i], L[i+1], \"...\")\n", + "print(\"\")\n", + "print(\"Time for sorting: \", t1-t0)\n", + "print(\"Time for searching:\", t2-t1)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sorting a list, fast version (with binary_search)\n", + "\n", + "The following code uses the function `binary_search()` above instead of `correct_position()` in our insertion sort algorithm.\n", + "\n", + "Complexity: $O(n\\log_2(n))$" + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running time: 0.03710268399998995\n" + ] + } + ], + "source": [ + "from random import randint\n", + "import time\n", + "\n", + "def binary_search(e, S, start, end):\n", + " if start == end:\n", + " return start\n", + " midpoint = (start+end) // 2\n", + " if e < S[midpoint]:\n", + " return binary_search(e, S, start, midpoint)\n", + " else:\n", + " return binary_search(e, S, midpoint+1, end)\n", + " \n", + "def sort_list(L):\n", + " S = []\n", + " for e in L:\n", + " cp = binary_search(e, S, 0, len(S)) # Changed here\n", + " S.insert(cp, e)\n", + " return S\n", + " \n", + "N = 10000\n", + "L = [randint(0,10**9) for i in range(N)]\n", + "\n", + "t0 = time.process_time()\n", + "sort_list(L)\n", + "t1 = time.process_time()\n", + "\n", + "print(\"Running time:\", t1-t0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fast exponentiation\n", + "\n", + "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", + "\n", + "Complexity: $O(n)$ for the slow algorithm, $O(\\log_2(n))$ for the other two." + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.71828179834636\n", + "2.7182817863957984\n", + "2.7182817983473577\n", + "Time for slow_power(): 3.234879998000004\n", + "Time for fast_power(): 9.059099999575437e-05\n", + "Time for Python's **: 0.00010159500000384014\n" + ] + } + ], + "source": [ + "import time\n", + "\n", + "def slow_power(a, n):\n", + " r = 1\n", + " for i in range(n):\n", + " r = r * a\n", + " return r\n", + "\n", + "def fast_power(a, n):\n", + " if n == 0:\n", + " return 1\n", + " if n%2 == 0:\n", + " return fast_power(a*a, n//2)\n", + " else:\n", + " return a * fast_power(a, n-1)\n", + "\n", + "a = 1.00000001\n", + "n = 100000000\n", + "\n", + "t0 = time.process_time()\n", + "print(slow_power(a, n))\n", + "t1 = time.process_time()\n", + "print(fast_power(a, n))\n", + "t2 = time.process_time()\n", + "print(a**n)\n", + "t3 = time.process_time()\n", + "\n", + "print(\"Time for slow_power():\", t1-t0)\n", + "print(\"Time for fast_power():\", t2-t1)\n", + "print(\"Time for Python's **: \", t3-t2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fast gcd\n", + "\n", + "Complexity: $O(\\log_2(n))$" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "126\n", + "Running time: 0.00017707599999994272\n" + ] + } + ], + "source": [ + "import time\n", + "\n", + "def gcd(a, b):\n", + " if b == 0:\n", + " return a\n", + " else:\n", + " return gcd(b, a%b)\n", + "\n", + "t0 = time.process_time()\n", + "print(gcd(155275387236018, 572335397352432))\n", + "t1 = time.process_time()\n", + "\n", + "print(\"Running time:\", t1-t0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fibonacci numbers\n", + "\n", + "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", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "9227465\n", + "9227465\n", + "Time for F_slow: 2.3301570169999906\n", + "Time for F_fast: 8.848800000293977e-05\n" + ] + } + ], + "source": [ + "import time\n", + "\n", + "F_memorized = [-1] * (10**6)\n", + "\n", + "def F_slow(n):\n", + " if n <= 1:\n", + " return n\n", + " else:\n", + " return F_slow(n-1) + F_slow(n-2)\n", + " \n", + "def F_fast(n):\n", + " if F_memorized[n] == -1:\n", + " if n <= 1:\n", + " F_memorized[n] = n\n", + " else:\n", + " F_memorized[n] = F_fast(n-1) + F_fast(n-2)\n", + " \n", + " return F_memorized[n]\n", + "\n", + "n = 35\n", + "\n", + "t0 = time.process_time()\n", + "print(F_slow(n))\n", + "t1 = time.process_time()\n", + "print(F_fast(n))\n", + "t2 = time.process_time()\n", + "\n", + "print(\"Time for F_slow:\", t1-t0)\n", + "print(\"Time for F_fast:\", t2-t1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} 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 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Diffie-Hellman key exchange\n", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Public key: p = 20747 and g = 13428 \n", + "\n", + "[[ Alice's secret key: a = 12403 ]]\n", + "[[ Bob's secret key: b = 17642 ]] \n", + "\n", + "Alice sends h1 = 14710 to Bob\n", + "Bob sends h2 = 10680 to Alice \n", + "\n", + "Alice computed 10455 using h2 and her secret a\n", + "Bob computed 10455 using h1 and his secret b\n" + ] + } + ], + "source": [ + "# Public information:\n", + "p = Primes()[10^3 + randint(1,10000)] # random prime\n", + "g = randint(2, p-1) # random integer\n", + "\n", + "print(\"Public key: p =\", p, \"and g =\", g, \"\\n\")\n", + "\n", + "a = randint(2, p-1) # Only Alice knows this\n", + "b = randint(2, p-1) # Only Bob knows this\n", + "\n", + "print(\"[[ Alice's secret key: a =\", a, \"]]\")\n", + "print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n", + "\n", + "h1 = (g^a) % p # Alice sends this to Bob\n", + "h2 = (g^b) % p # Bob sends this to Alice\n", + "\n", + "print(\"Alice sends h1 =\", h1, \"to Bob\")\n", + "print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n", + "\n", + "secret_a = (h2^a) % p # Alice can compute this because she knows a\n", + "secret_b = (h1^b) % p # Bob can compute this because he knows b\n", + "\n", + "print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n", + "print(\"Bob computed\", secret_b, \"using h1 and his secret b\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## General Diffie-Hellman\n", + "\n", + "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$." + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Public key:\n", + "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", + "g = (155 : 60 : 1) \n", + "\n", + "[[ Alice's secret key: a = 141 ]]\n", + "[[ Bob's secret key: b = 158 ]] \n", + "\n", + "Alice sends h1 = (29 : 125 : 1) to Bob\n", + "Bob sends h2 = (60 : 59 : 1) to Alice \n", + "\n", + "Alice computed (109 : 94 : 1) using h2 and her secret a\n", + "Bob computed (109 : 94 : 1) using h1 and his secret b\n" + ] + } + ], + "source": [ + "def genericDH(G):\n", + " if G.cardinality() == 1:\n", + " print(\"Group is trivial, can't do anything\")\n", + " return\n", + " g = G.random_element()\n", + " while g == G.identity(): # Make sure g is not trivial\n", + " g = G.random_element()\n", + " \n", + " print(\"Public key:\\nG =\", G, \"\\ng =\", g, \"\\n\")\n", + " \n", + " a = randint(2, G.exponent()-1) # Only Alice knows this\n", + " b = randint(2, G.exponent()-1) # Only Bob knows this\n", + "\n", + " print(\"[[ Alice's secret key: a =\", a, \"]]\")\n", + " print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n", + " \n", + " # \"Ternary operator\", I did not explain this\n", + " # https://docs.python.org/3/reference/expressions.html#conditional-expressions\n", + " h1 = g^a if G.is_multiplicative() else a*g # Alice sends this to Bob\n", + " h2 = g^b if G.is_multiplicative() else b*g # Bob sends this to Alice\n", + "\n", + " print(\"Alice sends h1 =\", h1, \"to Bob\")\n", + " print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n", + " \n", + " secret_a = h2^a if G.is_multiplicative() else a*h2 # Alice can compute this because she knows a\n", + " secret_b = h1^b if G.is_multiplicative() else b*h1 # Bob can compute this because he knows b\n", + "\n", + " print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n", + " print(\"Bob computed\", secret_b, \"using h1 and his secret b\")\n", + " \n", + "E = EllipticCurve(GF(157), [1,-1])\n", + "G = E.abelian_group()\n", + "genericDH(G)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Numerical methods for PDEs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.0", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} 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 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Diffie-Hellman key exchange\n", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Public key: p = 75521 and g = 58258 \n", + "\n", + "[[ Alice's secret key: a = 22794 ]]\n", + "[[ Bob's secret key: b = 69773 ]] \n", + "\n", + "Alice sends h1 = 31067 to Bob\n", + "Bob sends h2 = 54398 to Alice \n", + "\n", + "Alice computed 30031 using h2 and her secret a\n", + "Bob computed 30031 using h1 and his secret b\n" + ] + } + ], + "source": [ + "# Public information:\n", + "p = Primes()[10^3 + randint(1,10000)] # random prime\n", + "g = randint(2, p-1) # random integer\n", + "\n", + "print(\"Public key: p =\", p, \"and g =\", g, \"\\n\")\n", + "\n", + "a = randint(2, p-1) # Only Alice knows this\n", + "b = randint(2, p-1) # Only Bob knows this\n", + "\n", + "print(\"[[ Alice's secret key: a =\", a, \"]]\")\n", + "print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n", + "\n", + "h1 = (g^a) % p # Alice sends this to Bob\n", + "h2 = (g^b) % p # Bob sends this to Alice\n", + "\n", + "print(\"Alice sends h1 =\", h1, \"to Bob\")\n", + "print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n", + "\n", + "secret_a = (h2^a) % p # Alice can compute this because she knows a\n", + "secret_b = (h1^b) % p # Bob can compute this because he knows b\n", + "\n", + "print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n", + "print(\"Bob computed\", secret_b, \"using h1 and his secret b\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## General Diffie-Hellman\n", + "\n", + "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$." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Public key:\n", + "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", + "g = (53 : 90 : 1) \n", + "\n", + "[[ Alice's secret key: a = 145 ]]\n", + "[[ Bob's secret key: b = 65 ]] \n", + "\n", + "Alice sends h1 = (150 : 80 : 1) to Bob\n", + "Bob sends h2 = (4 : 58 : 1) to Alice \n", + "\n", + "Alice computed (28 : 28 : 1) using h2 and her secret a\n", + "Bob computed (28 : 28 : 1) using h1 and his secret b\n" + ] + } + ], + "source": [ + "def genericDH(G):\n", + " if G.cardinality() == 1:\n", + " print(\"Group is trivial, can't do anything\")\n", + " return\n", + " g = G.random_element()\n", + " while g == G.identity(): # Make sure g is not trivial\n", + " g = G.random_element()\n", + " \n", + " print(\"Public key:\\nG =\", G, \"\\ng =\", g, \"\\n\")\n", + " \n", + " a = randint(2, G.exponent()-1) # Only Alice knows this\n", + " b = randint(2, G.exponent()-1) # Only Bob knows this\n", + "\n", + " print(\"[[ Alice's secret key: a =\", a, \"]]\")\n", + " print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n", + " \n", + " # \"Ternary operator\", I did not explain this\n", + " # https://docs.python.org/3/reference/expressions.html#conditional-expressions\n", + " h1 = g^a if G.is_multiplicative() else a*g # Alice sends this to Bob\n", + " h2 = g^b if G.is_multiplicative() else b*g # Bob sends this to Alice\n", + "\n", + " print(\"Alice sends h1 =\", h1, \"to Bob\")\n", + " print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n", + " \n", + " secret_a = h2^a if G.is_multiplicative() else a*h2 # Alice can compute this because she knows a\n", + " secret_b = h1^b if G.is_multiplicative() else b*h1 # Bob can compute this because he knows b\n", + "\n", + " print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n", + " print(\"Bob computed\", secret_b, \"using h1 and his secret b\")\n", + " \n", + "E = EllipticCurve(GF(157), [1,-1])\n", + "G = E.abelian_group()\n", + "genericDH(G)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Numerical methods for differential equations\n", + "\n", + "## Euler's method (ODE)\n", + "\n", + "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:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "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", + "text/plain": [ + "Graphics object consisting of 2 graphics primitives" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "var('y')\n", + "\n", + "def euler_desolve(f, x0, y0, x1):\n", + " n = 5\n", + " h = (x1-x0)/n\n", + " S = []\n", + " Y = [y0]\n", + " for i in range(n+1):\n", + " S.append(x0 + i*h)\n", + " Y.append(N( Y[i] + h*f(S[i], Y[i]) ))\n", + " return S, Y\n", + "\n", + "f(x,y) = y\n", + "x0 = -1\n", + "x1 = 2\n", + "y0 = e^(-1)\n", + "\n", + "S, Y = euler_desolve(f, x0, y0, x1)\n", + "plot(e^x, -1, 2) + line([(S[i], Y[i]) for i in range(len(S))], color='red', marker='o', markersize=2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Sage also has an `eulers_method()` function \"for pedagogical purposes only\":" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " x y h*f(x,y)\n", + " -1 0.367879441171442 0.0367879441171442\n", + "-0.900000000000000 0.404667385288587 0.0404667385288587\n", + "-0.800000000000000 0.445134123817445 0.0445134123817445\n", + "-0.700000000000000 0.489647536199190 0.0489647536199190\n", + "-0.600000000000000 0.538612289819109 0.0538612289819109\n", + "-0.500000000000000 0.592473518801020 0.0592473518801020\n", + "-0.400000000000000 0.651720870681122 0.0651720870681122\n", + "-0.300000000000000 0.716892957749234 0.0716892957749234\n", + "-0.200000000000000 0.788582253524157 0.0788582253524157\n", + "-0.100000000000000 0.867440478876573 0.0867440478876573\n", + "-1.38777878078145e-16 0.954184526764230 0.0954184526764230\n", + "0.0999999999999999 1.04960297944065 0.104960297944065\n", + "0.200000000000000 1.15456327738472 0.115456327738472\n", + "0.300000000000000 1.27001960512319 0.127001960512319\n", + "0.400000000000000 1.39702156563551 0.139702156563551\n", + "0.500000000000000 1.53672372219906 0.153672372219906\n", + "0.600000000000000 1.69039609441897 0.169039609441897\n", + "0.700000000000000 1.85943570386086 0.185943570386086\n", + "0.800000000000000 2.04537927424695 0.204537927424695\n", + "0.900000000000000 2.24991720167165 0.224991720167165\n", + "1.00000000000000 2.47490892183881 0.247490892183881\n", + "1.10000000000000 2.72239981402269 0.272239981402269\n", + "1.20000000000000 2.99463979542496 0.299463979542496\n", + "1.30000000000000 3.29410377496746 0.329410377496746\n", + "1.40000000000000 3.62351415246420 0.362351415246420\n", + "1.50000000000000 3.98586556771062 0.398586556771062\n", + "1.60000000000000 4.38445212448168 0.438445212448168\n", + "1.70000000000000 4.82289733692985 0.482289733692985\n", + "1.80000000000000 5.30518707062284 0.530518707062284\n", + "1.90000000000000 5.83570577768512 0.583570577768512\n", + "2.00000000000000 6.41927635545363 0.641927635545363\n" + ] + } + ], + "source": [ + "# Usage: eulers_method(f, x0, y0, h, x1)\n", + "eulers_method(f, -1, N(e^(-1)), 0.1, 2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving the heat equation with a finite difference method" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def heat_fdm(u0j, u1j, ui0):\n", + " m, n = len(u0j)-1, len(ui0)-1\n", + " k, h = 1/m, 1/n\n", + " \n", + " u = [[0] * (m+1) for i in range(n+1)]\n", + " for j in range(m+1):\n", + " u[0][j] = u0j[j]\n", + " for j in range(m+1):\n", + " u[n][j] = u1j[j]\n", + " for i in range(n+1):\n", + " u[i][0] = ui0[i]\n", + " \n", + " for j in range(0,m):\n", + " for i in range(1,n):\n", + " u[i][j+1] = (k/(h*h)) * (u[i+1][j] - 2*u[i][j] + u[i-1][j]) + u[i][j]\n", + " \n", + " return u\n", + "\n", + "n, m = 20, 20\n", + "u0j = [10 - (j/m)*10 for j in range(m+1)] # One extreme goes from hot to cold\n", + "u1j = [(j/m)*10 for j in range(m+1)] # The other does the opposite\n", + "ui0 = [10 - (i/m)*10 for i in range(0,n+1)]\n", + "\n", + "u = heat_fdm(u0j, u1j, ui0)\n", + "for t in range(m+1):\n", + " show(line([(i/n, u[i][t]) for i in range(n+1)], ymin=-1, ymax =12))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.0", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} -- cgit v1.3