{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "48HONqnFpIvI"
   },
   "source": [
    "# Problem 7\n",
    "Let $A$ be the $20000 \\times 20000$ matrix whose entries are 0 everywhere except for the primes $2,3,5,7,\\cdots,224737$ along the main diagonal and the number 1 in all the positions $ a_{ij} $ with $ |i-j| = 1,2,4,8,\\cdots,16384$. What is the $(1,1)$ entry of $A^{-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JiZSIfm1qR-F"
   },
   "source": [
    "### Some introductory code "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "executionInfo": {
     "elapsed": 1558,
     "status": "ok",
     "timestamp": 1603893832461,
     "user": {
      "displayName": "Husnain Raza",
      "photoUrl": "",
      "userId": "16146658556003895350"
     },
     "user_tz": 300
    },
    "id": "wvpmKjZMqRiB",
    "outputId": "9ad6fbfd-f2b3-4ad4-eca6-97830b9d01b0"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(554466, 3)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def is_prime(N):\n",
    "    return N > 1 and all([N%i != 0 for i in range(2, int(N**0.5+1))])\n",
    "\n",
    "primes = [i for i in range(2,250000) if is_prime(i)][:20000]\n",
    "entries = [(i,i,primes[i]) for i in range(len(primes))]\n",
    "for diff in [2**j for j in range(15)]:\n",
    "    for i in range(20000):\n",
    "        j1, j2 = i+diff, i-diff\n",
    "        if 0<=j1<20000: entries.append((i,j1,1))\n",
    "        if 0<=j2<20000: entries.append((i,j2,1))\n",
    "entries = np.array(entries)\n",
    "print(entries.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vsdtHzxT3WQc"
   },
   "source": [
    "### Attempt 1: Using standard libraries\n",
    "\n",
    "Note that A is very sparse - with only $\\frac{554,466}{20000^2} \\approx 0.14\\% $ of its entries being non-zero\n",
    "\n",
    "![Plot of upper left 100x100 quadrant](media/prob7/plot1.png)\n",
    "![Plot of entire matrix](media/prob7/plot2.png)\n",
    "\n",
    "We therefore avoid calculating $A^{-1}$ explicitly (note that this is intractable as matrix inversion has a time complexity of $O(n^3)$) \n",
    "\n",
    "Instead, note that if we have $\\hat{v}$ be the first column of $ A^{-1} $, then we get that $ A \\hat{v} = (1,0,\\cdots, 0)^{T}$. We use [BIConjugate Gradient iteration in scipy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.bicg.html#scipy.sparse.linalg.bicg) to solve this sparse system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 229
    },
    "executionInfo": {
     "elapsed": 789,
     "status": "error",
     "timestamp": 1603893823785,
     "user": {
      "displayName": "Husnain Raza",
      "photoUrl": "",
      "userId": "16146658556003895350"
     },
     "user_tz": 300
    },
    "id": "3QZE1MQao2Zz",
    "outputId": "91fb2579-b63b-4656-d445-18ff6f771c42",
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7250783462684014"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.sparse import coo_matrix\n",
    "from scipy.sparse.linalg import bicg\n",
    "\n",
    "row, col, data = entries[:, 0], entries[:, 1], entries[:, 2]\n",
    "\n",
    "A = coo_matrix((data, (row, col))).tocsc()\n",
    "b = np.zeros(A.shape[0])\n",
    "b[0] = 1\n",
    "\n",
    "iterations = 0\n",
    "curr_xk = None\n",
    "change_iter = []\n",
    "def callback(xk):\n",
    "    global iterations, curr_xk\n",
    "    iterations += 1\n",
    "    if curr_xk is not None:\n",
    "        change_iter.append(xk[0] - curr_xk)\n",
    "    curr_xk = xk[0]\n",
    "x, exit_code = bicg(A,b,tol=1e-50,callback=callback)\n",
    "x[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEGCAYAAACdJRn3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAArm0lEQVR4nO3deXhU5fn/8fedhCTsCARUFhOQpVjcCLiiuIPQWq3FrbVa+rO2WpfaBbdabV1au9r6rfvS1mrRWqWI2uIC1iISKquAIIJElB0EgZDA/ftjJskkzExmksyczOTzuq5z5Zxnzpy55xhz85xnM3dHREQkGTlBByAiIplHyUNERJKm5CEiIklT8hARkaQpeYiISNLygg4gHbp37+7FxcVBhyEiklHmzJmzwd2Lor3WKpJHcXExZWVlQYchIpJRzGxVrNf02EpERJKm5CEiIklT8hARkaQpeYiISNKUPEREJGlKHiIikjQlDxERSZqSRwIqqvbw6pK1vL9+e9ChiIi0CK1ikGBTnfrr6azetBOAq04ZwPdOGxhwRCIiwVLNowHrt1XUJA6Ae15Zxp69zr2vLWfrzsoAIxMRCY6SRwNeXvTJPmWT533E3S8v5fYX3g0gIhGR4Cl5NCDaMr0fbgzVRCaVlQPw1oqNfPWhWVTt2ZvW2EREgqI2jwbc+s99axe/mfZeneOrn3qHtZ9WsHZbBb26tE1XaCIigVHNowFVe/etedSXnxe6jbc8v6hO+bptu1i3bVdK4hIRCZJqHk20cXttg/q0xWvrvDbi9lcAWHnX2LTHJSKSSqp5NMGQAzrx0ZadDZ8oIpJllDyaYHtFFbk5VnNc0r19gNGIiKRPRiYPM/ucmd1nZs+Y2beDimPbrkqumzSv5nj/ToVRz/vGY7PTFZKISFqkPXmY2SNmts7MFtYrH21mS81suZlNjHcNd1/s7pcD44HSVMYbz/aKKpZ8sq3meEflnpr93VW13XZfXbIurXGJiKRaEDWPx4DRkQVmlgvcC4wBhgAXmNkQMxtqZlPqbT3C7/ki8B/glfSGH9J7v7ZU7qntiWUGOyqqao6fmVNe5/zvPvlO2mITEUm1tCcPd58BbKpXPAJY7u4r3H038BRwlrsvcPdx9bZ14etMdvdjgYuifY6ZXWZmZWZWtn79+mb/HkUdC+p+HrBsXe3EiU++/WGd1/85bw0vLfy42eMQEQlCS2nz6AWsjjguD5dFZWajzOweM7sfmBrtHHd/wN1L3b20qKioeaMFuneomzyqh4N8sOEzAFZv3rHPe+au3trscYiIBKGljPOwKGUxR+e5++vA66kKJuJzYr5Wv+ZRbeP2Cg7sUkj3DgXs36mwTptIx8KWcrtFRJqmpdQ8yoE+Ece9gTUBxVIj3ujyono1j5vGfg6ACx58ix8/t4jl67bTqbBNnXM0C6+IZIuWkjxmAwPMrMTM8oHzgckBx0RFVeyJDrvXq3kM7dUZgMo9zt/KQk/gKvfuxSLqVHNWbY77ebvjfJ6ISEsSRFfdJ4GZwCAzKzezCe5eBVwJvAwsBia5+6J410mHioiut/VF1jwO6FxI+4J9H0m1y8/lhe+OpCA899WcVZspnvhC1OvNW72FgTe9yH+Xb2hi1CIiqZf2h/DufkGM8qnEaPwOSnXN4/azP89Jg3owe+Umlq3dzh9eW06/otrR5C9fewIfbd53mpJ2+XkMObATS382pk7S2F21l/y8HP4660M6tc3jsN5dOOveNwG48KFZ/O2yo9mwfTdjDz0gxd9QRKRx1IIbR3XyaJefy4Fd2nLW4b1wd6446eA6j6M6Fbah4/773srSg/ar2R9f2rtm/Y+du/eQn5fDDf9YAMDJg3vUed95D7wFwIiSU2M2zIuIBKmltHm0SBVVocdWBXm5NWVmRtv83JpHUZHlb994Sp2y/zeyX83+L849rGb/zHveqNOTK9YI9F1xHpuJiARJySOOispQzaN+ooBQsqivR8e6c1vl5NQ954qT+gPw0ZadzCtveMzHyF+8xv8+jN/ILiISBCWPOHbvqU4euQ2cmZjDenep2X+l3tofkc45snZ85NNlq2OeJyISFCWPOGpqHm2i36Y7zh7Ko5cOr1P2/h1nxrxeu/zadpHfv7o86jmvf38UF47oW3OssSEi0hKpwTyOI/p2YepVIzmoW7uor194VN99ynJzjKcvP4Zog9Pb5jdcgynu3p4du2vbOvYksAyuiEi6qeYRR/uCUFfbaGM44hle3JURJV2jXG/f5DH6kP0B6NWlLdN/MGqf8yITiYhIS6GaRxq1a7Pv7f7jV4/kL7M+5IxDetY0uEcmqzeWbeCwW//FE988isI2uRzYpbDO4y8RkSDor1Aa9excd8zGvRceiZnxtaMPqlNefwLFrTsrueTR2WzYXsHne3XimcuPpbBN8zTii4g0hh5bpVFBXi4Lbz2Dwft35O5zD405grwgL5e3bziFW74wpKZsw/YKABZ+9CmDb34pLfGKiMSi5JFmHQryeOmaE/hKaZ+45/XoVMhFRx0U8/UN2yuYPG/NPtPGT5m/hjeWrWfn7j3cP/19Xlr4CVt3VrJ83bYYVxIRSZ4eW7Vg+Xk5FLbJYVflvrPtlv5sGgBbduxm1geb+M34w9lZuYcr/xpa7rake/uahamGF+/H7JWbWfLT0XrcJSLNQjWPFm7eLafHff3Hzy/ihfkfM3f1Fm55fmFNeXXiAJi9MjRKPdY0KCIiyVLyaOEiR7ffdc5QZl5/ctTzbpuyiOfmxl8/683lG9hVuYeylaEl5N2dneoKLCKNoOSRQToU5nFA57YcFWUMycKPPgWg935ta8q6tKu7kmF+Xg5fuW8m5943k4feWMGwn03jsFv/xSdbd6U2cBHJOkoeGWDej0/nuycfXDOgMN6Y8++dNrBm/0uH96rz2qNvrmTBR6EJGX/2wmI2fbab3Xv28o4mXxSRJCWUPMysl5kda2YnVG+pDkxqdW7XhutOH0Rebug/V/+iDlHPy8/LYeSAIp674jgARg0qSuj6sz7Y1DyBikirYfW7eu5zgtnPgfOAd4HqB+Tu7l9McWzNprS01MvKyoIOo9lsr6jiH+98xIMzVvDhph1ce+pAunfMr9O1t3q1wu0VVbzz4WYK8nIZf//MmNdcedfYdIQuIhnEzOa4e2m01xKpeXwJGOTuZ7r7F8JbxiSObNShIK/OqPRDDuy0z5iQ/PAaJB0K8hg5oKjOXFuRgw+r3fniYi588C0tQCUiCUkkeawA2jR4VhqZ2Sgze8PM7jOzUUHHE5T7vzaMI/p24ej+3ZJ635eH9eaPFx3JnJtOZeSA7qFrTV/Bf9/fyOtL1Z1XRBqWSPLYAcw1s/vN7J7qrbEfaGaPmNk6M1tYr3y0mS01s+VmNrGByziwHSgEyhsbS6b73AGd+Md3jqNDgrP+PnP5MfzyK4fRqbANY4YeQLcOBdx5ztA65/xn+YZUhCoiWSaRvzqTw1tzeQz4A/Cn6gIzywXuBU4jlAxmm9lkIBe4s977vwG84e7Tzawn8GvgomaML2uVFneltLhuN9/e+9WuVTK0V2cWrfk03WGJSAZqMHm4++Nmlg9U9wFd6u6NXt7O3WeYWXG94hHAcndfAWBmTwFnufudwLg4l9sMFER7wcwuAy4D6Nt330WbpNZDF5eyZWclryxey+yVmxl7zxv8/MuH8vlenYMOTURaqAYfW4XbFJYRqhn8H/BeCrrq9gIiF+suD5fFiukcM7sf+DOhWsw+3P0Bdy9199KiosS6rLZWpw7pybnDepOfl8OG7RUsWvMpNz63kIqqPWzbpWVwRWRfiTy2+hVwursvBTCzgcCTwLBmjMOilMXsQ+zuzwLPNuPnC1CQV/tviXmrtzDoptDU7+rGKyL1JdJg3qY6cQC4+3s0f++rciByjvLeQPyJmqTZ/eCMwVHLtY66iNSXSPIoM7OHw91jR5nZg8CcZo5jNjDAzErC7Svn07yN9JKAoo61zUdtcmsrgys3hmbodXcee/MDtuzYnfbYRKRlSSR5fBtYBFwFXE1opPnljf1AM3sSmAkMMrNyM5vg7lXAlcDLwGJgkrsvauxnSOMtu30M7/1sDLNvPLWm7JRfTWftp7tYtOZTfvLPd7nqqbnBBSgiLUKD05Nkg2ybniRd3J2S66fuU94+P5dFt40OICIRSadGTU9iZpPCPxeY2fz6W6qClZbDzLjk2OJ9yj/TGiAirV683lZXh3/GG2chWW6/dvlRy6fMX8O4Qw9MczQi0lLErHm4+8fh3e+4+6rIDfhOesKToPXoVNuIPrBnB04Z3AOAK//6Dn95axXFE1+geOILmlBRpJVJZJzHacCP6pWNiVImWejsI3qRY6E1REqLu1K1Zy8H3/giADc9Vzs92fzyrYwo6craT3fRo2MBZtGG7ohItojX5vFtM1tAqFdUZHvHB4DaPFqJwja5nDe8b82cWHm5Odz/tX3Hh46/fyYfbtzBUXe8UiepiEh2itdV96/AFwiNt/hCxDbM3b+ahtikhTojvBwuwGlDetbsv7d2GwBPzPow7TGJSHrFfGzl7luBrcAFAGbWg9AU6B3MrIO76y9EK/bYpcOp3OOcNqQnZ//fm7zz4Ra++afa7tBbd1bSuW2LWgZGRJpRIhMjfsHMlgEfANOBlcCLKY5LWrhRg3rU1DoOjrKm+g3PLkh3SCKSRomMMP8ZcDTwnruXAKcAb6Y0KskoVRFzX5V0bw/ACws+jnW6iGSBRHpbVbr7RjPLMbMcd3/NzH6e8sgkY0wcM5hDe3dmRElXBu/fif43hEalb9heQfcOUZdbEZEMl0jNY4uZdQBmAE+Y2e+AqtSGJZmkZ6dCLj2uhEMO7ExujnH3uYcCcMmjb/P3Oa12lWCRrJZI8jiL0Drm1wIvAe8T6nUlEtWXj+wNwMKPPuW6p+cxZ9XmgCMSkeYW97FVeG3x5939VGAv8HhaopKMlpNTd4Dgl//4Xw7u0YFp3zsxoIhEpLnFrXm4+x5gh5lpMWtJSv22juXrtgcUiYikQiIN5ruABWb2b+Cz6kJ3vyplUUnG+/e1J7Cjcg8vL/yE26a8C0BF1R4K8nIDjkxEmkMibR4vADcTajCfE7GJxLRf+3x6dWnLpccVc0TfLgBMmr062KBEpNk0mDzc/XFgEvCWuz9evaU+NMkGZsZvxh8OwM3PL2LrjkqtiS6SBRIaYQ7MJdTTCjM73My0vrgkrLh7+5r10Q+77V/0v2EqT5fV1kJWrN/OD56eR0WVpnUXyRSJPLb6CTAC2ALg7nOBkpRFlAAzG2lm95nZQ2b23yBjkcQ8fumIOsc/eKZ2YuaTfzWdp+eUU7ZyM61hWWSRbJBI8qgKT5IYqdH/h5vZI2a2zswW1isfbWZLzWy5mU2Mdw13f8PdLwemoO7DGWHIgZ32Kbv3teVs+mx3zfFFD83ikTdXpjEqEWmsRJLHQjO7EMg1swFm9nugKf/afwwYHVkQHk9yL6FFpoYAF5jZEDMbamZT6m09It56IfBkE2KRNPrvxJP5zXmH8dI1IwG4++WlTHt3bZ1zfjrlXbWJiGSARJLHd4FDgApCa3xspXZ986S5+wxgU73iEcByd1/h7ruBp4Cz3H2Bu4+rt60DMLO+wFZ3/7SxsUh6HdilLWcf0ZvB+9fWQn749/l0LKzbY/yZOeqVJdLSJZI8xrr7je4+PLzdBHyxmePoBUT+xSgPl8UzAXg01otmdpmZlZlZ2fr165shRGlOz37n2Jr9bbuqmPGDk2qOr9d07iItXiLJ4/oEy5oi2oLXcZ9duPst7h7z8Zm7P+Dupe5eWlRU1OQApXkd3rsL7fNDAwbHDj2Avt3asfDWMwDY6zDi9mkAvLjgY+aXbwkqTBGJIeYIczMbA5wJ9DKzeyJe6kTzz6pbDvSJOO4NrGnmz5AWJCfHWHjrGZRv3kmPTqFuvB0K8rjznKFc/+wC1m2rYN7qLXz7if8BsPKusUGGKyL1xKt5rAHKCE1PEjmyfDJwRjPHMRsYYGYlZpYPnB/+HMliZkafru3qTFkyvrQPFx3VF4Cz7q1dc6x44gtpj09EYouZPNx9Xngkef/IkeXu/qy7N3qObTN7EpgJDDKzcjOb4O5VwJXAy8BiYJK7L2rsZ0jmys0xbj97aNTXIgcWikiwLNagLDOb5O7jzWwBUdof3P3QVAfXXEpLS72srCzoMCQJqzZ+xol3v75PuR5fiaSPmc1x99Jor8WbVbe6O+645g9JJL6DurWne4d8NmzfzZKfjmbwzS8B8NaKjRzdr1vA0YlIzJpHNlHNIzN9uquSyqq9dOtQwG+nvcdvpy0D4NFLh3PSoB4NvFtEmipezSORrroigehU2IZu4UWlvjmyX035pY/ODiokEQlT8pCM0KEgjwcvrv0H0OpNOwKMRkQSmZJ9n6lIopWJpNppQ3ryrRNCNZCRv3iNXZWawl0kKInUPL4epeySZo5DJCHXnjawZn/WB/WnSBORdImZPMzsAjP7J1BiZpMjtteBjWmLUCRCYZvaAYVff+RtXlm8Ns7ZIpIq8brq/hf4GOgO/CqifBswP+o7RNLgr988igsfmgXAhMfLmHPTqeTl5NC5XZuAIxNpPeKNMF/l7q8DpwJvuPt0QsmkN9EnMhRJi2MP7s6gnh1rjof9bBqH3fYvZryn2ZNF0iWRNo8ZQKGZ9QJeAS4ltKCTSGBuHPu5fcrKVjV61hwRSVIiycPcfQdwDvB7dz+b0Gp/IoE5YWARK+44k+siGtDveWUZry5RG4hIOiSUPMzsGOAioHpq03htJSJpkZNjfPeUASy/fUxNWfUodBFJrUSSwNWEFn/6h7svMrN+wGupDUskcXm5OQzs2YH31m5nfvlWbnpuAX9560PycoyvHn0QN439HHm5Gg8r0pw0t5VkhQ3bKzj119PZsqNyn9duHjeECceXBBCVSGZr0txWZlZkZneb2VQze7V6a/4wRRqve4cCXrhqZNTXZr6vYUkizS2RuvwTwBKgBLgVWElo5T+RFqVXl7bM/fFp+5TPfH+DpjIRaWaJtHl0c/eHzezq8FiP6WY2PdWBiTRGl3b5rLxrLAs/2kqnwja8s3ozVz81lw837WBgxNgQEWmaRGoe1Q+RPzazsWZ2BKGBgiIt1ud7daZvt3aUdG8PwNVPzWV31d6AoxLJHokkj5+ZWWfgOuD7wEPAtSmNSqSZVCePxR9/yvefnhdwNCLZo8Hk4e5T3H2ruy9095PcfZi7T05HcLGY2RAzm2RmfzSzc4OMRVq2joW1811NnreGj7bsDDAakewRs6uumf0eiNmP192vatQHmj1CaF30de7++Yjy0cDvgFzgIXe/K841rgPedvc3zGyyu38x3meqq27rVrVnL3dMXcIjb35QUzb9B6M4qFv7AKMSafka21W3DJgTZ2usx4DR9QLMBe4FxhCa+uSCcO1iqJlNqbf1AP4MnG9mdwPdmhCLtAJ5uTl868R+dcpOvPt1KveoDUSksWL2tnL3x1Pxge4+w8yK6xWPAJa7+woAM3sKOMvd7yRUS4nminDSeTbai2Z2GXAZQN++fZsjdMlgPTsV8up1J3Lyr2o7Cs5asYnjB3QPMCqRzNVS5mzoBayOOC4Pl0VlZsVm9gDwJ+DuaOe4+wPuXurupUVFRc0arGSmfkUdWHnXWF657kQAHnhjRcARiWSuljLBYbT1QeK1t6wkXKsQSVa/cA+sGe+tZ/bKTQwv7hpwRCKZp6XUPMqBPhHHvYE1AcUiWc7M2C+86uBX7ptJa5jfTaS5JTK3VYmZ/drMno1cy7yZ45gNDAh/Vj5wPhBod2DJbq9cN6pm/+JH3g4uEJEMlUjN4zlC81n9ntBa5tVbo5jZk8BMYJCZlZvZBHevAq4EXgYWA5PcfVFjP0OkIV3b59fMg/XGsg3s3avah0gyEmnz2OXu9zTXB7r7BTHKpwJTm+tzRBrSpV0+4w49gCnzP6bfDVOZ9r0TOLiH5r8SSUQiNY/fmdktZnaMmR1ZvaU8MpE0+P7pg2r2X1m8LsBIRDJLIjWPocDXgJOB6lFVHj4WyWh9urar2V+3rSLASEQySyI1j7OBfu5+Ynhuq5PcXYlDskJujjF4/9Cjqof/84F6XokkKJHkMQ/okuI4RALz0jUn1OyXXD+VhR9tDTAakcyQSPLoCSwxs5dT2FVXJFBXnzKgZn/c7/8TYCQimSGRNo9bUh6FSMCuOXUAry1dx/zyUK1j82e72a99fsBRibRciaznMZ3QGuYdw9vicJlI1jAzJl95PI9eMhyA0387Q+0fInEkMsJ8PPA28BVgPDBLCzBJthoZnmV3/bYKVm3cEXA0Ii1XIm0eNwLD3f3r7n4xoenTb05tWCLByMvN4ZJjiwGYu3pLoLGItGSJJI8cd48cPbUxwfeJZKQbzvwcANf8bW6wgYi0YIkkgZfCPa0uMbNLgBfQNCKSxfLzcvh8r04AvLTwk4CjEWmZ4iYPMzPgHuB+4FDgMOABd/9RGmITCcxvxh8OwOV/mcPG7Rp5LlJf3K667u5m9py7DyPGcq8i2ejgHh1q9qctXst5w7WUsUikRB5bvWVmw1MeiUgLYma8c3NoyvYf/X2BpmwXqSeR5HESMNPM3jez+Wa2wMzmpzowkaBFDhIcfvu0ACMRaXliPrYysxJ3/wAYk8Z4RFqUX48/jO9NmsfGz3azu2ov+XnqaCgC8Wsez4R/PuLuq+pv6QhOJGjnHNmbcYceAMDAm15kjx5fiQDxk0eOmd0CDDSz79Xf0hWgSNC+cXxJzX7/G6ZSuWdvnLNFWod4yeN8YBehR1sdo2wircLhvbtw7akDa45fXaIVB0WsocnfzGyMu7+YpniifX4/QlOkdHb3c2OVxVNaWuplZWWpDVSy3p/fWsXNzy0EYOVdYwOORiT1zGyOu5dGey2RWXUbnTjM7BEzW2dmC+uVjzazpWa23MwmNvD5K9x9QkNlIqn2taMPqtnftqsywEhEgpfqriOPAaMjC8wsF7iXUC+uIcAFZjbEzIaa2ZR6W48UxyeSlPu+OgyAI3/6b2795yKKJ77Ava8tDzgqkfRLafJw9xnApnrFI4Dl4drDbuAp4Cx3X+Du4+ptjX64bGaXmVmZmZWtX7++Cd9CpNYx/bsBULnHefTNlQD8/tVlAUYkEoxE1vNoZ2Y3m9mD4eMBZjauCZ/ZC1gdcVweLov1+d3M7D7gCDO7PlZZfe7+gLuXuntpUVFRE8IVqdW5bRu+M6p/nTJ31ANLWp1ElqF9FJgDHBM+LgeeBqY08jMtSlnMVnt33whc3lCZSLp061BQ57iiai/Pz13DucN6BxSRSPol8tiqv7v/AqgEcPedRE8AiSoH+kQc9wbWNOF6Imk1vrQ3Fx3Vl3m3nM6LV48E4PtPz2ODZt+VViSR5LHbzNoSrh2YWX+gKf+XzAYGmFmJmeUTGk8yuQnXE0mrjoVtuP3soXRu24bB+9cOeTrhF68FGJVIeiWSPG4BXgL6mNkTwCvADxO5uJk9CcwEBplZuZlNcPcq4ErgZWAxMMndFzUqepGAmRnPXB56ortj956AoxFJnwbbPNz932b2P+BoQo+rrnb3DYlc3N0viFE+Fa1GKFmitLgr5w/vw9//V87WnZV0btsm6JBEUi6R3lZHAgcBHxNqm+hrZv3NLJHGdpFW4bzhfajc40yavbrhk0WyQCKPrf4PeAt4AHiQ0GOop4D3zOz0FMYmkjEO690FgNunLmbNlp3BBiOSBokkj5XAEeExE8OAI4CFwKnAL1IYm0jGyMkxjuzbBYDJ89R5ULJfIsljcGSDtru/SyiZrEhdWCKZ5+nLj6WwTQ7zy7fwwIz3efg/HwQdkkjKJNJusdTM/kjoURXAeYQeWRUQHvshIpCbY4w79ECemVPO1AWfADAhYi0QkWySSM3jEmA5cA1wLbAiXFZJaH1zEQk7qGu7Osfzy7cEE4hIiiUyJftOd/+Vu5/t7l9y91+6+w533+vu29MRpEimuPjY4jrHX/zDm6zbtiuYYERSKJGuuseZ2b/N7D0zW1G9pSM4kUzTuW0bXrnuRO46Z2hN2YjbXwkwIpHUSKTN42FCj6vmABpCK9KA/kUdKOnWnqdmr2bu6i0ArN9WQVHHgvhvFMkgibR5bHX3F919nbtvrN5SHplIBsvJMZ674riaxaOG3z6N8s07Ao5KpPkkkjxeM7O7zewYMzuyekt5ZCJZ4Jh+3Wr275+up72SPRJJHkcBpcAdwK/C2y9TGZRItujcrg2XhBvR//zWKtxjLl0jklESmRhR3XFFmuAnXzyEddt2MXXBJ5RcP5VHLxnOSYN7BB2WSJMkNLmhmY0FDgEKq8vc/bZUBSWSbX5x7mE1Awd/O+09JQ/JeIl01b2P0Kjy7xKakv0rhGbZFZEEdSjI4z8/ClXi55Vv5f7p7wcckUjTJNLmcay7XwxsdvdbCa1l3qeB94hIPb33a8eXDj8QgDtfXBJwNCJNk0jyqJ5feoeZHUhoWhJN2CPSCNedPqhm/4lZqwKMRKRpEkkeU8ysC3A38D9CU7Q/Fe8NIhJdn67tePwbIwC48R8LqajSuFvJTInMbfVTd9/i7n8n1NYx2N1vTn1oIWbWz8weNrNnIso+Z2b3mdkzZvbtdMUi0hxOHFjET886BIDx97/Fnr3qviuZJ5GaB2Z2rJldSKjh/CwzuzjB9z1iZuvMbGG98tFmttTMlpvZxHjXcPcV7j6hXtlid78cGE9oDIpIRvnaMcW0z89l3uot9L9hKlV79gYdkkhSEult9WdCgwKPB4aHt0T/YD8GjK53vVzgXmAMMAS4wMyGmNlQM5tSb4vZn9HMvgj8B9Csc5KRLj+xf83+HVPVgC6ZJZFxHqXAEG/E0Fh3n2FmxfWKRwDLq1ciNLOngLPc/U5gXBLXngxMNrMXgL8mG5tI0K48+WDGDD2AU389nUfe/IBuHfK54qSDgw5LJCGJPLZaCOzfjJ/ZC1gdcVweLovKzLqFx5ocYWbXh8tGmdk9ZnY/MDXG+y4zszIzK1u/fn0zhi/SPMyMg3t04MYzPwfA3S8vDTgikcTFrHmY2T8BBzoC75rZ20BF9evu/sVGfqZFKYtZqwnP4Ht5vbLXgdfjfYi7PwA8AFBaWqoWSWmxvjmyhFkfbGLa4rXMWbWZYQftF3RIIg2K99gqVZMfllN3kGFvYE2KPkukxTMzrj9zMNMWr+WGZxfw3BXH0TY/N+iwROKK99jqI6DK3adHboRqCeVN+MzZwAAzKzGzfOB8YHITrieS8foXdeDH44awdO02xt8/M+hwRBoUL3n8FtgWpXxH+LUGmdmTwExgkJmVm9kEd68CrgReBhYDk9x9UTJBi2Sji47uC8CCj7by8qJPAo5GJL54yaPY3efXL3T3MqA4kYu7+wXufoC7t3H33u7+cLh8qrsPdPf+7n57oyIXyTIFebk8c/kxANw3/X0+q6gKOCKR2OIlj8I4r7Vt7kBEBEqLu/KVYb1558MtXPjQrKDDEYkpXvKYbWb/r36hmU0A5qQuJJHWrbqxfN7qLWzbVRlwNCLRxettdQ3wDzO7iNpkUQrkA2enOC6RVuvbo/rz6pJ1lG/eybzVWzl+QPegQxLZR8yah7uvdfdjgVsJzaS7ErjV3Y9xd7XmiaTIAZ3bMvXqkQB89eFZ7KrUzLvS8iSyhvlrwGtpiEVEwjoVtqnZnzxvDeNLtf6atCwJzaorIun3/BXHAfDDZ/bp9CgSOCUPkRbqsD5d6NGxgPYabS4tkJKHSAtWWrwfn+3ew+pNO4IORaQOJQ+RFmzM5w8AYPp7mhlaWhYlD5EWbNyhoeRx03ML2fTZ7oCjEaml5CHSgpnVrmBw4t3q9Cgth5KHSAv3xg9PAmDbriq27tCIc2kZlDxEWrg+Xdtx/9eGATC3fEuwwYiEKXmIZIARxV0BuG7S3GADEQlT8hDJAPu1zwdgw/bd3PL8woCjEVHyEMkYT4fX+nh85ip+/tISqvbsDTgiac2UPEQyxPDirow+ZH8A/vj6+zzy5gcBRyStmZKHSAa57axDavbvmLpEM+5KYJQ8RDJIj06FrLxrbM3xyb98PbhgpFVr8cnDzPqZ2cNm9kxE2Sgze8PM7jOzUcFFJxKMJT8dDcCarbsCjkRaq5QmDzN7xMzWmdnCeuWjzWypmS03s4nxruHuK9x9Qv1iYDuhddbLmzdqkZavsE0ulx5XDIC7BxuMtEqprnk8BoyOLDCzXOBeYAwwBLjAzIaY2VAzm1Jv6xHjum+4+xjgR4RWOhRpdXp2KgRg8cfbAo5EWqOUJg93nwFsqlc8AlgerlHsBp4CznL3Be4+rt62LsZ1q/sobgYKop1jZpeZWZmZla1frxlJJfucNqQnAG+t2BhwJNIaBdHm0QtYHXFcHi6Lysy6mdl9wBFmdn247Bwzux/4M/CHaO9z9wfcvdTdS4uKipovepEWon9RB3p1acttU95l2y7NeSXpFUTysChlMR/auvtGd7/c3fu7+53hsmfd/Vvufp67v56qQEVaurHhKdsf/+/KYAORVieI5FEO9Ik47g2sCSAOkYx33ekDAfjlv95jxfrtAUcjrUkQyWM2MMDMSswsHzgfmBxAHCIZryAvlytO6g/AC/M/DjgaaU1S3VX3SWAmMMjMys1sgrtXAVcCLwOLgUnuviiVcYhks6tOGQBA+eadAUcirUleKi/u7hfEKJ8KTE3lZ4u0FgV5ufQras/fylZzxzlDyc2J1qwo0rxa/AhzEWlYj46hHuvlm3cEHIm0FkoeIlngkmNLAFi2djt792rEuaSekodIFugaXizqm38q455XlwUcjbQGSh4iWWDYQfvx4MWldCrMY+2nFUGHI62AkodIFsjNMU4b0pOOhW2oqNIaH5J6Sh4iWaSgTQ67q7Q8raReSrvqikh6FeTl8tqSdZz26+ksWxcacb7o1jNoX6D/1aV5qeYhkkUmHF/CiYOKGNCzQ02ZZt2VVNA/R0SyyLnDenPusN4AFE98IeBoJJup5iEiIklT8hARkaQpeYiISNKUPEREJGlKHiIikjQlDxERSZqSh4iIJE3JQ0REkmbu2T/3v5mtB1Y18u3dgQ3NGE4m070I0X0I0X2ola334iB3L4r2QqtIHk1hZmXuXhp0HC2B7kWI7kOI7kOt1ngv9NhKRESSpuQhIiJJU/Jo2ANBB9CC6F6E6D6E6D7UanX3Qm0eIiKSNNU8REQkaUoeIiKSNCWPOMxstJktNbPlZjYx6HhSzcxWmtkCM5trZmXhsq5m9m8zWxb+uV/E+deH781SMzsjuMibxsweMbN1ZrYwoizp721mw8L3b7mZ3WNmlu7v0lQx7sVPzOyj8O/FXDM7M+K1rLwXZtbHzF4zs8VmtsjMrg6Xt8rfi6jcXVuUDcgF3gf6AfnAPGBI0HGl+DuvBLrXK/sFMDG8PxH4eXh/SPieFAAl4XuVG/R3aOT3PgE4EljYlO8NvA0cAxjwIjAm6O/WTPfiJ8D3o5ybtfcCOAA4MrzfEXgv/H1b5e9FtE01j9hGAMvdfYW77waeAs4KOKYgnAU8Ht5/HPhSRPlT7l7h7h8Aywnds4zj7jOATfWKk/reZnYA0MndZ3roL8afIt6TMWLci1iy9l64+8fu/r/w/jZgMdCLVvp7EY2SR2y9gNURx+XhsmzmwL/MbI6ZXRYu6+nuH0PofyigR7g82+9Pst+7V3i/fnm2uNLM5ocfa1U/qmkV98LMioEjgFno96KGkkds0Z5LZnu/5uPc/UhgDHCFmZ0Q59zWeH8g9vfO5vvxR6A/cDjwMfCrcHnW3wsz6wD8HbjG3T+Nd2qUsqy6F/UpecRWDvSJOO4NrAkolrRw9zXhn+uAfxB6DLU2XPUm/HNd+PRsvz/Jfu/y8H798ozn7mvdfY+77wUepPbxZFbfCzNrQyhxPOHuz4aL9XsRpuQR22xggJmVmFk+cD4wOeCYUsbM2ptZx+p94HRgIaHv/PXwaV8Hng/vTwbON7MCMysBBhBqGMwWSX3v8COMbWZ2dLg3zcUR78lo1X8sw84m9HsBWXwvwnE/DCx2919HvKTfi2pBt9i35A04k1Avi/eBG4OOJ8XftR+h3iLzgEXV3xfoBrwCLAv/7BrxnhvD92YpGdyDBHiS0OOYSkL/UpzQmO8NlBL6w/o+8AfCMzhk0hbjXvwZWADMJ/RH8oBsvxfA8YQeL80H5oa3M1vr70W0TdOTiIhI0vTYSkREkqbkISIiSVPyEBGRpCl5iIhI0pQ8REQkaUoeIgkws+3hn8VmdmEzX/uGesf/bc7ri6SCkodIcoqBpJKHmeU2cEqd5OHuxyYZk0jaKXmIJOcuYGR4XYtrzSzXzO42s9nhiQO/BWBmo8LrQfyV0AA7zOy58KSTi6onnjSzu4C24es9ES6rruVY+NoLw+tBnBdx7dfN7BkzW2JmT1SvEWFmd5nZu+FYfpn2uyOtRl7QAYhkmImE1rYYBxBOAlvdfbiZFQBvmtm/wueOAD7voSm6Ab7h7pvMrC0w28z+7u4TzexKdz88ymedQ2gywsOA7uH3zAi/dgRwCKF5kt4EjjOzdwlNHzLY3d3MujTvVxeppZqHSNOcDlxsZnMJTdndjdC8RhCa2+iDiHOvMrN5wFuEJtEbQHzHA096aFLCtcB0YHjEtcs9NFnhXEKP0z4FdgEPmdk5wI4mfjeRmJQ8RJrGgO+6++HhrcTdq2sen9WcZDYKOBU4xt0PA94BChO4diwVEft7gDx3ryJU2/k7oQWHXkrie4gkRclDJDnbCC1LWu1l4Nvh6bsxs4HhWYnr6wxsdvcdZjYYODritcrq99czAzgv3K5SRGiJ2JgzF4fXnujs7lOBawg98hJJCbV5iCRnPlAVfvz0GPA7Qo+M/hdutF5P9GVGXwIuN7P5hGZdfSvitQeA+Wb2P3e/KKL8H4TWvp5HaIbXH7r7J+HkE01H4HkzKyRUa7m2Ud9QJAGaVVdERJKmx1YiIpI0JQ8REUmakoeIiCRNyUNERJKm5CEiIklT8hARkaQpeYiISNL+Pwk2AGeXReA9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.yscale('log'); plt.xlabel(\"Iterations\"); plt.ylabel(\"Change from last iteration\")\n",
    "plt.plot(range(len(change_iter)), change_iter)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyPEr5R+vBWEUgLKt8e+FhsB",
   "name": "Problem 7",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}