q1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "565c999f-f82b-4927-83c4-c4d10dc6707a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import normal\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import linalg as la\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3aa8c1ab-de2c-41c2-a40e-513ad3da511c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1I0lEQVR4nO3deXRU9f3/8ddkD4GZkEASIgGigBDZFBSmoFSJBEgpAhaXKMFSrTQIQpXlV0QEaygqIt8KuBI8igtFrYUvKAbFogFkR9SAgA2ajYrJEFqyfn5/eDJfRxZhMskk1+fjnDnHe+/nfu778wn9zut7t7EZY4wAAAAsKsDfBQAAANQnwg4AALA0wg4AALA0wg4AALA0wg4AALA0wg4AALA0wg4AALC0IH8X0BjU1NQoPz9fLVq0kM1m83c5AADgPBhjdOLECcXHxysg4Oznbwg7kvLz85WQkODvMgAAgBeOHj2qtm3bnnU7YUdSixYtJH0/WXa73c/VAACA8+FyuZSQkOD+Hj8bwo7kvnRlt9sJOwAANDE/dQsKNygDAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLC/J3AQB+njrMWOuTfr6an+qTfgBYl9/P7HzzzTe67bbbFB0drfDwcHXv3l3bt293bzfGaPbs2WrTpo3Cw8OVnJysgwcPevRx/PhxpaWlyW63KzIyUuPHj1dZWVlDDwUAADRCfg073333nfr376/g4GCtW7dOn332mR5//HG1bNnS3WbBggVavHixli1bpq1btyoiIkIpKSk6deqUu01aWpr279+vDRs2aM2aNfrwww911113+WNIAACgkbEZY4y/Dj5jxgx99NFH+uc//3nG7cYYxcfH649//KPuu+8+SVJpaaliY2OVlZWlm2++WZ9//rmSkpL0ySefqE+fPpKk9evXa9iwYfr6668VHx//k3W4XC45HA6VlpbKbrf7boAAzorLWADq6ny/v/16Zuftt99Wnz599Jvf/EYxMTG6/PLL9eyzz7q3HzlyRIWFhUpOTnavczgc6tu3r3JyciRJOTk5ioyMdAcdSUpOTlZAQIC2bt16xuOWl5fL5XJ5fAAAgDX5NewcPnxYS5cuVadOnfTOO+9owoQJmjRpklasWCFJKiwslCTFxsZ67BcbG+veVlhYqJiYGI/tQUFBioqKcrf5sczMTDkcDvcnISHB10MDAACNhF/DTk1Nja644go98sgjuvzyy3XXXXfpzjvv1LJly+r1uDNnzlRpaan7c/To0Xo9HgAA8B+/hp02bdooKSnJY13Xrl2Vl5cnSYqLi5MkFRUVebQpKipyb4uLi1NxcbHH9qqqKh0/ftzd5sdCQ0Nlt9s9PgAAwJr8Gnb69++v3Nxcj3UHDhxQ+/btJUmJiYmKi4tTdna2e7vL5dLWrVvldDolSU6nUyUlJdqxY4e7zcaNG1VTU6O+ffs2wCgAAEBj5teXCk6ZMkW/+MUv9Mgjj2jMmDHatm2bnnnmGT3zzDOSJJvNpnvvvVcPP/ywOnXqpMTERD3wwAOKj4/XDTfcIOn7M0FDhgxxX/6qrKzUxIkTdfPNN5/Xk1gAAMDa/Bp2rrzySr355puaOXOm5s6dq8TERC1atEhpaWnuNtOmTdPJkyd11113qaSkRAMGDND69esVFhbmbvPyyy9r4sSJGjRokAICAjR69GgtXrzYH0MCAACNjF/fs9NY8J4doOHxnh0AddUk3rMDAABQ3wg7AADA0vjVcwBNGpfDAPwUzuwAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLI+wAAABLC/J3AQCalg4z1vq7BAC4IJzZAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlsbPRQCAfPczGF/NT/VJPwB8hzM7AADA0gg7AADA0gg7AADA0gg7AADA0gg7AADA0gg7AADA0vwadubMmSObzebx6dKli3v7qVOnlJGRoejoaDVv3lyjR49WUVGRRx95eXlKTU1Vs2bNFBMTo/vvv19VVVUNPRQAANBI+f09O5dddpnee+8993JQ0P+VNGXKFK1du1arVq2Sw+HQxIkTNWrUKH300UeSpOrqaqWmpiouLk4ff/yxCgoKNHbsWAUHB+uRRx5p8LEAAIDGx+9hJygoSHFxcaetLy0t1fPPP6+VK1fquuuukyQtX75cXbt21ZYtW9SvXz+9++67+uyzz/Tee+8pNjZWvXr10rx58zR9+nTNmTNHISEhDT0cAADQyPj9np2DBw8qPj5eF198sdLS0pSXlydJ2rFjhyorK5WcnOxu26VLF7Vr1045OTmSpJycHHXv3l2xsbHuNikpKXK5XNq/f/9Zj1leXi6Xy+XxAQAA1uTXsNO3b19lZWVp/fr1Wrp0qY4cOaKrr75aJ06cUGFhoUJCQhQZGemxT2xsrAoLCyVJhYWFHkGndnvttrPJzMyUw+FwfxISEnw7MAAA0Gj49TLW0KFD3f/do0cP9e3bV+3bt9frr7+u8PDwejvuzJkzNXXqVPeyy+Ui8AAAYFF+v4z1Q5GRkercubO+/PJLxcXFqaKiQiUlJR5tioqK3Pf4xMXFnfZ0Vu3yme4DqhUaGiq73e7xAQAA1tSowk5ZWZkOHTqkNm3aqHfv3goODlZ2drZ7e25urvLy8uR0OiVJTqdT+/btU3FxsbvNhg0bZLfblZSU1OD1AwCAxsevl7Huu+8+DR8+XO3bt1d+fr4efPBBBQYG6pZbbpHD4dD48eM1depURUVFyW6365577pHT6VS/fv0kSYMHD1ZSUpJuv/12LViwQIWFhZo1a5YyMjIUGhrqz6EBAIBGwq9h5+uvv9Ytt9yib7/9Vq1bt9aAAQO0ZcsWtW7dWpL0xBNPKCAgQKNHj1Z5eblSUlK0ZMkS9/6BgYFas2aNJkyYIKfTqYiICKWnp2vu3Ln+GhIAAGhkbMYY4+8i/M3lcsnhcKi0tJT7d4Cf0GHGWn+X0Kh9NT/V3yUAPxvn+/3dqO7ZAQAA8DXCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDS/vkEZAKzGVy9d5OWEgO9wZgcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFgaYQcAAFhakL8LANAwOsxY6+8SAMAvOLMDAAAsjbADAAAsjbADAAAsrdGEnfnz58tms+nee+91rzt16pQyMjIUHR2t5s2ba/To0SoqKvLYLy8vT6mpqWrWrJliYmJ0//33q6qqqoGrBwAAjVWjCDuffPKJnn76afXo0cNj/ZQpU/SPf/xDq1at0qZNm5Sfn69Ro0a5t1dXVys1NVUVFRX6+OOPtWLFCmVlZWn27NkNPQQAANBI+T3slJWVKS0tTc8++6xatmzpXl9aWqrnn39eCxcu1HXXXafevXtr+fLl+vjjj7VlyxZJ0rvvvqvPPvtML730knr16qWhQ4dq3rx5euqpp1RRUeGvIQEAgEbE72EnIyNDqampSk5O9li/Y8cOVVZWeqzv0qWL2rVrp5ycHElSTk6OunfvrtjYWHeblJQUuVwu7d+//6zHLC8vl8vl8vgAAABr8ut7dl599VXt3LlTn3zyyWnbCgsLFRISosjISI/1sbGxKiwsdLf5YdCp3V677WwyMzP10EMP1bF6AADQFPjtzM7Ro0c1efJkvfzyywoLC2vQY8+cOVOlpaXuz9GjRxv0+AAAoOH4Lezs2LFDxcXFuuKKKxQUFKSgoCBt2rRJixcvVlBQkGJjY1VRUaGSkhKP/YqKihQXFydJiouLO+3prNrl2jZnEhoaKrvd7vEBAADW5LewM2jQIO3bt0+7d+92f/r06aO0tDT3fwcHBys7O9u9T25urvLy8uR0OiVJTqdT+/btU3FxsbvNhg0bZLfblZSU1OBjAgAAjY/f7tlp0aKFunXr5rEuIiJC0dHR7vXjx4/X1KlTFRUVJbvdrnvuuUdOp1P9+vWTJA0ePFhJSUm6/fbbtWDBAhUWFmrWrFnKyMhQaGhog48JAAA0Po36h0CfeOIJBQQEaPTo0SovL1dKSoqWLFni3h4YGKg1a9ZowoQJcjqdioiIUHp6uubOnevHqgEAQGNiM8YYfxfhby6XSw6HQ6Wlpdy/A8viV8+blq/mp/q7BKDRO9/vb7+/ZwcAAKA+EXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAICleRV2Dh8+7Os6AAAA6oVXYadjx4669tpr9dJLL+nUqVO+rgkAAMBnvAo7O3fuVI8ePTR16lTFxcXp97//vbZt2+br2gAAAOrMq7DTq1cvPfnkk8rPz9cLL7yggoICDRgwQN26ddPChQt17NgxX9cJAADglTrdoBwUFKRRo0Zp1apV+stf/qIvv/xS9913nxISEjR27FgVFBT4qk4AAACv1CnsbN++XX/4wx/Upk0bLVy4UPfdd58OHTqkDRs2KD8/XyNGjPBVnQAAAF4J8manhQsXavny5crNzdWwYcP04osvatiwYQoI+D47JSYmKisrSx06dPBlrQAAABfMq7CzdOlS/fa3v9W4cePUpk2bM7aJiYnR888/X6fiAAAA6sqrsHPw4MGfbBMSEqL09HRvugcAAPAZr+7ZWb58uVatWnXa+lWrVmnFihV1LgoAAMBXvAo7mZmZatWq1WnrY2Ji9Mgjj9S5KAAAAF/xKuzk5eUpMTHxtPXt27dXXl5enYsCAADwFa/CTkxMjPbu3Xva+j179ig6OrrORQEAAPiKV2Hnlltu0aRJk/T++++rurpa1dXV2rhxoyZPnqybb77Z1zUCAAB4zaunsebNm6evvvpKgwYNUlDQ913U1NRo7Nix3LMDAAAaFa/CTkhIiF577TXNmzdPe/bsUXh4uLp376727dv7uj4AAIA68Srs1OrcubM6d+7sq1oA/EiHGWv9XQIANHlehZ3q6mplZWUpOztbxcXFqqmp8di+ceNGnxQHAABQV16FncmTJysrK0upqanq1q2bbDabr+sCAADwCa/CzquvvqrXX39dw4YN83U9AAAAPuXVo+chISHq2LGjr2sBAADwOa/Czh//+Ec9+eSTMsb4uh4AAACf8uoy1ubNm/X+++9r3bp1uuyyyxQcHOyx/Y033vBJcQAAAHXlVdiJjIzUyJEjfV0LAACAz3kVdpYvX+7rOgAAP+Crdyx9NT/VJ/0ATZlX9+xIUlVVld577z09/fTTOnHihCQpPz9fZWVlPisOAACgrrw6s/Ovf/1LQ4YMUV5ensrLy3X99derRYsW+stf/qLy8nItW7bM13UCAAB4xaszO5MnT1afPn303XffKTw83L1+5MiRys7O9llxAAAAdeXVmZ1//vOf+vjjjxUSEuKxvkOHDvrmm298UhgAAIAveHVmp6amRtXV1aet//rrr9WiRYs6FwUAAOArXoWdwYMHa9GiRe5lm82msrIyPfjgg/yEBAAAaFS8uoz1+OOPKyUlRUlJSTp16pRuvfVWHTx4UK1atdIrr7zi6xoBAAC85lXYadu2rfbs2aNXX31Ve/fuVVlZmcaPH6+0tDSPG5YBAAD8zauwI0lBQUG67bbbfFkLAACAz3kVdl588cVzbh87dqxXxQAAAPiaV2Fn8uTJHsuVlZX6z3/+o5CQEDVr1oywAwAAGg2vnsb67rvvPD5lZWXKzc3VgAEDuEEZAAA0Kl7/NtaPderUSfPnzz/trA8AAIA/+SzsSN/ftJyfn+/LLgEAAOrEq3t23n77bY9lY4wKCgr017/+Vf379/dJYQAAAL7g1ZmdG264weMzatQozZkzRz169NALL7xw3v0sXbpUPXr0kN1ul91ul9Pp1Lp169zbT506pYyMDEVHR6t58+YaPXq0ioqKPPrIy8tTamqqmjVrppiYGN1///2qqqryZlgAAMCCvDqzU1NT45ODt23bVvPnz1enTp1kjNGKFSs0YsQI7dq1S5dddpmmTJmitWvXatWqVXI4HJo4caJGjRqljz76SJJUXV2t1NRUxcXF6eOPP1ZBQYHGjh2r4OBgPfLIIz6pEQAANG02Y4zxdxE/FBUVpUcffVQ33nijWrdurZUrV+rGG2+UJH3xxRfq2rWrcnJy1K9fP61bt06/+tWvlJ+fr9jYWEnSsmXLNH36dB07duy0X2WvVV5ervLycveyy+VSQkKCSktLZbfb63+QwHnqMGOtv0tAE/fV/FR/lwDUG5fLJYfD8ZPf316d2Zk6dep5t124cOF5tauurtaqVat08uRJOZ1O7dixQ5WVlUpOTna36dKli9q1a+cOOzk5Oerevbs76EhSSkqKJkyYoP379+vyyy8/47EyMzP10EMPnfcYAABA0+VV2Nm1a5d27dqlyspKXXrppZKkAwcOKDAwUFdccYW7nc1m+8m+9u3bJ6fTqVOnTql58+Z68803lZSUpN27dyskJESRkZEe7WNjY1VYWChJKiws9Ag6tdtrt53NzJkzPQJb7ZkdAABgPV6FneHDh6tFixZasWKFWrZsKen7Fw3ecccduvrqq/XHP/7xvPu69NJLtXv3bpWWlupvf/ub0tPTtWnTJm/KOm+hoaEKDQ2t12MAAIDGwaunsR5//HFlZma6g44ktWzZUg8//LAef/zxC+orJCREHTt2VO/evZWZmamePXvqySefVFxcnCoqKlRSUuLRvqioSHFxcZKkuLi4057Oql2ubQMAAH7evAo7LpdLx44dO239sWPHdOLEiToVVFNTo/LycvXu3VvBwcHKzs52b8vNzVVeXp6cTqckyel0at++fSouLna32bBhg+x2u5KSkupUBwAAsAavLmONHDlSd9xxhx5//HFdddVVkqStW7fq/vvv16hRo867n5kzZ2ro0KFq166dTpw4oZUrV+qDDz7QO++8I4fDofHjx2vq1KmKioqS3W7XPffcI6fTqX79+kmSBg8erKSkJN1+++1asGCBCgsLNWvWLGVkZHCZCgAASPIy7Cxbtkz33Xefbr31VlVWVn7fUVCQxo8fr0cfffS8+ykuLtbYsWNVUFAgh8OhHj166J133tH1118vSXriiScUEBCg0aNHq7y8XCkpKVqyZIl7/8DAQK1Zs0YTJkyQ0+lURESE0tPTNXfuXG+GBQAALKhO79k5efKkDh06JEm65JJLFBER4bPCGtL5PqcPNDTes4O64j07sLLz/f6u0w+BFhQUqKCgQJ06dVJERIQa2fsJAQAAvAs73377rQYNGqTOnTtr2LBhKigokCSNHz/+gh47BwAAqG9ehZ0pU6YoODhYeXl5atasmXv9TTfdpPXr1/usOAAAgLry6gbld999V++8847atm3rsb5Tp07617/+5ZPCAAAAfMGrMzsnT570OKNT6/jx4zzyDQAAGhWvws7VV1+tF1980b1ss9lUU1OjBQsW6Nprr/VZcQAAAHXl1WWsBQsWaNCgQdq+fbsqKio0bdo07d+/X8ePH9dHH33k6xoBAAC85tWZnW7duunAgQMaMGCARowYoZMnT2rUqFHatWuXLrnkEl/XCAAA4LULPrNTWVmpIUOGaNmyZfrTn/5UHzUBAAD4zAWf2QkODtbevXvroxYAAACf8+oy1m233abnn3/e17UAAAD4nFc3KFdVVemFF17Qe++9p969e5/2m1gLFy70SXEAAAB1dUFh5/Dhw+rQoYM+/fRTXXHFFZKkAwcOeLSx2Wy+qw4AAKCOLijsdOrUSQUFBXr//fclff/zEIsXL1ZsbGy9FAcAAFBXF3TPzo9/1XzdunU6efKkTwsCAADwJa9uUK714/ADAADQ2FxQ2LHZbKfdk8M9OgAAoDG7oHt2jDEaN26c+8c+T506pbvvvvu0p7HeeOMN31UIAABQBxcUdtLT0z2Wb7vtNp8WAwAA4GsXFHaWL19eX3UAAADUizrdoAwAANDYefUGZQBA09Bhxlqf9PPV/FSf9AP4A2d2AACApRF2AACApRF2AACApRF2AACApXGDMlAPfHVTKACg7jizAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALM2vYSczM1NXXnmlWrRooZiYGN1www3Kzc31aHPq1CllZGQoOjpazZs31+jRo1VUVOTRJi8vT6mpqWrWrJliYmJ0//33q6qqqiGHAgAAGim/hp1NmzYpIyNDW7Zs0YYNG1RZWanBgwfr5MmT7jZTpkzRP/7xD61atUqbNm1Sfn6+Ro0a5d5eXV2t1NRUVVRU6OOPP9aKFSuUlZWl2bNn+2NIAACgkbEZY4y/i6h17NgxxcTEaNOmTbrmmmtUWlqq1q1ba+XKlbrxxhslSV988YW6du2qnJwc9evXT+vWrdOvfvUr5efnKzY2VpK0bNkyTZ8+XceOHVNISMhpxykvL1d5ebl72eVyKSEhQaWlpbLb7Q0zWFhahxlr/V0C4FNfzU/1dwnAaVwulxwOx09+fzeqe3ZKS0slSVFRUZKkHTt2qLKyUsnJye42Xbp0Ubt27ZSTkyNJysnJUffu3d1BR5JSUlLkcrm0f//+Mx4nMzNTDofD/UlISKivIQEAAD9rNGGnpqZG9957r/r3769u3bpJkgoLCxUSEqLIyEiPtrGxsSosLHS3+WHQqd1eu+1MZs6cqdLSUvfn6NGjPh4NAABoLIL8XUCtjIwMffrpp9q8eXO9Hys0NFShoaH1fhwAAOB/jeLMzsSJE7VmzRq9//77atu2rXt9XFycKioqVFJS4tG+qKhIcXFx7jY/fjqrdrm2DQAA+Pnya9gxxmjixIl68803tXHjRiUmJnps7927t4KDg5Wdne1el5ubq7y8PDmdTkmS0+nUvn37VFxc7G6zYcMG2e12JSUlNcxAAABAo+XXy1gZGRlauXKl/v73v6tFixbue2wcDofCw8PlcDg0fvx4TZ06VVFRUbLb7brnnnvkdDrVr18/SdLgwYOVlJSk22+/XQsWLFBhYaFmzZqljIwMLlUBAAD/hp2lS5dKkn75y196rF++fLnGjRsnSXriiScUEBCg0aNHq7y8XCkpKVqyZIm7bWBgoNasWaMJEybI6XQqIiJC6enpmjt3bkMNAwAANGKN6j07/nK+z+kD54v37MBqeM8OGqMm+Z4dAAAAXyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAAS/PrD4ECAJoGX/3eG7+xBX/gzA4AALA0wg4AALA0wg4AALA0wg4AALA0wg4AALA0wg4AALA0wg4AALA0wg4AALA0wg4AALA03qAM/ICv3hILAGg8OLMDAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsjbADAAAsza9h58MPP9Tw4cMVHx8vm82mt956y2O7MUazZ89WmzZtFB4eruTkZB08eNCjzfHjx5WWlia73a7IyEiNHz9eZWVlDTgKAADQmPk17Jw8eVI9e/bUU089dcbtCxYs0OLFi7Vs2TJt3bpVERERSklJ0alTp9xt0tLStH//fm3YsEFr1qzRhx9+qLvuuquhhgAAABq5IH8efOjQoRo6dOgZtxljtGjRIs2aNUsjRoyQJL344ouKjY3VW2+9pZtvvlmff/651q9fr08++UR9+vSRJP3P//yPhg0bpscee0zx8fENNhYAANA4Ndp7do4cOaLCwkIlJye71zkcDvXt21c5OTmSpJycHEVGRrqDjiQlJycrICBAW7duPWvf5eXlcrlcHh8AAGBNjTbsFBYWSpJiY2M91sfGxrq3FRYWKiYmxmN7UFCQoqKi3G3OJDMzUw6Hw/1JSEjwcfUAAKCx8OtlLH+ZOXOmpk6d6l52uVwEniauw4y1/i4BANBINdqwExcXJ0kqKipSmzZt3OuLiorUq1cvd5vi4mKP/aqqqnT8+HH3/mcSGhqq0NBQ3xcNADgnX/0/Jl/NT/VJP/h5aLSXsRITExUXF6fs7Gz3OpfLpa1bt8rpdEqSnE6nSkpKtGPHDnebjRs3qqamRn379m3wmgEAQOPj1zM7ZWVl+vLLL93LR44c0e7duxUVFaV27drp3nvv1cMPP6xOnTopMTFRDzzwgOLj43XDDTdIkrp27aohQ4bozjvv1LJly1RZWamJEyfq5ptv5kksAAAgyc9hZ/v27br22mvdy7X30aSnpysrK0vTpk3TyZMnddddd6mkpEQDBgzQ+vXrFRYW5t7n5Zdf1sSJEzVo0CAFBARo9OjRWrx4cYOPBQAANE42Y4zxdxH+5nK55HA4VFpaKrvd7u9y4AVuUAZ+XrhnB9L5f3832nt2AAAAfIGwAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALI2wAwAALC3I3wXg563DjLX+LgEAYHGc2QEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJbGz0UAAJocX/3UzFfzU33SDxo3wg4A4GeL0PTzwGUsAABgaYQdAABgaYQdAABgadyzA6/46jo3AAD1jTM7AADA0gg7AADA0gg7AADA0gg7AADA0gg7AADA0gg7AADA0gg7AADA0njPzs8M78cBAPzccGYHAABYGmd2AACoI349vXGzzJmdp556Sh06dFBYWJj69u2rbdu2+bskAADQCFjizM5rr72mqVOnatmyZerbt68WLVqklJQU5ebmKiYmxq+1kfYBAOeL74z6YYmws3DhQt1555264447JEnLli3T2rVr9cILL2jGjBl+rg4AgIZFaPLU5MNORUWFduzYoZkzZ7rXBQQEKDk5WTk5OWfcp7y8XOXl5e7l0tJSSZLL5fJ5fTXl//FJP+2mrPJJPwAAnC9fffd8+lCKT/r5sdrvbWPMOds1+bDz73//W9XV1YqNjfVYHxsbqy+++OKM+2RmZuqhhx46bX1CQkK91AgAwM+ZY1H99n/ixAk5HI6zbm/yYccbM2fO1NSpU93LNTU1On78uKKjo2Wz2fxYWcNyuVxKSEjQ0aNHZbfb/V1Ok8CcXTjm7MIxZxeOObtwVpgzY4xOnDih+Pj4c7Zr8mGnVatWCgwMVFFRkcf6oqIixcXFnXGf0NBQhYaGeqyLjIysrxIbPbvd3mT/ofsLc3bhmLMLx5xdOObswjX1OTvXGZ1aTf7R85CQEPXu3VvZ2dnudTU1NcrOzpbT6fRjZQAAoDFo8md2JGnq1KlKT09Xnz59dNVVV2nRokU6efKk++ksAADw82WJsHPTTTfp2LFjmj17tgoLC9WrVy+tX7/+tJuW4Sk0NFQPPvjgaZf0cHbM2YVjzi4cc3bhmLML93OaM5v5qee1AAAAmrAmf88OAADAuRB2AACApRF2AACApRF2AACApRF2LCAzM1NXXnmlWrRooZiYGN1www3Kzc095z6VlZWaO3euLrnkEoWFhalnz55av359nfttKuprzn5o/vz5stlsuvfee31cvX/U55x98803uu222xQdHa3w8HB1795d27dvr6+hNJj6mrPq6mo98MADSkxMVHh4uC655BLNmzfvJ38fqKlYunSpevTo4X7ZndPp1Lp16865z6pVq9SlSxeFhYWpe/fu+t///V+P7cYYzZ49W23atFF4eLiSk5N18ODB+hxGg/L1nFVWVmr69Onq3r27IiIiFB8fr7Fjxyo/P7++h1I/DJq8lJQUs3z5cvPpp5+a3bt3m2HDhpl27dqZsrKys+4zbdo0Ex8fb9auXWsOHTpklixZYsLCwszOnTvr1G9TUV9zVmvbtm2mQ4cOpkePHmby5Mn1OJKGU19zdvz4cdO+fXszbtw4s3XrVnP48GHzzjvvmC+//LIhhlWv6mvO/vznP5vo6GizZs0ac+TIEbNq1SrTvHlz8+STTzbEsOrd22+/bdauXWsOHDhgcnNzzf/7f//PBAcHm08//fSM7T/66CMTGBhoFixYYD777DMza9YsExwcbPbt2+duM3/+fONwOMxbb71l9uzZY37961+bxMRE89///rehhlWvfD1nJSUlJjk52bz22mvmiy++MDk5Oeaqq64yvXv3bshh+Qxhx4KKi4uNJLNp06aztmnTpo3561//6rFu1KhRJi0trU79NlW+nLMTJ06YTp06mQ0bNpiBAwdaJuz8mK/mbPr06WbAgAH1Vmdj4qs5S01NNb/97W/P2cZqWrZsaZ577rkzbhszZoxJTU31WNe3b1/z+9//3hhjTE1NjYmLizOPPvqoe3tJSYkJDQ01r7zySv0V7Wd1mbMz2bZtm5Fk/vWvf/m0zobAZSwLKi0tlSRFRUWdtU15ebnCwsI81oWHh2vz5s116rep8uWcZWRkKDU1VcnJyb4vtBHx1Zy9/fbb6tOnj37zm98oJiZGl19+uZ599tn6KdrPfDVnv/jFL5Sdna0DBw5Ikvbs2aPNmzdr6NCh9VC1f1VXV+vVV1/VyZMnz/oTQDk5Oaf97y0lJUU5OTmSpCNHjqiwsNCjjcPhUN++fd1trMQXc3YmpaWlstlsTfO3JP2dtuBb1dXVJjU11fTv3/+c7W655RaTlJRkDhw4YKqrq827775rwsPDTUhISJ36bYp8OWevvPKK6datm/vUuFXP7PhyzkJDQ01oaKiZOXOm2blzp3n66adNWFiYycrKqu9hNChfzll1dbWZPn26sdlsJigoyNhsNvPII4/U9xAa1N69e01ERIQJDAw0DofDrF279qxtg4ODzcqVKz3WPfXUUyYmJsYY8/0lG0kmPz/fo81vfvMbM2bMGN8X7ye+nLMf++9//2uuuOIKc+utt/q05oZC2LGYu+++27Rv394cPXr0nO2Ki4vNiBEjTEBAgAkMDDSdO3c2f/jDH0xYWFid+m2KfDVneXl5JiYmxuzZs8e9j1XDji//nQUHBxun0+mx3z333GP69etXL7X7iy/n7JVXXjFt27Y1r7zyitm7d6958cUXTVRUlKUCYnl5uTl48KDZvn27mTFjhmnVqpXZv3//GdsSdr7nyzn7oYqKCjN8+HBz+eWXm9LS0nqpvb4RdiwkIyPDtG3b1hw+fPi89/nvf/9rvv76a1NTU2OmTZtmkpKSfNJvU+HLOXvzzTeNJBMYGOj+SDI2m80EBgaaqqqq+hpGg/L1v7N27dqZ8ePHe7RfsmSJiY+P91nN/ubrOWvbtu1p9/XMmzfPXHrppT6rubEZNGiQueuuu864LSEhwTzxxBMe62bPnm169OhhjDHm0KFDRpLZtWuXR5trrrnGTJo0qT7KbRTqMme1KioqzA033GB69Ohh/v3vf9dXqfWOe3YswBijiRMn6s0339TGjRuVmJh43vuGhYXpoosuUlVVlVavXq0RI0b4pN/Grj7mbNCgQdq3b592797t/vTp00dpaWnavXu3AgMD62s4DaK+/p3179//tMexDxw4oPbt2/usdn+przn7z3/+o4AAz//zHRgYqJqaGp/V3tjU1NSovLz8jNucTqeys7M91m3YsMF9v0piYqLi4uI82rhcLm3duvWs97RYQV3mTPr+8fMxY8bo4MGDeu+99xQdHV2v9dYrv0Yt+MSECROMw+EwH3zwgSkoKHB//vOf/7jb3H777WbGjBnu5S1btpjVq1ebQ4cOmQ8//NBcd911JjEx0Xz33XcX1G9TVV9z9mNWuoxVX3O2bds2ExQUZP785z+bgwcPmpdfftk0a9bMvPTSSw05vHpRX3OWnp5uLrroIvej52+88YZp1aqVmTZtWkMOr97MmDHDbNq0yRw5csTs3bvXzJgxw9hsNvPuu+8aY06fs48++sgEBQWZxx57zHz++efmwQcfPOOj55GRkebvf/+72bt3rxkxYoSlHj339ZxVVFSYX//616Zt27Zm9+7dHv9+y8vL/TLGuiDsWICkM36WL1/ubjNw4ECTnp7uXv7ggw9M165dTWhoqImOjja33367+eabby6436aqvubsx6wUdupzzv7xj3+Ybt26mdDQUNOlSxfzzDPPNMCI6l99zZnL5TKTJ0827dq1M2FhYebiiy82f/rTn5rkl9CZ/Pa3vzXt27c3ISEhpnXr1mbQoEHuL21jTp8zY4x5/fXXTefOnU1ISIi57LLLTrs5t6amxjzwwAMmNjbWhIaGmkGDBpnc3NyGGE6D8PWcHTly5Kz/ft9///0GGpXv2IyxyCs3AQAAzoB7dgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgA0qDlz5qhXr17+LgPAzwhhB2hCxo0bJ5vNdtpnyJAh593HBx98IJvNppKSkvor1CKysrIUGRl5Xu3O9HcJCwu7oOPZbDa99dZb3hUL4KyC/F0AgAszZMgQLV++3GNdaGion6pBLbvdftqvt9tsNp8fp6KiQiEhIT7vF7AyzuwATUxoaKji4uI8Pi1btnRvt9lseu655zRy5Eg1a9ZMnTp10ttvvy1J+uqrr3TttddKklq2bCmbzaZx48ZJktavX68BAwYoMjJS0dHR+tWvfqVDhw65+/3qq69ks9n0xhtv6Nprr1WzZs3Us2dP5eTkeNT37LPPKiEhQc2aNdPIkSO1cOHCnzw78txzz6lr164KCwtTly5dtGTJktOO+/rrr+vqq69WeHi4rrzySh04cECffPKJ+vTpo+bNm2vo0KE6duzYBfd7tvF88MEHuuOOO1RaWuo+UzNnzpyzjsFms532d4mNjXVv/+Uvf6lJkyZp2rRpioqKUlxcnEd/HTp0kCSNHDlSNpvNvVx72e+5555TYmKi+2xRXl6eRowYoebNm8tut2vMmDEqKipy91e739NPP+3+e4wZM0alpaWSpA8//FDBwcEqLCz0GMe9996rq6+++px/L6DJ8fcvkQI4f+np6WbEiBHnbCPJtG3b1qxcudIcPHjQTJo0yTRv3tx8++23pqqqyqxevdpIMrm5uaagoMCUlJQYY4z529/+ZlavXm0OHjxodu3aZYYPH266d+9uqqurjTH/9yvIXbp0MWvWrDG5ubnmxhtvNO3btzeVlZXGGGM2b95sAgICzKOPPmpyc3PNU089ZaKioozD4XDX9+CDD5qePXu6l1966SXTpk0bs3r1anP48GGzevVqExUVZbKysk477vr1681nn31m+vXrZ3r37m1++ctfms2bN5udO3eajh07mrvvvturfs80nvLycrNo0SJjt9tNQUGBKSgoMCdOnDjjnC9fvtxjjGcycOBAY7fbzZw5c8yBAwfMihUrjM1mc/8ydXFxsfsX0QsKCkxxcbF7viIiIsyQIUPMzp07zZ49e0x1dbXp1auXGTBggNm+fbvZsmWL6d27txk4cKDHPEdERJjrrrvO7Nq1y2zatMl07NjR3Hrrre42nTt3NgsWLHAvV1RUmFatWpkXXnjhnGMBmhrCDtCEpKenm8DAQBMREeHx+fOf/+xuI8nMmjXLvVxWVmYkmXXr1hljjHn//feNJPPdd9+d81jHjh0zksy+ffuMMf8XDp577jl3m/379xtJ5vPPPzfGGHPTTTeZ1NRUj37S0tLOGXYuueQSs3LlSo995s2bZ5xO51mP+8orrxhJJjs7270uMzPTXHrppXXq98fjOZ8QU9tO0ml/lyFDhrjbDBw40AwYMMBjvyuvvNJMnz7dvSzJvPnmmx5tHnzwQRMcHOwOP8YY8+6775rAwECTl5d3Wu3btm1z7xcYGGi+/vprd5t169aZgIAAU1BQYIwx5i9/+Yvp2rWre/vq1atN8+bNTVlZ2U+OGWhKuGcHaGKuvfZaLV261GNdVFSUx3KPHj3c/x0RESG73a7i4uJz9nvw4EHNnj1bW7du1b///W/V1NRI+v5ySbdu3c7Yd5s2bSRJxcXF6tKli3JzczVy5EiPfq+66iqtWbPmjMc8efKkDh06pPHjx+vOO+90r6+qqpLD4TjrmGovD3Xv3t1jXe0Yve33x+O5EC1atNDOnTs91oWHh5/1WLXH+6m/iyS1b99erVu3di9//vnnSkhIUEJCgntdUlKSIiMj9fnnn+vKK6+UJLVr104XXXSRu43T6VRNTY1yc3MVFxencePGadasWdqyZYv69eunrKwsjRkzRhEREec/cKAJIOwATUxERIQ6dux4zjbBwcEeyzabzR1ezmb48OFq3769nn32WcXHx6umpkbdunVTRUXFWfuuvQH3p/o+m7KyMknf3+fTt29fj22BgYE/edwfr6uto679ejOegICAevm7SKq38BETE6Phw4dr+fLlSkxM1Lp16/TBBx/Uy7EAfyLsAD8ztU/yVFdXu9d9++23ys3N1bPPPuu+OXXz5s0X3Pell16qTz75xGPdj5d/KDY2VvHx8Tp8+LDS0tIu+Hj13W9ISIjHPNW34ODg8zpe165ddfToUR09etR9duezzz5TSUmJkpKS3O3y8vKUn5+v+Ph4SdKWLVsUEBCgSy+91N3md7/7nW655Ra1bdtWl1xyifr37+/jUQH+R9gBmpjy8vLTnqAJCgpSq1atzmv/9u3by2azac2aNRo2bJjCw8PVsmVLRUdH65lnnlGbNm2Ul5enGTNmXHBt99xzj6655hotXLhQw4cP18aNG7Vu3bpzPoL90EMPadKkSXI4HBoyZIjKy8u1fft2fffdd5o6deoF1+DLfjt06KCysjJlZ2erZ8+eatasmZo1a3bGtsaY0/4u0vdnTwICzu/B1w4dOig7O1v9+/dXaGiox1N2P5ScnKzu3bsrLS1NixYtUlVVlf7whz9o4MCB6tOnj7tdWFiY0tPT9dhjj8nlcmnSpEkaM2aM4uLi3G1SUlJkt9v18MMPa+7cuedVJ9DU8Og50MSsX79ebdq08fgMGDDgvPe/6KKL9NBDD2nGjBmKjY3VxIkTFRAQoFdffVU7duxQt27dNGXKFD366KMXXFv//v21bNkyLVy4UD179tT69es1ZcqUc75c73e/+52ee+45LV++XN27d9fAgQOVlZWlxMTECz6+r/v9xS9+obvvvls33XSTWrdurQULFpy1rcvlOu3vcr735NR6/PHHtWHDBiUkJOjyyy8/azubzaa///3vatmypa655holJyfr4osv1muvvebRrmPHjho1apSGDRumwYMHq0ePHh6P30vfX34bN26cqqurNXbs2POuFWhKbMYY4+8iAFjXnXfeqS+++EL//Oc//V3Kz8qcOXP01ltvaffu3T/Zdvz48Tp27Jj7fUyA1XAZC4BPPfbYY7r++usVERGhdevWacWKFaedTUDjUFpaqn379mnlypUEHVgaYQeAT23btk0LFizQiRMndPHFF2vx4sX63e9+5++ycAYjRozQtm3bdPfdd+v666/3dzlAveEyFgAAsDRuUAYAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJb2/wHWznGQkPk9+AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def haar_state(dim=2):\n",
    "    state = np.array([normal()+1j*normal() for _ in range(dim)])\n",
    "    return state/la.norm(state)\n",
    "\n",
    "def half_entropy(state):\n",
    "    dim = len(state)\n",
    "    dimA = int(np.sqrt(dim))\n",
    "    assert dimA == np.sqrt(dim)\n",
    "    S = la.svd(state.reshape((dimA,dimA)),compute_uv=False)\n",
    "    return sum(-np.log(S**2)*S**2)\n",
    "\n",
    "def test_entanglement(samples,dim=2**10):\n",
    "    results = [half_entropy(haar_state(dim)) for _ in range(samples)]\n",
    "    return results\n",
    "\n",
    "\n",
    "res=test_entanglement(5000)\n",
    "plt.hist(res,bins=25)\n",
    "plt.xlabel(\"Entanglement Entropy\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1016a7e4-6951-4b5f-b302-222c90e21121",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}