{ "cells": [ { "cell_type": "markdown", "id": "ac717206", "metadata": {}, "source": [ "---\n", "format:\n", " html:\n", " other-links:\n", " - text: This notebook\n", " href: L22-23.ipynb\n", " - text: preamble.jl\n", " href: preamble.jl\n", "---\n", "\n", "\n", "# Pivotted LU Decompositon\n", "\n", "::: {.callout-note}\n", "\n", "These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes. \n", "\n", "::: " ] }, { "cell_type": "code", "execution_count": 1, "id": "193b2d8d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "✓ file included! \n", "\n", "using: Plots, LaTeXStrings, Polynomials, PrettyTables, LinearAlgebra \n", "\n", "functions included: \n", " Ch2: simple_iteration, Newton, orderOfConvergence, \n", " Ch3: ChebyshevNodes, ChebyshevInterpolant, Lagrange, \n", "\n", "Use @doc <> for help\n" ] } ], "source": [ "# | code-fold: true\n", "\n", "include(\"preamble.jl\")" ] }, { "cell_type": "markdown", "id": "c06048a5", "metadata": {}, "source": [ "::: {.callout-note}\n", "## Warm-up\n", "\n", "Use LU decomposition combined with forward and backwards substitution to solve the linear system of equations\n", "\n", "\\begin{align}\n", "\\begin{pmatrix}\n", "2 & 3 & 4 \\\\\n", "4 & 5 & 10 \\\\\n", "4 & 8 & 2\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "x\\\\y\\\\z\n", "\\end{pmatrix}\n", "= \\begin{pmatrix}\n", "1 \\\\ 0 \\\\ -1\n", "\\end{pmatrix}\n", "\\end{align}\n", "\n", ":::" ] }, { "cell_type": "code", "execution_count": 2, "id": "ca2fbaa7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LU{Float64, Matrix{Float64}, Vector{Int64}}\n", "L factor:\n", "3×3 Matrix{Float64}:\n", " 1.0 0.0 0.0\n", " 1.0 1.0 0.0\n", " 0.5 0.166667 1.0\n", "U factor:\n", "3×3 Matrix{Float64}:\n", " 4.0 5.0 10.0\n", " 0.0 3.0 -8.0\n", " 0.0 0.0 0.333333" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "A = [2 3 4 ; 4 5 10 ; 4 8 2]\n", "lu(A)" ] }, { "cell_type": "markdown", "id": "b2a2da6a", "metadata": {}, "source": [ "
Exercise. \n", "\n", "What happens when you apply LU decomposition to the following matrix? \n", "\n", "\\begin{align}\n", "A = \\begin{pmatrix} 2 & 0 & 4 & 3 \\\\ -2 & 0 & 2 & -13 \\\\ 1 & 15 & 2 & -4.5 \\\\ -4 & 5 & -7 & -10 \\end{pmatrix}.\n", "\\end{align}\n", "\n", "The first row of $U$ is $u_1 = (2, 0, 4, 3)$ and so the first row of $L$ is $\\ell_1 = (1, -1, 1/2, -2)^\\intercal$ (the first column of $A$ divided by $(u_1)_1$). Therefore, we have \n", "\n", "\\begin{align}\n", "A^{(2)} := A - \\ell_1 u_1 &= \\begin{pmatrix} 2 & 0 & 4 & 3 \\\\ -2 & 0 & 2 & -13 \\\\ 1 & 15 & 2 & -4.5 \\\\ -4 & 5 & -7 & -10 \\end{pmatrix}\n", "%\n", "- \\begin{pmatrix} 2 & 0 & 4 & 3 \\\\ -2 & 0 & -4 & -3 \\\\ 1 & 0 & 2 & 3/2 \\\\ -4 & 0 & -8 & -6 \\end{pmatrix} \\nonumber\\\\\n", "%\n", "&= \\begin{pmatrix} 0 & 0 & 0 & 0 \\\\ 0 & 0 & 6 & -10 \\\\ 0 & 15 & 0 & -6 \\\\ 0 & 5 & 1 & -4 \\end{pmatrix}\n", "\\end{align}\n", "\n", "Now we have a problem! In the algorithm we learnt last time, we would set $u_2 = (0, 0, 6, -10)$ and then let $\\ell_2$ be the second column of $A^{(2)}$ divided by $(u_2)_2$ but this is zero! Instead, we will use a pivoting strategy. This is what the built in function does:\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": 23, "id": "bb0f5362", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LU{Float64, Matrix{Float64}, Vector{Int64}}\n", "L factor:\n", "4×4 Matrix{Float64}:\n", " 1.0 0.0 0.0 0.0\n", " -0.25 1.0 0.0 0.0\n", " 0.5 -0.153846 1.0 0.0\n", " -0.5 0.153846 0.0833333 1.0\n", "U factor:\n", "4×4 Matrix{Float64}:\n", " -4.0 5.0 -7.0 -10.0\n", " 0.0 16.25 0.25 -7.0\n", " 0.0 0.0 5.53846 -9.07692\n", " 0.0 0.0 0.0 -0.166667" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "A = [ 2 0 4 3 ; -2 0 2 -13 ; 1 15 2 -9/2 ; -4 5 -7 -10 ]\n", "\n", "display( lu(A) )\n", "\n", "# other examples (we went through some of these)\n", "B = [ 1 2 -3 ; 0 1 -1 ; 1 1 0 ]\n", "C = [0 6 -10 ; 15 0 -6 ; 5 1 -4 ]\n", "D = [ 4 0 8 6 ; -8 10 -14 -20 ; 2 30 4 -9 ; -4 0 4 -26 ];\n", "E = [2 3 4 ; 4 5 10 ; 4 8 2 ]\n", "E \\ [6 ; 12 ;7] ;\n", "F = [-2 1 0 0 ; 1 -2 1 0 ; 0 1 -2 1 ; 0 0 1 -2]\n", "lu(F);" ] }, { "cell_type": "markdown", "id": "80ec624e", "metadata": {}, "source": [ "## Pivoting\n", "\n", "Before, we chose $u_k$ to be the first nonzero row of $A^{(k)}$. Now we instead choose the $p(k)$ row so that we maximise the absolute value of $[A^{(k)}]_{p(k),k}$. Let's consider an example:\n", "\n", "Use LU decomposition to solve $Ax = b$ where\n", "\n", "\\begin{align}\\nonumber\n", "A = \\begin{pmatrix}\n", "1 & 0 & 2 \\\\ 0 & 2 & 1 \\\\ 2 & 0 & 3\n", "\\end{pmatrix}.\n", "\\end{align}\n", "\n", "First, we choose $p(1)$ so that $A_{p(1),1}$ is maximised: that is, $p(1) = 3$. Then, we let $u_1$ be the $p(1)$ row of $A$: that is, $u_1 = (2, 0, 3)$ and $\\ell_1$ be the first column of $A$ divided by $(u_1)_1=2$: that is, $\\ell_1 = ( 1/2, 0, 1 )^\\intercal$. At this stage we have \n", "\n", "\\begin{align}\\nonumber\n", "L = \\begin{pmatrix}\n", "1/2 & \\bullet & \\bullet \\\\ 0 & \\bullet & \\bullet \\\\ 1 & \\bullet & \\bullet \n", "\\end{pmatrix}, \n", "\\qquad \n", "U = \\begin{pmatrix}\n", "2 & 0 & 3 \\\\ \\bullet & \\bullet & \\bullet \\\\ \\bullet & \\bullet & \\bullet \n", "\\end{pmatrix}, \n", "\\qquad \n", "p(1) = 3.\n", "\\end{align}\n", "\n", "Notice that $L$ is no longer going to be lower triangular, but we will fix that at the end (this is why we are tracking the permutation $p$).\n", "\n", "Then we consider the matrix $A^{(2)} := A - \\ell_1 u_1$ and repeat the same proceedure:\n", "\n", "\\begin{align}\\nonumber\n", "A^{(2)} = \\begin{pmatrix}\n", "1 & 0 & 2 \\\\ 0 & 2 & 1 \\\\ 2 & 0 & 3\n", "\\end{pmatrix} - \n", "\\begin{pmatrix}\n", "1 & 0 & 3/2 \\\\ 0 & 0 & 0 \\\\ 2 & 0 & 3\n", "\\end{pmatrix}\n", "%\n", "= \\begin{pmatrix}\n", "0 & 0 & 1/2 \\\\ 0 & 2 & 1 \\\\ 0 & 0 & 0\n", "\\end{pmatrix}.\n", "\\end{align}\n", "\n", "We choose $p(2)$ to maximise $A^{(2)}_{p(2),2}$ and so we set $p(2) = 2$ and $u_2 = (0, 2, 1)$. Then we set $\\ell_2$ to be the second column of $A^{(2)}$ divided by $(u_2)_2 = 2$ (that is, $\\ell_2 = (0, 1, 0)^\\intercal$) and we have \n", "\n", "\\begin{align}\\nonumber\n", "L = \\begin{pmatrix}\n", "1/2 & 0 & \\bullet \\\\ 0 & 1 & \\bullet \\\\ 1 & 0 & \\bullet \n", "\\end{pmatrix}, \n", "\\qquad \n", "U = \\begin{pmatrix}\n", "2 & 0 & 3 \\\\ 0 & 2 & 1 \\\\ \\bullet & \\bullet & \\bullet \n", "\\end{pmatrix}, \n", "\\qquad \n", "p(1) = 3, p(2) = 2\n", "\\end{align}\n", "\n", "Then we consider the matrix $A^{(3)} := A^{(2)} - \\ell_2 u_2$ and repeat the same proceedure:\n", "\n", "\\begin{align}\\nonumber\n", "A^{(3)} = \\begin{pmatrix}\n", "0 & 0 & 1/2 \\\\ 0 & 2 & 1 \\\\ 0 & 0 & 0\n", "\\end{pmatrix} - \n", "\\begin{pmatrix}\n", "0 & 0 & 0 \\\\ 0 & 2 & 1 \\\\ 0 & 0 & 0\n", "\\end{pmatrix}\n", "%\n", "= \\begin{pmatrix}\n", "0 & 0 & 1/2 \\\\ 0 & 0 & 0 \\\\ 0 & 0 & 0\n", "\\end{pmatrix}\n", "\\end{align}\n", "\n", "and conclude \n", "\n", "\\begin{align}\\nonumber\n", "L = \\begin{pmatrix}\n", "1/2 & 0 & 1 \\\\ 0 & 1 & 0 \\\\ 1 & 0 & 0 \n", "\\end{pmatrix}, \n", "\\quad \n", "U = \\begin{pmatrix}\n", "2 & 0 & 3 \\\\ 0 & 2 & 1 \\\\ 0 & 0 & 1/2 \n", "\\end{pmatrix}, \n", "\\quad \n", "p(1) = 3, p(2) = 2, p(3)=1\n", "\\end{align}\n", "\n", "Now, we want to solve $A x = b$. Notice that \n", "\n", "\\begin{align}\\nonumber\n", "PL = \\big( L_{p(i),j} \\big) = \\begin{pmatrix}\n", "1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 1/2 & 0 & 1 \n", "\\end{pmatrix}\n", "\\end{align}\n", "\n", "is lower triangular (where $P$ is the permutation matrix with a $1$ in the $(i,p(i))$ entry), and so we have \n", "\n", "\\begin{align}\n", "(PL) U = Pb\n", "\\end{align}\n", "\n", "which we may solve using forwards and backwards substitution as before.\n", "\n", "
Example. \n", "\n", "Use the PLU decomposition as above to solve \n", "\n", "\\begin{align}\n", "\\begin{pmatrix}\n", "1 & 0 & 2 \\\\ 0 & 2 & 1 \\\\ 2 & 0 & 3\n", "\\end{pmatrix} \\begin{pmatrix} x\\\\y\\\\z\\end{pmatrix} = \\begin{pmatrix}0\\\\1\\\\-1\\end{pmatrix}\n", "\\end{align}\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "44e7b0f3", "metadata": {}, "source": [ "## Why use the largest pivot?\n", "\n", "Answer: We get better conditioning by choosing the largest pivot. Later, we introduce the condition number of a matrix, written $\\kappa(A)$, which will make the following discussion rigorous. \n", "\n", "Suppose we have the matrix \n", "\n", "\\begin{align}\n", "M(\\epsilon )= \\begin{pmatrix}\n", "-\\epsilon & 1 \\\\\n", "1 & -1\n", "\\end{pmatrix}\n", "\\end{align}\n", "\n", "with $\\epsilon \\ll 1$. \n", "\n", "1. It turns out that this matrix is well condtitioned for all small $\\epsilon$. Here, we can see the condition number is bounded for all small $\\epsilon$:" ] }, { "cell_type": "code", "execution_count": 5, "id": "1710afb3", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dZ0AU194G8DOzSxekd+lGRFAQNKJgRaOJJbYYTTMaJSaWFK+aa4p6E29evSYaTUyMuYm9xU6isaOAKCg2QFGkSa9KZ3fmvB8m2ewFlCK7s+X5fWKHs7v/PTu7z56ZMzMMpZQAAADoK1bsAgAAAMSEIIQOU1ZWdv369djY2JSUFI7j1F/AkSNHIiMjL1++rFiyYsWKmTNn1tXVtXjf999//91331VldZrrxo0bkZGRu3fvFruQFkRHR3/++edz5syJjIxMS0sTuxzQIRTg6fA8v3Pnzj59+rDs37+rOnXq9Oqrr964cUOdlXzyySeEkB07diiWBAYGEkIePnwo3Kytrf3hhx8OHDjQ9L52dnadOnVSU6Ea5siRI4SQefPmiV3Ik3zxxRfCqmVhYWFlZRUTE9Nss6VLlwrNRo8e3fS/MpnM2dlZaHDw4EEVlwxaQ6rO0AXdU19f//rrr+/du5dl2YiIiH79+pmampaWlkZHR2/fvn3Xrl0PHjxwdHQUq7w+ffpYWlpKpX+u51VVVZGRkSEhIePHj2/UMiwsrL6+Xu0FQqtwHPfFF19YW1snJiZ6enq22J5hmOPHj+fn5zs5OSkvP3HiRF5eHsMwFHMjQAmCEJ7KggUL9u7d6+HhcejQoV69ein/Kz4+fvbs2TKZTKzaCCGbNm1qZcsDBw6otBJ4Gjk5OTU1NQMHDmxNChJCBgwYEBMTs2PHjoULFyov//nnn6VSaUhISHx8vGoqBa2EIIT2S0hI+OGHH0xMTI4dO+br69vov/369bt48aLy9lJCSG5u7pkzZwoKCqytrcPDw5955hnl/9bU1KSmplpaWnp7excVFR0/fryoqMjLy2vkyJGmpqaNHp9SeuHChaSkJCMjo/Dw8B49ejStMDU1taamJigoiGXZwsLCW7duEUKqq6uvXLkiNLCxsfHw8CCE3Lhxg+d5YVNqo9eYmJhYW1vbpUuX4cOHW1paNno5BQUFXl5eVlZWSUlJsbGxPM/36dMnNDS0xd4TauvVqxfDMCdOnLh9+7aZmVlERISXl5dyszt37lRVVfXs2dPAwECxsL6+/tatWxYWFl27dhWW5OXl5efne3h42NjYJCYmxsfHGxgYDB48uFu3bkKDoqKiP/74o6ioyNfXd+TIkRKJpNmq0tLSzp07V1NT07Nnz8GDBzd6+wQFBQWnT5/Oz8+3sLAICwvz8/NT/m9paWlmZqaTk5Ozs/Pdu3ejo6MrKiqmTZum2CbZFKX08uXLiYmJ9fX17u7uERERnTt3Vvz36tWrt2/fFl618MaZm5s3WnMaGTVqVHp6+i+//KIchGVlZUePHh05cqShoeET7gv6SORNs6DNZsyYQQiZNWtWaxpzHLdkyRLlb3OGYV555ZWamhpFm6tXrxJCxo8fv2XLFhMTE0VLd3f3tLQ05UerqKgYOnSo8pr89ttvC/uHHrePcM2aNU3X/9dff11o2XQfYW5ublhYmHJjCwuLzZs3K7dZvHgxIWTr1q0vv/yycsvJkyfL5fInd0hISAghJDExsWfPnoo7SqXSdevWKTfr378/IaSgoEB54b179wghI0aMUCz57LPPCCE//vjjpEmTFI/GsuxXX31FKd20aZORkZFieXh4eHV1teK+wj7CuXPnLly4kGEYRbMBAwYUFxcrP69MJvvggw+U30RCyMsvv6z8aD///DMh5NNPP503b57i0U6dOvW4fsjOzu7Xr5/yA1paWm7dulXRQLlyweDBgx/3aMI68OWXXwoRmJCQoPjX+vXrCSH79u2bMGECwT5CUIIghPYTtlMdOXKkNY0//vhjQoinp+evv/6anp5+8uTJ4OBgQsjEiRMVbYQgdHNzMzMzW7FiRVxc3JkzZ0aPHt30u++FF14QkuDSpUvZ2dnbt2+3s7MTxhyPC8KsrKxff/2VEPLMM8+c/MutW7eElo2CsLa21t/fnxAyfvz4+Pj4tLS077//3sLCghCyd+9eRTMhCD09PX18fLZv33716tXdu3e7uroSQr777rsnd4gQhJ6eniNHjoyKikpMTFy1apWRkZFEIklJSVE0a1MQuru7d+3adffu3VevXt24caOZmRnLsmvXrjUxMVm5cuWlS5eOHz8ubMFetmyZ4r5CEDo7O1tZWf3888/Z2dkJCQljxowhhAwaNIjneUVL4adPcHDwwYMHU1NTz5w589xzzxFCpk6dqmgjBKGbm5u9vf2aNWuio6OjoqLS09Ob7YSqqiphzPrSSy9dvnz5zp0769ev79SpE8MwivXq9OnTGzduFIJZeNcSExMf16uKIBRG/3PnzlX8Kzg42Nrauq6uDkEIjSAIoZ3kcrmw3ezu3bstNs7Pzzc0NJRKpcoDu4qKCjs7O0KIYgagEISNwqaurk5IuKKiImFJdHQ0IcTb27u+vl7R7MSJE8J9nzBrtLi4mBASEhLStMJGQfjtt98SQvr168dxnGLh/v37hbBRjPaEILS1tS0tLW1UybBhw57cJ0IQjhgxQjlphAdcuXKlYkmbgtDa2rqkpESxcNmyZUKf7NmzR7EwJSWFYRh/f3/FEiEICSH79+9XLJTJZMJPgaioKGHJhQsXCCE9e/asq6tTNOM4rk+fPoSQpKQkYYkQhAzDxMfHP7kH6F/D9MGDByt3wo4dO4TfK4qFik0FLT6gIggppSEhIULyUUpv3rypyEUEITSC4wihnSorK3meJ4SYm5u32DgqKqqhoWHSpEmKfVqEkM6dO7/zzjuEECFgFDw8PCZPnqy4aWRkNGjQIEJIRkaGsOTQoUOEkHnz5inv7Bk+fHjTPXztdvDgQULIokWLlHeSjR8/3tfXNysrS7GLUfDWW29ZW1srbg4dOlQqlSqqfbJGWyNHjBhBlF5pW73xxhs2NjaKmwMHDiSEuLi4KPdn9+7d7e3tmz6Ft7f3iy++qLgplUoXLFhA/uoKQsi2bdsIIYsWLVLeVsmybGRkJCHk2LFjyo82ePDgZ599tsWChTlKixcvVu6EKVOmeHl5paWlCenVbtOnTxf2CxJChHiePn360zwg6CpMloF2MjMzE6ah19bWttg4JSWFEBIUFNRoubB1NDk5WXmhYn6HgoODAyGksLBQ+dEaTVIlhAQGBl67dq31L6HFgnv37q28kGGY3r173759Ozk5uW/fvorljSZuSCQSW1vbgoKC1jxRoxcrvNJW3rfFRxMG3F27dlWOGWF5YWFhZWWl8o+YgICARlNjhB8WincnKSmJEHL69OlG+ZSVlUUIyczMVF7YaAbN4zTbzxKJJDAw8P79+8nJyco7UNtq6tSpH3744ZYtW1588cWdO3f26NFDWN8AGkEQQjsZGBjY29sXFhamp6cLEy+foKqqivz1La/M3t6eEFJZWam8sOkEUeELWhiAKh5N+JZX1nRJuz2uYGFJawqmrTtSrdF9hVfayvs2pTzDSPFoj+vPRs8ivBdNlyhebEVFBSHk4MGDjWKVEGJlZdVooa2tbWsKrqqqYlm2aeNm+7mtrK2tx44de/Dgwc2bNxcUFDQ6lAJAAZtGof3Cw8MJISdPnmyxpTDyaDrQEQZ5wiSU1uvUqRMhpKioqNHypkvaTShYMQZVEF5CWwt+GkLAKH4ECKqrqzv8iZq+2EbvjtDtZ8+eLWvO999/344nNTc353le2HerrKP6efr06XK5/IMPPpBKpa+++upTPhroKgQhtJ+wx+Wnn35q+kUmoJQKJx0Vpl002rVGCElISCCEBAQEtOl5hUcTttQpa/r4jQjz/uVyeYtPIRyVmJiYqLyQ53nhKdpa8NMQzo3SKKWELYod6/r1643iVpiiIvQ2+WvLdlxcXAc+qfDgjfqZ4zjhqZ++n0eMGOHs7FxbWztq1Kim43sAAYIQ2u/5558fPnx4SUnJpEmTysvLG/23uLh48uTJeXl5hJDRo0cbGxsfPHjwzp07igbl5eXCtHjlqRytIcz627Bhg/LZtI8dOybMmH+Czp07d+rUKS8vr9E3flPC0XirV69WPnv4r7/+mpaW5u3t3XRnp+oIx6j8/vvviiVyufw///lPhz9RZmamcHiJoKGhYd26dYSQiRMnCkvefPNNQsjq1atLS0sb3Vcul7fmzOZNCf28atUq5Xdkx44dWVlZfn5+rdzR+ARSqfTYsWMnT54UDiIEaBaCENqPYZgdO3YEBgaeP3/+mWeeWbx48YEDB06dOrVjx47IyEhvb2/FdFA7O7ulS5fK5fIRI0bs3LkzOTn56NGjgwYNKi0tfeWVV5QnnrRG//79x48fn5GRMWrUqLNnz965c2fz5s2vvPKKm5tbi/cNCQkpKip66aWX1q1bt2nTprNnzzbb7I033ggKCrp06dK4cePOnj2bnJz89ddfz5w5k2GYNWvWNHu+FRWZMmUKy7Kff/756tWrY2Njd+/e3b9//0ePHnX4E7m5uc2aNeu77767ffv2+fPnx4wZk5qaOnz4cGEiKyEkLCzs3XffzczMDAkJ2bBhQ0xMzLVr1w4ePLhkyRI3N7cWf4U0a9asWT169Dh//vzEiROjo6Nv3bq1evXqt99+WzgVQNOdke3Qs2fPiIgId3f3p38o0FWYLANPxc7O7sKFC1988cWGDRtWrVql/K/AwMClS5d26dJFuLl06VJK6cqVK1955RVhCcuys2fP/uabb9rxvFu3bp06dWpUVNS5c+cIIQzDvPfee2ZmZp9//vmT77hx48YZM2YcPHhQCOnXX399yJAhTZsZGhr+8ccfr7322m+//fbbb78JC21sbDZt2jRu3Lh2FNxuwcHBq1evXrx48aJFi4Ql4eHhW7dubc1Z3Npk3LhxlpaWwpF2wpJhw4bt3btXOY3Wr1/v5ua2cuXKefPmKRYyDNOnT5/2zVQyNjY+efLkq6++eujQIeGoGEKIvb39tm3bhEP1AdQAZ2GHjlFfXx8fH5+VlVVTU2NtbR0UFKR8yKCCcGGKoqKizp07h4WFKWJS8SC5ublmZmaNdueUlpY+fPjQ0dGx0QTIxMTEa9euGRoaDhgwwNvbu7y8vLy83N7eXpjWQQjJycmpq6vz9vZuNIaTyWQFBQUymaxTp07C3MiMjAye5729vRsVnJKSkpCQUFtb6+HhER4ebmZmpvzfsrKyiooK5WcUZGVl8Tz/5DNE5+bmCqfWVD7tp0wmy8nJMTU1bXTJjqysrOjo6Pr6+h49eoSGhsrl8pycHBMTE8XVFZq+dkJIQ0PDgwcPmj7agwcPGhoaPDw8hG6pqakpKCiwsLCwtbXNyso6f/58XV1dz549+/bt2+yYrLq6Oi4uLjMzUyqVOjo69urVS/k8opWVlcXFxVZWVlZWVk94+Y3cvHnz6tWrtbW1Xl5e4eHhjaa/Pm7FaEp4R6ytrRudFVZZYWFhdXW1g4NDo3cT9BaCEAAA9Br2EQIAgF5DEAIAgF5DEAIAgF5DEAIAgF5DEAIAgF5DEAIAgF5DEAIAgF5DEAIAgF5DEAIAgF5DEAIAgF7r+CC8ePHiiRMnOvxhxdLi9XrgCdB7TwOnP2w3dN3T0MOPbccHYUxMzKlTpzr8YcWiikuB6w/0XrtxHNe+K/wBIaSmpgZZ2G56+LHFplEAANBrCEIAANBrCEIAANBrCEIAANBrCEIAANBrCEIAANBrCEIAANBrCEIAANAOdRz540HHHyGKIAQAAO2w9ha/5W7Hn/hG2uGP2G7vvfdeTEyM2FU0xvM8y7bz58LGjRv79OnTsfUAAOinkjry1U0ubmzHx5YGBWF8fPy8efP8/f3FLqRjfPTRRxkZGQhCAIAO8fk1booX62PBdPgja1AQEkJ8fX2Dg4PFrqJjWFlZiV0CAICOSH9Ed9zjUyYZqOLBsY8QAAA03ZIE/sMAiZ2xSh4cQQgAABrtYhG9XEwX9FBVYCEIAQBAc1FCPozn/hXMmqhsVx6CEAAANNfudF7Gk1d9VJhWmjVZBgAAQKFWTj5K4LcNlrAdP1f0bxgRAgCAhvryOjfAgQl3VGUMYkTYbqmpqb///vv169f9/PyWLFkidjkAALomp5p+l8onvqjynMKIsJ2io6NTU1PLy8tjY2PFrgUAQActvMTP6yFx76Ta4SBBELbo8uXLP/30k+LmpUuX/vvf/xJC3n777c2bNw8aNEi80gAAdFZcIb1YSBcGqCOkEIQtcHNz++CDD8rKyoSby5cvb2hoELckAADdxlOy4CK3+lnWVC277zR6H+F/0/gvkjr+RONP8IIb802oRHmJo6PjiBEjtm3btmDBguzs7Li4uN27d6uzJAAAffPTHd5QQl7yUtNQTaODcIoXO9hJ5VuHlVkaNvN0c+bMmTt37vz583/88ceXXnrJwsJCnSUBAOiVShlZdpU/PEKitm9/jQ5CMynxMldrEDZr6NChLMueP3/+559/PnTokNjlAADosuVXuee7MCG26vvyb0MQFhcXr1mz5vr167a2tu+8805oaKjqytI0b7311ptvvung4BASEiJ2LQAAOiv9Ed1yl78xQSVXmXic1m6BraioCA0NLSsrmzt37siRIysqKlRalqaZPn16QUHB7NmzFUu2bdvGMMw//vGPqKgohmHeeustEcsDANAN78Vzi3tJnEzV+qStHRGuXr3az89v06ZNKq1GYxUUFBgZGU2bNk2x5LXXXnvttddELAkAQMeczqOpFeTXYeo+nKG1z3fhwoXhw4evXLkyMjJy586dlFKVlqVRvvnmm2nTpi1YsMDc3FzsWgAAdJOcJ+9d5Nb2kxhJWm7csVo7IszKyvr3v//9wQcfDB069NNPP7179+5nn33WbMu0tLQTJ04kJiYqlqxdu9bT07PFp+B5tR4p0Xosy3700UcTJkxo070opbW1tZWVlSqqSitUVVWJXYK24jiuoaFBLpeLXYhWqqmp4TiOYcSfaqeNxPrYfpcmtTdkBlnVduy3pqmpqUTSQrS2NghNTU0HDhy4cOFCQoilpeVrr732uCB0d3fv16/frFmz/nwCqdTPz08qbfmJWFZDj+6fO3duO+7FMIyJiQkGkeiB9hGC0MTEROxCtBLLsqampgjCdlP/xza/hqxJlZ0fLTU3V81F6J+otUHo5ubm7Ows/O3i4lJeXs5xXLMxa2Rk5O7uHhER0WE1AgCATpt/kXvHj/W1FOe3S2sHYdOmTTt+/LhMJiOEHD58OCQkpMXBJgAAQIv+eECvldIlvUTLlNaOCKdNm3bkyBFfX197e/uioqIDBw6otCwAANAHtXLybhz3TajEWLyxVWuD0MDAYP/+/ZmZmXV1dT4+Pq3Z5wcAAPBkK69zIbbM813E3KHbtjzz8PBQTRmEEGJkZDRx4kQjIyPVPUU7UErbt8u9sLBQ+bhDAABo5O5D+kMqnzRe5JGVBg3sDh8+rLjakeaorq42MzNrxx0NDAxcXV07vB4AAJ0xJ5b7OEjiYiby/F4NCkJLS0tLS0uxq2issrISBwAAAHS47ff44jryTnfxD5zToCAEAAA98UhGliTwvw6TSMXPQVyhHgAA1G7JZW6cO9PPXiNOeoARIQAAqFViCT2cRW9N1JQAwogQAADUh6MkMoZb/SxrpTGHCCAIAQBAfb5N4S0NyTRvDUofTRmZAgCAzsuuop8ncRfGaFb0aFAmAwCAbouM4T4MkHTrrBFzZBQQhAAAoA4/p/EFteSDAI3LHc0anwIAgE4qqCVLErhjz0kNNC4HMSIEAADVezeWe9uX7W2rWRtFBRgRAgCAau25z6dW0B1DNDRxNLQsAADQDaX15P147tdhUhGvOPhk2DQKAAAqNC+Oe9WH7e+giRtFBQhCAABQld9y6KUi+llvTR0MEkKwaRQAAFTkYQN5J5bbNlhiptlRgxEhAACoxIeXuDFuzEBHzd0oKtDsmAYAAO10Jo+ezqM3JmhBymBECAAAHaxKRt66wH0/QGJuIHYprYAgBACADrbwEjfEiXnOVdM3igq0YNAKAABa5I8H9EQuvaYNG0UFWlMoAABovpI6MuM8t3OIxEIbNooKsGkUAAA6zDux3GtdmUFO2rFRVIARIQAAdIxf0viUCrp1sJYli5aVCwAAmulBNV2SwJ0YpbnnFH0cbBoFAICnxVPy2jnuHz0lPa21aaOoAEEIAABPa81NXk7Je/5amSnYNAoAAE8luZyuucldGieVaN9okBCMCAEA4GnUc+SVs9yqvhL3TtoZgwhCAAB4Gp9c4TzNmde7anGaYNMoAAC004UCuuMeva49J5FplhZnOAAAiKi8nrx2jvsxXGJrLHYpTwdBCAAA7TEnlpvoyTzfRVt3DSpo93gWAABE8W0Kf++R9p1Eplm68BoAAECdbpXT5Ve5C2OkhjqxVVEnXgQAAKhLtZy8dJr7qp+kW2et3ygqQBACAEAbzI/j+tkzr/roTnxg0ygAALTW3vv8hUJ65UWdyg6dejEAAKA66Y/ovIvc8ZFSc+256G5r6M7YFgAAVEfGk1fPcZ8ESYJsdGTXoAKCEAAAWrYkgbMxIu/66WBqYNMoAAC04EgWvz+DXh0v1bXBICEEQQgAAE+WVUUjY7j9EVJrI7FLUQ0dHOQCAEBHqePIhFPcp70l/R10cjRICIIQAACe4N1YrltnZk53XQ4LbBoFAIDm/XSHjy2kCbp11GBTOv7yAACgfW6U0X8mcude0LWjBpvS5dEuAAC0T0UDmXCK29Bf0t1SZ3cNKiAIAQDgf1BCZpznxroxkz31IiOwaRQAAP7Hl9f5/Bq6e6i+BIS+vE4AAGiN6EJ2QzJ/eZxEN6412BqtDcKqqqr4+HjFTV9fX1dXV9WUBAAA4siopG/FG+waJnEx0/1dgwqtDcL79++PHj06PDxcuLlgwQIEIQCALqmSkXEnuX/4yYc4GYpdi1q1YdOonZ3dyZMnVVcKAACIhRIy8wIXbMvM7sqJXYu6tSEIGxoajh49amJi0qdPn86dO6uuJgAAULPlV7kH1fTM89KGGrFLUbs2BKGtre22bdsePHhw7969gwcPDhgwoNlmxcXF8fHxK1euFG4yDDNz5kwrK6sOKFYMMplMJpOJXYW2Qu+1G8dxMplMKsV0tvYQVjyG0aO9XE/pSA758Ta5OJqwvEzHPrZSqbTFNaG1H7OAgIDk5GTh7xUrVsyZM+fGjRvNtpTL5XV1deXl5cJNhmEaGhp4nm/lE2kanue1t3jRoffajf+L2IVoJaHrEIStdPsheSeOOTCEOhgRntfHj21rg1B5lZowYcKKFSs4jpNIJE1bOjk5DR48eNWqVR1ToNgaGhqMjHT00iOqh95rN47jGIZB77WPXC43MjJCELZGWT2ZdE6+NpQd4PLn0RJ6+LFtz3Eily5d6tKlS7MpCAAA2kLOk8mn5ZM8maneenPMYHNaOyJcvnz5vXv3nnnmmczMzL17927evFmlZQEAgKq9H88ZsOTzEH0f1bT2V8D06dPDw8M5juvdu3dSUtKUKVNUWhYAAKjUlrv8yVy6e6hUovebkFs7InR3d589e7ZKSwEAAPW4WET/cYk7N1pqqV+HzjdPr7cLAwDooYxKOukUt3Ww1E8PLrHUGghCAAA98khGxp3klvRiR7oiBf+EIAQA0Bcynkw4KR/uwszrgS//v6EvAAD0AiXkrQuciZSs6qvv00QbwQmcAAD0wvKrXEo5PTca00QbQxACAOi+Xen81rv04lipGb71m0CXAADouPMFdMFF7swLUgcTsUvRSNhHCACgy1Ir6OTT8l1Dpf5W2CTaPAQhAIDOKq4jY09wq/tKhjkjBR8LQQgAoJtq5WTsCflrXdnXu+Kr/knQOwAAOoijZOpZzteS+SQI3/MtwGQZAAAdNDeOq5bTvcOk2CTaIgQhAICu+eQKl1hMz7wgNcRosBUQhAAAOuW7FH53Oo0ZIzU3ELsULYEgBADQHbvS+S+v8+dHS3DIYOshCAEAdMSpXPpePHdilNTDHHsG2wDbjwEAdMHlYjrtrPxAhLSXNVKwbRCEAABaL7mcjj0h3zJIOsABKdhmCEIAAO32oJq+8Af3f30lo7ogBdsDQQgAoMVK6sjwY9wCf/YNnD6mvdBxAADaqlJGRh2XT/Zk3vfHl3n7oe8AALRStZy88Ie8rz2zIhhXnH8qCEIAAO0jnFDb24JZH4oUfFoIQgAALdPAk8mn5TZGzOZwCYv5MU8NQQgAoE1kPJl8mjOSMDuHSCRIwY6AIAQA0BocJa9HczKe7hwikeL7u4PgFGsAANqBp2R6NFdWRw+PkBphz2DHQRACAGgBSsicWC6/hh4dITVGCnYoBCEAgKajhMyN426V0z9GSk3wtd3R0KMAAJpuyWXuagk9MUraCZcYVAEEIQCA5qKEfBjPXSigJ5/HhXZVBUEIAKChKCHvXeRiC+mJUVJLQ7Gr0V0IQgAATUQJmRfHXSulZ16QWmAsqEoIQgAAjcNT8tYF7u4jemwktoiqHIIQAECzcJTMPM9lVNJjz2F2jDogCAEANAhHyfRorriOHseREuqCbgYA0BQynkw9y9XI6aHhOGpefRCEAAAaoYEnU05zMp4ejMAZ1NQKJ20FABBfrZyMPyk3YMnB4UhBdUMQAgCI7JGMjDwutzZidg6RGOBbWe3Q5QAAYiqrJyOOyf2tmS2DcGUlcaDXAQBEk11F+x+RD3Zivu2Pa82LBkEIACCO2xU0PIqb5ct+2Qd7BcWEWaMAACK4WkJHn5B/HiKZ8QwGJCJDEAIAqNv5AjrplHzjAMlET6Sg+BCEAABqFZVNZ5yX7xgiHe6CvYIaAUEIAKA+2+7xiy9zv4+UhtgiBTUFghAAQE1W3eC/S+FPPy/tbokU1CAIQgAAleMoWXCRi86nF8ZIupghBTULghAAQLXqOfJ6NFdSR2PGSDvjQvOaBxOWAABUqKyeRByTG7Hk2EikoIZCEAIAqEpGJe1/RN7bhtkyWGKIr1tNhVZJHXcAAB9aSURBVHcGAEAlEopp2FFugT+7LlSCvYKarG1BWFFR4erqGhQUpKJqAAB0w4lcOvaE/Pswdk53jDc0XdveoQULFgQGBqqoFAAA3fB9Kj89Wn54hHSMG1JQC7ThTfr999+Li4unTZumumoAALQaR8l78dw3yfyF0dK+dtggqh1ae/jEo0ePPvzww99///3ixYtPbimXyysqKu7fv69Y4unpyTBYIQBAx1XJyLSzXB1H48ZKLTFBVHu0NggXLlw4Z84cT0/PFoPw7t27Bw4cOHnypHDTwMBgz549Pj4+T1WmeKqqqsQuQYuh99qN47iGhga5XC52IVqppqaG4zg1//7OrSFTLhgG29D/PCuT1JPKenU+eUfSsY+tqampRNLCVa5aFYQ3b948dOjQ0KFD9+3bd+nSpYqKin379o0fP14qbebu3bt3nz179qpVq9pTskYyNzcXuwQtht5rHyEITUxMxC5EK7Esa2pqqs4gjC+iE09z83uwi3uxhBir7XlVRN8+tq0KQkNDw/Hjx589e5YQcu/evcrKylOnTo0dO7bZIAQA0Ct77vPzL3I/hUtHu2EfkFZqVZJ169bthx9+EP7euXPn6tWrFTcBAPQWJWT5VW7bXXr2BakfzqOttdo8pHNzcxs+fLgqSgEA0CLVcjI9miuspZfGSW21fmuoXmtzEIaFhYWFhamiFAAAbZFRScef5HrbMtsHS41amIoBmg4HewIAtE10Ph1wVD7Vm/3vQAlSUAdgtgsAQBtsus0vu8ptHywd6oydgjoCQQgA0Cp1HImM4W6W0Ytjpe6dkIK6A5tGAQBallNNw4/KOUpixiAFdQ2CEACgBecLaL/D3CRPdvtgiSm2o+kcvKUAAE+y6Tb/2RVu+xDpMOwU1FEIQgCA5lXKyKwL3L1H9PKL0i5mSEGdhU2jAADNSK2goUfkhiw5PxopqOMQhAAAjW27x4cflb/vz27FTkE9gHcYAOBvdRxZfJn74wE9+4I0wBoDQb2AIAQA+NPdh3Tyac7bgrk0TtoZV9bVG9g0CgBACCGHs/jwKPmbz7D7IyRIQb2CESEA6Ds5Tz6+wu1Op4eHS5+1x+ZQvYMgBAC9llFJp57lHE2YaxOklhgI6iVsGgUA/XUgk+9/RD7enT0QIUEK6i2MCAFAH9XKyZIE7o8H9PeR0iAbbA7VaxgRAoDeSamgzx6Rl9eTxBeRgoARIQDoma13+UWXuTXPSl7xwUgACEEQAoD+KKkjb13g8mto7BiptwUGgvAn/CACAL1wLIf2OiDv1pnEIAXhf2FECAA6ro4jy65yu9Pp9iGSIU6IQGgMQQgAuuxmGX3lHOfbmbk6XmptJHY1oJEQhACgm3hK1ifzX1zj1vSTvIZ5MfB4CEIA0EEZlfSN85whS66Mx9UEoQX4lQQAOoUS8tM9tt8R7kV39uQopCC0DCNCANAd+TVkdoy8sEYS/YKkuxV+6EOrYEUBAB2xL4MPPCgLtmVORch8LTEQhNbCiBAAtF5BLYmM4bIq6R8jpYE2THW12AWBVsGIEAC0274MPvCArLslSXhRGogTh0LbYUQIANqqsJa8HcOlV9JjuIIEPAWMCAFA+1BCfk7jex6QBdowV3AFCXg6GBECgJZJf0QjY7iHDeTEKGkva0QgPC2MCAFAa8h5su4WH3pEPtyFjR+HFISOgREhAGiHG2X0rQucqZTEjJE+0xkRCB0GQQgAmq5WTv7vBrcxlf9XsGSWL4sMhI6FIAQAjXahgM66wPW0Zm5OMLA3Ebsa0EUIQgDQUCV1ZNFl7kwe/XaA5IUuGAeCqmCyDABoHErIT3d4//0yS0Nyc6IUKQgqhREhAGiWG2V0TixXz5GjI6R97BCBoHIIQgDQFDVysuoG930q/1EvybweLGbFgHogCAFAIxzN5ufF8QMdmZsTDeyMxa4G9AmCEABEdu8RnRfH5daQ7YMlYY4YBoK6YbIMAIimWk7+mcCFHpFHuLBXXpQiBUEUGBECgDiOZvPzL/J9bJmr46VdzBCBIBoEIQCo27VSOv8iVy0n2wdLBjggAkFkCEIAUJ/yerLsKrf7Pv/PXpK5PVgJQhA0APYRAoA6yHnyXQrv+6uMYcidyQYL/JGCoCkwIgQAlTv+gC68xDmakDPPS3tYIQBBsyAIAUCFblfQf1zm7jwkX4Swkz2xCQo0EYIQAFSitJ6suMrtSuc/DJD8Oow1kohdEMBjIAgBoIM18GRjCv/5NW6CB5s8CaeJAU2HIASADkMJ2XufX5LAB9kwcWOkXXEdedAGCEIA6Bhn8+mSyxxHyS8DJYOcEIGgNVobhElJSatWrbp//76hoWFERMSiRYtMTHCtaAAghJCUCrrsCn+5mP4zkH2rG64aAVqmtUHI8/zYsWP9/PzKy8sXLVpUUlKyfv16lVYGAJovp5p+nsQfyuI/8JdsGyzBjBjQRq0NwuDg4ODgYOHvd999d8OGDSorCQC0QFk9WXWD23ybf8uXTZts0NlQ7IIA2qsN+whlMllOTk5RUdHmzZunTp2qupoAQJNVy8k3yfxXN7mXvNhbkwwcsZMEtFwbgjAvL2/SpEk5OTleXl5Tpkx5XLOUlJSDBw/u27dPuGloaLhnzx4fH5+nrVQk1dXVDIM9Hu2E3ms3juMaGho4jhO7kP9RKyc/pUvWpkrCHeipYXJvc0o4UlUldllN1NTU8DyPda99dOxja2pqyrItnMmBoZS26UF5nv/000+PHz+emJjYbIMvv/zy/v37S5Ys+fMJGMbDw0N7u7WystLc3FzsKrQVeq/dhCDUnClpDTz5JY1fkcR3tySr+kqCbDT6E11dXW1qaqq9Xzvi0sOPbZsPn2BZdurUqV9++SXHcRJJM3vGJRKJpaWll5dXR5QHACKT8WRXOr/8Ku9lQQ4Nl4TYIl1A17Th8Inu3bsbGxs3NDRs2rQpODi42RQEAJ3BU7I/k/9nAu9hTnYNlfS1QwSCbmptEB48eHDQoEFWVlalpaXBwcHbtm1TaVkAICIhAj9J5G2NyaZwyRAcHQ86rbVBuGLFio8//ri4uNjGxsbYGKcOBNBNlJCobP7TK7yJhGzoL4lwQQSC7mvDPkJDQ0MXFxfVlQIA4jqVSxcncAYsWRHMjnHDJZNAX+BcowD6ThgFrrjKc5QsRwSC/kEQAugvOU/23Of/7wYvZcjHQex4D5wlFPQRghBAH9VzZM99/otrvJ0x+SKEHe2GCAT9hSAE0C+VMvLfO/x/bvKBNuTngZL+DkhA0HcIQgB9UVxHvk3hNiTzoQ7M4eGS3jg0HoAQgiAE0AeZlfTrW/zWu/wYN/biWFw4HuB/IAgBdNnNMrr6Bv9bDv+qD5s8SepsiggEaAxBCKCbYgro/93gEorp293Z9CkGlrheIMBjIAgBdI0QgcnlZEEPdu9Q1gSfcoAnwkcEQEfwlPyWw6+4ytfzZGEAezCCleLIeIBWQBACaL0qGdlyl//qJu9iRpYHS0Z1wYX4ANoAQQigxbKq6IZk/pe7/GAndttgHBQI0B4IQgCtdKWErrv153TQhHFSD3NEIEA7IQgBtEkDTw5n8V/f5IvryFw/dmOYgRk+xABPB58hAO1QXEf+e4ffkMK7dyLvB7ATPFgJBoEAHQFBCKDp7jyk36Xw2+7xo7uwx0ZK/K0QgAAdCUEIoKE4Sg5n8etu8XcfkXf92LsvGdgYiV0TgC5CEAJonPwa8mMq3ZwmdTXjF/izEzxYAxwRCKAyCEIADRJTQL9J5k/n8RM9mL2DuH7OxmJXBKD7EIQA4qtoIHvv898k8zwlb3RlfwgzsJByDQ1U7LoA9AKCEEBMV0roptv8vgw+wpld208S4fLnRBiOE7cuAD2CIAQQwSMZ2Z3Of5fC13JkxjNs2mQDW2wEBRAJghBArZSHgP95VjLMBecFBRAZghBAHeo4cjSbX3eLz6wir/owNyZIXc2QgAAaAUEIoFo3yuiPt/md6Xy4I7s0UPKcK8MiAQE0CYIQQCXK6snOdP6XNL64jrz5DHMdQ0AATYUgBOhIPCVn8ujWu3xUDh/hzH7Wmx3ligvkAmg0BCFAx0h7SHem87+kUVtj8poP+3UozogGoB0QhABP5ZGMHMrkt93jk8vpZE/28AhJL2tsAgXQJghCgPbgKYkrpNvu8b9m8P0dmNm+7IvuOCMogFZCEAK0zd2HdEc6v+Uu7WxA3nyGvYNj4QG0HIIQoFUKasmedH5nOp9dRad4swciJEE22AQKoAsQhABPUiMnv+XwW+/y5/Ppc67sx0HsSFdsAgXQKQhCgGbUc+RELr/vPv0th+/vwEz2ZHcPZc3wcQHQRfhkA/xNmAKzL4PffZ/vasFM9mTX9DOwwy5AAJ2GIAQghJDkcrovg996l5pIyGQvJn6s1NMcuwAB9AKCEPRadhXdlU5/ucvXyMmL7sx+TIEB0D8IQtBHOdX0QAbddZ/PrKQvebE/D5Q8a4/LIQHoKQQh6JGsKro/g/6awac9pGPd2WW9JRHODE4ECqDnEISg+7Kq6KFMui+Dv11Bn+/CfhTIPufKGiL/AIAQgiAEHZZZSQ9n0X0Z/J2HdJQru7gXDgEEgGYgCEHX3K+kR/83/3AhJAB4AgQh6Ajh+IeobFpYS57vwiD/AKCVEISgxXhK4ovokWz+UCat48hET2Z9f0k/zP8EgLZAEIL2qZWTU3n8kSx6NJt3NGHGujPbh0iCbZF/ANAeCELQGqX15LdsPiqbnszle1gxY9zY6NHSbp0RfwDwVBCEoOmEyS9ROXxCMQ13ZMa4sRv6G9ibiF0WAOgKBCFoIp6SpFJ6NJuPyqY51fQ5F3Z+D3aEC2skEbsyANA5CELQILVycjKXP5JNj2bzTibMWHfm+zDs/AMA1UIQgvjuV9JjOfS3HD62gIbYMePc2Y8DpR64+AMAqAWCEMQh50l8EY3K4U/l0uwqOtiJfdmL3TmEtTQUuzIA0DMIQlCrolpy/AEflU1P5fFe5kyEC7O2n6S/A8Ni+AcAImltEHIcFxsbe/nyZY7jBg4cGBoaqtKyQJdwlFz7a+ZLVhUd4sRGuDDrQg2cTMWuDACg9UG4devWNWvWjBo1SiqVjhkzZtGiRYsWLVJpZaDtMirp0XTJhVLuTB7vac6M6sKs7y/pa8dIMPgDAE3S2iCcNGnS9OnTGYYhhPTt23fOnDkIQmjqkYyczeNP5NKTufRRAx1kz47xYNb3N3DEYX8AoKlaG4Tm5uaKvyUSiaEhpjTAn4Qtn6dy6ak8/nIR7WvPRDizu4awQTZMdVWdubmx2AUCADxJmyfL1NbWLl269IMPPnhcg7y8vNjY2IULFwo3GYaZN2+evb19+2sUVV1dnYGBgdhVaKKMKnImnzmTz5wpIFaGzFAn+qY33RlGzf/qrYZ69F77cRzX0NDA4BDKdqmrq2NZFr3XPjr2sTU0NGTZFi5D07YgbGhoeOmll/z9/efPn/+4NkZGRsbGxtbW1ooiTExMJBJtPSOIRCLR3uI7XLWcxBeRMwXkdB55UEPC7UmEC1nzLHExJYQwhDT+3kHvPQ30XrsJXYcgbB8dW/Fasxq0IQjlcvm0adOkUunWrVufELA2Njb9+vX75z//2fpH1mQGBga69OOoHRp4cqmIns7jTzygt8ppfwdmuAv7yyDG37rl9Qu9124sy1JK0XvtI6x4CML20cOPbRsOn3jjjTeqqqoOHz6sb32kh+Q8SSihZ/PouXw+voj6WjJDnJgVwZIwR8ZYd34pAgAQ0vog/Omnn3bu3BkWFjZ69GhhydGjR42NMQ9Cd/CUpFbQ2EJ6KpeezOWtjZgIF2aWL7tnKGtlJHZxAAAq09ogHDly5MmTJ5WXYFyoG+5X0lO59FQuPZ3HWxoyES7MZC9mY5iBDcIPAPRDa4PQzc3Nzc1NpaWA2gjhF1NAz+RTI5YMcGAiXJiv+kldzbBPBQD0Ds41qhc4Sm6U0QsF9HwBjSngTaXMECdmhCvz7z6sC8IPAPQbglBnyXlyvYzGFNDYwj83ew5wYEa4MKv7Sj1xhSMAgL8gCHVKtZwkldDYQnoqj48vom5mTJgjM9mL+W6AgS0mNgEANAdBqPUeycjlInoqj48poNdKqa8lE+HCzO/B7huGa/sBALQMQaiV0h/RuCIaU0BjCmheDR3gwIQ7sqv6siF2jGEL5xICAID/gSDUDnUcuVJC4wppXCG9WMQbsMwAB2aAA/OOHxtghavaAgC0H4JQcxXUkoRiPraQxhTQ62V/7vAb78F8HiLtYYXoAwDoGAhCDcJRcruCCsl3pYTm1tA+dswAB2ZZb0l/B8YU7xUAgArgy1VkD6rp5WIaX0QvF9OrJdStE9PfgRnqzCwNYrt1xrAPAEDlEITqJhzhcKWEXimhsYW0tP7PYd+HAWx/BxYnNgMAUDMEocoJGzyv/BV+SaW0uyUjnNVscS/WzwqXigEAEBOCUCWyq2hiCU0opvFF9GoJdTZj+toxz9oxr3Vle1kzBjjCAQBAYyAIO0ZeDU0spoklNLGEXimhUoYJsWNCbJklvdhn7Rkc2A4AoLEQhO1U0UBulf29q6+snvawYoJtmZe92A39GS+czBMAQEsgCFvrYQO5Wfb3rr7cGupvxQTbYlcfAIB2QxA+VkkdSSqlF3OlNyu5KyW0vJ4G2zLBtswYd2Z5MIsLOAAA6AYE4d+yqmhSKU0qoUmlJKmUVslooA0TYEEmeDArQ1ifzhjzAQDoIL0Owryavzd1Xi6mMp74WTLBtswkT+bffdjulgzLkMrKWnNzE7ErBQAAVdGjIKyUkVvl9EYZvV5Kk0rprXLqaMIE2TBBtsy7fuxmG8YBeQcAoH90Ngg5StIf0Rtl9EYZvVlGbpTRwlrqZ8X0tGZ6WjMve7OBNoyFgdhVAgCA2HQnCIXjGVIqaHI5vVJCr5VSC0MSbMv0sBKu2MD6WjIS7OUDAID/pa1BKONJ2sO/Yy+lnJTWUx8LRtjJN9mT7WXDmGPABwAALdGaIMyroSnl5M/Yq6B3HlJ7Y8bPigTbMq93Zf0sGRzJBwAA7aChQdjAk7sP/8y85PL/mdIZ5sjM9mWDbRkTDa0dAAC0iUaESQNPblfQlHJ6s5ymVpAbZTS/hna3ZHpaMwHWzHAXtqc1Y2csdpUAAKCLxA/CNTf5T65wXuZMDyumhxXzijf5dx/WxwITWwAAQB3ED8L5Pdj5PVhcmQgAAEQhfhAiAgEAQERIIQAA0GsIwhZcuHChsrJS7Cq0EqX0jz/+ELsKbZWXl3ft2jWxq9BW8fHxZWVlYlehrY4fPy52CeqGIGzBypUrk5KSxK5CK5WXl8+ePVvsKrTViRMnfvjhB7Gr0FZff/11bGys2FVoJZ7nX375ZbGrUDcEIYAmopSKXYJ2QwdC6yEIAQBAryEIAQBArzEdvgHhnXfe+f3337t27dqxDyuWxMREHx8fS0tLsQvRPnK5PDY2dtCgQWIXopXy8vIqKir8/PzELkQrXbt2zdXV1dbWVuxCtA+l9MyZM8OGDRO7kA4zfvz4d95558ltOj4I09LSkpKSbGxsOvZhxZKTk+Pk5CSVin/ApTbKyMjw9PQUuwqtVF1dXVVV5eDgIHYhWikvL8/GxsbIyEjsQrSSjn1sPT09vb29n9ym44MQAABAi2AfIQAA6DUEIQAA6DUEIQAA6DUEIQAA6DXJsmXLxK5BQ2VlZR0+fLigoMDT05NlH/uLITU1NTk52cPDQ42laYH09PQjR46UlJR4eHgwTONrS3IcFxcXd+7cucLCQnd39yd0r56Ij48/efIkIcTJyanZBtnZ2YcOHcrPz3/y2qiHKKXR0dFnz541NTVt9niJsrKy06dPx8fHSyQSe3t79VeoyXJyclqzXt2+ffvmzZu6NJW0MQrNOX36tLW19cyZM/v06fPcc89xHNdss/z8fAcHB0tLSzWXp+GOHDliY2Mza9aswMDAiRMnNm3Qt2/f4ODg6dOnBwUF+fv7l5eXq79IzbF06VIPD4+3337bxcXl66+/btrg3Llz1tbWM2bMePbZZyMiIuRyufqL1Fivv/56jx49IiMj7ezsdu/e3ei/qamp5ubmo0aNmj59uq2t7aJFi0QpUjNFR0dbWVkJ69XQoUMft14VFRU5OTmZmpqquTx1QhA2r3///t9++y2ltKamxsPD4/jx4802mzhx4qJFixCEjQQEBGzdupVS+ujRIwcHh7i4uEYN7t27J/zBcVxgYGCz3/56oqCgwNjYOD09nVJ65cqVzp07V1VVNWoTHh6+bt06SmldXZ23t3dUVJQIhWqk69evd+7cuaysjFJ66NAhLy+vRr9ZHz58WFRUJPx98+ZNQkhhYaEIhWqkQYMGrV27llJaV1fXtWvXo0ePNttsypQpS5Ys0e0gxDaWZpSWlsbFxU2aNIkQYmJi8sILL0RFRTVttnfvXpZlx40bp/YCNVpmZmZycvKECRMIIebm5s8991zT3lMc38qyrIODQ0NDg7qr1BgnTpzw9/f38vIihPTu3dva2vr8+fPKDR4+fHjhwoWJEycSQoyMjEaPHt3s2qifoqKihg4damVlRQh5/vnn8/PzU1JSlBtYWFjY2dkJfzs6OjIMo88rm7LKysrz58+3uF4dPXq0qqpKaKbDEITNyMvLMzAwUOxOcHFxyc3NbdSmtLT0k08++frrr9VenabLzc21srIyMzMTbjbbewqxsbGXLl2aMmWKuqrTOLm5ua6uroqbTbsrLy+PZVnFvsMn96e+Ue494TP7hM7517/+NWLECOXe1md5eXkMwzx5vXr48OHixYu///57tVenbvp75rD333//7NmzjRaGhoZu3LiR4zjl/cYSiUQulzdqOX/+/EWLFrm4uGRlZam8Vs0TGRl56dKlRguHDRu2Zs0ajuOUZ8c023uCu3fvTpky5fvvv3d3d1dhrZqtUXdJpdJG3SU0ULR5Qn/qIY7jlE9/2LT3FH744YcjR47ExMSoqzRN15r1asGCBfPmzXN1dS0oKFB7gWql10E4c+bMRgs7depECHF0dKyvr3/48GHnzp0JIYWFhY3m8j148GD//v1mZmaXL18uLCysqamJjIz87LPPnJ2d1Va/uD766KOqqqpGC4XucnJyqqiokMlkBgYGpLneE2RkZERERCxfvlyfh4OEECcnp6KiIsXNwsLCRmuRo6Mjx3FlZWXC+Xsf15/6ycnJ6fbt28LfPM8XFRU1+xncsmXLF198cfbsWRcXF/UWqLmE9aq0tFTYdNx0vSoqKtq1a5eBgUFkZGRxcXFDQ0NkZOTHH3/cpUsXkUpWJbF3Umoinue7d+++d+9e4e+AgIBdu3ZRSuVyeUlJCaW0qqpq71/+9a9/mZmZ7d27t6KiQuS6NYNcLndzczt27Jjwt5eX12+//UYplclkpaWlQpvs7GwvLy9hOpKeu3fvnrGxsdAzmZmZxsbGxcXFlNKamppHjx4JbQICAnbu3Ekp5Xk+KCho27ZtIhasUc6dOyf8bKWUxsTE2NvbC39XVVUp5hzt2bPH2dk5JSVFzEI1Uq9evbZv304p5Xk+ODhYmOCm+JarqalRfMv9+9//NjIy2rt3r65O8EYQNm/79u0ODg7/+c9/pkyZ4ufnJ3y6EhISCCEymUy5ZWxsLGaNNrJx40ZXV9c1a9a8+OKLISEhwkS+06dPGxsbCw2CgoJcXV1n/0X4NOqtadOmhYaGrl27NjAwcMGCBcLCZcuWjRgxQvh7165d9vb2q1evnjp1qq+vb21trXjFapyBAweOGjXqq6++8vb2Xr16tbBw5syZM2bMoJTevHlTIpEMHjxYsbLduXNH1Ho1yJ49e+zs7FavXj1t2rRu3boJ69W1a9cIIdXV1cotExISdHvWKA6ob17Pnj379Olz69atbt26rV+/Xpj6YWho2LVr15CQEOWdOoaGht7e3oGBgeIVq3FCQkICAgJSU1MDAgLWrl0rXA3HyMjI19c3KCiIENKpU6fQ0FDnv/j4+OjzbsJx48YZGhpmZGRMnjz5ww8/FNYuCwuLnj17+vj4EEL8/f379euXnJzs4+Ozfv16c3NzsUvWIC+//HJtbW1+fv7s2bOnT58uLLSysgoMDHR3d+c4rmvXrgEBAYqVLSAgAB0o6NGjR2ho6K1bt7y9vTds2CB0i6GhoY+PT58+fZTnSRgYGHh7ewsfXp2EyzABAIBew+ETAACg1xCEAACg1xCEAACg1xCEAACg1xCEAACg1xCEAACg1xCEAACg1xCEAACg1xCEAACg1xCEABqqsrKyoKAA534CUDUEIYDG2bdvn7+/v4WFhZOTk5GR0fLly8WuCECX6e/1CAE00+bNm2fNmjV58uT169c7Ojo+ePDA0NBQ7KIAdBlOug2gQWpqalxdXcPCwo4cOSJ2LQD6AptGATRIQkJCeXn5jBkzxC4EQI8gCAE0SHFxMSFk/PjxjJK5c+eKXReALsM+QgANYmlpSQjZtGnTsGHDFAstLCzEqwhA92EfIYAGefTokZOT06RJk7Zs2SJ2LQD6QrJs2TKxawCAPxkZGRFCvv766+rqah8fH2Nj48zMzNzcXAcHB7FLA9BZGBECaBZK6Zo1a1auXFleXi4smT9//rp168StCkCHIQgBNJFcLk9NTa2trXV2dnZ2dmZZzGsDUBUEIQAA6DX8zAQAAL2GIAQAAL2GIAQAAL2GIAQAAL2GIAQAAL2GIAQAAL32/05aUPq1M7IVAAAAAElFTkSuQmCC", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "M(ϵ) = [-ϵ 1; 1 -1]\n", "plot( ϵ->cond(M(ϵ)), -1/2, 1/2, title=\"Condition number of M\", xlabel=\"ϵ\")" ] }, { "cell_type": "markdown", "id": "d2263eed", "metadata": {}, "source": [ "2. If we did not use pivoting, the LU decompositon of $M(\\epsilon)$ is \n", "\n", "\\begin{align}\n", "\\begin{pmatrix}\n", "1 & 0 \\\\ -\\epsilon^{-1} & 1\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "-\\epsilon & 1 \\\\ 0 & \\epsilon^{-1}-1\n", "\\end{pmatrix}\n", "\\end{align}\n", "\n", "Both of these matrices are illconditioned. In the following, we can see that the condition numbers explode near $\\epsilon \\approx 0$. " ] }, { "cell_type": "code", "execution_count": 19, "id": "4a8de60c", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd5hU1f0/8M+5d3rZmdnZ3thCB+lFEGmCiqKIKFhD7ElMjDEIWKIm5mvBqFF/ahRLlBgTFSwREEURBOm97y67sCzb+/R2z++PC+O6tGWZ2dm7+349PD4zd87c+exl2Lfn3HPuZZxzAgAA6KqEWBcAAAAQSwhC6IgkSSosLFy/fv2WLVuqqqpiUsO99947f/788NO9e/fecccdH3zwwVnfuHbt2jvuuGPp0qXRrK7jWrBgwT333NPQ0BDrQgBaC0EIHUtZWdm9996bnJzco0eP0aNHDx8+PDk5uX///i+++KLP52vPSt56661///vf4aelpaXvvPPOunXrwls2btz45ptvFhcXt3hjfn7+O++8s2PHjnYqtIP54osv3nzzTZfLFetCAFpLFesCAH6ybt26adOm1dbWpqWl3XXXXTk5OaFQKD8/f/ny5Q888MCePXvefvvtWNVmt9vHjx/fq1ev8JbFixc/99xzH3/8cU5OTvOWqamp48ePz87Obu8SAaBNEITQURQXF1911VX19fX333//008/rdPpwi+53e4XXnihsLAwhuUNGzZs1apVrWk5ZcqUKVOmRLseAIgUBCF0FHPnzq2vr7/hhhtefPHFFi8ZDIZHH320tra2+cZQKLR27drdu3cHg8Hc3NxJkyYZDIbmDYqKiurr63v37m0wGNauXbt9+3ZRFC+++OIBAwac/Om1tbXLly+vrKzs1q3blClTjEZjiwZNTU0FBQVJSUmZmZlEtGvXroqKCvlTtm7dKrfp27evXq+vra09fPhwenp6SkpK8z00NDSsXLmypKREp9MNHz58+PDhLT5i69atWq22f//+Lpdr6dKlJSUlKSkpl19+eUJCwpkPncPhyM/PT0hI6Nat27Fjx1asWFFfX9+jR4/LL79co9GEmzmdzoMHD9rt9ha91SNHjtTU1PTs2dNsNstbtm3bplKpBgwY4HA4li1bVlZWlpWVNXXqVK1WS0Scc/l4ajSayy67rEWHOMzv969YsaKwsNBisVx22WXp6eknt5EkacOGDdu3b/f5fDk5OZMnTzaZTM0b7N271+v1Dh48OBgMfvPNN4WFhampqTNnzpTfu3HjxoKCgqqqKrvd3q1bt9GjRzf//yeA1uIAHUB5ebkoioyxwsLC1rTft29f//79m3+Tk5OT//e//zVvM2vWLCL66quvJk6cGG7GGPvDH/7QYm+LFy+2WCzhNikpKT/++KNGo8nMzAy3+eqrr4jo17/+tfw0OTn55H9Ne/bs4ZzL47d//etfm3/E22+/HRcX17zxxRdffOzYseZtiCgnJ2fNmjXNdx4XF/fNN9+c+WisXLmSiO68886XXnpJrVaH39unT5/S0tJwszVr1hDRL3/5yxZvv+uuu4ho5cqV4S1arTY1NfX7779vnsE9evQoKSmprq4eO3ZseKNGo/nggw+a7+2iiy4iom+//bZ79+7hZjqd7rXXXmvxudu3b2/xl5iQkPDFF180b9O7d28i2rRpU25urtxm2LBhnPOSkpLBgwe3OP5ms/nMBwrglBCE0CF89NFHRNS3b9/WNK6qqkpNTSWiu+++e/v27fv373/66ac1Go1Kpfrxxx/DzeQgzMnJGTZs2CeffLJt27aFCxfabDYiWrZsWbjZ1q1b1Wq1Xq9/+eWXDx8+vGPHjtmzZ6ekpIiieIYgXL169fXXX09Ejz322DcnOJ1Ofqog/PjjjxljFovlH//4R35+/oYNG6655hoi6tevn9vtDjcjIqvVGh8f/6tf/erbb7/dsGHDb3/7WznjPR7PGQ6IHITZ2dlms3nBggUbN2785ptvxo8fT0TXXnttuNk5BaHJZLLb7b/97W9Xr169Zs2aq6++moiuvvrqyZMnjxw58rPPPtuyZcvTTz+tUqlMJlN1dXX4vXIQpqWlTZ06dfPmzSUlJe+++67NZmOMffnll+Fm+fn5FotFo9HMmzdv/fr1e/bsef31161Wq1qt3rx5c7iZHIRZWVlXX331Rx999OOPP3766aec86lTp8p/HVu3bi0pKdm2bduiRYumTp16hqMEcDoIQugQnn76aSKaMWNGaxo/+OCDRDRz5szmG1944QUiGjNmTHiLHIR9+vTx+XzhjQsXLiSiO+64I7zliiuuIKJXX301vEWSpEmTJhHRGYIwXMbHH3/corwWQRgMBuWhyCVLloTbhEKhUaNGEdErr7wS3ih3ax5++OHme7vkkkuI6Ouvvz7DAZGDsEWYNTY2WiwWlUoV/vHPKQiJaN68eeEtbrfbbrcTUa9evZofz1/+8pdEtGjRovAWOQj79evn9/vDGz/99FMiuuCCC8JbrrzySiJ69913m1fy9ddfE9Hll18e3iIH4SWXXCJJUvOWRqMxKyvrDMcEoPWwfAI6hMbGRiIKn6M6M/m3avNFfkR0zz332Gy2tWvXVlZWNt9+3333NT9PNnnyZCIKr3lwuVxff/21zWa7/fbbw20YY3/84x/b+JOcZPv27YcPH+7du7fcC5QJgjB37lwiWrJkSYv2cr62KLioqOisHzR06FA5NWVxcXEjR44MBoNHjx5tW+UPPPBA+LFer5dPat57773Nj+e4ceOo2fEMu++++5oP0k6bNq1Hjx67d+8uKCigE2dkMzIy5BwNmzx5co8ePVatWuX1eptvf/DBBxljzbdYrdba2tr8/Py2/WgAzWGyDHQI8uQUt9t91pZ+v//QoUPyVI7m2w0GQ58+fX788cd9+/Y1P8fWfMEDEckvyfNciCg/Pz8YDPbs2bPFJItBgwa19Udpad++fUQ0ePDgFr/Khw4dSkR79+5tUZ7Vaj254BbpfkotftLm783LyzvXsi0WS1JSUvMtiYmJRNSzZ8+TN55c3sCBA5s/ZYwNHDiwoKBg7969PXr02LFjhyRJarW6xf/NEFEgEPD5fOXl5c3n4PTr169Fs9tvv/3JJ5+84IILJk2aNGnSpEsuueSUc6AAWgNBCB2CPBWzNf0e+Tyc3W4XRbHFS/LvfYfD0Xxji6mkgiAQET8xDul0OomoxW98IkpMTGyRW212uo9oTbXhgiVJOusHnc97W783vV7ffKN8lPhJlyyWA7I5+QjIP6983ZmysrI333zz5I+22WwtLp5w8rzZJ554Ij09/fXXX1+2bNmyZcuIqEePHn//+9/lgW6Ac4KhUegQxowZQ0Q7d+6srq4+c0uTycQYq62tDYVCLV6S+3ktJmeedW90qg5NVVXVyb/c20Ye7z35QnFytc1nq0abnGQnHzc5qiPr5J9XPsjy3458TMaPH193GvKpwTMQBOGee+7ZsWNHaWnp+++/P3369EOHDk2bNm3btm0R/1mg00MQQoeQl5c3duzYQCDw/PPPn65NMBgkIo1G07Nnz2Aw2OIaZi6Xa9++fYyxFjPyz6xXr15qtfrgwYMtRmVb8/tUPgd2cq60INezdevWFsm6efNmOtWgX/TI6xpPjqj9+/dH/LO2b9/e/KkkSfIW+WgMGjSIMbZ58+ZAIHCeH5Senn7rrbcuWbLk4YcfDgaDn3/++XnuELogBCF0FM8++6xarf7b3/72/vvvn/zq//73v/AElhkzZhDRM88807zBq6++2tjYOH78+LMuP2/OYDBMmTKlsbHxrbfeCm/knJ8hj8PS0tKI6KxTUQYOHNi9e/f8/Pzm82JCodCCBQuISF6D0T4yMzPVavX69eubXxF7+fLl0bgs6ssvv+z3+8NPlyxZUlRUNHjwYPlsZVJS0pVXXllXVycfhBbO2kMNBoMtZtMQkbyipp0vSAudA84RQkdx4YUXvvHGG3fffffs2bPfe++9GTNm5ObmBoPB/Pz8xYsX//jjjzfffLPccs6cOe+///4nn3zyy1/+8u677zYajUuWLHnqqac0Gs0pf7Ge2V//+tcVK1bMnTvX6/VOnTq1qanpxRdfzM/PP/kcZAvyLMrnnnuuqakpPT2dMXb99dfL6xSbEwThhRdemDZt2m233VZaWjpp0qTKysq//e1vmzZtGjRo0G233XauBbeZRqOZPn36Rx99dOWVV86fP99kMq1Zs2bBggXdu3eP+OXr3G731KlTH3nkkZSUlO+++27+/PmCIDz77LPhBq+88srGjRsfffTR/fv3z5gxIy8vr76+vrCwcMmSJcFgcMWKFWfYeWNjY48ePW699dYJEybk5uYKgrBp06bHH39cEAT5f5IAzk0s124AnGTt2rXyQrTmDAbDnXfe2fwiKQUFBS0uUZaZmdniCizyOsINGzY03yj3GPr06dN849KlS+VFcrKMjIwtW7ac+coysmeeeab5G89wZZl///vfzVsS0WWXXVZZWdm8DRHl5OS0OCDvvfceET3++ONnOGjhK8u02D579mwi+uGHH8JbKioqml+QRavVvvHGG6e7sswp97Z69eoWh46IfvOb34S3yH99a9as6dOnT/O/wX/+858tdlhYWNj8oj8ym83W/IeVTxa2uJ5AY2Pjydc0t9vt//73v89wlABOh3HcoR46ntLS0k2bNlVVVanV6uzs7JEjR7a4BCURcc63bt26e/fuQCDQo0ePiy66qPn6NiKqqqpyOp3p6eny8vDwu4qLi9VqtTxPNaypqWnlypXV1dWZmZkTJ07U6XTFxcWiKGZlZckN3G73sWPHrFbryfMhGxsb5eugZmRkaDSapqamysrKhISEFr1Dt9u9Zs2aw4cP6/X6YcOGnXx2sKioSKVShT9R5nQ6q6qqbDbbyX3NMI/HU15ebjabW9RWXV3tcDjS0tKaLw4JhUKrVq06dOhQXFzcpEmTEhMT5WapqanhGaHFxcWCIHTr1u2se3O73RUVFXFxceER6bKyMq/XK/8Uq1evDn/Q6YasCwsLN2/e3NjYGB8fn5mZOWzYsOYLEI8ePRoIBHJyck6exHv48OHdu3dXVlZqtdrs7Ozhw4fjQqPQNghCAADo0jBZBgAAujQEIQAAdGkIQgAA6NIQhAAA0KUhCAEAoEtDEAIAQJeGIAQAgC4NQQgAAF0aghAAALo0BCEAAHRpsQnCjz76aN++fa1v37ZbbHdxuHhe2+DL1gY4aG2D49YG0ThosQnCL774osV9O8/M5XJFr5jOyuPx4J9ZG+DL1gY4aG3AOW9xO2hojWh82TA0CgAAXRqCEAAAujQEIQAAdGkIQgAA6NIQhAAA0KUhCAEAoEtDEAIAQJemgCCUJL5x/Y5YVwEAAJ2TAoKwtr4p+ev/F+sqAACgc1LFugAAAGijo0ePXnfddaFQKNaFtIdnnnlm0qRJ0dizAoKQMUaEy2YCALRUUVHhdDrff//9WBcSdc8880xBQUHXDUIAADgdo9E4dOjQWFcRdYmJidHbuQLOERIjFusSAACgs1JCEAIAAESNAoKQEWO4tR4AAESHAoIQAAAgehCEAADQpSkgCBkxxjA0CgAAUaGAIAQAgE5g9uzZS5YsIaJt27b97ne/a/HqTTfd9L///S8WdSkhCJlA6BACACja3r17t2zZMn36dCIqKyv78ssvWzR46KGHHn744ViUpoQgBAAABfnTn/5UV1dHRHv37n3yySe9Xi8RLVy48IYbbmDstMvCL7jgAlEU165d236FnqCAIGQMyycAABTj+eefb2ho2L59+9SpU8eMGaPT6Yjoq6++Gjt27JnfOHbs2K+++qpdavwZXGINAKDzcAdp5OdBd7CdPk5g9PIocUpmy37exo0bH3300Y8++mj48OFEFAqFCgoKcnNzz7y33NzcH374IVq1np4CgpAxzBoFAGgVg4qWXS4GpPb7xHTDKUY7582b5/F44uPj5ac+n0+SJK1We+Zd6fV6t9sd+RLPRgFBCAAArZdpjP3lmRcvXrxmzZobb7xx7dq1Go3GYDCYzeba2tqkpKQzvKu6ujolJaXdigxTwDlCAABQFrvd/sADDyQmJj766KPyllGjRu3YsaN5G96MvGXHjh2jRo1q71oVEYQCY4TJMgAAisIYe+eddxYtWiSvlLjhhhvCywQFQSgtLdU0U1tb6/V6f/jhh2uuuab9S8XQKAAARFL4PF9ycnJ5ebn8+IYbbnjmmWcqKipSUlKuuOKKQCDQ4l1vv/329OnTzzx2GiUK6BEywv0IAQCUTa/Xv/baawUFBadr4Pf7//znP7dnSWGR7xFu3769sLCwf//+ffr0ifjOAQBAoS655JIzvPrrX/+63SppIcI9wueee+7BBx+sqqqaPXv2Bx98EJmdCowI5wgBAJTN7/efdXVEU1OTJLXj4g8iingQLl26dMGCBffee+9jjz22bNmyiOyTEWMIQgAAhZs9e/a3335LRMuXL//LX/7S/KVZs2YdOXKEiJ544om33nqrnQuLcBDecccdzz///H//+9+FCxfOnj07sjsHAACF2rVr1549e6ZOnUpEBQUFq1evbv7qJ5980tDQQERz5sz561//evJUmqiKcBAyxpxOZ2Fhod/vj1T3ljH0CAEAFCN80e3t27c/9dRTfr+fiBYuXHjjjTee4aLbsrS0tF69erXz/ZgiEISVlZXHjh2TH8+dO/ef//znI4888re//e2xxx47/50TEWOYNQoAoBjyRbe3bNly/fXXT5w4UaPRENGKFSsuvvji1rz94osvXrFiRZRr/JnWzhq96aab1q1bV1JSsmzZsilTpsgbA4HAzTffvGbNGo1G07Nnz88++ywxMXHfvn0XXXTR3r17ExMTI1Ii7j4BANBK3O+rfP63FGyvq24TWa/9la7fyBYb16xZ8+STT3766acDBgwgomAwWFhYmJOT05od5uTkyKcS201rg/CKK66YN2/e1Vdf3Xzjv/71r4MHDxYXF2s0mssuu+yll156++23582bFwgEzGbziy++GJESmYChUQCAVmEabeJvnuFBfzt9HBNFa8LJ2//85z87HA65L0hEfr+fcx5+KghCKBQKN5YkiXMuCMdHKHU6ncfjiXLhP9PaILzllluIKFyo7MMPP7ztttv0ej0R3XPPPX/+858feeSR1iT5/v37165d+/zzz8tPtVrtsmXL1Gr1KRv7A0FG5HQ6W1kqyNxudzAYbPFXBmflcrnOehoDWsBBawPO+fnfacHtdvOTBsxEi/08d3v+lixZsm7duhtuuGHDhg06nc5gMFgslurqavnCMenp6eErzhBReXk55zwjI0N+WlVVlZaWdvI+fT6f0+k81y+bTqdTqc6SdOe1oP7w4cPdu3eXH/fo0aO4uLiVb8zNzZ0+fXp4iFWn09lsttM19geCTUQmk+l8Su2CBEHQ6XQIwnPFOceX7VzhoLUB55wxdp7HzWAwdMz/BbFYLPfee+/XX389d+7cl19+mYhGjx69ffv2fv36EdFFF11UUVGxZMmSa6+9VpKk5557btSoUeEU2LZt2+jRo0/ep1arNZlM0fiynddvSafTKXcHichgMHg8nmDrBqa1Wm1OTs7QE+RDc9oScdFtAADlsFqtgiAwxt56663ly5evXLmSiG688cYvvvhCbpCUlPTxxx8/9NBDqampSUlJO3fufP/99+WXgsHgihUrZs6c2Z4Fn1ePMCkpSV75QUR1dXUJCQln7YG2AW7MCwCgIGVlZfKDpKSk8MVFZ86c+dRTT5WWlspDoJdeeunBgwddLpdarQ6fOySiTz/9dNy4cdnZ2e1Z8Hn1CAcNGrRhwwb58caNGwcNGhSJklrCrFEAAKXTarWvv/56YWFh841Go7F5ChKRw+F4+umn27e0VvcIv/3227q6OpfLtXbtWqfTOWnSJJvN9pvf/ObSSy8dM2aMyWR65pln3njjjWiUyBh1yDFwAAA4B+PHjz9rm9tvvz36hbTU2iDcvHlzcXHx9OnTa2pqVq5cOXLkSJvNNmLEiPfff/+1114LBAILFixosbgiUtAjBACA6GltEM6fP/+U26+++uoo5R8AAEA7UMbces7YyWtlAAAAzp9CgpBjAQUAAESFQoIQPUIAAIiOyC/7iwbOEYQAAC3p9frt27fHx8fHupCoc7vdUVqYQIoJQsZCnCujVgCA9tK/f/+amppI3fy1g7NYLFHaszLChRPjEnqEAAAtRS8eug6cIwQAgC5NIUFImDUKAABRoYwglEiQeJcYBAcAgHamjCAkYoRzhAAAEAXKCEJOhBgEAIBoUEoQCqGuMT8YAADamUKCELNGAQAgOpQShIS5MgAAEA3KCEKJBImQhAAAEHnKCEJcWQYAAKJEKUFIOEcIAADRoJAgZCyEHiEAAESBUoJQ4FIo1lUAAEAnpIwglIhJ6BECAEAUKCMIOTGOBfUAABAFCglCJoQwWQYAAKJAIUFIjGNFPQAARIEyglBiAs4RAgBANCgkCLGgHgAAokMZQcgFQcJkGQAAiAKFBCFuwwQAANGhlCDE0CgAAESFMoJQEgSsIwQAgGhQRhByYhLWEQIAQBQoJQhFCdcaBQCAKFBIEDJcaxQAAKJCGUEoMVxrFAAAokIZQcgJk2UAACAqFBKEohAKIggBACDyFBKEJOCi2wAAEA0KCUImhkKYNQoAAJGnjCCUcK1RAACIDoUEIRM5eoQAABAFyghChnWEAAAQHcoIQolhaBQAAKJCGUHImchxiTUAAIgChQShIEg4RwgAAFGglCDEZBkAAIgKZQQhCSJ6hAAAEA3KCEKcIwQAgChRRhCSgPsRAgBAVCgjCLkgUBBBCAAAkaeMIESPEAAAokQxQUiYLAMAAFGgmCDEZBkAAIgG5QRhKBjrIgAAoBNSSBCqMDQKAABRoYwgZIKaQoFYVwEAAJ2QQoJQxCXWAAAgKpQRhCSKJOEcIQAARJ4yglAQVQyTZQAAIAoUEoQqFXqEAAAQDapYF9AqgqjC8gkAAIiGqPQIDx06dODAAc55pHYoqERCEAIAQBREuEdYU1MzY8aM5ORkURSvvfba66+/PiK7FdWqEJZPAABAFEQ4CJ944olf/OIXd9xxR2R3K6o0ooQgBACAyIvw0OiyZcsYYzNnznz88cfdbnekdiuq1GLIH6m9AQAAhEU4CCsrK3fv3v3qq6+63e6HH344UrtVqVUCeoQAABAFEQ7ChISEu+66KzEx8Ve/+tXq1asjtVu1WiXiHCEAAERBa4Nw7ty5EyZMyMvL27hxY3gj5/zBBx+0Wq0Wi+W+++6TJGny5Ml79+4loj179mRnZ0eqSpVarUaPEAAAoqC1k2V0Ot1999131113eTye8Mb//Oc/n3/+eX5+vkqlGjdu3LvvvvuXv/zlF7/4xT/+8Q+Px/Puu+9GqkqNVh1EjxAAAKKAndNqv9TU1A8//HD8+PHy00svvXTKlCl/+MMfiOjNN99ctGjRDz/8QETBYFClOlPEXnTRRYFAIDMzU35qNptfeukltVp9uvZFh0uD/3yi5xNvtb5UcLvdOp1OEJRx8aCOw+l0mkymWFehMDhobcA5d7lcOG7n6ly/bBqN5sx5ROe5fCI/P3/OnDny4379+hUUFBzf6dk+NT4+PicnZ8yYMfJTnU5nNpvP0N5kMrlCAZ1Odz7VdjWSJCEI2yAQwDftnOGgtQHnPBQK4bidq3P9sjHGztrmvIKwoaEhnMxms7murq6Vb7RYLCNHjpw5c2Yr2xv02iD34Xf6ORFOiHUhCoOD1gY4aG3AOcdxa4NoHLTz2p3dbm9sbJQfNzQ0JCYmRqKkU9Bp1RopEMFrtgEAAMjOKwj79Omza9cu+fGuXbt69+4diZJOgTEWENReH+bLAABAhLU2CPft27d169ZAIJCfn79161av10tEd91116uvvnrw4MGioqK///3vd955Z/QK9TONx4OLywAAQIS19hzhk08+WVBQkJ2d/eabbxLRxx9/nJOTM23atIMHD15xxRWSJN1xxx033HBD9Ar1iVqvz0uEGVYAABBJrQ3CDz/88JTb586dO3fu3MjVc1o+QeP2+NrhgwAAoEtRzIQlv6j1ehGEAAAQYYoJwoBK5/N4Y10FAAB0NooJwqBK6/ciCAEAIMIUE4Qhlc6HIAQAgEhTUhAGfDhHCAAAEaaYIJQ0+pAnYre8BwAAkCkmCEmrD/owNAoAABGmnCDU6CWf5+zNAAAAzoViglDQ6jmCEAAAIk0xQajS6cmHc4QAABBhyglCvUFEEAIAQKQpJgi1JqPod8W6CgAA6GwUE4Q6g1HjR48QAAAirLV3n4g5o8noD6JHCAAAEaaYHqHZZNAF0CMEAIAIU0yPMM5idgedsa4CAAA6G8X0CONMeq3kDwZDsS4EAAA6FcUEIWPMJRoamnCaEAAAIkkxQUhETrWpoaEp1lUAAECnoqQg9KjNTU2OWFcBAACdimImyxCRVxcnNiIIAQAgkpTUIwzqzJ4mDI0CAEAkKSkIJb0l4GyMdRUAANCpKCkImTEu6EKPEAAAIklJQag2W5gLPUIAAIgkJQWh3mpVu+tjXQUAAHQqSpo1arZaQ96GWFcBAACdipKC0Ga3Bb0YGgUAgEhSUhAmJ9qDAQyNAgBAJCnpHKHBoJWY2NCEmzEBAEDEKCkIiaheba2qqo11FQAA0HkoLAgdentdDYIQAAAiRknnCInIZ0zgCEIAAIgchfUIpbgEb311rKsAAIDOQ2FBqLIl8MaaWFcBAACdh8KC0GRPVDehRwgAABGjsHOE9uQkclfFugoAAOg8FBaEaenJggc9QgAAiBiFDY3arGbGeW2DM9aFAABAJ6GwICSiakNSaWl5rKsAAIBOQnlB6DCn1JQhCAEAIDKUF4QhW5qrCkEIAACRobwg1CSm8Zpjsa4CAAA6CeUFoS01zdBQFusqAACgk1DY8gkiysrJqnKhRwgAAJGhvB5hSrJdI/lr6rGCAgAAIkB5QUhEFcaM4kNHYl0FAAB0BooMwiZ7Tu3hQ7GuAgAAOgNFBiFLy5OOFca6CgAA6AwUGYT2nO7mGvQIAQAgApQ3a5SIeq5zOqgAACAASURBVPXKbXAdlUJBQVRk/QAA0HEos0do1pZpk0swXwYAAM6bIoOQiGrj844WYnQUAADOl1KDUErJc5QgCAEA4HwpNQgTcrvrKhGEAABwvpQahH379khvKpIkHutCAABA2ZQahAk2k0Ntyi/GRUcBAOC8KDUIiagqsW/x7j2xrgIAAJRNwUEo5g4IFu2KdRUAAKBsCg7CvIED0isQhAAAcF6iEoR33nnnsGHDorHn5rrnZkjEjpTgJr0AANB2kQ/CpUuXqtXqpqamiO/5ZGXJ/Q9sR6cQAADaLsJB2NjY+OKLL86fPz+yuz0dbfeBvsKd7fNZAADQKUU4COfMmfPEE0/odLrI7vZ0+g0ZkFm5E4sJAQCgzSIZhDt37tywYcO+ffsWLVrU2Nj43//+N4I7P6WMzDRJUO0+WBLtDwIAgM6qVUFYVVX13HPP3XDDDZdffnnz7R6P5/e//33v3r3Hjh27atWqhISE3/3ud9Gp87QqMkcUbVrfzh8KAACdRqvu51deXp6fn5+Xl7dgwYLm2//0pz/t3r176dKlmzZtuuaaawoKCu6++24iqqysfPPNN2fNmhWVkn8uYcgo5/JFRO3xWQAA0Pkwzlt7hm337t1DhgwJBALy00AgkJycvGzZsgsvvJCIrrjiiokTJ86ZM4eIJElyOBwWi+V0u5o4caJer+/bt6/81Gw2z5kzRxTF07V3OBxms/mULwUCwSOP/9L4h1dSE22t/EG6CLfbrdPpBEHBS0Vj4gxfNjgdHLQ24Jy7XC6TyRTrQhTmXL9sarX6rL8G2/5bsqysrL6+fujQofLToUOH7tlz/IJngiCcIQWJSBRFo9FoPcFgMDDG2laGWq0qTh2y68dNbXs7AAB0ca0aGj2l6upqg8GgVqvlp1ardevWra18b3Jy8pQpU26++eZWtvf7/Vqt9nSvWgaMbtr0jVZ7dSv31kWEQiGtVose4bk685cNTgkHrQ0458FgEMftXEXjy9b235JWq9Xr9YZCIfmpw+GIj4+PUFXnZtiYEd3r95XXu2Py6QAAoGhtD8K0tDS1Wp2fny8/PXjwYG5uboSqOjd6g/5I8sCNq36IyacDAICitSoIOef19fXyVdPq6+sbGxuJyGAwzJgx4/nnn+ecHzhwYNmyZa0f6oy4uBGT9bu+idWnAwCAcrUqCGtra/Py8q666iqz2ZyXlzd69Gh5+3PPPXfw4MGEhITRo0c/9dRTvXr1imapZzJ8zIhUZ+n+IlyAGwAAzk2rJsskJCTU1dWdvD0tLe2HH35wuVw6ne4Mix/agahSlXUf61j1bZ/cW2NYBgAAKE4EphQajcbYpqCs14RJmQXf+oK48CgAAJyDzjO3PqdXz5DW9O33WFAIAADnoPMEIRGxMTPY2k9iXQUAAChJpwrCUZeMs3urt2zbH+tCAABAMTpVEDJBqB1yTcXKxbEuBAAAFKNTBSERjblySlbVnqLDx2JdCAAAKENnC0KzUVvS5/KdX3wa60IAAEAZOlsQEtGYa6/pWbLmcGVjrAsBAAAF6IRBaLVZj+VctO6LL2JdCAAAKEAnDEIiGjJtxgX5ywqrnbEuBAAAOrrOGYQJGRl1OSPW//fDWBcCAAAdXecMQiIaefNtg4+s3Li/NNaFAABAh9Zpg1BvsdYMv678k4W49igAAJxBpw1CIrr42ulp7qP/W7k51oUAAEDH1ZmDUFSp4qbdZVv5ZrUrGOtaAACgg+rMQUhEvS8cxayJn378ZawLAQCADqqTByER9bvlnlF7/rN6H25eDwAAp9D5g9CW0S04dmbgg6fr3YFY1wIAAB1O5w9CIhp81XRmTfj2nbdiXQgAAHQ4XSIIibER9/wxu3TDym/WxboUAADoWLpGEBKZ40yGWx5OWvFKcUl5rGsBAIAOpKsEIRH17d+rdOh1R956yuvHagoAADiuCwUhEU25YYbXkPD5m2/HuhAAAOgoulYQMsYu+vUfu5f++MVXOFkIAABEXS0IichsMVl/8VD2yld+OFAR61oAACD2ulwQElFe39508QzXoqf21mBlIQBAV9cVg5CIBlx9XUqCZeXbb5a6cHcKAIAurYsGITE24J55l/j3ffD6u5WeWBcDAACx01WDkEgwmPr+8dnL3Fvf+8f79b5YVwMAADHSdYOQiASDacADT13WuO6df/zLidOFAABdUpcOQiISTJb+c56dWLtm4RsfOJCFAABdT1cPQiISTdYLHnx2YtXqV1//sMYb62oAAKB9IQiJiFRmW/85z06t/e7//eO/xzCPFACgK0EQHifG2fr84elZdV+99NbnxQ5kIQBAV4Eg/IloTeh5/7N31Xz66sLFG6qQhQAAXQKC8GfE+KTcPzz3K8+3W9588YN8f6zLAQCAqEMQtiTGJ+XNeeGKBKdt0fwFP9aiYwgA0LkhCE+BafW59/xp0Mghl35+//2L85uwrAIAoPNCEJ4GY2lTb8mbdee9mx77/Ts/7KpDzxAAoHNCEJ6Jeci4nN/936OlCz95e9F7+aFYlwMAAJGHIDwLdUZe9pwX7+BbxI8X3PaNoxETaAAAOhcE4dmJFnvWH/52eXfjg6t/d+t7u9ZUYJgUAKDzQBC2ClNrEmbdl3vL71468uzmd197eKPPi4FSAIBOAUF4DnS9h3ab/9qN8dXTv7zv2g8OrceiewAA5UMQnhvBZEm754l+V1//2r6HPvnXp39cH3QFY10TAACcBwRhWxiGT8p44IU/St9fufqJCR+ULS6WYl0RAAC0EYKwjVRJGal/eGHY4L4fH7j/wOf/uexLz/4GjJQCACgPgrDtmKiKu+ymjHmv/Srh2Eub7n3kgw2PbQ1hpBQAQFkQhOdLtCbYb3kw55b7Xm549+Jv/jzhX8feOiiF0DkEAFAIBGFkaHsOSp//2qhRQz7J/6Pn638N/dj9MU4cAgAoAYIwYpioMo27Jn3eqzfHVy7ed+/SlRsnLw9uxBILAICOTRXrAjob0WKPv/lBw4GtTy9+rdy5/MHyXxgysv88RBySwGJdGgAAnAKCMCp0vYemzH/DvPbLRd8+dtTZ9/elNydmZT0xRBgQjzgEAOhYEITRIo+UGkdNiVu//L/fPVri6POrklvtGekPDRRHJyMOAQA6CgRhdDGN1jTuGuOoy80/fPHJdw9WuAfMqZztt6bOGyhMzRKQhwAAMYfJMu2BaXTmS2amPvZejz49/3Vwzt+PvfDyD2VDPw3+M1/y4eLdAAAxhSBsP0yrN18yM+Xht3pkJL6374EPGl/btu9I9n8C8zeHjrowuRQAIDYQhO1NMJgtU29LeWhheqr9oW2PbKz9U87RDSM+C8z8NrSuEnEIANDecI4wNgSTJe6ym82TZnl3r7/2xy+vqlm4k035XeWlGpP5t32F63IEnRjrEgEAugYEYSwxUaUfdLF+0MX+owUXrl/+5fa7anJGvOKfcf/6rBvzhDt6CYPsmE8DABBdER4a/eyzz0aNGjVkyJDLL7+8qKgosjvvxDSZPWwz70t+eGFOXvZf9jy+o2b+RdXrZn4T6PdJ8NmdUq0v1vUBAHReEQ7CzMzM5cuXb9u2bdasWQ888EBkd97piWab+ZKZKX/6p33CNRMPf7am4J5F4uKDFU09PgrMXh1aXc4lnEMEAIi0CA+NDh06VH6Qm5vr8Xgiu/Muovl4qWbNZ3/dfOf/XTB2uXDV79dn1vnoxjx2c3dcoQYAIGKico7Q4/HMmzfvqaeeisbOuw5NZo/4mx8MNda61n455euHr07Lqe878SN+4bUrdVqBrs9lt3YX8uKQiAAA5yXyQej3+6+77rrbb7994sSJEd95FyRa7HFXzjZfdpN3zwa25dvbDv3j131HFOVNeNc96MIvpD5WNjNXuDabpRmQiAAAbdHaICwqKvr+++937tw5cODA22+/Pby9srLy//7v/woKCkaMGDF//nyVSjVz5swJEybcfffd0Sm4i2IqtTxeKrmd3r0bum1e8six5/4yYMxO+4T3qvs+vjXUx8pm5AgzslmWCYkIAHAOWhuEixYt2rlzZ2VlZV1dXTgIOedXXHHFkCFDHnzwwWefffbXv/51enr65s2bk5OT77nnnqSkpCeffDJqlXdRgsFkGD7JMHxSqL7Kve37vt/9/VlB+PugcdsyJn5Yn/zUjlCumc3IEa7NZmlYGgMA0AqM83OYifinP/3p8OHDixYtkp9+//33s2bNKisrE0WxrKwsNzd33bp14cZarbZ///6n3M+VV16ZnJw8YsQI+anZbJ41axZjp+3KOBwOs9nc+jq7lGBliWfLd57NKwWLXTtk4rb08R9Vx31RwuNU/KosNjVLGJVEInqJrYYvWxvgoLUB59zlcplMplgXojDn+mUTBOEM4SI7r17Dli1bLrzwQlEUiSgtLS09Pb2urm7y5MlnfWNTU5PL5QoXp9Vqr7rqKpXqtMX4/X6fD4vpTsOarJl0o3rC9cGC7f4dP/T6+oO/5F3wzMAx+yx9v22Kf2gzP9hI41NoSrp0ZQa3qLEC4yzwZWsDHLQ24Jz7fD61Wh3rQhTmXL9sGo3mDOEiO68grKysjI+PDz+12+3l5eWteWO3bt2mTJly8803t/KDQqGQwWBoS4ldypBxNGSc5HV7dq71bP++W9Grv8rp98CA0XUjRixtiF9SIv1xCx+eyK7KEi7PZL0s6CSeGr5sbYCD1gacc845jtu5isaX7byC0Gg0lpaWhp+63W4Mj8ScoDMYR15qHHmpq6FeKD3g3bdZtXzRdKP5pkFj+biRqyh36VH+wh5JYDQ5nV2WziamCTZtrIsGAIid8wrCzMzMFStWyI+DwWBpaWlmZmYkqoIIYBqttu9Iff9RdN1vfYf3e/du9Hy4YGTAN673sFcuGFGWPnRlhfDfIn7X2kCumU1KZ5PShLGpTIP7kQBAF3NeQTht2rT77rtv165dAwYM+OSTTxITE8NXloEORBC0uf20uf0sV90erC337tnYtPJjXdULM3oMvKXfSHbhqPWNhq+PSXM3hQ47+dgUYUIqm5DGLog/2/llAIBOobVBuHDhwnnz5nk8HkmSli5deu+99z755JMJCQnPPffchAkTevXqVVhY+MEHH5x1cg7Elsqeahp3jWncNZKz0bt/s2fvRt/i1/unZo8YdPH/TbioVpv4XZm0qpy/tl+q8/HxqcL4VDYhjfW14q8VADqt1i6f8Pl8brc7/FSn0+n1evlxbW3tkSNHevXqZTQaW/mpt9xyyzlNlsHk7DZwu906nU4QzjLWyX0e74Gtnj0bvPs2qeJTdH2HaXsO1mT3KfOKq8r5qjK+qpx7gnx8mjAmmV2cwvrbmNCpYxFftjbAQWsDLJ9om2h82VrbI9RqtVrtqedU2O12u90euZKgXTGtXj9wjH7gGJIkX9Fe74EtjZ8tDFaXanP7Tes5eFbPweqLs4+46Ptyvqacv7xXqvLwi5LZmBTh4hQ2LIFpcQNhAFA4XH0EThAEbfcLtN0voKm3cZ/Hf+SA9+D2+g9fCNZXmboPuK7n4Jv7D1WNTa700NoK6YdK/vv10sFGPsTOLkpmFyaxC5OEJH2sfwQAgHOHIIRTYFq9tudgbc/BRBRqqvMX7fXmb29a8QFTabS9Bl/Rc/D0AYOEC82OAK2v4usr+ev7pdmrQ3adnIhsVBIbGM9UmIAKAEqAIISzEOPi5et9E+eBsmJv/nbXxq/r//OiKilD23PwuF6DJ1/Ql6k1nOhAA99QxddX8TcPSMUOPtjOLkxiIxLZ8ETWDZcCB4COCkEIrcaYOj1XnZ5rnjCDh4L+4v2+gu2NS98LlBWp03K1uf2yc/r1yul7W08LETUFaFMV31jN/1XI798g+UN8eCIblsCGJwrDE1kyBlEBoMNAEEJbMFEln1CMm/IL7vf5Swv9xXtdG76q//AFQW/S5PbV5vQbl9tv0sAsYgIRlbn5lmq+uYb/v32hLdXcoGLDE9nQBDYkgQ22IxcBIJYQhHC+mEYrL9g3X0LEeaCyxF+011e81/Htx9zvVWf11Ob2s+f0uyqr59Xdjl9fuMjBN1fzrTX8b7ukbbXcoGKD7TTEzgYnsMF2jKMCQLtCEEJEMaZO6aZO6WYcfQURhRpr/cX7fEV7Gj/9R6CqVB5B1eb2zc7pn5trmpV7/E2HHXxbLd9eyxcekHbUkifEB8WzgXY2MJ4NtLN+Nlz4DQCiCEEIUSRa7Mcn2hBJXpe/aJ+/eJ/ju8X+0mdViWnanH6arJ6arF7ZSRnZZuHa7OPvqvLQzjq+o5avLOPP75YOOXiPONbfxi6IZ/1s1M/Gcsy4ghEARAyCENqJoDPq+g7X9R1ORDwUDBwt8B/e792/pWnFvyVnozqzh6ZbL01WL01WzyRrwuR0Njn9eNj5QrSnnu+t53vr+ev7+b56qvXxPlZ2gY31ldPRSulGJCMAtBGCEGKAiSpNdh9Ndh/5qeR1B8qKA0cL3NtWNXzyChFTZ/bQZPbQZHbXZPfVGuOGJrChCT9FXVOAChr53nq+r4G/sleSo7F7HOtrZf1srK+NhiUIqbjLGwC0DoIQYk/QGeTpNqZx1xBRqL7KX3LQf+SgY9WSQOkCwWzTZPVSZ+RpMvLU6d0FgylOTS2isdZHu+v43nq+p54vPcr31Ic0AvW3sT5W1tvKellYTwtlYQ4OAJwKghA6HNGWpLcl6QdeTHRiGmpJfqC0sGn3+kBZkWAwq9Pz1Bl56vQ8TUaeaE0kIruWxqey8ak/Rd0xF9/bQPvq+Z56vrhYOtDIm/zU08J6WVlvC+tloV5W1tPCjPgXANDl4dcAdGwnpqHSiMnyhlBjbaC0wH+00LV+ecPRAu73qVKzNZndNZk91Jk91EmZJAhElG5k6Ua6NP1nA6oHG/jBRn6ggS8+TPmNUn4jT9KzXhbqZfmp42iLzc8JADGDIASFES120WLX9btQfhpqqgscOxQoPeTZu7FpxQchR4M6tZsmvbs6PVed2V2dks3UGrllnJqGJ7LhiT9Fo8TpsJPnN9KBhp86jg1+XZ452D2O5cVR9ziWF8e6x1GWiYkYWAXopBCEoGxiXLwYF6/rM1x+KnndgWOHAseK/If3O9ctDVYdFa2J6rRcdWo3dWo3dWqOKiGNTtyjUWCUa2a5Zro846eUq2lwNImmoiYqcvC99XzxYamoicrcPM3AcuPk9scf9LEyA/4BASgf/h1DpyLoDNq8C7R5Fxx/LoWC9VXBiiP+o4Xu7WuCyxcFaytU9hRVSjd5xFWVkqVOzqJm6xK14vF0JPppoy9ERQ5e2MQPNVFhE/++nBc20VEXTzWwPDN1j2PdLSzXTDlmlm1itlPfuBMAOigEIXRqgqiyp6rsqeGhVO73BipLAmXFwfIjrvXLA+WHecCvSummTu2mTstRp2Rzs51Ouv+1VqQ+VtbH+rPh0aBER5z8kIMKG/khB19bQYcdUrGDC4yyzSzHzLJN8gPKMbMcMybmAHRQ+KcJXQvT6DSZPTWZPcNbJLcjUFYcqDgSKCv2bPveX37EKQjq5CxVcqYqKUOdkqVKylTFJ9NJV7NRCZQXx/LifjYlh4jqfHTYwYsd/LCT8hv516W82EHFTm5S/Swgs0ysm4myTCxO3R4/OACcDoIQujrBYNZ2H6DtPkB+6nA4jCIL1pYHK44EKkqca78MVpSEmmrFOLsqJUud0k2dkqVK6aZOzmQa3Sl3GK+leC0bktAyOCs9VOzghx38sJN21vIvS6QjTipxcoFRppFlmynLxDKNLNNI3Uysm5lS9bi5MUB7QBACtCQYTBpDD01mj/AW7vcGq0oDlUeDFSWePRuD334crC0XLXZVUqbcZVQnZ6qSMgRj3Bl2m6ynZD27MKllQNb76KiLH3HyEieVOPmuOipxSkecVO3lSTrWzUxZRpZlokzj8U5kupHF4zQkQOQgCAHOjml06ozu6ozuP22SQsHaikDl0WDlUX/xXtf65cGqUmJMlZiuSkxTJWaoEtNUiemqxHRBd5arvdm0ZNOyAfEtAzIgUZmblzhJzsjd9XzZUemIk465uTdEmUaWZqBMI0s3UrqRZRkpzcjSDSzFQFjoAXBOEIQAbSKIcs5R/wvD2yS3M1hbHqqtCNaWe/dvDq5aHKgqZaKosqeK9hR5nqpoT1EnZTDt2W9GrBaom4l1M9HFJ0WbO0hHXbzMTUedvNRFBxr4N8fomEs65uJ1Pko1sAwjZRhZqoHSjSxFTxlGlqynDCMz43wkwEkQhAARI4+pUrMxVSIKNdQEq48Fq48Fa8rcW78LVB0L1VUIRosqIVVlT1UlpIr2FFVCmsqecuaR1eYMKuplYb0sRCdlpF+iMhcvdVGpi1d4qNTFd9TSMZckP+a8ZTSG/5tpwrxW6KLwxQeILtGaIFoTtD0G/rRJkoIN1aHa8mBNebC2wr9rXai2IlhTTpyLJ9JRZU9RJaSJ9hSVLZEEsfUfpxEo28yyf74OMswVpFIXr/Qc/+9RJ99aQ6UuSd7C2M+iMc3AkvSUomepBkrSsUQ9Bl2hc0IQArQ7QVDFJ6vik7U9BjXfLLkdwdryYE15qLbCX3LQvX11sKZcctSL1kS5+yjaU1QJqar4ZDE+WTC0XOzYGsbTdyWJyBn4qR9Z5aEyN99ZSxUeqdxNVV5e76NEHUsxUKqeEvUs3UBJepaip1QDM4RYroZwJQFQKAQhQEchGMwag7n5Gkci4qFgqK4yWFMWrK0I1pT7i/cF6ypDdZXEJTE+WbQlq+zJqvhkMT5FFZ8sxie1LSBlJjX1trLeVjplTAYlqvLyCjeVe6jKw8vcdKiJr6ukcrdU4dJUegPeECXqWZKOUgyUqGOJOko1sETd8fhM0lGinmmwIAQ6HgQhQIfGRNXxWTk/J3mcobrKYF1VqK4yWFfhO7Q7WFcZqq0k4mK8HI3JYnyyypooWhNFW6IYF3/yNQHOiUqgNANLOz4H9me7cjgcZrPZF6JqL6/wUJWHqr288kSfssorVbip2kvVXm5SU7KeJeooSc+S9ZSgpQQdS9Ad3yg/ViMsoX0hCAEUSdCbhHSTOj2vxfYTAVkZqqsK1lX6i/eF6qtDDdWSq0m0Joq2BNGWpLImibZE0ZooxiepbEmtmcLaGlqRMowswyg/O3Xo1nip2survVTp4VUeqvFSfiNfV0mVHqnGSzVeXuMjk4qS9CxBRwk6lqg73pVM0FGCliXoKFFHdh2mv0IkIQgBOpXTBSQPBkINNaGG6lB9VbCuyn+0ILT7x1B9VbC+mgmiaE043n08HpDJojVBtNiZGOFfEXK89SE6XVISUa2Pqj28xks1Pl7toSovlTj5thqq8UpVXqrxUq2XBySy6yhBx+xasutYoo7kB/J/E7Rk11Gijlk0kS0fOicEIUCXwFRqVUKqKiH15JcktzPUUB2qrwzWV4fqqwPlh0P1VcGGaqmxTjCahbh40ZKosiYIFrvKmiBaE4Q4u8qWGKl+5MnsWrJr5Zg8bVh6Q1Tr5bU+qvVRjZfXeKnGS0ecfFsN1fqkWi/J/Ut36Pje7DqK17J47fEL4MVrya4je7MtJnQxuzAEIUBXJxhMgsGkTstp+QLnIUdDqLFGaqoN1tdIjbXegp2hxppQY12ovooYU9kSRUuCYLFLhjhmTxatCaLZJloTBLMt4l3JFnQipRtZ+hmHYYkoIB1Pyjof1fnk/1Kdlx9xUq2P6rySvKXWx4MShcPSpmU2Ddm0FK9lNi3ZNGTTMpuW4rVk0zCbFmcxOxsEIQCcBmNinE2MsxH1OPlFyes+PtbaVOeuKA1WHvXlbw856kMNNZKjgelNotkqWhNFs1WOxvaMyTC1QCl6StGfpX9JRN7QT0lZ7+P1Pqr3U52PH2igeh/V+6U6H9XLL/lJJ1L8ibC0ydmpIauWWTVk05JVw6wasmrJqiGrBndvVgD8FQFAWwg6g5CSpU7JIiLJ4TD//CaOIUe95GgINdSEnA2h+upgVakvf3vI0RBqqP4pJi12wWwVzTYxLl4wWcW4eCHOJposgsna/j+OTjztnNiTNQV+Cst6H5czssHPy9zyA6nBRw1+avDzeh9xOp6I4Wi0acmqIYuaDExMMkkWDbNqyaIhi5osGgzSxgCCEAAiTzTbRLPtFMOtREQkORtCjoZQY63kqA85GoINNdLRwpCjTmqqDzkbuMclmCyC2SbGxYsmi2Cxi2abYLaKcfGC0SKaLIIx7jyXgpynODXFqVk3k/zsLJV4Q9TgpwYfb/BTg5/qTzxo8PN8J3PX8ka/1OCnxuN/uDdElhPBaVGTVcssmp9iMk5z/NU4NcVpKE5NcRrMoT1fCEIAaG+CySqYrOrU7FO+ykNBydkYaqqTY1JqqgvWlEnF+0JNdSFno+RskDxO0WgRTBbBZBHNNsFkEYwW0WwVzNbjSWm2CnpT+/5Mp6UTmw/PUjg4OecuV9BkanlXy6BEjQFq9PN6HzX4qdHPT2Qklbp4UwM1+aneLzX5qSlATX5qCnBXgKxasmhYnJosJ9IxTv3TRrOaTGqKUzOrlsxq+Q/i8ycIQgDoWJioEi120WI/bQspFHI2Sq6m40npbJRcjf6SfMnVGHI2yiFKwUCLpBSMcYLRLJqsgtEiGM2CMU4wxrXb2crWUwkn5s0eH2w+e99X4tTopwY/l6Ox0U9NAS4/aPDzUhc5A+QIkCMg1fvkB9wRIGeArJrjiWjWkFlNNg2TI1NOSpuWjCoyqcmkZlYNmdRkUpFJ3QkXpXS4LwEAwFkIohgXL8bF02n6lETEgwHJ2Sg5G0KOesnZJLkaJZcjcKzI52yQXA7J1SS5mkKuJkGjFUxWwRQnGOIEY5xojBOMcYJBjkmzYDj+J3prRSJCYMen7TTb1qqh43rf8VCU/zT4eZP/eGrW+XiRg5wBcgbJGZAa/SRnpzPAmwLHu5gmNYtTU5yaTGpmUJFFQyY1GVRkA084CAAADwRJREFUUjGLhgwqMqjIqmFGNRlUZFZTnJrJGzuajlcRAMB5Yyq1fN+PM4//SW6n5GoM56L8IFhfJbmaJLdDcjskl0NyOygUPB6QBnPzgPwpMvVmZjAJemPHGZJtjbbFJxE1BcgZ4HJkNvrJFeTuIDX6yRkgV5Dq/bzYQe4guUPU6JecAXIHyRmgRj93h8gTJJuWjCp2Ih3lyGQ27fHsjFMzs5oMKjKqj0+7lVPWqGI6kaKxdAVBCABdl7yGkk66lGsLPOA/novhPy6H5HY0i0yn5HFwt0vyugSDSdCbmd4o6E3CiXRkBpOgD/8xMoOJ6YwkEZGSgjNMni7UbMM5zF3iRA2+49kp56g7SO4gb/CTK0juINX7+RGnvJHq/ZL8oMlPziD3Bskv6Y7eSAktT62eFwQhAMBZMLXmLKctwziXPE7J45TcTu5xSR6H5HFJbqfkcQYbayW3U/K45Abc45S87oZgQNAZBb2J6Q2CzsC0RkF+oDMeT1Ddiac6A9MZBJ1R0OnP6RaVHQ07RU+UWh+lDofDrIvwPB8EIQBA5DAmD5nS2UKTc+5yuUwGveR1Sx4n97gkr1vyurnHJfnc3OuW3M5gbQX3uiSvW/K4uffEdq+biSqm0Qk6I9MbBI2OafVMqxf0JqbRCVrd8cdaHdPoBa2e6Y1MrWVqtaA3M7WGqTvdXJfzhiAEAIgdQTwenOeCB/zc55F8bu5xS34P93m4zyt5nNznkXxeye0I1lZwn0d+yr0uye+joF/yOHnAzwN+QW9iajVT65jewFQaptEJOgNTa5hWL2j1TKVhOj3T6pkgCnoTiaKgM5CoYhrd8UDVGkgQBb0xtqs5IwhBCACgMHLHTjBZ2vZ2ye3kQT8P+LjHxQN+HvBJHhcP+rnPK3ndPBiQ3E5eX8UlSXI7SJIkr5uCAR7wSX4fBQOSz01SSHI7SRAEnYGJaqbRESN5ohDT6EhUMVGUp9oKOgMJIlOpmVpLRILBRERMrWUqNQmioDMQEdPqmKgiQcW0OiKSI5apNEytDe82qhCEAABdi5xGESBJcnDygJc4SR4nEXG/l0JBHgpxn4eIJK+bpBAPBnjAR0SS20lEPFDLgwGSQpLXTUTc5+GhEIWC3O8lIsnjJE484ONBP0lc8rroRPYTEYlq4/w3BGNcZH4EIkIQAgBAGwlCxDL1bORBXSJyul2RTUFCEAIAQMcX7hGyEI/4znFbLQAA6NIUEITBYHDlypWxrkJ5Nm/eXF1dHesqFMbtdn///fexrkJ51q5d63A4Yl2FwtTX169fvz7WVSjPqlWrfD5fZPepgCAsLS29//77Y12F8rz88str1qyJdRUKs3fv3ieeeCLWVSjP008/vXXr1lhXoTDr169//vnnY12F8jz00EMFBQWR3acCghDajPPID6YDAHQyCEIAAOjSEIQAANClsZiMno0ePdrhcKSkpLSmsdfr3bp160UXXRTtqjqZXbt2JScnJycnx7oQJWlqajpw4MCIESNiXYjCbN26NTc312azxboQJamtrS0pKRk8eHCsC1GY9evXDxgwwGg0trL99OnTf/Ob35y5TWyCcO3atXV1dQaDoTWNOedHjhzJzs6OclGdTXl5uc1m0+kiereSzi4UCh07diwrKyvWhSjM0aNHU1NTVSqsSz4Hfr+/uro6Pf0sd4CCFo4cOZKZmSkIrR3OzMnJycvLO3Ob2AQhAABAB4FzhAAA0KUhCAEAoEtDEAIAQJeGIAQAgC5N7GgXlCoqKvr888+rq6tzcnLYqW5/7HA4vvjii927d2dkZGBKpIxzvmrVqtWrVxsMBrvdfso25eXlW7ZsUavVcXERvoOJcpWXl3/22WdHjhzJyckRRbHFq5zzbdu2rVy5srS0NCMjQ61Wx6TIjsbv9y9btmzTpk2JiYlmc8v7qkuStG/fvlWrVm3btk0QBKzeCdu8efPXX38dCoXS0tLO0GzdunUNDQ04bjKfz7d06dItW7YkJyebTC3v9xQKhb777ruiE0Kh0Ol++50d70iWLVtmt9vvvPPOIUOGTJs27eQGlZWVOTk5V1111XXXXZeRkXH06NH2L7IDuummm/r373/33XcnJCR88sknJzcYN26cwWAwGAyvvfZa+5fXMW3dujU+Pn727Nljx44dOXKk1+tt0eDGG2/s06fPrbfeOm7cuPT09KKiopjU2aH4fL7/397dhjTVv3EA/+mmPUi6+ZTbUqdFmqigKJnMnKWV+NCyUZalLcNR0ZsiLSk0JipBphAUKVkGWkgPzJllbGkJvRAcmEmRloZrzZlzMpfmOTv/F4f/8N6s+8Zi5+Suz6uzwwXny4/f2bWdx4SEBIFAcOTIEW9v797eXpuC0dHRsLCwvLy8/Px8Pz+/M2fOUJKTbmQyWWBgoFQqDQwMrKqq+llZR0eHm5tbTk6OI7PR1szMTGxsrFAozM/P9/Hx6e/vty9ACG3bti01NTU1NfXGjRtL3ha9GmFMTMytW7cIgjCZTDwer7u726agrKxMJBKRy4cPHz59+rSjI9IP+YVuMBgIgmhtbd24caPFYrGp+fTpE4ZhKSkp0Aitdu/eXV5eThAEhmExMTF37961KRgaGrIui0SiU6dOOTQfLTU3N0dFRc3PzxMEUVFRkZmZ+Yvivr4+BoNhMpkclY6myHumBwcHCYJ4+/ath4fH1NSUfZnRaIyMjJRKpdAISQ0NDfHx8RiGEQRRWlq6b98+mwKyEX7//v33t0Wjc4RjY2NqtXrv3r0IIQ8Pj127dikUCpsahUJBFiCExGKxfYETUigU27dvZ7FYCKGsrKyRkZH379/b1PD5fPtDf87MYrE8efJELBYjhBgMhkgksp9LC2/C5XA4f/zNL38jhUIhEonIG+fFYvHTp08xDPtZsdls9vT0hEPKSqUyNDR006ZNCKGIiIigoCCVSmVfVlJScuLEiXXr1jk8IE0pFIqcnBzyi+sX3/Y9PT1dXV1Go/F3tkWjRvjlyxdPT0/rGSwej6fRaGxqNBqN9UEMPB5vbGzMoRFpSaPRWHeeFStW+Pr62o8bsDE+Pj4/P79wLv1i0IaHh+/duyeRSByVjr5sdkAMw3Q6nX1ZXl6eUCg8ePDg48eP3d3dHZuRdhbuoegnk62rq6u/v18qlTo2Gq3ZTDaz2WwwGGxqOBxOXV1daWkpn8+Xy+VL3haNHomE4/jCq2MYDIb9j00cx61P1mEwGOQxwEWvqXEeNuPGZDJ/8SMdkHAcRwhZx23RyUbS6/XZ2dnFxcUJCQmOy0dXNjsgQmjRcTt+/Pjk5OSdO3cuXryoUqmc/GjEv+6hZrP55MmTra2t//2xYc7gXyfbqlWrxsbGyJrGxkaJRKLT6Zb2nD8ajXtAQMD09PTs7Cz5UafTcTgcmxoOhzM+Pr6wwMm7IPrnmOA4rtfrf31ZGkAI+fv7u7q66vV68qNOp1t00AwGw44dO/bs2XPu3DnHBqQpmx3QxcXFfidFCAkEguzs7Pv37w8MDPT09Dg2I+0sHDS02GR78OCB2Wyuq6uTSqVyuVytVp89e9bhMWnHZrK5u7v7+vouLHBxcbF2ytzc3MnJyZGRkaVti0aNMDg4OCQk5Pnz5wghHMeVSmVKSgpCCMOwyclJskYoFHZ2dpLLnZ2dQqGQorA0IhQKX7x4MT8/jxB69eoVi8UKDw9HCM3MzJAnk4E9Nzc3gUBgP5csFsu3b98IgkAIGY3GnTt3pqSkVFRUUBiVVoRC4bNnz8jlzs7OxMRE8sjnwp+wViaTaXZ2Fm7XEQgEAwMD5DFkrVY7ODgoEAgQQrOzs9PT02TB5cuXyUsf169f7+fnt3XrVopD04DNZEtOTib/9hiNRvsT9mq1mslkLv0/wO9fb/MHNTQ0cLncK1eu5OTkxMTEkBenvXz5kslkkgVDQ0NsNru4uPjChQteXl5v3ryhNC8tWCyWxMTEjIyMmpoaPp9fW1tLri8oKJBKpeRyfX19UVERl8sVCARFRUV9fX3U5aWLjo4ONptdVVV17NixoKAg8kI+8helRqMhCEIsFnt5eRX939WrV6mOTD2j0cjn8yUSSXV1tbe3d3t7O7l+8+bNNTU1BEE0NzcfOHCgqqqqrKwsIiIiOzsbx3FKI9NCYWFhfHx8bW1tXFycda+srq5OSkqyqZTJZHDVKGliYoLL5Uql0srKShaLpVKpyPVRUVHknRJNTU2HDh2qrq4uKSnx8fG5dOnSkrdFu7dPKJVKlUq1du1aiURC3q779evX9vb2wsJCsuDjx4/Nzc0Wi2X//v1hYWGUhqULs9nc2Nio0WiSkpLS09PJld3d3Uwmk3yPo1KpHB4ettanpaWFhIRQk5VOent75XK5p6dnQUGBv78/QshkMrW0tOTl5a1evbqtrU2r1VqLeTxeRkYGdWHpQq/X3759e3p6Oisry/riRrlcHhoaGhkZaTAY2traPnz44ObmFhcXl56eDicvEEI4jre0tAwMDERHR+fm5pIH9Pr7+z9//pyZmbmwUq1WT0xMpKWlUZSUXrRabVNT08zMjEgkio2NJVc+fPgwIiIiPDxcq9W2tbWNjIysWbMmOTk5MTFxyRuiXSMEAAAAHIlG5wgBAAAAx4NGCAAAwKlBIwQAAODUoBECAABwatAIAQAAODVohAAAAJwaNEIAAABODRohAAAApwaNEAAAgFODRgjAX2NqamrR9/8BAH4HNEIA/gL19fUbNmxgs9kBAQErV668fv061YkAWD5o9GJeAMCiZDJZWVlZYWFhY2Ojr6/v6Ogom82mOhQAywc8dBsAWtPpdMHBwfn5+Tdv3qQ6CwDLExwaBYDWurq65ubmjh49SnUQAJYtaIQA0Jper0cIbdmyxWUBmUxGdS4Alg84RwgArZGnAx89ehQdHW2zEgDwR0AjBIDWBAIBg8F4/fq1SCSiOgsAyxOjvLyc6gwAgJ9isVjj4+PXrl1zdXXl8/nu7u5DQ0MTExN+fn5URwNgmYCrRgGgOwzDysvL6+rqTCYTuaaysvL8+fPUpgJg2YBGCMDfYW5u7t27dz9+/ODxeFwul+o4ACwf0AgBAAA4Nbh9AgAAgFODRggAAMCpQSMEAADg1KARAgAAcGrQCAEAADg1aIQAAACc2v8AGWVWsFsNZxwAAAAASUVORK5CYII=", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot( ϵ->cond([1 0 ; -1/ϵ 1 ]), 0, 1/2, yaxis = :log10, title=\"Condition numbers\", xlabel=\"ϵ\", label=\"κ(L)\")\n", "plot!( ϵ->cond([-ϵ 1; 0 -1+1/ϵ]), 0, 1/2, label=\"κ(U)\")" ] }, { "cell_type": "markdown", "id": "ffa35135", "metadata": {}, "source": [ "Notice this is not a problem with $M(\\epsilon)$, but with our choice of LU decomposition. That is, without pivoting, our algorithm is unstable.\n", "\n", "3. Using PLU decomposition, we arrive at \n", "\n", "\\begin{align}\n", "LU = \\begin{pmatrix} -\\epsilon & 1 \\\\ 1 & 0 \\end{pmatrix} \\begin{pmatrix} 1 & -1 \\\\ 0 & 1-\\epsilon \\end{pmatrix} \n", "\\end{align}\n", "\n", "(with $p(1) = 2$ and $p(2) = 1$). These factors are again perfectly conditioned:\n" ] }, { "cell_type": "code", "execution_count": 22, "id": "638b19b4", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeXgT1d4H8HMme9MtTds0bSm07HvLVnaQHcHrRVQE5bqg4uUVl6vCxV1crldceK8r3AsvKrjghoqCQtkLZS8WylqW0i1tk27Zk5nz/jE15hYsBdtO0nw/j49PMj3N/Dok883MnHOGMsYIAABAqOKkLgAAAEBKCEIICOfPn8/Jydm3b19JSYkkBSxatGjevHm+p0VFRXPmzHnvvfeu+It5eXlz5sz55JNPWrK6wLVs2bK5c+cWFhZKXQjAtUMQgpTMZvPChQsTExNTU1OHDBmSmZmZlJTUqVOnl156yWq1tmYln3zyyX/+8x//wlauXLl582bfkl9++WX58uXHjh1r8IsXL15cuXLlnj17WqnQALNp06bly5dXVFRIXQjAtZNLXQCErry8vClTply8eDE2Nvauu+7q3LkzY+zcuXM//vjjM888k5OTs379eqlqCw8PHz16dO/evX1LNm3a9Pjjj7/99ts9e/b0b6nX60ePHt2lS5dWrxEAmgeCEKRRWVl5/fXXFxUVzZ49+913342IiPD9yO12f/DBBz///LOE5XXs2HHr1q1NaZmZmdnElgAQmBCEII0XXnihqKhozJgxq1at4rj/OkWvVCofeuih22+/3X8hY2zfvn0HDx50Op3t27cfN25cVFSUf4OLFy+Wl5d37NgxOjr6wIEDe/bsYYwNHjx40KBBl669trZ2w4YNRUVFiYmJkydPjo6ObtDA4XDk5+frdLq0tDRCSH5+flFRkbiWgwcPim26dOkSERFRW1t7+vTp+Pj4du3a+b+CzWbbtGnT+fPn5XJ53759hw0b1uDP/OWXX3iez8jIcLlcGzZsOHv2bExMzIQJExITExvfdGJtUVFRnTp1qqio2Lhxo8lkSk1NnTRpklar9TVzuVxHjx6NjIzs3Lmz/6+XlJSUlpampqbGxMSIS/Ly8txud//+/Z1O54YNG86fP280Gm+44Qbfq+3fv3/fvn2MsbFjx3bv3v2yVfE8n5WVdfz4cbVaPWHChNTU1EvbMMYOHTq0f/9+q9WakpIyfvx4nU7n3+DUqVN1dXW9e/eWyWRbt27Nz8+PiIi4++67xZ8ePHgwPz/fZDLpdLrk5OThw4f7/70A144BtDqHwyHuwrKzs5vSvrCwcPDgwf7v2+jo6I8++si/zSOPPEII+eyzz2666Sb/lnfccQfP8/4ts7Ky4uPjfQ2ioqJ+/PHHlJQUhULha5Obm0sImTZtmvi0b9++l352Nm/ezBj74YcfCCEPPvig/yq++eab2NhY/8bp6emnTp3yb5OQkKDRaI4cOdKhQwdfM7Va/fnnnze+NY4ePUoImTp16urVq8PCwny/265du+PHj/uanTp1ihAyadKkBr/+7LPPEkJWr17tW9K+fXtK6eHDh1NSUnyvZjQa8/Pz6+rq/vSnP/kWymSypUuX+r/a9OnTCSHffvttenq6r5lcLl+8eHGD9Z45cyYzM9N/m0RGRq5atcq/zejRowkhW7du7dOnj9gmOTmZMWY2m0eNGtVg+yuVSrPZ3Pi2AmgKBCFIYMeOHYQQvV4vCMIVG1ut1q5duxJCbr311n379p08efKdd94JDw+nlH777be+ZmIQpqamduvW7ZNPPjl06NDq1auNRiMhZOXKlb5mZ8+eDQ8Pl8vlL7/8ckFBQX5+/sMPPxwREREREdFIEO7Zs2fu3Lli4G36lbgXvjQIt2/fLpPJ1Gr1kiVLTpw4cfDgwbvuukvcp/vvuBMSEhQKRXJy8syZM3/66ad9+/Y9+eSTHMdptdry8vJGNogYhO3atdNqtc8991x2dva2bdtuvPFGQsjQoUN9za42CFNSUmbPnr158+Y9e/bMnj2bEDJ48ODbb7+9R48ea9euPXTo0Ntvv63RaORy+enTp32/KwZhYmLiiBEjsrOzCwsLv/jiC/GgdsWKFb5mZWVlRqNRJpPNmzdvx44dx48f//DDDxMSEiilGzZs8DUTgzAlJWX06NFr1qzZvXu3+LXgvvvuI4TMnDkzJyensLDwl19++fzzz6dNm4YghGaBIAQJ/N///R8hZNiwYU1p/MYbbxBCRo0a5Z+a4nCFzp07+472xCA0Go01NTW+Zt9//z0hZMqUKb4l99xzDyHkmWee8V/FHXfcQQhpJAgZY6+//joh5O23325Q3qVBKJ6Mfe+99/ybiUH11FNP+ZYkJCQQQu65555LK2lwsNuAGISEkDVr1vgWut1u8XiuqKhIXHJVQUgIufPOO31LeJ7v2LEjIcRgMPhvz7///e+EkCVLlviWiEFoNBqtVqtvodiH1mg0ejweccmcOXMIIa+88op/JYcPH5bJZOnp6b4lYhD26dPH7Xb7t+zYsaNGo/G9GkDzwvAJkEBNTQ0hxL+DTCO++eYbQsjChQsppb6FM2bMSEtLO336dF5enn/j+++/PzIy0vd03LhxlNJz586JTxlj69atUygUDz30kP9vPfbYY9f6pzRUXFy8b9+++Ph4MXF9xAj5+uuvG7R/4okn/J+OHz+eEOIruBHioaTvqUKhuO6665r4u5f1t7/9zfeY47gRI0YQQubMmeO/PUeOHHnZVTzwwAP+l+sGDx48bNiw0tJSMRE9Hs9nn32m0WgabOf09PShQ4fm5uaWlpY2qEShUPgv0el0Lpfr8OHD1/anATQOnWVAAuJO0263N6Vxfn4+IaRfv37+CzmOS09PP3v2bH5+vv8FPPEkqo9arY6KiiorKxOfmkwmi8XSvn37BhfwevXqJZc3z2dBHGjYs2dPlUrlvzwjI4PjuJMnT/I8L5PJxIUymaxTp07+zQwGAyHEV3AjunTp4v/NwPe7JpPp2ipv0KcmLi5OXMulCy9dxaXXUDMyMrKzs48dOzZixIhTp07ZbLb4+HjxYNSfxWIhhIjdc3wLGwxQIYTcfffdBw4cGDx48JgxY8aOHTtu3Lj+/fs3+PMBrhmCECQgdrA8e/ZsUxpbrVaO48RdsD9xv19XV+e/0L/ziIjjOPbrhLriIH3/njIiuVweExNTVVXV1D+g0WovuwqVShUdHW2xWKxWq6+/q1KpbBDAYs9S1oQZgC/7lxJCBEG4hrIppRqN5tJXu+zCS8u79F9H3ALi1qiuriaEVFVVLV++/NJV63Q6j8fjv6TB1xRCyLx583Q63dKlS7OysjZv3rxo0aLk5OTXXnvN/5gY4Jrh1ChIYPDgwXK5vKio6MSJE1dsHBERIQjCpXOXiIdN/ifurig8PJxc7oDG6/Wazeamv04jxPO95eXlDZY7nc6qqiqO48QaWoEYWjzPN1jeElP2XPr3ihtZ3Bri/3v06GH5HeIZ18bNnDlz7969paWln3322e23324ymW6//XZpB5tCm4EgBAlERUWJgxxeffXV32vj9XrFB+KJsgMHDvj/lOf5Q4cOEUJ69erV9PUaDIbY2NiSkpIGO+4jR45cGhgNKJVKcrlcaUCs5+jRo06n03/5gQMHGGPdu3f3nRdtaWJnnEsj6vjx482+rkuv3omjLcWt0bVr17CwsGPHjonXhv8Ig8EwY8aM1atXv/nmm4yxr7766g++IABBEIJUFi9eHB4e/uGHH77xxhuXnmrbvn37vffeKz6++eabCSGvvfaa/0m/Tz755MKFC927d7/0elIjKKXTpk3zer1Lly71Xy72CG2cOCTg4sWLjTczGo3Dhw+vqKhYsWKFbyFjTIz8W265penV/kFardZgMOTn51+4cMG38MCBAz/99FOzr2vZsmX+56h37dqVk5OTnJwsjv5UqVSzZs3yer1PP/30pb/blCPUS9uI/xwul+sP1Q1ACME1QpBK165d16xZc9tttz3++ONffPHFzJkzxc4aZ8+eXbduXVZW1tixY8WW991337Jly3bs2DF9+vRHHnlEr9dv3Ljxueee4zjuzTffvNoeE88888zatWv/+c9/ymSyW2+91e12L1++fOPGjZGRkQ6Ho5Ff7NevH8dx//73vxUKRUpKikwmmzJlSlJS0qUtlyxZMnLkyMcee6y6unrq1KlWq/Xtt9/+4YcfOnTo8Oijj15VtX/QjBkz/vWvf02dOvX555+Pj4/Pycl56aWXOnXqJI6saEYqlWr8+PGLFy/u0KHD3r17H3/8cULIq6++6jv8/cc//pGVlfXOO+9cuHBh9uzZXbp0qaurKygoWL9+fUFBgXhw34iUlJTp06dPmjQpLS1NpVLl5uY++eSThJBbb721ef8QCFHSjdwAYLm5uRMnTmwQZkqlcubMmf6jtktKSsaMGePfJj4+/osvvvB/KXEc4bp16xqsIiYmRqfT+S/Jzs72DzC9Xp+VldX4zDKi5cuX+3dubGRmmQ0bNjTIyCFDhpw7d86/jTizTINqxYteDzzwQCMbzTezTIPl4giNtWvX+pbU1NSII/NEMpns5Zdf/r0B9Q1ebeHChYSQBtPciCc8b7rpJt8ScRzhDz/84D+VnVKpfP311xu8YElJyY033tjg3zo8PNx/04nVNthQjLEePXo02HFFRES88847jWwlgKajDHeoB6mZTKacnByTySTOb5KZmXnp5J+EkLy8vEOHDjkcjrS0tBEjRjTo0Gg2m2tqahISEhp0pzx//jwhxH8aM0KIw+HYvHlzSUlJQkLC2LFjw8PDCwsLeZ73zZDpdrsvXLgQHh7un3wiq9UqXngzGo0ajcZutxcXF0dHRzfoOelyuXbt2nXmzBmVStWnT5+MjIwGGXDhwgVBEBrMyelwOEpLSyMjIy/tOenjdruLiorCwsLEq4A+FoulurraYDD4D+ljjInDGLRa7ejRo5OTk6uqqqqqquLj433ddgoLC71erzirauOv5nK5iouLxZOu4hKTyWSz2ZKSkhQKRXZ29okTJzQazZgxY35vxtTCwsKcnByLxRIZGZmSkjJgwAC1Wu37aUlJidPpTElJuXQ0S0lJyeHDh00mE8dx7du3HzhwYKt1O4I2D0EIAAAhDZ1lAAAgpCEIAQAgpCEIAQAgpCEIAQAgpCEIAQAgpCEIAQAgpCEIAQAgpCEIAQAgpCEIAQAgpCEIAQAgpDV/EO7Zs6ct3S3z2u73DSJsvT8C0x9eM2y6PyIEP7bNH4S7du3avHlzs7+sVGw2m9QlBDFsvWvG83yDW/tC09ntdmThNQvBjy1OjQIAQEhDEAIAQEhDEAIAQEhDEAIAQEhDEAIAQEhDEAIAQEhDEAIAQEhDEAIAQHBgbpfzxMFmf1l5s78iAABAS6j68h0qk6u79W/elw2gIHzkkUd27doldRUNCYLAcc1/3BwZGZmVlUUpbfZXBgBok2x7Nnguno5/dGmzv3IABWFOTs78+fN79eoldSGtITMzk+d5uTyAtj8AQMDyFJ2p+WFV/ENvUKW62V88sHbE3bp169+/mY95AxOOBQEAmkiw15n/72XdLfPl8ckt8froLAMAAAGMMcvq1zS9h2j6Dm+hNSAIAQAgcNVu/Ji5XVF/mtNyq0AQAgBAgHLk7bHt36y/60nCyVpuLQhCAAAIRN7youq1/9Lf9TQXHt2iK0IQAgBAwGEuh3nl4sipdylTurT0uhCEAAAQYBgzf/RPVed0bebEVlgbgvAa5efnDx48mDFGCHniiSd2797t/9OjR48OGzZM/CkAAFyV2h8/Yk5b1J/vb53VIQiv0ZNPPrlgwQJxOOCGDRsKCwv9f9qrVy+dTrd+/XqJqgMACFaOI7vsh7bq736KylpppDuC8AqysrK+/vprQoggCG+99dbOnTsJIcXFxbt27Zo6dWojv3jHHXd88MEHrVQlAECb4CkqqP7yXf2c51q6g4w/BOEVZGdnb9y4kef5e++9d//+/UOGDCGEbNq0aeDAgUqlspFfHDly5NatW10uV2tVCgAQ3PjaqsoVL0RPn6dITG3N9QbWFGsNrDwlvHxYaM01Tkmh/xrScLSKx+OZNWtWdHT0f/7zH3EC7uPHj3fs2LHxl0pMTCSEnD9/vmvXri1ULQBAm8E8bvPKF8KHXq9JH9HKqw7oIJyRxo02tuqcnNHKy6xuzZo1RqPx5ptv9t2Gwul0qlSqK76aWq222+3NXCIAQNvDWNVnb8ljEiLGzWj9lQd0EGrlJC1C+smpp06dunz58r59+44dO/b6668nhBgMhpMnTzb+W263u6amxmg0tkqNAABBrHbjam9FcdyDS4gUNyTANcIri42NjY2NXbVq1d13311SUkIIGTJkyOHDhxs0Y78Snx45ciQlJSUhIaG1ywUACCr2w9vtB7bE3r+YKq98pq0lIAibavz48XfcccesWbN4nh8xYkRVVdXZs2fFHykUitmzZyt/tWDBAkLI999/P3PmTElLBgAIdK6zx6q/ek9/b6t2E20goE+NBoJnn33W9/iNN97wPX7kkUeWL1/+6quvEkIuPTp0uVyff/75li1bWqdIAIBg5K0ssax6Sf+XRQpjBwnLuOojwq1bt27atKklSgkuDz74YFRU1O/NHXP27NnFixcnJSW1clUAAMFCsNZULnsm8vo7VV3Spa3k6o4I9+7de8MNN3Tu3PnSY6BQo1KpFi1a9Hs/7d69e/fu3VuzHgCAIMI87soVL4T1G6UdPEnqWq7miNDlcs2bN++RRx5puWqCS01NTeOzibrdbgyfAABoSBAsH78qjzFETpotdSmEXFUQvvjii3/+85979OjRctUEkdzc3EmT6r/IzJ49+/Tp074fffbZZ0uXLiWEmEymgQMHut1uaUoEAAhI1euWCQ6bbtZjkgyWuFRTT40eOXLkxx9/zMnJ+fLLLxtvaTKZduzY8cQTT4hPFQrFgw8+qNfrr7iK4LpXw1NPPbVo0SJx0u2vv/76wQcf9P3o6NGjBQUFhJB27dplZmZ+9NFH995776Wv4HK5eJ5vtYIl4XK5Gp+IDn4Pz/Nut9s3hwNcFZfLJZPJaGDsZINOS39s7Vu/dJ3Ji5r7stvLE2+L7wOVSuUV3wlN+ph5vd777rvvnXfeacrWoZSq1eroX2m12qB+O/om3eZ5fsmSJTk5OYSQwsLCffv2TZ48+Yq/fvvtty9fvrzFqwQACAbOA1nOfT9H3f0Mpw6TupbfNOmIMDs7u6CgYOnSpUuXLi0sLDx//vytt9768ccfX3aasfj4+MGDBz/11FNXW0pg5mV2dnZRUdGNN9541113aTSaxx57jBCyefPmQYMGKRSKK/76kCFDjhw5UllZGRsb2+BHKpVKLm/jw1fcbndT5qKDS/E8TynF1rs2Xq9XpVIF5i4l8LXcx9Z5bK/9p4/j/uc1eVxgTbnVpB1xnz59Pv/8c/Hxli1bKisr77///lbYidt2/1iX9UVLr8Wfqls/3S3zGyx0uVw33XRTamrqW2+9JX60Tp48mZaW1pQXDAsLi4uLO3Xq1KVBCAAQOtznT1g+fTP23uflhnZS19JQk8JMp9ONGzdOfFxeXr5hwwbf0xYVNnCcqmtGK6zIR3a5qQ3WrFmTlpaWnJzs+4LpdDr9zxJzHCcIv90lQ/wu73uqVqudTmeLlQwAEOg8pefNK16IuWOBskMgjiu76qO64cOHx8XFtUQpl6IKpVwv/RH0TTfdtGzZsvT09HHjxk2bNo0QYjQa8/PzfQ0SExPFOUhFxcXFycnJ4mPGWEVFhXhLJgCAEOQ1l1YuezrqpgfU3fpLXcvlXXWftJSUlPHjx7dEKQErOjpap9OtXr36vvvuu3DhAiFkyJAhhw4d8jWYMmXKu+++a7VaCSEnTpz47rvvpkyZIv6ooKBAoVB06dJFksoBAKTF15gr338ycvzMsIxRUtfyu9A5+wo0Gk1YWBghZMSIEfPmzZs/f74gCMOHD7fZbKdOnRLbPP/88ykpKcnJySkpKaNHj37ppZdGjar/J1+3bt3MmTPRCR4AQpBgq618/0ntkMnaYVOkrqUxbbzX4h/nGxBJCFm8eLHv8aOPPrps2TJxGu7IyMhVq1YJglBXVxcVFeVrw/P8hx9++M0337RmwQAAgUBw2is/eFrTe2jE2FulruUKcKRyjf7617/Gxsb695HhOM4/BQkhFy5ceOKJJzp16tTq1QEASIm5XeZ/P6fs0C1yyp1S13JlOCK8RgqFopFJt0VpaWlNHGUBANBmMI/bvHKxTJ8QfdNfpa6lSXBECAAAzYbxXvOqV6gqLOa2RwJkKtErQhACAEAzEQTL6teojNPf+XfCyaSupqkQhAAA0BwEwbL6Nea0x/xlURClIEEQAgBAMxAEy+rXBIdVP+c5Kr/yPMwBBUEIAAB/jCBYPnldsNfp73k26FKQBFSvUY1GM378+DZ/QwaRXC7HKHsAaAvEFLTW6Oc8RxVBef/RAEqdDRs2OBwOqatoyGq1hoeHN/vLKhQKBCEABD1BsKx5XbAFcQqSgApCtVqtVqulrqIhuVweEREhdRUAAAGH8V7L6teYy6G/9/lgPCPqE0BBCAAAwYLxXsuHrzCvV3/3M0GdggRBCAAAV4t5PZZVLxO5Qj/nSSoL+hzBZSoAALgKzO2sXPY0VYfp/7KoDaQgwREhAAA0nWC3Vi5/RmHsoLv1oWCZQe2KcEQIAABNIlirK95dqGzXuS2lIMERIQAANAVfVV7x3iJNn2FRN9wjdS3NDEEIAABX4C0vqvjgyYjR08NH3ih1Lc0PQQgAAI1xF54y/+f5qBvmhA0cK3UtLQJBCAAAv8t1+ojl41d1tz2q7jFI6lpaCoIQAAAuz354e803H+jvfkaZ2kPqWloQghAAAC7Duu3rum3fxP71HwpjB6lraVkIQgAA+G+M1az/P+exnPiH35Dp4qWupsUhCAEA4DeM91o+/qfXYoqb/zqnjZS6nNaAAfUAAFBPsNfZP3yR8Xzc//wzRFKQ4IgQAABEXnNp5fJnudRe+hltauKYK0IQAgAAcV84YV6xOGLCTNZ3dEilIEEQAgCA48iuqi/ejpn5N3XPzLq6OqnLaW0IQgCAEMZYXdZa6671cX/9hyIpTepqpIEgBAAIUczrqVr7L0/x2fhH3pRFx0ldjmQQhAAAoUiw1ZpXvshFRMc//CZVqqQuR0oIQgCAkOMpPmte8UJY5oTICbNCrWvMpRCEAAChxXFkZ9UX7+pu/h9N+gipawkICEIAgJDBWN2WL6y71sfOXaxs10XqagIFghAAICQwl8Oyeolgq43/2//KInRSlxNAMMUaAEDb5zVdNL0xX6aLi3vwn0jBBnBECADQxonj5aNumKPNnCB1LYEIQQgA0HYJQs0Pq+yHtsXev1iZ0lXqagIUghAAoG0SbLXmD/9BmGB47F9ceLTU5QQuBCEAQBvkOnvM8tGr2swJkZPuwEjBxiEIAQDaFnH60B3rdDP/pu4+UOpqggCCEACg7RCsNZY1rzO3M/6xd2RReqnLCQ4IQgCANsJ1Js+y+p+aPsOibryPyrB7bypsKQCA4CdOGbP9G93Mx9TdB0hdTZBBEAIABDfBWm1ZvYTxXpwOvTYIQgCAIOY6nWtZ87p28KTICbMIh8nCrgWCEAAgKDGvp/bHD+2HtsXcsVDVqbfU5QQxBCEAQPDxlF2wrH5NrjcanniP00ZKXU5wQxACAAQVxqw7vq3d9GnU9Xdqh14vdTVtAYIQACBo8NUVljWvM487/pG35LGJUpfTRiAIAQCCgyN3Z9WX74YPn4p+Mc0LQQgAEOgEp636y3c9RWfiHnhJkdxJ6nLaGgQhAEBAcx4/UPX5/2r6Do9/7B2qUEpdThuEIAQACFCC3Vq9bpm74GjM7Y+pOqdLXU6bhSAEAAhEzvx9VWvfVvcYaFjwHlVppC6nLUMQAgAEFsFurVm/0nXqcMztj6s695W6nLYPQQgAEEAcuTurv3pP3XuIYcH7VKmWupyQgCAEAAgIfF1V9Zfvek0X9fc9r0zpKnU5IQRBCAAgNcZse3+qWb8qfOjkmNkLqVwhdUGhBUEIACAlb0Vx1dq3mdMWO3exsl0XqcsJRQhCAABpMI+7Lmutddf6yPG3hY+8kVAqdUUhCkEIACABZ/6+6i/fVaR0MSx8Xxahk7qckIYgBABoVXyNuWb9Sve549G3PqTu1l/qcgBBCADQWhjvte38rnbTZ+Ej/qSb8Qg6xQQIBCEAQGtwnjxU880Hsui4+EfeksclSV0O/AZBCADQsvjqipofVrnPHY+aercmfYTU5UBDCEIAgJbC3K66LV9Yd34XPuJPulsfxr0jAhOCEACgRTiP5VR/9b4ipYvh8Xdkunipy4HfhSAEAGhm7ounq7/+gHg9MbMXKlN7SF0OXAGCEACg2fCW8pofVrnOHImc/Bdt5gSMkQ8KCEIAgGYgOKx1m9facjZqB08yLPo3pw6TuiJoKgQhAMAfwnivfe/PtRs/VnXtZ1j4gSwyRuqK4OogCAEArp3zWE71N8tlMYbYB15RJKZKXQ5cCwQhAMC1cJ3OrVm/igiCbsbDuI98UEMQAgBcHff54zU/fMhXV0ROnh2WMQo9YoIdghAAoKk8pedrf1rjPn8icsJM7eCJhJNJXRE0AwQhAMCVeU0Xazd/5jpxKHz0tJjbn8AcMW0JghAAoDF8VXntps8cv+wOHz5V9/RKqtJIXRE0MwQhAMDlec1ldVlrHUd2hY/4U8LTKzE0sK1CEAIANOStKK7b9Jnj2N7wYVMSnvwPp42UuiJoQQhCAIDfeMoK67I+dx7bpx0yOeGplVxYuNQVQYtDEAIAEEKIp/R83ZYvXCcOaYdNSXh2FafWSl0RtJKmBmFRUdEXX3xx7tw5pVI5duzYyZMnt2hZAACtxl14su7nT90XT0dcN113y0NUqZK6ImhVTQ3CvLy8wsLCbt261dXV3XPPPY899tjjjz/eopUBALQsxpwnDtZt+YI3l4WPuTnmzicxKCI0NTUIJ0+e7DsK1Ov1//73vxGEABCsGHPm76396RPm9URcN13TbzSV4TpR6BWJ3hQAACAASURBVLrqf3uXy7Vjx45+/fq1RDUAAC2KuZ22PRvrtn8ti4qLnDhL3SMTE6TBVQThmTNnBg0aZLVaMzIyNm/e/HvNLly48PPPPxcUFNSvQC5/6aWXkpKS/milEnE4HDIZZlG6Rth614znebfbzRiTupCgZLfbCSH0vxNOqDE7cza4DmxWdOytnfm4PKmTQIjd4ZCoxsDVxj62arWa47jG21xFEKamphYUFJjN5gULFtx9991ffvnlZZvpdLpu3brdeuutviUGg0GlCtaLz263O3iLlxy23jXjeZ5Siq13bbxer0ql8gWh5+Jp2/ZvXKcOa/pfF/vwm3K9UdryAlwb+9jSJhzxX0UQymQynU6n0+kWL16ckZHB8/xlvzVERkZ269ZtxowZV1FpAJPJZG3py1Erw9b7I7D1rpm46Sghzvy9ddu/5c2l2qHX626Zj0GBTRGCb7ymBmFtbW1kZP3cCrt3705JSQm1LQUAwYLZ6+pyfrRlfy+Pjg8f9WdN76HkSifHIJQ1NQjnz59//Pjx9u3bl5WVnThxYs2aNS1aFgDANXAXnrLt+t6etzus91D93c8o23WWuiIIAk0NwhUrVhw6dKioqCg2NrZfv37h4TjDAACBgnk9zqM5ddvXCXUW7ZDJUePf1cYamnJxCIA0PQjlcvmgQYMGDRrUotUAAFwVb2Wpbc8G296fFYmpEaP+rOkzjHCczWaTui4IJhhDCgDBh/FeZ/4+2+4N7ountZnj4x9dKtcnSF0UBCsEIQAEE29FsS3nJ/v+zfLYRO2Qyfp7nsG8aPAHIQgBIAiIVwGtu3/0FJ3RpI+I/esrCmMHqYuCNgJBCAABzVNWaN+/WbwKGD70enXvIZgXFJoX3k8AEIgEa4390Db7/s28rUabOdHw+Nuy6Dipi4K2CUEIAAGEeT3OYzn2/VmugqPqXplRN9yj6pyOebGhRSEIASAguC+etu/fbD+0XR6XpB04Nmb2QqrSSF0UhAQEIQBIyVteZD+0zX5wK+Fk2oFjDY+/I4uOlbooCC0IQgCQAF9dYT+03X5om1BXpUkfGTN7oTKli9RFQYhCEAJA6xHsdc5je237szxFZ9Q9B0VNvkPdfSDhMIM/SAlBCAAtTrDVOvL2OHJ3uC+cVPcYFDHqz6pu/TEKAgIE3ogA0FL4uipn3h577k5P4SlVt37azIn6Oc9hIhgINAhCAGhmgq3Wmb/PnrvTXZCn7NhbO3CsZs6z6AIKAQtBCADNw1tZ4vhltyNvj9dUqO6ZqR0yWX/XUzj+g8CHIASAP4Axd+FJR94e59E9gsOm7jU4csJMVZd0XP+DIII3KwBcNeZxu88dcxzNcRzZReVKdc9B0bc+rErtgSlgIBghCAGgqfjqSmf+Pufx/a7TvyiS09S9hsTNf10ea5S6LoA/BEEIAI0SBHdxgfPYXuexvV6LSdW5r7r7wOhbHpJF6qSuDKB5IAgB4DL42irXiQPOEwedJw/JouPUPQZG3/RXZftuhOOkLg2gmSEIAaAe83rcZ485TxxwnjzEV1equqSru/aPuvE+WZRe6tIAWhCCECC0MeYpPe86nes8ech99pjc2EHdrb/u1oeU7brg4A9CBIIQIBR5K0tcp3Jdp484Tx/hNFpV577azIkxsxdymnCpSwNobQhCgFDB11W5C446Tx12nTzE3C5lx17qLhmRN9wjjzFIXRqAlBCEAG0ZX13pKshzFeS5Co4K1mpVpz6qzukRo2+SxydLXRpAoEAQArQ13opi19mjroKj7oI8weVUdeyl6tg7fNgURWIaBrwDXApBCBD8BMFTet51Lt9dkOcqyCOcTNWptyqtV+TYW+WGdlIXBxDoEIQAQUmw1rgvHHedP+E+f9xdeEoWHadK7a7uMRDX/ACuFoIQIEgIgqf8oufiade5fPfZo3x1pdyYqkrrGTF6mrJDD04bKXV9AMEKQQgQuASnzVN4ynX2mPviaffZY1xYhDKthzK5s3bIZGVyJ1zwA2gWCEKAAMLXmD1FZ9xFZ9wXT7svniEep7J9d2WHbuEj/qScvYBTa6UuEKANQhACSImvMXuKTrsvnnFfPO25eJp53PKE9sp2ndS9h6rHzwpP6YzDPoCWhiAEaFX+yee+cJJynDyhvSIhJSx9hOKGexSGFDH5eJ53u91IQYBWgCAEaEGM93rLizxFBe6iM56iM57iAi4sQpHcSdGuU/iwqYoZj+BmRgCSQxACNB9B8JpLPaXnPaUXPGXnvaXnveYyWUyCMjFVkdxJM3GWIrkTFxYhdZUA8F8QhADXjq8xe02FntILHvH/xQWcOkw81anpMUg+5hZFQnuqUEpdJgA0BkEI0CSM93orS7zlRd6KEm95kafsvLeskKrDFAntFYmpyvbdtEMmKQwpVKmWulIAuDoIQoBLMMZXVXgriz3lxd6KIm95kbeimK8xy6Lj5HFJivhkZUqXsEHjFcb2uGkRQBuAIIRQJ9itXnMpby7zmks9ZRe8ZYWe8iIqk8n1RnlCiiKhvSq1p0yfgJOcAG0VghBChsB7qyt5i8lrLuMtJq+lzFtR7C0vJozJ4xLlccny+GR1j0GKUdPkcUlUpZG6XABoJQhCaHMY42stvKXMazZ5LSbeXOa1mHhLGV9j5sKj5foEWYxBrk9Qd06XDbleEZ/EhUdLXTEASAlBCMGKedx8dQVfY+arKviaSq/FxFtM4v+5sAgx7WQxBmVq97D+o2UxCTJdHJXhDQ8ADWG/AAFNsFv5mkq+uoKvsfDVFXx1JV8rJp+ZeVyyqFhZlF6mi5NF6RWJaZreQ2QxBnlMAi7mAUDTIQhBYszl4GvMvLVaqKvia6v4GjNfU8lXV/I1Zr66gsoUsii9LDpWFqWXRccpU7rIomNl0XGyKD1uPAQAzQJBCC1LcNqEumq+rlqwVvO1FqGumq+rEmotvLWGrzULddWEk8mi9LLwKC4iRhapk0XGKLpkyKJjZZF6mS6eKlVS/wUA0MYhCOGaCIJgr+VttYKtTrDVCnbxQY1gFR/XCvY6wVbH22rrlCouQieLiObCo2WRMVxEtDK5ExcZIwuP4iJjZBE6nMYEAGkhCKEe87gFh405bYLTJjhszGETHFbBYWMuuyA+ttXVB561RnDaOG0kp43kwiJl2ghOG8mFRXDaSHlcEhcWydUvibQJJDIak0oDQEBDELYtAi84HcxlZx634HIwp5153MztFJw25nYxj0tw2JjHxZyOX9POKjhsgtPGHDbCcZxaS9VhnEZL1VpOE85pwji1lmq0cr2R02jrk+/XzGtKObSurqX/YgCAP0j6IBRstV6LicoV4tUgKldSxW+XhTiNtq3dkk3gBaeDEEKYIDhthBDCiOCwEkIIY8xpY4wxh40QItjrCCGCw0YIExw2whhz2hkTmNPOBJ65nMzrYU6b4HYyj4s57czlZEzg1GFUqaEKJacOoyoNVaioSs2pw6hCRZVqTqOlEdFUFcZpwjhNONVoObVWTD4qV0i3UQAAJCN9EDqO7LLt2cC8HuZ2EUKY1808Lt9PxQDwPaUqDZXJfnuqUFG50u+pwj9EiUxGlb/ND0I5jqrDrrY8r9frlssJIYTnBZejwU+Zy0EE/r+WCDxz/tbM9+cwnmfir3McJ5ZBOU6tFZtxYeGEEEIpVWsppVSjJYSI9+vhNFpCqPiFgMbEU8pRdRjlZFSlpjIF1WipUkUVKk6tpSo1xskBAFwt6feb2qHXa4de38TGzOVg/G/Bwzwu5nX7PfX4hyjheeb2yyRBYE771ZbncDg0Gg0hhMhk3CXTblGVhnCy/1rCyajaL31lCvF2BFQmw6xdAAABSPogvCpUpfnv86QtPve/t65OE4E7qQIAtFmc1AUAAABICUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhDUEIAAAhTd7EdiaT6aOPPtq7dy/P86NGjfrrX/+qUqlatDIAAIBW0NQgzMrKOnHixIwZM+Ry+QsvvJCXl7dixYoWrQwAAKAVNDUIZ82aNWvWLPGxTqebNm0aghAAANqAa7lGWFBQkJyc3OylAAAAtL6mHhH6FBUVLVq0aOXKlb/X4NSpU999992mTZvqVyCXr1ixIjU19dprlJTVapW6hCCGrXfNeJ53u91er1fqQoKS3W7neZ5SKnUhQamNfWzDwsJkMlnjba4uCE0m0/jx4x9//PGpU6f+Xpv27dtPmjTpoYce8i3p0aOHXH7ViRs4IiIipC4hiGHrXRsxCDUajdSFBCWO48LCwhCE1yzUPrZXkU+VlZXjx4+/7bbbFixY0EgzlUplMBj69+//h2sDAABocU29RlhdXT1p0qSJEyc+99xzLVoQAABAa2pqEK5atergwYMrVqyI+ZXdbm/RygAAAFpBU0+NPvDAA3feeaf/Ely9AACANqCpQahWq9VqdYuWAgAA0Pow1ygAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0BCEAAIQ0udQFkFdyhZWnhEnJdHI7brSRaqWvCAAAAohXIDnl7KdiYeNFZtCQ9RObOSekj50n07kbUujGIvZmHj9zC8uMpxOTuYlJtHcMlbo0AACQTKGV/VTEfipmW0qE1Ag6MYm+MVg2NL75o0H6ICSE9I6hvWPoE304q4dsLRU2FrGbNgt2L5mQTCcm0XFJXKxa6hIBAKDl2b1keyn7qVj4qYhVudj4JO7P7em7QxUGTQuuNCCC0CdcQW5I4W5IIYSQs3VsczH7+jz7n92e1Ag6LomOS+RGGqkSlzUBANoWcYf/faGwtYR1i6ZTU+jq0bIMPeVa5cxgYAWhv7QIen83en83wjNZrpltLmb//IW/OYsNjKPjErkb2tMe0Th3CgAQrMwusqVE2FzMNlxkCo6MS6J/6cytHs1FKVu7ksANQh8ZJf1jaf9YurAvV+kkm4uFn4rZv34UwuRkYjKdkETHJHLhCqmrBACAK/EIZE85+7lI+LmYna5ho43cxGS6sC+XFiHlgU0QBKG/WDW5rSN3W0dCCMmzsJ+K2Tv5wh3b+HQ9HZfEjU+ig+KoDAeKAACBJL+abS5mm4qFHaWsSxSdmExfz5QNjafywLjUFWRB6E/sYvN4b87hJdkmtrlEmL9bOFXDMuPpuERuXBLtH4tIBACQRoWTbCsVNhezjUXMI7DhBu6WVG7VKE6vkrqySwRxEPpo5GRcEh2XJCMDf9v00zcLboENN3DjkujUFJoYhlAEAGhZvsOSzcXMd1jy9TguwA9L2kIQ+otTk1tSuVtSCfm1G9LmYrZwH58YRm9oT8clcsMTqFomdZUAAG3I2Tr2/QW2/qKwx8S6RdNxSfTVgbIg6uTf1oLQn6/fqUeQ5ZSzTcXCMwf5Y1VsqIGOS+LGJtK+Ma3UNxcAoI05V8eyStjmYralRNCr6fgk+lBP7utxXEQQdl1sy0Hoo+DIiAQ6IkG2uD+pdpOtJUJWCZu1Vah0sjGJ3NhEOi6JSttnCQAg8FU6yZYSQcw/u5eNTeImJtMlmfJ22uDef4ZEEPqLVpJpHbhpHQghpMhW/43m+UO8WkbHJtKxiXRMIhffklMYAAAEEauH7CxjWSXClhJ2to6NTODGJtL5PbleuuAOP38hF4T+krX0zs70zs6EEFl+NcsqZp+dZfN2e5K1YiLSUUacOgWAkOMWSE45Ew/+cs1sQCwdm8i9M5QbFBcoAx6aV0gHob8e0bRHNJ3fk/BMdrCSbSlh/zom3L6V7x6pHNeOH5PIDY2nGmwtAGijeEYOVbKtpeznQsV+s6drNB2bSJ/NkA03tP1dX1v/+66ejJJBcXRQHP17X87Jk6zztr3V8ucO8kcsbEAsvS6RG2OkmfFU0Ra/FgFASGGEHLWwLSVsSynbWSYkhdExifTeTvyXE9TRrT7PmYQQhI1Ry8jIeGFKR9ni/sTqIbtMbEuJ8GiOcKqGDTHQMYncdUbaLxZz2QBAMDlVw7aWsi0lbFupEKWk1xnprI7038MVYveIujohIpRSkCAImy5cQSYl00nJMkJIlYvsKBO2lLA5O4WLVjY8gV5n5EYbaXprzZUOAHBVCmrZtlK2tZRtLWFyjlxnpNe3o68Hf4fPZoEgvBY6FbmxPXdje0IIqXCS7aXCtlK24qRQ6mAjE7jrjHR0Iu2tQygCgJQuWNnWEra1lG0rZTwj1xnpdUa6uL/EM1wHIAThHxWnJjencjenEkKIyUG2lQrbStn7x4VKJxtp5EYb6SgjQhEAWsm5Ora9lG0rZdvLmNPLRidy1xnpU+lclyjsg34XgrA5GTRkRho3I40QQkrtZFupsL2UvX9cKHewEQncaCMdmUD76nFNEQCa05latr2UbS9l28uYVyCjjHSUkf69L9cNN21tGgRhSzGGkZkduZkdCSGk3EH2VgjZJnb/LuF0DRsUT8clcsMMdFB80MzFBwAB5Wwd21XGsk3s52Lm9LIRCdwwA324F9cvliL9rhaCsDXEa8gNKdwNKUS8P0ZOuZBtYo/k1N80apiBDjdwQTRBLQBIwhd+vhsbjUuiD/XkerahSV4kgSBsbXHq30Kx0kl2lgnbStnje/lzdWyogY5I4EYm0IFxVIVbZACEPIGRY1VsexnbWca2lwphcjrKSEcb6aK+XAd0eGk+CEIpxap/m/jU4iK7yoQdZeyxvUJ+NeunpyONdEQCNzSehgfhbO4AcG08AjlQyXaVsR1lQraJGTR0RAKdmkKXDJKnhCP8WgSCMFDEqMif2nN/ak8IITYvOVzJsk1syS98TjlL0dLhCXRcEr3OyMWqpS4UAJqb3UsOVbJsE9tlErJNzKihwxPoLance8Moxvm1AgRhINLKyfAEOjyBLuzLeQVyxMJ2lbEvzrIHdnmilXSYgQ5PoBOSKM6NAASvOg/ZW16ffOL9bIcZ6F86cx+N4nQqqYsLMQjCQCfnSP9Y2j+WPtyL8Ex2xMx2lrGfitgzB/kwOR1uqM/FHtEYqggQ6IpsLNvEdpWx7WXsQh0bHE9HJHDPZHCD4qga3QKkgyAMJjJK+sXSfrH04V6EENmJavFcCnvzqFDpZEPj6bAEbriBDsSHCiAwiL1ddpnq88/uZcMM3EgjvbMLlx7TNm9pFIwQhEGsWzTtFk3ndCWEEJOD7KsQsk1s0X4h11x/mmV4Ah1t5OJwWRGgFXkE8oulfpzDlhJBKSPDDdwwA30Ug/wCFYKwjTD4hir6XXj/6LQwdxev+/Wy4jAD7aHD5xCg+dV6yD6/C37tw+nwBDo1hb4xGLNaBwEEYRsU5utrQziekTwL21nGtpSwxYcFgbFhBm6YgQ410Aw97qoIcI0YISer2Z7y+iO/MgcbaqBD47mn07mBcTQMe9aggn+uNk5GSbqepuvp/J6EEHK+ju0ysT3lbNUpoaCOZejpUAMdGk+HGHAGFeAKbF6yv4LtNrE95cIeE4tQ0mEGOjSePtyL64WJ9YMZgjC0dIigHSLoHZ0IIcTqIblmlm1i/z4p3LOD913JGJ5AM3BjRQBCCCEldnawsr6ri+/q+y2p3NtDMH6p7UAQhq5wxW9nUInfNIbLTwhFNtY7pv6y4jADF4NRTRAyfCN3D1aynSbm9LIBcbR/LH2+n2x4Avpjt00IQqiXFkHTIuhfOhNCSIWT7DEJe8rZ678Is8x8WgQdZqCD42lmPO0She420NaUOcgek7DbxPaUs1wz6xpNh8bTicl0cX9M6RkSEIRwGXHq3+Z78wjksJntMbGNRey5Q0Ktmw2Op5nx3OB4mhlHo5RS1wpw9dwCOVTJ9pazvRVsTzmrc7PB8XSIgXtxADcwFrP7hhwEIVyBgiOD4uiguPrvxeKAxYOV7O1jwm0mJo7NEOe+yYxHN1QIXP5X+45YWIqW9o+lo430qXQOw4pCHIIQro7/gEWekRPV9TuX5SeEizbWJ4b2j6XDE+ioBC5eI3WtENrEAbUHKxte7VvYlxuRwEXjZAb8CkEI105GSU8d7amrv7JY5SI55WxvhbDypHD/Tl6vpv11iuGJQmY8TdfjtsPQ4gRGjlezfRUsu0R+sIo/XcP6xNDMeDolhb40gMM9jOD3IAih2ehUZHI7OrmdjBDCCDlRzbYXunOr2IqTwpla1iuGiqdYB8XRzuhxA82k2Mb2VbC9FWxfOTtYyRLC6KA4mh7F7u3O9Yvj8PULmgJBCC2CEtI9mibL+AciZMRv9sWfi9g/cuuHZ4gnUUcmcAacRIUmE8e/iic8s03M7GID4+gwA320Nzckvv6GnTabMyyM4tsWNBGCEFqD4tebSYlPyx1kXwXbVyGsOCnM3cVHKX87WOwXS7V4V4Ifr0Dyqtjecravgu2rYBesLF1PB8XRKSkY3gDNA7sckEC8hkxNoVNT6k+inq5h+yrY/gr25Tkhz8LSIumAWDowjg6IpX1xcTH0MEJO1bADFWx/JTtQwY5YWIdwOiiOZsbTh3pyvXS4exE0MwQhSIwS0iWKdomqn/jNI5A8CztQyfZXsGXHhdO1rKfut1zsoaMyHAC0RResbH9F/b/7wUqmV9GBcXRAHJ3WnusXSyMwsA9aEoIQAouCq7/58P3dCCHEI5BTNSzbxLaWsNd/+W2Ehvhf92jMiRqsqt3kQAXbZRIOVrL9FYyj9SfPH+nFZcZh7A20KgQhBDQFVz9CQ8zFWg85WMkOVLD1hey5Q0K1iw2IowNjqTg+rD36xwcwi4uImSce9jl5NiCWDoij93Xllg2niWH4twPJIAghmEQqyHVGep2xfqdZ6STijnXVKfbgbt7NE9/BYv9Y3BxAYpVOcrCSHfq1h6fFyfrF0oFx9LY0+kYml4p/HQgYCEIIYrFqceSiuEuVVbnIsSp2sJJ9e4E9e1AotrNeut9yEdNotbRqNzlqqY+9g5XMN0jmTyn0hX4czmNDwEIQQtuhU9XfWEp86svFzcXsn0eQi83v95JvXBJd2BfJB0EDQQhtVoNcLHMQcX/99Xn21AHB5mX9Y2mGnvbT036xtFMk9tpXVmRjh83ssLl+Dk+bl/XT0/6xdGoKfa4flxaJ7xYQlBCEECoSNGRKOzql3W+30ThUyQ6b2Zfn2ZMHhEonS9fTDD3tF0sz9LRHNAarEUZIQS0Tt5L4f46SDD3N0NPbO9E3B3NpuM4HbQKCEEKUQeN/fZFUuchhMztkZj8XsVePCIVW1lNHxZ1+up72iaFhIfBZ8QrkeDU7ZGaHzexwJcs1sxh1/UaY35PL0JMkLZIP2qAQ+HADNIFORcYk0jGJ9Tt6q4f8YmGHzexgJfvPSSG/mqWG03Q9zYitD4YYlbT1Ng/xzzxiYblmdtjMjlWxlPD6P/CGflyb+TMBGocgBLiMcAUZaqBDDfW56BFIfjU7XMlyLez7C8IRC4tS0r4xNF1P+sbQjFiaGhEcl8eKbSzXQo6YWa6Z5VpYiY310NF0PU2PoXd25vrE4ObsEIoQhABXpuBI3xjaN6Y+7Bgh5+pYrpnlmtmHp9nf9grVLtZXT9P1YjrSnjqqlklbMiGEeAVyooblmtkRM8u1sFwz4yhJj6HpejqtA32+P9c1ClPWATQ5CB0Ox4oVKw4cOFBcXPzRRx8ZjcYWLQsgkFFC0iJoWgS9qUP9EouLHDazI2a2o4y9fUw4VcvSImjfGOpLx9a51VStp/5oTzzbebyatdPSdD1N19PHenN9Y6gxrDXKAAguTQ1Cu92+d+/evn37fvjhhw6Ho0VrAgg6MSoyNpGO/fUSo1sg+VXsiIUdMbN/Fgm5ZqbgSJ+Y+n43fWJo92iq+MO9UsUD0yNm9ouF/GJhuWZW7mS9dDRdTwfG0fu6cr1icE8rgCtr6qdEr9d//PHHbrf7iSeeaNGCANoAJUfE4zDSuX5Jka0+rtYXsldyhfNW1jWK9tbRPnraN4b2adohY62HHLUwsXtLnoXlWVi0ivaJIX1i6G0d6T8GchgNCXAN8HURoDUka2myllz/62gNJ0+OVbEjFvaLhW24KBwxMxlHxIPF3jraJ4Z2jSSMkZM17BcLy7PUh2i5g/XU1R9T3pbG9YmhOvTqBPjDKGOs6a3dbrdKpSooKEhLS/u9Nnfddde6det0Op1vyVdffdW5c+ffax/grFZreHi41FUEK2y9pit1kGM1sqPV9Gg1PV7Dna4jhBCjhvSKZj2jhB5RQu9olhrOcMDXFHa7XaPR0ODoyRtw2tjHNiwsTCa7Qte15j8i7N69u0qlWrhwYf0K5PKUlJRmX0trioiIkLqEIIat10QREaRLPJn261Onh7c73TERuC/fteA4LiwsDEF4zULtY9v8QchxXFRUVCOHjABwRQqOaHDhAqBVXEXHterq6qqqKkJITU1NVVXVVZ1TBQAACExXEYT9+vXr3r27TqcbO3Zsx44dQ2QQxc6dO+vq6qSuIigxxn766SepqwhWJSUlubm5UlcRrHJyciwWi9RVBKuNGzdKXUJru4ogPHv2rMVPWFhIDM195ZVXDh8+LHUVQamqqur++++Xuopg9fPPPy9btkzqKoLVW2+9lZ2dLXUVQUkQhNtuu03qKlpbyN9pBiAg4dLDH4QNCE2HIAQAgJCGIAQAgJB2UU17XwAACbRJREFUdQPqm2LevHk//vhj8I6gb+DAgQOdOnWKjo6WupDg4/V6s7OzR40aJXUhQamkpKS6urpHjx5SFxKUcnNzk5OTY2NjpS4k+DDGtmzZMnbsWKkLaTbTpk2bN29e422aPwhPnTp1+PBhvV7fvC8rlYsXLxqNRrkcQ7quxblz51JTU6WuIijZbDar1WowGKQuJCiVlJTo9XqVChPQXYs29rFNTU3t2LFj422aPwgBAACCCK4RAgBASEMQAgBASEMQAgBASEMQAgBASJM9//zzUtcQoC5cuPDtt9+WlZWlpqZy3O9+Yzh+/PixY8c6dOjQiqUFgYKCgu+++66ysrJDhw6X3g2H5/ndu3dv27bNZDK1b9++kc0bInJycjZt2kQIMRqNl21QWFi4bt260tLSxt+NIYgxtn379q1bt4aFhV12vITFYsnKysrJyZHJZPHx8a1fYSC7ePFiU95XJ06cyMvLa0tdSRticDlZWVkxMTFz5swZOHDgxIkTeZ6/bLPS0lKDwRAdHd3K5QW47777Tq/X33fffenp6dOnT7+0waBBg/r373/XXXdlZGT06tVLvJlJyHrqqac6dOjwwAMPJCUlvfXWW5c22LZtW0xMzD333JOZmTlu3Div19v6RQasv/zlLz179pw7d25cXNxnn33W4KfHjx+PiIiYPHnyXXfdFRsbu2DBAkmKDEzbt2/X6XTi+2rMmDG/974qLy83Go1hYWGtXF5rQhBe3tChQ999913GmN1u79Chw8aNGy/bbPr06QsWLEAQNtC7d++PPvqIMVZbW2swGHbv3t2gwZkzZ8QHPM+np6dfdu8fIsrKytRqdUFBAWPs4MGDUVFRVqu1QZsRI0b87//+L2PM6XR27Nhx/fr1EhQakI4cORIVFWWxWBhj69atS0tLa/Cdtaampry8XHycl5dHCDGZTBIUGpBGjRq1dOlSxpjT6ezcufP3339/2WYzZsz4+9//3raDEOdYLsNsNu/evfvmm28mhGg0milTpqxfv/7SZmvXruU47sYbb2z1AgPa+fPnjx07dtNNNxFCIiIiJk6ceOnW841v5TjOYDC43e7WrjJg/Pzzz7169RJvZN2vX7+YmJgdO3b4N6ipqdm5c+f06dMJISqVaurUqZd9N4am9evXjxkzRqfTEUKuv/760tLS/Px8/waRkZFxcXHi44SEBEppKL/Z/NXV1e3YseOK76vvv//earWKzdowBOFllJSUKBQK3+WEpKSk4uLiBm3MZvMzzzzz1ltvtXp1ga64uFin02m1WvHpZbeeT3Z29t69e2fMmNFa1QWc4uLi5ORk39NLN1dJSQnHcb5rh41vz1Djv/XEz2wjG+fFF1+cMGGC/9YOZSUlJZTSxt9XNTU1Cxcu/OCDD1q9utYWujOHPfroo1u3bm2wcMiQIe+//z7P8/7XjWUymdfrbdDyoYceWrBgQVJS0oULF1q81sAzd+7cvXv3Nlg4duzYN954g+d5/94xl916otOnT8+YMeODDz5o3759C9Ya2BpsLrlc3mBziQ18bRrZniGI53n/6Q8v3Xo+y5Yt++6773bt2tVapQW6pryvHn744fnz5ycnJ5eVlbV6ga0qpINwzpw5DRaGh4cTQhISElwuV01NTVRUFCHEZDI16MtXVFT01VdfabXaffv2mUwmu90+d+7c5557LjExsdXql9aiRYusVmuDheLmMhqN1dXVHo9HoVCQy2090blz58aNG/fCCy+E8uEgIcRoNJaXl/uemkymBu+ihIQEnuctFos4f+/vbc/QZDQaT5w4IT4WBKG8vPyyn8EPP/zw5Zdf3rp1a1JSUusWGLjE95XZbBZPHV/6viovL//0008VCsXcuXMrKircbvfcuXOffvrpdu3aSVRyS5L6ImUgEgShe/fua9euFR/37t37008/ZYx5vd7KykrGmNVqXfurF198UavVrl27trq6WuK6A4PX601JSdmwYYP4OC0t7YcffmCMeTwes9kstiksLExLSxO7I4W4M2fOqNVqccucP39erVZXVFQwxux2e21trdimd+/en3zyCWNMEISMjIyPP/5YwoIDyrZt28SvrYyxXbt2xcfHi4+tVquvz9Hnn3+emJiYn58vZaEBqW/fvqtXr2aMCYLQv39/sYObby9nt9t9e7l//OMfKpVq7dq1bbWDN4Lw8lavXm0wGF5//fUZM2b06NFD/HTt37+fEOLxePxbZmdno9doA++//35ycvIbb7zx5z//ecCAAWJHvqysLLVaLTbIyMhITk6+/1fipzFkzZo1a8iQIUuXLk1PT3/44YfFhc8///yECRPEx59++ml8fPySJUtmzpzZrVs3h8MhXbEBZ+TIkZMnT37zzTc7duy4ZMkSceGcOXPuuecexlheXp5MJhs9erTvzXby5ElJ6w0gn3/+eVxc3JIlS2bNmtW1a1fxfZWbm0sIsdls/i3379/ftnuNYkD95fXp02fgwIFHjx7t2rXr22+/LXb9UCqVnTt3HjBggP9FHaVS2bFjx/T0dOmKDTgDBgzo3bv38ePHe/fuvXTpUvFuOCqVqlu3bhkZGYSQ8PDwIUOGJP6qU6dOoXyZ8Mb/b+/uVVQHAjAM5+zBaGPQSg2CRWJlIjaCgt1ewCIYsBRBK7HxBrSxE5FYiY21nbdhq1hqoShYCIraKJ7Cw7I3sCRk3qec6mNI+PIzmXx9ybK8Xq8ty2q1Wu+jS1GUdDqt67okSYZh5HK55XKp67pt28Fg0OnILlIul+/3+36/r9frlUrlPRgOhzOZTCKReD6fyWTSNM3vg800TSbwLZVK5fP5xWKhadpwOHxPiyzLuq5ns9mf6yR8Pp+mae+T15P4DRMAQGh8PgEAEBpFCAAQGkUIABAaRQgAEBpFCAAQGkUIABAaRQgAEBpFCAAQGkUIABAaRQi41OVyORwO7P0E/DaKEHCd6XRqGIaiKLFYzO/3dzodpxMBXibu/wgBdxqPx7VazbIs27aj0eh2u5Vl2elQgJex6TbgIrfbLR6PFwqF2WzmdBZAFDwaBVxkPp+fTqdqtep0EEAgFCHgIsfjUZKkYrH454dGo+F0LsDLeEcIuEgoFJIkaTQafX5+fg8qiuJcIsD7eEcIuMj5fI7FYqVSaTKZOJ0FEMXfdrvtdAYA//n9fkmS+v3+9XrVdT0QCGw2m91uF4lEnI4GeBZ3hIC7vF6vXq/X7XZPp9N7pNlsDgYDZ1MBHkYRAm70eDxWq9X9fldVVVXVjw/WtQG/hSIEAAiNy0wAgNAoQgCA0ChCAIDQKEIAgNAoQgCA0ChCAIDQ/gHYlT49kd5UcwAAAABJRU5ErkJggg==", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot( ϵ->cond([-ϵ 1 ; 1 0 ]), -1/2, 1/2, title=\"Condition numbers\", xlabel=\"ϵ\", label=\"κ(L)\")\n", "plot!( ϵ->cond([1 -1; 0 1-ϵ]), -1/2, 1/2, label=\"κ(U)\")" ] }, { "cell_type": "markdown", "id": "35702e56", "metadata": {}, "source": [ "## Condition numbers of linear systems of equations\n", "\n", "In order to fully understand the previous section, we need to make sense of the condition number $\\kappa(A)$. Recall that when we were considering functions $f : \\mathbb R \\to \\mathbb R$, the condition number was (approximately) the constant of proportionality between the relative error in the input to the relative error of the output. That is, if $x \\approx \\tilde{x}$ then, we have \n", "\n", "\\begin{align}\n", "\\frac{|f(x)- f(\\tilde{x})|}{|f(x)|} \\approx \\kappa_f(x) \\frac{|x- \\tilde{x}|}{|x|}.\n", "\\end{align}\n", "\n", "We now take this idea but the \"input\" is the vector $b$ and the \"output\" is $x$ solving the linear system $A x = b$. Again, the condition number measures how sensitive the solution $x$ is to errors in the data $b$. Roughly speaking, if the condition number is $10^\\rho$, we expect to lose $\\rho$ digits of accuracy when we go from $b$ to $x$. This made perfect sense on $\\mathbb R$ (with the absolute value) but now we have many choices for the norm. We shall fix the Euclidean norm for the remainder of the course:\n", "\n", "\\begin{align}\n", "|x| := \\sqrt{\\sum_{j=1}^n |x_j|^2}.\n", "\\end{align}\n", "\n", "We also need the corresponding *operator norm* on the space of matrices:\n", "\n", "\\begin{align}\n", "\\|A\\| := \\max_{x\\not=0} \\frac{|Ax|}{|x|}. \n", "\\end{align}\n", "\n", "Notice that, we have $|Ax| \\leq \\|A\\| |x|$ for all $x$. We are now ready to define the condition number of the matrix $A$ (or the corresponding linear system).\n", "\n", "Suppose we have the exact data $b$ and exact solution $x$ (solving $Ax=b$) and we have the approximate data $\\tilde{b}$ (for example, the floating point representation of $b$) and corresponding solution $\\tilde{x}$ solving $A\\tilde{x} = \\tilde{b}$. First, notice that $|b| = |Ax| \\leq \\|A\\| |x|$ and that $A(x - \\tilde{x}) = b - \\tilde{b}$ and so\n", "\n", "\\begin{align}\n", "|x - \\tilde{x}| \\leq \\|A^{-1}\\| |b - \\tilde{b}|.\n", "\\end{align}\n", "\n", "Therefore, we have \n", "\n", "\\begin{align}\n", "\\frac\n", "{|x-\\tilde{x}|/|x|}\n", "{|b-\\tilde{b}|/|b|}\n", "%\n", "\\leq \\frac\n", "{|x-\\tilde{x}|}\n", "{|b-\\tilde{b}|}\n", "\\frac\n", "{|b|}\n", "{|x|}\n", "%\n", "\\leq \\|A^{-1}\\| \\|A\\| \n", "\\end{align}\n", "\n", "We call the right hand side the condition number: $\\kappa(A) := \\|A^{-1}\\| \\|A\\|$." ] }, { "cell_type": "markdown", "id": "22fd3dbc", "metadata": {}, "source": [ "
Exercise. \n", "\n", "1. Use the definitions directly to show that the following matrices are ill-conditioned,\n", "\n", "\\begin{align}\n", "\\begin{pmatrix}\n", "1 & 0 \\\\ -\\epsilon^{-1} & 1\n", "\\end{pmatrix}, \n", "\\begin{pmatrix}\n", "-\\epsilon & 1 \\\\ 0 & \\epsilon^{-1}-1\n", "\\end{pmatrix}\n", "\\end{align}\n", "\n", "2. Show that $\\|A\\| = \\max_{|x|=1} |Ax|$, \n", "3. Show that $\\|AB\\| \\leq \\|A\\|\\|B\\|$ and hence $\\|A^k\\| \\leq \\|A\\|^k$, \n", "4. What is $\\|D\\|$ for diagonal $D$,\n", "5. $A$ is normal (i.e. $A^TA = AA^T$) if and only if there exists orthogonal $P$ (i.e. $P^{-1} = P^T$) and diagonal $\\Lambda$ such that $A = P\\Lambda P^T$. Show that the columns of $P$ are eigenvectors of $A$,\n", "6. Show that, for normal matrices $A$, the operator norm $\\|A\\|$ is equal to $\\max_i |\\lambda_i|$ where $(\\lambda_i)$ are the eigenvalues of $A$.\n", "\n", "
" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.11.6", "language": "julia", "name": "julia-1.11" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.11.6" } }, "nbformat": 4, "nbformat_minor": 5 }