{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "cb5df239",
      "metadata": {},
      "source": [
        "# Orthogonal and QR Iteration\n",
        "\n",
        "Follow a real $8\\times8$ matrix through orthogonal iteration, Hessenberg reduction, and **complex Wilkinson-shifted QR** to an upper-triangular Schur form. Cells marked `#clear` contain short algorithms to type in class; surrounding instructions are supplied."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "qr-demo-01",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Eigenvalues: [-1. +0.j  0.5+0.j  2. +0.j  3. -2.j  3. +2.j  4. +0.j  6. +0.j  8. +0.j]\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "import matplotlib.pyplot as pt\n",
        "\n",
        "np.set_printoptions(precision=3, suppress=True, linewidth=100)\n",
        "rng = np.random.default_rng(3)\n",
        "n = 8\n",
        "T0 = np.diag([8., 6., 4., 3., 3., 2., -1., .5])\n",
        "T0 += .35 * np.triu(rng.normal(size=(n, n)), 1)\n",
        "T0[3, 4], T0[4, 3] = -2., 2.             # One block with eigenvalues 3 +/- 2i.\n",
        "U, _ = la.qr(rng.normal(size=(n, n)))\n",
        "A = U @ T0 @ U.T\n",
        "print(\"Eigenvalues:\", np.sort_complex(la.eigvals(A)))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "qr-demo-02",
      "metadata": {},
      "source": [
        "The matrix is dense and nonsymmetric, with six real eigenvalues and one conjugate pair.\n",
        "The supplied plotter shows the numerical fill pattern (white means $|M_{ij}| \\le 10^{-12}\\|A\\|_F$)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "qr-demo-03",
      "metadata": {},
      "outputs": [],
      "source": [
        "def show_fill(frames):\n",
        "    fig, axes = pt.subplots(1, len(frames), figsize=(2.5 * len(frames), 2.7), squeeze=False)\n",
        "    for ax, (M, title) in zip(axes[0], frames):\n",
        "        ax.imshow(np.abs(M) > 1e-12 * la.norm(A), cmap=\"Greys\", vmin=0, vmax=1)\n",
        "        ax.set(title=title, xticks=range(n), yticks=range(n))\n",
        "    pt.tight_layout()\n",
        "    pt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "qr-demo-04",
      "metadata": {},
      "source": [
        "## 1. Orthogonal iteration and QR iteration\n",
        "\n",
        "Orthogonal iteration: $AQ_k=Q_{k+1}S_{k+1}$, starting with $Q_0=I$.\n",
        "QR iteration: $B_k=V_kR_k$, $B_{k+1}=R_kV_k$, starting with $B_0=A$.\n",
        "Track $W_{k+1}=W_kV_k$ so $B_k=W_k^TAW_k$.\n",
        "\n",
        "**Sign convention:** NumPy's `la.qr` does *not* require positive diagonal entries in $R$.\n",
        "The two orthonormal bases agree up to column signs: $W_k=Q_kD_k$, where\n",
        "$D_k$ is diagonal with entries $\\pm1$. Thus\n",
        "$B_k=D_k(Q_k^TAQ_k)D_k$; we adjust signs **only for the comparison**, not for either algorithm."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "qr-demo-05",
      "metadata": {},
      "outputs": [],
      "source": [
        "Q_orth, B, W, k = np.eye(n), A.copy(), np.eye(n), 0\n",
        "# Orthogonal iteration: QR-factor A @ Q_orth; keep its orthogonal factor.\n",
        "# QR iteration: QR-factor B, form R @ V; also retain V to track the similarity."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "qr-demo-06",
      "metadata": {},
      "outputs": [],
      "source": [
        "def orthogonal_iteration(A, Q):\n",
        "    Q, _ = la.qr(A @ Q)\n",
        "    return Q\n",
        "\n",
        "def qr_iteration(B):\n",
        "    V, R = la.qr(B)\n",
        "    return R @ V, V"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "qr-demo-07",
      "metadata": {},
      "source": [
        "Repeat the next cell with **Ctrl-Enter** to advance both algorithms.\n",
        "Compare the raw matrices and the comparison adjusted for column signs. Why can the first difference be large?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "qr-demo-08",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Step 1: raw difference 1.81e-16, after sign alignment 1.81e-16\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPAAAAEKCAYAAAAsK9hZAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIFJJREFUeJzt3XlUk1f+BvAHCIZFEEERFxBM3UBFrRitW1Fbp+5UrR7r1sVjp62nU22n2trjTHVcW9fROR2XurZuU3Fp0eNutQhKcRurIFpAUBF1ArKFkO/vD39kjGwJvkHfzvM5555jbt733pvgk3fNjZOICIhIlZyf9gCIqPoYYCIVY4CJVIwBJlIxBphIxRhgIhVjgIlUjAEmUjEG2A5JSUk4d+6cTcumpKQgMTGx3OdycnJw9uxZnDp1CiaTSckhlnHhwgVcvnzZoX3Q0+OkpjuxLl26hIKCAjz//PPlPn/69GnUrl0brVu3dkj/Q4cOxeXLl20KxIQJE3Dw4EHcuHHDqv6vf/0r5s2bh1atWkGr1WL//v1ISkpC/fr1ERwcrPiY27dvj4CAAOzbt8/mdQoLC3H9+nUYjUY0adIEfn5+5S539uxZFBYWAgCcnJzg7u6OkJAQeHl5KTL2ymRmZiItLQ0RERFwcXFxeH9VSU9PR0ZGBoKCgtCoUaOa61hUpFevXuLn51fh856entKnTx+H9T9kyBBp2bKlTct+8cUXMmjQIKu6y5cvi5OTk3zzzTdW9Z6envLee+8pNUwr4eHh0q9fP5uWvXHjhowZM0bc3d1Fp9NJu3btxN3dXbp06SJHjhwps7xOpxMPDw/R6/Wi1+ulWbNm4uLiIlFRUZKVlaXwK7E2d+5cASD37993aD+VycvLk0WLFolerxc/Pz8BIHPnzq3RMWhq7qPif8vnn39epu78+fMQkQr3IJ6mlJQU9OjRA/Xr18fp06cRFhYG4OHu/pQpU9C3b1+sW7cOY8aMsVqvefPmOHXqlOVxTEwMBg4ciMLCQvz44481+hpq2u3bt5Geno4lS5ZAo9EgIiKixsfwuw5wUlISCgoKEB4eDpPJhCtXrkCr1eK5554rd/mSkhKkp6cjNzcXOp0OHh4eFbZdVXspKSnIyclBhw4dAADnzp2zHBNfunQJeXl5qFOnDgwGA8xmM27fvm0JQp06dcocBuTn5+Pq1atwcXFBy5YtodGU/6e7fv068vLy0KJFC9SqVavqN+n/jRkzBgUFBdi/fz8CAgIs9d7e3li9ejUyMjIwceJE9OzZE0FBQRW288orr6BXr17Yv38/cnJy4O3tDQAwm82Ij49H3bp10bJlS5vGlJGRgezsbDRt2hQ+Pj5WrzE9PR0AcObMGdSuXRsA0KZNG8u/S/3222+4d+8emjRpAn9/f6vnbty4gRs3bkCv1wMAkpOTYTKZ0KpVKzg7V316KCQkBIsWLbKM46mo0e39E7J3F7p0lzchIUFat24toaGhotVqpWPHjnLz5k2rdbdv3y6NGzeWgIAACQ8PF29vbxk1apTcvn27Wu2NHz9eGjdubHk8cuRIadasmQCQ8PBw0ev18uGHH4perxdnZ2dp0KCBZVd08uTJlvVyc3Nl4sSJotVqpVmzZtKkSRPx8/OTtWvXWvV3/fp16dy5s2g0GmndurWEhITI4cOHbdqFjo2NFQDy0UcfVbjMTz/9JABk+vTpljqdTifh4eFllh0+fLgAkBs3bli9DgAyYMCASsciIhIfHy9hYWHi4+Mj4eHhUq9ePenVq5dcvHhRRETmz58vgYGBAkA6depked9KnxcR+eGHH6R58+bi5eUlYWFh4ubmJv3797fatZ85c6YAkF9//VU6duworVq1Ek9PT2natKkcP368ynE+6vTp009lF/p3H+BGjRrJmDFjLMdKSUlJ4uPjI2+99ZZlubt370qtWrVk2rRpYjabRUTEbDbL1q1bJT4+3u72RMoGWERk1apVlv8wj4+7omPgvn37SoMGDazGsXr1agEge/bsERERo9EooaGh0qpVK0lLSxMRkZycHBk2bJgEBQVVGeD58+dbtVee4uJi0Wq10q1bN0tdeQE2Go0SHBwsfn5+YjKZLPX5+fmi1+tl6tSplY6ltN2BAwdKUVGRpe6nn36SnTt3Wh5Xdgx88OBBcXFxkffff9/Sxs2bN6V9+/bStWtXy9+4NMBRUVFy9epVEXl4XNu/f3/x8vKyvJe2YIBtUJ0AA5CkpCSr5SZOnCh16tSxPD5z5owAkO3bt1fav63tiSgT4AMHDgiAMltbkYfvRffu3UVE5F//+pcAkOjoaKtlSk+aVRXgKVOmCACrD4nyNGrUSHQ6neWxTqeT5s2bS2xsrMTGxsquXbukf//+4uzsLOvXr6+0rYoYjUZxcnKSzz77rNLlKguwXq+X5s2bW32AiPz3/Tx69KiI/DfAmzZtslouNTVVNBqNfPLJJzaP+2kF+Hd9DAwA9erVQ/Pmza3qmjVrBoPBYDlGCw0NRUhICCZOnIjY2Fi8/PLL6NatW5njKVvbU8rRo0cBAH5+fjhz5gzk4QcuAKBhw4bYtWsXACA+Ph4A0L17d6v1W7ZsWea4rzxarRYAYDQaK12uqKioTHsZGRn405/+BBFBWloabt26hc8//xzjxo2r+gWWw9XVFf369cPChQuRmZmJAQMGoGfPnqhfv75N6+fn5yM+Ph6DBg1CYmKi5f0SEZjNZgBAQkICevXqZVnn8fctKCgITZs2RVxcXLVeQ01SVYBdXV0rvfGhuLgYrq6uVnXlXcN0c3MD8PCP7e3tDXd3d8TFxWHp0qWIiYnB0qVL4ezsjKioKKxYsQL16tWzqz2l5OTkAAD+9re/wcnJqczz7dq1g9lsRl5eHgCU23edOnWq7KdFixYAgCtXrqBbt27lLnP37l3cvXsXPXr0sKp/9Cy0iGDu3Ln47LPPEBgYiIkTJ1bZd3l27tyJlStXYteuXRg7diwKCwvRs2dP/OMf/6jyGv+DBw8gIkhISMD7779f5nm9Xl/mPSnvPapTp47lfX2WqSrAjRo1Qk5OTrlbuqysLBiNRjRu3LhabdevXx+zZ8/G7NmzkZubi+3bt+Odd95B7dq1sWbNGiWGX6HywgnAcrZ3w4YNlZ65bdKkCQAgLS0NOp3OUl9SUoIbN24gJCSk0v5feeUVaLVabNu2DW+++Wa5y2zduhUA8Nprr1X6Oj799FMcOHAAU6ZMwYABA6p1U4ObmxumTJmCKVOmwGg04uDBg5gwYQLGjBmDhIQES1/l8fX1haenJyIjI7Fx40ab+ktNTbU6y126N/H4h9WzSFW3Ug4aNAgigm3btpV57rvvvgMADB482O528/LyLLtaAODl5YU333wTOp0OKSkp1R+wjXx8fMr9tB8xYgRq1aqF5cuXl7tebm4ugIfvi5OTE9atW2f1/NatW1FUVFRl/w0aNMCf//xn7N+/H5s2bSrz/PXr1/GXv/wF7du3x7Bhw6psb+HChXjw4AG++OILS53ZbMapU6dw5cqVStc1Go1WY65Vqxb69++PPn36WP0tSgP3+Pum0WgwatQo7N69G2lpaWXaLy4uLvOePP6a9+zZg+zsbAwZMqTyF/oMUNUWeNiwYRg1ahTee+89XL16FT179oSI4PDhw1i2bBnGjx9frQDHxsbigw8+wPjx49GmTRtotVrExMTgypUr+PTTTx3wSqx1794d+/fvR3R0NAICAizXgZs2bYo1a9bgrbfewu3bt/Haa6/Bz88P169fxw8//AA/Pz+sWrUKoaGhmDp1KubMmYOSkhJERkbiwoUL+Pnnn9GmTRubxjBz5kxkZ2djwoQJOH78OIYOHQoPDw/Ex8dj4cKFaNiwIfbs2WPTteVOnTph6NChWLNmDT7++GPodDrk5+eja9euGDBgAPbu3Vvhunfv3kXXrl0xduxYdOzYEXXr1sWZM2ewa9cuvPvuu5blXnjhBTg7O2PevHl47bXX4OrqarkO/OWXX+L8+fPo2rUrPvroI7Rt2xa5ubk4d+4cNmzYgMOHD1vdtqrRaDB16lT069cPV65cwYwZM9C3b98yN62U58yZMzCZTJbba9PT0y2HFG3btoWnp2eVbTwJVQXYyckJ3333HaKjo7Fr1y4sWLAAwMML6t9//z0GDRpktXzLli0tx6ePatiwIfR6veU/Y9++fREdHY21a9dixYoVMJlM0Ol0OHXqFDp37mx3ewCg0+lw7949q+X8/f2h1+vh7u5uVb9s2TLMmTMHS5cuRWFhISIiIrBs2TIAD2+wiIiIwJo1a7Bq1SqUlJRAp9NhwoQJVq934cKFaN++PbZs2YK4uDj06NEDGzduxKRJkyq8l/lRLi4uWLlyJSZMmIBvv/0WS5YswX/+8x+cPn0aHTt2RFxcXJmbRzp06FDmtZSaNWsWbt68ie+//x4ff/wxXFxcoNfr0apVq0rH0bBhQ8THx2PNmjXYuHEjDAYDAgMDsW3bNgwcONCyXNu2bbFlyxZ8++23mDZtGkpKSrBmzRqEhYXBx8cHJ06cwObNm7Fv3z7s3LkT/v7+aN++PU6ePImGDRta9Tlz5kxs3LgRixcvhslkwueff47JkyfbdI/19OnTLXtCer0eCQkJlt38devWVfl6n1iNnvMm1Zk1a5YAkBkzZjztoSiu9DJSQUHB0x5KtalqC0w1b8aMGcjPz8fBgwfRo0cPvPzyy097SPQIBpiqNGfOnKc9BKqAqs5CEympSZMm0Ov1Nn1x4Vmlqi/0E5E19X70EBEDTKRmNX4Sy2w2IzMzE15eXhXeDkf0v0xEkJubi0aNGlV5fF7jAc7MzERgYGBNd0ukOunp6Zb73CtS47vQNTFjIdHvgS1ZqfEAc7eZyDa2ZIUnsYhUjAEmUjEGmEjFGGAiFat2gEsnCCOip8fuAK9atQpBQUFwdXXFc889hx07djhiXERkC3u+PLx7925xdXWVrVu3Sl5enqxYsUJcXFzk1KlTNrdhMBgEAAsLSxXFYDBUmSe7vo3Uu3dv+Pr6Wm11u3TpAp1Oh82bN9vURk5Ojk1TnRL9rzMYDFVOU2zzLrSIIC4uzmpCbACIjIxEbGxs9UZIRE/E5nuhc3NzkZ+fX2aGfH9/f9y+fbvC9YqKiqym8SydrJyInpzdJ7EeP/tsNpsrveVr7ty5qFOnjqXwiwxEyrE5wF5eXvD09ERWVpZVfVZWltXvyT5u+vTpMBgMllL6u65E9ORsDrCTkxO6deuGI0eOWNUfOnSowt/TAR7+cJa3t7dVISKF2HMZ6cCBA+Li4iJff/21ZGZmypw5c8TV1VXOnj3Ly0gsLAoXWy4j2T2x+5YtW6RNmzbi7e0tzz//vOzbt8+u9RlgFhbbiuLXgZXA68BEtlH0OjARPXsYYCIVY4CJVIwBJlIxBphIxRhgIhVjgIlUjAEmUjEGmEjFGGAiFWOAiVSMASZSMQaYSMUYYCIVsyvAIoIDBw7g1VdfRXBwMKKjox00LCKyhV0BXrJkCebPn4/Ro0cjNTUVDx48cNS4iMgGNk8rCwCTJ0/Ghx9+6KixEJGd7NoCazR25Z2IHMzhieTE7kSO4/Cz0JzYnchxHB5gTuxO5DgO34XWarXQarWO7obofxJv5CBSMbsCfPLkSQQHByM4OBgA8OGHHyI4OBjTpk1zxNiIqAp2TexeWFiIW7dulan38vKCn5+fTW1wYnci29gysbtdx8Bubm6WrS8RPX08BiZSMQaYSMUYYCIVY4CJVIwBJlIxBphIxRhgIhVjgIlUjAEmUjEGmEjFGGAiFWOAiVSMASZSMbu+jWQ0GrF582YcO3YMxcXFiIiIwKRJk+Du7u6o8RFRJez6PnCnTp0QHh6OHj16wMXFBV9++SVcXFxw8uRJm0PM7wMT2caW7wND7HDnzh2rx5mZmeLk5CQ7duywuQ2DwSAAWFhYqigGg6HKPNl1DFyvXj2rx3Xr1oVGo+FPrBA9JU80K+WyZcvg7OyMPn36VLgMJ3YnciB7dqEfFRMTIxqNRr7++utKl5s5c+ZT3xVhYVFjsWUXuloBPnTokLi7u8u8efOqXLawsFAMBoOlpKenP/U3hoVFDcUhAT58+LB4eHjInDlz7F1VRHgSi4XF1qJ4gI8ePfpE4WWAWVhsL7YE2K7rwL6+viguLka3bt2s6seMGYMxY8bY1AavAxPZRvF5obdu3YqSkpIy9c8995x9IyMiRdi1BVYCt8BEtrFlC8wvMxCpGANMpGIMMJGKMcBEKsYAE6kYA0ykYgwwkYoxwEQqxgATqRgDTKRiDDCRijHARCrGABOpmN2T2l27dg0//vgjsrOz0aJFC7z66qtwc3NzxNiIqAp2bYGXLFmCqKgopKSkQESwYMEChIaG4ubNm44aHxFVxp7pcP7973+L2Wy2PM7PzxdfX19ZsGABp9RhYVG42DKljl270KGhoVaPTSYTSkpK4OPjY08zRKQQu4+Bb926hS+//BJ5eXn46aef8MYbb2DChAkVLs+J3Ykcx+6z0K6urggICICfnx+cnZ3xyy+/4P79+xUuP3fuXNSpU8dSAgMDn2jARPRfTzQnVlFREdq3b4/IyEisXLmywmUe3wIzxERVU3xWysdptVp07NgRly5dqnQZrVb7JN0QUQXs2oU+duyY1eN79+7h+PHj6NChg6KDIiLb2LULPXLkSKSmpiI8PBwFBQWIiYlBWFgYvv/+e/j6+trUBqeVJbKNLbvQdh8DX7hwAXFxcXB1dUVYWBg6depk16AYYCLbOCTAT4oBJrINJ3Yn+p1jgIlUjAEmUjEGmEjFGGAiFWOAiVSMASZSMQaYSMUYYCIVY4CJVIwBJlIxBphIxRhgIhWrdoCzs7OxevVqHDp0SMnxEJEdqjWljohg3LhxOHLkCAYMGIA+ffooPS4iskG1tsCLFi2CiKB3795Kj4eI7GB3gM+cOYPFixfjm2++gZOTkyPGREQ2smsXOjc3F6NGjcLKlSsREBBg0zqc2J3IcezaAk+aNAl9+/bF4MGDbV6HE7sTOY7Nc2IdPXoU/fr1w8KFC+Hh4QEAWL58OVxdXfHOO+9g1KhRqF27dpn1OLE7UfUoOrG7j48Pxo4di/Pnz1vq7t+/DxcXF5w6dQpRUVHlBpgTuxM5zhPNSjlw4EC4ublhx44dNq/DWSmJbMNZKYl+557ot5H69esHV1dXpcZCRHbixO5EzyjuQhP9zjHARCrGABOpGANMpGIMMJGKMcBEKsYAE6kYA0ykYgwwkYoxwEQqxgATqRgDTKRiDDCRitn1dcL09HSkpqZaN6DRoEuXLooOiohsY1eA169fjwULFqBdu3aWOi8vL8TExCg+MCKqmt1f6A8NDcWJEyccMRYispPdx8Amkwnnzp3DlStXYDKZHDEmIrKR3QFOTEzE6NGj0adPHwQEBGDjxo2VLl9UVIScnByrQkQKETvs27dPfvvtN8vjxYsXi7Ozs8TGxla4zsyZMwUACwuLncVgMFSZySeeE6tVq1YYMGAAvvrqq3Kf58TuRNWj6MTuFWnQoAEyMjIqfJ4TuxM5jl3HwHl5eVaP79y5g7Nnz6JNmzaKDoqIbGPXFjgyMhJDhw5Fhw4dkJ2djQULFsDf3x/vvvuuo8ZHRJWwawu8f/9+lJSUYMWKFdizZw/Gjx+Pc+fOwdfX11HjI6JKcGJ3omcUJ3Yn+p1jgIlUjAEmUjEGmEjFGGAiFWOAiVSMASZSMQaYSMUYYCIVY4CJVIwBJlIxBphIxRhgIhWrdoBNJhOMRqOSYyEiO9kd4OTkZPTv3x+enp4ICAjA66+/jnv37jlibERUBbsCfOvWLXTv3h3+/v7IysrC3bt3MXToUMTFxTlqfERUGXumlZ08ebIEBQWJ0Wi0ZzUrBoPhqU/XycKihmLLtLJ2bYH37t2LoUOHwtXVtcwEd0RU8+wKcFpaGjQaDTp37oz69eujdu3aGD9+PO7fv1/hOvxlBiIHsmf318XFRbRarRw6dEhERFJSUqRFixYyatSoCtfhLzOwsFSv2LILbVeAg4KCJCoqyqpu6dKl4uXlJWazudx1CgsLxWAwWEp6evpTf2NYWNRQbAmwXfNC9+rVC1lZWVZ1+fn50Gq1cHJyKncd/jIDkQPZswW+cOGCeHh4yIoVKyQ1NVV+/PFH8ff3l08++cTmNngWmoXFtqL4LrSIyIkTJ6R3797SuHFj6dSpk3z11VdSXFzMALOwKFxq5NcJ7cWJ3Ylsw4ndiX7nGGAiFWOAiVSMASZSMQaYSMUYYCIVY4CJVIwBJlIxBphIxRhgIhVjgIlUjAEmUjEGmEjFGGAiFbMrwF27dkW9evXKlLFjxzpqfERUCbum1ImJiYHJZLI8/u233xAREYHIyEjFB0ZEVbMrwD4+PlaPly9fDi8vL4wcOVLJMRGRjap9DGw2m7Fu3TqMHj0anp6eSo6JiGxk1xb4UQcOHEBaWhomTpxY6XJFRUUoKiqyPObE7kQKsndSu1IjRoyQDh06VLkcJ3ZnYaleccislCIi2dnZUqtWLVm5cmWVy3JidxaW6hXFJ3YvtWHDBmg0GowePbrKZTmxO5HjVOsk1tq1azFy5EhOD0v0lNkd4Li4OFy8eLHKk1dE5Hh270KHh4fjzp07qFevniPGQ0R2sDvAbm5ucHNzc8RYiMhO/DIDkYoxwEQqxgATqRgDTKRiNR5gqdlfMyVSLVuyUuMBzs3NrekuiVTJlqzU+A98m81mZGZmwsvLC05OTlUun5OTg8DAQKSnp1f5Y8fVVRN91FQ/7EP9fYgIcnNz0ahRIzg7V76NrfbXCavL2dkZTZo0sXs9b29vh4arpvqoqX7Yh7r7sPU2ZZ7EIlIxBphIxZ75AGu1WsycOdOhX0msiT5qqh/28b/VR42fxCIi5TzzW2AiqhgDTKRiDDCRitX4dWB7mEwm/Pzzz7h//z46duyIwMBAh/Rz6dIlnD9/Ht27d6/WNeqqFBQUIDExEQaDAWFhYQgKClK8DwD49ddfkZycjICAAHTq1KnKmwCexOHDh5GVlYUhQ4bA3d1dsXZ3796N/Px8q7rQ0FC0a9dOsT5K3b59G6dPn0b9+vURERGh6Pu1c+dOq+mUSwUHB6NLly6K9fPMBjgzMxMvvfQSCgoKEBISgtjYWMyaNQtTp05VrI/Y2FhMnz4dN2/eRFJSErZv347hw4cr1j4ArFy5EvPnz0eTJk3g7e2N48eP480338Ty5csV6yM5ORlvvPEGcnJyEBISgvPnz0Oj0WDv3r1o2bKlYv2UOnLkCAYOHIiCggKkp6cr+qH37rvvonHjxggJCbHUiYiiARYRzJgxA0uWLEG3bt1QUlICEcHu3btRu3ZtRfqIiYmxmgO9sLAQu3btwvTp0xUNcLXnhXa0YcOGSefOnaWwsFBERLZv3y5OTk5y/vx5xfrYt2+fHD58WAoKCgSAbN++XbG2S23YsEGys7Mtj3/55RfRaDSydetWxfpITEyUc+fOWR6bTCZ54YUXZPDgwYr1UerOnTvStGlTmTNnjgCQ9PR0Rdtv3LixfPPNN4q2+bilS5eKp6enJCQkWOpOnDhh9XdS2saNG8XJyUmSk5MVbfeZDHBOTo5oNBpZv369VX3Tpk1l2rRpivfnyACXp3nz5jJ9+nSH9jF27FiJjIxUtE2z2SwDBgyQ2bNny4EDBxwW4OnTp8vOnTvlzJkzYjQaFW3fbDZLw4YNZcqUKYq2W5UXX3xR8b+HSDXnhXa0y5cvw2QyoU2bNlb1bdq0wYULF57SqJRx7do1XL9+vcxrU8Lu3bthMBiQmJiIY8eOYdu2bYq2v3jxYhgMBkybNg1HjhxRtO1H7dixA+fPn0diYiK8vb2xZcsWhIeHK9J2cnIybt68iZdffhmXLl1CcnIygoODFWu/PCkpKTh27Bg2b96seNvPZIANBgMAwNfX16rez88PKSkpT2NIiigsLMTrr7+O8PBwjBgxQvH29+/fj4yMDPzyyy9o164dGjZsqFjbCQkJmDdvHuLj4+Hi4qJYu49btWoVXnnlFQAP36/hw4dj5MiRuHjxIjSaJ//vmpWVBeDh3OaJiYlo2bIl4uLi0KpVK+zdu9chX2pYu3Yt6tati1dffVXxtp/Jy0ilt509ePDAqv7BgweqnRHTaDRi+PDhyM7Oxu7du+Hq6qp4HytWrEB0dLTlQ07JD4lJkyahd+/eOHXqFLZs2YKjR48CeLjVT0xMVKyf0vACD2dA/eyzz3DlyhUkJSUp0n7p/5/79+/j0qVL2LNnDy5fvoxr165h1qxZivTxqJKSEqxfvx7jxo1zyO2Uz2SAdTodACAtLc2qPjU1Fc2aNXsaQ3oiRqMRI0aMwOXLl3HkyBE0atTIof25urpi5MiRiI+PL/dSRnV07doVABAdHY3o6GicPHkSwMOtviMPa0q3iHfu3FGkvdL/P0OGDLFs0X19fREZGYmEhARF+njUvn37kJGRgbffflvxtgE8u2ehw8PDZcKECZbHycnJ4uTkJNHR0Yr35ciTWEajUYYMGSI6nU7xEz6lMjMzy9R98MEH4u/v75D+RMQhJ7Hu3LkjxcXFVnVffPGFaLVauXfvnmL9dOvWTd5//32ruk6dOsm4ceMU66NUVFSUvPDCC4q3W+qZPAYGgEWLFuEPf/gD3N3d0bp1a/z9739Hnz59MHjwYMX6yMzMxPHjx1FcXAwAOHnyJEwmE1q0aIGOHTsq0sdbb72FH374AXPnzsWJEycs9SEhIdDr9Yr0sWjRIly9ehUvvvgiPDw8cOLECWzZsgVr165VpP2acunSJXzwwQeIiopC48aN8fPPP2PTpk1YuHAh6tatq1g/ixYtwksvvQSNRoO2bdsiJiYGSUlJ2LRpk2J9AA+Pt/fu3Yt//vOfirb7qGf620hnz57F+vXrcf/+fURERGDixImoVauWYu0nJiZi/vz5Zer79euHN954Q5E+pk6dioyMjDL1vXr1wh//+EdF+gCAgwcPWm4eCAkJweuvv46mTZsq1v7jLl68iNmzZ2PlypVlTjY+iWvXrmHTpk1ITU1FYGAgRowYgbCwMMXaL5WUlIQ1a9YgKysLOp0Ob7/9NgICAhTt48iRI1i1ahVWr14NDw8PRdsu9UwHmIgq90yexCIi2zDARCrGABOpGANMpGIMMJGKMcBEKsYAE6kYA0ykYgwwkYoxwEQqxgATqRgDTKRi/wcVBpIr4Jai3AAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 250x270 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "Q_orth = orthogonal_iteration(A, Q_orth)\n",
        "B, V = qr_iteration(B)\n",
        "W = W @ V\n",
        "k += 1\n",
        "signs = np.where(np.sum(Q_orth * W, axis=0) < 0, -1., 1.)\n",
        "B_orth = Q_orth.T @ A @ Q_orth\n",
        "B_aligned = signs[:, None] * B_orth * signs[None, :]\n",
        "print(f\"Step {k}: raw difference {la.norm(B-B_orth)/la.norm(A):.2e}, \"\n",
        "      f\"after sign alignment {la.norm(B-B_aligned)/la.norm(A):.2e}\")\n",
        "show_fill([(B, f\"Unshifted QR: step {k}\")])"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "qr-demo-09",
      "metadata": {},
      "source": [
        "## 2. Reduce to upper Hessenberg form\n",
        "\n",
        "For column $j$, use a Householder reflector on rows $j+1,\\ldots,n-1$.\n",
        "Apply it on both sides to preserve $H=Z^TAZ$. The supplied `House(x)` returns a unit\n",
        "vector $v$ defining $I-2vv^T$ (or $v=0$ for a zero input).\n",
        "Why do we preserve the first subdiagonal rather than eliminate it too?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "qr-demo-10",
      "metadata": {},
      "outputs": [],
      "source": [
        "def House(x):\n",
        "    \"\"\"Unit vector v for the Householder reflector I - 2vv^T; zero if x is zero.\"\"\"\n",
        "    v = np.array(x, dtype=float, copy=True)\n",
        "    normx = la.norm(v)\n",
        "    if normx != 0:\n",
        "        v[0] += np.copysign(normx, v[0])   # Choose sign to avoid cancellation.\n",
        "        v /= la.norm(v)\n",
        "    return v\n",
        "\n",
        "H, Z, j = A.copy(), np.eye(n), 0\n",
        "hess_frames = [(H.copy(), \"Dense A\")]\n",
        "# Use House(H[j+1:, j]) and apply I - 2vv.T from both sides; also update Z."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "qr-demo-11",
      "metadata": {},
      "outputs": [],
      "source": [
        "def hessenberg_step(H, Z, j):\n",
        "    v = House(H[j+1:, j])\n",
        "    H[j+1:, :] -= 2 * np.outer(v, v @ H[j+1:, :])\n",
        "    H[:, j+1:] -= 2 * np.outer(H[:, j+1:] @ v, v)\n",
        "    Z[:, j+1:] -= 2 * np.outer(Z[:, j+1:] @ v, v)\n",
        "    H[j+2:, j] = 0.                     # Clear computed roundoff below the subdiagonal.\n",
        "    return H, Z"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "qr-demo-12",
      "metadata": {},
      "source": [
        "Run the next cell in place to eliminate one column at a time (six reflectors).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "qr-demo-13",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPAAAAEKCAYAAAAsK9hZAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHYpJREFUeJzt3XtU1GX+B/D3yEyAXEUh1FQIEhUUb4RJpugSbpaXyrXVTN3qtNtm1z272c3z09Zra6RJFyw1daP0JBGFrhnelbyMFzQExZRAA4wdQLnP5/eHP+bnyG2+8J2Br75f5zznOM98v8/zMPLme39GJyICItKkDm09ACJqOQaYSMMYYCINY4CJNIwBJtIwBphIwxhgIg1jgIk0TN/WA9AiEUFeXh7y8/Ph5eWFkJCQth5Sq5hMJuTk5KCystJS5+rqivDw8DYcFdmCAb7OlStXcOLECbi7uyMsLKzBZYqKivDQQw8hOzsbQUFBGDFiBF577TVkZWWhf//+cHNzc/CoW+fNN9/EO++8g759++K2224DAGRlZcHPzw+ZmZl26TMrKwvl5eVt9geisrISZ8+ehaurKwIDA9tkDKoRsli2bJkAEBcXFykuLm5wmTlz5oinp6f897//tdRt3LhRAMjBgwcdNFJ1ZGRkCAD5/PPPreonTJggISEhduvX3u03pqKiQv72t7+Jp6enhISESGhoqPTp00f+85//OHwsauEx8HU++eQT3H333aiursb69esbXObYsWMIDg6Gl5eXg0envqNHjwIAhgwZ0rYDcQARwSOPPIKUlBSkp6cjMzMTGRkZSElJwc8//9zWw2sx7kL/n3379uHkyZP47rvv8MEHHyAhIQHPPfec5f2ysjJkZGTgl19+QW1tLQ4cOAAA8Pb2RnZ2NgAgIyMDNTU1AIDg4GB06dLFqo/8/HxcunQJfn5+uOOOO6zeKyoqwpkzZxAeHg5XV1fk5eUhLy8PgwcPhl5f/7/JluVra2tx5swZlJeXIygoCB4eHpb1jUYjjEajZdyXL1+Gr68vgoKCmvycmmrzemazGWfPnkVlZSV69+5t2T0/ceIEiouLUVFRYfkMAWDYsGFW65eXlyM7OxtOTk7o3bs3DAaD1fuZmZmorq5G//79UVNTg+zsbOh0OvTp06fB8axbtw7ffvstjEaj1TJBQUHN/sztWlvvArQXM2fOlDvvvFPMZrOkpqYKAElPT7e8n5GRIZGRkeLh4SHu7u4SGRkpkZGR8tZbb8ldd90lACQsLMxSv2XLFsu6Bw8elKFDh4qrq6uEhoZa1s/KyrIss27dOgEgu3fvlgceeECCg4PFxcVFCgsLGxxvc8vHxcVJ586dxc/PT0JCQsTFxUVmz54tVVVVIiIyceJECQgIEAAycOBAiYyMlLlz54pI47u4zbV5/XK+vr7i5eUloaGh4u3tLa+99pqYzWaZPn26dOrUSVxcXCyfVWRkpGXdqqoqefnll8XFxUV69eol/v7+4unpKUuXLrXqIzY2VsLDw2Xnzp0SHBwswcHBMnbs2Eb/f0eMGCEDBgwQEZGLFy/K4cOH5eLFi40urxUMsIiYTCZxc3OTJUuWiIiI2WyWoKAgefrpp+stGxUVZfULJ9L0MfBPP/0k7u7uMnHiRMtxc2lpqYwdO1YCAgKkrKxMRP4/kNHR0bJ//34RETl9+rSYTKYGx9zU8gsXLhSdTierVq2yLH/48GHx8fGRF1980VK3YsUKASDnzp2zaruhANva5vz580Wn08ny5cvFbDaLiEh5ebn885//lMrKykbbrzN79mzR6/WSlJRkqVu+fLkAkLi4OEtdbGys9OzZUyZPnmz5jBo7B2E2m8XZ2VnGjBkjMTEx4unpKX379hW9Xi8jR46Un3/+ucH1tIABFpEPPvhAXFxcpKioyFK3dOlScXd3l9LSUqtllQZ42rRp4uHhYXXSS0Tk3LlzAsASiLpAvvnmmzaNubHlS0tLxc3NTaZOnVpvnfnz58ttt90mJSUlImJ7gG1t02Qyiaurq/zhD39ocuyNBfjy5ctiMBhkxowZ9d677777xNfX1/JHITY2VgBIZmZmk32JXPsDDUAAyPjx4+Xq1asiIpKZmSldu3aVPn361NuL0AoeAwNISEjAsGHDkJ2dbTme7devH65cuYLExEQ89dRTLW57x44dCA4OtrQr1/5oQkTg4eGBw4cP48knn7QsP2LECEXt37j8oUOHcOXKFQQGBuLQoUOQ/5uvQUTQsWNHVFVV4cSJExg+fLjNfdjaZkVFBcrLy3H//fcr+hnqHDlyBNXV1Rg1alS990aNGoVdu3YhJyfHcszq5+dn0zV4Jycny7/ff/99uLq6AgBCQkLwj3/8Ay+++CK2bNmChx56qEXjbku3fICNRiOOHDmCyMhIvPjii1bvBQYGYtWqVa0KcElJCSoqKqxOiNXp168fbr/9dqu6G18358blS0pKAACbNm3C999/X2/5yMhIdOig7OKDrW2WlZUBADp16qSo/TpXr14FgAbP8Ht7ewO4dq2+jq2flZubGzw9PeHk5IQePXpYvTdw4EAAwOnTpxlgLUpISEB4eLjVGdE6hw4dQkREBE6cOIH+/fs32oZOp2v0vZ49e8LDwwP79++3aTxKw3Xj8j179gQAzJ49G3/9618VtdUYW9s8duwYAODMmTNNttfY59W9e3cAwPnz5+u9d+7cOQCwOnuv5LOKiIjA7t27YTabrdar++OktRtw6tzS14HLy8vx73//G2PHjm3w/SFDhsDPzw8JCQlNttPQ1qHO448/joMHDyI9Pb3eeyLS4DqtER4ejrCwMHz88ceoqqqq935paand2hwwYABCQ0Px4YcfoqKiwmqZiooKy663t7d3gz/3wIED0atXL6xevRq1tbWW+rKyMiQmJmLEiBHw8fFRPH4AmDVrFqqqqvDdd99Z1W/evBkdOnRATExMi9pta7d0gL/88kuYTKZGA6zT6RAbG4v169fX+4W83sCBA+Hu7o4VK1Zg586dOHDgAIqKigAAr7zyCsaNG4fY2FjMmzcPW7ZsQUpKCpYsWYLQ0FAcPHhQ1Z9Jp9Ph888/R2FhIYYNG4bVq1cjLS0Nn3/+OWbPnl3vequabep0OmzYsAGlpaWIiIjAmjVrsH37drz33nsYNGgQqqurAQD33nsvfvnlF8THx2Pfvn2WvR8nJyckJCQgKysLDzzwAL7++mt8+eWXGDlyJMxmMz744IMWfy5Tp07FpEmTMGPGDMTHx2Pbtm144YUXsGbNGsybNw/BwcEtbrst3dK70EajESNGjEBUVFSjy0yZMgVZWVn48ccfcd999yE0NNRq6wAAnTt3RnJyMj766CO8/vrrqKmpwf/8z/8gNjYWBoMBSUlJ2Lx5M5KSkrB9+3Z4e3sjNDQUycnJll8cX19fREZGomPHjjaNvanlw8LCcPLkSaxatQqbN2+GyWRCjx49EBUVhUWLFlmW8/f3R2RkJJydna3WDwkJgYuLS4vaDA8PR0ZGBj788EMkJiaitrYW4eHh2Lp1q+VmjlmzZqG4uBjJyclYt24dRMQS4piYGBiNRqxcuRJxcXFwcnJCTEwMZs+ebdnFBoC+ffvCz8/Pps8KuPbHZdOmTVizZg2Sk5ORmJiInj17YuvWrZrd+gKAToTTyhJp1S29C02kdQwwkYYxwEQaxgATaRgDTKRhDDCRhjn8OrDZbEZ+fj48PDyavAWR6FYlIigtLUW3bt2avV3U4QHOz8+vd0M5EdWXm5tbb+aWGzl8F7qxKViIyJotWXF4gLnbTGQbW7LCk1hEGsYAE2kYA0ykYQwwkYa1OMBms1nNcRBRCygOcEJCAnr27AmDwYDg4GBs2rTJHuMiIlsomYM2OTlZDAaDfPHFF3LlyhVZuXKlODk5yYEDB2xu4/o5ellYWBovjU3qfz1FM3KMHj0aPj4+VlvdYcOGISgoCBs2bLCpjZKSkpvii8GI7M1kMsHT07PJZWzehRYRpKenY+TIkVb10dHRNk+ZSkTqsvle6NLSUly9ehW+vr5W9X5+fvj1118bXa+ystLqm9/r5uElotZTfBLrxrPPZrO5yVu+Fi5cCC8vL0vhgwxE6rE5wB4eHnBzc0NBQYFVfUFBAfz9/Rtdb86cOTCZTJaSm5vb8tESkRWbA6zT6RAVFYW0tDSr+u3btzc5r7KzszM8PT2tChGpRMllpG3btomTk5N89NFHkp+fLwsWLBCDwSBHjx7lZSQWFpWLLZeRFH8/cGJiooSFhYmnp6cMGTLE6pvoGWAWFvWK6teB1cDrwES2UfU6MBG1PwwwkYYxwEQaxgATaRgDTKRhDDCRhjHARBrGABNpGANMpGEMMJGGMcBEGsYAE2kYA0ykYQwwkYYpCrCIYNu2bXj44YcREBCApKQkOw2LiGyhKMBxcXFYvHgxpk6divPnz6OsrMxe4yIiG9g8rSwAzJ49Gy+99JK9xkJECinaAuv1ivJORHZm90RyYnci+7H7WWhO7E5kP3YPMCd2J7Ifu+9COzs7w9nZ2d7dEN2SeCMHkYYpCvDevXsREBCAgIAAAMBLL72EgIAAvPrqq/YYGxE1Q9HE7hUVFbh06VK9eg8PD3Tu3NmmNjixO5FtbJnYXdExsIuLi2XrS0Rtj8fARBrGABNpGANMpGEMMJGGMcBEGtZmjxfZcoq8NXQ6nd3aJmovuAUm0jAGmEjDGGAiDWOAiTSMASbSMAaYSMMYYCINU3QduKqqChs2bMDOnTtRXV2NiIgIPPPMM3B1dbXX+IioCYq2wMOHD8eePXswatQojB07FqtXr0ZUVBTKy8vtNT4iaoooUFhYaPU6Pz9fdDqdbNq0yeY2TCaTABCTyaSka8UAsLBoutiSEUVb4C5duli97tSpE/R6Pb9ihaiNtOpe6OXLl6NDhw4YM2ZMo8twYnciO2rpLmpqaqro9Xr56KOPmlxu7ty5Ld49aI2G+mRh0VKxJSMtCvD27dvF1dVVFi1a1OyyFRUVYjKZLCU3N9fmwbVGW3/4LCytLXYJ8A8//CAdO3aUBQsWtChYPInFwmJbUT3AO3bsaFV4RRhgFhZbiy0ZUTQvtI+PD6qrqxEVFWVV//jjj+Pxxx+3qY26eaH5QD9R01SfF/qLL75AbW1tvfrg4GBlIyMiVSgKcExMjL3GQUQtwIcZiDSMASbSMAaYSMMYYCINY4CJNIwBJtKwNvtmBntTcH9Ki/FmEWpr3AITaRgDTKRhDDCRhjHARBrGABNpGANMpGGKLyPl5OTgu+++Q1FREXr37o2HH34YLi4u9hgbETVD0RY4Li4OkyZNwtmzZyEiWLJkCfr164eLFy/aa3xE1BQl09ScPHlSzGaz5fXVq1fFx8dHlixZYnMbjppSxxHQDqZdYbl5iy0ZUbQL3a9fP6vXNTU1qK2thbe3t5JmiEglio+BL126hHfeeQdXrlzB7t27MWvWLMycObPR5TmxO5H9KD4LbTAY4O/vj86dO6NDhw44cuQIiouLG11+4cKF8PLyspQePXq0asBE9P8UzUp5o8rKSgwcOBDR0dGIj49vdJkbt8A9evSw+6yUjsCHGcieVJ+V8kbOzs4YPHgwTp061eQyzs7OremGiBqhaBd6586dVq9/++037Nq1C4MGDVJ1UERkG0W70FOmTMH58+cRHh6O8vJypKamIjQ0FF999RV8fHxsasNRE7s7AnehyZ5syYjiY+ATJ04gPT0dBoMBoaGhGDp0qKJBMcBEtrFLgFuLASayjS0Z4cMMRBrGABNpGANMpGEMMJGGMcBEGsYAE2nYTTuxuyM44gocL1VRU7gFJtIwBphIwxhgIg1jgIk0jAEm0jAGmEjDWhzgoqIirFq1Ctu3b1dzPESkQIuuA4sInnjiCaSlpWHcuHEYM2aM2uMiIhu0aAu8bNkyiAhGjx6t9niISAHFAT506BDeffddrF69mncJEbUxRbvQpaWleOyxxxAfHw9/f3+b1uHE7kT2o2gL/Mwzz+B3v/sdxo8fb/M6nNidyH5snhNrx44diI2NxdKlS9GxY0cAwIoVK2AwGPDnP/8Zjz32GNzd3eutdzNP7O4IPEy5dak6sbu3tzemT5+O48ePW+qKi4vh5OSEAwcOYNKkSQ0GmBO7E9lPq2alfPDBB+Hi4oJNmzbZvM7NNCulI3ALfOvirJREN7lWPdAfGxsLg8Gg1liISCFO7N7OcRf61sVdaKKbHANMpGEMMJGGMcBEGsYAE2kYA0ykYZzYvZ3j5PHUFG6BiTSMASbSMAaYSMMYYCINY4CJNIwBJtIwRZeRcnNzcf78eesG9HoMGzZM1UERkW0UBXjt2rVYsmQJBgwYYKnz8PBAamqq6gMjouYpvpGjX79+2LNnjz3GQkQKKT4GrqmpwbFjx3D69GnU1NTYY0xEZCPFATYajZg6dSrGjBkDf39/rFu3rsnlKysrUVJSYlWISB2KAhwREYGcnBycPHkSv/zyC9544w3MnDkTBw4caHQdTuxOZD+tnhOrT58+GDduHP71r381+D4ndm//+DBD+6TqxO6Nuf3225GXl9fo+5zYnch+FO1CX7lyxep1YWEhjh49irCwMFUHRUS2UbQFjo6OxsSJEzFo0CAUFRVhyZIl8PPzw7PPPmuv8RFRExRtgbdu3Yra2lqsXLkS33zzDWbMmIFjx47Bx8fHXuMjoiZwYnfiSax2ihO7E93kGGAiDWOAiTSMASbSMAaYSMMYYCIN48Tu5JDJ4wFerrIHboGJNIwBJtIwBphIwxhgIg1jgIk0jAEm0rAWB7impgZVVVVqjoWIFFIc4OzsbDzwwANwc3ODv78/pk2bht9++80eYyOiZigK8KVLl3DvvffCz88PBQUFuHz5MiZOnIj09HR7jY+ImqDogf7nn38eX3/9Nc6cOQODwdCiDvlA/62Ld2Ipo/oD/SkpKZg4cSIMBkO9Ce6IyPEUBfjChQvQ6/W4++674evrC3d3d8yYMQPFxcWNrsNvZiCyI1HAyclJnJ2dZfv27SIicvbsWendu7c89thjja4zd+5cAVCvmEwmJV3TTaCh3wOWxostGVF0DNyrVy8MGTIEX331laVu+fLleOONN2AymRo8xuE3M1AdHgMro/o3M4wcORIFBQVWdVevXoWzs3Oj/zn8ZgYi+1F0DPz3v/8du3fvRnx8PC5cuIDU1FS8++67ePLJJ+01PiJqguJ5offu3Yu33noLp0+fRteuXfHHP/4Rzz//PPR62zbmvIx06+IutDK2ZIQTu5PDMMDKcGJ3opscA0ykYQwwkYYxwEQaxgATaRgDTKRhnNidHMYRVyxvtUtV3AITaRgDTKRhDDCRhjHARBrGABNpGANMpGEMMJGGKQrwPffcgy5dutQr06dPt9f4iKgJim7kSE1NRU1NjeX1zz//jIiICERHR6s+MCJqnqIAe3t7W71esWIFPDw8MGXKFDXHREQ2avExsNlsxpo1azB16lS4ubmpOSYislGL74Xetm0bLly4gKeffrrJ5RqaVpaI1NHiLfAnn3yCQYMGYciQIU0ut3DhQnh5eVlKjx49WtolEd2gRZPaXb58Gd26dUNcXBz+8pe/NLksJ3YnR7qZnkZSfWL3Op999hn0ej2mTp3a7LKc2J3Iflq0C/3pp59iypQp8PLyUns8RKSA4gCnp6cjIyOj2ZNXRGR/inehw8PDUVhYiC5duthjPESkgOIAu7i4wMXFxR5jISKF+DADkYYxwEQaxgATaRgDTKRhDp8Xuu7GL94TTdQ0W26SdHiAS0tLAYD3RBM1o7S0tNmbpRz+Bd9msxn5+fnw8PCw6b7Vununc3Nz7XbvtCP6cFQ/7EP7fYgISktL0a1bN3To0PRRrsO3wB06dMAdd9yheD1PT0+7P/zgiD4c1Q/70HYftt6mzJNYRBrGABNpWLsPsLOzM+bOnWvXRxId0Yej+mEft1YfDj+JRUTqafdbYCJqHANMpGEMMJGGOfw6sBI1NTXYt28fiouLMXjwYLvdvXXq1CkcP34c9957b4uuUTenvLwcRqMRJpMJoaGh6Nmzp+p9AMBPP/2E7Oxs+Pv7Y+jQoc3eBNAaP/zwAwoKCjBhwgS4urqq1m5ycjKuXr1qVdevXz8MGDBAtT7q/Prrrzh48CB8fX0RERGh6ue1efNmq8kc6wQEBGDYsGGq9dNuA5yfn4+YmBiUl5cjMDAQ+/fvx/z58/HKK6+o1sf+/fsxZ84cXLx4EVlZWdi4cSMeffRR1doHgPj4eCxevBh33HEHPD09sWvXLvzpT3/CihUrVOsjOzsbs2bNQklJCQIDA3H8+HHo9XqkpKQgJCREtX7qpKWl4cEHH0R5eTlyc3NV/aP37LPPonv37ggMDLTUiYiqARYRvPHGG4iLi0NUVBRqa2shIkhOToa7u7sqfaSmplrd719RUYGvv/4ac+bMUTXAkHbqkUcekbvvvlsqKipERGTjxo2i0+nk+PHjqvWxZcsW+eGHH6S8vFwAyMaNG1Vru85nn30mRUVFltdHjhwRvV4vX3zxhWp9GI1GOXbsmOV1TU2NDB8+XMaPH69aH3UKCwulV69esmDBAgEgubm5qrbfvXt3Wb16tapt3ui9994TNzc3OXz4sKVuz549Vv9Palu3bp3odDrJzs5Wtd12GeCSkhLR6/Wydu1aq/pevXrJq6++qnp/9gxwQ+666y6ZM2eOXfuYPn26REdHq9qm2WyWcePGydtvvy3btm2zW4DnzJkjmzdvlkOHDklVVZWq7ZvNZunatau8/PLLqrbbnFGjRqn+/yEi0i53oTMzM1FTU4OwsDCr+rCwMJw4caKNRqWOnJwcnDt3rt7Ppobk5GSYTCYYjUbs3LkTX375partv/vuuzCZTHj11VeRlpamatvX27RpE44fPw6j0QhPT08kJiYiPDxclbazs7Nx8eJF3H///Th16hSys7MREBCgWvsNOXv2LHbu3IkNGzao3na7DLDJZAIA+Pj4WNV37twZZ8+ebYshqaKiogLTpk1DeHg4Jk+erHr7W7duRV5eHo4cOYIBAwaga9euqrV9+PBhLFq0CD/++COcnJxUa/dGCQkJ+P3vfw/g2uf16KOPYsqUKcjIyIBe3/pf14KCAgDX5jY3Go0ICQlBeno6+vTpg5SUFLs81PDpp5+iU6dOePjhh1Vvu11eRqq77aysrMyqvqysTLMzYlZVVeHRRx9FUVERkpOTYTAYVO9j5cqVSEpKsvyRU/OPxDPPPIPRo0fjwIEDSExMxI4dOwBc2+objUbV+qkLL3BtBtTXX38dp0+fRlZWlirt1/3+FBcX49SpU/jmm2+QmZmJnJwczJ8/X5U+rldbW4u1a9fiiSeesMvtlO0ywEFBQQCACxcuWNWfP38ed955Z1sMqVWqqqowefJkZGZmIi0tDd26dbNrfwaDAVOmTMGPP/7Y4KWMlrjnnnsAAElJSUhKSsLevXsBXNvq2/Owpm6LWFhYqEp7db8/EyZMsGzRfXx8EB0djcOHD6vSx/W2bNmCvLw8PPXUU6q3DaD9noUODw+XmTNnWl5nZ2eLTqeTpKQk1fuy50msqqoqmTBhggQFBal+wqdOfn5+vboXXnhB/Pz87NKfiNjlJFZhYaFUV1db1c2bN0+cnZ3lt99+U62fqKgoee6556zqhg4dKk888YRqfdSZNGmSDB8+XPV267TLY2AAWLZsGcaOHQtXV1f07dsX77//PsaMGYPx48er1kd+fj527dqF6upqAMDevXtRU1OD3r17Y/Dgwar08eSTT+Lbb7/FwoULsWfPHkt9YGAgIiMjVelj2bJlOHPmDEaNGoWOHTtiz549SExMxKeffqpK+45y6tQpvPDCC5g0aRK6d++Offv2Yf369Vi6dCk6deqkWj/Lli1DTEwM9Ho9+vfvj9TUVGRlZWH9+vWq9QFcO95OSUnBxx9/rGq712vXTyMdPXoUa9euRXFxMSIiIvD000/jtttuU619o9GIxYsX16uPjY3FrFmzVOnjlVdeQV5eXr36kSNHNvvVrEp8//33lpsHAgMDMW3aNPTq1Uu19m+UkZGBt99+G/Hx8fVONrZGTk4O1q9fj/Pnz6NHjx6YPHkyQkNDVWu/TlZWFj755BMUFBQgKCgITz31FPz9/VXtIy0tDQkJCVi1ahU6duyoatt12nWAiahp7fIkFhHZhgEm0jAGmEjDGGAiDWOAiTSMASbSMAaYSMMYYCINY4CJNIwBJtIwBphIwxhgIg37X3o5Ulg/eq88AAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 250x270 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "if j < n - 2:\n",
        "    H, Z = hessenberg_step(H, Z, j)\n",
        "    j += 1\n",
        "    hess_frames.append((H.copy(), f\"After reflector {j}\"))\n",
        "show_fill([hess_frames[-1]])\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "qr-demo-14",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAEVCAYAAAAW3T11AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYltJREFUeJzt3XdYU3f7P/A3S0RkKOJGBcG96satKIiz7tE6qtU6Wmetta3ap+0j2ro6nTjrqgvcG7UqpSioKKiAgIgDFwkzQPj8/vBHvkaGCZxDgOf9uq7PdZmTcz73ncid5E7OMBJCCBARERERERGRLIwNnQARERERERFRacbGm4iIiIiIiEhGbLyJiIiIiIiIZMTGm4iIiIiIiEhGbLyJiIiIiIiIZMTGm4iIiIiIiEhGbLyJiIiIiIiIZMTGm4iIiIiIiEhGpoZOgIioOIuKisL58+cRHx+Pli1bolevXoZOqcDS0tJw9OhRREdHIyMjQ7PcyckJw4cPN2BmREVLqVTi2LFjiI2NRcWKFTFx4kRDp1QoZ8+exe3bt5GSkqJZZmpqis8//9yAWREVraysLJw4cQL37t1DWloaZs6cCQsLC0OnVWAhISHw9/fHy5cvtZaPHTsW1atXN1BWVBhsvIuhsLAw+Pr6am6XKVMGNjY2qFevHlq1aoVy5coZMDt5PHjwADt37oSdnR0mTZpk6HTof0BISAiOHj2ab9N55swZ9O/fH3379oWzszNSU1Px77//4ty5cyXuDT01NRUtW7ZEmTJl0KtXL5iavn7537x5M9577z3ZGu/NmzfD0tLSII19eHg4jh8/jpSUFEydOhU2NjZFngMVLYVCgTVr1qBMmTKYPXs2jIyMcqwTHx+PFi1aoHbt2ujcuTPMzMzw5MkTbNmyBYMGDUL9+vUNkHnBjR8/HsePH8fIkSM1r0lnz55FSEiIbI332bNnERwcXOSNfVZWFvz9/XHr1i0YGxvD2dkZXbt2hbExd+As7dasWQOFQpFv0+nh4YHw8HAMGjQI5ubmEEJg5cqVaNq0aYn70nzTpk2YNm0aRo8ejcqVKwN4/Z524MABdOvWTZbGOzw8HPv37zdYY//w4UOcO3cOCoUCLVq0QOfOnYs8B9kJKnZ27dolAIj3339fzJ8/X3z++edi/PjxolGjRqJ8+fLis88+EwqFwtBpSmrevHkCgDAyMhKRkZGGTof+BwwbNkwAEJaWlkKpVOa6jru7u2jXrp3Wsp9++kkAEM+ePSuKNCWzbds2AUBERUVpLW/evLnw8PCQLa7c8+fm2LFjokWLFsLFxUU0bdo018dNpdPvv/8uAAgA4ty5c7mus2TJEmFqaioSExM1ywIDAwUAsXfv3qJKVRIxMTECgNi8ebPW8pkzZwpzc3PZ4so9f27OnDkj6tevL7p16yamT58uxo0bJypWrChcXFxEcHBwkeZCRSs0NFRT14sWLcp1nStXrggA4vjx41rLLS0txfTp04siTUk5OTmJcePGaS07ePCgACD8/f1liSn3/HlRq9Vi3rx5wtLSUgwdOlTMnj1buLu7i44dO4r09PQizUVu/IqwGPvggw+wdOlS/PTTT9i8eTNu376N3bt3Y9u2bejRo4fWLmUlWUZGBrZt24bRo0fD0tIS3t7ehk6JSrnnz5/D19cXH3/8MVJSUrB79+5c14uMjEStWrWKODt5REREAECpeTz5sbGxwaZNm3Dv3j24u7sbOh0qQhs3bkSvXr3g6OiIjRs35rpOZGQkKlSogPLlyxdxdtKLjIwE8L9R1xUrVsSFCxfg5+eH3377DVu2bEFoaKjmV1AqvTZs2IAKFSpg2LBh2Lx5M7KysnKsU5pqITMzE9HR0aXisejiq6++gre3N/7991/s3bsXK1euxMmTJ/HTTz9BCGHo9CTFXc1LmL59++KXX37BuHHj8PPPP2PBggVa92fvpvH8+XPUqFEDnp6esLa21tz/zz//4Pz585g9ezYSExPh6+uLlJQUdO/eHU2aNMkR78mTJzh37hyePXuGWrVqoUuXLrCzs8ux3rvi5ufQoUN4+vQp5s+fD2tra2zevBn/+c9/NLvCEklt69atMDIywpIlS/Dw4UNs2LBB6xCH7MM9nj17hrt372Lp0qUAAGdnZ/j5+QEAfvnlF81hHx9++CFq1qyp2T42Nhbnzp3DixcvULNmTXh6esLKykpz/5UrV3Dx4kXMnTsXr169wtGjR5GQkIDZs2fnmq8u6/v7+yMoKAiZmZlo2bKl1i5aS5cu1eT9448/AgA6duz4zt248pvzTf/++y+CgoIghEDbtm3RqlUrAMDq1asRHx+P1NRUzXNoY2ODqVOnarbNyMjA6dOnce/ePVhYWKBTp05o3Lix1vwbNmyAnZ0dBg8ejMDAQPzzzz+oV68ePDw8cs2nQ4cO+T4uKp2uXbuG4OBg+Pr6IiwsDIsXL8avv/6KihUrAgCSk5Px66+/4tq1a1p/k/Xq1dPUh4+Pj+ZLqi5dumj9LSmVSpw+fRoPHjyAjY0N3N3dter+4cOH+PPPPzFs2DDUqFEDR44cQVRUFMaMGYOqVavmyFeX9cPDw3Hx4kUolUo4OTnBw8MDZcuWBQDs3LkTp06dAgDs3r0b//77LxwcHPDBBx/k+zzlN+eb7t+/j4sXL0KhUKBhw4Zwc3ODiYkJDhw4gKtXr0KtVmueQwA5DucICAjQrNe0aVN069ZNa9f/kydPIiwsDLNmzcKDBw9w4sQJGBkZ5Xm42XvvvZdjWZUqVeDq6oojR45ACJHroQVUsqWnp2Pbtm0YP348RowYgfbt2+PEiRPo06ePZp0NGzbgxIkTAABvb2/Y29vD0dERUVFRyMjIQFBQkOZvtU6dOhg5cqRmW7VaDT8/P4SGhsLExAQdOnTI8bf2448/ok2bNujevTv8/Pxw48YNdOjQAW3bts0153et//z5c5w5cwaPHj2Cvb09evfuDXt7ewDApUuXcO7cOWRlZeHKlSuavOfPn5/v85TfnG96+fIlzp49i4cPH6JWrVpwc3ODra0trl69igMHDgAAtm/fjvPnzwN43Xc0bdpUs/39+/fh5+cHpVKJOnXqwMPDQ+sQ2OzPTxMmTICFhQUOHz6Mhw8fYvr06bC0tMyRT0xMDJYvX47ly5ejUaNGWve5urrm+5hLIv7iXQKNHj0a1tbW2LVrl9byb7/9Fk5OTti+fTvu37+PFStWoG7durh06ZJmnfPnz2PBggW4cOEC+vbti6tXr+LgwYNo1qwZNm3apDXfn3/+iTp16mDbtm2Ijo7Grl270L59e/j4+OgdNz8bNmxAx44d0axZM0ybNg2PHz/G0aNHC/bkEOlg48aNGDp0KOzt7TFt2jQEBgbi5s2bmvszMjKQkJAAIYTm3wkJCUhJSUFqaiqA18eSZi9Xq9WabRcuXAgnJyfs2LED9+/fx48//ghnZ2f4+/tr1jl37hwWLFiAU6dOwdPTE//88w+2bt2aZ775rf/06VN07twZffr0gb+/P0JCQjBkyBD06tULSUlJAICEhASoVCrNvxMSEjSPIze6zAkAjx8/RufOneHm5qY5pnTq1KmaDzUKhQJqtVrrOVQqlZrt79y5g0aNGmHy5Mm4ffs2Tpw4gebNm2PixIlaz+mKFSuwbds2zJkzB9988w2CgoJw/PjxfP6H6X/Rhg0bUKtWLfTt21dzsrTt27dr7hdCICEhARkZGZp/Z9d1YmIigNfNefby7JoBXjfktWvXxvfff4+IiAgcPHgQzs7OWLNmjWad6OhoLFiwAOfOnYObmxsOHz6s+dCZm/zWz8zMxCeffIJGjRrh0KFDiIiIwIIFC9CoUSPcuXMHAJCUlITk5GTNvxMSEjSPIze6zJm93uTJk1GvXj3s2rULERERWLVqFVq3bo3nz58jOTkZKpVK6zlMSEjQ/AqZlJQET09P9OjRA5cuXUJQUBCGDBkCV1dXPH36VBPn4MGD+Pbbb7Fv3z4MHToU165dy/G55l2SkpIQHBwMV1dXNt2l1IEDB/Dy5UtMnToV7dq1Q8uWLbFhwwatdRITEzXvaYmJiUhISNDUshACKpVK83eaXTPA6/egJk2aYNy4cbhx4wb++ecfdO7cGePHj9f6Vf2rr77C4cOHMWzYMPz6668ICAjI9zNufuuvX78etWvXxs8//4z79+9j8+bNcHR0xP79+wG8PgmqQqHQ/Ds77/x++X3XnNnWrl0LBwcHeHl5ITw8HLt370bbtm1x6dIlqFQqzXOT/RwmJCQgPT1ds/2CBQtQr1497N27F3fv3sUXX3yBunXr4p9//tGsc+PGDc3rWo8ePXDu3Dns2rUrz9emffv2Qa1Wo3///rh48SJ+/vlnrF27FqGhoXk+3hLNYDu5U56yj/HO71gzV1dXYWxsLNRqtRBCiI0bNwoAYsuWLZp1srKyxAcffCAqV66sOZbNy8tLABDDhw8XKSkpmnUHDRokKleuLFJTUzXLqlatKiZPnqwVV6lUiqtXr2pu6xo3LzExMcLY2Fj8+eefmmUdO3YU/fr1y3c7ooL6+++/BQBx6dIlIcTrY4tq164tPv300xzr1qhRQ4wYMUJrWX7HeK9du1YAENu3b9csU6vVYuTIkaJq1aoiOTlZCCHE999/LwCIwYMHa+owPj4+z5zzW79z586ievXqIi4uTrN+XFycqFSpkpg2bZpm2dy5c0VuL/m5HYOt65yurq6iWrVqOc7LcOHChXznz35eGjZsKJydnbWey3379gkAYtmyZZpl9evXF5UrVxYrV67ULMvv+XpT9uPmMd6lW3JysrC2thY//PCDZtkHH3wgmjRpkmPdDz74QFSpUkVrWX7HeN+8eVOUKVNGTJ48WfOeK8Tr48mNjY1FYGCgEOL/XlsaNmwooqOjhRBCZGZmihcvXuSac37rf/PNN8LY2FicPXtWs75KpRLdunUTTZo0EVlZWUIIIQ4fPiwACD8/P625czsGW9c5v/zyS2FkZCSOHj2qtf2tW7fE8+fP85w/26RJk4SpqanmeRFCiIiICGFrays8PT01yz755BNhbm4uPvroI5GZmSmE0K2u9+7dK7y8vMTcuXNF/fr1Rf/+/UVsbOw7t6OSyc3NTfTs2VNze8OGDcLU1FQ8fvxYa70NGzYIACIsLExreV7HeKempgpHR0fRsmVLrfMmXb9+XZiZmWm935iYmIjKlSuLQ4cOaZbl97ea1/qnT58WRkZG4ttvv9Vaf8GCBcLCwkLExMQIIYRITEwUAMTixYu11svtGGxd5zx+/LgAIL788ktNrQshxPPnz8WtW7fynD/bjh07BACxevVqzbKkpCTRvn17UblyZc25crJ7GFdXV817e1JSkubzz9vGjh0rAIju3bsLR0dHMWXKFOHp6SmMjIy0Pm+UFmy8iyFdGm8PDw8BQCQkJAghhGjYsKFo2bJljvXCwsIEALFz504hxP813idPntRa76+//hIANCcoycrKEpaWlmL48OFaHzTepmvcvCxatEjY29sLlUqlWbZjxw5hYmIiHj58mO+2RAUxduxY0bx5c61l//3vf4Wtra3Wl1FC6N94u7i4iLZt2+ZYHhISIgCIv/76Swjxf430kSNHdMo5r/UDAgIEAPHLL7/k2Gb+/PnC0tJSZGRkCCF0b7x1nfPy5csCgPj555/zzT2vxtvPz08AEGvXrs1xX9u2bUWdOnU0t+vXry/s7e0LdJIVNt7/GzZt2iTKlCkjnjx5olmW/Tf69odIfRvv7Eby7ZMwqtVqUblyZc2Hw+xG+vPPP9cp57zWz8jIENbW1mLw4ME5tsn+8HzlyhUhhO6Nt65zqlQqYWlpKQYNGpRv7nk13qmpqcLc3FyMHDkyx31ffPGFAKD5kuGTTz4RAMT169fzjfW2LVu2iC+++EKMGTNGVKpUSXTv3l3cuXNHrzmoZIiMjBRGRkbiwIEDmmXJycnC1tZWeHl5aa2rb+Od3UgeO3Ysx31DhgwRjRo10tw2MTERrVu31jnvvNb39PQUlStX1nzRlC0pKUmYmJiIJUuWCCH0a7x1nbNXr16iSpUqWp+335Zf4921a1dRo0aNHHGOHTum9QNcdg/zZoOen4EDBwoAolGjRlqvscuWLRMAhLe3t07zlBQ8iLaESkpKgpGREaysrJCWloawsDC0adMGy5cv1+yOIoRAZmYmAODevXta2799DGW1atUAAHFxcWjRogWMjIwwf/58LFq0CIGBgfD09ETXrl3h7u4OW1tbAChQ3DdlZWVh8+bNcHBwwMqVKzXL09PTkZWVhU2bNmHhwoWFeJaItCkUCuzduxcdOnTQOjbxxYsXSEhIwL59+zBmzJgCzZ2UlITw8HC0a9cuRz1k76r1dj28edyULt5ePygoCMDrY6qyY2bHvXfvHpKTkxEXF4fatWvrHEPXOW/cuAEAaN26tV6PIVv2bmTNmzfPcV+LFi2wfv16pKamai6P1KBBA5iZmRUoFpV+GzZsQI0aNbB582at5eXLl8eGDRvQvn37As8dFBQEW1tbrFu3DgA0NSGEQJkyZSSv64iICCiVSiQlJeV4LYmPjwfwuhb1Of5R1zltbW2RnJxc4Lq+f/8+VCpVnnUNvK797NckY2PjHMd1vsu4ceM0/3716hVcXV3h7u6OO3fulKhLPNK7bdy4EcbGxggNDcXdu3c1y2vVqoWNGzdi/vz5BT7EIPu9zt/fH7dv39Z6r3v27JnkdZ0d08rKCqtWrQKg/VpSrly5fD8z5/c4dJnzxo0baN68OcqUKaN3DOB13bZs2RImJiZay9+s6zfp+nxl1+z06dO1zoUze/ZsfPvtt9i5cycmTJhQoJyLIzbeJVBWVhZu376NRo0awdjYWHMspBACz58/z7H+/Pnz0a5dO61lb5/gIPtEZhkZGZplCxcuRP/+/eHr64tLly5h06ZNMDMzw/r16zFy5MgCxX3TiRMn8PTpU4wePRoJCQla9w0YMACbNm3C119/zetzkmR27NiB8uXLo3Xr1lp/c2ZmZujWrRs2btxY4MY7ux6ysrLyrIc2bdpoLcs+6ZOu3l4/O2ZaWlqOmPXq1cP8+fO1TnqiC13nzD7+raD1mf0BJ7fts5dlrwPo/1zR/47Q0FD4+/tj5syZOd5LBg8ejD179mD16tVaH+r0oVarYWJikmtdf/DBB6hTp47WMqnqWq1W54hpamqK+fPno0GDBnrF0HXOoq5rKyurQn2hVqFCBUyaNAmff/45AgMD0aVLlwLPRcVLZmYmtmzZggEDBuQ4PtjT0xMrVqzA+fPn0b179wLNn10TSUlJWscxA69P6uXq6oqsrCzN321h6zo7prGxca6vJdOmTcv1C6t30XXONx9LQQghdK5rQPfnK/v18+0zuJuZmaFatWp49OhRAbItvth4l0D79+/XOqOxpaUlateuDVtbW61f8aTQokULzbdZCQkJ6NGjB2bMmIGRI0cWOm72SdVy2zYqKgpOTk44c+YMLwdEktmwYQP69euX69+cn58fevTogbt376J+/fp5zpHXG5eNjQ1q1KgBOzs7yeswL9l7rri5ub3zTMZSz5l9FYTr16/n+2tiXs9XvXr1AAC3b9/OcWbYW7duoUaNGnp/aUD/mzZs2ABHR0esXr06x32pqan466+/sGvXLkyePDnPOfL7QNq4cWNERETghx9+KJKrbTg5OaFs2bJwdHSU7LVE1znT0tJgYWGB69ev5ztfXs+Xo6MjzMzMcPv27Rz33bp1C8D/1b5Usr8s4MnVSpejR4/i8ePHWL16da6X1fL398eGDRve2Xjn9bea/V43YsSIfH8kklLjxo0RFxcn6WcEXeds0qQJbt26hczMzDxfx/J7HaxXrx7CwsJyXD2gsHXdtWtXLF26NEeDrVar8eTJk0LtrVQc8afEEuby5cuYNm0amjRponUpodmzZ+Ps2bM4fPhwjm3CwsJy/SYsPykpKbh27ZrWMltbWzg5OWl2Iy9M3CdPnuDIkSPo3bt3rvc7Ojqifv36Oc5cSVRQV69exfXr1/P8m+vUqROsrKzyvPZvtipVqgB4vSva22bPno2TJ0/mesbt0NBQvHz5sgCZ561Lly5o2bIlvvvuuxxzZ1+KRK45O3fujBYtWmDp0qV48uSJ1npvfmivUqVKrs9Vt27dUKdOHaxatUrrTOlnz57FpUuXNGelJsqPSqXC9u3b86xrCwsLdO3a9Z3vJfnV9WeffYbU1NRcD3169uyZ1hnBpWBhYYEpU6Zg+/btOd6HgdeX73vzjOtSzlm2bFlMmjQJ+/bty3HW5piYGM0eBVWqVIFKpcrxS2S5cuUwfPhw/PXXX1q7Bj9+/Bjr1q1Dly5d4OzsrFfu2d48c3I2hUIBb29v1KhRI89LO1HJtGHDBjRs2DDPa1l7enpqznien7zeg0aMGIFq1arhyy+/zHGVD5VKhX///bfgyedh1qxZiIiIwG+//ZbjvtjYWMTExMg254wZMxAXF4effvpJa53k5GTNJRTzex386KOPEB0djW3btmmWZWRkYMmSJbC1tcWQIUP0zh0AevbsiYYNG2LdunVar2sbN25ESkoKRo8eXaB5iyv+4l2MZV9PVK1W48WLF5prYo4YMQIrV67U2m1u5syZiIuLw+DBg+Hh4YHmzZsjOTkZN2/exLNnzzTXN9SVWq3G1KlTYWpqitatW8POzg7BwcE4ceIE1q5dW+i4W7ZsQWZmZp4floDXL6q///47nj17luu1CIn0sXHjRpiYmKBXr1653m9mZgY3Nzds27YNS5YsyXP3x169esHGxgYTJ06Ep6cnzMzMNNfxnjNnDh49eoQBAwbA09MTzZo1Q2JiIm7evIkXL17g5MmTkj4mY2Nj+Pj4YPDgwXBxccHgwYNRo0YNPHjwAJcvX0avXr30vqa1rnMaGxvj4MGD6N+/P5o0aYLBgwfD3t4eV69ehampqeaSgMOGDcPEiRMxefJkODo6wtbWFlOnToWZmRn27t2L/v37o1WrVhg4cCCeP3+OnTt3YuDAgfj6668L/LxERUVhz549AKD58LRmzRpUqFABFSpUwCeffFLgual4OXjwIF68ePHO95JZs2ZpjnHMTY0aNdChQwd4eXnh0aNHsLS01FzHu127dti+fTs++eQTnDlzBl27doWZmRnu3LmD4OBgeHt7673r97ssW7YML168gKurK95//300aNAAL168QFBQEIQQmmvsyjHnsmXLEBcXh+7du2PIkCFwcXFBREQEbty4obne+cCBA7Fw4UIMHToUXbp0gYmJieY63j///DMiIiLQvn17jBo1CmXKlMHevXthZ2en9aFdX97e3pgyZQratWun2QX1wIEDsLKygq+vL8zNzQs8NxUvcXFxOHHiBGbMmJHnOp6enliwYAG2b9+OmTNn5rnesGHDsGrVKsyaNQtVq1bVXMfbysoKx48fx5AhQ1C/fn30798fdnZ2iIqKwqVLlzBt2jTJv8x5//33sXr1anz++efYt28f2rdvD7VajdDQUNy7dw979+7V65ws+sw5aNAgeHl5YeHChTh16hRcXV3x/PlzXLhwAX/88QecnZ3RokULODs7Y/78+bh27RrKli2ruY73xIkT8e+//2LChAk4ceIE6tSpg5MnTyIyMhJ79+5FhQoVCvScmJqaYu/evejduzfatGmDfv364cGDB9i1axemTJlSqo7vBgAj8fZO+WRwN2/exM6dOzW3y5QpAxsbG7i4uKBjx46ws7PLc9vo6GicPHkScXFxsLe3R/PmzdG5c2fNbiEXLlzA8ePH8e2336Js2bKa7WJjY/H7779j7NixWic6+ffffxEQEIBXr16hZs2a6NevHypXrqx33LetXbsWjx49wnfffffO52HEiBF477338n7CiHTg5eUFc3NzzJkzJ891sutjypQpqFOnDry8vFCnTh2MGjVKa73o6GgcOnQIT58+1XxJ9eabZVRUFE6ePIlHjx6hcuXKmnrI5ufnh5MnT+K7777T6UQn71pfCIGLFy8iICAAKSkpcHR0RLdu3bRyOnbsGC5evJhjd7RffvkFVlZW+Oijj/SeE3j9Jd3JkycRHBwMU1NTtGnTBj169NBa59y5c/D390dSUhJsbW0xf/58zX1JSUnw9fXFvXv3YGFhgU6dOqFTp05a269cuRL29vY6H39/7949bNq0Kdf7KlWqhM8//1yneaj427NnD4KDg/HNN9+gfPnyua6T/f7m7u6OHj16YM+ePYiIiMjx5Y5SqcT+/ftx//59ZGRkwMPDQ2s31oSEBBw9ehQREREoV64cGjRoAHd3d03DFxMTgzVr1mD8+PE6NeK6rB8aGoqzZ8/i2bNnqF69Olq3bq114rM7d+5gy5YtmDx5MpycnDTLjxw5goCAAHz//fd6z5nt6tWrOH/+PFJSUtCwYUMMGDBAq7kNCwvDsWPH8OLFC2RlZWH+/PmaD99ZWVk4efIkrl27BrVajaZNm6J///5aX2j6+Pjgxo0bWLx48Tufq2xRUVHw8/NDTEwMLC0t0bx5c/To0YMnXixlAgMDsX//fowePRrNmjXLc71FixahVq1a+PjjjzXbzJ07V+vHmqysLBw+fBjXr19HWloanJ2dtfaoyszMxKlTp3D9+nWo1WrUrVsXvXr10prj66+/RseOHdGnTx+d8n/X+k+fPsWxY8cQHR0NW1tbNGrUCD179tScuCw9PR2LFi1Cz5490bNnT812oaGh2LZtG6ZPnw4HBwe95swWExODY8eO4cmTJ6hduzb69u2r+aUbeP1r94EDBxAbG4vMzEwMGzYMrVq10twfFBSEc+fOQalUok6dOhg4cKBWX5L92X3mzJmaEzfrIikpCT4+PggPD0f58uXh5uaGli1b6rx9ScHGm4iIiIiIiEhGPMabiIiIiIiISEZsvImIiIiIiIhkxMabiIiIiIiISEZsvImIiIiIiIhkxMabiIiIiIiISEZsvImIiIiIiIhkZFrUAbOysvDo0SNYWVnleY1nItKNEAKJiYmoXr06jI0N9z0a65pIOqxrotKpONQ265pIWvrUdZE33o8ePcpx0XciKpzY2FjUrFnTYPFZ10TSY10TlU6GrG3WNZE8dKnrIv+6zcrKqqhDEpV6hq4rQ8cnKo0MXVeGjk9UWhmytljXRPLQpbaKvPHmbi1E0jN0XRk6PlFpZOi6MnR8otLKkLXFuiaShy61xZOrEREREREREcmIjTcRERERERGRjNh4ExEREREREcmIjTcRERERERGRjNh4ExEREREREclI7+t4v3r1Ctu3b0dMTAxcXFwwduxYlCtXTo7ciKiIqNVq7N69G8HBwbC3t8fo0aN5nU+iUuDChQs4ceIETE1NMWDAALRp08bQKRFRIYWHh2PPnj1QKBTo0KEDBg0aZOiUiEgHev3i/fTpU7Rs2RJ79+6Fra0t1q9fj/bt2yMpKUmu/IhIZkIIDBgwAAsXLkT58uVx5coVNGnSBDdv3jR0akRUCMuWLUPfvn2RmZmJhIQEdOrUCZs3bzZ0WkRUCBcuXECzZs1w584dlC1bFlOnTsVHH31k6LSISBdCD59++qlo0KCBSEtLE0IIkZCQIOzt7cV///tfnedQKBQCAAcHh4RDoVDoU8pa9u3bJ0xMTERkZKRmmbu7u+jVqxfrmoPDgKMwdf3o0SNRpkwZsXXrVs2y//73v8LGxkYkJyezrjk4DDgKU9tNmjQRH330keb25cuXBQBx+fJl1jUHhwGHLnWt1y/evr6+GDZsGMzNzQEANjY26N+/P3x9ffWZhoiKEV9fX3To0AFOTk6aZR9++CHOnj3LvVmISqjjx48DAIYOHapZNmbMGCgUCly4cMFQaRFRIdy/fx+3bt3Chx9+qFnWoUMH1K1bl5/FiUoAnRtvlUqF2NhYODo6ai13cnJCeHh4vtsplUqtQUTFR3h4eK51nZWVhfv37+e6DeuaqHgLDw9H1apVtc7B4uDgADMzszzfs1nXRMVbdu2+/Z7t6OjIuiYqAXRuvFNTUwEAVlZWWsutra2RkpKS53ZeXl6wsbHRDJ6wiah4SU1NzbWuAeRZ26xrouItt7oGXr+Hs66JSqaCfBZnXRMVHzo33paWljAyMkJCQoLW8levXmk+pOdmwYIFUCgUmhEbG1vgZIlIeuXLl8+1rgHkWdusa6LiLbe6FkJAoVCwrolKqPLlywOAXp/FWddExYfOlxMzMzODi4sLwsLCtJaHhYWhUaNGeW5nbm6uOSaciIqfxo0b4+rVq1rLwsLCYG5urnXc95tY10TFW+PGjfHkyRO8evUKFSpUAADcu3cParU6z/ds1jVR8ZZdu2FhYXB2dgbw+gu1u3fvokuXLrluw7omKj70OrnaiBEjsGfPHrx8+RIAEBMTgyNHjmDEiBGyJEdE8hsxYgSCgoIQEBAAAEhPT8fGjRsxcOBAlC1b1sDZEVFBeHp6oly5cli/fr1m2e+//44aNWqgY8eOBsyMiAqqevXq6Ny5M9auXQshBADg4MGDePLkCYYNG2bg7IjonfS5hEFSUpJwdXUVderUEaNHjxbVqlUTffv2Fenp6TrPwcsYcHBIPwpzaRIhhJgxY4awsbERI0eOFE2bNhVOTk4iNjaWdc3BYcBR2LretWuXsLCwEP369RM9e/YU5cuXF6dOnWJdc3AYeBSmtm/fvi2qVq0q2rZtK4YNGyYsLS3FDz/8wLrm4DDw0KWujYT4/1+Z6UitVuPMmTN48OABXFxc0LVrVxgZGem8vVKphI2NjT4hiegd8jtuU1fXrl1DcHAw7Ozs0Lt3b1hYWOi8LeuaSHpS1PWDBw9w7tw5mJqaomfPnqhatarO27KuieRR2NpWKpU4efIklEol2rdvj8aNG+u1LeuaSHq61LXejXdhseCJpCfFB/TCYF0TSY91TVQ6GbK2WddE8tClrvU6xpuIiIiIiIiI9MPGm4iIiIiIiEhGbLyJiIiIiIiIZMTGm4iIiIiIiEhGbLyJiIiIiIiIZMTGm4iIiIiIiEhGbLyJiIiIiIiIZMTGm4iIiIiIiEhGbLyJiIiIiIiIZMTGm4iIiIiIiEhGbLyJiIiIiIiIZMTGm4iIiIiIiEhGbLyJiIiIiIiIZMTGm4iIiIiIiEhGBW68MzMzIYSQMhciMjC1Wo2srCxDp0FEEsrKyoJarTZ0GkQkISEEMjMzDZ0GEelBr8Y7MTERa9asQbNmzWBmZoYdO3bIlRcRFRG1Wo1Dhw6hd+/eMDMzw+TJkw2dEhFJ4OrVq5g4cSLKly+P+vXrGzodIpJAbGwsFi5ciOrVq8PMzAwPHz40dEpEpCO9Gu/t27fjxo0b2LZtm1z5EFERCw4OxsaNGzFjxgx06tTJ0OkQkURmz54NV1dXTJ8+3dCpEJFEvv/+e5iZmWHFihWGToWI9GSqz8rTpk2TKw8iMpDWrVvj0KFDAIAff/zRwNkQkVT+/vtvAMC3335r2ESISDLr168HAJw/f96wiRCR3nhyNSIiIiIiIiIZ6fWLd0GoVCqoVCrNbaVSKXdIIpIZ65qo9GFdE5U+rGui4kP2X7y9vLxgY2OjGQ4ODnKHJCKZsa6JSh/WNVHpw7omKj5kb7wXLFgAhUKhGbGxsXKHJCKZsa6JSh/WNVHpw7omKj5k39Xc3Nwc5ubmcochoiLEuiYqfVjXRKUP65qo+NCr8RZCQK1Wa25nZWUhMzMTxsbGMDbmedqISqrMzEwAr2tcCIHMzEwYGRnBxMTEwJkRUUGp1WoIIZCVlQXg/+rc1FT279yJSCZZWVnIysrSfB5Xq9XIzMyEiYkJjIyMDJwdEeVL6OHs2bPCxMQkx5g8ebLOcygUCgGAg4NDwqFQKPQpZS2ZmZm51nX9+vVZ1xwcBhyFqWshhOjatWuutf3w4UPWNQeHAUdhanvRokW51vWOHTtY1xwcBhy61LWREEKgCCmVStjY2BRlSKJST6FQwNra2mDxWddE0mNdE5VOhqxt1jWRPHSpa+4fTkRERERERCQjNt5EREREREREMmLjTURERERERCQjNt5EREREREREMmLjTURERERERCQjNt5EREREREREMmLjTURERERERCQjU0MnQESlR1Fdm9TIyEj2GET0GuuaqPRhXRMVPf7iTURERERERCQjNt5EREREREREMmLjTURERERERCQjNt5EREREREREMmLjTURERERERCQjNt5EREREREREMmLjTURERERERCQjva/jHR0djUuXLiEjIwOtW7dG06ZN5ciLiIpQQkICzp8/jydPnqB+/fro1q0br71JVMJlZmbiwoULuHfvHqpXrw53d3dYWFgYOi0iKqRr164hKCgIVlZW6NKlC6pXr27olIhIB3r94j1p0iT06tULJ06cwPnz5+Hq6opPP/1UrtyIqAisX78eLVq0wNatW3H9+nWMGzcOHTt2RFJSkqFTI6IC8vf3R6NGjbB06VLcvHkTP/zwA5ydnRESEmLo1IiogF68eIGOHTtiypQpCAwMxI4dO+Di4oI///zT0KkRkS6EHo4cOSLUarXm9pUrVwQAce7cOZ3nUCgUAgAHB4eEQ6FQ6FPKWvz8/IRSqdTcfvnypahSpYpYtGiR3nVdmDz0Yejnm4OjKEZh6unGjRsiOjpaczsrK0t4eHiILl266DwH65qDQ55R0Jp69uyZuHLlitayJUuWiLJly4qUlBSd5mBdc3DIM3SpKb1+8e7bty+Mjf9vk/bt26NMmTIIDw/XZxoiKka6desGKysrze0KFSqgWbNmrGuiEqxZs2aoXbu25raRkRG6dOnCuiYqwSpVqgRXV1etZV27dkVaWhpiY2MNlBUR6UrvY7zfdPToUaSnp6NNmzZ5rqNSqaBSqTS3lUplYUISkcyePHmCK1eu4IcffshzHdY1UcmSlZUFHx8ftG3bNs91WNdEJc+BAwdQsWJF1K1bN9f7WddExUhBdx2Jjo4WVapUERMmTMh3vcWLFxv8p38OjtI+pNplTKVSia5du4omTZrku9taXnXNXdc4OKQbUtbTl19+KSwtLUVoaGie67CuOTiKZkhVU0ePHhUmJiZi27Ztea7DuubgKJqhS00VqPGOi4sTzs7Oom/fvkKlUuW7blpamlAoFJoRGxtr8CeGg6O0DSneQNPT08XAgQOFk5OTiI2NLVBd842cg0O6IVU9LVu2TFhYWIgzZ87kux7rmoOjaIYUNeXn5yfKlSsn/vvf/+a7Huuag6Nohi41pfeu5o8fP0b37t1Rr1497N+/H2XKlMl3fXNzc5ibm+sbhoiKUEZGBkaMGIGbN2/i/PnzqFmzZr7rs66JSoYVK1Zg8eLFOHToENzc3PJdl3VNVDJcvHgR/fr1w5dffomvvvoq33VZ10TFh14nV3vy5Am6d+8OZ2dnHDhwgIVMVApkZmZi5MiRuH79Os6fP49atWoZOiUiksCqVauwcOFC+Pr6olevXoZOh4gk8Pfff6NPnz6YP38+Fi5caOh0iEgPev3i7e7ujri4OIwZMwa//vqrZnmHDh3QoUMHyZMjIvnNmTMHBw4cwGeffYa//vpLs9zBwQEjRowwYGZEVFD79+/HnDlz0LdvX9y8eRM3b97U3Pf5558bMDMiKqiYmBj06dMHDg4OsLCwwPLlyzX3jRgxAg4ODgbMjojeRa/Gu2/fvsjIyMCrV6+0liclJUmaFBEVnWbNmmHu3LkAXu/Vks3S0tJQKRFRIVWsWDHXuiaiksvIyAiffPIJgJx1nZ6eboiUiEgPRv//xAdFRqlUwsbGpihDEpV6CoUC1tbWBoufXddFlYeRkZHsMYgMjXVNVDoZsrZZ10Ty0KWm9DrGm4iIiIiIiIj0w8abiIiIiIiISEZsvImIiIiIiIhkxMabiIiIiIiISEZsvImIiIiIiIhkxMabiIiIiIiISEZsvImIiIiIiIhkZGroBIiIiIiEEEUSh9cVJio6RVXXAGubij/+4k1EREREREQkIzbeRERERERERDJi401EREREREQkIzbeRERERERERDJi401EREREREQkIzbeRERERERERDJi401EREREREQkI70a79TUVKxatQrdunVDkyZNMHjwYFy4cEGu3IioCAghcODAAfTv3x+NGjVCz549sWnTJkOnRUSFFBQUhLFjx6Jp06bo0KEDvv32WyQmJho6LSIqhMePH+Pzzz9H69at0apVK0yePBn37983dFpEpAO9Gu+vv/4aL1++xH//+1/s2rULDRs2RM+ePeHv7y9XfkQks127dsHX1xfTp0/H/v37MWbMGMyYMQM//fSToVMjogKKjY3F7Nmz0adPH+zZsweLFi3Cnj17MHLkSEOnRkSFMG7cONSpUwcbN27E77//jqdPn6JTp05QKpWGTo2I3sFICCF0XVmtVsPExERrmYODAyZNmoRFixbpNIdSqYSNjY1+WRJRvhQKBaytrQu0bW51PWbMGMTFxeHcuXM6zZFd14XJQx9GRkayxyAytMLUU1ZWFoyNtb9b9/b2xpQpU6BSqXLcl5uiruuiwtcPMjQp37MjIyPh7OyMCxcuoEuXLu/cvrTWNcDaJsPSpab0+sX77Q/nFy5cQHx8PDp16qR/dkRULLxd1w8fPsSVK1fQtWtXA2VERIX1dmOdnJyMI0eOoFOnTjo13URUPL35np2VlYW9e/eiYsWKaNy4sQGzIiJdmOq7wb179+Du7o7k5GQkJyfD29sbPXr0yHN9lUoFlUqluc1dYYiKp5YtW+LJkyeIj4/HtGnT8t2LhXVNVDJMmzYNR44cQXx8PNq2bYtDhw7luS7rmqhk2LZtGxYtWoSXL1/C1tYWZ8+ehZ2dXa7rsq6JihGhp/T0dBEVFSWCg4PFggULRLly5cSlS5fyXH/x4sUCAAcHh4xDoVDoW8o5PHjwQISGhoqNGzcKa2trsWLFCr3rWoo8dGHo55uDoyiGFPUUHx8vIiIixOHDh0XDhg3FgAED8lzX0HVdVAz9/8rBUdiaUiqVIioqSly6dEn0799fODo6iufPn+e67v9KXQvB2uYw7NClpvQ6xjs3Hh4esLCwgI+PT6735/ZNm4ODQ2FCEtFbpD5Wy8vLCz/++CNevXqV6/151TWP8SaSjtT15O/vjw4dOiA4OBgtWrTIcb+h67qo8PWDDE3KmkpPT4ednR2+//57zJo1K8f9/yt1DbC2ybB0qSm9dzV/W/ny5fO9PIm5uTnMzc0LG4aIilD58uWRkpICIUSub2Ssa6KSp3z58gCAlJSUXO9nXROVPGXKlIG5uTnrmqgE0OsMK3PmzMHDhw81tw8ePIgjR45gyJAhkidGREXjp59+wrVr15C980tISAhWrVqFwYMH89tjohJq//79OHjwINLT0wEA8fHx+PLLL1GrVi20bt3awNkRUUHcvXsXP/74o+YHL5VKhYULF0KhUGDAgAEGzo6I3kWvxrt169Zwc3ODnZ0dLC0tMXv2bKxYsQKffPKJXPkRkcy6du2KuXPnokKFCqhYsSI6d+6MgQMHYsOGDYZOjYgKqH379vDx8UGlSpVQpUoV1K5dG+bm5jhz5gzKlClj6PSIqABq166N1NRU1K1bF1WqVIGNjQ1OnjyJY8eOoUmTJoZOj4jeoUDHeCckJMDMzAyWlpZ6B+R1vImkJ8WxWqmpqUhNTUXFihX13pbX8SaSnhT1pFar8eLFC9jb2+tdN6X1er98/SBDk6Kmnj17BltbW5iZmem1XWmta4C1TYYl2zHetra2BdmMiIoxCwsLWFhYGDoNIpKQiYkJKleubOg0iEhi9vb2hk6BiPSk167mRERERERERKQfNt5EREREREREMmLjTURERERERCQjNt5EREREREREMmLjTURERERERCQjNt5EREREREREMirQ5cSIiIiISiIhRJHF4nWFiYpOUdU265oKir94ExEREREREcmIjTcRERERERGRjNh4ExEREREREcmIjTcRERERERGRjNh4ExEREREREcmIjTcRERERERGRjNh4ExEREREREcmowI23UqmEj48PAgMDpcyHiAzs+PHjOH78uKHTICIJ/fvvv/Dx8UFiYqKhUyEiidy/fx8+Pj6IiooydCpEpIMCN96TJ0/G8OHDsWzZMinzISIDWrNmDQYOHIgxY8YYOhUiksj169fh5uaGQYMGISYmxtDpEJEEEhMT0adPHwwZMgRHjx41dDpEpIMCNd4bNmxAXFwc3NzcpM6HiAwkJCQES5YswaeffmroVIhIIsnJyRg1ahTmz59v6FSISEJTp07F4MGDYWFhYehUiEhHejfeoaGhWLRoEbZv3w4TExM5ciKiIpaSkoKRI0fil19+QfXq1Q2dDhFJZPr06XB3d0fv3r0NnQoRSWTLli24c+cO/vOf/xg6FSLSg6k+K6elpWHkyJH46aefUKdOHZ22UalUUKlUmttKpVKvBIlIfjNmzICrqysGDRqE5cuXv3N91jVR8bdjxw5cu3YNgYGBuHXr1jvXZ10TFX93797F/PnzcfHiRZiZmb1zfdY1UfGh1y/es2bNQuPGjfHhhx/qvI2XlxdsbGw0w8HBQe8kiUg+e/bswcWLF/Hzzz/rvA3rmqh4i4iIwKxZs7Bjxw6ULVtWp21Y10TFm0qlwsiRI/H999+jfv36Om3DuiYqPnRuvK9cuYLNmzfDw8MDPj4+8PHxwdOnT/Ho0SP4+PggJSUl1+0WLFgAhUKhGbGxsZIlT0SFk5aWhsmTJ2PYsGE4ffo0fHx8cPv2baSnp8PHxwcPHjzIdTvWNVHxNm3aNHTo0EFz1uMLFy4AAM6dO4egoKBct2FdExVvK1asgFKpROXKlTWfxdVqNW7evIlTp07lug3rmqj40HlXcxMTE3h6esLHx0ezLDY2FiYmJtiyZQs6dOiAcuXK5djO3Nwc5ubmkiRLRNJSq9Xo3r07bt++jdu3bwN4/UuZSqXCli1bYG1tjVq1auXYjnVNVLw1btwYUVFR2LJlCwAgISEBAODr6wsAaNmyZY5tWNdExVuVKlXQtGlTTV0DQEZGBgICAgAA7u7uObZhXRMVH0ZCCFHQjfv164eyZcti3759Om+jVCphY2NT0JBElAuFQgFra2tJ5lq+fDmWLl2K58+f67xNdl1LmUd+jIyMZI9BZGhS1tPVq1fRpk0bhISEoEmTJjptU9R1XRrxtYpyI2VNlS9fHkuXLtX5iiSs68JjXVNudKmpAl/Hm4iIiIiIiIjeTa+zmr/N1dVVpzMqElHJ4eLigj59+hg6DSKSUIUKFTBw4ED+wkVUyvTv3x9OTk6GToOIdFCoXc0LgruaE0nP0LuMcVdzIun9r9V1acTXKsqNIWuKdV14rGvKDXc1JyIiIiIiIjIwNt5EREREREREMmLjTURERERERCQjNt5EREREREREMmLjTURERERERCQjNt5EREREREREMmLjTURERERERCQjU0MnQEREREREVBIIIYosFq8ZXrrwF28iIiIiIiIiGbHxJiIiIiIiIpIRG28iIiIiIiIiGbHxJiIiIiIiIpIRG28iIiIiIiIiGbHxJiIiIiIiIpIRG28iIiIiIiIiGel1He/t27fD29tba1n58uVx5MgRSZMioqITEhKCzz77LMfyjRs3wtnZ2QAZEZEU3NzcoFartZZ9+umnGDp0qIEyIqLCmjt3Lq5du6a1rFOnTvjhhx8MlBER6UqvxjsmJgbx8fH4448//m8CU72mIKJiRqFQ4MKFCzh58iTKlCmjWV61alUDZkVEhXXhwgUsWbIEbdu21Szjl2lEJduNGzdQvXp1TJ48WbOsUqVKBsyIiHSld9dsbW2Nbt26yZAKERlSly5dULZsWUOnQUQSatasGd+ziUqZWrVqsa6JSiC9G++oqCj0798fZcuWRdu2bfHpp5/CwsJCjtyIqAiNHj0aGRkZaNSoEWbMmIEaNWoYOiUiKqSlS5di9erVcHR0xLhx49C+fXtDp0REhXT06FFcv34dVapUQe/evTFq1ChDp0REOtCr8TYxMcH777+P3r17IyEhAT/99BO2bt2KwMDAPJtvlUoFlUqlua1UKguXMRFJzt3dHYMHD4a5uTm2bduGBg0aICAgAI0aNcp1fdY1UfHXoEEDDBs2DLVq1cLZs2fRpUsXbNiwAePGjct1fdY1UfFXoUIFdOrUCa1bt0ZYWBimT5+O06dPY9OmTbmuz7omKkaEHpKTk7VuP336VFhZWYkVK1bkuc3ixYsFAA4ODhmHQqHQp5S1pKSkaN3OysoS7du3F/3799e7rguThz4M/XxzcBTFKGw9vf2ePX/+fGFrayvUanWu6xu6rksjQ/8NcRTPUZiaeruujx8/LgCIa9eu5bo+67pkM/TfKoe0dW30//9TC6xr166oXbs2tm3bluv9uX3T5uDgUJiQRPQWhUIBa2tryeZbvHgxNm/ejAcPHuR6f151LXUeeTEyMpI9BpGhSV1Pfn5+6NGjByIjI+Hk5JTjfkPXdWnE1yrKjZQ1lZmZibJly2L9+vWYMGFCjvtZ1yUbX0NKDl1qqtCnJH/27Bnq16+f5/3m5uYwNzcvbBgiKkLPnj3L90RrrGuikufZs2cAkGdts66JSp6XL19CrVazrolKAGN9Vl69ejXS0tI0t3/77TeEhYVh2LBhkidGREVj27ZtePTokeZ2QEAAtm7diuHDhxswKyIqjLNnzyIwMFBz+/Hjx/j+++/RsWNHVK9e3YCZEVFBPXjwALt27dLcTktLw5w5c2BjYwMPDw8DZkZEutDrF+/09HQ4OjqievXqePHiBZKTk+Ht7Y1evXrJlR8RyaxSpUro3r07TExMYGxsjMjISEydOhWLFy82dGpEVEA1a9bElClTcP/+fdjb2yM0NBRubm5Yu3atoVMjogKqUKECTp8+jZkzZ6JOnTq4f/8+qlWrhmPHjsHOzs7Q6RHRO+h9jLdKpUJYWBjKlSsHR0dHmJmZ6RVQqVTCxsZGr22IKH+FPVYrKysL4eHhSEtLg7OzMywtLfXaPruueYw3kXSkqKe4uDg8fvwYjo6Oen8wL+q6Lo34WkW5KWxNKRQK3Lt3D1WrVkXNmjX1+jtjXZcsfA0pOWQ5xtvc3BwtWrQoaE5EVAwZGxvne64GIiqZatSogRo1ahg6DSKSkI2NDdq0aWPoNIhIT3od401ERERERERE+mHjTURERERERCQjNt5EREREREREMmLjTURERERERCQjNt5EREREREREMmLjTURERERERCQjvS8nRkRERERERPISQhRJHF4vvGjwF28iIiIiIiIiGbHxJiIiIiIiIpIRG28iIiIiIiIiGbHxJiIiIiIiIpIRG28iIiIiIiIiGbHxJiIiIiIiIpIRG28iIiIiIiIiGRW48X748CHi4uKkzIWIDOzly5eIiIgosutGEpH8UlJScPfuXahUKkOnQkQSUavViIiIQEJCgqFTISId6d14BwQEoGnTpmjevDm6deuGbt264eHDh3LkRkRF5OnTp+jXrx9q1aqF/v37o379+vj7778NnRYRFUJmZibmzJkDe3t7DBgwAE5OTli7dq2h0yKiQtqyZQuqVauGHj16oGnTpvjkk0+Qnp5u6LSI6B30arwjIyPRs2dP9O3bF/Hx8QgPD8eSJUsQGRkpV35EJLP09HR4eHhApVLh4cOHCAsLw7lz5xAdHW3o1IioED777DPs3bsX165dw927dxEREYG0tDRDp0VEhbBnzx5MnjwZ69evx4MHD/DgwQO0bdsWCoXC0KkR0TsYCT32KZ04cSIuX76M0NBQGBsXbC91pVIJGxubAm1LRLlTKBSwtrYu0Lbbt2/H+PHjER0dDQcHhwLNkV3XhclDH0ZGRrLHIDK0wtRTTEwMnJycsH37dowePbpAcxR1XZdGfK2i3BSmpurWrQs3NzesX7++QNuzrik3fK0qPF1qSq/u+dSpU+jfvz+ysrJw9+5dvHr1qlAJEpHhnTp1Ci1btoSDgwOioqJ47gaiUuDMmTMAgP79+yM+Ph6RkZHIzMw0cFZEVBiRkZG4f/8+Bg4ciKSkJNy9exepqamGTouIdKRz4y2EwKNHj/Dq1SvUq1cP/fr1g4ODA9zc3PL9oK5SqaBUKrUGERUfcXFxsLGxQdeuXdG9e3e0aNECzs7O+R7jzbomKt7i4uJgZ2eHL7/8Ek2aNIGbmxvs7OzyPcabdU1UvGV/3v7777/h6OiIfv36oUKFCpgyZUqeX6yxromKD50bbyMjI5iYmGDv3r04cuQIwsPDERsbC6VSicmTJ+e5nZeXF2xsbDSjoLuyEpE8TExMcPbsWUyaNAnR0dF48uQJ3NzcMGTIECQnJ+e6DeuaqHgzMTHBs2fPkJmZicePHyM6Ohq///47pk2bhoCAgFy3YV0TFW8mJiYAgLNnzyI8PBzh4eG4du0adu3ahVWrVuW6DeuaqPjQa1dzBwcHeHh4oFGjRgCAChUq4KOPPoKfn1+elx9asGABFAqFZsTGxhY+ayKSTK1atWBvb48PP/wQwOs39pkzZ+LZs2e4detWrtuwromKt1q1agEAZs6cqfmw/uGHH6JSpUo4f/58rtuwromKt+y6/uijj2BrawsAaNy4MTw8PHDu3Llct2FdExUfpvqs3KtXL9y5c0drWXx8PCpUqJDnQfnm5uYwNzcveIZEJKtevXph165dSElJQbly5QC8rmsAqFixYq7bsK6JircePXrA2NgY8fHxmi/LU1JSkJSUxLomKqEcHBzQoEEDzXt0tvj4eNSoUSPXbVjXRMWHXo33l19+iZYtW+Krr75Cv379cPv2baxcuRILFy6UKz8iktnQoUOxYsUKjBw5EnPmzIFSqcT8+fPRt29fuLi4GDo9IiqAGjVqYMaMGZg6dSqWLVsGa2trrFy5EnZ2dhgyZIih0yOiAvLy8sL48eNRo0YNNG7cGEeOHMHly5dx8eJFQ6dGRO+g1+XEAODOnTvw8vLC3bt3Ua1aNYwaNQrDhw/XeXteToxIeoW9LEhCQgKWLFmCS5cuwcrKCm5ubpgxYwbKli2r0/a8nBiR9ApbT1lZWfjjjz9w4MABqNVqvPfee5g/fz6qVaum0/a87FDh8bWKclPYmjpx4gT++OMPPHv2DM7Ozpg1axZatWql07asa8oNX6sKT5ea0rvxLiw23kTSM/QbKBtvIun9r9V1acTXKsqNIWuKdU254WtV4Ul+HW8iIiIiIiIi0g8bbyIiIiIiIiIZsfEmIiIiIiIikhEbbyIiIiIiIiIZsfEmIiIiIiIikhEbbyIiIiIiIiIZsfEmIiIiIiIikpGpoRMgIiIiIiIiwxBCFFms/+VrhvMXbyIiIiIiIiIZsfEmIiIiIiIikhEbbyIiIiIiIiIZsfEmIiIiIiIikhEbbyIiIiIiIiIZsfEmIiIiIiIikhEbbyIiIiIiIiIZ6XUd7wULFuDx48c5ljdr1gxz5syRLCkiKjpbt26Fn59fjuVlypTB+vXrDZARERVWQEAA1qxZk+t9ixYtgpOTUxFnRESFpVAoMHPmzFzvGz58OPr06VPEGRGRPvRqvNu2bQuFQqG5nZSUhM8++wyLFy+WPDEiKhr16tWDEEJr2cKFC1G3bl0DZUREhVW5cmV069ZNa9mOHTvwzz//4LfffjNMUkRUKGXKlMlR19euXcNvv/2G8ePHGyQnItKdXo33oEGDtG5v3LgRxsbGmDBhgqRJEVHRcXV1haurq+Z2REQE4uLisHTpUgNmRUSF4ejoCEdHR83trKwsLF68GCNHjkT58uUNmBkRFZSFhUWOBvvMmTNwcXFB165dDZMUEemsUMd4e3t7w8PDA7Vq1ZIqHyIyMG9vb9ja2mLIkCGGToWIJHL69Gk8ePAAkyZNMnQqRCSRhIQEHDhwAB9//DGMjIwMnQ4RvYNev3i/KTQ0FP/88w8OHDiQ73oqlQoqlUpzW6lUFjQkEclMrVZj69atGDNmDMqWLZvneqxropLF29sbzZo1Q9u2bfNch3VNVLLs3LkTmZmZ+e5mzromKj4K/Iu3t7c3qlativ79++e7npeXF2xsbDTDwcGhoCGJSGbHjh3D48eP3/mrGOuaqOR48eIFfH19WddEpYy3tzcGDBiAypUr57kO65qo+DASb59VSQcZGRmoUaMGJk6cCC8vr3zXze2bNhY9kbQUCgWsra0LPc/777+Pp0+fwt/fP9/18qprqfJ4F+5SR/8LpKqn1atX46uvvsKjR49ga2ub53qGruvSiK9VlBspaur69et47733cOLECXh4eOS5HuuaipvS+rqoS00VaFfzQ4cO4fnz5/j444/fua65uTnMzc0LEoaIitDTp09x9OhRrFu37p3rsq6JSg5vb28MHTo036YbYF0TlSTe3t6oXbs2evXqle96rGui4qNAu5p7e3ujR48evNwQUSmydetWlCtXDiNGjDB0KkQkkX///Re3bt3C5MmTDZ0KEUlEpVJhx44d+Pjjj2FsXKjzJBNREdL7F++HDx/i5MmT2Llzpxz5EJGBbNq0CaNHj4alpaWhUyEiiXh7e6NBgwbo1KmToVMhIokcOHAASqUSH330kaFTISI96N14q1QqeHt757imNxGVXMnJyfjyyy/Ro0cPQ6dCRBJyd3fHxIkTDZ0GEUmoVq1aOHDgAGrUqGHoVIhID3o33nXr1uUu5kSljKWlZb6XIyGikmnIkCGGToGIJNaxY0dDp0BEBcADQ4iIiIiIiIhkxMabiIiIiIiISEZsvImIiIiIiIhkxMabiIiIiIiISEZsvImIiIiIiIhkpPdZzQtLCFHUIYlKPUPXVXZ8pVJp0DyIShPWNVHpZMjaZl0TyUOXui7yxjsxMbGoQxKVeomJibCxsTFofABwcHAwWA5EpQ3rmqh0MmRts66J5KFLXRuJIv7aLSsrC48ePYKVlRWMjIx02kapVMLBwQGxsbGwtraWNb+iisXHVDJiFffHJIRAYmIiqlevDmNjwx05wrou2jilNRYf02slua6B4v3cFvdYfEwlI1ZB4xSH2mZdF32s0viYijJWcX9M+tR1kf/ibWxsjJo1axZoW2tra9mf8KKOxcdUMmIV58dkyF/EsrGuDROntMbiYyr5dQ0U3+e2JMTiYyoZsQoSx9C1zbo2XKzS+JiKMlZxfky61jVPrkZEREREREQkIzbeRERERERERDIqEY23ubk5Fi9eDHNz81ITi4+pZMQqjY+puCiNz21pfExFGYuPqXQojc8tH1Pxj1OUsVjXJT9OUcYqjY+pKGOVpsdU5CdXIyIiIiIiIvpfUiJ+8SYiIiIiIiIqqdh4ExEREREREcmIjTcRERERERGRjIr8Ot76SkhIwOXLl2FqaoqOHTuifPnyssWKjo5GQEAA3nvvPdSrV0+2ODExMQgJCUGlSpXQsmVLlClTRpY4QgiEhIQgOjoaNWvWxHvvvQcjIyNZYmULDw/HtWvX0LJlS8mfw8TERBw9ejTH8u7du6NKlSqSxsqmUCjg7+8PCwsLuLq6Sv5/FRUVhYCAgFzv6927N2xtbSWNV1ykp6fjypUrUCgUaNOmDapXry5bLIVCgTNnzqBKlSro1KmTbHFevnyJa9euwdTUFC1atECFChVki5X9GmJra4tWrVrBwsJCtljA6+fw+PHjqFmzpizPoa+vL1JTU7WWNW7cGE2bNpU8FgBkZGTgn3/+QVJSElxdXSWvs+TkZBw+fDjX++R4bSxOgoKCEBMTAxcXFzRp0kS2OJmZmThz5gxSUlIwePBg2eKkpaUhKCgICQkJaNy4MWrXri1brOzXECMjIzRv3hz29vayxcrm4+ODrKwsWZ5Df39/xMTEaC2zt7eHm5ub5LGy3bhxAzExMWjRogVq1aol+fyHDx9GcnJyjuVyvTYWF1FRUbhx4wbs7e3Rvn17mJiYyBbr33//xf3792X9DJT9+TgmJgZ16tSR7b0G+L/XkFevXqFBgwaoW7eubLGyXb58GbGxsbI8h5GRkQgMDMyxfMSIEbL1GA8ePMD169dRq1YttGjRQvL5AwICEBUVlWN5mTJlJH9tLNaN9/HjxzFy5Eg0atQIaWlpiIuLw6FDh9C+fXtJ44SGhmLevHkIDQ3Fo0ePsGzZMlk+GMXFxWHSpEkICwtD48aNERkZidTUVOzduxdt2rSRNNb169cxceJECCHg4OCAq1evws7ODseOHUPNmjUljZUtMTER/fr1Q2RkJJYvXy75cxgXF4dRo0ZhwIABWo1GkyZNZGm8165di3nz5uG9995D+fLlER8fj4MHD8LBwUGyGA8ePICPj4/WsuvXryMyMhJxcXGSxSlOoqOj0atXLwCvP6wEBARg5cqVmDJliqRxEhMTMXfuXBw5cgRZWVlo3769LB+MhBCYNm0aDh06hMaNGyMlJQU3b97E6tWrMWHCBEljvXr1ChMmTEBISAgaN26M+/fv4+nTp9ixY4fmOZXDpEmT4OvrC09PT1mew08++QSOjo5aTY2xsbEsH4YuX76MUaNGwdraGs7Ozpg9ezZ+++039OzZU7IYKSkpOer60aNH+Pvvv7F///5S2XirVCoMGTIEAQEBaNmyJQICAjBw4EBs3rwZxsbS7ly3bNky/PHHHzAxMcHLly9la7zXrVuHJUuWoEaNGrCxscHFixcxduxY/PHHH5J/wPzyyy+xc+dONGnSBElJSbh69Sq+/fZbfPHFF5LGedPvv/+OOXPmwMLCQpbn8Oeff0ZgYKDW55sGDRrI0ng/efIEQ4YMQVRUFNq1a4d79+5h/PjxmDdvnqRxTp06hWfPnmluZ2ZmYv/+/fjss89KbeP9/fffY+nSpXB1dcW9e/dQqVIlnDp1CpUqVZI0zuHDh7Fo0SIkJycjPDwcwcHBsjRZf//9N6ZPnw4jIyPUqlULgYGBcHFxga+vLypWrChprN27d+Prr7+Gg4MDLCwscOnSJXh6emLHjh0wMzOTNFa24OBgeHp6IjExUZbn8PTp05g3bx769u2rtXz48OGSvy5mZmZi2rRp2LlzJ7p27YqEhATY29vj4MGDksYKCgrChQsXtJadOnUKNWvWlP61URRTiYmJws7OTnz99deaZePHjxd169YVarVa0lhXrlwRhw8fFmq1WtjZ2YlVq1ZJOn+2sLAwcfz4cc1ttVotRo4cKVxcXCSPdeXKFXH37l3N7dTUVNGoUSMxYcIEyWNlGzVqlPj666+FjY2NLM9hWFiYACBiY2Mln/ttR48eFUZGRuLw4cNa8d98TuXSvHlz8f7778sex1A8PDxEt27dREZGhhBCiM2bNwtTU1MRHh4uaZynT5+KdevWiaSkJDFkyBAxcOBASefPplarxZo1a4RKpdIs++OPP4SpqamIiYmRNFZcXJw4ffq01rIJEyYIJycnSeO8ad26daJz586if//+sj2HVapUEdu3b5dl7jfFxcUJW1tbMW/ePM2yV69eiQsXLsgee+bMmcLe3l6kp6fLHssQvLy8ROXKlUVcXJwQQojQ0FBhYWEhtmzZInmsX3/9VTx8+FD8+uuvwsbGRvL5s23fvl08e/ZMc/v69evCzMxM7NixQ/JYO3bs0LwmCiHEtm3bBADJX0Oy3bhxQ9SsWVN88cUXsj2HI0aMEBMnTpRl7jdlZWWJdu3aia5du4rk5GQhxOvX5Tffv+Vy8OBBAUBcv35d9liGcOXKFQFA876TlJQkGjduLMaPHy95rD179oigoCAREhIiAIjg4GDJYwghxMmTJ8Xt27c1txUKhXBxcRGTJk2SPJavr694+fKl5vb9+/eFqamp2LZtm+SxhHjdO9WvX18sW7ZMtudwzZo1om7dupLPm5t58+aJKlWqiPv372uWZfdrcoqPjxdmZmZi9erVks9dbBvvvXv3CmNjYxEfH69ZdvPmTQFAXLp0Sba4cjbeudmzZ48AIJKSkmSP5enpKYYNGybL3N7e3qJ169YiPT1d9sZ7z5494tChQ7I2wR07dhQDBgyQbf68XL16VQAQR48eLfLYReHp06fCyMhI7N+/X7MsMzNT2Nvbix9++EG2uHI23rl5+vRpkf0/Llu2TNjb28sy961bt0TVqlVFdHS0GDhwoKyN9zfffCMOHjworl27ptWASOnrr78WlSpV0vqSpCioVCphZ2en1fCXNo0bNxafffaZ1rLBgweLnj17yhZT7sY7Nw0bNiyS/8eAgAABQKtBkEpycrJo2LChOHjwoFi1apWsjfeAAQPEwYMHxd9//y2USqUscU6fPi0AiKCgIFnmz0+/fv1EmzZtijxuUZk+fbpo2rSp1rJffvlFlCtXTrbXUbkb79xMnTq1SP4fMzIyhLW1tfj9999lmX/MmDFi1qxZIjg4WNbG28HBQRw7dkycPHlSPHr0SPIYQgihVCqFhYWFWLlypSzz52f58uXC3NxcvHjxQvK5i+2u5iEhIahSpYrWMU6NGzeGsbExQkJC0LFjRwNmJ50zZ87A0dERlpaWssy/d+9epKamwt/fH+Hh4Tl2f5TCnTt3sGDBAly6dEm2XWeyGRsbY8WKFbC1tcWVK1fQrVs37NixA9bW1pLFSE9PR0BAAFavXo2YmBhcv34d1apVQ6tWrWQ9rgkAvL294eDggN69e8sax1Bu374NIYTWsZ8mJiZo1KgRQkJCDJiZtM6cOQMjIyM0atRIlvkvXLiAuLg4hIeHw9vbG7/99pvkMVJTUzFixAisWLFC1uNas/31118IDg5GcHAwbG1tsXv3bsl3Nb948SK6deuG1NRUnDp1ClZWVnjvvfckff3IjY+PD168eIGPP/5Y1jiGkpmZibCwMMyYMUNredOmTbF27VoDZSW9mJgYREREyHbselhYGG7cuIGnT59i/fr1mDt3riyvIZ9++ik6deqE999/H6tXr5Z8/jcFBAQgMzMT0dHRePLkCdauXYthw4ZJGuPixYuoVKkSmjRpAj8/P6hUKjRr1kzWc4cAwOPHj3H8+HGsWbNG1jiGFBISkuPvvWnTpkhJScH9+/fRoEEDA2UmnczMTJw/f17yw1izvXjxAqdPn0ZSUhL27NmDtm3bYuzYsZLH2b59O4KDg3H16lWEhYVJPv+bXr58iRUrViAtLQ2BgYGYM2cOvLy8JI0RGBiI1NRUuLu7IygoCA8fPoSLiwsaNmwoaZzcbNq0CUOGDJH80AOgGB/jrVAocjxgY2Nj2NraIiEhwTBJSezw4cPw9vbGrl27ZI2RkJCAq1evonPnzpKfrEWlUmHEiBH44Ycf4OLiIuncb6tQoQKuXr2K9957DwDw8OFDtG/fHvPmzcO6deski/Py5UtkZmbi7NmzWL58OZo2bYobN27AysoKR44cQZ06dSSL9abU1FTs2rULM2fOlPyYyOJCoVAAQI7atrOzKzV1HR0djVmzZmHq1Kmy/a34+/vj6tWruH37NqpVqwZHR0fJY8yaNQvNmzfH6NGjJZ/7bZs3b4anpyeA13UwePBgjBw5EiEhIZLWQnx8PExMTNCyZUs0aNAAjx8/xoMHD/Dnn3/K+mWXt7c3unbtWiqP7QaApKQkZGVlleq6VqlU+OCDD9C4cWOMHDlSlhjh4eE4ePAgYmNjkZycjObNm0seY/fu3bh8+TKCgoIkn/ttkyZNwrZt2zQnJf32228xbtw4tGrVCk5OTpLFiY+Ph5WVFTp16gQLCwuYmpri8uXLWLRoERYsWCBZnLdt2bIFFhYWsv09FAcKhSLHF6B2dnYAUGpqe8GCBXj48CG+/vprWeZ/9eoVfHx88PLlS9y4cQMff/wxypYtK2mM8PBwzJkzB+fOnYO5ubmkc7+tTZs2iIqK0vQTZ8+ehbu7O5o1a4ZRo0ZJFic+Ph4A8N133yE0NBS1a9fGpUuX0L17d/z111+y/dDn7++P0NBQ/P7777LMX2x3Nf/iiy9E7dq1cyw3NzeXdVfwotrV3M/PT1hYWIglS5bIHkuI18d4d+jQQfTp00fSeb28vESdOnXErl27NKNcuXJi7NixwtfXV9JYufnuu+9EtWrVJJ3z1atXAoBo0KCB5hCAtLQ00bp1a1mPvd6+fbswNjaW7Zi+4uDYsWMCgIiKitJa7u7uLuuu4EW1q3lcXJxwdnYW/fr1K7JjeefMmSPs7e01xzZKwd/fX5iZmYm1a9dq6rp169aidevWYteuXSIxMVGyWLm5ePGiACD54STNmzcX5ubmIjQ0VLPss88+E5UqVRKZmZmSxsoWExMjjI2NxZ9//inL/MVBSkqKACA2b96stXzJkiWy7gpeVLuap6eniwEDBggnJ6ciOceIEK+PDQUgAgICJJtTpVIJa2tr8fnnn2vqeuzYsaJcuXJi165dIjIyUrJYucnIyBDm5uZi3bp1ks47c+ZMAUB4e3trlu3fv18YGRnJuruys7Oz+Pjjj2Wbvzho27atGDdunNay7OO+5Xpui3JX82XLlgkLCwtx9uxZ2WMJIURUVJSoWLGi+PHHHyWd18PDQ/Tp00dT10uXLhUAxNKlS4W/v7+ksXLTo0cPMWrUKEnnzD5/wtixY0VWVpYQQojo6GhhZWUly7HX2SZOnChcXFw0MaVWbH9Wq1u3Lp4+fQqVSqVZ9uTJE6hUKkm/KTWECxcuoF+/fvjqq69k/Tb2TWXLlsXQoUPx999/Szqvg4MD2rVrBx8fH83IyMhASEgITp48KWms3FhbW2udYVQKtra2sLOzQ+/evTWHAJibm6Nfv364du2apLHe5O3tDQ8PD1kugVJcZF9G48GDB1rLY2JiSnxdP3r0CN27d0e9evWwb98+2Q+7yPbhhx/i2bNnuHPnjmRzmpubY/DgwfDz89PU9ePHj/H48WP4+PjkejkdKWXv+i11bdetWxcNGzbU2lVtyJAheP78eY7LHkll8+bNsLGxwZAhQ2SZvziwsLBAtWrVSmVdZ2RkYMSIEQgJCYGfn59sVwV524ABA2BtbY1Lly5JNqcQAp6enoiNjdXUdUhICDIyMuDj4yNbDWQzNTWFhYWFLHUNQKvGBg4cCGNjY9l+2b9w4QIiIiIwadIkWeYvLurWrZtrXRsZGcmyp1VRWr58ORYvXgxfX1/06NGjSGLWqVMHnTp1kvyzeKtWrWBlZaWpaz8/PwCAn5+frJ9bs8nxWTy7rgcPHqw5g3nt2rXRqlUr2R5TcnIy9uzZg48//li+yy/L0s5LIDo6WpiYmIh9+/Zplq1evVqUK1dOKBQK2eLK/Yv3xYsXhaWlpfjuu+9kiyGEyPVkBx988IFo1KiRrHGFELKdXO3tx5SVlSW6du0qXF1dJY81btw40bt3b61lI0aMEB06dJA8lhBCRERECCMjI3HgwAFZ5i9OnJ2dxaeffqq5ff36dQFAnDlzRraYcv/i/ejRI1GvXj3Rp08fkZaWJlucx48f51i2YcMGYWxsLJ48eSJbXCGEbCdXi4+Pz3EytUWLFomyZcuKhIQESWNt3rxZ2NnZae0dsGbNGmFqaipSUlIkjSXE69eo2rVrixkzZkg+d3EzceJE0bx5c82eAykpKaJ69erim2++kS2m3L94Z2RkiMGDBwtHR0dZ90R69eqVSE1N1Vp27949AUD89ddfssUVQsh2crW0tDStszkLIcS5c+cEAHHy5ElJY0VGRgpjY2Pxzz//aJZln4z11KlTksbKNmbMGNGsWTNZ5i5Otm/fLszMzLQ+f73//vuiU6dOssUsil+8V6xYISwsLCT/W3zb259b09LSRN26dcW0adNkjSvnydXefkzPnj0TFStWFAsWLJA0TlZWlnB0dBRLly7VLEtPTxcODg7iq6++kjRWNm9vb2FmZibr56lie4x37dq1MXfuXEyaNAkRERFIS0vD0qVL4eXlJfmJcBISEnDixAkAr0+sFRwcjN27d6NmzZqSXpcxNDQUffr0QbNmzeDi4oLdu3dr7uvXrx/Kly8vWaxvvvkGKSkp6NixI8zMzHD27FkcPnwY+/fvlyxGUfP29sbly5fh7u6OsmXLYu/evbh58yaOHTsmeaz//Oc/aNeuHT766CN07twZgYGBOHjwII4fPy55LOD1iRyqVKmC/v37yzJ/cbJq1SoMGjQIpqamcHR0xOrVqzFw4EBZru164MABpKen4+HDh1Cr1di9ezfKlCkj6XUZVSoV3NzckJCQgOHDh+PgwYOa+9q2bSvpL37Hjx/Hli1b0KdPH9jb2+PmzZvYsGEDvvrqK1muZV8UQkJC8Pnnn2PQoEGoXr06Ll++jB07dmDlypWwsbGRNNaHH36o2bNk7NixePz4MVasWIFFixbBwsJC0ljA65PsxcTElPpfxQBg0aJFaN26Nd5//3306dMHe/bsgbm5OWbPni15rEuXLuHhw4cICgpCRkaG5r20d+/esLW1lSzO5MmT4evriyVLluDKlSu4cuUKgNe/Wkl5Iqa4uDiMHDkSgwYNgpOTE2JjY7Fu3Tp0794dgwYNkixOUUpNTYWrqysGDBiABg0aIDIyEr/++itGjRoFd3d3SWM5OTlh3rx5GD58OGbPng1TU1P88ssv6NWrlyzvKwqFAvv27cOPP/4o+dzFzciRI7FmzRq4u7tjypQpCAoKwvHjx3H+/HnJY927d09zEi0AOHnyJO7cuYOWLVtKen6Mbdu2Ye7cuZg4cSJevnypef0oV64cBgwYIFkcAPDw8ED37t3RtGlTJCcnY8eOHVCpVPjiiy8kjVOUJk6ciGrVqqFt27ZITEzEmjVrUL16dcydO1fSOEZGRprXjMTERDg6OmLPnj1Qq9U5TuQpFW9vbwwYMEDWz1NGQggh2+wS2LdvH06cOAFTU1MMGjQIHh4ekseIiYnB/Pnzcyxv06aNpH9IAQEBWLVqVa73/fzzz5L/Rx8+fBhnzpxBWloa6tatizFjxqBatWqSxsjNxIkTMXToUM3JkqR08eJFHD58GEqlEvXr18f48eNlOesg8PqMpevWrUN0dDQcHBwwduxY2U4gN3fuXDRs2LDUnvX4bYGBgfjzzz+hVCrh6uqKCRMmwNRU+u8BJ06cmGO3aEtLS3h7e0sWIzExMc/G6pNPPkH37t0liwUAwcHB2LdvHx4/foyaNWti8ODBaNGihaQxcrNixQoAkPzNFQAiIyPx559/4sGDB6hVqxaGDx8u25lLVSoVNm7ciMDAQNjZ2aFv376y7Wbo7e2N27dvY+XKlbLMX9zExcVh7dq1ePDgAVxcXDB16lTNiZiktGLFCgQGBuZYvmzZMknPwD9v3jzExsbmWN65c2dMnz5dsjjA60Pptm7dirt376JSpUro3Lkz+vXrJ9/ujv/f8ePHsW/fPklfE7MpFAps2bIFISEhqFSpEnr06CF50/0mX19fHD16FKampmjXrh0+/PBDWa5E4u/vj19++QVr1qyR9Iue4iotLQ3r1q1DcHAw7OzsMHHiRFnOtn/8+HFs3bo1x/Jx48ZJ+nly8+bNuR4KaWdnJ/kJtdLS0rB9+3YEBgaibNmyaNasGT744ANZvuh9U3ZfI/VrIgCo1Wr89ddf+Pvvv2FmZoaWLVvigw8+kOUzHPD6M8+2bdugVCrRoEEDTJo0SZa6UygUmDJlCj777DN06NBB8vmzFfvGm4iIiIiIiKgkK7YnVyMiIiIiIiIqDdh4ExEREREREcmIjTcRERERERGRjNh4ExEREREREcmIjTcRERERERGRjNh4ExEREREREcmIjTcRERERERGRjNh4ExEREREREcmIjTcRERERERGRjNh4ExEREREREcmIjTcRERERERGRjNh4ExEREREREcno/wE9SHwt4lujwwAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 1000x270 with 4 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Hessenberg similarity error: 5.009328107237119e-16\n"
          ]
        }
      ],
      "source": [
        "# Finish any remaining columns; this also makes Run All work.\n",
        "while j < n - 2:\n",
        "    H, Z = hessenberg_step(H, Z, j)\n",
        "    j += 1\n",
        "    hess_frames.append((H.copy(), f\"After reflector {j}\"))\n",
        "H0, Z0 = H.copy(), Z.copy()\n",
        "show_fill([hess_frames[i] for i in (0, 1, 3, 6)])\n",
        "print(\"Hessenberg similarity error:\", la.norm(H - Z.T @ A @ Z) / la.norm(A))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "722745b3",
      "metadata": {},
      "source": [
        "## 3. Complex Wilkinson-shifted QR\n",
        "\n",
        "For the active leading $m\\times m$ Hessenberg matrix $X$, choose **Wilkinson's shift** $\\mu$: the eigenvalue of the trailing $2\\times2$ block closest to $x_{mm}$.\n",
        "With complex arithmetic, one QR step suffices even when $\\mu$ is nonreal; no real double shift is needed."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "1f1ae210",
      "metadata": {},
      "outputs": [],
      "source": [
        "def wilkinson_shift(X):\n",
        "    \"\"\"Trailing 2 x 2 eigenvalue closest to the bottom-right entry.\"\"\"\n",
        "    C = X[-2:, -2:]\n",
        "    roots = la.eigvals(C)\n",
        "    return roots[np.argmin(np.abs(roots - C[-1, -1]))]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8c0378e4",
      "metadata": {},
      "source": [
        "Factor $X-\\mu I=QR$ and update\n",
        "\\[\n",
        "X^+=RQ+\\mu I=Q^*XQ.\n",
        "\\]\n",
        "To maintain $A=ZHZ^*$ after a deflation, also transform the upper-right coupling and accumulate $Z$. The explicit dense QR below is for teaching ($O(m^3)$ per step); an implicit Hessenberg QR step costs $O(m^2)$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "3e89a126",
      "metadata": {},
      "outputs": [],
      "source": [
        "def shifted_qr_step(H, Z, m):\n",
        "    mu = wilkinson_shift(H[:m, :m])\n",
        "    Q, R = la.qr(H[:m, :m] - mu * np.eye(m))\n",
        "    H[:m, :m] = np.triu(R @ Q + mu * np.eye(m), -1)  # Clear roundoff below the subdiagonal.\n",
        "    H[:m, m:] = Q.conj().T @ H[:m, m:]\n",
        "    Z[:, :m] = Z[:, :m] @ Q\n",
        "    return mu"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5d3de062",
      "metadata": {},
      "source": [
        "### Stop when the trailing subdiagonal is small\n",
        "\n",
        "For active size $m$, stop iterating on its bottom eigenvalue when\n",
        "\\[\n",
        "|h_{m,m-1}|\\le\\max\\left\\{\\tau\\bigl(|h_{m-1,m-1}|+|h_{mm}|\\bigr),\\;\\epsilon_{\\rm mach}\\|H\\|_\\infty\\right\\},\\qquad\\tau=10^{-12}.\n",
        "\\]\n",
        "Then set that subdiagonal entry to zero, **deflate one eigenvalue** ($m\\leftarrow m-1$), and continue. Stop when $m=1$; there are no remaining subdiagonal entries. A step cap guards against failure to converge."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "5cb0cd7a",
      "metadata": {},
      "outputs": [],
      "source": [
        "tol = 1e-12\n",
        "\n",
        "def deflation_threshold(H, m):\n",
        "    scale = abs(H[m-2, m-2]) + abs(H[m-1, m-1])\n",
        "    return max(tol * scale, np.finfo(float).eps * la.norm(H, np.inf))\n",
        "\n",
        "def advance(H, Z, m):\n",
        "    if m <= 1:\n",
        "        return m, \"complete\"\n",
        "    if abs(H[m-1, m-2]) <= deflation_threshold(H, m):\n",
        "        H[m-1, m-2] = 0.\n",
        "        return m-1, \"deflate 1 x 1\"\n",
        "    mu = shifted_qr_step(H, Z, m)\n",
        "    return m, f\"QR with shift {mu:.4g}\"\n",
        "\n",
        "H, Z, m = H0.astype(complex), Z0.astype(complex), n\n",
        "qr_steps = 0\n",
        "qr_frames = [(H.copy(), \"Hessenberg\")]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2b2f8418",
      "metadata": {},
      "source": [
        "Repeat the next cell with **Ctrl-Enter**. It either takes one Wilkinson-shifted QR step or deflates the bottom eigenvalue when its subdiagonal is small enough. Watch the magnitude and cutoff."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "36ac5cbe",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "QR with shift 2.074+0j; trailing |h|: 2.59e-01 -> 7.10e-03 (cutoff 2.12e-12)\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOUAAAESCAYAAAACKRIDAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJ+xJREFUeJzt3XlYVPX+B/A3iwzbOAgoiwsqiyFqIq65oYASmiSa+6XMpZuamj4/LZe4pmGlqGVp5cXlunVDy11DRU2RVBDlpmaoiAjmAjIgMTAMn98fxjyOs8BhZvAon9fznOdxvuecz/kMzme+Z/2OBRERGGOiYfmsE2CMaeKiZExkuCgZExkuSsZEhouSMZHhomRMZLgoGRMZLkrGRIaLkjGR4aJkTGS4KA1YsWIFDh48WOPlY2NjcfToUTNmJGwb+vJ/8OABNm/ejNjYWGzcuNHEGbIqubm5WLJkCa5duyZoPYtnce9reXk5Tp06hUuXLqG8vBwtWrRA//794eLiorXsmjVrUFBQ8DhZCwvY2dnBx8cHYWFhsLOzM2uerq6uGD58OL755ht12+eff46XX34ZAwcO1Fre2toaM2fOxPLly82Wk5Bt6Mo/OzsbQUFB6NChA7p37w4vLy+Eh4dj8+bNGDNmDFq3bl2jPCorK5GSkoLffvsNKpUKL730EoKDg2FpWfPv+ZKSEhw4cAA3b96Em5sbBg0apPUZOHv2LBITE/XG6NmzJ/r166dznlKpxIoVK6BUKjFu3Di0bNnSYD7nz5/HoUOHMG/ePJPk/uuvv6JHjx7Yu3cvBg8eXG3MKnXeUx48eBC+vr54++23cenSJeTl5WH16tVo1qwZYmJiUFlZqbH8ihUr8NVXX0GhUKC0tBRXr17FxIkT0apVK5w5c8asuc6ePRsREREabR9//DH27t1r1u2aiq78v/32WygUChw6dAixsbF45513kJWVhYULF+KPP/6oUdwDBw7A29sbkZGRSE1Nxf/+9z+MHTsWHTt2xG+//VajGJcvX4a/vz+WLFmCu3fvYsOGDfDx8UFSUpLGciqVCgqFQmvavXt3tTnHxMTggw8+wMKFC2vUW509exYfffSRyXJv1qwZ5s+fD19f32pjaqA6tH//frKysqK33nqLysrKNOZt27aNrKysaMqUKRrt3t7e9PLLL2u05ebmUsOGDcnf39/cKWtxcHCgqVOn6pxnZWVFs2fPNuv2jd3GyJEjqVWrVhptx44dIwB08ODBGsWIiYmhsWPHklwuV7cVFBRQq1atyNfXlyorKw2uX1FRQf7+/tS5c2eNz8GIESPI2dmZCgoKDK5fWVlJ3t7eZG9vr5HDk5KTk8na2ppGjBhBAOjw4cPVvq+1a9eSlZWVWXOviToryrKyMmratCn5+flpFWSVGTNmEAA6ffq0uk1XURIRvfbaawSA7t69a3C7mzdvpo0bN2q0HT58mBYvXkzp6eka7cuWLaMDBw6oX8fFxalf//XXX7R48WKysbGhbt260eLFi2nx4sW0ZcsW9fJVBVNZWUk///wzLVu2jLZt20Z//fWXwRyrqFQqSkxMpNWrV9OaNWs0/g612caT+RcUFNDixYvJ39+fnJ2d1fknJCTQ22+/TQAoOjpa3X7mzBm9ef7xxx862+fOnUsA6PLlywbf58GDBwkAbd26VaP9/PnzBIC++OILg+tXfYmMHz9e5/xHjx6Rt7c3TZo0idauXWvSohSS++3bt2nx4sWUmZlZ7bafVGe7r4mJicjNzcWECRNgY2Ojc5l3330XALBhw4Zq49Hfh8IWFhYGl7tw4QKmTZsGpVKpbluyZAkWLlyocax17do1/N///R/y8vLUbbGxsdi9e7d6ewqFAkSksUv1ZFzg8bHWuHHjEB8fjxs3bmD69Ono3Lkz/vrrL4N5yuVydO7cGRMmTMDly5dx5coVzJkzB8HBwVrL1nQbuvJXqVTqf1flX15eDuDxMVhVu0ql0purvt2x/Px8ANX/n5w4cQIA0KNHD432jh07wsHBQT1fn/Xr1wMAJk+erHP+7NmzUVpaimXLlhmMUxtCcs/JycHChQvx+++/C9uIoBI2wr/+9S8CQEeOHDG4nKOjIwUGBqpf6+op79y5QzKZjNq2bVvtdg8dOkQA6Pjx40T0+FvUxsaGgoKCqHXr1url1qxZQwDo5s2b6jYXFxd65513NOJVt/vq5uam8a1c9Q26cuVKg3l++eWXZGlpqdXz//LLL7Xehq78IyMjydvbW6NN6O6rLtnZ2eTg4EB+fn7V7r6OGjWKAFB5ebnWPF9fX43//6fJ5XKys7Oj9u3b65x/4MABAkC7du0iItLbU1ZUVKj3CqqmIUOGkKWlpVb7vn37apV7SkoKAaC9e/ca/Hs8zbp23xfCFRYWAgCcnJwMLieTydTLVrl79y6WLFkCALh37x4SEhLg6uqKrVu3VrvdPn36QCKRIDExEX379sWJEyegUqkQGxuLgQMH4vr16/D29sbhw4fh5+cHLy+v2rw9tTZt2iA0NFT9OjAwEG3atMGpU6cwc+ZMves9evQIRIR79+6hSZMm6vbevXubbBvmolAoMGLECJSWlmLNmjXV9pR//fUXLCws0KBBA615tra2Bvcqtm/fjtLSUp29ZH5+PiZMmIARI0YgMjLSYA5KpRILFy7UOe/p9gkTJmDQoEFG515TdVaUMpkMAFBUVGRwuaKiIrRt21ajjZ7YdczPz0dhYSFCQ0Ph7+9f7Xbt7OzwyiuvIDExEZ988gkOHz6MLl26IDQ0FM7Ozjh8+DBatmyJY8eOYcyYMbV/g397OncA8PDwQG5ursH1oqOj8c0336BTp07o168fgoODERoaii5duphsG+agVCrxxhtv4MyZM/j6668REhJS7Tp2dnYgIiiVSq0Pt0KhgKOjo9514+PjYWdnh3HjxmnNe/fdd1FWVobVq1dXm4O1tTXmz5+v0abvkkjnzp1NkntN1VlRBgYGAgDS0tL0XlfKzMxEcXGx1ofO3d1d3VMCwPHjxxESEgJHR0esXbu22m2HhYVhwYIFyM/Px+HDhxEVFQVLS0v0798fiYmJ6NixIwoLCxEWFmbEO3zMwcFBq83Kykrr2PNpTZs2xZUrV7B7924kJSVh48aNmDdvHvr27Yv9+/drxK3tNkxNpVJh7Nix2LdvH1auXIkpU6bUaL2qvZG8vDyNPRMiQl5ent7/h99++w3nzp3Dm2++qXOPKyEhAX369MF3332nbjt37hwAYOvWrfj111/x4YcfwsrKCtbW1hqfKQD45ptvkJiYqNVuityFqLMTPQMHDoSbmxvWr1+PiooKnctU/THHjx9vMFZwcDAmTZqEb7/9Vv1HNyQsLAyVlZXYvHkzLl26pP7DhYWFISkpCQcPHoSVlZXeL4snCbk4LpS9vT1Gjx6NdevW4erVq9i0aRNOnDiBLVu2mG2bQO3eU2VlJd58800kJCQgLi5O0G5znz59AAApKSka7RkZGSgpKVHPf1p8fDwAYNKkSTrnz58/H71799a4nln1RVV1EouMvFemtrkLIugI1Eg//vgjWVhY0DvvvENKpVJj3k8//UTW1tYUHR2t0a7vkkhubi7Z2dlRaGhotdtVqVTk4uJCLi4uJJVK1QfpWVlZBIBcXFyoR48eWuvpOlHi4+NDb7zxhs7t6LuGGBISQkFBQQZzPHv2LBUXF2u0VZ3AWbNmTa22UdMTPb///jsBoE2bNhnMsUplZSVNnDiRANDy5curXfbpS0fl5eXk4+ND3bt31/gcjB07lmQyGd2/f18rTllZGbm6ulJAQECNcqxi6ksiQnIX/YkeABg6dCgSEhIwdepUJCUlYfDgwbC3t8fZs2eRlJSEd999F3FxcTWK5enpiSlTpiAuLg5JSUno37+/3mWrdlUTEhLw2muvqY8FWrZsCR8fH1y7dq3Gux2jRo3C8uXLMW3aNLi7u6NVq1YYO3ZsjdY15MKFCxg+fDi6desGHx8fFBcXY8eOHejZsyf+8Y9/GB3fEF9fX3Tq1AkxMTG4evUq7OzsMGDAAHTt2lXn8itWrMC///1vvPTSSygtLdXa3Xvydj2VSoWFCxciJCRE/Xdq0KABfvjhB7z66qvo0aMHwsPDce7cOSQnJ+P777+Hq6ur1jb37NmDBw8eYMGCBSZ+98LUJneh6rQoAWDYsGGIiIjA0aNHcfnyZRQWFiIlJQXt2rXDp59+qnUNc+rUqXqva3744YewtbXF/fv3q93u9OnT4efnp3HWEnh8Le/ixYs6T/LMnj0bAQEBGm0ff/wxunfvjgsXLmhdp5w/f77W9SsAGDduHEpLSw3mN2nSJIwcORKJiYnIzMyEh4cHEhIS0KtXL43lhGxDV/4jRoxA3759NdosLS2RlJSEn376CdnZ2dVep2zTpo36JIlCodCa/+S6lpaWmD9/Pnx8fDSWCQwMxNWrV7F7925kZ2dj6NCh2LhxI9zd3XVu08rKCvPnz0d0dLTevHQJCgrC/Pnz0apVqxovWx2huQv1TG5If9qRI0cwaNAg9OnTB/v27YNEInnWKTFmtFOnTqF3795ITEwUdAJIFI9uhYaGYvfu3ejWrZvZbzJnrK7cunULANCiRQtB64mip2TsRfLw4UMsX74cmzdvRps2bXD48GFB64uip2TsRUJEkEgk+Oyzz3DgwAHB63NPyZjIcE/JmMhwUYrEnTt3cPv27Wedhk5izu1FxLuvNVBWVoacnBzY2dmhadOmtY6Tm5sLIkKzZs205oWHh+PBgwdITU01JlWzeNa5FRQUoKKiQuPpmRcZ95Q18N1338HX1xdt27Y16tGcsWPHYvjw4TrneXp6onnz5rWObU7PIjelUomYmBg0btwYrVu3RkBAABwdHTFt2jTI5fI6zaXOCbopr57q0KEDtWvXjgDQhg0bah2nb9++1K1bN9Ml9gKbNWsWAaD4+Hh1W2JiIkkkEhoyZMgzzMz8ePe1GmfPnkW3bt2wZ88erFy5EmVlZUhOTja4Tn5+PlQqlcbuVnZ2NkaMGIHy8nIkJCQAeHzrWNXtX3fu3IFKpVLv2t6+fRtEpLOHksvluH//Plq0aKF1C+L9+/dRXl4ODw8PQU9/yOVyKBQKuLm5ac17OreioiLcu3dPZxwbGxudF8uF5uXn5wcLCwtcvXpVoz0yMhL79++HUqms9mHq59Yz/lIQvYkTJ5KXlxepVCpKSEggAHTp0iWdy3733XfUunVrsra2piZNmlDTpk3pm2++IaLHvaStrS1JJBLy9vYmb29vjaEjBg4cqPGUx5QpU0gikVB+fr7WdsaNG0dOTk5UWlqqbvv222/Jy8uLbG1tydnZmVxcXCguLq7a95eQkECtW7cmR0dH8vT0pEaNGtG8efOopKREb24//PCD+j08OVlaWmo9gVLbvPr370/u7u5UUVGh0d6zZ09q0aJFtes/z7goDSguLiapVEqxsbFERKRUKsnT05Pef/99rWU//vhjsrCwoGXLlqlH68vLy6OZM2eqlzG0+/r0B7/qsa0vv/xSY7nCwkKys7PTGCcoNjaWLCws6Ouvv1Z/iPfs2UMSicRgAWRlZZG1tTV99NFH6vXkcjktXbpUY0S6p3PTpWqMo0mTJhmdFxFRWloaubm50ejRoyklJYXS0tJo9uzZ5ODgQPv37ze47vOOi9KAdevWkY2NjcZgVjExMeTi4kIKhULddv/+fZJIJDRu3DiD8YQUJRFRp06dtJ4lrfrwnz9/noiI8vPzydbWVudwi9OmTSNnZ2et3qbKvn37CAD9+uuvBvOurij37dtHVlZW1L9/f/WzqsbkVWXTpk0kk8nI0dGRnJ2dqUGDBvTJJ59Uu97zjs++GrBu3Tr069cPRUVFuHbtGq5du4b+/fujoKAAP/30k3q55ORklJWVYciQISbd/oQJE3Dx4kWcP39e3bZ+/XoEBgaqh1dJTk6GQqFAYGAgbt68iaysLNy4cQM3btxAy5YtUVBQgMzMTJ3xg4KCIJVKMWbMGKxcuRIZGRmCn8xPT0/HqFGj0KZNG+zcuVP9rKoxeQHAvHnzMH78eKxevRrFxcXqoVw+++wzjBo1SlCOz51n/a0gVhkZGQRA57FT48aNKSQkRL3spk2bCAAlJSUZjCm0p6zaVa0aNb4qp6+++kq9zMaNGwkANWvWTGeu3t7edOHCBb05Xbp0iaKjo8nd3Z0AUOPGjWnu3LkaewL6espbt26Rp6cnubm5UVZWlsY8Y/IqLi4ma2trGjp0qNa8hQsXEgCD7+l5V+cPOT8v1q1bhw4dOuDixYta806ePIm+ffuqh6esuqHg5s2bJs1BJpNh2LBh2LZtG+Li4hAfHw9bW1uNB7Krtr18+XKMHDlS8Dbatm2LTZs2AQCuX7+O+Ph4LF26FA4ODnqHYAQen4EdNGgQCgsLcfz4ca0fzzEmrzt37qCiokLncJ9VbTk5OXj55ZcFxX1uPOtvBTEqLS2lRo0a0Zw5c3TOVyqVJJPJ6MMPPyQiIoVCQe7u7tStWzetgYifHMdlyJAh1KFDB50x9fVGx48fV18fdXFxoTFjxmjl6uHhQWFhYToHQdY1aHB18xo1aqQxDtHTuSmVShowYABZWlrSTz/9pDOGMXmVlpaSg4ODzr2Kt956iwDQtWvX9K7/vOOi1GHLli3V7o4OHz6cPDw81EV36NAhsrW1pbCwMDpw4AClp6fThg0bND7MS5YsoQYNGtCePXvojz/+oBs3bqjnGTqZ4uvrS05OTgSAjh49qjX/yJEj5ODgQBEREbR//37KyMiggwcP0rx586h79+5638P69espLCyMNm/eTKmpqZSWlkYffPABWVhY0I4dO/TmtmjRIgJA06dPp8zMTI0pOzvb6LyIHv+uC/7+fZPk5GQ6d+6cOrfJkycbXPd5x0Wpw/jx4ykgIEDvDxEREW3fvp28vb3pxIkT6rbLly/T5MmTKTAwkDp16kQTJkzQ+DGckpISmjlzJnXq1Il8fHw0rlOOHz+eXn/9dZ3b+vLLL8nb25t69uyp9ycBrl+/Tu+//z5169aNAgICKCIigj799NNqfwXq6NGjFB0dTUFBQdS5c2caO3as1k8lPJ3bwoUL9R4n9u/f3yR5ERGdPHmSxowZQ127dqXAwECKiorS2zO/SPiOHsZEhi+JMCYyXJSMiQwXJWMiw0XJmMhwUTImMlyUjIlMnd9mV1lZiby8PEil0hf3IVXG8Hj81+LiYnh6egp64LzOizIvL0+0Y9EwZg45OTk6B0vTp853X6VSaV1vkrFnSuhnvs6LkndZWX0j9DPPJ3oYExkuSsZEhouSMZHhomRMZAQXZUpKCqKiohAUFIRRo0bht99+M0dejNVbgooyPT0d/fr1g4+PD1atWgVHR0f06tXL5GPTMFavCXkiOioqioKDg9WvKysryc/PT2Ng4OrI5XICwBNP9WaSy+VCykzYuK/Hjh1DRESE+rWFhQUiIiJw7NgxIWEYYwbU+Da7kpISPHz4EB4eHhrtHh4eyMnJ0bteWVkZysrK1K+LiopqkSZj9UeNe8qKigoA0PqVJ4lEAqVSqXe9pUuXQiaTqSe+75Uxw2pclFKpFDY2NsjPz9doz8/Ph6urq971PvzwQ8jlcvVkqFdljAkoSktLS3Ts2BFnzpzRaD99+jSCgoL0rieRSNCwYUONiTFmgJCzQvHx8eTg4ECpqalERPTzzz+TlZWVoJ8m47OvPNW3SejZV8GDMc+ZM4ckEgl5enqSnZ0dff7554LW56Lkqb5NQouyVoMxFxcX488//0TTpk1hb28vaN2ioiLIZDKhm2TsuSWXywUdttVq5AGpVMoPKzNmJnxDOmMiw0XJmMhwUTImMlyUjIkMFyVjIsNFyZjIcFEyJjJclIyJDBclYyLDRcmYyHBRMiYytSrKX375BatWrcLvv/9u6nwYq/cE3ZB+4sQJTJ06FU5OTkhOToarqyteeuklc+XGWL0kqCglEgm+//57tGvXjn89izEzEVSU3bt3N1cejLG/mf2XnHmIScaEMfvZVx5ikjFhzF6UPMQkY8KYffdVIpFAIpGYezOMvTD45gHGREZQT5mTk4OdO3eqX//888948OAB2rVrh9DQUJMnx1h9JKgoS0tL1b9FOWPGDADAzZs3tX70hzFWe7Ua99UYPO4rq2+EjvvKx5SMiQwXJWMiw0XJmMhwUTImMlyUjImM2e/o0UfoGama4kfK2POOe0rGRIaLkjGR4aJkTGS4KBkTGS5KxkSGi5IxkRFUlKdPn0ZkZCQaNWoER0dHBAcH48yZM+bKjbF6SVBRxsXFYeLEicjKykJeXh4CAgIwYMAAHuKDMVMiI5SVlZG1tTWtX7++xuvI5XICQHK53JhN6wWAJ55ENQn9rBt1THn//n2oVCo0atTImDCMsSfU+jY7IsJ7772Hli1bYuDAgXqX43FfGROotruJs2bNIplMRunp6QaXi4mJMUmXXlO6tsUTT89yEvpZr1VRzp07lxo2bEhnzpypdlmFQkFyuVw95eTk1CrRmnrW/wE88fT0JPSzLnj3dd68eVi7di1+/vlndO3atdrledxXxoQRVJQfffQRVq9ejX379qFjx45QKBSPg1hbw9r6mT0FxtiLRUi3KpPJSCKRaE2LFy+ucQy+JMJTfZvMuvtaWFgoZHHGWC3wva+MiQwXJWMiw0XJmMhwUTImMlyUjInMC3dxkcz8e0U8hCUzN+4pGRMZLkrGRIaLkjGR4aJkTGS4KBkTGS5KxkRG8CWRs2fPYs+ePXjw4AH8/PwQHR0NV1dXc+TGWL0kqKdctGgR5s6dC6lUipdffhn79+9H27ZtkZWVZa78GKt/hDzndfv2bY3XSqWSmjRpQrGxsTWOYe7nKc0NIng+j6fnazLrEJNNmzbVeF1QUICSkhI0b95cSBjGmAGCjylv376NDz74ACUlJTh37hwWLFiAMWPG6F2eh5hkTBjBZ1+lUinCw8MREhICX19fbNq0yeDPFixduhQymUw9ca/KmGEWfx8n1YpKpULnzp3RsWNHbNiwQecyunrK5s2bQy6Xo2HDhrXd9DPDN6QzoYR+1o16SsTKygr+/v64efOm3mV4iEnGhBG0+5qQkKDxaNStW7dw5MgRvPLKKyZPjLH6SlBR/vLLL/Dz80NkZCQGDBgAf39/9O/fH/PmzTNXfozVO4KPKe/cuYO0tDQ0aNAA/v7+aNGihaANFhUVQSaT8TElqzeEftaNOtFTG1yUrL4R+lnnG9IZExkuSsZEhouSMZHhomRMZLgoGROZF27cV3Mz58lqPrPLAO4pGRMdLkrGRIaLkjGR4aJkTGS4KBkTmVoXpUKhwKlTp/D777+bMh/G6r1aF+X06dPRt29fLFiwwJT5MFbv1aooExISkJqaipCQEFPnw1i9J7gos7KyMGPGDGzduhU2NjbmyImxek3QHT0VFRUYPXo0FixYAH9//xqtw0NMMiaMoJ5y/vz5cHV1xZQpU2q8Dg8xyZgwNR55IDU1Fb169cK2bdvQpEkTAMDcuXNhY2ODxYsXo0uXLjpHrXvRhpg0J7739cVktiEmS0tL0blzZ6xYsULdduXKFVhaWuKDDz7Ajh074O7urrUeDzHJmDBGjdEzePBg2NraYseOHTVe53kfo8ecuKd8MfEYPYw954x6njIgIIAvizBmYjzEpIjw7uuLiXdfGXvOcVEyJjJclIyJDBclYyLDRcmYyPAQkyLCw1cygHtKxkSHi5IxkeGiZExkuCgZExkuSsZEhouSMZERdEnku+++w5dffqnR1rBhQ5w+fdqkSTFWnwkqynv37sHKygpbt25Vt1lZWZk8KcbqM8E3D9jZ2aFdu3bmyIUxhloUZWZmJnr06AFbW1t07doVc+fOhbOzszlyY6xeElSUEokE//znPxEeHo7CwkJ88skn+O9//4uMjAy9D3HyuK+MCUQCKJVKjdcPHz6kRo0a0dKlS/WuExMTQwC0JrlcLmTTzEi6/g94qptJ6Gfd6OFAgoOD0axZM2zZskXnfB73VRz4hvRnx2zjvuqTm5uLtm3b6p3P474yJoygmwc++ugjPHjwAACgUqmwaNEi3LhxA2PHjjVLcozVR4J6yhYtWiAwMBCWlpZ4+PAhmjRpgp07d6Jnz57myo+xeqdWx5S3b9+Gvb19rS6F8BCTzwYfUz47dXJM2axZs9qsxhirAb4hnTGR4aJkTGS4KBkTGS5KxkSGi5IxkeFxX+sJI++mrBZfcjEd7ikZExkuSsZEhouSMZHhomRMZLgoGROZWp19zcvLQ1JSEiwsLDBw4EC4urqaOi/G6i3BPeVXX30FX19f7NixA8ePH0dwcDDS0tLMkRtj9ZKgnvLYsWOYPn06du3ahSFDhgAAHj58iLt375olOcbqI0HPUw4ePBiPHj3C8ePHa71Bfp7yxcQ3D+hn1ucpU1JSMH36dFy+fBknT55EkyZNEBwcjEaNGuldh4eYZEyYGh9TEhEKCgqQlJSEqKgopKamYvny5fD19UVycrLe9ZYuXQqZTKaemjdvbpLEGXtRCdp9tbe3h5ubGy5dugR7e3sAwIgRI/DHH3/gwoULOtfhISbrB9591c+su69t2rSBv7+/uiABIDQ0FLt27QIR6fyP4SEmGRNG0CWRqKgoXLhwASqVSt2WmpoKHx8f/qZkzEQE7b4+evQIffr0gaOjI8LDw3Hp0iXs2rULP/74IwYOHFijGHz29cXEX8r6Cf2sC+opHR0dcfr0aYwfPx5FRUXo2rUrrly5UuOCZIxVz+jfEhGKe8oXE/eU+pm1p2SMmR8XJWMiw0XJmMhwUTImMlyUjIkMDzHJTMKcJ/Hr25ld7ikZExkuSsZEhouSMZHhomRMZLgoGRMZQWdfy8vLUVlZqdVuZWWFBg0amCwpxuozQT1lly5d4OTkpDHZ2dlh1KhR5sqPsXpHUFFevHgRCoVCPaWnpwN4PCQIY8w0jDqmjI+Ph6urK4YOHWqqfBir92pdlEqlEps3b0Z0dDRsbGxMmRNj9Vqtb7Pbu3cv7t27h4kTJxpcjsd9ZUyYWveU8fHx6NWrF/z9/Q0ux+O+MiZMrYYDyc3NhZeXF9avX4/o6GiDy/K4r8xYz/sN6WYd97XKxo0bIZVK8cYbb1S7LI/7ypgwgndfiQjr16/HuHHjYGdnZ46cGKvXBBflyZMnkZubi0mTJpkjH8bqPcG7r3369IFCoTBHLowx8A3pjIkOFyVjIsNFyZjIcFEyJjJ1Pppd1b0KfLsdqy+E3p9T50VZXFwMAHy7Has3iouLIZPJarx8nf/qVmVlJfLy8iCVSqu9farqlrycnByz3JJnzvice93HNnd8obGJCMXFxfD09ISlZc2PFOu8p7S0tESzZs0ErdOwYUOz3idrzvice93HNnd8IbGF9JBV+EQPYyLDRcmYyIi6KCUSCWJiYsz2lIk543PudR/b3PHNnXuVOj/RwxgzTNQ9JWP1ERclYyLDRcmYyIj2R2OJCGfPnsWdO3fQtm1b+Pn5mTS+XC7HkSNH4Obmhl69epk0dm5uLi5evAiZTIbAwEDY29ubLLZKpUJGRgZycnLQsmVLdOjQwWSxn5ScnIycnByEh4fDycnJ6Hjp6em4evWqRptUKsWgQYOMjl2lsrISaWlp+PPPP9GlSxe4u7sbHbOkpAR79+7VOa9Tp04m/1wCAEiEioqKqHfv3uTh4UFhYWHk4OBAM2fONFnsSZMmkYeHB7m5uVFkZKRJ4hIR5efnU1RUFDVv3pwiIiKoffv25O7uTkePHjVJ/JMnT1K7du2oS5cu9Nprr1Hjxo2pe/fu9ODBA5PEr3L+/HmSSqUEgNLT000Sc8aMGeTp6UkjR45UTzNmzDBJbCKizMxMat++PXl5edHrr79Obdq0ofj4eKPj3rt3TyPnkSNHUu/evQkA7dy50wSZaxNlUc6ePZtatWpF+fn5RER05swZsrS0pH379hkd++7du/Ttt9/So0ePaNiwYSYtylu3btHOnTupsrJS3TZ16lRydXUllUpldPykpCTKzs5Wv5bL5dSsWTOaM2eO0bGrFBcXU5s2beizzz4zeVEOGjTIJLGeplAoyM/Pj6Kioqi8vJyIiMrKyujgwYNm2d6MGTOocePG6m2ZmiiPKbds2YK33noLzs7OAICuXbuiV69e2LJli9GxmzRpgsmTJ8PBwcHoWE9r3rw5oqKiNO7pjYqKwoMHD3D79m2j4/fr1w8tWrRQv27YsCE8PDzUN/mbwpQpU/Dqq69iwIABJotZRS6XY8+ePTh+/DgePnxosrgJCQm4du0avvjiC/Wvv9nY2CA8PNxk26hSXl6u/nya65fmRHdMef/+fdy9exft2rXTaG/fvj2OHz/+bJIywpEjR+Dk5ISmTZuaJF5FRQV27NgBhUKBpKQklJaWYs6cOSaJvXnzZqSnpyM1NRVXrlwxScwnXb58GWvXrsXdu3eRmZmJuLg4TJ482ei4v/zyC9q2bQtXV1ckJibCwsICHTt2ROPGjU2QtaZdu3YhPz+/2l8GMIboilIulwOAupes4uLigsLCwmeQUe2dPHkScXFxWLVqFaysrEwSU6VSYdeuXXj06BHOnTuHQYMGoVGjRkbHzczMxKxZs5CUlGSWO1aGDh2K2NhY9UmvtWvXYsqUKejUqRM6d+5sVOx79+7B0tISXbp0gYeHBxQKBdLS0rBy5UqTFP2T4uPj0bdvX/Oc4Klilp1iI9y6dYsA0N69ezXaZ82aRd7e3ibdlqmPKZ+UmppKTk5ONGvWLLPEJ3p8TBkQEEDR0dFGxxo4cCBFRETQ9u3bafv27fTpp58SAPr0008pJSXFBNlqa9q0KS1atMjoOCNHjiQAtH//fnXb119/TQ0aNKCcnByj41fJzs4mS0tL2rJli8li6iK6Y0pPT0/Y2tri1q1bGu3Z2dlo3br1M8pKmPPnzyMsLAxvvfUW4uLizLadhg0bIjIyEidPnjQ6VlBQEKRSKXbt2oVdu3bh2LFjAIBjx44hLS3N6Pi6SKVS3L9/3+g43t7esLOzQ0REhLpt2LBhUCqVyMjIMDp+lQ0bNkAmk2HYsGEmi6mTWUu+liIjIyk4OFj9+uHDhySVSmnVqlUm3Y45esr09HRydnY26en+Knl5eVpt4eHh1KdPH5NvKz093aRnX5/OPSMjg6ysrGjTpk1Gx05JSSEAlJWVpW47duwYAaCMjAyj4xMRVVZWkpeXF02fPt0k8QwR3TElAHzyySfo0aMHRo8ejd69e2Pjxo3w8vIy2fHBjz/+iPLycty+fRsqlQrff/89bGxsEBUVZVTcW7duITQ0FJ6enujevTu+//579bywsDC4uLgYFX/KlClwcnJCt27dQETYt28fTp06hUOHDhkVty6EhIRgwIABaNeuHfLy8rB69WqEhIRg9OjRRsfu3r073nzzTbz66qt47733oFAosGLFCowZMwbt27c3QfaPT9hlZ2fXyS8DiPYpkevXr2PdunW4c+cOAgIC8O6770IqlZok9oQJE1BSUqLR5uDggPj4eKPiXr58GR9//LHOeUuWLIGPj49R8SsrK7Fz506cOHECKpUKfn5+iI6ONrrYdcnOzsbcuXPx2WefwcvLy+h4paWl+M9//oO0tDTIZDL07NkTr7/+uvGJ/o2IsG3bNiQlJcHe3h59+vTB8OHDTfaLXfHx8bh06RJWrFhhkniGiLYoGauvRHeih7H6jouSMZHhomRMZLgoGRMZLkrGRIaLkjGR4aJkTGS4KBkTGS5KxkSGi5IxkeGiZExkuCgZE5n/ByrIqdIAaOy+AAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 250x270 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "if m > 1:\n",
        "    old_m = m\n",
        "    before = abs(H[m-1, m-2])\n",
        "    cutoff = deflation_threshold(H, m)\n",
        "    m, action = advance(H, Z, m)\n",
        "    qr_steps += (m == old_m)\n",
        "    qr_frames.append((H.copy(), f\"{action}\\nActive size {m}\"))\n",
        "    print(f\"{action}; trailing |h|: {before:.2e} -> \"\n",
        "          f\"{abs(H[old_m-1, old_m-2]):.2e} (cutoff {cutoff:.2e})\")\n",
        "show_fill([qr_frames[-1]])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "eb101455",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAABdIAAAEsCAYAAADHIUHdAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAdFpJREFUeJzt3Xd4U/X7//FXSksLpQ1lr0LZyka2iOypTAFBBRVUhopbFAc4UVTciDLEgaJ+QBQcUASRvQUEZe+y6aBA9/v3h7/mS2gb2iZp2uPzcV3nuujJybnvhNy5k3fOeR+bMcYIAAAAAAAAAABkys/XCQAAAAAAAAAAkJ8xkA4AAAAAAAAAgAsMpAMAAAAAAAAA4AID6QAAAAAAAAAAuMBAOgAAAAAAAAAALjCQDgAAAAAAAACACwykAwAAAAAAAADgAgPpAAAAAAAAAAC4wEA6AAAAAAAAAAAu+Ps6AQAoyNLS0hQXF6fg4GAFBATkej8pKSmKj493WhcSEqJChQq5myKAXLp48aJSU1MVEhLi1n7i4+OVkpLi+DsgIEDBwcHupgdYUnpfLVq0qAoXLnzV7VNTU3X+/Hm3+7CnYlwtf2OM4uLiVLhwYRUpUsQr+SJ/SkpK0sWLF2W322Wz2XydTr5w8eJFpaSkKDQ0NFvbX7hwQWlpaW73ZU/FuFr+ly5dUmJiIv/nFhQXF6eAgADex/MQ37uRX3BEOgC4YefOnQoLC9M333yTq/vv2bNHbdq0UWhoqCIiIhQREaGKFSsqLCxMGzZsyPH+kpOTFRMTo9TU1Fzlkx1paWmKiYnRxYsXvRYjP8bGf8/QoUPVuHHjXN//tddeU+nSpVW+fHlHfZcsWVK33HJLrvZ3/vx5XbhwIdf5ZMfFixcVExOjtLQ0r8bJb7H/a5KSkpScnJzl7cYYxcTEOJbY2FinH4O86fDhwwoLC9Pnn3/uWJf+3p+UlJRh++3btyssLExz5871Wk45iZFZ/pJ0+vRp9e7dW8WKFVOVKlU0fvx4JSYmKiYmRsaYHOeUlpbm9v/JpUuXsrWdq9eKO9tenseVgxqupKSkKCYmRomJiTmOJf37w0h273u1WKmpqU61cvly5ev1k08+UVhYmE6ePJmrvHPrap+ZLl265JS3t/vM5UaPHq06deo4rbt48aLi4uIy3f72229Xs2bNvJpTTmJklr8kffnll6pcubJKlSqliIgIRUdHKzY2Nts1l5n0/6fc7uPSpUu5eq9JFx8fr9jY2BzFy06+ua3j3EhKSsr2Z4yrPV/VqlXTI4884qnUkA18785bfO/OGgPpAOBD999/v5KSknTmzBnHF5jZs2fnen8//fRTrj8MXM3KlSs1ZswYVa5cWWFhYRo9erTHY+TH2EBurF69Wk8//bQ++ugjnT9/3lHf9evXz/U+27Rpk+tBeFfOnj2rjz/+WB07dpTdbldYWJgOHz7s8Tj5LfZ/zYkTJzRmzBiFh4crJCREwcHBqlatmsaOHavo6GinbWNjYxUWFqayZcsqIiJClStXVmBgoOrWrauPPvrIq3n6+fnJbrcrMDDQsW737t0KCwvTV1995dXYnpBZ/pL0/PPPa8OGDTpw4IBiYmI0adIkvf/++woLC9PZs2ezte/4+Hh98MEHatq0qex2u4oWLaqqVatq/Pjx2R6Ijo6O1vDhwxUSEqKwsDBVqFBBEydOzDBgtHfvXo0dO1Y1a9ZUaGioihUrplatWumnn37KsM+cbJuZBx98UJ06dXK5TWpqqn799VcNGzZMZcuWVVhYmN54441s7T/dypUr1aZNGwUGBqpYsWKqWbOmpk+f7lasLVu2KCwsTBUrVnQMzKQvVw72BAYGym63y8/P+1/B//e//6lt27YKCgpSqVKlVLRoUXXp0kWLFy/OsO3TTz+tsLAwValSRRERESpRooTsdrsGDBigAwcOeDXPokWLym63O60bNmyYGjRo4NW4npJZ/keOHNHw4cM1YsQIXbhwQTExMSpRooQqVqyoJ554Ildx4uLidO211yosLExPP/10tu+XkJCgJ598UmXKlFFoaKiKFi2qAQMG6OjRo9m6/9GjR/XWW2+pRYsWstvtKlmypEfy3bBhg4YOHaoKFSrIbrcrNDRUXbt21R9//JHtx5YTP/zwg+rVq6eQkBAFBQWpc+fO+vvvvzNsl5PnK/09GAUH37vzf+yCgoF0APChVatWqX379vn+g1hqaqqeeuop1ahRQ8uWLfvPxAZya+XKlZKkm2++2ceZXN0XX3yhzZs3a9y4cXrggQf+M7H/S7Zt26ZGjRpp5cqV+vrrr5WQkKBLly7pk08+0ffff6+mTZvq4MGDGe43ZMgQxxHpZ8+eVatWrTR69GhNmTLFa7lWrlxZMTExGjJkiNdieFNW+a9evVrNmzdXmTJlcr3vpUuXKioqSp9++qmio6MVHx+vV199Va+//rr69+9/1funpaXp5ptv1vLly7Vx40YlJCRo+vTpeuGFF/T88887bfvwww9r3rx5evvttxUfH6+oqCi1atVKN998c4bB55xsm1sHDhzQO++8oxtuuEH/+9//cnz/devWqWPHjipTpoxOnDih+Ph4Pf7447rvvvv04Ycfuh1rxowZGY5Iv/I1cO+99yomJsat10B2jB49WoMGDVLnzp117NgxXbx4UQcOHFCjRo3UtWtXvfjii5neb8eOHY4jiBcsWKDff/9d7du39+rRiFOmTNGOHTu8tn9vyyz/jRs3KikpSTfddJPH4owZMyZXR5Pfeuutmjp1qmbOnKmEhATt2LFDhw4dUtu2bXX+/Pmr3v/dd99VVFSU3n//ffXs2dNj+U6ZMkV9+vTRtm3blJCQoK1bt8rPz08dOnTQqlWrsh0nO37++Wf169dPAwYMUHx8vE6ePKkiRYqobdu2On78uNO2OXm+9u3bp8mTJ3s0V3gX37vzd+yChIF0AMihxMTEHH2YvXjxYobtExISFBUV5fhykv6lKyuXnzocGxub6WmQiYmJjtNx4+Pjszy12FVeWSlUqJDj1+ny5ctn6z7p0qcIyOxU4aSkJJc5uhsbyAljTI5OMU5JScl0+5iYGMfRSwkJCYqJiVFCQkKW+0k/Yi0mJkbnz5/P9LTj2NhYpaamOqYZiImJyfLU96zyysrDDz/sOCo8p/NDpp+6ndlprXFxcVnm6InYyJ6kpCT169dPQUFBWrZsmW644QbZbDYVKlRInTp10u+//66YmBjdeuutLvdTvHhxffjhhypatKg++eSTq8bN7HWf/lq//DWe3iPSX7NXTuOSkpLiGMBIn/7najXl6jZXkpKSst0Xs4pxZf7pp36fOXNGhQoVcpo+I30fcXFxjvWuThHv1auXXn31VdWvX1/+/v4qXLiwBg8erEGDBmnRokU6deqUy5znzp2r1atX6+2331bt2rUlST169ND999+v119/XSdOnHBs27lzZ23evFk333yzChUqpNDQUE2ePFl169bVyy+/7LTfnGybWzVq1HAcJR4WFpbj+7/wwgsKCgrSp59+qlKlSikwMFAjRoxQ//79NW7cOKfPKO7Gykr6Zx53pte4munTp+ujjz7SW2+9pWeffdZxBHHZsmU1adIkjR07VuPHj9fPP/+c5T78/Px04403asyYMTp06FCmR7FfLr3PXS79c+uVg/CXLl1ymh7kymlc4uPjlZyc7Kijq30+Tk5OzvW0Cq4+e2Y3xpX5nz9/3tH/02s//X0vPWZOp9D5/vvvNXv27Aw/+FzN+vXr9eOPP+rRRx911Ga1atU0a9Ys7d+/X++8885V9/HGG2/orbfeUvPmzbMdNzv5fvrpp+rXr59KlSolSapatao++ugjpaamujxCODExMcdTWj3xxBNq1KiRxo8fr4CAAIWFhWnmzJmKi4vTxIkTHdvl9PmKi4tza6oeZA/fu/nenR8xkA4A2bR27Vo1a9ZMQUFBstvtevDBB7OcCzQqKkp33nmnihcvruLFiyskJERDhw51nL79zjvvOOZUfOeddxynAWf1QX3nzp2ObapUqaIiRYqocuXKmjBhgmNAYubMmbrvvvskSX379nVsf/nRVFfLyxtsNptGjBihMmXK6K+//nKsT05OVvv27XXNNdd4NT5wNcnJyXrsscccNVGvXj3HEeWZ+fLLL9WwYUPH9AB169Z1qrOIiAjHQGN6Hb711ltZ7u++++5zbFe2bFkVLVpUnTp10p9//unYpkGDBtq5c6eWL1/u2PaGG27IUV7e8Pfff6tcuXK6++67ndYvXLhQxYsXd/m4kTfmzJmjffv2ady4cRmmIJCkChUqaMyYMVq/fr2WLFnicl+BgYEqV65chqP4MtOkSZMMr4umTZsqLCzMKc7q1asVFhamyMhISRnnGF++fLk6duwoSXr88ccdr//XX389Q8wpU6aofPnyjmlrFi5ceNU8U1NTNW7cOFWoUMEx3Unnzp21du3aTLe/Wowr81+wYIEiIiJ0/Phxx78jIiI0ZswYvfLKK5Kkxo0bO9bndp7W7FzEeP78+SpcuLC6dOnitL53795KTk52eiwPPfRQphdbrFSpUoZpDnKyra9s3LhRdevWzXBRyNatWysuLu6qg8XZkX7hyazkxRzpr776qipVqpTlGT7PPvusihUrpldfffWq+4qIiJCkq9b7d999p7CwMO3atcux7n//+5/CwsIyTNnTq1cvRz1LGecY79u3rxYsWKCjR486TZNz5UDW0aNH1blzZ4WEhKho0aK67bbbsjUwvW7dOrVr105FihRR8eLFFRERoRdeeCHTo+6zE+PK/Dt27KixY8c6/h0REaHhw4crIiJCFy5c0GeffeZ4TLfffvtV8z116pRGjBjhGAzOiY0bN0qSrr/+eqf1derUUfHixb3y+cCdfNO/U13541VaWpomT56sGjVqKDg4WEFBQWrZsmW2poHZs2ePdu7cmeFo+lKlSun666/XvHnzHOty+nwxR7p38b07d/jenTcYSAeAbNi1a5fjlOCoqCidOXNGzZo1y3AqtPTvh8iWLVtqx44dWr58uRITE7V582bt2LFDnTt3VnJysp566inHl8unnnrK8St2Vkdk1q9f3+nInMTERL311lt68803HfN2jho1ynEUR2RkpGPb2267Ldt5ecv06dMVHh6uW265xXHkzhNPPKF169Zpzpw5/OINnxo1apSmTp2q6dOn68KFC/r55581adIkHTt2LMO2EydO1J133qk777xTsbGxOn/+vEaNGqWBAwfqu+++k/TvkS7pgxjpdfjMM89kGX/27NmO7S5evKhdu3YpJCREN910k2Pu6kOHDql+/frq2LGjY9tt27blKC9vuO666/T+++/riy++0Mcffyzp32kRhg4dqk6dOmn8+PFei43sSR+g7tq1a5bbpN+Wvm1W0s+2qFat2lXjdu7cWUuWLHEcgXX06FH9888/CgkJcRq0jIyMVEBAgNq2bZvpfjp27Kj169dL+ncQO/31f+Vra/bs2YqNjdWePXsUHR2thg0batCgQTpz5ozLPF977TW9//77+uqrr5SQkKATJ05o3Lhxmc7HnpsY/fr1U0xMjCpUqOD4d0xMjGbMmKGXXnpJ0r9TBKSvb9mypct8pf87E+TgwYN677339N1332nixIlXHUj/66+/VKVKlQzzt6cfnX75l+7MnD59WitWrFDdunWvmmNW22Z2cc6kpKQs13uKzWbLdJA7/bPP9u3b3dr/3XffreLFi6tw4cJq1KiRPvvsM7f2lxv79u3TgQMH1KlTpyw/TxYrVkytW7fW2rVrrzqvfvrr4Wr1nj5YfmVdh4SEaP369Y4j0BMSErRy5Up17tw5y31FRkaqX79+jimS0pfLX7OJiYl67LHH9Oabb+rixYtauHChvv/+e0c9ZSU+Pl7du3dXeHi44+jUlStXymazad26dU7b5jbG+vXr9d577zn+HRMTozlz5igmJkbBwcGO6X1iYmI0f/58l/uS/p0OKCwsLNPvG1djs9kkKdPXfUpKinbu3OnxiyTmJN/0I2fPnDmj1atXa9iwYapVq5YefPDBDPt8/vnnNWHCBMXHxys6Olo33nijOnXq5OgNWUl/DdeqVSvDbbVr19axY8ccRyb74vlC5vje7R6+d3sfA+kAkA0TJ06Un5+fvvzyS5UvX16FCxfW0KFDHUfrXO6VV17RyZMnNXfuXDVs2FA2m021atXSZ599pi1btrh9BEhKSoouXryozp07q0+fPo6j3q7G23m5EhISorlz5+ro0aMaNmyYvvnmG7377rt65ZVX1K5dO6/FBa5m//79mjlzph5++GENGDBAAQEBqly5sqZOnaotW7Y4bXvixAlNmDBB9957rx599FEVK1ZMQUFBeuCBB9SvXz89++yzbudz4cIF2e12vfTSS4qKirrqEcJ5lZcr9957r4YOHaqHHnpIK1euVP/+/RUcHKyvvvoqTy6qB9fSvzxWqlQpy23Cw8Ml/XuRvMulnwYcHR2tP//8U4MHD1ZqaqrLH4bSderUSWfOnNHmzZsl/d/A2n333ZdhwK1FixaZHtGcE2lpaXr66adVrFgxhYaGatKkSbpw4cJVf0hau3atrrnmGrVr1042m01BQUFq3769YzDMEzE87ZlnnlFERISqVaumhx9+WI888ogeeuihq94vOjo607MS0tddedHZK40YMUIXL17M1v9/VtuuWbNGYWFhTssXX3yhjRs3Zlg/c+bMq8bJrpYtW+qvv/7KMP3Nb7/9Jkm5PkIvMDBQzz33nP7++28lJSXpyJEjatmype666y6PTWuTXem1nl7PWQkPD1dqaqqioqKc1qdPMXTs2DHNmjVLH3zwga6//nqXA9+SVL58edWtWzdDXT/00EPy8/NzPMcrVqxQQkLCVS8sezUHDx7U448/roYNG8rPz0+dO3dWv3799Omnn7q8365duxQdHa1BgwY5jnquVKmSnn/+ebVv394jMTxpxowZWrBggaZPn66goKAst7t8eriYmBjHDyTpP8qlP//p1q1bp/j4eKWkpFx1+jVv5Jvu9OnTjotZt27dWidPntTXX3+tcuXKObbZsGGDZs6cqQkTJuiOO+5QUFCQQkJC9Prrr6tevXp64YUXJGX9HKS/p2XnfS+vny9kje/d7uF7t/fx7QYAsmH58uVq0aJFhtMNe/TokWHbX375RfXr11fJkiUVHx+v8+fPKy4uTuHh4SpevLjLKSOyYozRm2++qdq1ayswMFAVK1Z0nD62b9++bO3DG3nlRN26dfXxxx9r7ty5uv3229WzZ089+eSTXo0JXM2KFStkjFH37t2d1leoUEENGzZ0Wrd06VIlJSWpS5cuTjUUFxen5s2ba/fu3Tp9+nSOc9i5c6d69+6t0NBQhYWFqUqVKo5pW7JT397KKyc++ugj1apVS+3atdP27dv17bffOuY+hW+l/5jham7O9KPsChcu7LT+m2++UUREhMLDw9W4cWMtWbJEc+bMydZFdDt27Cg/Pz/H4FpkZKTat2+vHj16aPv27Tpx4oTi4uK0fv36qw7UZceV+6hRo4YCAgJ04MABl/dr3ry5Nm7cqBEjRmj58uUu51fPbQxPmzx5suMMlm+++UYfffSRevTocdX5V202W6bbpK9LPyIzM+PGjdP333+vESNGaODAgS7juNrW399fdrvdaSlcuLAKFSqUYf2VR867Y/z48TLGqH///tq0aZMOHDigZ599VsuXL3fklRv169fXiy++6BjgqVixoqZOnaoOHTropZdeuuqPE56UXuuZXWfjclnVe6tWrRxTGdx9991q1KiRIiMjs/WDaOfOnfX7778rOTlZu3fv1uHDh9WvXz+1bNnScaZLZGSkihQpkmFaspyy2+1q1qyZ07prr71Wp06dcnlh1Bo1aigsLEwPPfSQpk2bpkOHDnk8hqccOHBAjzzyiEaNGqU2bdq43PbWW291+gGqW7dukv6dMmrgwIH68MMPNWXKFB05ckTLli3TkCFDVKRIEUm5f927k2+6MmXKON7H9u/fr3r16ql169ZOBxD88ssvkv79Yfbyzzfnz59X8+bNHd9dsnoO0t/TsvO+l5fPF1zje7f7+N7tXQykA0A2xMTEqHTp0hnWZ7bu9OnT2rZtmypVqqRKlSopPDxclStXVpUqVXJ8McN0EydO1NixY/XEE08oNjbWcVGT4cOHZ/uiO97IK6c6deqkkJAQpaamasSIES6/tAN5If2U3uzUd/pg9F133eVUQ5UrV9arr74qu93udBG17IiNjVW7du107tw5rV+/XomJiY6pI6TMTzG+kjfyyqmiRYuqa9euSk1NVePGjbM1PQXyRtWqVSXJ5WDvwYMHJUlVqlRxWj9kyBDH0X0HDhxQixYtNHLkSO3du/eqccPCwtSkSRMtXrxYxhj99ttv6tKli2644QYVLVpUkZGRWrZsmVJSUjwykF62bFmnv202m4oWLXrVIwifeuopTZo0SatXr1a7du0cgzBXTvXgTgxvCQoK0oABA/T8888rMjLyqlPzlCpVKtOB3XPnzknK/H1Q+nfe7YkTJ+qOO+7QlClTXMa42rYtW7bMMIXLkCFD1LRp0wzrr5xj3x3XXXed1qxZo9KlS6t3795q3bq1du7c6Tgq8GpHcedUjx49lJSU5HStC2/LTq1L/9a7n5+fKlas6LR+x44divn/F/KdPn261q1bp3HjxmUrdufOnRUfH6/Vq1crMjJSpUuXVqNGjdSlSxenH9PatGnj9g8kV9ah9O+UNZJc1qLdbtfvv/+uRo0a6bHHHnOc1TF+/PgMP6DlNoanPPvssypatKjTNBTpcRMTExUTE+OYGqJYsWJOP0BdfnbP7Nmz9cILL+iTTz7RddddpzFjxuihhx7SzTffnGHbvMo3M1WrVtXXX3+t4OBgp+k70j/ftGvXLsPnmzlz5sjPz09paWlZPgfpP+i7et9LvyCvlHfPF1zje7dn8L3be/hJDQCyoXTp0jpx4kSG9ZldgKls2bIqXbq0VqxY4bH43377rW688Ubdc889Tuv379+f7X14I6+cSE1N1aBBgxQYGKgaNWrovvvu05YtW1SmTBmf5ANI//eh/MSJE455gtNdWd/pX6y///57p4uluWPZsmU6ffq0vvvuO11zzTWO9TmtbU/nlVOLFy/W5MmT1alTJy1ZskSvvvpqtqZ/gPd1795dM2bM0I8//qjHHnss020WLFggSS6PNI+IiHC8TocPH+44kteVzp07680339SqVat06tQpde7cWYULF9aNN96oyMhIhYaGym63q3nz5rl7cJfJ7RfEgIAAPfHEE3riiSd09uxZLV26VBMmTFDHjh21d+9ep2kG8uuX0PT3jiun5rlSw4YN9cUXX+jChQtO86nv3LnTcfuVJk+erGeeeUaDBw/WrFmzXB6dnJNtfaFRo0aaO3eu07r0qTqunNrDXenzu2c1B683VKpUSfXq1dPixYuVmJiY6YD12bNntXr1anXr1k0BAQGZ7sff31/Dhw/XP//8ozfffFM9e/a8am9p27atChcurEWLFmnnzp3q1KmTbDabOnfurOeee05r167V1q1bNWnSJLcfpzt12KBBA/3vf/9TWlqatm7dqq+//lovvviiLly4oDfffNMjMTwhNTVVCQkJql+/vmNd+hHUM2fO1Ndff63PP/9cvXr10pw5c7Lcj7+/v8aOHeu4AGq6N954w6Ov+Zzkm5XAwEBFREQ4vY+lf77ZunVrhh96L5fVc5D+npb+Hne5nTt3qnr16k6D43n1fME1vne7j+/d3pW/Pt0AQD7VuXNnrVmzJkMDz2x+s759+2rDhg3avXt3pvu62qnXmcnsy86hQ4cyzOOX/mEws1PTvZFXTowbN04rVqzQ7Nmz9f333yshIUG33XbbVU9BBrypXbt28vf3zzC4snfvXm3dutVpXefOnRUcHKwvv/wy0315qrYlZTo3cEhISKa17Y28cuLIkSO6/fbb1aZNG/3666965JFHNH78eC1dutSrcZE9ffr0UePGjfXaa6855lC+3N9//60pU6aoQ4cOGaYyuFL58uX1+OOP648//tAPP/xw1didO3dWUlKSnnnmGVWpUsVxwbcuXbo4jqBu3779VQcbXfU2d11eHyVLltSAAQM0fvx4XbhwIdPBF0/K6ePKqpb/+OMPSf8OEqZLv3jn5fseMGCAUlJSMlzk8JtvvlFwcLBuuukmp/UfffSRHnvsMQ0aNEhffPGFy/+nnGzrbZk99sykpKTogw8+UKdOnTL9ESE7Mvs/SU1N1ffff69ixYrpuuuuy9V+c+vFF1/UmTNnNGHChExvf/zxx5WcnJytU/yfffZZlShRQk8++eRV+0hwcLBatWqlX375RcuWLVOXLl0kSc2aNVNYWJjGjh0rY0y2zj7Jqte56/LH4Ofnp8aNG2vSpEm69tprr3rRSk/IyeNKv0Dp5cuOHTsk/XuRw5iYGJeD0q7MmzdPhw8f1qOPPuq0Pn2e8dzISb5ZvZZOnz6tf/75x+l9rHfv3rLZbLn+fFO5cmU1b95c33//vdOFQg8dOqR169bp1ltvvepjy+r5gvfwvdt9fO/2LgbSASAbnn76aYWEhKh3795au3atDhw4oJdffjnTZjRu3DjVrVtX3bt315w5c3Tw4EHt27dPCxYsUJ8+fbRw4cIcx+/fv7+WL1+ujz76SMeOHdPSpUs1YMAA9ezZ02m7OnXqKCAgQHPnztWJEycUExPjOCrK3bzST2tLnyIi/SJ0VztdU5J+/PFHvfHGGxo/fry6dOmiKlWq6PPPP3cc+Xc17sQGXKlQoYIef/xxffTRR5o8ebIOHjyoFStWaNSoUWrbtq3TtiVLltSUKVP0xRdfaOTIkdq0aZNOnDih9evX6/XXX8/Vl9rWrVurXLlyeuqpp7R161bt3btXzz//fKav60aNGunPP//U5s2bFR0d7Tht2t28UlJSHPWUfqpp+kXnzp8/7/K+ycnJGjhwoPz9/TVnzhwVKlRIkyZNUosWLXTbbbdlevSQp2IjewoVKqR58+YpLCxMrVu31uzZs3Xs2DEdOXJEM2bMULt27VS1alV988032drfI488otKlS+vZZ5+96hey66+/XsHBwfrjjz+cBtA6d+6sEydOaPfu3dkaWCtfvrzKlCmjhQsX6siRI9kaJM2unj176rXXXtOGDRt08uRJbd68WVOnTlXZsmXVpEkTj8TISqNGjST9Owh1+vRpxcTEOA32XGno0KGaNGmSNmzYoBMnTmj79u167rnnNHnyZN11111OP4Rs2LBBYWFhTj22e/fu6tmzpx599FEtWrRIJ06c0JQpU/Tpp5/q1VdfdZqP9vPPP9f999+vm266SR988IHOnz/vNEh2+SBATrZ1x5XvDQkJCY4Ylz9vmT323bt3q3///lq5cqVOnjypNWvW6Oabb1Z0dHSmP1xmN9bdd9+tiRMnav369Tp+/LjWrFmjvn37auPGjZo8ebJjOpC80rdvX7388st6/fXXde+992rz5s06ffq01q9fr1tvvVVffPGF3n333Qz9LTN2u11PPPGENm/enK0L43Xu3Fl//vmn4uLiHHXt5+enDh066I8//lCZMmWcBkmz0qhRI506dUq//vqrzp07l+vB3SstXLhQvXv31o8//qj9+/fr2LFjmj59uvbs2ZPp3Mue1qhRI61YsUL//POPYmJidOHCBa/HfOSRRzR16lTt27dPhw4d0ocffqg777xTL774Yoa5zG+55RaFhYU5vbdm9nk7/e/c5r9mzRr17NlTCxcu1L59+3T48GEtXLhQXbt2lb+/vyZOnOjYtn79+ho/frxeeOEFPffcc9q+fbuioqK0atUqjRs3TiNGjLhqvLffflsHDx7UyJEjdejQIW3btk2DBg1S5cqV9cQTT+T6+YL38L2b7935ngEAZMvu3btNv379TLly5UytWrXMa6+9Znbs2GHsdrv55ptvnLa9ePGimThxomnevLkpU6aMqVWrlunTp4/58ccfTVpamjHGmPj4eGO3283EiROd7rtw4UJjt9vNhg0bHOtSUlLMxIkTTb169Uz58uVNt27dzObNm82zzz5rSpYs6XT/zz77zDRq1MiULFnS2O12M3v27BzllZVOnToZu92e6XLl47/c0aNHTUREhOnZs6dJTU11uu3ZZ581YWFhZvny5V6JDWRHWlqaeeedd0y9evVMuXLlTPfu3c3OnTvN0KFDTePGjTNsv3btWjNo0CBTtWpVU7ZsWdOyZUszbtw4ExUV5djm2WefNXa7PcN9b7jhBnPLLbc4rdu+fbu56aabTMWKFU3t2rXNCy+8YGJjY43dbjdvvPGGY7tTp06ZAQMGmEqVKpnixYub+vXr5zivzCxZsiTL+srs8V8uqxpOr/tevXq5fG9xJzZyJj4+3kyaNMlcf/31pnz58qZo0aJGkhk0aJBJTEx02jYmJsbY7Xbz4IMPZrqvd955x9jtdjN37tyrxr3lllsy3bZWrVrGbreb3bt3O60/dOiQsdvt5vPPP3daHxkZaVq1amVKly5t7Ha7mTBhgjHGmG3btmWZS+XKlc0jjzziMr+oqCjz7LPPmubNm5uyZcuaunXrmtGjR5sDBw44tslJjKzyr1Onjhk2bFiG+0+cONFce+21pkSJEsZut5s1a9ZkmeuZM2fMyy+/bFq1amXKli1rKleubLp3725mz56doc42bNjg9Dylu3TpkpkwYYKpU6eOKVOmjGnZsqX58ssvM8Ry1Xftdrs5e/ZsrrbNzIMPPmg6derkchtjjKlWrVqWMS7/zJTVY//hhx9M+/btTbly5Uy9evXMU089Zc6dO+dWrDNnzpiXXnrJ8X8SERFh+vfvn+nnmvfff99IMsePH7/qY3XXmjVrzB133GFq1KhhSpQoYfz8/ExgYKBZv359hm2feuopY7fbzbFjxzLcduHCBVO9enXTokWLDJ/hrpT+vDdr1sxp/YwZM4zdbs/09T9q1ChTp04dp3UJCQlmxIgRpmrVqqZ48eLGbrebhIQEY4wxt912W4b9G2PMlClTjN1uNydPnswyv7S0NLNgwQLTp08fU61aNVO5cmXTvn37DLWakxiZ5f/5559n+t62a9cu061bN1OuXDlTvHhx07t37yxzzcyxY8eM3W43Tz31VLbvc+rUKfPQQw+Z2rVrm4oVK5qbbrrJ/PLLL5lum/5enf5cG2PMJ598kmUdXPlZJif5LlmyxAwcONDUqlXLlC5d2jRo0MCMGTPGHDp0KNN9/fzzz6Znz54mPDzcVKhQwbRp08a8/vrrJjo6OlvPw9q1a02PHj1M+fLlTZUqVczw4cMz/WyUk+erZMmSZsSIEdmKj5zjezffu/MzmzFePqcAAAAAQL5y4cIFdejQQTt37tRvv/3mkXnKAWTtgw8+0IMPPqgTJ05kejFLb/rmm2902223aeDAgZo9e3a+m78eKGhKlSqlAQMG6KOPPvJ1KgDyGB0UAAAA+I8JDg7Wzz//rBo1amjw4MHau3evr1MCLC0uLk5+fn4qWbJknse+9dZb9f777+uXX37R008/nefxAas5f/6844L1AP5bOCIdAAAAAAAvSE5O1uHDh3XrrbcqKChIK1eu9HVKAHIpOjpac+fO1b333qvIyEh16tTJ1ykByGP+vk4AAAAAAAArWrBggR588EE1btxYb7/9tq/TAeCGbt266cKFC3rjjTcYRAf+ozgiHQAAAAAAAAAAF5gjHQAAAAAAAAAAFxhIBwAAAAAAAADABQbS4VH79+9Xt27dtHz5cl+nAhRot956q95//31fp5Gp/JwbkN8tWLBA3bp1U3R0tK9TySA/5wYUBPm5P+bn3ID8Lj/3x/ycG1AQ5Of+mJ9z+y9jIN1LLly4oG7duumdd97J9Pb169erW7duioyMzNvEvCwuLk6LFi3S8ePHfZ0K4DWfffaZunXrpldeecWt/fTt21dTpkzJ9LbffvtNO3bscGv/3pLXuR0/flzdunXLchk5cmSe5QJrS0hIUK9evdStWzcdOnQo1/uZO3euunXrpri4uAy3HTp0SIsWLVJiYqI7qXqFL3O7dOmSPvnkEw0ZMkQDBw7U5MmTdeHChTzPA9ZF78773MaPH59p3x4yZEie5gFro3f7Prf0sY30C3ECnkLv9k1uR48e1dtvv62ePXvqpptuyvP4+Z2/rxOwquTkZC1atEiVKlXK9PZTp05p0aJFGjRoUB5nBsBdr7zyig4ePKiVK1dqzJgxCgkJydV+IiMjVbFixUxv+/bbb1WuXDl30vSavM7Nbrfr4YcfzrB+5cqVeuWVVzRs2LA8ywXWNnfuXC1YsEABAQGaMWOGXnzxxVzt58CBA1q0aJGSkpIy3NarVy/VqFFDJUqUcDddj/NVbvv371fXrl1VpkwZDRs2TKVLl9batWvVoUMHrVu3Lk9zgXXRu/M+t02bNmnHjh2aNm2a0/qgoKA8zQPWRu/2bW4XL17UkCFDdODAASUnJys5OdknecCa6N15n9sTTzyhb7/9VrfccotOnjypzZs352n8goCBdADIgd9//1179uzR1KlTNXr0aH311VcaMWKEx+N06NDB4/v0lLzOrWjRourWrVuG9bNmzZIk3XfffXmaD6xr2rRpatWqlWrXrq2ZM2dq/PjxKlSokEdjVK5cWZUrV/boPj3FF7klJyerT58+ql69un766SfH892rVy9OU4fH0Lt9l1twcHCmPRzwFHq3b3N78sknlZaWpgEDBuirr77yWR6wHnq3b3J76KGHNGnSJNlsNvXp04eB9EwwkJ4PrVu3Tl9//bUOHDigkJAQde/eXYMHD5af3//NxLNlyxZ98cUXOnTokMLCwtSqVSsNHTpUAQEBOdomO/H279+v0aNH6+mnn1bFihX13nvv6fDhw6pTp44eeeQRlS5dOtPHsXPnTn3wwQeKiopS/fr19dBDD6lUqVI5fryXxy9fvrw+/PBD7du3T6+++qoaNGigM2fO6L333tO2bdtUrlw53X///YqPj9cLL7ygd955R9dcc43b/ydAuunTp6tmzZq677779OOPP2r69OlZNvQDBw5o1qxZ2rFjh4KCgnTjjTfq7rvv1qVLlzRw4EAlJCRowYIF2rt3rySpcePGmjhxoqR/50O74YYb9OCDD0qShgwZovDwcL366qsZ4owbN07Hjh3TZ5995li3f/9+zZw5Uzt37lShQoXUunVr3XvvvQoODnb5+I4cOaLp06dr165dKlSokBo1aqR7771XxYsXd2xzZW7vv/++fvrpp0z3165dOz311FNu53Wlc+fOaf78+WrQoIFatGiRo/sCmdmzZ4/++OMPff7557rmmms0a9Ys/frrr5mezpiamqq5c+dq8eLFOnPmjGrUqKG7775bdevW1ZQpUzR9+nRJ/9ZKes996623VLduXS1YsEAffvihvv76a4WFhenXX3/VO++8o8mTJ6tOnTpOcXbu3KlHH31Ujz32mDp37uyI/d1332nRokU6e/asKlWqpCFDhqhVq1YuH19aWpq+//57LVq0SKdPn1blypXVu3dvpw/oV+YWHR2twYMHZ7nP6dOnO868y21e3377rbZv366ZM2dmGPgICwtzeV8gu+jd+aN3A55G7/ZN7063ePFiffTRR/rtt9/05ZdfZus+QHbRu33Tu7OaVQP/hznS85kXX3xR119/vS5duqTbb79dTZs21SOPPKJ+/fopLS1NkvTLL7+oadOmunjxogYPHqy2bdtq/fr16tevn2M/2dkmu/HS5z1ftmyZHnjgATVq1Eg333yzPv/8c7Vv316pqakZHsfWrVv1wAMPqFmzZurevbu+/fZbXXfddTp27Fiu4y9evFgjR45UgwYN1Lx5c8XExCgqKkrXXXed5syZo27duqlVq1Z67LHH9NNPP2nRokWKiYnx5H8P/uOio6M1d+5cjRw5UjabTaNHj9bGjRv1559/Ztj2u+++U506dbRq1Sp17dpVPXr00JYtWzR48GAVKVJEDz/8sAICAtSwYUM9/PDDevjhh3XLLbc47n/lfGilSpXSm2++qdOnTzvFOX36tN58802VLVvWsW7evHmqW7eu/vzzT/Xt21fdunXTjBkz1KJFC5dHeO7bt08NGjTQhg0b1KtXL/Xu3Vvnzp1Ts2bNnLa7MreOHTs6HkP6UrduXS1atEj79+93O6/MzJ49W4mJibr33ntzdD8gK9OnT1fJkiU1YMAANW3aVE2bNs0wHYH07zVQOnTooOHDh6tChQq64447VLFiRd1xxx3asGGD2rdvrzZt2kiSRowY4aiJ9NNJr5zLtEmTJlq2bFmmsT755BP9/vvvatKkiaR/T5/u2LGjRo8erZo1a2ro0KGOLwsff/yxy8f30EMPafjw4apRo4buuusu1apVS6+99prTBYyuzC04ODhDbd9///3asWOHIiMjZYxxO6+ffvpJoaGhCgsL09ixY9W7d28NHTpUs2fPdnwOANxB7848t7zq3adOndJtt92mXr16adSoUVkOAAC5Qe/2Te+W/n1vHTZsmO655x61a9fuqtsDOUHvzjw3X3zvRiYMvCI6OtpIMpUqVTJdu3bNsFx33XVGkvn0008d91m6dKmRZN58802nfW3atMnYbDYza9YsY4wxt9xyi7nuuusyxDx79qzj39nZJrvxtmzZYiSZa6+91iQkJDi2i4yMNJLM3LlzHevSt61Ro4a5ePGiY/3JkydNsWLFzKBBg3Idv3r16k77TExMNHfccYcJDg42J06ccFpfr149I8msWbMmw3MA5Na7775rihQpYs6dO2eMMSY1NdVUrVrVjB492mm7Y8eOmSJFiph+/fqZtLQ0p9vS72uMMcHBweb+++/PNFbJkiXNiBEjHH//9ddfRpJ56623nLZ78803jSTz999/G2OMOXHihClatKhTrRnz73tSyZIlzQMPPJDl45s4caLx9/c3iYmJTusvf9/ILLcrbd++3YSGhpqaNWuaM2fOuJ1XZho1amSKFClioqOjc3Q/IDNJSUmmbNmyZuzYsY51M2fONP7+/iYqKspp2zFjxhg/Pz+zefNmp/UpKSkmLi7OGGPMG2+8YSSZ06dPZ4j1/vvvG0nm+PHjjnX9+/c3pUqVcqq9xMREU7JkSaeaeeyxx0xAQIDZtm2b0z5feuklU7hw4Qy5pktNTTVFixY1zz//fIbbLq/vzHK70rBhw4wk8/7777udlzHGNG3a1ISEhJiQkBAzcuRI8+2335pnnnnGBAQEmJ49e5qUlJQs7wtkB70789yu5I3e3bdvX3P33XebTz/91HzxxRdm6NChxmazmVtvvZXahtvo3VnndiVP925jjBk8eLCpUKGCiYmJMcYYM3z4cCOJz+bwCHp35rldydvfu3v37m0KFSqU7e3/KxhI95L0gfSuXbuaX375JcPywgsvZBhIv/POO01AQIC5dOlShv3Vr1/f9OnTxxhjzF133WVCQ0PN0qVLs4yfnW2yGy99IPuFF15w2iYpKcnYbDYzYcIEx7r0bZ977rlM4xUpUsQkJyfnKv64ceOctklNTTUhISFm8ODBGe7/yiuvMJAOj6tfv765++67nda99tprpnjx4k4/8rzzzjtGklm/fr3L/eWkoRtjTIsWLUydOnWc1tWpU8e0bt3a8Xf6h+m1a9dm2Ofdd99tKlWqlGU+H3zwgZFkpk+f7qjT7OaWLioqylSuXNmUKlXK7N271yN5XWnjxo1GkrnzzjuzfR/Alblz5xo/Pz9z4MABx7qLFy+aEiVKmFdffdWxLi0tzYSFhZkePXq43F9Ov4z/8ssvRpL59ttvHeu++eYbI8lERkY6YpcqVcp069Ytwz4PHTrkqN2slC5d2rRs2dLpMWYnt8u99NJLRpJ5+OGHHevczatu3bpGknnwwQed1r/33ntGkvniiy+yvC+QHfTurHNL563enT7Adrn0gYiPP/7Y5X2Bq6F3Z53b5bzRu+fMmWMkmR9++MGxjoF0eBK9O+vc0uXF924G0jPHHOleVqlSpUwvsJOSkpJh3d69e+Xv76/+/ftLksy/P3TIGKOoqCgFBgZKksaPH69t27apQ4cOKl++vG688UZ169ZNt912mwoXLpztbbIbL13VqlWd/g4ICFBISIhOnjyZ4bFcua0kVatWTZcuXdKpU6dUoUKFHMevVq2a098xMTE6f/58prEiIiIyrAPcsW7dOm3fvl02m82ppuPi4hQTE6PvvvtOQ4cOlSTHaVW1atXyaA733HOP7r33Xq1bt04tWrTQ2rVrtXPnTs2cOdOxTfq8b88//7wCAgIcNSVJu3btUlRUlIwxstlsGfY/bNgwRUZG6p577tGjjz6qG264QZ06ddJdd92VrbmKL1y4oJtvvlmnTp3Sb7/9purVq3skryulP14uMgpPmTZtmoKCgjRy5Ein9YUKFdL06dP11FNPyWazKSYmRtHR0R6v7S5duig8PFwzZ87UgAEDJP37Oo+IiFDHjh0lSbGxsTpz5oz++ecf3XzzzY76McY4pkA5evRoljGmT5+u4cOHq2rVqqpbt67atm2r/v37q3379tnK8auvvtLzzz+vPn366K233nKsdzev0NBQSdKgQYOc1t92220aM2aMFi1apDvuuCNbOQJXonf7tnfb7fYM6x588EE99dRT+uGHH+jjcAu9++q80bvPnDmj0aNHa+DAgerVq1e28gBygt6df753I3MMpOcjhQsXVkhIiB544IFMb0//shkREaFNmzZpx44d+v3337V06VLdd999eu+997RmzRoFBgZma5vsxkt35UVKJclms2U6h2l8fHyW64oUKZKjx5vuyosipA+0u4oFeMq0adPUsmVLjR8/PsNts2bN0rRp0xwNPf01HhcXl+mXyNwaNGiQHnnkEcf8ZjNmzFBISIgGDhzo2Cb9h7Lhw4dnqKGrKVKkiObPn68jR45o6dKl+uOPPzRhwgRNnDhRGzZsUJUqVbK8b2pqqgYNGqQtW7Zozpw5uv76651udyevyyUkJOirr75S3bp1M8QAcuPw4cNavHix3nvvPacPoZKUmJiofv36aenSperYsaOCgoIk/VvbnuTn56e7775bL7/8so4ePSpjjCIjIzV+/HjHh9z0Gqpfv75Gjx6dYR8PP/xwhvwv16tXL0VFRWnVqlVauXKlFi5cqClTpuiJJ57QpEmTXOa3YsUK3X333WratKlmz57tdPFzd/OqXbu21qxZk+FLg91ul81m8/hzjf8Wenf+6N1X7jM0NFSxsbEe2R/+m+jdvuvdBw8e1Llz53Tw4EGnQc70eZz79++vGjVqaOrUqS7zA7JC785/vRtX8Orx7v9h6VO7DB8+PNPbFyxYkGFql2effdbYbDZz5MiRHMdLPzVkxYoV2d4mu/HSp1b5+uuvM9xmt9udTjVJ3/b222/PsG2bNm1MtWrVHH97In7t2rVNq1atMqy/5557mNoFHnP+/HlTrFgx88Ybb2R6+6+//uo0X9oPP/xgJJnPPvvM5X6LFy9uRo4cmeltWZ3GNWzYMBMaGmpOnTplQkJCzD333ON0+/z58zOcaumOtWvXGknm7bffdpnb/fffbySZ119/PdP9eCqv2bNnG0nmnXfecWs/QLrx48cbu92e5WmVLVu2dJpjsHHjxqZmzZou5/dNP8308ut3pMvqFOyDBw8aPz8/89JLL5kXX3zR+Pn5mUOHDjlt06hRI3PddddlmAMyt7p162bsdrvL3Hbt2mVKlChhIiIiMn087uaVfhr8//73P6f127ZtM5LM448/nuN9AsbQu/NT777c3r17jaQMp+wDOUHvtrvMzZu9OzY2NtOpa7t27ero5ytXrszVYwPo3fmrdzO1S+YYSPeS3Ayknzp1ypQvX960a9cuw8U9Fi1aZH7++WdjzL8XXti3b5/T7c8884yx2Wxm9+7d2d4mu/FyM5Bev359x9xwxhjz2WefZbjAiSfif/TRR0aS48KkxhizYsUK07hxYwbS4TGffPKJkWS2b9+e6e2XLl0yRYoUMY899pgx5t/5+5s3b24qVqzodFGj6OhoM3PmTMffDRs2NF26dMl0n1k19FWrVhlJpmPHjpm+xlNTU831119vIiIizNatW51u+/PPP53iX2n27NkZ5lKbN2+ekWTmzJmTZW6TJ082ksx9992X5b7dyetyHTp0MEFBQU4XjwFyKzU11YSHh5tbbrkly21eeOEFExgY6LiAz48//mhsNpt56KGHTFJSkmO7RYsWmb/++ssYY8z3339vJJnVq1dn2J+ruUw7d+5sqlWrZqpWrWq6du2a4faffvrJ+Pn5maefftrp4kQXL140H330kTl8+HCmj+HMmTNm8uTJjguqGfPvBdGaN29urrnmmixzO336tKlevbopXry42blzZ5bPUW7zMsaY5ORkU79+fdOwYUPHvLQXL1403bt3N0WLFjV79uzJ8r6AK/Ru3/bu3bt3m88++8xpoPP48ePmhhtuMIGBgebPP//M8r6AK/Ru3/fuzDBHOjyB3p1/vncbw0B6VhhI95LcDKQbY8yePXscA0WNGzc2rVq1MqVLlzbdunVzFMGsWbNM7dq1TbVq1Uy7du1MzZo1TdmyZZ0uCJKdbbIbLzcD6Z9++qnp27evadiwoalTp44pVKiQeeKJJzL84u1u/LS0NDN27FhTqFAhc+2115rGjRub3r17m88//9xIMhs2bMjqvwjItmbNml31ohw9evQwpUuXdnwQPXPmjLnllluMv7+/qVu3rmnWrJmpWLGiU+P66quvjL+/v2nYsKHp2rWreeqppxy3ubqwyLXXXmskmbp162Z6e3R0tBk6dKgpXLiwqVOnjmnTpo2pVKmSadKkifnpp5+yfAzLli0z119/valQoYJp06aNadSokQkJCTFPPvmkU+1emVuTJk0cHzK6du3qtEycONHtvNIdOHDA2Gw2c8cdd1x1WyA7fvrpJyPJTJs2Lctt1q9fbySZyZMnO9Z9/fXXpnz58qZkyZKmdevWpnLlyqZ79+6OL7GJiYmmcePGJiwszHTo0MF07drV8UXd1Zfx9KOzJZnvvvsu03zmz59vqlat6ojdoEEDU7p0aTN69GgTGxub6X0uXbpkHn30UVO2bFlTv359065dO1OqVClz3XXXOX3puDK39M8qERERGWq7a9euTmeU5SavdIcPHzatW7c2drvdtG7d2pQpU8ZUqVLF5QXTgauhd/u2d589e9bcc889pkyZMqZp06amadOmpkiRIubaa6+ltuEWenf+6N1XYiAdnkDv9v337h9//NGxzzJlyhibzeb4+5lnnnF53/8KmzH/f9Z5eFRKSoqWLFmi8PBw1a1bN8Ptp0+f1qZNm1S/fn1VrFgxw+0nTpzQrl27FBQUpFq1amWYO9QYo4MHD+rQoUMqVaqUateunWEO8+xsk51458+f16pVq9SoUSOVK1fO6X6//fabypcvrzp16mS67a5du3T8+HFde+21Klu2bJbPV27jpzt58qR27typ8uXL65prrtGHH36oBx54QHv37nU57xxwNWlpaYqMjFSpUqXUpEmTLLfbs2eP9u3bp9atWyskJMSx/tSpU/r7778VEhKiOnXqOOZpTHf27Fn99ddfSkhIUMmSJdW0aVNJ0tKlS1WuXDlHbV1u586dOnz4sCpXrpzp7eliY2O1Y8cOpaamqnr16qpQoUK2HvOJEye0Z88eFStWTDVr1lSxYsWcbr8ytzVr1mQ512n58uXVsGFDj+QVFRWlbdu2Zfm+CeTU7t27tX//frVq1SrLeRWNMVq8eLHKlCmjxo0bO9anpqbqr7/+UnR0tGrWrJnhNZmWlqa//vpLJ0+eVGpqqlq2bKnixYvr8OHD2rlzpzp06OCYwzBdUlKSli5dKkmZ3n55TukXCypTpoxq1KiR4b0lM4mJidq9e7eio6MVHh6e4WLdV+aW/lklK23atHG6fklu80q3d+9eHT58WCVLllT9+vWd5nMFcoLenX9698WLF/XPP/8oNjZWVapUUbVq1bJ1PyAr9O781bvT7dixQ0eOHFGnTp3k78+l+JBz9O780bvT31MyU7p0aZf/N/8VDKTDkm6//Xb98ssvOnPmDF/EAQAAAAAAALiFEUYUeDNmzFBKSorj7/nz5+ubb77RyJEjGUQHAAAAAAAA4DbOuUGBt3XrVlWsWFE1a9bUyZMndeDAAQ0bNkwTJkzwdWoAAAAAAAAALICpXWAJcXFx+vvvv5WUlKS6deuqRIkSvk4JAAAAAAAAgEUwkA4AAAAAAAAAgAtMIA0AAAAAAAAAgAsMpAMAAAAAAAAA4EKeX2w0LS1NUVFRCgkJkc1my+vwgOUYY3T+/HlVqFBBfn6++22M2gY8j/oGrInaBqyL+gasKz/UN7UNeF5OajvPB9KjoqIUHh6e12EByzty5IgqVarks/jUNuA91DdgTdQ2YF3UN2BdvqxvahvwnuzUdp7/hBYSEpLXIYH/BF/Xlq/jA1bm6/rydXzAqnxdW76OD1iZr+vL1/EBK/NlfVHbgPdkp77yfCCdU08A7/B1bfk6PmBlvq4vX8cHrMrXteXr+ICV+bq+fB0fsDJf1he1DXhPduqLi40CAAAAAAAAAOACA+kAAAAAAAAAALjAQDoAAAAAAAAAAC4wkA4AAAAAAAAAgAsMpAMAAAAAAAAA4EKOB9J37typ++67T127dtUDDzygQ4cOeSMvAHksNTVVH374oXr27Km+ffvqs88+kzHG12kB8AB6N2BN9G7AuujdgDXRu4GCLUcD6f/8849atmyp1NRUjRw5UsePH1fz5s11/Phxb+UHII+MHDlSr776qvr3769u3brpoYce0vjx432dFgA30bsB66J3A9ZE7wasi94NFHAmB2677TbTsmVLx9/JyckmIiLCPPbYY9neR2xsrJHEwsLi4SU2NjYn5exk165dxmazmV9//dWxbvr06aZw4cLmzJkz1DYLi48Xd+qb3s3Ckn8XejcLi3UXejcLi3WX3NY3vZuFJX8v2antHB2RvnjxYvXq1cvxt7+/v26++WYtXrw4J7sBkM9ERkaqSJEi6tixo2Ndnz59lJSUpOXLl/swMwDuoncD1kTvBqyL3g1YE70bKPj8s7vhxYsXdebMGVWsWNFpfcWKFV3O15aYmKjExETH33FxcblIE4A3HTp0SGXLlpW///+9JZQsWVJBQUFZ1je1DeR/9G7AuujdgDXRuwHroncDBV+2j0hPSkqSJBUpUsRpfdGiRR23ZWbixImy2+2OJTw8PJepAvCWpKSkDLUt/VvvWdU3tQ3kf/RuwLro3YA10bsB66J3AwVftgfSQ0JC5O/vr7NnzzqtP3v2rMLCwrK839NPP63Y2FjHcuTIkdxnC8ArwsLCMtR2amqqYmJisqxvahvI/+jdgHXRuwFroncD1kXvBgq+bE/tUqhQIdWvX1+bN292Wr9hwwY1btw4y/sFBgYqMDAw9xkC8LrGjRvr1KlTOnbsmOM00k2bNskYk2V9U9tA/kfvBqyL3g1YE70bsC56N2ABObnC8HvvvWeKFy9udu3aZYwxZt26dSYgIMB899132d4HVxhmYfHOktsrhxtjzKVLl0zFihXNqFGjjDHGpKWlmT59+ph69eqZtLQ0apuFxceLO/VN72Zhyb8LvZuFxboLvZuFxbpLbuub3s3Ckr+X7NR2jgbSU1JSzLBhw0xQUJBp0KCBCQwMNE888UROdkHRs7B4aXHnw7oxxqxevdqUL1/eREREmAoVKpjq1aubHTt2UNssLPlgcae+6d0sLPl3oXezsFh3oXezsFh3cae+6d0sLPl3yU5t24wxRjl09OhRHTlyRNWqVVPZsmVzdN+4uDjZ7fachgRwFbGxsQoNDXVrH8nJydq6dav8/f1Vv359FSpUKNv3pbYB7/FEfdO7gfyH3g1YF70bsC5365veDeRP2antbM+RfrlKlSqpUqVKuUoKQP4VEBCgpk2b+joNAF5A7wasid4NWBe9G7AmejdQcPn5OgEAAAAAAAAAAPIzBtIBAAAAAAAAAHCBgXQAAAAAAAAAAFxgIB0AAAAAAAAAABcYSAcAAAAAAAAAwAUG0gEAAAAAAAAAcMHf1wkAAAAAAAAAyJ7Y2FiFhoZ6PY7NZvN6DKAg4Yh0AAAAAAAAAABcYCAdAAAAAAAAAAAXGEgHAAAAAAAAAMAFBtIBAAAAAAAAAHCBgXQAAAAAAAAAAFxgIB0AAAAAAAAAABcYSAcAAAAAAAAAwIVcD6QnJCQoLS3Nk7kAyAdSUlKUlJTk6zQAeAG9G7AmejdgXfRuwJro3UDBlKOB9NjYWL3//vuqW7euihQpoq+++spbeQHIY2vXrtWdd96pkJAQ1alTx9fpAPAQejdgXfRuwJro3YB10buBgi1HA+lfffWVdu3apW+++cZb+QDwkaefflrt27fXmDFjfJ0KAA+idwPWRe8GrIneDVgXvRso2PxzsvGoUaO8lQcAH1u2bJkkacKECb5NBIBH0bsB66J3A9ZE7wasi94NFGxcbBQAAAAAAAAAABdydER6biQmJioxMdHxd1xcnLdDAsgD1DZgXdQ3YE3UNmBd1DdgTdQ2kL94/Yj0iRMnym63O5bw8HBvhwSQB6htwLqob8CaqG3AuqhvwJqobSB/8fpA+tNPP63Y2FjHcuTIEW+HBJAHqG3AuqhvwJqobcC6qG/AmqhtIH/x+tQugYGBCgwM9HYYAHmM2gasi/oGrInaBqyL+gasidoG8pccDaSnpaUpKSnJ8XdycrISEhLk7+8vf3+vj8kD8KKkpCSlpaUpJSVFxhglJCRIkoKCgnycGQB30LsB66J3A9ZE7wasi94NFHAmB5YuXWoCAwMzLCNHjsz2PmJjY40kFhYWDy+xsbE5KecMOnTokGl9Hz16lNpmYfHx4k5907tZWPLvQu9mYbHuQu9mYbHu4k59e6p3u/sZIrt8/VyzsOTlkp26sv3/wsgzcXFxstvteRkS+E+IjY1VaGioz+JT24D3UN+ANVHbgHVR34B1+bK+02s7r3Kw2WxejwHkF9mpK69fbBQAAAAAAAAAgIKMgXQAAAAAAAAAAFxgIB0AAAAAAAAAABcYSAcAAAAAAAAAwAUG0gEAAAAAAAAAcIGBdAAAAAAAAAAAXGAgHQAAAAAAAAAAF/x9FTg2NlahoaFejWGz2by6fwAZ5UVtS9Q34Av0bsCa6N2AddG7AbjDGJNnsXgvQUHAEekAAAAAAAAAALjAQDoAAAAAAAAAAC4wkA4AAAAAAAAAgAsMpAMAAAAAAAAA4AID6QAAAAAAAAAAuMBAOgAAAAAAAAAALjCQDgAAAAAAAACAC/45vcOlS5e0detWJScnq27duipRooQ38gLgAwcPHtSePXtUoUIF1a1b19fpAPAQejdgXfRuwJro3YB10buBgitHA+mvvPKKpkyZovDwcBUqVEhbt27VSy+9pEceecRb+QHIAzt37tTo0aN16NAh1ahRQzt27FClSpU0b948VapUydfpAXADvRuwJno3YF30bsCa6N1AwZejgXS73a6dO3fKbrdLkubOnav+/furTZs2atq0qVcSBOB9Z8+e1Ysvvqgbb7xR0r9HwLRr106jRo3SggULfJwdAHfQuwFroncD1kXvBqyJ3g0UfDmaI/2BBx5wNHNJ6tevnwoXLqzNmzd7PDEAeadNmzaOZi5JRYoUUe/evbVp0yYfZgXAE+jdgDXRuwHroncD1kTvBgq+HM+Rfrm1a9cqKSlJ1157bZbbJCYmKjEx0fF3XFycOyEB5JHly5dT24AF0bsB66J3A9ZE7wasi94NFDAml6Kjo03t2rVN9+7dXW43fvx4IynDEhsbm9vQ2ZZZXBYWqy6erKkpU6aYQoUKmZUrV2a5jS9r2xjqm+W/tXiqrujdLCz5a6F3s7BYd6F3s7BYd/FUXRWE3p2XfP3/ysKSnbqy/f8Xa47Ex8erS5cuSkxM1G+//abixYtnuW1mv56Fh4crNjZWoaGhOQ2dIzabzav7B/ITT9XUt99+qzvuuEPTpk3TnXfemeV2vqxtifrGf4sn6oreDeQ/9G7AuujdgHV5oq4KSu/OS7yXwNeyU1c5ntolPj5e3bt3V0JCgpYsWeKymUtSYGCgAgMDcxoGgA/873//05AhQzR16lSXzVyitoGChN4NWBe9G7AmejdgXfRuoODK0cVGL1y4oB49eujChQtasmSJSpQo4a28AOSxefPm6fbbb9eUKVM0bNgwX6cDwEPo3YB10bsBa6J3A9ZF7wYKthwdkd6rVy9t2bJF7777rv744w/H+muuuUbXXHONx5MDkDeWLl2qQYMGqWfPnipZsqTmz5/vuK1Pnz4+ywuA++jdgDXRuwHroncD1kTvBgq+HA2klyxZUh07dtSPP/7otH7QoEE0dKAAi4mJUY8ePZSamqpZs2Y53UZDBwo2ejdgTfRuwLro3YA10buBgi9XFxt1R1xcnOx2Oxc9ATzM1xcbycvalqhv/Lf8l+qb2sZ/yX+ptiXqG/8t/6X6prbxX+PL+s7r3p2XeC+Br2WnrnI0RzoAAAAAAAAAAP81DKQDAAAAAAAAAOACA+kAAAAAAAAAALjAQDoAAAAAAAAAAC4wkA4AAAAAAAAAgAsMpAMAAAAAAAAA4AID6QAAAAAAAAAAuODv6wS8yRiTZ7FsNluexQJAfQNWRW0D1kV9A9aUl7Wdl3gfAfJWXr2XUNtwB0ekAwAAAAAAAADgAgPpAAAAAAAAAAC4wEA6AAAAAAAAAAAuMJAOAAAAAAAAAIALDKQDAAAAAAAAAOACA+kAAAAAAAAAALjAQDoAAAAAAAAAAC745/QO//zzjxYuXKgzZ86oVq1aGjhwoIoVK+aN3ADkoeTkZH3//ffasmWLQkJC1LlzZzVr1szXaQHwAHo3YE30bsC66N2ANdG7gYItR0ekv/HGG7rjjjt09uxZ2e12TZ06VXXq1NGxY8e8lR+APHD27Fk1atRIP/74o8LCwhQVFaW2bdvqhRde8HVqANxE7wasid4NWBe9G7AmejdgASYHdu/e7fR3QkKCKVmypHn99dezvY/Y2FgjycTGxuYkdL4niYXFp4s7NXX+/HkTFRXltO7NN980wcHBJjU1NVv7sGptG0N9s/h+caeu6N1Z8/X/KwsLvdt7fP1/y8JC7y74fP0aYsm/S27rit6dP/j69cOSf5fs1FWOpnapWbOm09+XLl1ScnKySpYsmZPdAMhnihUrluFU0djYWIWFhcnPj0spAAUZvRuwJno3YF30bsCa6N1AwZfjOdKPHz+uV155RRcuXNDq1as1evRo3XnnnVlun5iYqMTERMffcXFxucsUgNfNmDFDmzZt0v79+3Xu3DnNmzcvy22pbaDgoHcD1kXvBqyJ3g1YF70bKLhy/JNXYGCgrrnmGlWvXl0hISFaunSpTp8+neX2EydOlN1udyzh4eFuJQzAeypWrKjatWsrIiJC+/fv17p167LcltoGCg56N2Bd9G7AmujdgHXRu4GCy/b/5wfKleTkZDVq1Eg33HCDPv7440y3yezXs/DwcMXGxio0NDS3ofMdm83m6xTwH+fpmvrqq680dOhQHThwINNm/V+pbYn6hu95sq7o3f+H2oav0bu9h/qGr9G7Cz7eR5AVT9YVvTvvUdvISnbqKsdTu1wuICBAjRo10u7du7PcJjAwUIGBge6EAeADzZo1U2pqqvbv359pQ6e2gYKJ3g1YF70bsCZ6N2Bd9G6gYMnR1C6LFy92+vv06dNatmyZmjRp4tGkAOStLVu26MyZM07rvvnmGwUFBalu3bo+ygqAJ9C7AWuidwPWRe8GrIneDRR8OToi/csvv9TYsWPVoEEDXbp0SZGRkWrRooWee+45b+UHIA9cuHBBbdu2VUREhMqVK6cdO3Zoz549mjFjhkqVKuXr9AC4gd4NWBO9G7AuejdgTfRuoODL8Rzpu3fv1vr16xUQEKC6deuqXr16OQoYFxcnu91uufmcmGMJvuZuTV26dEkrVqzQ0aNHVbFiRbVu3VrFihXL9v2tWtsS9Q3fc7eu6N2Zo7bha/Ru76G+4Wv07oKP9xFkxZ26onf7HrWNrHhljvRatWqpVq1auU4KQP5UpEgRdenSxddpAPACejdgTfRuwLro3YA10buBgi1Hc6QDAAAAAAAAAPBfw0A6AAAAAAAAAAAuMJAOAAAAAAAAAIALDKQDAAAAAAAAAOACA+kAAAAAAAAAALjAQDoAAAAAAAAAAC4wkA4AAAAAAAAAgAv+vk7AKowxeRbLZrPlWSwAeVff1DaQt+jdgHXRuwG4i88JgDVR23AHR6QDAAAAAAAAAOACA+kAAAAAAAAAALjAQDoAAAAAAAAAAC4wkA4AAAAAAAAAgAsMpAMAAAAAAAAA4AID6QAAAAAAAAAAuMBAOgAAAAAAAAAALuR6IH3u3LmqUaOGHnzwQU/mA8DHkpKS1KVLF9WoUUN79uzxdToAPIjeDVgTvRuwLno3YE30bqBgytVA+qFDh/Twww8rKChIx48f93ROAHxo7NixiouL0759+5SYmOjrdAB4CL0bsC56N2BN9G7AuujdQMGU44H0lJQUDR48WBMmTFBERIQXUgLgKz/99JN++eUXTZo0ydepAPAgejdgXfRuwJro3YB10buBgivHA+nPPfecypcvr+HDh3sjHwA+EhUVpXvvvVezZ89W0aJFfZ0OAA+idwPWRO8GrIveDVgTvRso2PxzsvGSJUv05Zdf6s8//8z2fRITE51OU4mLi8tJSAB5IC0tTbfffrvGjBmjJk2aaOPGjVe9D7UNFAz0bsCa6N2AddG7AWuidwMFX7aPSD916pSGDh2qmTNnqmTJktkOMHHiRNntdscSHh6eq0QBeM8rr7yitLQ0Pfnkk9m+D7UN5H/0bsC66N2ANdG7AeuidwMFn80YY7Kz4fz58zVgwABVqVLFse748eOy2WwqV66cli9frooVK2a4X2a/noWHhys2NlahoaEeeAj/PTabzdcpIB9yp6YaNWqk48ePKyQkRNK/dXv06FGFh4dr8ODBev311zPch9r2PGobWcltXdG78w/qG5mhdxd81DayQu9GTvBeUrDktq7o3f891HbBkp26yvbULp07d9bff//ttO6+++5TYGCg3n//fZUtWzbT+wUGBiowMDC7YQD4wIIFC5ya819//aW+ffvq448/VpMmTTK9D7UN5H/0bsC66N2ANdG7AeuidwMFX7YH0oODg1WjRg2ndUWLFlVQUFCG9QAKlitPD4uJiXGsL1OmjA8yAuAJ9G7AuujdgDXRuwHroncDBV+250gHAAAAAAAAAOC/KNtHpGdm2rRpzPcDWFD9+vW1Z88eVa5c2depAPAwejdgTfRuwLro3YA10buBgsetgfTy5ct7Kg8A+UhgYCCnjgIWRe8GrIneDVgXvRuwJno3UPAwtQsAAAAAAAAAAC4wkA4AAAAAAAAAgAsMpAMAAAAAAAAA4AID6QAAAAAAAAAAuMBAOgAAAAAAAAAALjCQDgAAAAAAAACAC/6+TgA5Z4zJs1g2my3PYgH/ddQ2YF3UN2BNeVnbeYn3ESBv5dV7CbUN5C2+A1gPR6QDAAAAAAAAAOACA+kAAAAAAAAAALjAQDoAAAAAAAAAAC4wkA4AAAAAAAAAgAsMpAMAAAAAAAAA4AID6QAAAAAAAAAAuMBAOgAAAAAAAAAALvjnZOPo6GidPXvWaV2hQoVUtWpVjyYFIO/t27dPxhindaVLl5bdbvdRRgA8gd4NWBe9G7AmejdgXfRuoGDL0UD6hx9+qJdfflmVKlVyrAsNDdXmzZs9nhiAvFW7dm2VKVNGRYsWdax79tlnddddd/kuKQBuo3cD1kXvBqyJ3g1YF70bKNhyNJAuSY0aNdLatWu9kQsAH5s5c6a6devm6zQAeBi9G7AuejdgTfRuwLro3UDBlas50k+ePKnY2FhP5wLAxxISEnT06FGlpaX5OhUAHkbvBqyJ3g1YF70bsCZ6N1Bw5Xggfd26dapbt67KlSun2rVr69dff/VGXgB8YMCAAWrSpImCg4N13333KS4uztcpAfAAejdgXfRuwJro3YB10buBgitHA+nXXnut1q9frzNnziguLk79+vVT7969tX379izvk5iYqLi4OKcFQP4zYcIERUdH6+TJk9qwYYMiIyM1cuTILLentoGCgd4NWBe9G7AmejdgXfRuoIAzbkhLSzPVqlUzTz75ZJbbjB8/3kjKsMTGxroTGnkks/87lvy5eLqmZs2aZWw2m4mPj8/0dmq7YPP165XFd/VN77Y+X79eWXxT28bQu+E5vq4NKyz0buRHvq4LqyyerCt6NzzF13VhhSU7dZWrOdLT2Ww2Va5cWYcOHcpym6efflqxsbGO5ciRI+6EBJBHIiIiZIzJsmapbaBgoncD1kXvBqyJ3g1YF70bKFj8c7JxSkqK/P3/7y5xcXHatm2bRo0aleV9AgMDFRgYmPsMAXjdlbUtSStXrlRgYKDCw8MzvQ+1DRQM9G7AmujdgHXRuwFroncDBV+OBtLbt2+vu+++W40bN9aZM2f00ksvKSAgQKNHj/ZWfgDywKxZs7Rx40bdcsstKl26tJYsWaJXXnlFTz75pIKDg32dHgA30LsBa6J3A9ZF7wasid4NFHw5Gkj/8ssvNWnSJE2ZMkVFixZV8+bNNXfuXJUuXdpb+QHIA8OHD1dQUJAmT56s48ePq2rVqvr666/Vu3dvX6cGwE30bsCa6N2AddG7AWuidwMFn+3/T0ifZ+Li4mS32xUbG6vQ0NC8DI1csNlsvk4B2eTrmqK2CxZqu2DxdV1R3wUL9V1w+LqmqG1khfcR9/m6rqhvZIba9gxf1hW1jaxQ3+7LTl25dbFRAAAAAAAAAACsjoF0AAAAAAAAAABcYCAdAAAAAAAAAAAXGEgHAAAAAAAAAMAFBtIBAAAAAAAAAHCBgXQAAAAAAAAAAFxgIB0AAAAAAAAAABf8fZ0A8jdjTJ7EsdlseRIHwL/yqrYl6hvIa3lZ33mF9xEgb/E5AbAmahuwLuo7b3BEOgAAAAAAAAAALjCQDgAAAAAAAACACwykAwAAAAAAAADgAgPpAAAAAAAAAAC4wEA6AAAAAAAAAAAuMJAOAAAAAAAAAIALDKQDAAAAAAAAAOCCf27udPLkSS1btkw2m02dO3dWiRIlPJ0XAB9IS0vTqlWrtHv3btWrV08tWrTwdUoAPITeDVgTvRuwLno3YE30bqDgyvER6VOnTlX16tX15Zdf6tdff1WbNm20ZcsWb+QGIA+dPHlS119/vYYMGaLVq1dr3Lhxuueee3ydFgAPoHcD1kTvBqyL3g1YE70bKNhydET68uXLNXr0aM2dO1d9+/aVJJ09e1anTp3ySnIA8s7gwYNls9n0999/q0iRIpKklStX+jgrAO6idwPWRe8GrIneDVgXvRso2HI0kP7WW2/phhtucDRzSSpZsqRKlizp8cQA5J3Nmzdr2bJl+u233xzNXJJuuOEGH2YFwBPo3YA10bsB66J3A9ZE7wYKvhwNpK9atUpjxozRrl27tHLlSpUpU0Y33nij7HZ7lvdJTExUYmKi4++4uLjcZwvAK1atWqXChQurWbNmmj9/vmJjY9WgQQM1btw4y/tQ20DBQO8GrIneDVgXvRuwJno3UPBle450Y4zOnTun5cuXq2fPnlq5cqVefvll1axZU2vXrs3yfhMnTpTdbncs4eHhHkkcgOecPXtWwcHB6ty5s2bMmKHIyEi1bdtWQ4cOzfI+1DaQ/9G7AeuidwPWRO8GrIveDRR8NmOMye7GRYsWVenSpbVz504FBwdLkvr376/9+/dr8+bNmd4ns1/PwsPDFRsbq9DQUDfTh1XYbDZfp1DguVNTkyZN0tixY/X+++/rgQcekCRt27ZNjRo10vz589WrV68M96G2kV3Ut/vcqSt6N6zAqu8j9G6A+s4MvRtWYNXalnJf3/RuWIVV6zs7dZXtI9IlqVatWmrdurWjmUtS586d9ddffymr8fjAwECFhoY6LQDyl1q1akmSunTp4ljXoEEDlS1bVtu3b8/0PtQ2UDDQuwFroncD1kXvBqyJ3g0UfDkaSO/bt6+2bdum1NRUx7otW7aoWrVqlv01Avgv6Nixo0JCQpyOcDl69KhOnz6t6tWr+zAzAO6idwPWRO8GrIveDVgTvRso+HJ0sdFHH31U8+bNU6dOndS9e3ft2LFD3333nebOneut/ADkgZCQEL333nsaOXKkNm/erNDQUM2aNUs33nij+vfv7+v0ALiB3g1YE70bsC56N2BN9G6g4MvREekhISFau3atbrvtNp06dUoNGjTQjh071L17d2/lByCP3HXXXVq6dKmCgoJ06dIlvfbaa1qyZIn8/XP0exuAfIbeDVgXvRuwJno3YF30bqBgy9HFRj0hLi5OdrudCyPACacous/XNUVtIyvUt/t8XVfUN3zNqu8jvq4pahv5AfXtHdQ3fM2qtS35tr6pbeQHVq1vj19sFAAAAAAAAACA/xoG0gEAAAAAAAAAcIGBdAAAAAAAAAAAXGAgHQAAAAAAAAAAFxhIBwAAAAAAAADABQbSAQAAAAAAAABwgYF0AAAAAAAAAABc8Pd1AoAkGWPyLJbNZsuzWADytr7zCu8jQN7icwJgXXlV39Q2kLfo3YB1/Zd7N0ekAwAAAAAAAADgAgPpAAAAAAAAAAC4wEA6AAAAAAAAAAAuMJAOAAAAAAAAAIALDKQDAAAAAAAAAOACA+kAAAAAAAAAALjAQDoAAAAAAAAAAC7452TjUaNG6ciRIxnWN23aVBMmTPBUTgDy2IoVK/T6669netvkyZNVq1atPM4IgKfQuwFroncD1kXvBqyJ3g0UfDkaSL/11lsVHx/v+DsmJkZDhgzRDTfc4PHEAOSd6tWra+TIkU7rpkyZojVr1ig8PNxHWQHwBHo3YE30bsC66N2ANdG7gYIvRwPp7dq1c/r7ww8/lL+/v+666y4PpgQgr1WoUEEVKlRw/J2SkqJ7771Xd9xxh4oUKeLDzAC4i94NWBO9G7AuejdgTfRuoOBza470GTNmqGfPnipXrpyn8gGQD/z88886ceKE7r33Xl+nAsDD6N2ANdG7AeuidwPWRO8GCp4cHZF+uS1btmjLli165ZVXXG6XmJioxMREx99xcXG5DQkgj0yfPl3NmzdXgwYNstyG2gYKHno3YF30bsCa6N2AddG7gYIn10ekT58+XZUrV1bXrl1dbjdx4kTZ7XbHwrxPQP52/Phx/fLLL1f9VZzaBgoeejdgTfRuwLro3YA10buBgilXA+kJCQn66quvNGzYMPn5ud7F008/rdjYWMeS2dXHAeQfn332mYKCgjRo0CCX21HbQMFC7wasi94NWBO9G7AuejdQMOVqape5c+cqLi5Ow4YNu+q2gYGBCgwMzE0YAD4wc+ZMDR48WMWKFXO5HbUNFCz0bsC66N2ANdG7AeuidwMFU64G0mfMmKFu3bpxSglgMX/88Yf27Nmj2bNn+zoVAB5G7wasid4NWBe9G7AmejdQcOV4IH3//v36/fffNW/ePG/kA8CHZsyYoYYNG6pZs2a+TgWAB9G7AeuidwPWRO8GrIveDRRcOR5IL1y4sH788Ud169bNG/kA8KE777xT5cqV83UaADyM3g1YF70bsCZ6N2Bd9G6g4MrxQHqlSpVUqVIlb+QCwMc6dOjg6xQAeAG9G7AuejdgTfRuwLro3UDB5frS3wAAAAAAAAAA/McxkA4AAAAAAAAAgAsMpAMAAAAAAAAA4AID6QAAAAAAAAAAuMBAOgAAAAAAAAAALvjndUBjjCQpLi4ur0MDlpZeW76OT20Dnkd9A9ZEbQPWRX0D1uXL+qa2Ae/JTm3n+UD6+fPnJUnh4eF5HRqwtPPnz8tut/s0vkRtA95AfQPWRG0D1kV9A9bly/qmtgHvyU5t20we/5SWlpamqKgohYSEyGazZes+cXFxCg8P15EjRxQaGurV/KwYy4qPKS9j5ffHZIzR+fPnVaFCBfn5+W62ptzUtsT/Y0GJZcXHVBBiFeT6zu/PbX6PZcXHlJex8vtjKsi1LfH/WFBiWfExFYRYBbm+8/tzm99jWfEx5WWsgvCY8kN907utHcuKj6kgxMpJbef5Eel+fn6qVKlSru4bGhrq9SfcyrGs+JjyMlZ+fky+PNolnTu1LfH/WFBiWfEx5fdYBb2+8/NzWxBiWfEx5WWs/PyYCnptS/w/FpRYVnxM+T1WQa/v/PzcFoRYVnxMeRkrvz8mX9c3vfu/EcuKjym/x8pubXOxUQAAAAAAAAAAXGAgHQAAAAAAAAAAFwrEQHpgYKDGjx+vwMBAYuXjOFaNZcXHlJ/w/1gwYlnxMVk5Vn5g1eeWmisYsaz4mPIT/h8LRiwrPiYrx8oPrPrcUnMFI5YVH1N+wv9jwYhlxcdktVh5frFRAAAAAAAAAAAKkgJxRDoAAAAAAAAAAL7CQDoAAAAAAAAAAC4wkA4AAAAAAAAAgAv+vk7gauLj47Vy5UqlpqaqdevWKl68uNdixcTE6LffflO5cuXUunVrr8WJj4/Xpk2blJiYqAYNGqhcuXJeixUTE6ONGzfKGKP69et7NVa65cuX6/jx47r55ptVrFgxj+47KipKf/zxR4b1PXv2VHBwsEdjpYuNjdXatWsVGBio66+/XoULF/bo/qOjo7Vo0aJMb2vVqpWqVKni0Xj5yfbt27V3715FRESocePGXo21bt06HThwQN27d5fdbvdanN27d2vXrl0qX768rrvuOvn5ee/3yh07dmjfvn2qUKGC12NJ0tmzZxUZGamIiAi1bNnS4/tfvHixzp0757SuatWqatGihcdjpdu2bZsOHjyohg0beqXWVqxYoWPHjmVYHxISoptuusnj8fILerd76N3uo3d7D73bPfRu99G7vYPe7R56t/vo3d5D73YPvdt9Vurd+XogfdWqVerdu7eqVKmigIAA/fPPP/r222/VpUsXj8aJjY3VY489pp9//lmpqalq06aN1xr6yy+/rKlTp6p69eoqXLiwVq1apbFjx2r8+PEej/Xiiy9q+vTpqlOnjhISErR+/Xo98cQTeuGFFzweK93q1avVo0cPXbx4UXv27FGNGjU8uv/Nmzfr9ttv14ABA5zWd+zY0SsNfdq0aXrsscfUsGFDhYSE6MSJE/r+++89WvjR0dGaP3++07oDBw5o/fr1Wr58uSUbempqqu6880799NNPat68uTZt2qTWrVvru+++8/gHpvnz52vChAm6cOGC9u7dq+3bt3ulof/999+69957dfr0adWqVUt//fWXihYtqh9++MHjdbBz504NGzZMCQkJqlKlirZs2aKiRYvqp59+UvXq1T0aK50xRkOHDtWSJUs0YMAArzT0cePGKTk5Wddee61jXdu2bb3S0E+ePKn+/ftrz549atmypfbu3avbb79dTz/9tEfjrFmzRps3b3Zat2DBArVo0cKyX8bp3e6hd7uP3u0d9G730LvdR+/2Hnq3e+jd7qN3ewe92z30bvdZsnebfCo5OdlUqVLFjBgxwrHuscceM2XKlDEXLlzwaKyoqCgzbdo0Ex8fb3r37m1uueUWj+7/ctOnTzdxcXGOvxcvXmwkmeXLl3s81tdff20SExMdf8+bN89IMjt37vR4LGOMiY6ONtWqVTMTJ040ksyePXs8HmPBggUmMDDQ4/vNzKJFi4zNZjPz5893rNu1a5f5+++/vR779ttvNzVr1jRpaWlej+UL06ZNM8WKFTO7d+82xhhz6NAhU6JECfPmm296PNacOXPMli1bzJYtW4wks337do/HMMaYdevWmVWrVjn+TkpKMu3atTPt27f3eKwNGzaYv/76yylWkyZNzMCBAz0eK91bb71lunbtajp06GBuv/12r8Ro0qSJmThxolf2faXWrVubG264wcTHxxtjjElNTTULFizwetw9e/YYSebLL7/0eixfoHe7j97tHnq399C73UPvdh+92zvo3e6jd7uH3u099G730LvdZ8XenW8H0pctW2YkmV27djnWHT16NMMbrKd5u6FnJiQkxLz77rtej7Njxw4jyaxbt84r++/Xr5959tlnzYoVK7za0AsXLmwWL15sfv75Z3Po0CGPx0jXtm1b06NHD6/tPyvR0dGmSJEi5vXXX8/z2Hmlbdu2ZvDgwU7r7rvvPtOoUSOvxfR2Q8/MW2+9ZcLCwvIk1oABA0z37t29su+NGzeaChUqmOPHj5uOHTt6taGPGjXKzJs3z6xbt85cunTJK3HS+8v69eu9sn9XnnrqKRMWFua1x+Zr9G7Po3fnDL3be+jdnkfvzj56t/fQuz2P3p0z9G7voXd7Hr07+6zau/Pt1C7bt29X4cKFVatWLce6ihUrqkSJEtq+fbt69+7tw+w8Z/369Tp//rzq1avnlf3v2bNHmzZt0pkzZzRt2jSNGjVKzZs393icKVOm6MiRI/rmm2+0du1aj+//cjabTS+//LL8/Py0evVqDRkyRB9//LEKFSrksRgpKSlas2aNJk2apMOHD2vLli2Oubf8/b1bNl999ZVSUlJ05513ejWOL23fvj3DqaL169fXp59+KmOMbDabjzLzrCVLlnittiXp+++/1/nz57Vx40Zt2LBBc+fO9XiM8+fPa9CgQZoyZUqezPX4888/6/Dhw9qxY4fS0tL05Zdfqk2bNh6N8ccffygsLEyNGjXS77//rkuXLqlBgwaqWLGiR+NcKTU1VZ999pmGDBmioKAgr8byFXq3Z9C7c4fe7V30bs+gd+cOvdt76N2eQe/OHXq3d9G7PYPenTtW7d35diA9NjZWYWFhGdaXLFlSMTExeZ+QF8TExOjOO+9Uly5d1KFDB6/EOHDggObPn69jx47p3LlzXrmwxPbt2zV+/HitXr3a682uevXq+ueffxQRESFJ2rp1q66//nrVrl1bTzzxhMfixMTEKCkpScuXL9fbb7+tBg0aaPv27SpSpIgWLlyoatWqeSzWlWbMmKFevXqpbNmyXovha7GxsSpRooTTupIlSyo5OVkXLlzw+MVyfOHjjz/W4sWLtWzZMq/F+OWXX3Ty5Elt2rRJTZs29cprZuTIkerQoUOefIl6+eWX1bVrV9lsNqWkpOiee+7RwIEDtXv3boWEhHgszqlTpxQaGqobb7xRhQsXVuHChbVy5UqNGzdOzz33nMfiXOnnn3/W8ePHdc8993gthq/Ruz2D3p079G7vond7Br07d+jd3kPv9gx6d+7Qu72L3u0Z9O7csWzv9vgx7h4yadIkExwcnGF9hQoVzDPPPOO1uHl1itn58+fN9ddfb6677joTHR3t9XjGGLNkyRIjySxdutSj+23durXp16+f+frrr83XX39tJkyYYCSZd955x2zYsMGjsTIzdOhQ06ZNG4/u8/z580aSqVmzpjl//rwxxpjExETTokULc/PNN3s01uXST4P69ddfvRYjPyhatKh54403nNZ98sknRpJJTk72Ssy8PMXsm2++MQEBAWbWrFlej2XMv6/NTp06mRtvvNGj+42MjDRBQUFmxowZjvquV6+eueGGG8zXX3/t9VOc9+/fbySZ3377zaP7feyxx4wk8/HHHzvW/fDDD0aS2bhxo0djXa53796mZcuWXtt/fkDv9jx6d/bRu72L3u1Z9O6coXd7D73b8+jd2Ufv9i56t2fRu3PGqr3bzzvD8+6rXr26Lly4oHPnzjnWXbp0SadPn/bqr5J5IT4+Xt27d9elS5cUGRmp4sWL50ncjh07qkKFClqxYoVH99uiRQsFBARo/vz5mj9/vmP/kZGR2rp1q0djZSY0NFSnT5/26D6LFSumsmXLqmvXro5faQsXLqyePXtq06ZNHo11uRkzZqhKlSrq3Lmz12LkB9WrV9fhw4ed1h06dEjh4eFeP7rC2/73v/9pyJAhmjp1ap6dJli4cGENHDhQq1atUlpamsf2GxwcrN69e2vx4sWO+j59+rSOHDmi+fPnKzEx0WOxMhMaGipJHq/v9Cus33LLLY51PXv2VEBAgNfq++TJk/rpp5907733emX/+QW92/Po3dlH7/Yuerdn0btzht7tPfRuz6N3Zx+927vo3Z5F784Zy/Zurw3Ru+ncuXMmKCjITJ061bHuyy+/NP7+/ubo0aNei+vtX8bj4+NNmzZtTOPGjc3Zs2e9FicuLs5xVdx0R44cMf7+/mbmzJlei2uM8epFT6Kiopz+vnTpkqlevboZMmSIx2MNHz7cdOrUyWnd7bffbpo3b+7xWMYYk5CQYMLCwsyLL77olf3nJ2PHjjVVqlQxCQkJxhhjkpOTzbXXXmtGjhzptZh58cv43LlzTeHChc306dO9FsOYjHVgjDEjRowwlStX9mpcY4zXLnoSHR2d4Zf2Tz75xNhsNrN3716Pxjp48KApVKiQWblypWPd3r17jSTz888/ezRWutdff92EhIRkeF+2Gnq3e+jd7qN3ew+92z30bvfQu72H3u0eerf76N3eQ+92D73bPVbt3fn2J6iwsDC9+OKLevTRR3Xq1CkFBARo4sSJGjt2rFcmpp87d66Sk5N17NgxFSpUSHPmzFFgYKD69u3r0Th9+/bVxo0b9eabb2rx4sWO9XXq1FGDBg08FufMmTPq2bOn+vTpo+rVq+v48eOaNm2amjVrpsGDB3ssTl575plndPHiRbVp00YpKSn69NNPlZycrBdeeMHjsSZMmKBmzZpp6NChateunTZt2qTvvvtOCxYs8HgsSZo3b57i4uI0bNgwr+w/P3n88cc1Z84c9ejRQ/3799eCBQt07tw5PfPMMx6PtWvXLm3ZskWHDh2SJP3666/666+/1KRJE9WsWdNjcZYvX65Bgwapa9euCg4O1pw5cxy3DRw4UH5+njsBaOLEiTp+/LhuvPFGBQYG6o8//tB3332n2bNneyxGXjt27JhuvfVW9evXTxEREdqyZYumTZumcePGOX7J9pQqVaroqaee0uDBg/Xoo48qICBA77//vtq3b5/hYjyeMnPmTN12220KDg72yv7zC3q3e+jd7qN3ew+92z30bvfQu72H3u0eerf76N3eQ+92D73bPVbt3TZjjPHa3j3g559/1g8//KC0tDT16NHD4w023V133aWEhASndaGhofrkk088Gmf48OG6cOFChvW9e/f2eKM9c+aMZs2apb///lthYWG6/vrr1adPH4++sWRm165dGj9+vN59912vXIThxx9/1JIlS5Samqp69erprrvuUpEiRTweR/r3tJCpU6fqwIEDqlSpkoYMGaLatWt7JdY777yj6Ohor3w4yY/OnDmjKVOmaN++fapSpYpGjRql8uXLezzOwoUL9eWXX2ZYP2zYMI++ef/yyy/67LPPMr3tiy++UEBAgMdiSf9+MFm8eLHi4+NVtWpV3XHHHQoPD/dojMy88sorKl68uO6//36P7/v48eP67LPPtGfPHlWoUEG9evVSs2bNPB4n3YIFC7Rw4UL5+fmpRYsWuuOOO7xyiuPRo0f1+OOP65lnnlH9+vU9vv/8iN6de/Ru99G7vYfe7R56t/vo3d5D7849erf76N3eQ+92D73bfVbr3fl+IB0AAAAAAAAAAF/KtxcbBQAAAAAAAAAgP2AgHQAAAAAAAAAAFxhIBwAAAAAAAADABQbSAQAAAAAAAABwgYF0AAAAAAAAAABcYCAdAAAAAAAAAAAXGEgHAAAAAAAAAMAFBtIBAAAAAAAAAHCBgXQAAAAAAAAAAFxgIB0AAAAAAAAAABcYSAcAAAAAAAAAwAUG0gEAAAAAAAAAcOH/AaRwjZ9kauHWAAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 1500x270 with 6 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Converged after 20 Wilkinson-shifted QR steps (plus deflations).\n"
          ]
        }
      ],
      "source": [
        "# Finish the QR iteration if the repeatable cell has not yet converged.\n",
        "for _ in range(100):\n",
        "    if m == 1:\n",
        "        break\n",
        "    old_m = m\n",
        "    m, action = advance(H, Z, m)\n",
        "    qr_steps += (m == old_m)\n",
        "    qr_frames.append((H.copy(), f\"{action}\\nActive size {m}\"))\n",
        "else:\n",
        "    raise RuntimeError(\"QR iteration did not converge in 100 further advances.\")\n",
        "\n",
        "show_fill([qr_frames[i] for i in np.linspace(0, len(qr_frames)-1, 6, dtype=int)])\n",
        "print(f\"Converged after {qr_steps} Wilkinson-shifted QR steps (plus deflations).\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a90fe336",
      "metadata": {},
      "source": [
        "## 4. Read the complex Schur form\n",
        "\n",
        "The result $A=ZHZ^*$ has a **triangular** $H$; its diagonal entries are the eigenvalues, including the nonreal conjugate pair. Complex Schur vectors are generally not eigenvectors, but their leading groups span invariant subspaces."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "id": "2de1bbcd",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Complex Schur form:\n",
            " [[ 8.   -0.j    -0.828+0.337j -0.132+0.124j -0.14 -0.036j  0.01 +0.179j -0.545-0.416j\n",
            "  -0.043-0.033j  0.158+0.12j ]\n",
            " [ 0.   +0.j     6.   +0.j    -0.105+0.04j  -0.063-0.047j -0.036+0.106j -0.166-0.281j\n",
            "  -0.048-0.081j  0.156+0.265j]\n",
            " [ 0.   +0.j     0.   +0.j     3.   -2.j     0.032+0.038j  0.   -0.001j -0.112-0.36j\n",
            "  -0.324-0.306j -0.02 -0.244j]\n",
            " [ 0.   +0.j     0.   +0.j     0.   +0.j     4.   +0.j     0.024-0.043j -0.091-0.033j\n",
            "   0.238+0.113j  0.486+0.225j]\n",
            " [ 0.   +0.j     0.   +0.j     0.   +0.j     0.   +0.j     3.   +2.j     0.284-0.247j\n",
            "   0.442-0.086j  0.157-0.192j]\n",
            " [ 0.   +0.j     0.   +0.j     0.   +0.j     0.   +0.j     0.   +0.j    -1.   +0.j\n",
            "  -0.24 +0.j     0.531+0.j   ]\n",
            " [ 0.   +0.j     0.   +0.j     0.   +0.j     0.   +0.j     0.   +0.j     0.   +0.j\n",
            "   0.5  +0.j     0.214+0.j   ]\n",
            " [ 0.   +0.j     0.   +0.j     0.   +0.j     0.   +0.j     0.   +0.j     0.   +0.j\n",
            "   0.   +0.j     2.   +0.j   ]]\n",
            "Eigenvalues from Schur form:  [-1. +0.j  0.5+0.j  2. +0.j  3. +2.j  3. -2.j  4. +0.j  6. +0.j  8. -0.j]\n",
            "Eigenvalues of A:           [-1. +0.j  0.5+0.j  2. +0.j  3. -2.j  3. +2.j  4. +0.j  6. +0.j  8. +0.j]\n",
            "Largest subdiagonal entry: 0.0\n",
            "Relative reconstruction error: 6.149743752757245e-15\n",
            "Unitarity error: 4.610546013904431e-15\n"
          ]
        }
      ],
      "source": [
        "print(\"Complex Schur form:\\n\", H)\n",
        "print(\"Eigenvalues from Schur form: \", np.sort_complex(np.diag(H)))\n",
        "print(\"Eigenvalues of A:          \", np.sort_complex(la.eigvals(A)))\n",
        "print(\"Largest subdiagonal entry:\", np.max(abs(np.diag(H, -1))))\n",
        "print(\"Relative reconstruction error:\", la.norm(A - Z @ H @ Z.conj().T) / la.norm(A))\n",
        "print(\"Unitarity error:\", la.norm(Z.conj().T @ Z - np.eye(n)))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4f965e1e",
      "metadata": {},
      "source": [
        "Observe the two nonreal eigenvalues on the diagonal; no $2\\times2$ real Schur block is required. The matrix may remain dense **above** the diagonal, but all entries below it converge to zero (to the chosen tolerance)."
      ]
    }
  ],
  "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.13.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}