Examples of LU decomposition using gaussian elimination and cholesky factorization.

LU_decomp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chris Sullivan | LU Decomposition | 10/30/217"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below are two methods for LU decomposition of a matrix, Gaussian Elimination (w/o pivoting) and Cholesky Factorization for symmetric positive definite matrices. Evaluate each cell using iPython to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg\n",
    "from scipy.linalg import hilbert\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def backward(R,b):\n",
    "    nrows,ncols = R.shape\n",
    "    if nrows != ncols: raise RuntimeError(\"This function only supports square upper triangular matrices.\")\n",
    "    x = np.zeros(ncols)\n",
    "    for i,row in enumerate(np.flip(R,axis=0)):\n",
    "        pos = ncols-(i+1)\n",
    "        x[pos] = (b[pos] - row[pos+1:].dot(x[pos+1:]))/row[pos]\n",
    "    return x\n",
    "\n",
    "def forward(L,b):\n",
    "    nrows,ncols = L.shape\n",
    "    if nrows != ncols: raise RuntimeError(\"This function only supports square lower triangular matrices.\")\n",
    "    x = np.zeros(ncols)\n",
    "    for i,row in enumerate(L):\n",
    "        x[i] = (b[i] - row[:i].dot(x[:i]))/row[i]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gaussian_elimination_outer(A):\n",
    "    L,U = [np.eye(A.shape[0]),A]\n",
    "    for k in range(A.shape[0]-1):\n",
    "        v = A[k+1:,k]/A[k,k]\n",
    "        A[k+1:,k:] -= np.outer(v,A[k,k:])\n",
    "        L[k+1:,k] = A[k+1:,k]/A[k,k]\n",
    "    return L,U\n",
    "\n",
    "def gaussian_elimination(A):\n",
    "    U=A\n",
    "    L=np.eye(A.shape[0])\n",
    "    for k in range(A.shape[0]-1):\n",
    "        for j in range(k+1,A.shape[0]):\n",
    "            L[j,k]=U[j,k]/U[k,k]\n",
    "            U[j,k:] -= L[j,k]*U[k,k:]\n",
    "    return L,U\n",
    "\n",
    "def linear_system_lu(A,b):\n",
    "    L,U = gaussian_elimination(A)\n",
    "    return backward(U,forward(L,b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probelm 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shown below is the vector two-norm of the residual between the true solution to the linear hilbert system, and the estimated solution via LU decomposition via Gaussian Elimination. For Hilbert matrices = 2, the error is on the order of machine epsilon, but increases with increasing matrix size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEKCAYAAADenhiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUVNW5/vHvK4o4GwU1AQkaWShRo9yOoD8HnFFRCI5g\nHCKIeoNDFBXjEOKERo25IJELioBGkCCajrYjioCiDAIKtiTIdWg0gjMIyvT+/tgHKcuu7q7uU3Wq\nup7PWrW66tSpqqeLot4+e++zt7k7IiIicdok6QAiItL4qLiIiEjsVFxERCR2Ki4iIhI7FRcREYmd\niouIiMROxUVERGKn4iIiIrFTcRERkdhtmnSApDRv3tzbtGmTdAwRkaIye/bsT9y9RW37lWxxadOm\nDbNmzUo6hohIUTGz9+qyn5rFREQkdiouIiISOxUXERGJnYqLiIjETsVFRERip+IiIiKxU3EREZHY\nlex5LiIijdKaNfDll9+/fPXV92//7GfQs2dOY6i4iIgUim+++WFhqK441HT/qlW1v06PHioudWVm\n3YETgG2B+9392YQjiUgxW78evv0WVq/+/s+GXF+5subisHp17bm23hq23Ra22y5cdtgBdttt4+3U\nS+p+qduaNs3521fQxcXMRgJdgaXuvnfK9i7A/wBNgPvc/TZ3fxx43Mx+BNwJqLiIlKr16+Gdd2DO\nnHB5551wVJBNMVi7Nt5Mm20GW2zx/S/6XXaBdu1qLwip9zVpEm+uHCno4gKMAu4BxmzYYGZNgKHA\n0UAVMNPMyt39rWiX66L7RaQUrF4Nb721sZDMmQPz5sHy5eH+TTcNfQxbbhn+Yt988/DX/w47hOsb\nttXlen3v32wz2KS0xk8VdHFx9ylm1iZt8wHAIndfDGBm44BuZlYJ3AY85e6v5zWoiOTHihWhcKQW\nkgULNjYnbbUV/OIXcPbZsP/+4fLzn4cveMmrgi4uGbQEPki5XQV0BC4GjgK2M7M93H1Y+gPNrC/Q\nF6B169Z5iCoi9bZs2feLyJw58O9/g3u4v3nzUDwuu2xjIdljj6JpNmrsirG4VMvdBwODa9lnODAc\noKyszPORS0Rq4Q7vvffDQrJkycZ9fvrTUDzOPHNjIWnZEsySyy01KsbisgTYNeV2q2ibiBS6tWth\n4cLvF5G5c+Hzz8P9m2wCe+4JnTtvLCL77Rf6R6SoFGNxmQm0NbPdCEXlDKBXspFE5Afc4fXXYdas\njYXkjTfCqC2AZs1gn33g1FM3FpJ99gkd71L0Crq4mNlYoDPQ3MyqgD+4+/1m1g94hjAUeaS7L0gw\npoikmz0bLr8cpkwJt7ffPhyBXHTRxkKy555hJJc0SgX9L+vu1Z5C6u4VQEWe44hIbT78EH7/exgz\nBnbcEYYMgRNOgDZt1D9SYgq6uIhIkVi5Eu68E26/PfSr9O8P114bTvyTkqTiIiL1t349PPwwXHMN\nVFXBKaeEArP77kknk4SV1imjIhKfl1+GTp3grLNg551D/8rf/67CIoCKi4hk6//+D047DQ4+OJyL\nMno0zJgBhxySdDIpIGoWE5G6+eoruPVWuPvuMMpr4MDQt7LVVkknkwKk4iIiNVu7Fu6/H66/PkzJ\ncvbZoci0bJl0MilgKi4iktlzz4XzVebPD81eFRVQVpZ0KikC6nMRkR96+23o2hWOOQa+/homTICX\nXlJhkTpTcRGRjT79FC6+GPbeG6ZOhT/9CSor4eSTdRKkZEXNYiIS1kMZOhRuvDF03F9wAfzxj9Ci\nRdLJpEipuIiUMncoL4crrwxrpRx7LNx1V1hgS6QB1CwmUqrmzoUjj4Tu3cPQ4ooKePppFRaJhYqL\nSKn56CPo0wc6dAhT4A8dGn4ed1zSyaQRUbOYSKlYtQr+/GcYNCj0sVx+OVx3XZgOXyRmKi4ijZ07\njBsHAwbA++9Djx5hcsk99kg6mTRiahYTacxefRUOOgh69Qrrq0yeDI8+qsIiOafiItIYvfdeKCgH\nHhiuP/BAWG74sMOSTiYlQs1iIo3NY4/Br38dmsNuuCEMM95666RTSYlRcRFpLNzDjMX9+0PHjjB+\nPOy6a9KppESpWUykMVi7Fvr1gyuuCFO1vPCCCoskSsVFpNgtXw7dusFf/wpXXQWPPAJbbJF0Kilx\nahYTKWZLloTZi998E4YNC3OCiRQAFReRYjVvHpxwQpho8oknoEuXpBOJfEfNYiLF6Kmnwhr2ZjBt\nmgqLFJxGU1zMbHczu9/MJiSdRSSnhg2DE0+Etm3htddg332TTiTyAwVRXMxspJktNbP5adu7mNlC\nM1tkZgNqeg53X+zuvXObVCRB69eHYcYXXRSOVKZMgZ/8JOlUItUqlD6XUcA9wJgNG8ysCTAUOBqo\nAmaaWTnQBBiU9vjz3H1pfqKKJGDlSjjrLJg4EX77W/jLX8I0+SIFqiA+ne4+xczapG0+AFjk7osB\nzGwc0M3dBwFd85tQJEEffxyGGs+YEU6SvPRSLTksBa8gmsUyaAl8kHK7KtpWLTPb0cyGAfub2TUZ\n9ulrZrPMbNayZcviTSuSC5WV0KlTWG9l4kS47DIVFikKBXHkEgd3/xS4sJZ9hgPDAcrKyjwfuUTq\n7cUXw/T4m28OL70Ev/xl0olE6qyQj1yWAKnzV7SKtok0fqNHwzHHhA77V19VYZGiU8jFZSbQ1sx2\nM7OmwBlAecKZRHJrw0zG554LnTvDyy9DmzYJhxLJXkEUFzMbC0wH2plZlZn1dve1QD/gGaASGO/u\nC5LMKZJT334bRoTddBOcdx5UVGgJYilaBdHn4u49M2yvACryHEck/z77DH71q3Duyi23wDXXqONe\nilpBFBeRkrZoUZgj7N134eGHoWe1f2uJFBUVF5EkvfJKOIfFHSZNCvOFiTQCBdHnIlKSxo+HI46A\nH/0Ipk9XYZFGRcVFJN/c4bbb4PTTwxDj6dPDJJQijYiKi0g+rVkDffuGDvuePeG552DHHZNOJRI7\nFReRfPnyy9Bxf999cO218NBD0KxZ0qlEckId+iL58P77obC8/TaMHAm/+U3SiURySsVFJNdmzw7r\n3K9aBU8/DUcemXQikZxTs5hILpWXw6GHhsknX3lFhUVKhoqLSK4MHgzdu8PPfx4mn2zfPulEInmj\n4iISt3Xr4JJLwqJe3brB5Mmwyy5JpxLJKxUXkTitWBHmCBsyBC6/HCZMgC23TDqVSN6pQ18kLitW\nhD6VWbPgnnvCWvciJUrFRSQOa9eGM+5nz4ZHHw19LSIlTMVFpKHc4b//O6y/Mny4CosI6nMRabhb\nb4URI8JZ9+efn3QakYKg4iLSEA8+CNddt3EFSREBVFxE6m/SpLAc8ZFHhvnCtHKkyHdUXETq4803\noUcP2Guv0IHftGnSiUQKioqLSLaqquC442CbbUIn/nbbJZ1IpOBotJhINr78Eo4/HpYvh6lToVWr\npBOJFCQVF5G6Wr0aTj4ZKivD7Mb77pt0IpGCpeIiUhfuYZjxpEkwerRmNxaphfpcROrihhtgzJgw\n3Pjss5NOI1Lwsi4uZna0mY0ws/2i233jj1U/ZraVmc0ys65JZ5FGZMQIuPlm6NMnnCgpIrWqz5HL\necCVwK/N7Ahgv4aGMLORZrbUzOanbe9iZgvNbJGZDajDU10NjG9oHpHvPPUUXHRRGB127706l0Wk\njurT57Lc3b8A+pvZbcAvY8gxCrgHGLNhg5k1AYYCRwNVwEwzKweaAIPSHn8e8AvgLaBZDHlEwiSU\np54Kv/gFjB8Pm6qLUqSu6vO/5ckNV9x9gJld3NAQ7j7FzNqkbT4AWOTuiwHMbBzQzd0HAT9o9jKz\nzsBWQHtglZlVuPv6hmaTEvXuu3DCCdC8OTz5JGy9ddKJRIpKrcXFzFqnbZqTtu0fG267+/sxZmsJ\nfJByuwromGlnd78WwMzOBT6prrBE/UN9AVq3Tv+1RCKffRaawVavhhdf1CqSIvVQlyOX0YADmRqb\nN9znwBEx5ao3dx9Vw33DgeEAZWVlnq9MUkS++SZMmb94MTz/fJjeRUSyVmtxcffD8xGkGkuAXVNu\nt4q2ieTG+vVwzjnhzPtx4+CQQ5JOJFK0Cvk8l5lAWzPbzcyaAmcA5Qlnksbs6qtDx/0dd4RVJUWk\n3urT55JRfftczGws0BlobmZVwB/c/X4z6wc8QxghNtLdF9Tn+UVqdc89cOed0K8fXHFF0mlEil4c\nfS4b1LvPxd17ZtheAVTU5zlF6uzxx+GSS0Jfy1/+onNZRGJQyH0uIrn36qvQsyd07Ah/+xs0aZJ0\nIpFGoZD7XERya9EiOPHEMG1+eTlsuWXSiUQaDRUXKU3LlkGXLuH6U09BixbJ5hFpZDSfhZSelSvD\nEcuHH8ILL8AeeySdSKTRUXGR0rJuHfTqBTNmwMSJ0KlT0olEGqW6DEVeThgJ9oO7AHf3bWNPJZIL\n7nDZZfCPf8CQIWF0mIjkRF1Gi22TjyAiOXfXXeF8lv79w/ksIpIzWTWLmdmPgLakTGvv7lPiDiUS\nu0cegSuvDGfe33570mlEGr06Fxcz6wNcSpjjay7QCZhOAUxWKVKjKVPC0sSHHgqjRsEmGiQpkmvZ\n/C+7lLAw2HvRiZX7A1/kJJVIXCoroVs32H13eOwxaKa15ETyIZvi8o27fwNgZpu7+9tAu9zEEonB\nRx+FdVmaNQvnsuywQ9KJREpGNn0uVWa2PfA48JyZfQ68l5tYIg20YgV07QqffBKaxdq0STqRSEmp\nc3Fx919FVwea2YvAdsBTOUkl0hBr18Jpp8G8efDPf0KHDkknEik52XTo31DN5v2AG+OLI9JA7nDR\nRaEZbMSI0CwmInmXTbPY1ynXmwFdgcp444g00C23wH33wXXXQZ8+SacRKVnZNIvdlXrbzO4kLOQl\nUhjGjIHrrw/Djm/UAbVIkhoy4H9LwjkvIsmbNAl694ajjgrNYVrwSyRR2fS5vMnGOcaaAC2Am3IR\nSiQrixfDKafAXnvBhAnQtGnSiURKXjZ9Ll1Trq8FPnb3tTHnEcnOypXQo0c4Unn8cdhuu6QTiQh1\nmxX58hruw93/HG8kkTpyhwsugDfegIqKcBa+iBSEuhy5bJgVuR1h+pfy6PaJwIxchBKpk6FD4aGH\n4KabNq4qKSIFoS5T7v8RwMymAB3cfXl0eyDwZE7TiWTy8svwu9+FFSV///uk04hImmxGi+0MrE65\nvTraJpJfH30UOvDbtAnDjzXLsUjByaZDfwwww8wei253B0bFnkikJmvWhKldvvoKnnsOtt8+6UQi\nUo1sTqK8xcyeBg6ONv3G3efkJlb2zOwQ4EzC79Te3Q9KOJLkQv/+MG0ajBsHe++ddBoRySCrlSjd\nfTYwO+4QZjaSMNR5qbvvnbK9C/A/hPNq7nP322rINhWYambdgZlxZ5QC8NBDMHhw6Gs5/fSk04hI\nDeoyFHmaux9sZsvZeBIlgAHu7tvGkGMUcA+h6W3D6zYBhgJHA1XATDMrJxSaQWmPP8/dl0bXewG9\nY8gkhWTePOjbFw47TMsUixSBuowWOzj6uU1t+9aXu08xszZpmw8AFrn7YgAzGwd0c/dBfP+Ezu+Y\nWWvgyw0j2qq5vy/QF6B169bxhJfc+/zzcKLkDjvAI4/AZpslnUhEalHnYTZmdqqZbRNdv87MJprZ\n/rmLRkvgg5TbVdG2mvQGHsh0p7sPd/cydy9r0aJFDBEl59avh1//Gj74IEztsrMGKIoUg2zGcF7v\n7svN7GDgKOB+YFhuYtWPu//B3V9JOofE6MYbw9n3gwdDp05JpxGROsqmuKyLfp4ADHf3J4FczhC4\nBNg15XaraJuUiieegD/+Ec49N0zzIiJFI5vissTM/hc4A6gws82zfHy2ZgJtzWw3M2savW55LY+R\nxmLRotAc1qED/PWvmkJfpMhkUxxOIywOdoy7fwHsAFwZRwgzGwtMB9qZWZWZ9Y5mXO4XvWYlMN7d\nF8TxelLgvv46dOA3aQKPPgpbbJF0IhHJUjbnuawCtgJ6AjcCmwFfxBHC3Xtm2F4BVMTxGlIk3MOQ\n4/nz4emnwxQvIlJ0sjly+SvQiVBcAJYTzkMRic/gwfDww3DzzXDMMUmnEZF6yubIpaO7dzCzOQDu\n/nnUFyISjylTwvQu3bvDgAFJpxGRBsjmyGVNdNa8A5hZC2B9TlJJ6fnwwzAh5e67w+jRmulYpMhl\nc+QyGHgM2MnMbgFOAa7LSSopLatXhyn0V6yAF16AbeOYUUhEkpTNrMh/M7PZwJGEecW6u3tlzpJJ\n6bj8cpg+HcaPh/btk04jIjGoU3ExMwNaufvbwNu5jSQlZcyYsFxx//5w6qlJpxGRmNSpYdvdHQ0J\nlrjNmRPOvD/8cBiUPtG1iBSzbHpNXzezX+YsiZSWzz4LJ0o2bx4W/to0q6WFRKTAZTUUGTjTzN4D\nvmbjei775iSZNF7r1kGvXmGE2NSpsNNOSScSkZhlU1yOzVkKKS0DB8Izz8D//i8ccEDSaUQkB7IZ\nLfZeLoNIiSgvD2ff9+4N55+fdBoRyRGdqSb5869/wVlnQVkZ3HOPZjoWacRUXCQ/VqwIHfhNm4aZ\njps1SzqRiORQrc1i0br0deLu7zcsjjRK7tCnD1RWwrPPQus6f6REpEjVpc9ldB2fy4EjGpBFGqu7\n74ZHHoHbboMjj0w6jYjkQa3Fxd0Pz0cQaaQmT4arroKTTw4/RaQkqFlMcqeqCk4/Hdq2hQceUAe+\nSAlRs5jkxrffhpmOV66El16CbbZJOpGI5JGaxSQ3LrsMXnstjAzbc8+k04hInmkossTvgQdg2DC4\n+uow/FhESo6Ki8Rr9my46KIwKuzmm5NOIyIJUXGR+HzySRgVtvPOMHasZjoWKWH63y/xWLcOevaE\n//wHpk2DFi2STiQiCVJxkXhcfz08/zzcd1+YO0xESlpRNouZ2e5mdr+ZTUjZtpWZjTazEWZ2ZpL5\nSs5jj4WVJPv2DbMdi0jJy3txMbORZrbUzOanbe9iZgvNbJGZDajpOdx9sbunf4v1ACa4+/nASTHH\nlkwqK+Gcc8K6LIMHJ51GRApEEkcuo4AuqRvMrAkwFDgOaA/0NLP2ZraPmT2Rdsm0bGEr4IPo+roc\nZZdUkyfDwQfDFlvAhAmw+eZJJxKRApH34uLuU4DP0jYfACyKjkhWA+OAbu7+prt3TbsszfDUVYQC\nAxl+LzPra2azzGzWsmXL4vh1Stfw4XD00WFk2CuvwK67Jp1IRApIofS5tGTjUQeEQtEy085mtqOZ\nDQP2N7Nros0TgZPN7F7gn9U9zt2Hu3uZu5e10Gim+lm7Fi65BC64AI45BqZPh5/9LOlUIlJginK0\nmLt/ClyYtu1r4DfJJCoRn38Op50WRoVdcQXcfjs0aZJ0KhEpQIVSXJYAqe0qraJtUigWLoQTT4R3\n34WRI+E3quMiklmhFJeZQFsz241QVM4AeiUbSb7z7LPhiKVpU3jhhdCJLyJSgySGIo8FpgPtzKzK\nzHq7+1qgH/AMUAmMd/cF+c4madxhyBA4/viwNPGMGSosIlIneT9ycfeeGbZXABV5jiOZrF4NF18c\nRoV16wYPPQRbb510KhEpEoUyWkwKySefhJFgw4fDNdfAxIkqLCKSlULpc5FCsWBB6Lj/8MNwtHKm\nZtIRkeypuMhGTz4ZZjbeaquwNHHHjkknEpEipWYxCR33d94ZjljatoWZM1VYRKRBVFxK3bffhnNW\nrrwSTjkFpk6FVq1qf5yISA1UXErZxx/DEUfA6NEwcCCMGwdbbpl0KhFpBNTnUqrmzYOTToJly2D8\neDj11KQTiUgjoiOXUvTYY3DQQWFp4mnTVFhEJHYqLqXEHW69FXr0gL33Dh33HToknUpEGiE1i5WK\nVavCEsRjx4ZzV0aMCIt8iYjkgI5cSsGHH8Jhh4XCcuut8OCDKiwiklM6cmnsZs0Kc4N9+SU8/ni4\nLiKSYzpyaczGj4dDD4VNNw1LEauwiEieqLg0RuvXwx/+AKefHjrsZ86EffdNOpWIlBA1izU2X38N\n55wDjz4azry/917YfPOkU4lIiVFxaUw++CCcGPnGG3DXXfC734FZ0qlEpASpuDQWr74K3buHIcf/\n/GdYPVJEJCHqc2kMHnwQOncOU+VPn67CIiKJU3EpZuvXw4ABcPbZcOCBYY379u2TTiUiomaxorVi\nBfTqFZrALrgAhgyBzTZLOpWICKDiUpyWLg1NX3PnhqLy29+q415ECoqKS7FZvBiOPRaWLAln3Hft\nmnQiEZEfUHEpJnPmwHHHwZo1MGlS6GcRESlARdmhb2a7m9n9ZjYhZdteZjbMzP5uZn2SzJcTL7wQ\nJp9s2jSswaLCIiIFLO/FxcxGmtlSM5uftr2LmS00s0VmNqCm53D3xe7eO21bpbtfCJwOHBt/8gSN\nHw9dukDr1mGOsL32SjqRiEiNkjhyGQV0Sd1gZk2AocBxQHugp5m1N7N9zOyJtMtOmZ7YzE4CKoBx\nuYufZ0OGwBlnQMeOMHUqtGqVdCIRkVrlvc/F3aeYWZu0zQcAi9x9MYCZjQO6ufsgoM491u5eDpSb\nWTnwaDyJE+IO110X1l/p1i2sxaI1WESkSBRKn0tL4IOU21XRtmqZ2Y5mNgzY38yuibZ1NrPBZjYc\nmJzhcX3NbJaZzVq2bFl86eO2di306RMKS9++MGGCCouIFJWiHC3m7p8CF6Ztm0yGopKyz3BgOEBZ\nWZnnKF7DrFwZpsp/4gm44QYYOFDnsIhI0SmU4rIE2DXldqtoW2n59FM48cQwCeW998KFF9b+GBGR\nAlQoxWUm0NbMdiMUlTOAXslGyrP33w8jwhYvDs1gPXoknUhEpN6SGIo8FpgOtDOzKjPr7e5rgX7A\nM0AlMN7dF+Q7W2Lmz4eDDgpn3T/zjAqLiBS9JEaL9cywvYIwjLi0TJsWmsK22CIMNdZyxCLSCBTK\naLHS9I9/wNFHw047hZMjVVhEpJFQcUnKiBGh+WvffeHll6FNm6QTiYjERsUl39zhppvC+SvHHhvm\nDGvePOlUIiKxKpTRYqVh3Tq4+OIwzPjss+G++7TAl4g0SjpyyZdvvgknR957L1x1FYwapcIiIo2W\njlzy4YsvoHt3eOkluPtuuOyypBOJiOSUikuuffhhWOCrshIefhh6VjsSW0SkUVFxyaWFC0On/aef\nwpNPhmHHIiIlQMUlV2bMgOOPh002gcmT4b/+K+lEIiJ5ow79XHjqKTj8cNhuu3BypAqLiJQYFZe4\njRkDJ50E7dqFkyP32CPpRCIieafiEhd3uOMOOOccOOyw0BS2yy5JpxIRSYSKSxzWr4crrgjnr5x2\nWui833bbpFOJiCRGxaWhVq+Gs84K569ccklY637zzZNOJSKSKI0Wa4jly+Hkk+G552DQILj6ai1J\nLCKCikv9LV0ahhrPnQsPPADnnpt0IhGRgqHiUh+LF4eTI5csCWuynHBC0olERAqKiku25s+Ho46C\nNWtg0iQ48MCkE4mIFBx16Gdrl13CAl/TpqmwiIhkoCOXbDVvDs8+m3QKEZGCpiMXERGJnYqLiIjE\nTsVFRERip+IiIiKxK8oOfTPbHbgW2M7dT4m2bQLcBGwLzHL30QlGFBEpaXk/cjGzkWa21Mzmp23v\nYmYLzWyRmQ2o6TncfbG7907b3A1oBawBquJNLSIi2UiiWWwU0CV1g5k1AYYCxwHtgZ5m1t7M9jGz\nJ9IuO2V43nbAK+5+OXBRDvOLiEgt8t4s5u5TzKxN2uYDgEXuvhjAzMYB3dx9ENC1jk9dBayOrq+P\nIaqIiNRTofS5tAQ+SLldBXTMtLOZ7QjcAuxvZtdERWgiMMTMDgFeyvC4vkDf6OYKM1sYR/gMmgOf\n5PD541IsOaF4sipnvIolJxRP1obk/GlddiqU4pIVd/8UuDBt20ogvR8m/XHDgeE5jPYdM5vl7mX5\neK2GKJacUDxZlTNexZITiidrPnIWylDkJcCuKbdbRdtERKQIFUpxmQm0NbPdzKwpcAZQnnAmERGp\npySGIo8FpgPtzKzKzHq7+1qgH/AMUAmMd/cF+c4Ws7w0v8WgWHJC8WRVzngVS04onqw5z2nunuvX\nEBGRElMozWIiItKIqLg0gJntamYvmtlbZrbAzC6tZp/OZvalmc2NLjcklPVdM3szyjCrmvvNzAZH\nMyS8YWYdEsjYLuV9mmtmX5nZZWn7JPZ+Vje7hJntYGbPmdm/o58/yvDYOs9AkaOcd5jZ29G/7WNm\ntn2Gx9b4OclDzoFmtiTl3/f4DI/N2/tZQ9ZHUnK+a2ZzMzw2L+9ppu+jxD6j7q5LPS/Aj4EO0fVt\ngH8B7dP26Qw8UQBZ3wWa13D/8cBTgAGdgNcSztsE+A/w00J5P4FDgQ7A/JRtfwIGRNcHALdn+F3e\nAXYHmgLz0j8nech5DLBpdP326nLW5XOSh5wDgf51+Gzk7f3MlDXt/ruAG5J8TzN9HyX1GdWRSwO4\n+0fu/np0fTlhMELLZFPVWzdgjAevAtub2Y8TzHMk8I67v5dghu9x9ynAZ2mbuwEbJkkdDXSv5qHf\nzUDh7quBcdHj8pbT3Z/1MHAG4FXCcP9EZXg/6yKv7yfUnNXMDDgNGJvLDLWp4fsokc+oiktMoilt\n9gdeq+bug6LmiKfM7Od5DbaRA8+b2exopoJ01c2SkGShPIPM/1kL4f3cYGd3/yi6/h9g52r2KbT3\n9jzCUWp1avuc5MPF0b/vyAxNOIX2fh4CfOzu/85wf97f07Tvo0Q+oyouMTCzrYFHgcvc/au0u18H\nWrv7vsAQ4PF854sc7O77ESYH/a2ZHZpQjlpF5zqdBPy9mrsL5f38AQ/tCwU9/NLMrgXWAn/LsEvS\nn5N7CU0z+wEfEZqbCl1Paj5qyet7WtP3UT4/oyouDWRmmxH+If/m7hPT73f3r9x9RXS9AtjMzJrn\nOSbuviT6uRR4jHAYnKqQZkk4Dnjd3T9Ov6NQ3s8UH29oPox+Lq1mn4J4b83sXMJEsGdGXzI/UIfP\nSU65+8fuvs7d1wMjMrx+QbyfAGa2KdADeCTTPvl8TzN8HyXyGVVxaYCorfV+oNLd/5xhn12i/TCz\nAwjv+afO7gkuAAAEwElEQVT5SwlmtpWZbbPhOqFzd37abuXA2dGosU7AlymH0vmW8S/BQng/05QD\n50TXzwH+Uc0+ic9AYWZdgKuAkzzMw1fdPnX5nORUWj/frzK8fuLvZ4qjgLfdvdo1pPL5ntbwfZTM\nZzTXIxga8wU4mHCI+QYwN7ocT5hU88Jon37AAsLoi1eBgxLIuXv0+vOiLNdG21NzGmFNnXeAN4Gy\nhN7TrQjFYruUbQXxfhIK3kdsXJCuN7AjMAn4N/A8sEO070+AipTHHk8YvfPOhvc/zzkXEdrUN3xO\nh6XnzPQ5yXPOB6PP3xuEL7cfJ/1+ZsoabR+14bOZsm8i72kN30eJfEZ1hr6IiMROzWIiIhI7FRcR\nEYmdiouIiMROxUVERGKn4iIiIrFTcZFGz8xWpN0+18zuia5faGZnR9dHmdkp0fV3G3JyppntV8OM\nvp3NzM2sT9r+bmb9a3ne7mbWvob7v/t96pizk5m9Fs3YW2lmA6PtJ+VjtmFpvDZNOoBIktx9WNzP\nGZ21vR9QBlRk2G0+YbLD+6LbPQnnQtSmO/AE8FZ1r1uP32c0cJq7zzOzJkA7AHcvR0uNSwPoyEVK\nmoX1QzIdLVwVrcMxw8z2iPZvYWaPmtnM6PL/Up7nQTN7mXAi4I3A6dERwenVPPd7QDMz2zk6s7oL\nKZNJmtn50fPPi15vSzM7iDDn2h3R8/7MzCab2V8srBNy6Ybfx8w2jR7fOXq+QWZ2SzU5diKcHIiH\naVfeivZPPbpLXWNnlZkdFp15PjJ6b+aYWU5nJZbioyMXKQVb2PcXctqBuv1V/qW77xM1M/2FMC/X\n/wB3u/s0M2sNPAPsFe3fnjBJ4apoHq8yd+9Xw/NPAE4F5hAm5Pw25b6J7j4CwMxuJpwRPsTMygnr\n2UyI7gNo6u5l0e2BAO6+NsowwcwuJhSvjtVkuBtYaGaTgaeB0e7+TeoOHiZdxMxOJEwh8wrwR+AF\ndz/PwsJjM8zseXf/uobfV0qIiouUglUbviDhuwkcy+rwuLEpP++Orh8FtI++1AG2tTALLUC5u6/K\nItd4woSHe0avcVDKfXtHRWV7YGtCEcuk2kkT3X2BmT1IaEY70MM6Hen73GhmfyPMedWL0DzXOX0/\nM2sL3AEc7u5rzOwY4KSUo75mQGvCGiIiKi4iNfBqrm8CdEr/6z4qNln91e7u/zGzNcDRwKV8v7iM\nArpHfSHnUs0XfoqaXncf4AtC81emHO8A95rZCGCZme2Yen9UPMcD5/vGyUwNONndF9bw2lLC1Oci\nktnpKT+nR9efBS7esIOZ7Zf+oMhywlKztbkBuNrd16Vt3wb4yMIU6mfW43kxsx6EJsBDgSFR81X6\nPifYxsOwtsA6QjFKNRJ4wN2npmx7hrCo14YZqvevSyYpHSouIpn9yMzeIBxV/C7adglQZmGlxLcI\nMzZX50VC81mmDn0A3P0Vd69uwbPrCasIvgy8nbJ9HHBl1In+s0zPGw2jvg3o4+7/Au4h9BelO4vQ\n5zKXMBDhzNRCZ2Y/BU4Bzkvp1C8DbgI2A94wswXRbZHvaFZkERGJnY5cREQkdiouIiISOxUXERGJ\nnYqLiIjETsVFRERip+IiIiKxU3EREZHYqbiIiEjs/j/CieyjV8faswAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8591fab7f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "norms = []\n",
    "sizes = [2,4,6,8,10,12,14,16,18,20]\n",
    "for n in sizes:\n",
    "    A=hilbert(n)\n",
    "    x = np.ones(n)\n",
    "    b = A.dot(x)\n",
    "    result = linear_system_lu(A,b)\n",
    "    norm = np.linalg.norm(result-x,ord=2)\n",
    "    norms.append(norm)\n",
    "fig,axes = plt.subplots()\n",
    "axes.plot(sizes,norms,color='red')\n",
    "axes.set_ylabel(r\"$\\vert\\vert$ residual $\\vert\\vert_2$\")\n",
    "axes.set_xlabel(r\"Hilbert Matrix Size\")\n",
    "axes.set_yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cholesky_factorization(A):\n",
    "    R=A\n",
    "    for k in range(len(R)):\n",
    "        for j in range(k+1,len(R)):\n",
    "            R[j,j:] -= R[k,j:]*R[k,j]/R[k,k]\n",
    "        R[k,k:] /= np.sqrt(R[k,k])\n",
    "    return R\n",
    "\n",
    "def linear_system_cholesky(A,b):\n",
    "    R = cholesky_factorization(A)\n",
    "    #R = numpy.linalg.cholesky(A)\n",
    "    return backward(R,forward(R.transpose(),b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEKCAYAAADenhiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xmc1WX5//HXO8ytNDf0myKB/YjCJbVR3EqzVCwVd0ES\nUxRRMXOH1DSXMDW3UBA33IIQUVER1EIBQQRcQVAJFwY1yAVxQ4Hr98d9qHFiYAbOOZ+zvJ+Pxzzm\nfO4553Ou4wPnmnu7bkUEZmZm+fS1rAMwM7PK4+RiZmZ55+RiZmZ55+RiZmZ55+RiZmZ55+RiZmZ5\n5+RiZmZ55+RiZmZ55+RiZmZ5t1rWAWRlo402ilatWmUdhplZWZkyZcq/I6L5ip5XtcmlVatWTJ48\nOeswzMzKiqQ3G/M8D4uZmVneObmYmVneObmYmVneObmYmVneObmYmVneObmYmVneObmYmVneObmY\nmVWLxYuhb194/PGCv5WTi5lZNZg8Gdq3h1NOgXvuKfjbObmYmVWy+fNTQtlxR5gzBwYPhv79C/62\nTi5mZpUoIiWS738fbrgBevaEGTPgiCNAKvjbV21tMTOzivXaa3DyyfDYY/CjH8GDD0JNTVFDcM/F\nzKxSLFwIF10EW28NEyemyfuJE4ueWMA9FzOzyvD3v8NJJ8Grr0KnTnDVVfDtb2cWjnsuZmbl7N13\noUsX+PnP01LjUaNg0KBMEws4uZiZlafFi6FfvzRhP3QoXHABTJ0Ke++ddWRABQ2LSdoCOBf4VkQc\nmnU8ZmYF8+yz0KMHTJoEP/tZWg32ve9lHdVXlETPRdKtkuZKmlqvvYOkVyTNlNRrefeIiFkR0a2w\nkZqZZeijj+C3v4UddoA334S7704rwkossUDp9FwGAn2BO5Y2SGoGXA/sBdQCkyQNB5oBfeq9/tiI\nmFucUM3MiiwC7r0XTj0V3nkHTjwRLr0U1lsv68gaVBLJJSLGSGpVr3lHYGZEzAKQNBjoGBF9gP1W\n5n0kdQe6A7Rs2XKl4zUzK5p//jNtgBw5ErbbDu67L+22L3ElMSzWgM2A2XWua3NtyyRpQ0n9ge0k\n9V7WcyJiQETURERN8+bN8xutmVk+LVyYeidbbQVPPQXXXAPPPFMWiQVKpOeSDxHxHtAj6zjMzFbZ\nE0+koa8ZM+Cww+Dqq2GzBv+2Lkml3HOZA2xe57pFrs3MrDLNnQtdu8JPf5p6LiNGwJAhZZdYoLST\nyySgjaTWklYHOgHDM47JzCz/liyBAQOgbdtUbPLcc2HaNNh336wjW2klkVwkDQImAG0l1UrqFhGL\ngJ7AKGA6MCQipmUZp5lZ3r3wAuy6K5xwAmy7bbq+5BJYa62sI1slJTHnEhGdG2gfAYwocjhmZoW3\nYAFceCFcey1ssAHccQf86ldFKYdfDCWRXMzMqkYE3H8//OY3UFsL3btDnz4pwVSQkhgWMzOrCq+/\nDvvvDwcfnJLJ+PFw440Vl1jAycXMrPC++AIuuwy23DItM/7zn2HKFNh556wjKxgPi5mZFdKYMWnP\nyssvpx7LNdfA5puv+HVlzj0XM7NCmDcPjjkGdt8dPvkkHTV8771VkVjAycXMLL+WLIFbbknnrNx1\nF/TqlXot+61UScSy5WExM7N8eemlNAT21FPw4x+nw7y23DLrqDLhnouZ2ar65BM4+2zYfvtUD+y2\n2+DJJ6s2sYB7LmZmq2b4cDjlFHjrLejWDf70J9hww6yjypx7LmZmK+PNN6Fjx/S17rowbhzcfLMT\nS46Ti5lZU3z5JVxxBbRrB48/Dpdfns6033XXrCMrKR4WMzNrrKeegh49YOpUOOAAuO46+M53so6q\nJLnnYma2Iu+9B8cdB7vtBvPnp9pgDzzgxLIcTi5mZg2JgIED0zkrAwfCWWelPSsdO2YdWcnzsJiZ\n2bK8/HLaszJmDOyyC/TvD1tvnXVUZaNiei6SfiCpv6R7JB2XdTxmVqY+/RR694Yf/jDNrdx8M4wd\n68TSRCWRXCTdKmmupKn12jtIekXSTEm9lnePiJgeET2AI4B9ChmvmVWohx9OGx8vuywd3DVjRtq7\n8rWS+FVZVkrlv9hAoEPdBknNgOuBfYF2QGdJ7SRtLemhel8b515zAOnkysHFDd/MylptLRxySKr/\ntdZaqSz+bbdB8+ZZR1a2SmLOJSLGSGpVr3lHYGZEzAKQNBjoGBF9gGVWgIuI4cBwScOBewsXsZlV\nhEWL0nLiCy6AxYvTiZCnnw6rr551ZGWvJJJLAzYDZte5rgXaN/RkSXsABwNrAk808JzuQHeAli1b\n5ilMMytLo0fDqaemYpO/+AX07QutW2cdVcUo5eTSJBHxBA0klTrPGQAMAKipqYnCR2VmJef11+HM\nM2HYsLRP5d574aCDQMo6sopSKnMuyzIHqHuqTotcm5lZ0338MZx3HvzgBzByJFx8MUyfnk6HdGLJ\nu1LuuUwC2khqTUoqnYAjsw3JzMpOBNx9N5xzDrz9NnTpklaDtWiRdWQVrSR6LpIGAROAtpJqJXWL\niEVAT2AUMB0YEhHTsozTzMrMpEmpoORRR8Gmm6baYHfd5cRSBCXRc4mIzg20jyAtLTYza7x33oHf\n/S6VbNlkk7SsuGtX71cpopJILmZmebFwIVxzDVxySXp89tlw7rnpvBUrKicXMyt/EfDgg2mPyj//\nmcrhX3kltGmTdWRVy31EMytvL78M++yTKhWvvnpaCfbAA04sGXNyMbPy9MEHaRPkNtukiftrr4UX\nXkiJxjLnYTEzKy+LFsFNN8H556cEc8IJcNFFsNFGWUdmdbjnYmblY/Ro2H57OOmkVAL/uefghhuc\nWEqQk4uZlb7XX09Vi/fcEz76CIYOhX/8Iw2JWUlycjGz0tVQyZZDDnHJlhLnORczKz0u2VL23HMx\ns9Liki0VwcnFzErDu+/CMcfAjjvCrFmpZMvEibDLLllHZivBw2Jmlq2FC9MelYsvdsmWCuLkYmbZ\ncMmWiuZhMTMrPpdsqXhOLmZWPC7ZUjU8LGZmheeSLVWnYnoukvaQNFZSf0l7ZB2PmeWMHg0/+pFL\ntlSZkkgukm6VNFfS1HrtHSS9ImmmpF4ruE0AHwNrArWFitXMGun11+HQQ1PJlvnzXbKlypTKsNhA\noC9wx9IGSc2A64G9SMlikqThQDOgT73XHwuMjYgnJW0CXAV0KULcZlbfxx+n3fRXXgnNmqUlxmec\nAWutlXVkVkQlkVwiYoykVvWadwRmRsQsAEmDgY4R0QfYbzm3+wBYoxBxmtlyuGSL1VESyaUBmwGz\n61zXAu0berKkg4F9gPVIvaBlPac70B2gZcuWeQvUrOpNmpRWgU2YADU1cM893llf5Uo5uTRJRAwD\nhq3gOQOAAQA1NTVRjLjMKtq770Lv3jBwIGyySSrZ0rUrfK0kpnMtQ6WcXOYAm9e5bpFrM7OsuWSL\nrUApJ5dJQBtJrUlJpRNwZLYhmVU5l2yxRiqJvqukQcAEoK2kWkndImIR0BMYBUwHhkTEtCzjNKtq\nLtliTVASPZeI6NxA+whgRJHDMbO6PvgALrwQrr8e1lknDYedeCJ8/etZR2YlrMk9F0l7SbpJ0ra5\n6+75D8vMMrdoEfTrl3omfftC9+7w2mvwm984sdgKrUzP5VjgROA8SRsA2+Y3JDPL3OjR8Nvfwosv\nwh57pN6Kd9ZbE6zMnMuCiPgwIs4E9gZ2yHNMZpYVl2yxPFmZ5PLw0gcR0Ys6JVvMrEwtWQJ9+sAP\nfgCPPJKWGE+fDoccAlLW0VkZWuGwmKT6W9mfq9f2wNLriHgrn8GZWRH8+9/wq1/BqFGp13L11S7Z\nYqusMXMut5MqDjf058vSnwWwZ57iMrNiGD8ejjgC5s6F/v3TpL17KpYHK0wuEfHTYgRiZkUUkXoo\n55wDLVummmDbb591VFZBSmKfi5kV0YcfwjHHwP33w0EHwa23wnrrZR2VVZiVmXNpkOdczErclClw\n2GEwezZcdVVabuxhMCuAfMy5LOU5F7NSFQE33pjK4m+8MYwZAzvvnHVUVsE852JW6T7+GE44Af76\nV+jQAe680+fXW8GVROFKMyuQqVNhhx1g8GC45BJ4+GEnFisKT+ibVao77oAePdIZK48/Dj/1IIQV\nj3suZpXms8/guOPg6KOhfXt4/nknFis6JxezSvLqq7DTTnDLLelkyMceg//7v6yjsirUmKXIC0gr\nwf7nR0BEhM81NSsF99wD3bqlg7xGjIB99806IqtijVkttk4xAllVkn4MdCF9pnYRsUvGIZkVx8KF\ncOaZ6cyVnXaCIUNg882zjsqqXJMm9CWtD7QB1lzaFhFjVjUISbcC+wFzI2KrOu0dgGuBZsDNEXFZ\nQ/eIiLHAWEkHApNWNSazsvDGG3D44TBpEpx2Glx2Weq5mGWs0clF0nHAqUAL4HlgJ9K59/nYODkQ\n6Eud8v2SmgHXA3sBtcAkScNJiaZPvdcfGxFzc4+PBLrlISaz0vbQQ9C1KyxeDMOGpVIuZiWiKRP6\np5IOBnszt7FyO+DDfASR6/28X695R2BmRMyKiC+AwUDHiHgpIvar9zUX/lOqZn5ELMhHXGYladGi\nVHBy//2hVSt49lknFis5TUkun0fE5wCS1oiIGUDbwoQFwGbA7DrXtbm25ekG3NbQDyV1lzRZ0uR5\n8+blIUSzIpszJ50Sefnladf9+PHw3e9mHZXZ/2jKnEutpPWA+4HHJH0AvFmYsFZORFywgp8PAAYA\n1NTULGsFnFnpevxxOPJI+PRTuOsu6NIl64jMGtTo5BIRS/vdF0oaDXwLeKQgUSVzgLpLXlrk2syq\ny+LFqXTLH/6QjiEeOjR9NythTZnQ//0ymrcFLspfOF8xCWgjqTUpqXQiTdabVY+5c9MRxI89Bkcd\nBf36wTe+kXVUZivUlDmXT+p8LQb2BVrlIwhJg0grz9pKqpXULSIWAT2BUcB0YEhETMvH+5mVhXHj\nYLvtYOxYuPlmuP12JxYrG00ZFvtz3WtJV5J+8a+yiOjcQPsIYEQ+3sOsbETAlVdC797QunXabf/D\nH2YdlVmTrEpV5LVJ8yBmli/vvw+//jU8+CAcemiqEbauKyxZ+WnKnMtL/LfGWDOgOXBxIYIyq0qT\nJqUjiN9+G667Dnr29BHEVraa0nPZr87jRcC/cvMiZrYqIuD66+H00+Hb305zLO3bZx2V2SppTFXk\n05fzMyLiqvyGZFZFPvoIjj8+FZv85S/TAV8bbJB1VGarrDE9l6VVkduSyr8Mz13vDzxTiKDMqsKL\nL6Z5lVmzUsHJs86Cr/mIJasMjSm5/wcASWOA7ZfW7ZJ0IfBwQaMzq1S33gonnwzrrw//+Af85CdZ\nR2SWV035M2kT4Is611/k2syssT79FI45Jh3qteuu8NxzTixWkZoyoX8H8Iyk+3LXB5JK5ZtZY8yY\nkVaDTZsGF1wA558PzZplHZVZQTRlE+WlkkYCu+WajomI5woTllmFGTw4TdyvuSaMHAl77511RGYF\n1aRNlBExBZhSoFjMKs/ChemEyH790jDY3/4Gm63o5Aiz8rfCORdJ43LfF0j6qM7XAkkfFT5EszI1\na1ZKKP36pZVgo0c7sVjVaMxqsd1y39dZ0XPNLOeBB+Doo9MO+wcegAMOyDois6Jq9GoxSYdJWif3\n+DxJwyRtV7jQzMrQl1/CmWfCgQdCmzbpCGInFqtCTVmKfH5ELJC0G/Bz4Bagf2HCMitDs2fD7rvD\nn/+c9rCMG5eqGptVoaYkl8W5778EBkTEw8Dq+Q/JrAyNGpXOXnnppbQyrG9fWGONrKMyy0xTkssc\nSTeSToQcIWmNJr7erPIsXpz2q+y7L2y6KUyeDEcckXVUZplrSnI4nHQ42N4R8SGwAXBWQaJaCZLa\nSRoiqZ+kQ7OOx6rAu+/CXnul8+2POQaefhrats06KrOS0JTk8hnwDWDpqZFfBz7MRxCSbpU0V9LU\neu0dJL0iaaakXiu4zb7AXyLiRKBrPuIya9CTT6ZhsKefhttuS4d6rb121lGZlYymJJcbgJ34b3JZ\nAFyfpzgGAh3qNkhqlrv/vkA7oHOud7K1pIfqfW0M3Al0knQFsGGe4jL7qiVLoE8f2HNP+Na3YOLE\ndHKkmX1FU3bot4+I7SU9BxARH0jKy4R+RIyR1Kpe847AzIiYBSBpMNAxIvrw1YPL6jo5l5SG5SMu\ns6947z3o2jWdad+pEwwYAOt4+5fZsjQluXyZ+8UdAJKaA0sKElWyGTC7znUt0ODxfLnk9DvS0N0V\nDTynO9AdoGXLlnkK06rCxIlw+OFpnuWGG6BHDx9BbLYcTRkWuw64D9hY0qXAOOCPBYlqJUTEGxHR\nPSK6RMS4Bp4zICJqIqKmefPmxQ7RylEEXHst/PjH6SCvp56CE090YjFbgaZURb5b0hTgZ4CAAyNi\nesEigznA5nWuW+TazIpj/vx07sq990LHjmnifv31s47KrCw0KrlIEtAiImYAMwob0n9MAtpIak1K\nKp2AI4v03lbtJkxI8yuvvw5XXgmnn+7eilkTNGpYLCICGFGoICQNAiYAbSXVSuoWEYuAnqS9NdOB\nIRExrVAxmAHwxhtpsn6XXeDzz9OS4zPOcGIxa6KmTOg/K2mHiJiU7yAionMD7SMoYFIz+4+PPkpL\njK++Os2t/P73qUz+N7+ZdWRmZalJS5GBLpLeBD4hzbtERGxTkMjMimHRorQB8vzzYd48OOoo+OMf\noUWLrCMzK2tNSS77FCwKsyw8+mga8po6Na0GGzECamqyjsqsIjRltdibhQzErGhefjmdufLII7DF\nFmk12EEHeV7FLI9c1diqx7x5cNJJsM02MH58Onfl5Zfh4IOdWMzyrCnDYmbl6fPP4brr4NJL4ZNP\n0ibICy6AjTbKOjKzirXC5CKp0XVSIuKtVQvHLI8iYOhQOOectF9lv/3giivg+9/POjKziteYnsvt\njbxXAHuuQixm+fPMM3DaaWn4a5tt4LHH4Oc/zzoqs6qxwuQSET8tRiBmefHWW9C7N/z1r7DJJnDT\nTekgr2bNso7MrKp4WMwqw4IF8Kc/pUl6gHPPTcNhLolvlgkPi1l5W7w4FZQ87zz417+gS5e0CdJH\nKphlysNiVr4efzxtgnzxxVQL7IEHoH2DR/6YWRF5n4uVnxkzYP/9Ya+9Uk2wIUNg3DgnFrMS4uRi\n5ePf/4ZTToGttoIxY9Icy/TpcNhh3gRpVmK8idJK38KF0LcvXHxxmrg/4QS48ELYeOOsIzOzBji5\nWOmKgPvug7PPhn/+E/bdN22C3HLLrCMzsxXwsJiVpsmTYffd4ZBDYM01YeTIVLXYicWsLJRlcpG0\nhaRbJA1dXpuVodradLzwDjvAK6/AjTfC88/DPj7xwaycFD25SLpV0lxJU+u1d5D0iqSZknot7x4R\nMSsiuq2ozcrIxx+nYpLf+15a/dWrF7z2GnTvDqt59Nas3GTxf+1AoC9wx9IGSc2A64G9gFpgkqTh\nQDOgT73XHxsRc4sTqhXc4sVwxx1pR/0776Tz6/v0gVatso7MzFZB0ZNLRIyR1Kpe847AzIiYBSBp\nMNAxIvoA+xU3Qiua0aPh9NPTsFf79unQrp13zjoqM8uDUplz2QyYXee6Nte2TJI2lNQf2E5S74ba\nlvG67pImS5o8b968PIZvTfLqq9CxI+y5J7z/PgwaBBMmOLGYVZCyHMyOiPeAHitqW8brBgADAGpq\naqJgAdqyvf8+XHQRXH89rLVWGv469dT02MwqSqkklznA5nWuW+TarBJ88QXccENKLPPnw3HHpceb\nbJJ1ZGZWIKUyLDYJaCOptaTVgU7A8IxjslUVAfffn/amnHYa1NSk+ZUbb3RiMatwWSxFHgRMANpK\nqpXULSIWAT2BUcB0YEhETCt2bJZHzz6b5lQOOgi+/vW0AXLUKNh666wjM7MiyGK1WOcG2kcAI4oc\njuXb22+nZcW33w4bbpjmV7xXxazq+P94y49PPoErr4TLL4dFi+DMM+F3v4P11ss6MjPLgJOLrZol\nS+Cuu1IimTMHDj00lcLfYousIzOzDJXKhL6VoyefTDXAjj4aNt0Uxo6Fe+5xYjEzJxdbCTNnwsEH\nwx57wNy5qefy9NOw225ZR2ZmJcLJxRrvgw9SuZZ27eDRR+GSS1Ll4i5d4Gv+p2Rm/+U5F1uxRYug\nX790+uMHH8Cxx6ZTIb/97awjM7MS5eRiy/fGG3Dkkan21557wlVXwQ9/mHVUZlbinFysYffcA8cf\nn1aE3X03dO4MUtZRmVkZ8EC5/a9PP00bHw8/HNq2TSVbjjzSicXMGs3Jxb7qpZdSDbCbboJzzoFx\n47y02MyazMnFkohUuXiHHVJp/EcfhcsuS3XBzMyayMnFUjI55BA4+eS0d+WFF2CvvbKOyszKmJNL\ntRs7Nq3+euihVBtsxAiXwzezVebkUq0WL4Y//CH1VNZYA8aPhzPO8GZIM8sLL0WuRrNnw69+BWPG\npN31N9wA666bdVRmVkGcXKrNAw+kHfYLF6YzV7p2zToiM6tAZTkGImkLSbdIGlqn7QeS+ku6R9Jx\nWcZXkj7/HHr2hAMPhO98J50U6cRiZgWSxTHHt0qaK2lqvfYOkl6RNFNSr+XdIyJmRUS3em3TI6IH\ncASwT/4jL2PTp0P79ulUyNNOS6Vcvve9rKMyswqWRc9lINChboOkZsD1wL5AO6CzpHaStpb0UL2v\njRu6saQDSEclDy5c+GUkAm65JW2KfPttePjhVBtsjTWyjszMKlzR51wiYoykVvWadwRmRsQsAEmD\ngY4R0QfYrwn3Hg4MlzQcuDc/EZep+fPhhBPgb39LBSfvvDMd6GVmVgSlMueyGTC7znVtrm2ZJG0o\nqT+wnaTeubY9JF0naQDwRAOv6y5psqTJ8+bNy1/0pebpp2HbbWHoUPjjH9NueycWMyuislwtFhHv\nAT3qtT1BA0mlznMGAAMAampqokDhZWfJErj8cjjvPGjRIm2Q3HnnrKMysypUKsllDrB5nesWuTZr\nrHfegaOOgr//PVUzvvFGWG+9rKMysypVKsNik4A2klpLWh3oBAzPOKby8cgjqYTL+PGpmvHgwU4s\nZpapLJYiDwImAG0l1UrqFhGLgJ7AKGA6MCQiphU7trKzcGE60/4Xv0hHDk+ZAscd53NXzCxzWawW\n69xA+wjSMmJrjNdeg06d0mbIk09ORSfXXDPrqMzMgNKZc7GmuPNOOOmkdNbKffelXfdmZiWkVOZc\nrDEWLEiT9l27wvbbp3NXnFjMrAQ5uZSLKVNSQvnrX+HCC+Ef/4DNN1/hy8zMsuDkUuqWLIE//znt\nV/n8c3jiCbjgAmjWLOvIzMwa5DmXUjZ3Lhx9NIwcmYa/brkFNtgg66jMzFbIPZdS9fjjae/K6NGp\nmvGwYU4sZlY2nFxKzZdfQq9esPfesP76MGlSWhnmvStmVkY8LFZKXn8dOneGiRPh+OPhmmtg7bWz\njsrMrMmcXErF4MGpRL4EQ4bAYYdlHZGZ2UrzsFjWPvkEunVLPZYtt4Tnn3diMbOy5+SSpeefhx/9\nCG67DX73O3jySWjVKuuozMxWmZNLFiKgb990rv1HH6WVYZdemsq5mJlVAM+5FNt778Gxx8Lw4fDL\nX6ZeS/PmWUdlZpZX7rkU05NPpr0rI0emlWAPPujEYmYVycmlGBYtSiVb9twzLS2eMAFOPdV7V8ys\nYnlYrNDeegu6dIFx41Ipl7594ZvfzDoqM7OCKsuei6QtJN0iaWidtj0kjZXUX9IeGYb3X/fdB9tu\nm0rj33UXDBzoxGJmVSGLY45vlTRX0tR67R0kvSJppqRey7tHRMyKiG71m4GPgTWB2vxG3USffZZK\nthx8MHz3u/Dcc6n3YmZWJbIYFhsI9AXuWNogqRlwPbAXKTFMkjQcaAb0qff6YyNi7jLuOzYinpS0\nCXAVkM1v82nT0vHDU6fCWWfBJZfA6qtnEoqZWVaKnlwiYoykVvWadwRmRsQsAEmDgY4R0QfYr5H3\nXZJ7+AGwRn6ibYIIuOkm+O1vYZ110oqwffYpehhmZqWgVOZcNgNm17muzbUtk6QNJfUHtpPUO9d2\nsKQbgTtJPaNlva67pMmSJs+bNy9/0X/wARx+eKoNtttuaY7FicXMqlhZrhaLiPeAHvXahgHDVvC6\nAcAAgJqamshLMOPHp7pgb78Nl18OZ5wBXyuVnG1mlo1S+S04B6h7IHyLXFvpWrw4lWz5yU9gtdXg\nqafSHIsTi5lZyfRcJgFtJLUmJZVOwJHZhrQcc+bAUUelUyI7d4b+/WHddbOOysysZGSxFHkQMAFo\nK6lWUreIWAT0BEYB04EhETGt2LE1ylNPpRIuEyemumB33+3EYmZWTxarxTo30D4CGFHkcJqudWvY\nfnv4y1+gbdusozEzK0mlMixWPjbdFB59NOsozMxKmmefzcws75xczMws75xczMws75xczMws75xc\nzMws75xczMws75xczMws75xczMws7xSRn+LA5UbSPODNlXz5RsC/8xhOOfBnrg7+zNVhVT7zdyKi\n+YqeVLXJZVVImhwRNVnHUUz+zNXBn7k6FOMze1jMzMzyzsnFzMzyzsll5QzIOoAM+DNXB3/m6lDw\nz+w5FzMzyzv3XMzMLO+cXJpA0uaSRkt6WdI0SadmHVOxSGom6TlJD2UdSzFIWk/SUEkzJE2XtHPW\nMRWapN65f9tTJQ2StGbWMeWbpFslzZU0tU7bBpIek/Ra7vv6WcaYbw185ity/7ZflHSfpPXy/b5O\nLk2zCDgjItoBOwEnS2qXcUzFcirpCOpqcS0wMiK+D/yQCv/skloB3YEfRcRWQDOgU5YxFchAoEO9\ntl7A3yOiDfD33HUlGcj/fubHgK0iYhvgVaB3vt/UyaUJIuKdiHg293gB6RfOZtlGVXiSWgC/BG7O\nOpZikPQt4CfALQAR8UVEfJhtVAX3EfAlsJak1YC1gbezDSn/ImIM8H695o7A7bnHtwMHFjWoAlvW\nZ46IRyNiUe7yaaBFvt/XyWUl5f7S2w6YmG0kRXENcDawJOtAiqQ1MA+4LTcUeLOkb2QdVCFFxPvA\nlcBbwDvyr3s9AAAFL0lEQVTA/IiolvO8N4mId3KP3wU2yTKYDBwLPJLvmzq5rARJ3wTuBX4bER9l\nHU8hSdoPmBsRU7KOpYhWA7YH+kXEdsAnVN5QyVdI+i5wGimxbgp8Q9Kvso2q+CItn62aJbSSziUN\n99+d73s7uTSRpK+TEsvdETEs63iKYFfgAElvAIOBPSXdlW1IBVcL1EbE0l7pUFKyqWQ1wPiImBcR\nXwLDgF0yjqlY/iXp2wC573MzjqcoJP0a2A/oEgXYk+Lk0gSSRBqHnx4RV2UdTzFERO+IaBERrUgT\nvP+IiIr+izYi3gVmS2qba/oZ8HKGIRXDK8BOktbO/Tv/GRW+iKGO4cDRucdHAw9kGEtRSOpAGuo+\nICI+LcR7OLk0za7AUaS/3p/Pff0i66CsIE4B7pb0IrAt8MeM4ymoiHgeuAOYDLxE+t1QcTvXJQ0C\nJgBtJdVK6gZcBuwl6TXg57nritHAZ+4LrAM8lvs91j/v7+sd+mZmlm/uuZiZWd45uZiZWd45uZiZ\nWd45uZiZWd45uZiZWd45uVjFk/RxvetfS+qbe9xDUtfc44GSDs09fkPSRqvwnts2tExd0h6SQtJx\n9Z4fks5cwX0PXF6x1Lqfp5Fx7iRpYm456nRJF+baD5BU0VUJrLBWyzoAsyxFRP7X96fCj9uSdr2P\naOBpU4HD+W8x0M7AC424/YHAQyxjU6ek1Vbi89wOHB4RL0hqBrQFiIjhpM2FZivFPRerapIuXE5v\n4WxJL0l6RtL/yz2/uaR7JU3Kfe1a5z53SnoKuBO4CDgi1yM4Yhn3fhNYU9ImuR3xHahTPFDS8bn7\nv5B7v7Ul7QIcAFyRu+93JT0h6RpJk4FTl34eSavlXr9H7n59JF26jDg2JhWqJCIWR8TLuefX7d09\nX+frM0m7S/pG7pyQZ3LFPTs27b+8VTr3XKwarCXp+TrXG9C4v8rnR8TWuWGma0h1mK4Fro6IcZJa\nAqOAH+Se3w7YLSI+y9VtqomInsu5/1DgMOA54FlgYZ2fDYuImwAkXQJ0i4i/SBoOPBQRQ3M/A1g9\nImpy1xcCRMSiXAxDJZ1CSl7tlxHD1cArkp4ARgK3R8TndZ8QEdvm7r0/qWTIeOAPpFJAxyodNPWM\npMcj4pPlfF6rIk4uVg0+W/oLEv5TsK+mEa8bVOf71bnHPwfa5X6pA6ybq5INMDwiPmtCXEOAvwHf\nz71H3UKRW+WSynrAN0lJrCF/W1ZjREyTdCdpGG3niPhiGc+5SNLdwN7AkaThuT3qP09SG+AK4KcR\n8aWkvUkFTZf2+tYEWlI99chsBZxczBoWy3j8NWCn+n/d55JNk/5qj4h3JX0J7EU66bNuchkIHJib\nC/k1y/iFX8fy3ndr4EPS8FdDcfwT6CfpJmCepA3r/jyXPIcAx9c590TAIRHxynLe26qY51zMGnZE\nne8Tco8fJRW1BNIqrwZeu4BUGHBFfg+cExGL67WvA7yjdMRDl5W4L5IOJg0B/gT4i5ZxTrqkX+q/\n3bA2wGJSMqrrVuC2iBhbp20UcMrS10rarjExWfVwcjFr2Pq5qsinkg7SAvgNUCPpRUkvAz0aeO1o\n0vBZQxP6AETE+Ii4fxk/Op90yulTwIw67YOBs3KT6N9t6L65ZdSXAcdFxKukKrjXLuOpR5HmXJ4n\nLUToUjfRSfoOcChwbJ1J/RrgYuDrwIuSpuWuzf7DVZHNzCzv3HMxM7O8c3IxM7O8c3IxM7O8c3Ix\nM7O8c3IxM7O8c3IxM7O8c3IxM7O8c3IxM7O8+//znkALU+xOAAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8571325400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "norms = []\n",
    "sizes = [2,4,6,8,10,12]\n",
    "for n in sizes:\n",
    "    A=hilbert(n)\n",
    "    x = np.ones(n)\n",
    "    b = A.dot(x)\n",
    "    result = linear_system_cholesky(A,b)\n",
    "    norm = np.linalg.norm(result-x,ord=2)\n",
    "    norms.append(norm)\n",
    "fig,axes = plt.subplots()\n",
    "axes.plot(sizes,norms,color='red')\n",
    "axes.set_ylabel(r\"$\\vert\\vert$ residual $\\vert\\vert_2$\")\n",
    "axes.set_xlabel(r\"Hilbert Matrix Size\")\n",
    "axes.set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f8571231278>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEKCAYAAADenhiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmczfX3wPHXyTYqaUG/pEIhSraxRAtKIVtFWZKtpCgU\nhVIqSyJaiEb2RHaSJdmF7GQnkfFVBtHYzTi/P96XrmkGM+7cz52Z83w85uHe9/3ce8+dpjnz3s5b\nVBVjjDEmkK7yOgBjjDGpjyUXY4wxAWfJxRhjTMBZcjHGGBNwllyMMcYEnCUXY4wxAWfJxRhjTMBZ\ncjHGGBNwllyMMcYEXHqvA/BKtmzZNHfu3F6HYYwxKcqqVasOqGr2S12XZpNL7ty5WblypddhGGNM\niiIiuy/nOhsWM8YYE3CWXIwxxgScJRdjjDEBl2bnXOJz5swZIiMjOXnypNehpEhhYWHkypWLDBky\neB2KMcZjllz8REZGkiVLFnLnzo2IeB1OiqKqHDx4kMjISPLkyeN1OMYYj9mwmJ+TJ09y0003WWJJ\nAhHhpptusl6fMQaw5PIflliSzr53xphzLLkYY0xaceYMe7sO5c/R85L9rSy5hKC//vqL+vXrkzdv\nXkqUKMH999/PpEmTmD9/PlmzZqVo0aLnv3766SevwzXGhDpVYsZO5LOcPbm7c23avXt1sr+lTeiH\nGFWlVq1aNGrUiG+//RaA3bt3M3XqVG644QYefPBBpk2b5nGUxpgUY8kSfnlpCC02tGQtT1ElfD8f\njCmV7G9rPZcQM3fuXDJmzEiLFi3Ot91xxx28+uqrHkZljElxtm3j7+rP83K5ddy/IYL91+dn/Hex\n/LA8B3nvTP75Ueu5JKRNG1i7NrCvWbQofPrpRS/ZuHEjxYsXT/DxRYsWUbRo0fP3J0yYwJ133hmw\nEI0xKdz+/WiX9xk1MJo36M0ByUbrV2L5oMc1ZMkSvDAsuYS4li1bsnjxYjJmzEivXr1sWMwYE7/j\nx6FvX7Z0n8grx3szjwqULn6amV9fRbFiwR+ksuSSkEv0MJLLPffcw4QJE87f79+/PwcOHCA8PNyT\neIwxIS42FoYN40Tn7nTb14SPZRnXXHcVAz+GF1/MyFUeTX7YnEuIqVixIidPnmTAgAHn244fP+5h\nRMaYkKQK06dD0aLMeGE89x5aQDfeoe5zGdiyLR0vvYRniQUsuYQcEWHy5MksWLCAPHnyUKpUKRo1\nakTPnj2Bf+dczn2NHz/e44iNMUG3ejU8+ih7n3iROr9/TFVmkDHPrcydCyNGwM03ex1gKhoWE5G8\nwNtAVlWt7XU8V+KWW25hzJgx8T525MiRIEdjjAkZu3bBO+8QM2oM/a7pQOdMM4iJzUC3btCunZAx\no9cB/iskei4iMkRE9ovIhjjtlUVkq4jsEJEOF3sNVd2pqs2SN1JjjPHA339D+/ZQoAC/jPuDkjf/\nQdtjXXmwYkY2bhQ6dSKkEguETs9lGNAPGHGuQUTSAf2BSkAksEJEpgLpgB5xnt9UVfcHJ1RjjAmS\nU6egf3/o2pW//4ZOBWbw1bYK5EwvjB8PTz0FoVrSLySSi6ouFJHccZpLATtUdSeAiIwBaqpqD6Ba\ncCM0xpggOnsWvvsOOnVCd+1iVOGevCFtObgjA23awPvvE9Q9K0kREsNiCbgV2ON3P9LXFi8RuUlE\nBgLFRKRjAtc0F5GVIrIyKioqsNEaY0wgzJ8PpUtD/fpsCSvKI0UO0vDXN8mTLwMrV0KfPqGfWCC0\nk0uiqOpBVW2hqnf6ejfxXROhquGqGp49e/Zgh2iMMQnbuBGqVYMKFTix7zDv1FjPfb9NZM3uGxk4\nEJYscUU+UopQTi57gdv87ufytRljTOrxv//Biy/CfffB4sVMbzKOezJspdvUwtSrJ2zdiud7VpIi\nlMNdAeQTkTwikhGoC0z1OKag+PPPP6lbty533nknJUqUoGrVqkRERFCtWuKmmsqXL8/KlSsT/f6N\nGze2/TPGJLfoaHj3XciXD4YPJ7JJZ2o/9BdPDK1NprCrmDcPhg+HHDm8DjRpQmJCX0RGA+WBbCIS\nCbynqoNFpBUwC7dCbIiqbvQwzKBQVZ588kkaNWp0fq/LunXrmDo1TeRVY1K/M2fg66+hSxfYv5+Y\nOvX4It/nvPt5NmJi8O1ZCb2lxYkVEj0XVa2nqreoagZVzaWqg33t01U1v28epZvXcQbDvHnzyJAh\nwwUl94sUKcKDDz7I0aNHqV27NnfffTcNGjRAVQGYM2cOxYoVo3DhwjRt2pRTp07953V//PFH7r//\nfooXL06dOnU4evQoAB06dKBQoULcd999tGvX7j/P69y5M40bN2bOnDnUqlXrfPvs2bN58sknA/3x\njUm9VGHyZChcGF55BQoUYNnXGwjf9i2vd8/GQw+5aZdQ3LOSFCHRcwlFHlXcZ8OGDZQoUSLex9as\nWcPGjRvJmTMn5cqV4+effyY8PPz8L//8+fPz/PPPM2DAANq0aXP+eQcOHKBr16789NNPXHPNNfTs\n2ZM+ffrQsmVLJk2axJYtWxARDh8+fMH7tW/fnujoaIYOHQq4Cs1RUVFkz56doUOH0rRp0yv7hhiT\nVixb5jZBLl4Md9/N36Om03FBZSJeFHLmhAkT4MknQ3fPSlKERM/FXJ5SpUqRK1currrqKooWLcqu\nXbvYunUrefLkIX/+/AA0atSIhQsXXvC8ZcuWsWnTJsqVK0fRokUZPnw4u3fvJmvWrISFhdGsWTMm\nTpzI1Vf/e/Tphx9+yJEjRxg4cCAigojQsGFDvvnmGw4fPszSpUupUqVKUD+/MSnOjh1Qpw7cfz9s\n344OGMiINzdQoE0Vvh4stG0LmzeH9mbIpLKeSwI8qrjPPffck+BkeqZMmc7fTpcuHTExMZf1mqpK\npUqVGD169H8eW758OXPmzGH8+PH069ePuXPnAlCyZElWrVrFoUOHuPHGGwFo0qQJ1atXJywsjDp1\n6pA+vf34GBOvAwfggw9gwADIlAm6dGHzE+14pf01zJ8PZcrA7NlQpIjXgSYf67mEmIoVK3Lq1Cki\nIiLOt61fv55FixbFe32BAgXYtWsXO3bsAGDkyJE8/PDDF1xTpkwZfv755/PXHDt2jG3btnH06FGO\nHDlC1apV6du3L+vWrTv/nMqVK9OhQweeeOIJoqOjAciZMyc5c+aka9euNGnSJKCf25hU4cQJ6NED\n7rzTlW1p2pTj67bz9un3KFL2Gtatg4gI+Pnn1J1YwHouIUdEmDRpEm3atKFnz56EhYWRO3fuCybT\n/YWFhTF06FDq1KlDTEwMJUuWvGAxAED27NkZNmwY9erVOz/Z37VrV7JkyULNmjU5efIkqkqfPn0u\neF6dOnWIjo6mRo0aTJ8+ncyZM9OgQQOioqIoWLBg8nwDjEmJYmNh5Ejo3BkiI6F6dejZk+m/F6RV\nJfj9d2jUCD7+OOUuLU4sObfiKK0JDw/XuHtANm/ebL80L6FVq1YUK1aMZs3iL0Bt30OT5syaBW++\nCevXQ8mS0KsXkXc+TJs2bqK+YEE3OhZnQCHFEpFVqnrJo3FtWMxcthIlSrB+/Xqee+45r0Mxxntr\n18Jjj0Hlym5D5JgxxCxeRt/VD1OwIPzwA3Tv7i5LLYklMWxYzFy2VatWeR2CMd774w83/DVyJNxw\nA/TtCy+/zLI1mWhRCtatg6pVoV8/yJPH62C9Yz2XONLqMGEg2PfOpGqHD8Nbb0H+/K4cfvv28Ntv\nHHq+DS+9lomyZeHgQTcUNm1a2k4sYMnlAmFhYRw8eNB+SSaBqnLw4EHCwsK8DsWYwDp9Gj77DO66\nC3r1gmeega1b0Y96MmLq9dx9NwweDK+/nnr3rCSFDYv5yZUrF5GRkdhZL0kTFhZGrly5vA7DmMBQ\nhbFjXT2WnTvhkUdccilWjM2b4eUKsGCB2x85cKAramz+ZcnFT4YMGciT1vuyxhhYuNANey1f7mqB\nzZwJjz3G8RNC107Quzdce63bs9KsWcorhx8M9i0xxphzNm+GmjXd8q69e2HIEFizBh5/nB+mC/fc\n4/ZI1q8PW7e6Y1gsscTPvi3GGPPnn9CiheulzJvn1hBv2wZNmhC5Lx1PP+0Oicyc2Z1CPGwY2GG2\nF2fDYsaYtOvoUfjkEzeXcuqUK4XfuTNkz05MDHzeB957z23A79HDTdqnhnL4wWDJxRiTNi1Z4sa3\ndu+Gp5922SNfPgCWLnUdmfXrXY/l889taXFipZphMREpKCIDRWSciLzgdTzGmBAVGwsffggPPeQm\nTBYuhPHjIV8+Dh2C5s2hbFk4dAgmTYKpUy2xJEVIJBcRGSIi+0VkQ5z2yiKyVUR2iEiHi72Gqm5W\n1RbAs8DjyRmvMSaF2rMHKlZ0Z9c/84ybrH/wQVTdefUFCrg5/Hbt3Nx+rVq2ZyWpQiK5AMOAyv4N\nIpIO6A9UAQoB9USkkIgUFpFpcb5y+J5TA5gOjAlu+MaYkDdxoqtzv2qVm5EfNQqyZmXTJqhQARo3\ndqNiq1e7KZhrr/U64JQtJJKLqi4EDsVpLgXsUNWdqnoalzBqquqvqlotztd+3+tMVdXKQKPgfgJj\nTMg6ftxNoDz9NOTN63orjRoRfVTo0MHlm/XrYdAgdwqxbYYMjFCe0L8V2ON3PxIondDFIlIeeAoI\nA+YncE1zoDnA7bffHqAwjTEha/16qFcPNm1ymyK7duVs+oyMHA4dOrgVyI0bu3NWbGlxYIVyckkU\nVZ1PAknF75oIIALceS7JH5UxxhOq7iTIdu3g+uvdmSuPPcbSpdC6NaxYAaVLw5QpUKqU18GmTiEx\nLJaAvcBtfvdz+dqMMSZhBw64Xfavvurqga1fz957HqNhQ7cKbO9eVy1/yRJLLMkplJPLCiCfiOQR\nkYxAXWCqxzEZY0LZnDlu0mTWLPj0U06Mm0bXiBzkzw/jxsHbb7uyLc89Z2VbkltIfHtFZDSwFCgg\nIpEi0kxVY4BWwCxgMzBWVTd6GacxJkSdOQMdO0KlSnDddeiyXxh/a2sKFhI6d4YqVdzS4q5dbRVY\nsITEnIuq1kugfTpuabExxsTvt9/cTvvly+GFF1jX7HNat83MggWuEzN3rltqbIIrJHouxhiTJKNG\nQbFisHUrUV9PpkW6QRQvl5kNG2DAALelxRKLN0Ki52KMMYkSHQ0tW8LIkZy5/yH6PzqJLm/cyLFj\n8NprbgP+DTd4HWTaZsnFGJOyrFjh9q78/jsz64+g7ern2PKh8Pjj0LcvFCzodYAGbFjMGJNSnD3r\ndjuWLcu2Y7dSrfR+qnzbkNhYYdo0mDHDEksosZ6LMSb07dsHzz/PkZ+W82G+CXz2e3WuPi707u22\ns9gZK6HHei7GmND2ww/EFi7KoAX5yZflT/rsqE7jxsK2bfDGG5ZYQpX1XIwxoenkSXjrLRZ+vobW\nYQtYe+ZuHigNMz+D4sW9Ds5civVcjDGhZ/Nmdhd/kmc/L8vDLORgtgKMGePO9bLEkjJYz8UYEzpU\nOdZ/GB+33cfHMRORjBnp0gnatxeuvtrr4ExiWHIxxoQEPfQ3Y6oM583lTxPJbdSteYKen6fDTsdI\nmWxYzBjjuVWD1/DgLTuov7wNOW7NwKIFZxk9ObMllhTMkosxxjN/RsbQrNhqSr5QhO2xeRj89k6W\n7/4/HnjIfjWldPZf0BgTdKdOQa9Of5M/9ylGrr2XdoVmsH13Jpp2zUu6dF5HZwLB5lyMMUGjCt9/\nD683P8pvf91A9XTT+aTXGfK1q+l1aCbArOdijAmKjRvh8UdjqVkTMv71BzPzv8bUrQUssaRS1nMx\nxiSrQ4egSxf48kslix7lM97l5fbXkqHbJ5Ahg9fhmWSSanouIlJeRBaJyEARKe91PMakdTEx8OWX\nkC+f0r/fWZprBNuzl+O1n2qQ4eNullhSuZBILiIyRET2i8iGOO2VRWSriOwQkQ6XeBkFjgJhQGRy\nxWqMubQ5c9wZXi1bQhFZzxotypdVp5Ftw3x45BGvwzNBEBLJBRgGVPZvEJF0QH+gClAIqCcihUSk\nsIhMi/OVA1ikqlWAt4D3gxy/MQbYuROeegoefRSORp1gwvXNmBNdmvu+aA5Tp0K2bF6HaIIkJOZc\nVHWhiOSO01wK2KGqOwFEZAxQU1V7ANUu8nJ/A5nie0BEmgPNAW633VnGBEx0NHTvDn36QIYMSvcH\nZ9J20ZOEFcwLY5a7w+xNmhIqPZf43Ars8bsf6WuLl4g8JSJfASOBfvFdo6oRqhququHZs2cPaLDG\npEVnz8Lw4ZA/P3z0EdR94h+23VmVjouqEvZSY1i50hJLGhUSPZdAUNWJwESv4zAmrVi2zJ1Xv2IF\nlCoFk1+aQelPnnET9RMmuPExk2aFcs9lL3Cb3/1cvjZjjIf27oWGDeH++yEyEkZ8dYKldzWk9PtV\nXT38desssZiQ7rmsAPKJSB5cUqkL1Pc2JGPSrhMn3JxK9+5umXGnTtDx0RVc+0Jd2L0bPvjANVr9\nFkOI9FxEZDSwFCggIpEi0kxVY4BWwCxgMzBWVTd6GacxaZGqG+UqVAjeeQcqV4bNG8/SLctHXPtY\nWZdpFiyAzp0tsZjzEt1zEZFKwDNAf1VdKyLNVTXiSoJQ1XoJtE8Hpl/Jaxtjkm7dOmjTBubPh3vv\ndftXKt79PzcuNncuPPMMfPUVXH+916GaEJOUnktToD3wnIhUBIoGNiRjjNeiouDll90Uyvr1bqf9\nmjVQ8dj3bvXXsmUweDCMGWOJxcQrKcklWlUPq2o74DGgZIBjMsZ45MwZ+Owzt7R40CBo1Qq2b4eX\nm5wkfdtXoUYNuO02WL0amjYFEa9DNiEqKcnlh3M3VLUDMCJw4RhjvDJ/vuuUtGkDJUu6Hstnn8GN\nf25ya4379YO2bV2vpUABr8M1Ie6Scy4iEncr+5o4bVPO3VfVPwIZnDEm+cXGQrdurnJx3ryuSku1\naiAoDPzKJZQsWWD6dKhSxetwTQpxORP6w3FFIRPq/557TIGKAYrLGBMEUVHw3HPw449Qv76bm7/2\nWlyd/BdegEmT4LHH3Db8//s/r8M1Kcglk4uqVghGIMaY4Fq8GOrWhQMHXFJ58UXfFMqCBS7j/PUX\n9O7tei5XhcSuBZOC2E+MMWnM2bPQqxeULw9hYbB0KTRvDhIbA+++CxUrQubM7oE33rDEYpIkKXMu\nCbI5F2NC26FD0LixO8f+6afdauKsWXHjY3XquF5L48bwxRe+8TFjkiYQcy7n2JyLMSFs+XK35/F/\n/3OrwF591TcMtno1PPkk7N8PI0e6ITFjrpDNuRiTyqm6VcRvvAG33AKLFkHp0r4HR4+GZs3cIV6L\nF0OJEp7GalIPG0w1JhU7csT1Vl57zS36WrPGl1hiY+Gtt9wSsfBwd+6KJRYTQJZcjEml1q51eWPS\nJOjZ0+1fufFG4O+/3UaWjz92NV5++gly5PA6XJPKhHLJfWNMEqjC11+7OZWbboJ58+DBB30PbtoE\nNWu6EvkREW79sTHJwHouxqQiR4/C88+7pcUPPuiGwc4nlilT3JhYdLTLOJZYTDK6ZHIRkWgR+See\nr2gR+ScYQRpjLm2TrwTYqFGulMvMmb7RrrNn3UFetWpBwYJufqVcOa/DNanc5awWyxKMQIwxSffN\nN/DSS3DNNa6Uy6OP+h6IjoZGjdzEy/PPu634YWGexmrShkTNuYjIDUA+4PxPp6ouDHRQSSEiDwIN\ncJ+pkKqW9TgkY5LdiRPQurUrj//QQ25lcc6cvgd/+83Nr2zZAn37ugutRL4JkstOLiLyAtAayAWs\nBcrgjia+4o2TIjIEqAbsV9V7/dorA58B6YCvVfWjhF5DVRcBi0SkFrDiSmMyJtRt3+421a9bBx07\nupGv9Of+j549G5591iWTWbPgkUc8jdWkPYmZ0G+NOxhst29jZTHgcIDiGAZU9m8QkXRAf6AKUAio\nJyKFRKSwiEyL8+W/jrI+8G2A4jImJI0b57al/PEHTJsG3bv7EosqfPKJO+g+Vy5YscISi/FEYobF\nTqrqSRFBRDKp6hYRCciJQaq6UERyx2kuBexQ1Z0AIjIGqKmqPXC9nP/w1UE7oqrRCTzeHGgOcPvt\nl10yzZiQceoUtG/vSn+VLg3ffQd33OF78MQJt0zsm29c4bBhw6w+mPFMYnoukSJyPTAZmC0iU4Dd\nyRMWALcCe/zf39d2Mc2AoQk9qKoRqhququHZs2cPQIjGBM+uXW5Z8RdfuNMiFy70Syx79rgHR42C\nrl1d18YSi/HQZfdcVPVJ380uIjIPyArMSJaokkhV3/M6BmOSw/ffu8VeZ8/ChAnw1FN+Dy5e7Hoq\nJ064vSzVq3sWpzHnJGZC/914mosCHwQunAvsBW7zu5/L12ZMmnHmDLz9tjt/pVgx1yG5806/C776\nym3Fz50b5s93+1iMCQGJGRY75vcVi5toz50MMZ2zAsgnInlEJCNQF5iajO9nTEiJjIQKFVxiadEC\nlizxSyynT7vGFi3cppblyy2xmJCSmGGxT/zvi0hvYFYgghCR0UB5IJuIRALvqepgEWnle490wBBV\n3RiI9zMm1P34IzRo4Ea6Ro1yxYvP++svqF3bDYd16ODmWNKl8yxWY+JzJYUrr8YNVV0xVa2XQPt0\nYHog3sOYlCA2Ft5/3+WLQoVg/Hi4+26/C1audAd7HTwIY8a4vSzGhKDEzLn8ijttElxPIjvwYXIE\nZUxa9Ndfrocyd66r2NK/vyvnct4337hikzff7MbIihb1LFZjLiUxPRf/vSUxwF+qGhPgeIxJk+bP\nh3r14PBhGDIEmjTxezAmxh3s1acPlC8PY8eCLaU3Ie6SyUVEXr/IY6hqn8CGZEzacfYsfPQRdO4M\nd93lKrXcd5/fBYcOQd26rpzLq6+63fcZMngWrzGX63J6LueqIhfAlX85t2KrOrA8OYIyJi04eBAa\nNoQZM9zUyaBBkMW/BvmGDa7wZGRkPN0ZY0Lb5ZTcfx9ARBYCxc+VVhGRLsAPyRqdManU0qUuofz1\nl5tbefnlOAWLJ050uyavuw4WLIAyZTyL1ZikSMw+l5uB0373T/vajDGXSdVVv3/oIVdocskSeOUV\nv8Ry9iy8+67bcX/vvW51mCUWkwIlZkJ/BLBcRCb57tfCVTM2xlyGw4fdyNbkye5QyKFD4frr/S74\n5x83TjZ1qrvwyy/tYC+TYiVmE2U3EZkJPOBraqKqa5InLGNSl1Wr3Nkre/a4Ofm2beMMg23f7uZX\ntm1zlSlbtrSDvUyKlqhNlKq6CliVTLEYk+qowsCBropxjhxu+qRs3DNSZ850K8LSp3erwipU8CRW\nYwLpknMuIrLY92+0iPzj9xUtIv8kf4jGpEzR0W5T5CuvQMWKsGZNnMSiCj17QtWqrvDkypWWWEyq\ncTmrxR7w/ZvlUtcaY5xff3Xlv3bscKdEvvUWXOX/p9zx49Cs2b8lXAYPjrMd35iU7bJXi4lIHRHJ\n4rv9johMFJFiyReaMSnT0KHulMh//oE5c9z59hcklt27oVw5d4xkjx4werQlFpPqJGYpcmdVjRaR\nB4BHgcHAwOQJy5iU5/hxt8iraVO3enjNGlet5QILFkB4OPz+O0yb5qoa28S9SYUSk1xiff8+AUSo\n6g9AxsCHZEzKs3Wr660MH+5KucyeDf/3f34XqLrdko8+CtmyufNXqlb1LF5jkltiksteEfkKd2jX\ndBHJlMjnG5MqjR7tOiN//ulKuXzwQZzjVU6dctWMW7WCKlVg2TLIn9+zeI0JhsQkh2dwB3c9pqqH\ngRuB9skSVRKISCERGSsiA0SkttfxmNTv5ElXtqV+fShSxA2DPf54nIv27XMrwAYPhnfecTsos2b1\nJF5jgikxyeUEcA1w7mCvDMDhQAQhIkNEZL+IbIjTXllEtorIDhHpcImXqQJ8oaovA88HIi5jErJz\np5uTHzgQ2reHefMgV9yj8375xXVp1q+HcePgww/jzOwbk3ol5if9S6AM/yaXaKB/gOIYBlT2bxCR\ndL7XrwIUAur5eieFRWRanK8cwEigroj0Am4KUFzG/MekSVC8uEswU6bAxx/HUwV/2DBXQCxTJlel\nsrZ1pk3akpgd+qVVtbiIrAFQ1b9FJCAT+qq6UERyx2kuBexQ1Z0AIjIGqKmqPbjw4DJ/LX1JaWIg\n4jLG3+nTbnFX376uQzJunNv7eIEzZ6BdO/j8c3jkEbfc+Cb7W8ekPYlJLmd8v7gVQESyA2eTJSrn\nVmCP3/1IoHRCF/uSUyfc0F2vBK5pDjQHuP322wMUpkkL/vjD7XVctsyd2dWrl+uUXODAAXjmGTdG\n1rat69KkT1SFJWNSjcT85H8OTAJyiEg3oDbwTrJElQSqugtf4rjINRFABEB4eLgGISyTCkyf7ooV\nnznjThiuUyeei9atc6WO9+1z65Gft2k/k7YlpiryKBFZBTwCCFBLVTcnW2SwF7jN734uX5sxQRET\n445W6dHDrQYbNw7y5YvnwrFj3e7JG26ARYugZMmgx2pMqLms5CIiAuRS1S3AluQN6bwVQD4RyYNL\nKnWB+kF6b5PG7dnjeisLFsALL7gplMyZ41wUG+t2TPbo4SpSTpgQZ+ekMWnXZa0WU1UFpidXECIy\nGlgKFBCRSBFppqoxQCvc3prNwFhV3ZhcMRgDcPSo660UKAArVrgRrkGD4kksR46481d69HAbJOfO\ntcRijJ/EzLmsFpGSqroi0EGoar0E2qeTjEnNmHNiY10iefttt9O+bl2XN/6zGgxgyxaXWHbudKdF\ntmhh9cGMiSNRS5GBBiKyGziGm3dRVb0vWSIzJkjmzIE33nBz8mXKwMSJcP/9CVw8bRo0aOCWis2Z\n4/ayGGP+IzHJJW5hC2NStK1b3e767793PZQxY9xK4ng7IaquK/POO1CsmNtJacvZjUlQYlaL7U7O\nQIwJloMH4f33YcAAN5fSsye89hqEhSXwhKNH3Wqw8eNdIbFBg+Dqq4MaszEpje3wMmnGqVPQrx90\n7eoO8nrpJejSxZ1tn6Dff3f7VzZscDsn33jD5leMuQyWXEyqp+rmUd58083BV6ni8sQ991ziiXPn\nunGy2FiA4ihCAAAV6klEQVS3k/I/JY+NMQm5ZHIRkcseWFbVP64sHGMCa8UKeP11WLzYJZOZMy8j\nRxw+7MbKevVya5KnTIG77gpKvMakFpfTcxl+ma+lQMUriMWYgNmzBzp1gm++ccNeX33ljh++aKmv\n48fdbsmePV2CadjQnR6ZJUvQ4jYmtbhkclHVCsEIxJhAOHrU5Ybevd1wWKdO8NZbcN11F3nS6dPw\n9dfuvJU//4QnnoBu3VzNF2NMktiwmEkVYmPdESrvvOPyQ/360L073HHHJZ40ejS8956bjHnwQVdA\n7IEHghW2MamWDYuZFO+nn9wirvXrXYmvyZOhdIKHM+C6NN9/77bjb9gARYu6CfvKlW0lmDEBYsNi\nJsXavNltgvzhB7cJcuxYd+DjRfPD/PlurGzpUlfieMwYV0Pfjh82JqDs/yiT4hw4AK1aQeHCrsL9\nxx+7RFOnzkUSy6pVbplYhQru5K+ICNi40Z0AZonFmICzfS4mxTh1Cr74wm2CPHr0302Q2bNf5Elb\ntriy+OPHu+OGe/eGV16Jp8yxMSaQLLmYkKfqjkp58023Yb5qVbcFpVChizzpjz9cjZdhw1yplnff\ndRMzF102ZowJFEsuJqQtX+42Qf78sxsG+/FHqFTpIk+IinLLxL780t1/7TXo2PESNV6MMYFmycWE\npD/+cPPuo0bBzTe7WpFNmkC6dAk84Z9/4JNPoE8ftxmycWO3xNgqFxvjiRSZXEQkL/A2kFVVayfU\nZlKe6Gj46COXI8CtFn7rrYtskj950vVSund35Y5r13abIe++O2gxG2P+K+jLZERkiIjsF5ENcdor\ni8hWEdkhIh0u9hqqulNVm12qzaQcsbGud5Ivn8sTTz/tzlvp2jWBxBIT43bV58vn5lJKlHCFxMaN\ns8RiTAjwoucyDOgHjDjXICLpgP5AJSASWCEiU4F0QI84z2+qqvuDE6oJhtmzXX749VcoVw6mToVS\npRK4+OxZt/LrnXdg+3Z3dOTIkVC+fDBDNsZcQtCTi6ouFJHccZpLATtUdSeAiIwBaqpqD6BacCM0\nwbJpk9sEOX065MnjOh1PP32RkyBnzXITMWvWuBLHkydDjRq2q96YEBQqu8duBfb43Y/0tcVLRG4S\nkYFAMRHpmFBbPM9rLiIrRWRlVFRUAMM3iREVBS1bwn33uVVgvXu7TZAJ7q5fssT1TKpUcdWKR4xw\nB97XrGmJxZgQlSIn9FX1INDiUm3xPC8CiAAIDw/XZAvQxOvkyX83QR47Bi+/7BZ0ZcuWwBPWr3cz\n+tOmuSVj/frBiy9CxoxBjdsYk3ihklz2Arf53c/lazOpgKqbJnnrLbcJslo1V7KlYMEEnvDbb27T\n4+jRkDWrm+F/7TW45pqgxm2MSbpQSS4rgHwikgeXVOoC9b0NyQTCL7+4TZBLlrhhsNmz4dFHE7j4\nf/9zy4i//hoyZHDZ6M034YYbghqzMebKebEUeTSwFCggIpEi0kxVY4BWwCxgMzBWVTcGOzYTOLt3\nuzNVypRxR6V8/TWsXp1AYjl0yCWSu+5yFzZv7novPXpYYjEmhfJitVi9BNqnA9ODHI4JsH/++XcT\npIhbMfzmmwnsVTl6FD77zBUK++cfaNDA1QPLmzfocRtjAitUhsVMChcTA0OGuALE+/fDc8+5qZLb\nbovn4lOnXMn7rl3dxTVquNuFCwc9bmNM8rDkYq7Yjz+6TZAbNriTgqdNg5Il47kwNha++cYtEdu9\n2y0vnjwZ7r8/2CEbY5JZqOxzMSnQpk2u/P3jj7takRMmwIIF8SQWVZg0yc3oN27s1h7PmgVz51pi\nMSaVsuRiEm3/frdH5b773CqwTz5xieapp+LZ0zhnjpvVf+op13MZN87VAHvsMdsAaUwqZsnFXLaT\nJ6FnT7eoa9Agd6Djjh1uqXGmTHEuXr7cLQ179FHYtw8GD3bjZpc85N4YkxrYnIu5JFUYO9atFt69\nG6pXd5sg4y0+vGmTWyI2aZIb/urbF1q0gLCwoMdtjPGO9VzMRUVFuWRSty5cfz389JOrWvyfxLJr\nl5tPKVzYXfT++26DS5s2lliMSYOs52ISNHeuW1J88CB8+im0ahXPSZB//QXdusHAgXDVVdC2LXTo\ncJGCYcaYtMCSi/mPM2egSxe3QT5/fpgxA4oUiXPRkSNu8+Onn7rJmCZNXD2weDe2GGPSGksu5gK7\ndkG9erBsGTRr5jbQX1Av8vhxV534o4/g77/h2Wfhgw9cFjLGGB9LLua8sWNdWS9VGDPG5Y3zzpxx\nK74++MCt/qpc2W3BL1bMs3iNMaHLJvQNx465Y1KefdaVwV+71i+xqLrdkQULus0tefK4nZIzZlhi\nMcYkyJJLGrduHYSHu05Jp06wcKHLH4DLMhUquL0pmTPD99/D4sXw0EOexmyMCX2WXNIoVTd1Urq0\nm5ufPdst+sqQAbcFv3lzKF7cbXwcMMCdW1+tmm2ANMZcFptzSYMOHICmTV1H5IknYOhQyJ4dOH0a\nPv/cHdh1/Di0bu1WgNmZKsaYRLLkksbMn++OTTlwwK0ifu01EBSmfu9KG+/Y4apRfvJJAlvwjTHm\n0lLksJiI5BWRwSIy3q+toIgMFJFxIvKCl/GFopgYd9ZKxYpw7bVuqXHr1iAbN7gikjVrQvr0bqL+\nhx8ssRhjrogXxxwPEZH9IrIhTntlEdkqIjtEpMPFXkNVd6pqszhtm1W1BfAs8HjgI0+5du+Ghx92\n53E1bgyrVkGx2w9Cy5Zud+TKlW5Dy/r1bomxMcZcIS+GxYYB/YAR5xpEJB3QH6gERAIrRGQqkA7o\nEef5TVV1f3wvLCI1gFeAQYEPO2UaPx5eeAHOnoVvv4V6tc/Al1+6LfjR0W558fvvw003eR2qMSYV\nCXpyUdWFIpI7TnMpYIeq7gQQkTFATVXtAVRLxGtPBab6EtOEwEScMh0/7sp8RURAqVIwejTk3TYT\n7msLW7ZApUquYvE993gdqjEmFQqVOZdbgT1+9yN9bfESkZtEZCBQTEQ6+trKi8jnIhIBzE/gec1F\nZKWIrIyKigpc9CFm/Xq3dyUiwpXJXzx4K3lffQKqVHGTL1OnupMgLbEYY5JJilwtpqoHgRZx2uaT\nQFLxuyYCiAAIDw/XZArPM6puxOuNN9zq4R8nRFNp0btQrB9cfTX07g2vvgoZM3odqjEmlQuV5LIX\n8C+nm8vXZi7TwYOu0OSUKVCl8lmGVRhBjubt4NAhN+nStSvkyOF1mMaYNCJUhsVWAPlEJI+IZATq\nAlM9jinFWLAAihaF6dOhT4ttTNtTlBxvNYF774XVq934mCUWY0wQebEUeTSwFCggIpEi0kxVY4BW\nwCxgMzBWVTcGO7aUJiYG3nvP7V3JnP40y8q+TtuBBbjqWLRbJjZvnss6xhgTZF6sFquXQPt0YHqQ\nw0mx/vjD7bRfvBga3bOSL7Y9TpaoU64Mftu2drSwMcZToTIsZhJhwgQoUkRZu+I031z3CsM2liRL\ngxqwfTt07GiJxRjjOUsuKciJE9CihauAn+/0RtaeKkiDe9fBihWu+uQtt3gdojHGAKGzWsxcwoYN\nUPfp02zclpH2fEzXGwaScXB3d6qXlcE3xoQY67mEOFUY8OlJShY7w4Fth5iVsTofdzlBxm0boG5d\nSyzGmJBkPZcQdujAWV6ovIdJq+7gcWYy/Mkp3PzZl3DbbZd+sjHGeMh6LiFq0YANFL11P9NW3ULv\n2z5l+uKs3DxxgCUWY0yKYD2XEBOzK5JuNX7hg19rkTfdbpZ0WUV459fgKvs7wBiTclhyCRXHj7On\ncwQNPg1n0dmnee7etXz5411kuSWv15EZY0yiWXLxmip89x2TWs2h2cGenEmfmRG9omj4uu2sN8ak\nXDbW4qWVKzlR9hFeqXeIpw4OIm+BDKzelJmGr2f3OjJjjLki1nPxwr590KkTG4etoG66cWygIG+0\nPUv3j7JYNXxjTKpgPZdgOnkSevRA8+UnYmQYJdOvZv+NBZgxA3r3ucoSizEm1bCeSzCowsSJ0L49\nf//+Ny/eMosJx8pSqRKMGAH/939eB2iMMYFlPZfktnatq4lfuzaLeYAi2fcxJaosH38MM2daYjHG\npE6WXJLL/v3w0ktQvDix6zfyYbVfeHj3cDJkCePnn6F9e9u6YoxJvezXW6CdPg2ffAL58sGQIUQ2\nfZdHCu7l3WmlqFtXWLMGSpXyOkhjjEleKTK5iEheERksIuP92sqLyCIRGSgi5YMelCp8/707Wrhd\nO3jgAaZ8+jtFJnVh5doMDBsG33wD110X9MiMMSbovDjmeIiI7BeRDXHaK4vIVhHZISIdLvYaqrpT\nVZvFbQaOAmFAZGCjvoSNG+Hxx6FGDUiXjpOTZ/Jq3h+o1SoXd9zhjrFv1MgKGBtj0g4vVosNA/oB\nI841iEg6oD9QCZcYVojIVCAd0CPO85uq6v54XneRqi4QkZuBPkCDZIj9QgcPukPsBw6ELFng00/Z\nXOEVnn0uA7/+6k4b7tEDMmVK9kiMMSakBD25qOpCEckdp7kUsENVdwKIyBigpqr2AKpd5uue9d38\nG0jeX+dnzsCAAdClCxw5Ai1aoF3eZ/CUbLxWBq69Fn74AapWTdYojDEmZIXKnMutwB6/+5G+tniJ\nyE0iMhAoJiIdfW1PichXwEhczyi+5zUXkZUisjIqKippka5cCUWKQOvWUKIErFvH4W79ebZlNl58\nEcqWhXXrLLEYY9K2FLmJUlUPAi3itE0EJl7ieRFABEB4eLgm6c1vvNGtIZ4yBapXZ8lSoX5R2LsX\nPvrIlhgbYwyETnLZC/ifgpXL1xZ68uaFX38l9qzwUXc35XL77bBoEZQp43VwxhgTGkIluawA8olI\nHlxSqQvU9zakhO39n9CwIcyb546xHzgQsmb1OipjjAkdXixFHg0sBQqISKSINFPVGKAVMAvYDIxV\n1Y3Bju1yLF7splx++QWGDIFvv7XEYowxcXmxWqxeAu3TgelBDifR8uaF4sXhiy+gQAGvozHGmNAU\nKsNiKUbOnPDjj15HYYwxoc3WNRljjAk4Sy7GGGMCzpKLMcaYgLPkYowxJuAsuRhjjAk4Sy7GGGMC\nzpKLMcaYgLPkYowxJuBENWnFgVM6EYkCdifx6dmAAwEMJyWwz5w22GdOG67kM9+hqtkvdVGaTS5X\nQkRWqmq413EEk33mtME+c9oQjM9sw2LGGGMCzpKLMcaYgLPkkjQRXgfgAfvMaYN95rQh2T+zzbkY\nY4wJOOu5GGOMCThLLokgIreJyDwR2SQiG0WktdcxBYuIpBORNSIyzetYgkFErheR8SKyRUQ2i8j9\nXseU3ESko+9ne4OIjBaRMK9jCjQRGSIi+0Vkg1/bjSIyW0S2+/69wcsYAy2Bz9zL97O9XkQmicj1\ngX5fSy6JEwO8oaqFgDJASxEp5HFMwdIadwR1WvEZMFNV7waKkMo/u4jkBpoDJVT1XiAdUNfLmJLJ\nMKBynLYOwBxVzQfM8d1PTYbx3888G7hXVe8DtgEdA/2mllwSQVX3qepq3+1o3C+cW72NKvmJSC7g\nCeBrr2MJBhHJCjwEDAZQ1dOqetjbqJLdP8AZILOIpAeuBv7nbUiBp6oLgUNxmmsCw323hwO1ghpU\nMovvM6vqj6oa47u7DMgV6Pe15JJEvr/0igG/eBtJUHwKvAmc9TqQIMkDRAFDfUOBX4vINV4HlZxU\n9RDQG/gD2AccUdW0cqD3zaq6z3f7T+BmL4PxQFNgRqBf1JJLEojItcAEoI2q/uN1PMlJRKoB+1V1\nldexBFF6oDgwQFWLAcdIfUMlFxCRO4G2uMSaE7hGRJ7zNqrgU7d8Ns0soRWRt3HD/aMC/dqWXBJJ\nRDLgEssoVZ3odTxBUA6oISK7gDFARRH5xtuQkl0kEKmq53ql43HJJjULB5aoapSqngEmAmU9jilY\n/hKRWwB8/+73OJ6gEJHGQDWggSbDnhRLLokgIoIbh9+sqn28jicYVLWjquZS1dy4Cd65qpqq/6JV\n1T+BPSJSwNf0CLDJw5CCYStQRkSu9v2cP0IqX8TgZyrQyHe7ETDFw1iCQkQq44a6a6jq8eR4D0su\niVMOaIj7632t76uq10GZZPEqMEpE1gNFge4ex5OsVHUtMAJYCfyK+92Q6naui8hoYClQQEQiRaQZ\n8BFQSUS2A4/67qcaCXzmfkAWYLbv99jAgL+v7dA3xhgTaNZzMcYYE3CWXIwxxgScJRdjjDEBZ8nF\nGGNMwFlyMcYYE3CWXEyqJyJH49xvLCL9fLdbiMjzvtvDRKS27/YuEcl2Be9ZNKFl6iJSXkRURF6I\nc72KSLtLvG6tixVL9f88lxlnGRH5xbccdbOIdPG11xCRVF2VwCSv9F4HYIyXVDXw6/td4ceiuF3v\n0xO4bAPwDP8WA60HrLuMl68FTCOeTZ0ikj4Jn2c48IyqrhORdEABAFWdittcaEySWM/FpGki0uUi\nvYU3ReRXEVkuInf5rs8uIhNEZIXvq5zf64wUkZ+BkcAHwLO+HsGz8bz2biBMRG727YivjF/xQBF5\n0ff663zvd7WIlAVqAL18r3uniMwXkU9FZCXQ+tznEZH0vueX971eDxHpFk8cOXCFKlHVWFXd5Lve\nv3e31u/rhIg8LCLX+M4JWe4r7lkzcd95k9pZz8WkBZlFZK3f/Ru5vL/Kj6hqYd8w06e4OkyfAX1V\ndbGI3A7MAgr6ri8EPKCqJ3x1m8JVtdVFXn88UAdYA6wGTvk9NlFVBwGISFegmap+ISJTgWmqOt73\nGEBGVQ333e8CoKoxvhjGi8iruORVOp4Y+gJbRWQ+MBMYrqon/S9Q1aK+166OKxmyBHgfVwqoqbiD\nppaLyE+qeuwin9ekIZZcTFpw4twvSDhfsC/8Mp432u/fvr7bjwKFfL/UAa7zVckGmKqqJxIR11jg\nO+Bu33v4F4q815dUrgeuxSWxhHwXX6OqbhSRkbhhtPtV9XQ813wgIqOAx4D6uOG58nGvE5F8QC+g\ngqqeEZHHcAVNz/X6woDbSTv1yMwlWHIxJmEaz+2rgDJx/7r3JZtE/dWuqn+KyBmgEu6kT//kMgyo\n5ZsLaUw8v/D9XOx9CwOHccNfCcXxGzBARAYBUSJyk//jvuQ5FnjR79wTAZ5W1a0XeW+ThtmcizEJ\ne9bv36W+2z/iiloCbpVXAs+NxhUGvJR3gbdUNTZOexZgn7gjHhok4XURkadwQ4APAV9IPOeki8gT\n8m83LB8Qi0tG/oYAQ1V1kV/bLODVc88VkWKXE5NJOyy5GJOwG3xVkVvjDtICeA0IF5H1IrIJaJHA\nc+fhhs8SmtAHQFWXqOrkeB7qjDvl9Gdgi1/7GKC9bxL9zoRe17eM+iPgBVXdhquC+1k8lzbEzbms\nxS1EaOCf6ETkDqA20NRvUj8c+BDIAKwXkY2++8acZ1WRjTHGBJz1XIwxxgScJRdjjDEBZ8nFGGNM\nwFlyMcYYE3CWXIwxxgScJRdjjDEBZ8nFGGNMwFlyMcYYE3D/D5/3fL2Y6cbKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f857139b828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "norms = []\n",
    "sizes = [2,4,6,8,10,12]\n",
    "for n in sizes:\n",
    "    A=hilbert(n)\n",
    "    x = np.ones(n)\n",
    "    b = A.dot(x)\n",
    "    result = linear_system_lu(A,b)\n",
    "    norm = np.linalg.norm(result-x,ord=2)\n",
    "    norms.append(norm)\n",
    "fig,axes = plt.subplots()\n",
    "axes.plot(sizes,norms,color='red',label='GE')\n",
    "axes.set_ylabel(r\"$\\vert\\vert$ residual $\\vert\\vert_2$\")\n",
    "axes.set_xlabel(r\"Hilbert Matrix Size\")\n",
    "axes.set_yscale('log')\n",
    "\n",
    "norms = []\n",
    "for n in sizes:\n",
    "    A=hilbert(n)\n",
    "    x = np.ones(n)\n",
    "    b = A.dot(x)\n",
    "    result = linear_system_cholesky(A,b)\n",
    "    norm = np.linalg.norm(result-x,ord=2)\n",
    "    norms.append(norm)\n",
    "axes.plot(sizes,norms,color='blue',label='Cholesky')\n",
    "axes.set_ylabel(r\"$\\vert\\vert$ residual $\\vert\\vert_2$\")\n",
    "axes.set_xlabel(r\"Hilbert Matrix Size\")\n",
    "axes.set_yscale('log')\n",
    "axes.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see above, the two-norm of the residual between the true answer and the calculated result is comparable cholesky factorization and gaussian elimination indicating that both methods are comparably stable, and that the conditioning of the problem will be the dominating factor in accuracy for the tested LU decomposition methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probelm 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now compare with QR decompositions from previous work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def classic_gram_schmidt(A):\n",
    "    \"\"\" The jth column of Q is the jth column \n",
    "    of A with the projected subtractions of \n",
    "    every preceeding column in A, then normalized \"\"\"\n",
    "    \n",
    "    nrows,ncols = A.shape\n",
    "    \n",
    "    Q = np.copy(A,order='F') # columns in contiguous memory\n",
    "    R = np.zeros(shape=(ncols,ncols))\n",
    "    for colj in range(ncols):\n",
    "        for coli in range(colj):\n",
    "            R[coli,colj] = Q[:,coli].dot(Q[:,colj])\n",
    "            Q[:,colj] -= R[coli,colj]*Q[:,coli]\n",
    "        norm2 = np.linalg.norm(Q[:,colj],ord=2)\n",
    "        if  norm2:\n",
    "            Q[:,colj] /= norm2\n",
    "        else:\n",
    "            raise RuntimeError(\"One of the vectors of A is the 0 vector, or A is not of full rank\")\n",
    "            \n",
    "        R[colj,colj] = norm2\n",
    "        \n",
    "    return Q,R\n",
    "def modified_gram_schmidt(A):\n",
    "    \"\"\" Subtract the ith normalized column of A from every jth column of A such that j>i \"\"\"\n",
    "    nrows,ncols = A.shape\n",
    "    \n",
    "    Q = np.copy(A,order='F') # columns in contiguous memory\n",
    "    R = np.zeros(shape=(ncols,ncols))\n",
    "    \n",
    "    for coli in range(ncols):\n",
    "        R[coli,coli] = np.linalg.norm(Q[:,coli],ord=2)\n",
    "        if R[coli,coli]:\n",
    "            Q[:,coli] /= R[coli,coli]\n",
    "        else: \n",
    "            raise RuntimeError(\"One of the columns of A is the 0 vector\")\n",
    "        for colj in range(coli+1,ncols):\n",
    "            R[coli,colj] = Q[:,coli].dot(Q[:,colj])\n",
    "            Q[:,colj] -= R[coli,colj]*Q[:,coli]     \n",
    "    return Q,R\n",
    "def householder_reflector(x):\n",
    "    v = np.copy(x)\n",
    "    v[0] += np.sign(x[0])*np.linalg.norm(x,ord=2)\n",
    "    return v\n",
    "\n",
    "def householder_qr(A):\n",
    "    R = np.copy(A,order='F')\n",
    "    nrows,ncols = R.shape\n",
    "    \n",
    "    hvectors = []\n",
    "    for icol in range(ncols):\n",
    "        #print(R[icol:,icol])\n",
    "        v = householder_reflector(R[icol:,icol])\n",
    "        norm2 = np.linalg.norm(v,ord=2)\n",
    "        if norm2:\n",
    "            v /= norm2\n",
    "        else:\n",
    "            raise RuntimeError(\"Norm of zero vector\")\n",
    "        \n",
    "        R[icol:,icol:] -= 2*np.outer(v,v.dot(R[icol:,icol:]))\n",
    "        hvectors.append(v)\n",
    "        \n",
    "    # calculate Q according to 10.3\n",
    "    Q = np.eye(nrows)\n",
    "    for icol in range(ncols):\n",
    "        for k in reversed(range(ncols)):    \n",
    "            Q[k:,icol] -= 2*hvectors[k]*(hvectors[k].dot(Q[k:,icol]))\n",
    "\n",
    "    return Q[:nrows,:ncols], R[:ncols,:ncols]\n",
    "\n",
    "def back_substitution(R,b):\n",
    "    nrows,ncols = R.shape\n",
    "    if nrows != ncols: raise RuntimeError(\"This function only supports square upper triangular matrices.\")\n",
    "    x = np.zeros(ncols)\n",
    "    for i,row in enumerate(np.flip(R,axis=0)):\n",
    "        pos = ncols-(i+1)\n",
    "        x[pos] = (b[pos] - row[pos+1:].dot(x[pos+1:]))/row[pos]\n",
    "    return x\n",
    "\n",
    "def least_squares_via_qr(A,b,qr_decomp_method):\n",
    "    # Step 1, compute the reduced QR factorization A = QR\n",
    "    Q,R = qr_decomp_method(A)\n",
    "    # Step 2, Compute the vector [Q*]b\n",
    "    rhs = Q.transpose().dot(b)\n",
    "    # Step 3, Solve Rx = [Q*]b\n",
    "    return back_substitution(R,rhs)\n",
    "\n",
    "def least_squares_given_qr(Q,R,b,qr_decomp_method):\n",
    "    # Step 1, Compute the vector [Q*]b\n",
    "    rhs = Q.transpose().dot(b)\n",
    "    # Step 2, Solve Rx = [Q*]b\n",
    "    return back_substitution(R,rhs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEKCAYAAADenhiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlcVOX+wPHPwyYiIC4sKipIIoIsCiKCKKLgrrmVpqZW\nem23W/6y7WZdrbzVtTJvZotaWZalpmkqKoR7Qq7grqioIIsgyCLMPL8/BidQUVBgBnjerxcvOWfO\nnPOdAefLc55zvl8hpURRFEVRqpKJoQNQFEVR6h6VXBRFUZQqp5KLoiiKUuVUclEURVGqnEouiqIo\nSpVTyUVRFEWpciq5KIqiKFVOJRdFURSlyqnkoiiKolQ5M0MHYCjNmzeXLi4uhg5DURSlVomPj0+X\nUtrfbbt6m1xcXFyIi4szdBiKoii1ihDibEW2U6fFFEVRlCqnkouiKIpS5VRyURRFUapcvZ1zuZ2i\noiKSk5MpKCgwdCiKUoalpSXOzs6Ym5sbOhRFqRCVXEpJTk7GxsYGFxcXhBCGDkdRAJBSkpGRQXJy\nMq6uroYOR1EqRJ0WK6WgoIBmzZqpxKIYFSEEzZo1UyNqpVZRyeUmKrEoxkj9Xiq1jUouiqIo9YRW\nW8zps+tr5FgquRih1NRUHnnkEdq1a4e/vz/du3dn1apVxMTE0LhxY/z8/PRfmzdvNnS4iqLUAqfP\n/s7h75vSJHYwmVknqv14akLfyEgpefDBB5k4cSLff/89AGfPnmXNmjU0adKE0NBQfvvtNwNHqShK\nbVFcXMD2jQ8SdGUj+VKQ0OYJQmzdqv24tXLkIoRoJ4T4Sgjxc6l1YUKIbUKIhUKIMAOGd1+2bt2K\nhYUF06ZN069r27Ytzz77rAGjUhSlNnPMjGWfaSuKBuynR69FCJPq/+g3mpGLEOJrYDBwWUrZqdT6\n/sDHgCnwpZTyPSnlaeDx0skFkEAuYAkk33dA06fD/v33vZsy/Pzgo4/uuElCQgJdunQp9/Ft27bh\n5+enX/7ll19wc6v+v0IURak9CguvsitqJL6hC2nS2I1Ww4/T0dq5RmMwppHLEqB/6RVCCFNgATAA\n8ATGCiE8y3n+NinlAOBl4K1qjLNGPf300/j6+tK1a1cAQkND2b9/v/5LJRZFUUrbk7yHJ5d6E3Z1\nM4fidB+FtjWcWMCIRi5SylghhMtNqwOBkyUjFYQQy4FhQOJtnq8t+fYK0OB2xxBCTAWmArRp0+bO\nAd1lhFFdvLy8+OWXX/TLCxYsID09nYCAAIPEoyhK7ZCXn853m//Bk/tW09KmJdu6zaen7zMGi8eY\nRi630wo4X2o5GWglhGgmhFgIdBZCvAIghBghhPgc+Bb49HY7k1IuklIGSCkD7O3v2o7AIMLDwyko\nKOCzzz7Tr8vLyzNgRIqiGL3UGHJWt2Ni9kpm+I0j4akEQg2YWMCIRi6VIaXMAKbdtG4lsNIwEVUd\nIQSrV6/mhRde4D//+Q/29vY0atSIuXPnArfOubz++uuMGjXKUOEqimJAV3OTMT/4Gg2TvqGJVRsS\nXd/mPd/phg4LMP7kcgFoXWrZuWRdndaiRQuWL19+28eys7NrOBpFUYzR3r1v0erov3E01YDHi1j4\nvI2fmZWhw9Iz9uSyF2gvhHBFl1TGAI8YNiRFURTDycw6wZGN/QjRnOEkDTji9zmdPB83dFi3MJo5\nFyHED8AuoIMQIlkI8biUshh4BtgIHAF+klImGDJORVEUQzl46H9o1nYgsPgMMdahtH7oslEmFjCi\nkYuUcmw569cDNVMMR1EUxYjZOwRyXtiRGbiQsPYPGTqcOzKa5KIoiqKUJaVk+x/TMLu0jqCHz9HC\nMYAW4zMNHVaFGM1pMUVRFOVv57LPMfD7gcQeXkQjzVWycs4YOqRKUSMXRVEUI6LVFrN9ywQ+OLya\nbfmmDA6fR6fAZzAxqV0f12rkYoRSUlIYM2YMbm5u+Pv7M3DgQBYtWsTgwYMrtZ+wsDDi4uIqffxJ\nkybx888/331DRVGq1tUTaKJ60TNtOU/ZN+HwU4d5Omh6rUssoEYuRkdKyfDhw5k4caL+XpcDBw6w\nZs0aA0emKEp10RRf56+YCQRkrMHcpAFpnd6jX6cZNVK9uLrU3sjrqOjoaMzNzcuU3Pf19SU0NJTc\n3FxGjRqFh4cH48aNQ0oJwJYtW+jcuTPe3t489thjFBYW3rLfTZs20b17d7p06cLo0aPJzc0FYObM\nmXh6euLj48NLL710y/PeeOMNJk2axJYtW3jwwQf166Oiohg+fHhVv3xFqXdOnF7F0R+a0vXyT1xq\n5A2DErH3eblWJxZQI5dyGajiPocPH8bf3/+2j+3bt4+EhARatmxJSEgIO3bsICAgQP/h7+7uzqOP\nPspnn33G9Ol/l4BIT09n9uzZbN68WV9K5r///S9PP/00q1at4ujRowghyMrKKnO8GTNmkJOTw+LF\niwFdhea0tDTs7e1ZvHgxjz322P29IYpSj12/nsvOjUMJvhrNVQQ72zxH9+B5UMuTyg1141XUE4GB\ngTg7O2NiYoKfnx9JSUkcO3YMV1dX3N3dAZg4cSKxsbFlnrd7924SExMJCQnBz8+PpUuXcvbsWRo3\nboylpSWPP/44K1euxMrq79IR//73v8nOzmbhwoUIIRBCMGHCBL777juysrLYtWsXAwYMqNHXryh1\nyf6f3QjLiSbOtC1icCLBPT6u9aOV0tTIpRwGqriPl5dXuZPpDRr83UnA1NSU4uLiCu1TSklERAQ/\n/PDDLY/9+eefbNmyhZ9//plPP/2UrVu3AtC1a1fi4+PJzMykadOmAEyePJkhQ4ZgaWnJ6NGjMTNT\nvz6KUhn5BZmYCDMaNLDF3PP/+LMwi+Bu/zZ0WNWi7qTJOiI8PJzCwkIWLVqkX3fw4EG2bdt22+07\ndOhAUlISJ0+eBODbb7+lV69eZbYJCgpix44d+m2uXbvG8ePHyc3NJTs7m4EDBzJv3jwOHDigf07/\n/v2ZOXMmgwYNIicnB4CWLVvSsmVLZs+ezeTJk6v0dStKXbfn+CpSVjixa8MgADr7vUhgHU0soEYu\nRkcIwapVq5g+fTpz587F0tISFxeXMpPppVlaWrJ48WJGjx5NcXExXbt2LXMxAIC9vT1Llixh7Nix\n+sn+2bNnY2Njw7BhwygoKEBKyX//+98yzxs9ejQ5OTkMHTqU9evX07BhQ8aNG0daWhodO3asnjdA\nUeqYnIKrvLL1Vf63dwE/tbamXZshhg6pRogbVxzVNwEBAfLme0COHDmiPjTv4plnnqFz5848/rhx\nFsury9TvZ+1zaP+HmB54hcjzRYwKeJ454XNoZNHI0GHdFyFEvJTyrq1x1chFqTB/f38aNWrEhx9+\naOhQFMW4Xb8Cf72I9+nFnDaxYO3IJXTuONHQUdUolVyUCouPjzd0CIpi9PbsfgXPpE+xkfng+Qpt\nvV7D1Lx2j1buhUouiqIoVSAtM4ETmwYQrD3PMU1DXPtvx8I+CFNDB2YgtfJqMSFEOyHEV0KIn++0\nTlEUpbpJrZYd257GZJ03/przxNj2pd2YdCzsgwwdmkEZTXIRQnwthLgshDh80/r+QohjQoiTQoiZ\nAFLK01LKMjPKt1unKIpS3c5djCHw3P+4KBpxPvhXwgZHYW5uPL3sDcVokguwBOhfeoUQwhRYAAwA\nPIGxQgjPmg9NURTlb1KrZd9+3aX7bZ3DOer3OZ5jMnjAdaiBIzMeRpNcpJSxwM0t1gKBkyWjkuvA\ncmBYjQenKIpS4lTmKb5c1p7OiS9yKOFzALy9pmJqZmHgyIyL0SSXcrQCzpdaTgZaCSGaCSEWAp2F\nEK8A3G7dzYQQU4UQcUKIuLS0tGoPvqrMmjWLDz74oMr2FxwcXGX7Mjb32sPmhn/9619s3rz5lvUx\nMTH6fjoxMTHs3Lnzno+h1E6a4ut8vv1tvD/zZs7Fy/zh9CidOk4xdFhGq1ZeLSalzACm3W3dbZ63\nCFgEupsoqy1AI1ddH4zFxcW1vt7Y22+/fddtYmJisLa2rtNJWilLZiVybF13/Auv0td1AP8bvAhn\nW2dDh2XUjH3kcgFoXWrZuWRdjQhbEnbL1wc7P7jnxyvqm2++wcfHB19fXyZMmFDmsS+++IKuXbvi\n6+vLyJEjycvLA2DFihV06tQJX19fevbsCUBCQgKBgYH4+fnh4+PDiRMnALC2ttbvb+7cuXh7e+Pr\n68vMmTPLjWnv3r34+Pjg5+fHjBkz6NSpEwBLlixh6NChhIeH06dPH3Jzc+nTpw9dunTB29ubX3/9\nFYCkpCQ8PDyYNGkS7u7uPPLII2zatIng4GDat2/Pn3/+We6x//jjD/z8/PDz86Nz5876Wmflxb5i\nxQoCAwNxd3fX12RbsmQJDz74IBEREbi4uDB//nw++OADOnfuTFBQEJmZujOypbtwbtiwAQ8PD7p0\n6cLKlSv1r2PhwoXMmzcPPz+/cmu+KXVDcXE+8vBsxIbOtDPVUPjAk/w65jeVWCrA2P/M3Au0F0K4\noksqY4BHDBtS9UpISGD27Nns3LmT5s2bk5mZySeffKJ/fMSIEUyZohuKv/7663z11Vc8++yzvP32\n22zcuJFWrVrp+7IsXLiQ559/nnHjxnH9+nU0Gk2ZY/3+++/8+uuv7NmzBysrK/0H7O1MnjyZL774\ngu7du9+ShP766y8OHjxI06ZNKS4uZtWqVdja2pKenk5QUBBDh+omOU+ePMmKFSv4+uuv6dq1Kz/8\n8AM7duxgzZo1vPPOO6xevfq2x/7ggw9YsGABISEh5ObmYmlpecfYi4uL+fPPP1m/fj1vvfWW/jTX\n4cOH2bdvHwUFBbi5ufGf//yHffv28cILL/DNN9+U6YFTUFDAlClT2Lp1Kw888AAPP/wwAC4uLkyb\nNg1ra+vbNldT6o6jx3+APY/jYZoPbUZj6T+fkIaOhg6r1jCa5CKE+AEIA5oLIZKBN6WUXwkhngE2\nAqbA11LKhJqKKWZSTLU+fjtbt25l9OjRNG/eHEBf7v6Gw4cP8/rrr5OVlUVubi79+vUDICQkhEmT\nJvHQQw8xYsQIALp3786cOXNITk5mxIgRtG/fvsy+Nm/ezOTJk/V9XG4+1g1ZWVnk5OTQvXt3AB55\n5BF+++03/eMRERH650opefXVV4mNjcXExIQLFy6QmpoKgKurK97e3oCutUDfvn0RQuDt7U1SUlK5\n70lISAj//Oc/GTduHCNGjMDZ2fmOsd94/f7+/mX227t3b2xsbLCxscHOzo4hQ3QFBL29vTl48GCZ\nYx49ehRXV1f9ezZ+/PgylaqVuqugMIvdGwbRI3cnGZgQ5/YqAd3mGDqsWsdoTotJKcdKKVtIKc2l\nlM5Syq9K1q+XUrpLKd2klPX+Jzxp0iQ+/fRTDh06xJtvvklBQQGgG6XMnj2b8+fP4+/vT0ZGBo88\n8ghr1qyhYcOGDBw4UN+rpao1avR3aYtly5aRlpZGfHw8+/fvx9HRUR9j6X40JiYm+mUTE5M79qaZ\nOXMmX375Jfn5+YSEhHD06NE7xnNjvzf3vLnX4yv1x6GERVz4yZGwazvZZf4ADYaeVInlHhlNclF0\nwsPDWbFiBRkZGQC3nKrKycmhRYsWFBUVsWzZMv36U6dO0a1bN95++23s7e05f/48p0+fpl27djz3\n3HMMGzbslr/OIyIiWLx4sX7eprzTYnZ2dtjY2LBnzx4Ali9fXm782dnZODg4YG5uTnR0NGfPnq38\nm3CTU6dO4e3tzcsvv0zXrl05evRohWO/Vx4eHiQlJXHq1CmAMo3WbGxs9PM+St1yJWU7DdAQ5z6b\n0IdOYNfY1dAh1VpGc1pM0fHy8uK1116jV69emJqa0rlzZ1xcXPSP//vf/6Zbt27Y29vTrVs3/Yfc\njBkzOHHiBFJK+vTpg6+vL3PnzuXbb7/F3NwcJycnXn311TLH6t+/P/v37ycgIAALCwsGDhzIO++8\nc9u4vvrqK6ZMmYKJiQm9evWicePGt91u3LhxDBkyBG9vbwICAvDw8Ljv9+Sjjz4iOjoaExMTvLy8\nGDBgAA0aNKhw7PfC0tKSRYsWMWjQIKysrAgNDdW/10OGDGHUqFH8+uuvzJ8/n9DQ0Co7rlLz9u3/\nkPxryQSHzCO09xKuFfwHZysnQ4dV66l+LqWofhnly83N1V9l9t5773Hp0iU+/vhjA0dVv6jfz6p1\nKecSr22ZyYtXvgGzRniOu1qnethXF9XPRalS69at491336W4uJi2bduyZMkSQ4ekKPckLz+dXVse\nYcKhHaQXF+HRdSrPhr2rEksVU8lFKePpp59mx44dZdY9//zzTJ48WX85bnVZvHjxLaOhkJAQFixY\nUK3HVeoJqYWk7xFxL9CnKJ2Zrv4MivwRt6Zuho6sTlKnxUpRpx0UY6Z+P+/dgYOf0vLkh9gXJCGb\n+nOi7TTcOz5h6LBqJXVaTFGUeu9scjQXt02gu7xAisaM4uCvMWs3EXehToFVN5VcFEWpc7KvnmXf\nllF0z4ujuYSYxuEE9v0Rs4bNDR1avaGSi6Iodc7Zs+vomRfHDvP2tO+9nDD7LoYOqd5RyUVRlFpP\nSklc3NtcS9tN2MDf8fF+ivNNOxHaqqehQ6u31InHOszFxYX09HSgbA+XGTNm4OXlxYwZM1i4cCHf\nfPNNpfZbuqpybVW6P8u9uHjxIqNGjbrtY6V7ylTljZ3K7R1KPUS/7/qRsm8WzhnRFBZeBaC1SiwG\npUYu9UTpHi6LFi0iMzMTU1NTg8Wj0WgMevz71bJlS31p/jt55513bqmMoFSNtPRDJEY/xJSTx0g3\nsSMldA79Ap7DwqL2//FTF6iRy51sDrv715EPym5/83LprwqoSN+TzMxMHnzwQXx8fAgKCtLXDMvI\nyCAyMhIvLy+eeOIJSl9mfmO0MXToUHJzc/H39+fHH38s0+Xy1KlT9O/fH39/f0JDQ/UFIs+cOUP3\n7t3x9vbm9ddfv2P8Wq2Wp556Cg8PDyIiIhg4cKD+Q9jFxYWXX36ZLl26sGLFinJ700yaNIknn3yS\noKAg2rVrR3R0NBMnTqRjx45MmjSp3GNrNBomTZpEp06d8Pb2Zt68eYCu1H/fvn3x9fWlS5cu+nph\nubm5jBo1Cg8PD8aNG6d/v1xcXHjllVfw8/PD39+f+Ph4IiMjcXNzY+HChfqf042eNvn5+YwZM4aO\nHTsyfPhw8vPzAV3Bzfz8fPz8/Bg3blxFfvxKRRTnQ8K7NN0SRPD1o8zyjODkcyeZEvyqSizGREpZ\nL7/8/f3lzRITE8uuiOp196/E98tuf/Ny6a8KOHPmjDQ1NZUHDx6UGo1GdunSRU6aNElqtVq5evVq\nOWzYMPnMM8/IWbNmSSml3LJli/T19ZVSSvnss8/Kt956S0op5W+//SYBmZaWJqWUslGjRvpjlP7+\nzTfflO+/r4s5PDxcHj9+XEop5e7du2Xv3r2llFIOGTJELl26VEop5aefflrm+TdbsWKFHDBggNRo\nNPLSpUvSzs5OrlixQkopZdu2beXcuXP126anp+u/f+211+Qnn3wipZRy4sSJ8uGHH9a/Zmtr6zLv\nx759+2577Li4ONm3b1/98pUrV6SUUgYGBsqVK1dKKaXMz8+X165dk9HR0dLW1laeP39eajQaGRQU\nJLdt26aP83//+5+UUsrp06fLTp06yatXr8rLly9LBwcH/c/Jy8tLSinlhx9+KCdPniyllPLAgQPS\n1NRU7t2795b3+n7d8vtZz2g1Grkj9hmZ+YONlMuQ8o/hMiN1t6HDqneAOFmBz1h1WuxO+sbc3/aV\nfX6Ju/U9OXv2LL/88gugq6KckZHB1atXiY2N1XdMHDRoEE2aNKnwMXNzc9m5cyejR4/WryssLARg\nx44d+uNNmDCBl19+udz9bN++ndGjR2NiYoKTkxO9e/cu83jpu/zL600DuuKQN16zk5NTmfcjKSkJ\nPz+/W47drl07Tp8+zbPPPsugQYOIjIwkJyeHCxcuMHz4cEBXkPKGwMBAnJ11HQX9/PxISkqiR48e\nAPoGZ97e3uTm5ur7wDRo0EDfjO2G2NhYnnvuOQB8fHzw8fEp9/1R7s2hhM8Rf71IsOk1jmoaYhW+\ngQat+nH7DkSKMagzyUUIEQqMQ/eaPKWUtbbB+d36jpibm1f5MbVaLXZ2duzfv/+2jwshquQ4pXu/\nTJo0idWrV+Pr68uSJUuIiYnRP1b6Nd/8fpTXe6VJkyYcOHCAjRs3snDhQn766ac7Ftcsvd/yer9U\n5vhK1buUupcz0SMJ1p7nEiZsb/kEwT0/w8Skznx01VlGPecihPhaCHFZCHH4pvX9hRDHhBAnhRAz\nAaSU26SU04DfgKWGiLemhIaG6nu5xMTE0Lx5c2xtbenZsyfff/89oGthfOXKlQrv09bWFldXV1as\nWAHoTpceOHAA0NX3utHDpXQPmdsJCQnhl19+QavVkpqaWiZh3Ky83jT3Kj09Ha1Wy8iRI5k9ezZ/\n/fUXNjY2ODs761soFxYW6ud2qkrp9/3w4cNl+uaYm5tTVFRUpcerTzSaQjyKkomxCcN25CV6hH2h\nEkstYdTJBVgC9C+9QghhCiwABgCewFghhGepTR4Bvq+pAA1h1qxZxMfH4+Pjw8yZM1m6VJdL33zz\nTWJjY/Hy8mLlypW0adOmUvtdtmwZX331Fb6+vnh5efHrr78C8PHHH7NgwQK8vb25cOHCHfcxcuRI\nnJ2d8fT0ZPz48XTp0qXc3i83etOEhIRUSd+XCxcuEBYWhp+fH+PHj+fdd98F4Ntvv+WTTz7Bx8eH\n4OBgUlJS7vtYpT355JPk5ubSsWNH/vWvf+Hv769/bOrUqfj4+KgJ/Qoq1hYTs3k8239sB4Bzyx5Y\njrpM2JBoGlk5GDg6pTKMvnClEMIF+E1K2alkuTswS0rZr2T5FQAp5btCiDbAG1LKKXfbrypcWX1u\n9H7JyMggMDCQHTt24OSkmi/dr7r++/n7id95KeolHpeJhDdujNuDidhYtzR0WMpN6nLhylbA+VLL\nyUC3ku8fBxaX90QhxFRgKlDpv+qVihs8eDBZWVlcv36dN954QyUW5Y4unN/IpdhxvH8hg0JLN9r1\n+RHfjqNUf5VarjYml3JJKd+8y+OLgEWgG7nUSFB11KFDh5gwYUKZdQ0aNGDPnj13nGepKt26ddNf\nzXbDt99+q7+qTDF+Mi8FcXgWLU9+gbWElzqPp2/4V1iYWhg6NKUK1MbkcgFoXWrZuWSdUoO8vb3L\nvbKsJuzZs8dgx1buT0FhFrs3j6VrVhRWJgLh/gyNvF5lYENHQ4emVKHamFz2Au2FEK7oksoYdJP4\niqIYManVsmvnC7Q5s4AwUw17hBOefddiYx9QKz+IlDsz6p+pEOIHIAxoLoRIBt6UUn4lhHgG2AiY\nAl9LKRMMGKaiKHdx/kIsWdGDCDbJ5RiW7POcSze/Fw0dllKNjDq5SCnHlrN+PbC+hsNRFKWSNMXX\nMTWzoEkTd3Klhm0tJhEc+jmmZmpepa5Tl2MYmdIFEWtSZcvoly54WZqh4leMy9XCq/y+JpyTPzSm\nqCgPaysnPMblEtp7sUosBiSl5FByNkdTrlb7sVRyUQzqfkupqFIsxqW4uIAv4hbSfn57fj0dTZpF\nS/Ly0wDUpcUGUqzRsvNkOrPWJBDy3laGfLqdz2JOVftx1U/bCGk0GqZMmYKXlxeRkZHk5+ezf/9+\ngoKC8PHxYfjw4frSLqUbU6Wnp+Pi4gJAQkICgYGB+Pn54ePjw4kTJwD47rvv9Ov/8Y9/oNFo9Md9\n7bXX8PX1JSgoiNTUVEA3EgkPD8fHx4c+ffpw7ty5W+KNj4/H19cXX19fFixYUOZ1zJgxg65du+Lj\n48Pnn38O6ErWhIaGMnToUDw9PW/Z3w1z5szB3d2dHj16MHbsWP1IKSwsjOnTpxMQEHDH2mFKzUo5\n/g1nfrDjz21P4t7MnSce2kuPh0/R2LatoUOrd/Kva9iYkMI/f9pPwJzNPPLlHn748xxerRrz/igf\n3hziVe0xqORyB0vClrAkbAnpx3TdHHd+sJMlYUvY+YGu8Vb6sXT9NjesnbqWJWFLOLb2GADH1h5j\nSdgS1k5dW+HjnjhxgqeffpqEhATs7Oz45ZdfePTRR5k7dy4HDx7E29ubt9566477WLhwIc8//zz7\n9+8nLi4OZ2dnjhw5wo8//siOHTvYv38/pqam+ppe165dIygoiAMHDtCzZ0+++OILAJ599lkmTpzI\nwYMHGTdunL76b2mTJ09m/vz5+lpkN3z11Vc0btyYvXv3snfvXr744gvOnDkDwF9//cXHH3/M8ePH\nbxt/fHw8y5cvZ//+/axfv569e/eWefz69evExcXx4otqUtjQ8jPiIXogTnETaWRiwqSgV4idFEtA\ny7vexK1UoSvXrvNzfDJTv4mj87838Y9v49ly5DLhHRxYON6fff+K4ItHAxgd0Jqmjar/1KRRT+jX\nV66urvqS8v7+/pw6dYqsrCx69eoFwMSJE8uUxr+d7t27M2fOHJKTkxkxYgTt27dny5YtxMfH07Vr\nV0DX5MrBQVevycLCQt/219/fn6ioKAB27dqlL+M/YcIE/u///q/McbKyssjKyqJnz576bX7//XcA\nNm3axMGDB/XNwrKzszlx4gQWFhYEBgbi6upabvzbtm1j+PDhWFlZAX+XwL+hdOl+xTDSM4+QuHU0\nwYUJSIvGiM4f0NL9GVqaNrj7k5UqcSErn6iEFDYmpPJnUiYaraRFY0seDmhNpJcTga5NMTc1zBhC\nJZc7mBQzqcxy8EvBBL/0dyX/5h2a37LNkEVDyix3GNKBDkM6VOq4N5eCv7l/SGlmZmZotVoACgoK\n9OsfeeQRunXrxrp16xg4cCCff/45UkomTpyoL+hYmrm5ub6s/s3l5++VlJL58+eX6dMCutNipUvv\n34v7fb5y7woLr7Jry1g6Z64nWMCOBt50jliNbeN2hg6tzpNScuJyLhsPp7ApMZVDF7IBaO9gzZO9\n3Ij0csS7VeMqa5FxP9RpsVqgcePGNGnShG3btgG6Mic3RjEuLi7Ex8cDlOnpfvr0adq1a8dzzz3H\nsGHDOHghxxi9AAAgAElEQVTwIH369OHnn3/m8uXLAGRmZnL27Nk7Hjs4OLhMuf3Q0NAyj9vZ2WFn\nZ8f27dv129zQr18/PvvsM33J+ePHj3Pt2rUKveaePXuyevVq8vPzycnJYe3aip9WVKrPhZRdpP7U\njLCs9RwzceBcj3X0GnVQJZZqpNVK4s9m8u76I/T+IIbIebF8GHUcM1PBzAEebH2xF1H/7MVL/Trg\n42xnFIkF7mHkIoSIAB4CFkgp9wshppbU7FKq0dKlS5k2bRp5eXm0a9eOxYt19TlfeuklHnroIRYt\nWsSgQYP02//00098++23mJub4+TkxKuvvkrTpk2ZPXs2kZGRaLVazM3NWbBgAW3blj/hOn/+fCZP\nnsz777+Pvb29/rilLV68mMceewwhBJGRkfr1TzzxBElJSXTp0gUpJfb29vq+KnfTpUsXHn74YXx9\nfXFwcNCfylMMIzPrBE3t2tPCoSs7LdqQ7v4PAjv/392fqNyTwmINu05lsDEhlajEVNJzCzE3FXR3\na86Unu2I6OiIg63l3XdkQJUuuV9y1/yTwOvobmQcJaV8qhpiq1aq5H7tMmvWLKytrXnppZcMHYrB\nGOL381z2OeI2DKL39QTkkGM0tWtfo8evT3IKiog5lsamxFSij14mt7CYRhamhHk4EOnpSG8PB2wt\nq74LbWVVZ8n9HCllFvCSEOI9QP1JqSh1TE7uRebtnse7ez6lt6WGpg90J8DCxtBh1TlpOYVsPpLK\nxoQUdp7M4LpGS7NGFgz2aUE/Lye6uzXD0tzU0GHek3tJLutufCOlnCmEeLYK41HqmYyMDPr06XPL\n+i1bttCsWTP98qxZs2owqnpMq+FK4n8p2jcTy2wtIzuO450+79Cmsep/VFWS0q+xKTGFTQmpxJ+7\ngpTQpqkVE4PbEunlRJc2TTA1MY55k/tx1+RS0t2xtH03rfv1xrKU8tY77BTlDpo1a2bQ0v3K3zKS\nVtIscRZNsg5xpkFLBvWZzf91nGzosGo9KSUJF6+yMUGXUI6l5gDg1dKW6X3c6dfJkQ6ONkYzEV9V\nKjJyWQpIoLxXfuMxCYRXUVyKotSQnNyL7NsQSc/rCRQ1dMa8x0+4th4FdezDriYVa7T8mZTJppIJ\n+QtZ+ZgI6OrSlH8N9iTC05HWTa0MHWa1umtykVL2rolAFEWpefF/zcUx4TV6mGiIadiFwP4bMW/Y\n3NBh1Ur51zVsO5HGxoRUthxNJSuviAZmJoS2t+f5vu3p29GxRu6MNxbqJkpFqYc0xdfZsdKHnsXH\nOIM5CX4LCPP6h6HDqnWy8q6z5chlNiWmEHs8nfwiDbaWZvTp6Eg/L0d6uttjZVE/P2bvehOlEKJN\nRb9qIuC6LiUlhTFjxuDm5oa/vz8DBw4st/5WZcTExOjLu6xZs4b33nsPgNWrV5OYmKjf7l//+heb\nN2++7+PdydixY/Hx8WHevHll1t+ujL+Liwvp6brabsHBuuoIpcv6L1myhGeeeea+4vnoo4/Iy8u7\n7WNhYWG0adOG0pfsP/jgg3dtUZCVlcX//ve/O25z4/UYgqmZBSayiJiGXXEanYK3SiwVdjErn6U7\nk3jki934z97MiysOcOB8NqMDnPnu8W7EvxHBvIf96N+pRb1NLFA1cy43GHTORQjREXgeaAZslFJ+\naahY7pWUkuHDhzNx4kT9XfEHDhwgNTUVd3f3KjvO0KFD9bW6Vq9ezeDBg/XVid9+++0qO87tpKSk\nsHfvXk6ePFnp5+7cubPK49FoNHz00UeMHz9eX8fsZnZ2duzYsYMePXqQlZXFpUuX7rrfG8nlqadu\nvQWsuLgYMzOzank9d4wp+wyHNvbH0f8d3N1GEjL6hCqDXwFSSk5eztVNyCemcjBZV3LlAQdrpvVq\nR6SnEz7OxlFyxZjc9TdLStlbShle8u+dvqo8sQghvhZCXBZCHL5pfX8hxDEhxEkhxMySOI9IKacB\nDwP9brc/YxcdHY25uTnTpk3Tr/P19SU0NBQpJTNmzKBTp054e3vz448/AroRSVhYGKNGjcLDw4Nx\n48bp/8resGEDHh4edOnSRV98Ev7+a3/nzp2sWbOGGTNm4Ofnx6lTp5g0aZK+jMyWLVvo3Lkz3t7e\nPPbYYxQWFgK60cSbb75Jly5d8Pb25ujRo7e8loKCAiZPnoy3tzedO3cmOjoagMjISC5cuICfn5++\nnE1FlTdaOH/+PGFhYbRv375Mtejy2gtYW1vz4osv4uvry5w5c7h48SK9e/emd+/bTy+OGTNGn+xX\nrlzJiBEj9I/l5ubSp08f/Xvx66+/AjBz5kxOnTqFn58fM2bMuG2bgRuvZ9WqVfTp0wcpJZcuXcLd\n3Z2UlJRKvTd3s+bYGgYv7obf9eOkJK0CVH+VO9GVXLnCu78fIfzDP4iYF8sHm45jaiJ4ub8HW17s\nxeZ/9mJGPw98WxtPyRWjIqU02i+gJ9AFOFxqnSlwCmgHWAAHAM+Sx4YCG4CRd9u3v7+/vFliYuIt\n62rSxx9/LKdPn37bx37++WfZt29fWVxcLFNSUmTr1q3lxYsXZXR0tLS1tZXnz5+XGo1GBgUFyW3b\ntsn8/Hzp7Owsjx8/LrVarRw9erQcNGiQlFLKxYsXy6efflpKKeXEiRPlihUr9Me5sXzj+ceOHZNS\nSjlhwgQ5b948KaWUbdu2lZ988omUUsoFCxbIxx9//JZ4P/jgAzl58mQppZRHjhyRrVu3lvn5+fLM\nmTPSy8vrtq/xzTfflC1btpS+vr76L3Nzc5mWliallLJRo0ZSSllmH4sXL5ZOTk4yPT1d5uXlSS8v\nL7l3716ZmJgoBw8eLK9fvy6llPLJJ5+US5culVJKCcgff/xRf9y2bdvqj3GzXr16yd27d0tvb29Z\nXFwsIyIi5JkzZ/SxFBUVyezsbCmllGlpadLNzU1qtdpbXmd0dLS0srKSp0+f1q+7sQ8ppRw3bpyc\nP3++HDRokPz+++9vG8u9/H5mXDkul/7oL5mF9PnMR+4/u7XS+6gvCos0MubYZfnKyoMyYHaUbPvy\nb9LtlXVy/Je75be7kmRKdr6hQzQKQJyswOe3UZ8QlFLGCiFcblodCJyUUp4GEEIsB4YBiVLKNcAa\nIcQa4Jeb9yeEmApMBWjT5s5TRG+tTSDxYtW2AvVsaXvPTXq2b9/O2LFjMTU1xdHRkV69erF3715s\nbW0JDAzE2dkZAD8/P5KSkrC2tsbV1ZX27XXlOsaPH8+iRRUvAXfs2DFcXV31p+MmTpzIggULmD59\nOoD+r3d/f/8yo6LS8T77rO7+Wg8PD9q2bcvx48extbW943FfeOGFMiVebjQ/u5OIiAj9DZcjRoxg\n+/btmJmZldtewNTUlJEjR951vzeYmprSo0cPli9fTn5+fpmYpJS8+uqrxMbGYmJiwoULF/SN1m52\npzYD8+fPp1OnTgQFBTF27NgKx3ZH51divn0CY7V55AQ/zZTw/2JhWn+uVqqI3MJiYo5dZlOCruRK\nTmExVham9O7gQKSXI2EdHGjc0PAlV2ojo04u5WgFnC+1nAx0E0KEASMASyDmdk+UugKbi0BXW6xa\no7wHXl5eZSobV9TNJfprovXvjWPW1PHu5OZTEkKIO7YXsLS0xNS0ciU1xowZw/Dhw2+pFLBs2TLS\n0tKIj4/H3NwcFxeXMq0PSrtTm4Dk5GRMTExITU1Fq9Vich+nrNIyE7HY9yKNUzdgZefHKffXePqB\nUfe8v7rmRsmVTQkp7ChVcmWgdwv6dXIk2K15rS25YkxqY3K5LSllDOUklXtRE21AbxYeHs6rr77K\nokWLmDp1KgAHDx4kOzub0NBQPv/8cyZOnEhmZiaxsbG8//77t53vAN1oISkpiVOnTuHm5sYPP/xw\n2+1sbGzIycm5ZX2HDh1ISkri5MmTPPDAA2XK/FdEaGgoy5YtIzw8nOPHj3Pu3Dk6dOhQocnwyoqK\niiIzM5OGDRuyevVqvv76a6ysrBg2bBgvvPACDg4OZGZmkpOTc9sK0Dfeg+bNy7+/IzQ0lFdeeeWW\nUUV2djYODg6Ym5sTHR2tb2FQ3vt6O8XFxTz22GP88MMPLF26lP/+97/3VKBTarXs2vkC7mfm09AE\npO9sTD3/D3cT9Zf37UqutG7akEe760qu+LetGyVXjElFyr/koLsS7JaHACmlvPN5jqp3AWhdatm5\nZF2tJ4Rg1apVTJ8+nblz52JpaYmLiwsfffQRPXr0YNeuXfj6+iKE4D//+Q9OTk7lJhdLS0t9GX4r\nKytCQ0Nv+2E3ZswYpkyZwieffFJm1GRpacnixYsZPXo0xcXFdO3atcyFBnfz1FNP8eSTT+Lt7Y2Z\nmRlLliwpM8KqSoGBgYwcOZLk5GTGjx9PQICuYGtF2wtMnTqV/v3707JlS/2FBzcTQtz2A3/cuHEM\nGTIEb29vAgIC8PDwAHRlbUJCQujUqRMDBgwo0w7hZu+88w6hoaH06NEDX19funbtyqBBgypVAfly\n+kFOb+pPMJc4TCOyui/jgXbDKvz8ukZKyeELV9mUmMLGhBSOp+YCdb/kijGpdMn9mlYy5/KblLJT\nybIZcBzogy6p7AUekVImVGa/quS+UtuU9/t55PgynPZMoKGQ7G46kB6Rv2BmZty9PqpDkUbL3jOZ\nbExIISoxlYvZBZgICHRtSqSnE5Fejjg3qdslV2pCtZTcF0I0Adqjm9cAdJPulQ+vwsf7AQgDmgsh\nkoE3pZRfCSGeATaiu3Ls68omFkWpC7TaYkxMzHBtO4i4fa60DPqUsLYDDB1Wjcq7Xkzs8TQ2JaSy\n5ehlsvN1JVd6utvzz8gOhHs41KuSK8akwslFCPEEupsUnYH9QBCwi2q8cVJKedvLZqSU69E1KlOU\nekdKScwf02hx/lucR57G2sqJHg+fMnRYNSbz2vWSCflUtp1Io7BYi52VOX07OhLp5Uho++b1+s54\nY1GZn8Dz6BqD7ZZS9hZCeADvVE9YiqLczrnsc0xdOxW71I3MdGxM7rWLWFs5GTqsanc+M49Nibor\nvPYmZaKV0MquIWMD2xDp5UigS1PMTNVNocakMsmlQEpZIIRACNFASnlUCNGh2iIzECmlmuRTjI5W\nq6UgP41F33Rme44Zc/vMx6frNExM6uZf6FJKjlzKKZmQT+XIJd09Zx5ONjzT+wEivZzwammr/q8a\nscr8ZiYLIeyA1UCUEOIKcLZ6wjIMS0tLMjIyaNasmfqlVYyGLC4g4+Ix7PL281ATG54btwfXpu0M\nHVaV02glcUmZbExIZVNiCslX8hECAto24bWBHYn0cqRts/LvFVKMS4WTi5RyeMm3s4QQ0UBj4Pdq\nicpAnJ2dSU5OJi0tzdChKApSSgoL0rHUXsOy8Ax2TQQuwal1qiZYQZGGbSfS2ZSQwpajl8m8dh0L\nMxN6PNCcZ3o/QJ+OjtjbVM8l7Er1qsyE/r9us9oPqN4yujXI3Ny83PIcilKTzp7fTNYfI+lscpUz\nlu649ouCRnWjq0VW3nW2HtWVXPnjeBr5RRpsLM0I93Cgn5cTPd3tsW5QN0/31SeV+QleK/W9JTAY\nOFK14ShK/aYpvs62zQ8RmP4rdsD2FpMJ6fUl1PLRysWsfKISdae7dp/ORKOVONo2YJS/M5FejnRz\nbYaFWe1+jUpZlTkt9mHpZSHEB+juNVEUpYps/zWQsMID7DVxwLnPOno43vVeNaMkpeTE5Vw2Jegm\n5A9d0PVAcbNvxD96tiPSywmfVo0xUSVX6qz7GXtaobvnRVGU+1BcXMC1/Ms0tmmDS+CHbD+9gpDQ\n/9W6uRWtVrLv/BU2JaSyMSGFpAxdd0+/1na83N+DCE9HHnC4cwdPpe6ozJzLIf6uMWYK2AP/ro6g\nFKW+SEjZhzYqlDxTGwLHXKBt6z60bd3H0GFVWGGxhp0nM9iUmEJU4mXScwsxNxV0d2vOE6HtiPB0\nxNG2/pWiUSo3chlc6vtiIFVKadha64pSSxUVX2fuzv/w9h9v8469GcGdHjJ0SBV2taCI6KOX2ZSY\nSszRy1y7rqGRhSlhHg5EejrS28MBW0tVibm+q0hV5H/e4TGklP+t2pAUpW5LOrOGa9vHsu5iHsM7\nPsTEAZ9i38je0GHdUerVgpIJ+VR2nUqnSCNpbm3BUL+WRHo6EfxAMxqYqR4oyt8qMnKxKfm3A7ry\nL2tKlocAf1ZHUIpSJ2muQ8I7tE2YQ7qQvB/2Bj26Ge+V/KfScvXzJ/vPZwHQtpkVk0Nc6efliF9r\n1QNFKd9dk4uU8i0AIUQs0EVKmVOyPAtYV63RKUodkXhsGbb7nsdZm4FwGU/Tzh/So6GDocMqQ6uV\nHLyQzcaEFDYlpHAqTXf3gXerxrwU6U6klxPtHaxV9YpaKuXPc+x8cyPpl7W4D3AjbHbfaj1eZeZc\nHIHrpZavl6xTFKUchYVX2bVhAD1yd5KmNSE18AscOzyBsZxAul6sZffpGxPyqaReLcTURNDNtSmP\ndnchwtORlnYNDR2mci+0Wgp3xnP6l7/YsF5y9LgJ9qTxBz1xKbAibHb1Hr4yyeUb4E8hxKqS5QeB\nJVUekaLUEQlHFtMg7knCTAvZZt4en/5RNLa9tRNmTcspKCLmWJp+Qj6nsJiG5qb0crcnwtORPh0d\nsLNSPVBqpfPnKd6wmV3Lz7IhtiEmxUVE0Zc9BOHf5BQRoQXMfdaVbr2rv2laZW6inCOE2AD0KFk1\nWUq5r3rCqjwhRDvgNaCxlHKUoeNR6reCwiyaxz2BBkGc+78JDXjdoPFcvlpAVEkPlJ0lE/LNGlkw\n0LsFEZ6O9GjfHEtzYxlPKRWWnQ0xMVxdt411a4rZmerGGoZwDhfCxVbcm6Tx0NjmfP8itGvnVqOh\nVeomSillPBBf1UEIIb5Gd6nz5RvtjEvW9wc+RndfzZdSyvfuENtp4HEhxM/lbaMo1e3Ise9wd3sI\nywZ2nAr4ktatIwiwNsy9xicv5+pLruw7p5uQb9PUiknBLkR6OdGljZqQr3WKimDPHoiK4szaw6zZ\n15rfGMQfzGEqX2JPOoO9zxP2elsi+4XTuLHhQq3IpcjbpZQ9hBA5/H0TJYAApJTStgriWAJ8iu7U\n243jmgILgAggGdgrhFiDLtG8e9PzH5NSXq6COBTlnly7fo2lGx9l2tWVbDu7hl6RP+HVcXKNxqDV\nSvYnZ7GppGT96VIT8i9G6Cbk3R3VhHytIiUcPQpRUWiitrJnay5r8/rwGyOxI4LO7ONqiw5MH2tK\ncMv+dO1tTasuxjEVXpGrxXqU/Gtzt23vlZQyVgjhctPqQOBkyYgEIcRyYJiU8l3K3tCpKAa1/eRv\nTFo/naQrp2jpHUSf4I9q7NiFxRp2ncpgU2IqUYmppOUUYmYiCGrXjEnBLvTtqCbka53UVNi8GTZv\nJmfjTjZd6sRahvC7+AIrmct50YaeIRp6n/6a5s42zFxgRssAM6BmT3vdTWXKv4wGNkgpc4QQrwNd\ngH9X47xLK+B8qeVkoNsd4msGzAE6CyFeKUlCN28zFZgK0KZN3ShfrhhQUS571/fBLftPGou2bJ4Y\nTZhLWLUf9sYd8lGJqcQcSyO3sBgrC1PCOtgT6elE7w4ONLZSd8jXGnl5EBurSyhRUZw7eIW1DGGt\n+XiiNYu4jjkONteYolmEed5VHt70OB59nbl+7TEsGhnvhReVmXN5Q0q5QgjRA+gLvA8s5A4f+DVJ\nSpkBTLvLNouARQABAQHyTtsqyp3IS1GIP6cQcO0csY38iB21kUZW1XffSkr2jQn5FHafztDfIT/E\ntwWRnk50d2umJuRrC40G9u2DqCiIikK7fSd7i3xZa/Iga61WchA3TNAQ1vwYLzVZR8Snw+gR2oit\nr3jhHORM+54tAIw6sUDlkoum5N9BwCIp5TohRHVeKX0BaF1q2blknaIYTHbOOQ5s6EfPoqNg0x4R\nsY1e9iFVfhwpJafSckta/qZyoOQOeZdmVjwW4kqkukO+djlzRp9M2LqVa5kFRBHB2qbPss58NalF\ntpgg6eVTyPvDoU/Xa6wN/xlrE2t8XbIwM2tC5PuRhn4VlVKZ5HJBCPE5EAnMFUI0AKqzJvheoL0Q\nwhVdUhkDPFKNx1OUO4qLm03LI7MIMdEQ26gbPQdEg1nVzWfoS9YnphKVkMrpdN2EvK9zY2b060Bk\nScl6NSFfC1y5Alu3/p1QTp8mmVasbTyBtdZvsNXci8IiU2yLYcAQ6OefjmbNOq4cu8wLz/0TUwtb\n7Hc+Rkv/lpjU0iZqlUkuDwH9gQ+klFlCiBbAjKoIQgjxAxAGNBdCJANvSim/EkI8g64hmSnwtZQy\noSqOpyiVkZV9hkMbIwktPskpLDjS+XN6ej5eJfsuKLoxIf93yXozE0F3t2ZM7uFKREdHnBqrkvVG\nr7AQdu36O5nEx6PVSv6yCmWt82zWOvdhX7IDZEO7ZjDtKQj3vkwby8v4jetEfmYjFi+5RtD0IDRF\nGkwtTHHuVrvbZVUmueQDjYCxwNuAOZBVFUFIKceWs349sL4qjqEo9+rc+Q10KzpJTKPuBA1Yj2UD\nu/vaX3Z+ETHHdD3kY47pStZbNzCjVwd7Ij0dCevgQOOGakLeqEkJhw/rEsnmzfDHH5CXR56JNVvc\nn2St32f8dtabSxkWmJyE7t3hvWdgyBDo2BGO/3aM5UOXc7RZQ7xHd6Rh04Y8efjJOjUqFVJWbF5b\nCPEZoAXCpZQdhRBNgE1Syq7VGWB1CQgIkHFxcYYOQzFSmVknOHLoI0JCFwBwKTWOFvfRcvhStq6H\nfFRiKrtOZVCsldjbNCDC05FIT0e6u6mS9Ubv4sW/k8nmzZCSolvtFspvrZ9k7bXebD7kSEGBwMYG\n+vXTJZOBA0GmpbH7o91IrWToF0MpyisiflE8PuN9sGpe/aVYqpIQIl5Kedf/DJUZuXSTUnYRQuwD\nkFJeEUIY9+UKilJJWqllyf4lWMY9xUirQi6lTqaFY0ClE0vpHvKbElM5mKzrId/OvhFPhLbTTcg7\n26ke8sYsJ0c3Iim5RJjERABkc3v2+z/OmsDRrD3jRfyhBnAKXFxgyhRdQunVC7R5BeRl5NG0eVOO\nbs/g4LcH8Rnvg5QScytzgqYHGfb1VbPKJJeikrvmJYAQwh7dSEZR6oTjp1byeuwcVpz7iwdbBxAQ\nNBP3SiQVjVay75xuQn5TqR7ynduoHvK1QnExxMX9PW+ya5dunaUlBcHhbO32FmtzwvhtdzOSNwqE\ngG7dYM4cXULp1AlunNWK+zyOjS9spG3PtozfMB73we68ePFFLO3qz/xZZZLLJ8AqwEEIMQcYBRi2\nGp+iVIWiXPL2zaTdiQWM01gwYOjXTPSbiIm4+1U6BUUadpxMZ1NCKluOppKeex1zU0GwW3Om9GxH\n346qh7xRS0qC9et1ySQ6WlcIUgjo3JmUf7zJOssRrD3WgaitpuRthUaNIDIS3n4bBg0Ch5Jbm/LS\n89g97wBN3JrgMcwDe097fCb44D/FHwATM5N6lVigEnMuAEIID6APurpiW6SUR6orsOqm5lwUqdVy\n/OB/6JC0APKSOdOsL427/Y+mdu3v+LzsvCK2HtNVGP7jeBp51zXYNDDT95AP62CPjeohb7zS0+Gn\nn2DZMti5U7eubVtk3wgOuo9i7ZUerN3aiD9L+uy2bq0bmQwZAmFhYFmSI6RWoinSYNbAjHVPrSPu\nszj8/+HP4IV1uzpVRedcKpRchO4SBmcp5fm7blxLqORSv527EEPqHw/RlTRyrdywDvkG7IPL3f5i\nVr6+wvDu05lotBJH2xsT8k4EtWuGRS29H6FeuHYNfv1Vl1A2bdKd7vLyovDhR4lpNY618S1Z+5vg\n3Dnd5l27/p1QfH3/Pt11w97P9rLz/Z0EPhNI939258qZKxRdK8Khk3F1F60OVTqhL6WUQoj1gPd9\nR6YoBlRYeJWdUSMIytpCEyCm+TB69F0OZmVPWUgpOZaaQ1TJHfKHLugm5B9wsOYfPdsR6eWET6vG\nakLemBUV6U53LVsGq1frani1bs3px+fwe9NxbDjUiq3v6VY3bAgREfDGG7rTXS1alN2VtljLyQ0n\ncershG0rW7LPZWPX1o7mHs0BaOLaxAAv0LhV5lLkpcCnUsq91RtSzVAjl/pHqy3m2DJbOprms0u0\nwiV8FS0c/76SvlijJe7sFf0I5XxmPkJAlzZNiPR0JMLTkXb2akLeqEmpm4hftkx36is9nXy7FsQE\nzWRDo5H8frAlJ07o/iBwdYUBA3SXCoeH6xJMeZb0WsLZ2LP0nt2bnq/1RGolop7+YVGlp8VKdngU\neAA4C1zj734uPvcTqKGo5FJ/ZGadoImtG8LEhG1bJ9HQ9gECSjpD5l0vJva4ruVv9NHLXMkrwsLM\nhBC3ZkR4OtHX0wEHm/o1EVsrJSbqEsr33yOTkjhu4c3vni+yQQzgjyP2FBQILC2hd2/o31+XVB54\n4NbTXaAbtR5ZeYS/vviLyA8jcfByIGFFAiamJrgPcce0nhcIrY77XPrdRzyKUuM0Wg2rt82g7/l5\n7HKdTnDIPELDl5CWU8jyP88RlZjK9pPpFBZradzQnD4eDkR4OtLT3Z5GDSrVpFUxhORkWL4cli0j\nd/8Jtoq+bGj1Cb83DyMp3Qb2Q4cOMG2aLqH07Hnn0UnmqUyatGuCEIJtc7aRl57H1fNXcfBywGu0\nV829rjqiwv+DpJRnqzMQRalK8ee384/fp3M0JZ5V7VrQwLo/C/84xaaEFPadz0JKcG7SkEe6tSHC\n05FAl6aYmaoJeaN35Qr88gvyu2Uk/JHO7/Rng+0itpl2oUhjSqMr0KcP/F9/XUJxdb37LosLi1nW\nfxlJMUk8tuMxWge3ZuyasVi3sMZE/U7cM/XnmVK3XM9i97o+WF/ZR3pmMCNdVvJuahNOH74GHKVT\nK1te6OtOhKcjHk42daqWU51VUAC//Ub2klVs3ljMhuK+bDD7nmR0s+6d2sD0AbpkEhICDRrceXc5\nF6yNxF4AACAASURBVHOInRNL6oFUJsVMwqyBGXYudoS/E07TB5oCYOtcFd3b6zeVXJQ6QWq15J/8\njl3bvmNTRgi/X50BxTbsuCro7tZQtfytbTQatFuiOfDpNn7fZMqGwjB2shQNZtg2KiaivylvloxO\nnMspHiylRAjB8XXH2fPxHmxa2vDgkgcxbWDKoe8O0TKgJXkZeVg7WjNs8bCafX31wF2TixCiwv2A\npfz/9u48ruoqf/z4633ZEQSRRVHAXcFdUXEBUQFRS802S1vUcqxsb1qmaWqmmWxqlvxVU98yc01t\nbNGpXHHf9yXFFTeUxRUB2e/5/XGuyThqisC9wHk+Hjy498O9n8/7GN0353zOeR91/PbCMYxbcz63\nkNlrlrBi62Z2ZbclTz2Ll6sQG1GPeFNhuGpRirNJO1jyj90sWO7OovwYMogDoGOzi7xyj4XEgRAV\n5YzLVf9Ji/KKcHZ3RkRY/NJikr9Jptdrveg8tjOFOYXkZuYS3CUYAM+6nrxy/pUaO9urstxMz2Xq\nTZ5LAX1vIxbDuCnHz15i8d50Fv2cytZjWVixEOjcmJ4NT/Jw3P1ENQ0wCxqriJIS2Pr9CRZ8dJiF\nG3zYlN8OKx3xc8kmoed5BjxaSMIdrtSrd2WYSlkVuZmXqBVYi0tnLzEjYQYZuzJ4JuUZfEJ8UEoR\nHBmMT5gPAG3ub0Ob+9v813VNYql4v5pclFJ9KiOQ2yUiscDbwB5gtlJqhV0DMsqNUordJ7N+KVm/\nLz0bgCZuR3kqcAMBdU6QOPgjAuu2ufGJDIeQkQGL/53FwilpLNpZj7PFIQgN6Oq9jzfid5I4vhld\n+tXGyckbgNzMXLLTrHjX9yb5u2TmjZqHT4gPT+x+Ag8/D2o3rE2ThCa/3Hzv/3czsdUROMSwmIhM\nBu4AMpVSbUodTwQmoneinKSUevdGlwdyAHcgtSxxGI6jsNjKhpSzvySU9Iv5WAS6NPLj94PC8c//\nkrbp75Hd7j06tJ9u73CNGyguho0bYcH3BSycm83Wo/6AD4HkM7DOGgYMgPhXOuHfLoKivCLSt6dj\nLfDAydOFeaPmsWPKDnr8tgfx78Xj18yPNsPb0LC7vtEiIgyfN9y+DTSu6VcXUYrI8ps8l1JKlWlY\nTERi0Ilh2uXkYivvfwCIRyeLzehdMJ2ACVedYjRwRillFZEg4B9KqRE3uqZZROl4LuYXsXxfJkv2\nZrBy/2myC4rxcHEipoU/fVv64Zn2W7zdSohN/B5ltVJcko+LS9XaaKmmOHkSFi2ChT+VsGRhCRdy\nXXGimO6sJ7HORgbc5UH7Z2M57xrE+SPnaT6gOdZiK3+t81cKcwp5eNnDNO7TmL1z93Lh2AWaxDWh\nXvt69m6WQTkuoqyMYTGl1CoRaXTV4a7AIaVUCoCIzAaGKKUmoHs513Me+JXJiIajOHUhj6XJ/71D\no7+XK4Pa1Sc+IoiezfzZfGodT/54J5Pc91LsGYqyWhGLBReLSSyOorBQFxhesAAWLlTs2qXvaQRL\nBnern0j0XkePwf5kt++Je9d7COvdiEMLDzFzwMe4eLrwatarWJwtxP01Du8G3tTvqKcZR9wTYc9m\nGbfBkaciNwBKV2FOBbpd78UiMgxdRcAX+Og6rxkLjAUIDb3p0T6jHCml2JeezeI9GSxJTufnkxcB\nvUPjmOjGJEQE0SGkDk4W4cy5ZNZ+15vhB/bj6RXKmdhZDIowQyCO4tgxWLhQJ5SkJMjJARdLMb1c\nNvJX5tHVdSeurZsTPr4fdR76Pxa8uJRNL2+i1V35hPVuRMOohgyePJiG3RoiTjoZdXmySu6ablyD\nIyeXW6KU+hb49lde8xnwGehhscqIy9AFITcdPffL/ZPU87ogZMeQ6+zQqKxw6HN8tr5ETHE2/+ww\nlLsSZlDLtZb9GmGQnw+rV1/unUCybTenMN8LPOT8E5EsIVCdYVA/hYwcwYdvNODc9gvUcm5CHRcX\nIp+IpPV9ranfSfdK3H3d6Tiqox1bZFQkR04uJ4GQUs8b2o4ZVUBugS4IuWRvBsv2Z3LBVhAyupk/\n4/s0o194EAHe/zt6uf/gHHx2v0a9/CO4BPYmI/xNRjaoEhMWq6VDh670TpYvh7w8cHNTDKy/nYd8\nlhCRtY7BF+ZzvlMcH27rRaZXC/pOfRZPf0/uaZmGd7A3XvX0Hw4B4QF2bo1RmRw5uWwGmotIY3RS\nGQ48aN+QjBvJzM4nKTmTxXvSWXv4LIXFVnw9Xehr26Exuvn1C0JezEll+6I76ZW/g/NWCwVRn+PW\nbAxBpjxLpbp0SSeRhQv116FD+nhiwBae89lN25BjDDn8D3YfDWcRiRS17YzM+4A6jRrxZPIZ/Fv5\n/7KG5HIPxaiZHCK5iMgsIBbwF5FU4E2l1BciMh5YhJ4hNlkptceOYRpXUUpx+HQOi23DXTtsBSFD\n/DwY2S2M+IggujSqc8OCkMpqZf36F2mc8v+ItlhZ4xpBu4QfcPO5iYqDRrkoKIAff4SpU/UML0tB\nHgnOSdzleZKQx6IYcHoaG384w+kSPyIKduD5wjg6DL2PTu3b4VLLFdD7bwREmJ6JccVN7+dS3Zip\nyGVTYlVsP355Q60MjpzJBaBtAx+9oVbrIFoG3XxByJRjCwhbM5ADVg+skR/TOnxURYZv2CgFmzfr\nhPLjzAs0y9qCv8cl6o+9k8TQPez6/b+pV3ySO4q+oY4vFA+7F+eHR0B0NFhM9YOarCL2czFqqLzC\nEtYeOsPivekkJWdyNrcQFychqkldRvdsRFxEEPV9br4gZH7BBXZuf49uUe/QJGwAOy58QNvWT+Dk\n7FqBrTBArz+ZPh2+/iKbU4dyyXKvx4NR5whdtY7mIYrh81/AcuQwCW7uyNA7YcQUSEzE+ddKDRvG\nVUxyMa4pMzufZcmZLE3WG2rlF1nxdnMm1nb/pHfLAGq733pByMWHF3NqxX086pFFSv1omoQNoEP7\nZyugBcZlly7Bd9/pXsrSpdBVbWAwi3BrFsKTfwzCe+Yn5KmV1DqYqzdDefMN5K67oLYpO2+UnUku\nBnBl/UlScgZLkjPZeeICAA18Pbg/MoS4iCC6Na5b5oKQaRlbmLDqTT7c+xO9/JvQruNrdAobUJ5N\nMEpRSk8bnjoVNsw6QkTeVix+YbzxRhcGNXDj0vQS2u5/D58R+yE4mFqvPwdjxkCjRvYO3agmTHKp\nwQqLrWw8cpalezNYmpzJyQt5ALQP8eWlhBbERdza/ZNrshZTkPwPvLe/yp35EBD7J17u+TJuzmaY\npSKkpMC0qYqfPj9JcpoPeHnzeNND1D2eQuzAM/Rc/zosWaLvmwwcCI+/p787m48Co3yZ36ga5nxu\nIcv3Z5KUnMnKA6fJKSjG3cVCr2YBPN23GX1bBRJY271crnXs0GzCDryL24WdnK3dnhax/yS+oVmz\nUt4uXoR//1v3UlavhpHMZBCHeejuPoz+fX2cvlyCy1fTcJqRAaGh8Mc/wujR199lyzDKgUkuNcDh\n0zkkJWewdG8mW46dw6og0NuNO9vXJy48iB5N/fFwdSq3653POszPi+8gumgf+S51ce81l+CQYWDW\nrJSbkhJdcmXqFMWxuZsJL9pFfthg/vKXQHq4heN7IIfwvX/CrWMSODnB4MHw+OOQkKCfG0YFM8ml\nGiousbL12HmWJmeQlJxJim26cHj92jzVpxlx4UG0beCDpZw3TFJWK2tXjaXVicl0tyhWeHSmc8J8\n3L2Cy/U6NVlyMkydVMSaKYdYe64Vvr7CU7W34ecHr7+eTKMdE2DaNDh3Dpo0gXfegUcfhfpmQaNR\nuUxyqSay84tYeeA0ScmZLLeVW7k8XfjRno3o2yqQhnUqroqwslrZPLs+vchkN16c6/Ylsc3uqbDr\n1SRnz8Ls2XrY6+fNeTzHB8RTyNh3HuO+39SB2a64fzUZHl0LLi5w1126l9K3r1mTYtiNSS5V2Ilz\nl0hKziBpXyYbUs5SVKKoYyu3EhceRHRzf7zLMF34VlzKO4OHmx9isVBQL5HVCD1jJ2GxmF+t21FU\npOt5zfw0m4uL1tPQepyCtqP58989CD/Ti1ZNCmi0fQLSZAZkZUGLFvD++/DIIxBgVsob9mc+AaoQ\nq1WxM/UCSbb1J5e3+20aUIvRPRsTFxFEp1Bdrr6iKaVI2j6R8D0vsT3scXrGfEJ036kVft3qTCnY\nsQOmfZjFqm/PsC2rKY3qwiOymfqxTZn45Rk8l86H2Z/p5fVubnDPPTB2rF45b+5pGQ7EJBcHl1dY\nwppDZ1i6V/dQzuQU4GQRIsPq8PrAcPqFB9IkwOvXT1SOjpw7xNMLn2PZoR/5Mcwb/7qmbPrtSE+H\nGTMU06YJRbv3MZw5DHT15K3vXyRxgBdqQz9cZ0yGto/rTVNat4YPPoCHHgI/P3uHbxjXZJKLA8q4\nqKsLJ9lWxxcU69XxMS0DiA8PIrZlAL6ediiVUlLIhiX34JX+AxvSPHg77m/06vYMLk4VO/RWHeXn\nw/z58M0HJ5D168nFE4+ud/Cbv4XR9Gws3e5vTJ21n0K3z3V3xsMD7r9f91KiokwvxXB4Jrk4AKUU\nyWnZttldGexMzQKgYR0PHugaSlx4EF0b+5V5dXy5xJi+DNnyFFEX97HBLZidj/1Eg4D2dounKlIK\n1q9TzHrvBIuWuXAwpz4xdc7Tx/M4rR/qxL2fKNiwAz6bAt3n6M1TOnSAf/0LHnwQfHzs3QTDuGkm\nudhJQXEJG1LO2dafZHAqKx8R6BDiy2/7tyQuPIgWQV63tzq+HBw9voS0taPork5Crcao3j8Q1WCQ\nXWOqao4fh2nTFNOnCy0O/EgkWxkc1ob+395N7+jWOF8MxjJrJrR9BvbsAS8vGDlS91I6dza9FKNK\nMsmlEp3LLWT5Pn0zftWB0+QWluDuYiG6eQDPxbWgT6vAa+7OaA9nziWzZ/lweuTvwl/BWr94eibM\nQ5xvvvpxTZaTA99+C4ve2YrX/q2spzv1YtoycFg7OjUJpePwlrhuWwmPfQ5z5+pNVbp0gc8/18Nf\n3t72boJh3BaTXCqQ3kwr95fhrq3Hzv+yOn5whwbERwTSo6k/7i6OtWJ69fJRtD85hZ4C61zDadVn\nNj3929k7LIdntcKyhYV8N2Efs7c159wlDx7yOkFwffjwLRd6jwVOe8DUryHyXjhwQFcefuwxvS6l\nvRlmNKqPapNcRCQCeAs4CyQppebaI47iEitbjp23FYPM4OjZSwBE1K/N+L7NiQsPpE1w+a+Ov10l\nxYUUlVzC3c0XV89gki318O/+OTGN77B3aA7vwH7F9BnC9GmKu47/i0CyGBkzhHv/0oFukXfg4mqB\nZcvg/jd17fuiIujRA373O7j3XvCsuMWthmEvDrETpYhMBu4AMpVSbUodTwQmorc5nqSUevcG53gR\n2KSUWi0i85VSg290zfLcifJifhEr958mKTmD5ftPk5VXhKuThaimdYkPD6RveBANfB1zOEkpxfJ9\nc2i4ZRSnfLoQO3iVvUOqEi5cgNnTCtj07nLqpO3lcxlLjwQv7mq+mz5DfGjeNwTJSIcvv4QvvtDl\niv384OGHdU+ldWt7N8EwyqSq7UQ5BfgImHb5gIg4AR8D8UAqsFlE5qMTzYSr3j8amA68KSKDgboV\nHfCJc5dYmqx7JxtTzlFs1avj48KDiAsPJLpFAF5ujvLPe207jyXx4soJJB1J4rsQL+oFRNk7JIdW\nXAz/mXaehZ8eZequjhQWuPCiywFqt2nIhn8VEBHtBSUReiP6u5+F//xHV5iMjYW334Zhw8C9fCpO\nG4ajc4hPP6XUKhFpdNXhrsAhpVQKgIjMBoYopSageznX8pQtKX17rR+KyFhgLEBoaGiZYt167By/\n+/Zn9mfo1fHNAr0YE92Y+PAgOlbS6vjblnuM3UuG0jh7B6cy6zAxcSIDI8fh6mS2Gb6W3TutTJth\n4acvM7jv7KfUB8aOasbDT3rTsf14nFwscOIE/PEj3Us5cUKXYHnhBd1LadHC3k0wjErnEMnlOhoA\nJ0o9TwW6Xe/FtuT0O6AW8P61XqOU+gz4DPSwWFmC8vdyo04tF34/KJy48CAa+dcqy2nsIuviUdz2\n/wP3w58RgWK1dxTrhs7B17tsibY6S02Frz9M58D/raAkK4cpzo8xcEAg9QMHMPSVlgQ299ZdmR9/\n1DO8FizQd/Tj4+Hvf4chQ8DVJGuj5nLk5HJLlFJHsfVKKlJY3VrMHtu9oi9TrgoKLrJ++cO0Oz0f\nVycFjR/Bqd3bxNYKsXdoDuPSJVi5UrFsehpb115i+fFmBCKMdj6JT5+2nJhRQr1gJ6ArHDkCv38f\nJk+GtDRdzv7VV/U2wU2a2LsphuEQHDm5nARKf/o1tB0zbpLVWsyGdS/Q8MgnxDoVs8VSF+8u/6Jl\n8/vsHZrdKQW7d+vbI4sX6x0c2xdsYiALifQKZOD7zRg4MIhWLZ/H4mSBwkKY+53upSxZok8yYAB8\n8gkMGmS2CTaMqzjy/xGbgeYi0hidVIYDD9o3pKpl9Xdd6V2wnf24s7XV20R2etXeIdlVZqbOC4sX\n6y+39KNEsRG/gECeeqoPvTu0wv+CE50eaoO7r+1NBw/DpEkwZYo+QcOG8Ic/6G2Cy3jfzjBqAodI\nLiIyC4gF/EUkFXhTKfWFiIwHFqFniE1WSu2xY5hVwqGUeXh4BtKgXneCO77BmhML6NHrYyw1sLhk\nQQGsW3eld7J9OwRwGi9fZ2IT6xCl0ihYdoIeLzag1ysAPqA6w/798NUyvXJ++XK9LfAdd+hyLP37\nm22CDeMmOMQ6F3soz3UujuBU9in+mfQS7+TOYoNLc6LvO2DvkCqdUnrR++VksmIF5ObqEasePSAu\n7wesm7fS+YlI7vjXIIryirA4W3A6eUIvcrz8lZamT9i0KYwapb+CzVbNhgFVb52LUUbZOaf4cdWz\njNnxE0UlRbRr15+BMRPtHValOX8ekpJ0Mlm0SBeJBGjWDB7vl0LjC9uJfy2S8MQw9n3fjPNH6tIu\nLgi++gqXy8nkyBH9psBAvTVwv376e+PGpmikYZSRSS5VlbUIDk9CbXmRe615rGk+iOf7TaSpX1N7\nR1ahioth06YrvZNNm/QM4Nq1dU54+fHz9L+3Ns1aOjEjcR0nd52kMKUefLeNVklJOpm8kKxP5uur\nFzg+/7xOJhERJpkYRjkxyaWKUVYrmza9Tof02bhdOoqrX1eSQ8fyUcQYe4dWYY4evZJMkpL0lvEW\niy4i/PvfQ0ICdO2imHPnTA6/cZjiZkMgJY07w3ZRK2M5zuNf02Nmnp4QE6OHufr21XulmPsnhlEh\nTHKpQn7eO4mSbS/SzXKRNIsv9WPm4d7gTtpUs7+2s7P1/ZLLQ10HD+rjISG6zmP//tCnjyJ33wl2\nzdhFl1YxuKzdSuDFg4SFHqbByD5QcgEfV1fo3h3eekt3a7p0MQsbDaOSmORSBRw7kcSpNY/QXZ0k\nU1lYFTSCHr0ngXP1qFNlteqZXJd7J+vW6cLBnp561Gr8eN07adkSCi7m417LmcLVG/gsMQlVUkL7\nL54lpCiFhMvdmRFP6J5Jjx6m4rBh2IlJLlXA0Y3P07nkJCt8YonsN4sYz3r2Dum2nTp1Zb3JkiVw\n5ow+3rGjLsnVv7/ODW5ugNVK9tKNfHVXEscOFvK8xyd45JxhBKHUa10X14QhOplER5utgA3DQZjk\n4oAu5Z9lU9KD+IUOoV3bJ2kb9w2XivOIrcIbduXl6VXwl3snP/+sjwcF6YXu/ftDXJx+jlKofftI\nf3MxqYv20OX4N3ieyyKHx4jyO4safC8M7ENobKwuEGkYhsMxycWBlFhLmLZzGhNX/o7V/ulsLcqG\ntk/i59vc3qHdMqX0dvCXk8mqVZCfr295xMTobU0SEqBdO9sErWPH4KdlqKRlyPJlnDjlzJeMxkUC\naffAENwS+/CbPn30CnnDMByeSS6OQCm2bXuXo7vfZ/Sx83Rt0JU9XT8itsXd9o7slpw+DUuXXhnu\nOnVKH4+IgHHjdO8kJsZ2GyQ9Xa9+/0ivNUlPyWUVvcl3qc3Dd8fQMLYPd54LI/w30bj5mfsmhlHV\nmORib+e2w46X6ZS+lLpOzswb8gl3tv8NUgVmgBUWwvr1V3on27bpHoufn648n5Cgv4eEoFc7rlgB\nr9gWLu7dyxn8uVQrgNC4thQPieHY1ELaP9oR63vxWJwsdLJ3Aw3DKDNT/sVOUtPWcnTVw/QsOYK4\n1qEw4jVoNg5XVy+7xfRrlIJDh64kk+XLISdHl1fp3l0nk/79oVMncMrLgTVrrpRUuZx5PD0hOpqN\nbjEsnF9E/U71GLv1NyilsBZbcXIx604Mw5GZ8i8OKuviMbYn3UtU7mb8ga1+/YiMm4urq++vvtce\nSkp0bvjmG51Ujh7Vx5s2vXLfpE8fqO2aDxs2wH+WwfPLYONGvZzexQW6d+fE6DfZcqIeQX0i6PFq\nNE33nyE+5gDtRupJCiJiEothVCMmuVSi1cseofXJ6cRYFGtdmtKk93Qi6znexmNK6Y7GzJkwa5a+\nPeLtbSuv8rJOKE3DimHrVp15PkyCtWv1HXuLBSIj4aWXOB/eA+ee3fBuGkjKn1ayf+56/HpZAfBv\n6Y9/S387t9QwjIpikksFs1qLAbBYnLEWZZNi8cO768dEN7/fzpH9r5QUnVBmztRV511d9T5YI0bA\noAFW3A/u1snkuWWwcqVeSg96yte4cXqtSUwM+Pjw45M/suWRLfT4rTvx78UT9VwUPX7bAxePmlf6\n3zBqIpNcKtCGfXPw3jyGcw3vI7rPZKLjv8Zicax/8tOn4euvdUJZv14f690bXnoJ7u6VQZ2NC2Hu\nQhi39MpKx+bN4cEHdTLp0wcCAkjdkMr2ydvpXD+X4EgfQnqE4FXfiw6PdADArbabnVpoGIY9ONYn\n3U0SkSbA64CPUuoe27Fw4FmgLrBIKTXJXvElp2/l5eVvseTgD6wMdcPFXQ//OEpiyc2FefN0Qlm8\nWN8aadsW/vpOCQ8020zItnnw8UJ4fId+w+WVjnFxOpmE6N2ns9Oy8arrhQBLfruEtG1phPQMITgy\n+Jd7KYZh1EyVPltMRCYDdwCZSqk2pY4nAhPRu05OUkq9exPnmns5uZQ6ZgHmKKXuvdF7K2S2WF4a\nWxYOwu/idqLSvHmx1+s80+0ZPFw8yvc6ZVBcrNegzJwJ332nE0xICDw4KIsRAYtpu/srXXI4O1tP\n/+rZExIT9Ve7dvpeio21xMrXw77mwA8HGLloJE3imnD24Fm86nnh5m16KIZRnTnybLEpwEfAtMsH\nRMQJ+BiIB1KBzSIyH51oJlz1/tFKqcxrnVhEBgNPAp+Xf9jXl5ObhvvBj3De/wEdrQWs8epA8hPf\nUtencWWG8T+Ugs2bdUKZPVtvAe/rq3gwJpURteYRvftfWD617W0SGqqHuhIT9XBX7dqlzqPI2JnO\n7q920/ftvji5OuEZ4EnPV3pSt0VdAOo2r2uPJhqG4aAqPbkopVaJSKOrDncFDimlUgBEZDYwRCk1\nAd3Ludlzzwfm2xLTN+UT8fUVF+exbsUYWqXNwcvJCqH34dT+HXp723fDroMHdUL56iv92M3Vyp3h\nhxnRYBYD9v4dtwUXdUXI2FgYN1YnlJYtf9koKzstm0PfbCf/fD7dX+iOtdjKpK56lDF8WDgNuzVk\n8KTBdmyhYRiOzjFuAkAD4ESp56lAt+u9WETqAn8BOorIa0qpCSISCwwD3IEV13nfWGAsQGhoaJmD\nvbxhV8DBfxDjVMhOqc2Zjh8QET6qzOe8XRkZMGeOTiqbNoGIok/wfl6tO5m7z/4fPjsv6gQybpRO\nJjExKA8PRIQ9/97Dgb98T2hMKJ0f70zGrgzmj55PrcBaRD0fhZOLE3fPupuw3mF41jWlWAzD+HWO\nklxuiVLqLDDuqmMruE5SKfWaz4DPQN9zKev118wNJ7r4AIdxZWPT39G1y9tIqXsSlSUnB77/HmbM\nUCxdCiUlQgevg7xvmcQD1hk0yLoI/fphTXiXix2iqd2jDdYSK7OHzCZ9zKcMnz+c4M7BHF9znCPL\nj+AfricehPYKZfz+8dRpWueXMjThw8IrvX2GYVRdjpJcTgIhpZ43tB1zSJ7NRrHq3G56xH5B00re\nsKuoSM/wmvllAfN+cOJSgTNhTqm8UjKNEcykaSNPMtsn4Nzn/+ChBFJWpTJr8Czcai/kpfQ2WJws\niAiN+zXG2U3/5+//9/4MmDjgl2u41nL95V6KYRhGWThKctkMNBeRxuikMhx40L4hXV/nTq9W6vWU\ngg1rS5g58TRzfvLmzKVa+JHDw8whwWM1dZrUpdPoDtQevpR5r29mx5QddHApYcgYnSQ6/6Yz9TrU\nQ1kVYhEe+M8D/3V+i3Pl97oMw6jeKj25iMgsIBbwF5FU4E2l1BciMh5YhJ4hNlkptaeyY3M0+1af\nZubf0vhqWRApOUHUp4QH+IJm3pmMe9YN10Hx/G3oJXL35BLcMobawcFEPhFJq2GtqN+pPgA+oT4k\n/jPRzi0xDKOmMVWRHUlRESkz1/P91Cy+2tgUp7xsOrONYosrPbtb6dHXnX///ThB7esxatUoLM4W\nTu89jU+oD65ervaO3jCMGsCR17kY6IWI5w6e48KmAwSe38+3Uy+yb3sBHuQxnRH4epUwLOIAAV4W\nOjzWgU6Pd0ZZFa++qbA4XRnGCogw2/wahuF4THKpJJk/Z3Js1TFqB7rTsnYahz9dxlffeWBF+Bsv\ncAkvhrr+RNtWRcz+kweRQxoALf/rHGIRBMffRMwwDMMklwqy5t01nNpyiujf9aK+Vw47nv2W9csK\nCLRk8oG1Dd/zDPXIJM/dl0eHFjLyaUVU94FUgQ0oDcMwfpVJLmVUXFBMYXYhnv6eZO7JZMH4BWSd\nyOLpg08jOTns+mQtJVk55K75f+zNOMdKHuR7p+HsLWmBh5uVoYMVIx4JJiFB76dlGIZRnZjk/Ajt\n7AAACXJJREFUUgYr/riC1X9eTev7WjNs5jDcvFwpOptF47oXKY6Nx2X9Ku4sCmSO6yPcyzR20ASL\nRREfJ7w6AoYOdcLb296tMAzDqDgmuZRBgy4N6DG+I2EemTB6ND4LF/JYWhpZ1GZ6w+eYGfwpy483\nRRUKXdrDxJFw//1CUJC9IzcMw6gcJrncqo0baf7n52m+cSNYrRT4BrGg9UvMaHAPP+wOoyBVaNYM\n/vAHXWS4RQt7B2wYhlH5THK5VT4+WEsUqx+axMwLA/n3ykAurBUCAmDsWL0lcNeumBvzhmHUaCa5\n3KJVma0YmbaeE5ugVi0YOhRGjtSbNDqbf03DMAzAJJdb1qSJ3hL43XdhyBCdYAzDMIz/ZpLLLWrY\nEH780d5RGIZhODZTDtcwDMModya5GIZhGOXOJBfDMAyj3JnkYhiGYZS7KnlDX0SaAK8DPkqpe2zH\nooER6DZFKKV62DFEwzCMGq3Sey4iMllEMkXk56uOJ4rIfhE5JCI33EdYKZWilBpz1bHVSqlxwA/A\n1PKP3DAMw7hZ9ui5TAE+AqZdPiAiTsDHQDyQCmwWkfnoLY8nXPX+0UqpzBuc/0FgzA1+bhiGYVSw\nSk8uSqlVItLoqsNdgUNKqRQAEZkNDFFKTQDuuNlzi0gokKWUyi6ncA3DMIwycJR7Lg2AE6WepwLd\nrvdiEakL/AXoKCKv2ZIQ6B7Llzd431hgrO1pjojsL2O8/sCZMr63qjJtrhlMm2uG22lz2M28yFGS\nyy1RSp0Fxl3j+Ju/8r7PgM9u9/oiskUpFXm756lKTJtrBtPmmqEy2uwoU5FPAiGlnje0HTMMwzCq\nIEdJLpuB5iLSWERcgeHAfDvHZBiGYZSRPaYizwLWAy1FJFVExiilioHxwCIgGfhaKbWnsmO7Bbc9\ntFYFmTbXDKbNNUOFt1mUUhV9DcMwDKOGcZRhMcMwDKMaMcnlFohIiIgsF5G9IrJHRJ61d0yVRUSc\nRGS7iPxg71gqg4j4ishcEdknIski0t3eMVU0EXnN9rv9s4jMEhF3e8dU3q5VIURE/ERkiYgctH2v\nY88Yy9t12vy+7Xd7l4h8JyK+5X1dk1xuTTHwolIqAogCnhKRCDvHVFmeRd8PqykmAguVUq2A9lTz\nttsWNo8FOiul2qCrYwy3Z0wVZAqQeNWxV4EkpVRzIMn2vDqZwv+2eQnQRinVDjgAvFbeFzXJ5RYo\npdKUUttsj7PRHzgN7BtVxRORhsAgYJK9Y6kMIuIDxABfACilCpVSF+wbVYW7CBQBHiLiDHgCp+wb\nUvlTSq0Czl11eAhX6hFOBYZWalAV7FptVkottk2kAtiAXv5RrkxyKSPbX3odgY32jaRSfAC8DFjt\nHUglaQycBr60DQVOEpFa9g6qIimlzgF/A44DaegySovtG1WlCVJKpdkepwNB9gzGDkYDC8r7pCa5\nlIGIeAHfAM8ppS7aO56KJCJ3AJlKqa32jqUSOQOdgE+UUh2BXKrfUMl/EZGmwPPoxBoM1BKRkfaN\nqvIpPX22xkyhFZHX0cP9M8v73Ca53CIRcUEnlplKqW/tHU8l6AkMFpGjwGygr4jMsG9IFS4VSFVK\nXe6VzkUnm+osElinlDqtlCoCvgVqyp5IGSJSH8D2/UZV16sNEXkUXRh4hKqANSkmudwCERH0OHyy\nUuof9o6nMiilXlNKNVRKNULf4F2mlKrWf9EqpdKBEyLS0naoH7DXjiFVhv1AlIh42n7P+1HNJzGU\nMh94xPb4EWCeHWOpFCKSiB7qHqyUulQR1zDJ5db0BB5C//W+w/Y10N5BGRXiaWCmiOwCOgDv2Dme\nCqWU2oHeY2kLsBv92VDtVq5fq0II8C4QLyIHgTjb82rjOm3+CPAGltg+xz4t9+uaFfqGYRhGeTM9\nF8MwDKPcmeRiGIZhlDuTXAzDMIxyZ5KLYRiGUe5McjEMwzDKnUkuRrUnIjlXPX9URD6yPR4nIg/b\nHk8RkXtsj4+KiP9tXLPD9aapi0isiCgReeyq1ysReelXzjv0RsVSS7fnJuOMEpGNtumoySLylu34\nYBGp1lUJjIrlbO8ADMOelFLlP79fF37sgF71/tN1XvYzcB9XioE+AOy8idMPBX7gGos6RcS5DO2Z\nCtynlNopIk5ASwCl1HzMVuPGbTA9F6NGE5G3btBbeFlEdovIJhFpZnt9gIh8IyKbbV89S51nuois\nBaYDfwLut/UI7r/GuY8B7iISZFsRn0ip4oEi8rjt/Dtt1/MUkR7AYOB923mbisgKEflARLYAz15u\nj4g4294fazvfBBH5yzXiCEQXqkQpVaKU2mt7fene3Y5SX3ki0ltEatn2CdlkK+455Nb+5Y3qzvRc\njJrAQ0R2lHrux839VZ6llGprG2b6AF2HaSLwT6XUGhEJBRYB4bbXRwC9lFJ5trpNkUqp8Tc4/1zg\nXmA7sA0oKPWzb5VSnwOIyJ+BMUqpD0VkPvCDUmqu7WcArkqpSNvztwCUUsW2GOaKyNPo5NXtGjH8\nE9gvIiuAhcBUpVR+6RcopTrYzn0numTIOuCP6FJAo0VvNLVJRJYqpXJv0F6jBjHJxagJ8i5/QMIv\nBfsib+J9s0p9/6ftcRwQYftQB6htq5INMF8plXcLcX0NzAFa2a5RulBkG1tS8QW80EnseuZc66BS\nao+ITEcPo3VXShVe4zV/EpGZQALwIHp4Lvbq14lIc+B9oI9SqkhEEtAFTS/3+tyBUGpOPTLjV5jk\nYhjXp67x2AJEXf3XvS3Z3NJf7UqpdBEpAuLRO32WTi5TgKG2eyGPco0P/FJudN22wAX08Nf14jgM\nfCIinwOnRaRu6Z/bkufXwOOl9j0R4G6l1P4bXNuowcw9F8O4vvtLfV9ve7wYXdQS0LO8rvPebHRh\nwF/zB+AVpVTJVce9gTTRWzyMKMN5EZFh6CHAGOBDucY+6SIySK50w5oDJehkVNpk4Eul1OpSxxYB\nT19+r4h0vJmYjJrDJBfDuL46tqrIz6I30gJ4BogUkV0ishcYd533LkcPn13vhj4ASql1Sqnvr/Gj\nN9C7nK4F9pU6Phv4re0metPrndc2jfpd4DGl1AF0FdyJ13jpQ+h7LjvQExFGlE50IhIG3AOMLnVT\nPxJ4G3ABdonIHttzw/iFqYpsGIZhlDvTczEMwzDKnUkuhmEYRrkzycUwDMModya5GIZhGOXOJBfD\nMAyj3JnkYhiGYZQ7k1wMwzCMcmeSi2EYhlHu/j+Hjt2E0L4azQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8571295ba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# LU factorization via Gaussian Elimination\n",
    "norms = []\n",
    "sizes = [2,4,6,8,10,12]\n",
    "for n in sizes:\n",
    "    A=hilbert(n)\n",
    "    x = np.ones(n)\n",
    "    b = A.dot(x)\n",
    "    result = linear_system_lu(A,b)\n",
    "    norm = np.linalg.norm(result-x,ord=2)\n",
    "    norms.append(norm)\n",
    "fig,axes = plt.subplots()\n",
    "axes.plot(sizes,norms,color='red',label='GE')\n",
    "\n",
    "\n",
    "# LU factorization via Cholesky Factorization\n",
    "norms = []\n",
    "for n in sizes:\n",
    "    A=hilbert(n)\n",
    "    x = np.ones(n)\n",
    "    b = A.dot(x)\n",
    "    result = linear_system_cholesky(A,b)\n",
    "    norm = np.linalg.norm(result-x,ord=2)\n",
    "    norms.append(norm)\n",
    "axes.plot(sizes,norms,color='blue',label='Cholesky')\n",
    "\n",
    "\n",
    "# qr decompositions\n",
    "for method,color,style in zip([classic_gram_schmidt,modified_gram_schmidt, householder_qr],['green','orange','purple'],['--','-.',':']):\n",
    "    norms = []\n",
    "    for n in sizes:\n",
    "        A=hilbert(n)\n",
    "        x = np.ones(n)\n",
    "        b = A.dot(x)\n",
    "        result = least_squares_via_qr(A,b,qr_decomp_method=method)\n",
    "        norm = np.linalg.norm(result-x,ord=2)\n",
    "        norms.append(norm)\n",
    "    axes.plot(sizes,norms,color=color,linestyle=style,label=method.__name__)\n",
    "\n",
    "hilbert_condition = lambda n:np.power((1 + np.sqrt(2)),4*n)/np.sqrt(n)\n",
    "axes.plot(sizes,[1e-12*hilbert_condition(n) for n in sizes],label='Condition of Hilbert Matrix')\n",
    "axes.legend()\n",
    "axes.set_ylabel(r\"$\\vert\\vert$ residual $\\vert\\vert_2$\")\n",
    "axes.set_xlabel(r\"Hilbert Matrix Size\")\n",
    "axes.set_yscale('log')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I find that the errors induced by the gram-schmit qr decomposition increase much more rapidly relative to the householder qr decomp. as well as gaussian elimination and cholesky LU factorizations. The LU factorizations and the Householder QR decomposition appear to scale equally well in terms of accuracy for small n in the linear Hilbert system.\n",
    "\n",
    "Furthermore, the householder, cholesky factorization, and gaussian elimination methods all scale with the condition number of the hilbert matrix, indicating that the accuracy of the solution is dominated by the condition of the input matrix. This is not the case for the gram schmidt qr decomposition."
   ]
  }
 ],
 "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.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}