diff --git a/.ipynb_checkpoints/Apollonian Circle Packings-checkpoint.ipynb b/.ipynb_checkpoints/Apollonian Circle Packings-checkpoint.ipynb index 791e551..40bfa46 100644 --- a/.ipynb_checkpoints/Apollonian Circle Packings-checkpoint.ipynb +++ b/.ipynb_checkpoints/Apollonian Circle Packings-checkpoint.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Calculations to find the quadratic form of the octahedron packing" + "# Calculations to find the quadratic form of the octahedral packing" ] }, { @@ -55,7 +55,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "For an Apollonian packing, you can make a matrix, $W$, where the rows are the coordinates of each of the circles in a quadruple. In an octahedral packing, this is impossible, as there are six circles in a unit. You can try creating a $6\\times 4$ matrix, but, later, we end up needing to invert $WPW^T$, the product of a $6\\times 4$, $4\\times 4$, and $4\\times 6$ matrix, which is singular. Fortunately, the sextuples come in three pairs of circles. Each pair consists of circles that aren't tangent to each other. For whatever reason (we aren't sure yet why, but the Guettler and Mallows says this is the case), the average of the coordinates of the circles in a pair is the same across the sextuple. So, we can make a matrix where the first three rows are the coordinates of circles from different pairs, i.e. three mutually tangent circles, and the fourth is the average of the coordinates in a pair. From this, you can recover the coordinates for all the circles in a sextuple.\n", + "For an Apollonian packing, you can make a matrix, $W$, where the rows are the coordinates of each of the circles in a quadruple. In an octahedral packing, this is impossible, as there are six circles in a unit. You can try creating a $6\\times 4$ matrix, but, later, we end up needing to invert $WPW^T$, the product of a $6\\times 4$, $4\\times 4$, and $4\\times 6$ matrix, which is singular. Fortunately, the sextuples come in three pairs of circles. Each pair consists of circles that aren't tangent to each other. It turns out that the average of the coordinates of the circles in a pair is the same across the sextuple. So, we can make a matrix where the first three rows are the coordinates of circles from different pairs, i.e. three mutually tangent circles, and the fourth is the average of the coordinates in a pair. From this, you can recover the coordinates for all the circles in a sextuple.\n", "\n", "Here $W$ is such a matrix computed for the $(0, 0, 1, 1, 2, 2)$ root sextuple." ] @@ -70,7 +70,7 @@ " [2, 0, 0, 1],\n", " [2, 0, 0, -1],\n", " [-1, 1, 0, 0],\n", - " [3, 1, sqrt(2), 0]\n", + " [6, 2, 2*sqrt(2), 0]\n", "])" ] }, @@ -107,10 +107,10 @@ { "data": { "text/plain": [ - "[ 1 -1 -1 -1]\n", - "[-1 1 -1 -1]\n", - "[-1 -1 1 -1]\n", - "[-1 -1 -1 -1]" + "[ 1 -1 -1 -2]\n", + "[-1 1 -1 -2]\n", + "[-1 -1 1 -2]\n", + "[-2 -2 -2 -4]" ] }, "execution_count": 4, @@ -123,6 +123,21 @@ "M" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$M^{-1}$ is the matrix of the quadratic form, although it is nice to normalize it to have ones along the diagonal, which is fine since it is equal to zero.\n", + "$$\n", + "\\left(\\begin{matrix}\n", + " 1 & 0 & 0 & -1/2\\\\\n", + " 0 & 1 & 0 & -1/2\\\\\n", + " 0 & 0 & 1 & -1/2\\\\\n", + " -1/2 & -1/2 & -1/2 & 1/4\n", + "\\end{matrix}\\right)\n", + "$$" + ] + }, { "cell_type": "code", "execution_count": 5, @@ -133,10 +148,10 @@ { "data": { "text/plain": [ - "[ 1/2 0 0 -1/2]\n", - "[ 0 1/2 0 -1/2]\n", - "[ 0 0 1/2 -1/2]\n", - "[-1/2 -1/2 -1/2 1/2]" + "[ 1/2 0 0 -1/4]\n", + "[ 0 1/2 0 -1/4]\n", + "[ 0 0 1/2 -1/4]\n", + "[-1/4 -1/4 -1/4 1/8]" ] }, "execution_count": 5, @@ -220,7 +235,7 @@ { "data": { "text/plain": [ - "[b_avg == b1 + b2 + b3 - sqrt(2*b1*b2 + 2*(b1 + b2)*b3), b_avg == b1 + b2 + b3 + sqrt(2*b1*b2 + 2*(b1 + b2)*b3)]" + "b1^2 + b2^2 + b3^2 - b1*b_avg - b2*b_avg - b3*b_avg + 1/4*b_avg^2" ] }, "execution_count": 9, @@ -229,28 +244,7 @@ } ], "source": [ - "quad = 2 * factor(simplify(D[1][1]))\n", - "solve(quad == 0, b_avg)" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "b1^2 + b2^2 + b3^2 - 2*b1*b_avg - 2*b2*b_avg - 2*b3*b_avg + b_avg^2" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "quad" + "2 * factor(simplify(D[1][1]))" ] }, { @@ -258,32 +252,12 @@ "metadata": {}, "source": [ "So, we end up deriving the quadratic form $$\n", - " b_1^2 + b_2^2 + b_3^2 + b_{\\text{avg}}^2 - 2b_{\\text{avg}}(b_1 + b_2 + b_3) = 0\n", + " b_1^2 + b_2^2 + b_3^2 + b_{\\text{sum}}^2/4 - b_{\\text{sum}}(b_1 + b_2 + b_3) = 0\n", ".$$\n", "\n", "This means that, given three mutually tangent circles with curvatures $b_1,b_2,b_3$, there are two solutions for $b_{\\text{avg}}$, allowing us to derive two new sets of three mutually tangent circles with curvatures $b_1' = 2b_{\\text{avg}} - b_1$ etc." ] }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "1/2*b1*h11 - 1/2*b_avg*h11 - 1/2*b1*h1_avg - 1/2*b2*h1_avg - 1/2*b3*h1_avg + 1/2*b_avg*h1_avg + 1/2*b2*h21 - 1/2*b_avg*h21 + 1/2*b3*h31 - 1/2*b_avg*h31" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "factor(simplify(D[2][1]))" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -295,38 +269,51 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ - "S1 = matrix([\n", - " [-1, 2, 2, 2],\n", - " [0, 1, 0, 0],\n", - " [0, 0, 1, 0],\n", - " [0, 0, 0, 1],\n", - "])\n", - "\n", - "S2 = matrix([\n", - " [1, 0, 0, 0],\n", - " [2, -1, 2, 2],\n", - " [0, 0, 1, 0],\n", - " [0, 0, 0, 1]\n", - "])\n", - "\n", - "S3 = matrix([\n", - " [1, 0, 0, 0],\n", - " [0, 1, 0, 0],\n", - " [2, 2, -1, 2],\n", - " [0, 0, 0, 1]\n", - "])\n", - "\n", - "# this is the only one that differs from the generators of the Apollonian group\n", - "S4 = matrix([\n", - " [1, 0, 0, 0],\n", - " [0, 1, 0, 0],\n", - " [0, 0, 1, 0],\n", - " [2, 2, 2, 3]\n", - "])" + "def weyl_generators(matrix, alphas):\n", + " retval = []\n", + " for alpha in alphas:\n", + " retval.append(identity_matrix(len(alphas)) - 2 * alpha * alpha.transpose() * matrix)\n", + " return retval" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "def standard_basis(dim):\n", + " return [ matrix(dim, 1, [0] * i + [1] + [0] * (dim - i - 1)) for i in range(dim) ]" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[\n", + "[-1 0 0 1] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 0 -1 0 1] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 0 0 -1 1] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 1 1 1 1/2]\n", + "]" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "S_i = weyl_generators(2 * M.inverse(), standard_basis(4))\n", + "S_i" ] }, { @@ -335,49 +322,38 @@ "metadata": {}, "outputs": [], "source": [ - "W = matrix([\n", - " [2, 0, 0, 1],\n", - " [2, 0, 0, -1],\n", - " [-1, 1, 0, 0],\n", - " [3, 1, sqrt(2), 0]\n", - "])" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [], - "source": [ - "def draw_circles_from_coords(coords, color=(0,0,1)):\n", - " drawing = None\n", - " averages = coords[-1]\n", - " avgb = averages[1]\n", - " for row in coords[:-1]:\n", - " b = row[1]\n", - " if not b == 0:\n", - " x = row[2] / b\n", - " y = row[3] / b\n", - " if drawing is not None:\n", - " drawing += circle((x, y), 1 / b, rgbcolor=color)\n", - " else:\n", - " drawing = circle((x, y), 1 / b, rgbcolor=color)\n", - " newb = 2 * avgb - b\n", - " if not newb == 0:\n", - " newx = (2 * averages[2] - row[2]) / newb\n", - " newy = (2 * averages[3] - row[3]) / newb\n", - " if drawing is not None:\n", - " drawing += circle((newx, newy), 1 / newb, rgbcolor=color)\n", - " else:\n", - " drawing = circle((newx, newy), 1 / newb, rgbcolor=color)\n", - " return drawing" + "S1 = S_i[0]\n", + "S2 = S_i[1]\n", + "S3 = S_i[2]\n", + "S4 = S_i[3]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Hmm. Doesn't look like it worked. And my crude code to draw the circles doesn't seem to be the problem either. I calculated a few of these circles by hand using the matrices and it comes out exactly where it's drawn here, so my code isn't the problem. Unless, of course, I computed the generators of the Weyl group wrong." + "We can test this out on the packing on page 3 of the Guettler and Mallows. The root sextuple is (-1, 2, 2, 4, 4 7), so the coordinates would be (-1, 2, 4, 6). After applying $s_1$, it should swap out $s_1$ for its pair, resulting in (7, 2, 4, 6). Likewise for $s_2$ and $s_3$. Then, for $s_4$, it should give the average between the triple (-1, 2, 4) and the other tangent triple, i.e. (-1, 2, 4, 14)." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(7, 2, 4, 6)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "root = vector([-1, 2, 4, 6])\n", + "S1 * root" ] }, { @@ -387,9 +363,8 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAEECAYAAADah2tfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbr0lEQVR4nO3dd3iT1RcH8G+6aaFl0zJlbxkFBVkiywEFBWXJENmgoqgMRRBRVByIFGQoQ0D4KVuUPUQ2BWSjbBDKbksLXcn9/fE1tIWWpm3evHmT83mePNE2eXNo1nnPvfdck1IKQgghhBAiYx56ByCEEEII4ewkYRJCCCGEyIQkTEIIIYQQmZCESQghhBAiE5IwCSGEEEJkQhImIYQQQohMSMIkhBBCCJEJSZiEEEIIITIhCZN4gIkCTSaTSe9YhBBCCGfglYXbSktwNxEdHY2goCBER0frHYoQQgjhCJkWCKTC5OT++OMPtGnTBkWLFoXJZMKyZcsyvc/mzZtRu3Zt+Pr6oly5cpg9e7bmcQohhBCuTBImJxcXF4caNWogPDzcptufOXMGzz33HJo2bYoDBw5gyJAh6N27N9asWaNxpEIIIYTrMmVh810ZktOZyWTC0qVL0a5duwxvM2zYMKxatQqHDx++97NOnTohKioKq1evtulxYmJi7g3JBQYG5jRsIYQQwtllOiSXlTlMwgB27NiB5s2bp/lZq1at8Prr72LvXuDoUeDECeDsWeDiReDaNSAx0YLkZAWz2QSzGUhODgBwGRUr5oa3N+DhAXh6Al5eQFAQULQoUKoUUK4cUKkSUKMGULiwLv9cYUB37/L1d/lyyuXSJeDqVSAhAUhKAsxmvt68vYGAACAkhJeiRVOuS5Tg61KIzCgF3LyZ8lpLfR0bCyQn85KUxGuzOeUzz/o69PICcudO/7WYPz8gS2Rcn1SYDORhFSaLBdi+HWjTZiby5m0JH5+SuHWLHwbx8Qp8mtO+oz08AF9fwGxOQGLiXQBmAJb/LgoFChSGUh6wWHh8sxlITOSHSnp8fAB/fyZVBQsymXr6aSAsDJBClXu6cwf46y8gIiLlcvQoX0tWgYH84ilSBPDzw70k3foldvt2SmKVmJhyv4AAoFYtIDQUqFOH1xUqSBLlrm7dAvbtA/bvT0nIrUnR/a8dgElOSAhff/cnRp6e+O/kMW0iFRPDY928mfZYPj5AcHDaJKpUqZTXZ758DvsziOzLNOWVhMlArAlTWFg7/PknsHw5k6S//+aHhfWpNJks8Pf3QJ48/FDw8rqOgwdXYeTILqhWzRtVqzKZ8fHh7RMSEpCQkHDvcWJiYlCiRIkMh+QsFn4QHTwIHD8OnDwJnD8PREayYhUdzS/K1IlVrlz8AKldG2jVCmjXTpIoV3X2LLByJbBiBbB5M79ovL2B6tVTEptKlVK+XAICbDuuUnydX7oE/PsvX3979zIJO3WKt8mXD3juOSbprVrJa8xVWZMjaxK+dy9w+jR/FxAAlCnzYBUo9XVwMJPz7EpI4Ofd/dWq1Ana6dM8YQWA0qVTXvuhofwczJ8/538HYVeSMLmC8+eByZOBCRO2I3fuOoiL87mXHHl4cDiscmWgcWNgxYqeaNw4LyZOnHjv/rNmzcKQIUNsbhNgrzlMN28CS5YAa9cCBw4AFy4A8fEpv0+dRPXsCbRoke2HEjo7fhyYP59J0sGDTJCefBJo3Rpo0ACoVo3VTK1Yv0A3b34whrZtgS5d5CzfqJKTgZ07eXJoTZCtyVHu3ClVHGulsXx556gyWizAP/+kxBwRwddo6iTKGvcTT/DiJZNk9CQJk1H9+ScwcSKwaVPq8m8y8uZNRq1afmjcGHjhBeDRR9Peb9iwYfjtt99w6NChez/r0qULbt686RSTvh+WRPn4cD5Uz55A794pFTDhnJKSmJxMmQJs3Ohc1R1rlWvlSr6HvL2Bzp2BgQP5BSWc2+3bwJo1fP5WrQJu3GByVLt2SpJhHYL1MNBab2sSlXqIet8+/nud6f3jpjKfhaaUsvUiNJSUpNTs2Uo1bKiUn59SHIBQKn9+s2rR4qaaNeuEAqC++uortX//fnXu3DmllFLDhw9X3bp1u3ec06dPK39/f/XOO++oY8eOqfDwcOXp6alWr15tcyzR0dEKgIqOjrb7vzM9V64oNXKkUhUqKGUy8d9tMilVtqxS77yj1OXLDglD2OjKFaXGjFGqaFE+Vw0aKDV/vlLx8XpHlr7ISKU+/lipEiUY72OP8b2WmKh3ZCK18+eVCg9XqlUrpXx8+FxVr67Ue+8ptWuXUmaz3hFqw2xWavdupUaNUqpGDf67vb2VatlSqW+/VersWb0jdBuZ5kGSMOnoxg2+SSpVUsrDIyVRKF1aqTffVOrCBaU2bdqkwOpemkuPHj2UUkr16NFDNWnSJM1xN23apGrWrKl8fHxUmTJl1KxZs7IUl6MTptSSkpSaMUOpJ55Qytc3JXEsWFCpTp34wSL0ER3N12tAAC/9+il14IDeUdkuKUmp5cv5hQwwIf/pJ9f9InZ2FotSe/cq9cEHStWqxefEy0upZs2U+uYbpU6f1jtCfZw9y0SpRQsmTgATqVGjlNqzh383oQlJmJzR2rVKhYamVFO8vfn/kyYpdfeu3tHpmzDdb9MmpZ5/Xql8+VKSp0KFlBo9WqmEBL2jcw/x8Up99ZVSBQqw+vnuu0z2jezgQaXatOHrqVYtpVavli8iR4mOVmryZKWqVOHfP29epbp0UWrhQqVu3dI7OucSFaXUokVKde2a8hlYuTK/K6Ki9I7O5UjC5Czu3lVqxAh+6VgrSXXrKrVypfOd4TpTwpTa6dNK9eihVK5cKWejTz+t1NGjekfmupYuVapkSaU8PZXq00epixf1jsi+tm7lkCKg1FNPKXX8uN4Rua6DB5Xq35/VSU9Ppdq3V2rNGhkatVVSklLr1in14ov87AsIUKpvX2NVeZ2cJEx6O3RIqebN+QEB8EXep49zn6E7a8JkZTYrNXOmUmXKpFSdypThUJ6zJZ9Gdf26Up0782/bpo1rJxIWC09cypdnBe2LL5RKTtY7KteQkKDUggWcmwkoFRLC6rCrJd6O9u+/Sn34oXHmERqEJEx6MJs5Bl2yZMoXeoUKSv34o96R2cbZE6bUjh9X6tlnecYFsPrUvTsnJovsWbJEqcKFOQQwf777DFXFxXHuoMnEOXSunCRq7dw5LuQoXJjvy6ZNlfr5Z6km2VtiolKLF3Pel3W6wvDhSp05o3dkhiQJkyOZzTx7sg4ZeXsr1a6d8SYvGilhskpK4hmX9QPa+iEtK+xsFxPDuSSAUm3buu/fbutWpcqVY7Vp0iT3SRjt4Z9/uDjDw0OpPHmUGjxYqSNH9I7KPRw7ptTrrysVGMik/8UXlTpxQu+oDEUSJkf57ju+UAGlgoKUGj+eX+JGMnnyZFW5cmVVoUIFwyVMqW3dmrLqxmRSqkMHpW7f1jsq53b6tFLVqvFLbt48SRLi4pR67TW+hnr2lKGOzFy6pNSAAaz0FivG9gDyntNHbCy/j0qU4FSQvn05hCcyJQmT1pYsUSo4mH9JPz9WOYw+j8aIFab0bN/OSoF1gviAAcZLYh1h40YuRihbVqoB9/vxR7a3qF/ffStuDxMVxaE3f38O4X7+uVJ37ugdlVCKC42+/FKp/Pk56jFsmFI3b+odlVOThEkrf/7JLxjrl/HAga7zZewqCZPVsmVpk9oxY4yf1NpLeDjPQps1c+6FCHratYuTlYsVYx8cwS/jCRNSvoxHjJCWAM4qKorNP/392cLhs88kqc2AJEz2dvx42uGeF190vdKzqyVMVtOmcbjUOmw6bZreEenHYmEXdYDzHlwl2dfKv/+yDUiuXOzZ5K6SkpT6/nulihdnot2/vwz3GMXlyzyx9/Li6rrp0+V9fx9JmOzl7l2lnnkmZULxU0+5boneVRMmpVhZ+vDDlO1nQkLcr3u42czJuIBSEyfqHY1x3Lmj1HPPcduO5cv1jsbxNmxg00RAqY4dlfr7b70jEtlx8mRKy5AKFdhIWSilbMiDDLRtoX5WrAAKFgR+/52bwx4/DmzYAAQH6x2ZyCoPD+CDD7jZ5YABwJUrwGOPAa+8wo0xXZ1SwGuvAZMnA999B7zxht4RGUeuXNw4unVroH17bgzrDmJjuWlxs2ZAoULA3r3AwoVA+fJ6Ryayo2xZYMECYP9+oFgxoGVLoF8/ICZG78icnyRMDxEfDzz9NNC2LXdmnzEDOHAAqFhR78hETnl5AVOmACdPAuXKAbNnA0WKAHv26B2ZdpQC3n6b/+7p0/khKbLGxwdYtIg7ynfoAKxZo3dE2tq4EaheHZgzh0n2pk1AaKjeUQl7qFmTJ/5TpzKBql4dWLdO76icmyRMGbBWldasAerWBS5fBnr31jsqYW+lSwP//AOMHg3cvMlqU69erlltmjAB+Oor4NtvgT599I7GuLy8gJ9+Alq0AJ5/HoiI0Dsi+4uNBQYNYlWpVCng0CH+v4d8Y7gUkwno35/Pb7lyUm3KjEkpZettbb6hkcXH80Nw9WqeTYaHu1+iFBMTg6CgIERHRyMwMFDvcBzmzBl+CZ46xaGHVauYLLuCVauANm2A4cOBTz7ROxrXcPcu0LgxEBnJyqSrDNFv2sSThqtXgc8+43CcJEquz2IBpk0D3nkHKFAA+P57oHlzvaNyKFNmN5C3QSrWqtLq1VJVckelS3OIbvRo4MYN16k2HTsGdOnCuTfjxukdjevIlQtYtgxITuacpoQEvSPKGWtV6amnWFU6eBAYPFiSJXfh4cF5nYcPs9rUooVUm+4nFSbwCzEsjGfh7lpVSs1dK0yp3V9t2rgRqFZN76iy7tYtJn6+vsD27YCbPp2a2rULaNKESen333OYw2j+/BPo1k2qSoLurzbNncvXuIuTClNmrl7l2dSqVUCdOlJVEnR/talmTc5bMRKl+CV+8yawfLkkS1p5/HFOop81ixPqjWb6dKBpU6B4cakqCUpdbSpdmnPZpkzhZ4o7c+u3xY4dwCOPABcvMpPeswfIn1/vqPQTHh6OKlWqoK6rTNyxgzFjgJ07WaHp0gV48029I7LdtGkcXl6wgEuJhXa6d+dw1jvvcBGBESQlMeZ+/XjZuFFeJyKtRx4B1q9nEj1oECeIJybqHZV+3HZIbsqUlDOphQu5RFiQDMk96OZNoHZt4Nw5oEEDYPNmrpZyVmfPcplwly5MnIT24uLYpy0kBNiyxbmrNNevAy++yKG48HCgb1+9IxLO7ocfmDDVqwcsXsypCi5GhuTS06MHs+XAQC6nlGRJZCZ/fuD0afbl2rYNKFECuHRJ76jSpxTw6quMecIEvaNxHwEB/FL5809g0iS9o8nYoUNc1HL4MPvwSLIkbNGrF1dQnjjB189ff+kdkeO5VcIUH88zwLlzgSpV+IVXubLeUQmj8PBgt/dRo7iUvEwZVpqczbRpHF6ZOVPmLTla48bA668DI0c659DcsmVA/fp8Xezdy3iFsFWDBnzdFCgAPPEEK03uxG0Spn/+AYoW5aTGTp2AI0cAf3+9oxJGNHYsJ1FbLFyC/eWXekeU4upV4N132ZiyRQu9o3FPn3zCz5pBg/SOJIVSwEcfscfc009zxWSpUnpHJYyoRAlg61b2devQgfM8jd56xVZukTD9+isrSlFRwMSJxlvtJJxPWBj7G+XLx+1GOnXSOyL6+GMuax8/Xu9I3FdAAPDFF9xmYsMGvaNhg82OHbmH4ocfAv/7H2MUIrv8/fk9+vHHPIHs0AG4c0fvqLTn8pO+Fy0COndmf6W1a6UEbQuZ9G27+HhOgvzrL3bF1XMvpjNnuM/hmDEcEhL6UYpDFsnJwO7d+vVmiotjJWDnTmDePOCFF/SJQ7iuFSu4uCQ0lMWJPHn0jijb3HvS95w5TJZy5eIQnCRLwt78/IB9+zg0t3490KiRfuXp0aM5t+CNN/R5fJHCZAI+/ZTzPfSa5xETA7RqxXYpa9ZIsiS0ERbGYsSBA3y9RUfrHZF2XDZhmjoV6NmTpedjx6S/iNCOhweHXp55hiukHnvM8UnTwYOsIHzwgQy3OIsmTThf6L33WGlypFu3OIft8GFWPRs1cuzjC/fyxBP8DDx+nE0ub9zQOyJtuGTCNHEiW/sHBXEJZMmSekck3MFvv/EsPiKCncEdmTR98gk78kqXeucyfjzw99+cN+QoN26w4nnyJFdL1qvnuMcW7qtOHbYdOHeOneOvXdM7IvtzuYRp6lR2Y86fnx8YRYvqHZFwJ4sXA127stdNnTqOSZoiI/m4r78OeHtr/3jCdjVrMnlx1JYpUVGsLP37L1te1K7tmMcVAmDbni1buFq3RQs2/HUlLpUwzZrFylLevKwsFSyod0TGIluj2Me8eeyivH8/+5ZonTTNnMlEqUcPbR9HZM/AgWx2qnWjv9u3OQR49iyH4apX1/bxhEhPlSqcz3nxIl+PrjSnyWVWyVlXw+XOzRJ4cLDeERmXrJKzj7ZtuYKkaVMOjWghOZlDcU8/DcyYoc1jiJxJSuKeXG3aAN99p81jxMVxDt1ff3EuSZ062jyOELbav5/V1SpVuOggd269I8qUe6ySW7EiZTXc4cOSLAnnsHw5y9KbNgGtW2vzGL/+yjO5AQO0Ob7IOW9vbj8yb542Z9sJCUC7dlytuXq1JEvCOdSqxUTp0CGeLNy9q3dEOWf4hOnECaB9e+4mf/CgTPAWzmXtWg7LrVrFDtz29v33XJUnc1WcW58+7Nm1aJF9j6sUk+WtW/kaq1/fvscXIicee4zbSe3axZMG2we0nJOhE6Y7d7gCxGzmmKm0DhDO6I8/gOLFuRHuL7/Y77hxcZyr0rGj/Y4ptFG0KNCwIavh9jRpEuduzpjBNgZCOJsGDbgp9bx5zrWNVHYYOmGqV4+rQiZN4pMihDPy8GCrgVy5uIXKsWP2Oe66dRyOadPGPscT2goL44ldXJx9jrd2LfDWW8A77wDdutnnmEJooVMn7j7w7rtsv2JUhk2YrEu3u3UDBg/WOxohHq5wYU7GtVg4bBIbm/NjrlgBVKoElC+f82MJ7bVpwwTXHtvn/PMPK4utWsm+gcIYPvqIczk7d2aDSyMyZML07bfAggVcNjt3rt7RCGGb+vWB8HBO/H388Zwdy2zmhO+wMPvEJrRXvjwT3JUrc3ac6Gg+70WKcANUT0/7xCeEljw8OCxXvDhfv7du6R1R1hkuYdq2jXtl5c3LDSWFMJIBA9gv6ehRnmll15497KQrw3HGEhbGRDe7k1/NZlbXL19mhTEoyL7xCaGlwEC+bq9f5+efo7cMyilDJUyRkdwR3tOTO4D7++sdkRBZN3s2O+IuXMhtfLJjxw6uDM1ppUo4VpMm7IJ8+nT27v/ee1x1tGgRUKGCfWMTwhHKlgV+/pnz+YYN0zuarDFMwpSczP4i8fHcl0nmbQgj27mT2/e89RZX0WWVdb862QrFWEJDeR0RkfX7zp8PfPYZ8MUXnLskhFE1a8aTxa++4gmkURgmYWrdmvsjjRgBPP+83tG4JtkaxXH8/Fgl9fTkl19MTNbuHxGR8uUrjKNIEc7hyGrCdPAg8OqrHM4dMkST0IRwqEGDuFl4v35sumoEhtga5ZdfuDdX/frA9u16ReE+ZGsUx8nOa/v2bc5dmTkT6NVL2/iE/bVrx1WS69fbdvukJDYANJuZZPv5aRqeEA6TmMj2QElJwN69nGagI+NvjXLnDtC9O/+Qq1frHY0Q9tWhA+fl7dgBzJlj230OHOCkYakwGVNoKCtMtp6rfvIJW6jMni3JknAtPj58XZ84wbYDzs7pE6awMO5BM3MmZ9gL4WpWruQChn79bBuaO3qUS3SrVNE+NmF/jz7KhruXL2d+27/+AsaNY9M/2f5GuKJHHwVGjQI+/TR7c/scyakTpl9+YbO/Bg2Al1/WOxohtOHnx/4kCQnA009nfvtLlzgXRiZ8G1OxYrzOLGFKSgJ69gQqVwbef1/zsITQzfDhTJx69uTnoLNy2oQp9VCckVupC2GL55+3fWju8mXuTSaMyfrcZZYwpR6K8/HRPCwhdOPtbYyhOadNmNq0kaE44V5SD81FRWV8u8uXgZAQh4Ul7KxwYcBkenjCdOCADMUJ95J6aG7vXr2jSZ9TJky//AJs3ChDccK9pB6ae/bZjG936ZJUmIzMy4tDqpcupf97GYoT7so6NPfKK845NOd0CZMMxT0oPDwcjzzyCPz8/PD4449j9+7dGd529uzZMJlMaS5+srTGMGwZmouMlAqT0YWEZFxh+uQT4PBhGYoT7sfZh+acLmFq21aG4lJbtGgR3nrrLYwePRr79u1DjRo10KpVK1y9ejXD+wQGBuLy5cv3LufOnXNgxCKnUg/NxcY++Pu4OCB3bsfHJewnd26eHN7v4EEZihPuLfXQnLM1tHSqhOnAATZzq1tXhuKsvvrqK/Tp0wevvPIKqlSpgu+++w7+/v744YcfMryPyWRCcHDwvUuRIkUcGLHIKT8/VpcSEliavl9yMod1hHF5eaW/8ejbbwPlyslQnHBvw4dzSHro0OxvVK0Fp0qYXn6ZkyF//lnvSJxDYmIiIiIi0Lx583s/8/DwQPPmzbFjx44M7xcbG4tSpUqhRIkSaNu2LY4cOfLQx0lISEBMTEyai9BXhw784ly8mENwqUnCZHzpJUwbNgDr1nFITobihDvz9gbGjwc2bwbWrNE7mhROkzBt2AAcOQI88wxQqpTe0TiH69evw2w2P1AhKlKkCCLv/xb9T8WKFfHDDz9g+fLlmDdvHiwWC5544glcvHgxw8cZP348goKC7l1KlChh13+HyJ65c3l21bVr2p97egIWiz4xCfuwWPg8WinFs+p69bh1ihDu7rnngIYNuX+ss3zeOU3C1KsXuxf/+KPekRhb/fr10b17d9SsWRNNmjTBkiVLUKhQIUybNi3D+4wYMQLR0dH3LhcuXHBgxCIj9etzG42NG4Fjx1J+7uXFlVTCuJKS0lYJFy/mUupPP2WVXQh3ZzLx/XDgALBokd7RkFMkTPPnA+fPc0guf369o3EeBQsWhKenJ65cuZLm51euXEFwcLBNx/D29katWrVw8uTJDG/j6+uLwMDANBfhHObP53XqKpOvLxdGCOOKj0/ZaDQ5GXjvPXZ5b9JE37iEcCYNGrAn4/vvc6NevTlFwvTGGxyzf0gRxC35+PggNDQUGzZsuPczi8WCDRs2oH79+jYdw2w249ChQwiRdeiGVLEi0KwZsH8/sG0bf1akCHBfDi0MJjKSzyMAzJoF/P0352wIIdL65BPgzBmunNeb7gnTV18BN24Ar78uO3Gn56233sKMGTMwZ84cHDt2DAMGDEBcXBxe+W/5VPfu3TFixIh7tx87dizWrl2L06dPY9++fXj55Zdx7tw59O7dW69/gsihefM4XN2jB///YT18hPNTKqVb+507wJgxQOfOQM2aekcmhPOpVg3o1g0YOzb9NiuOpOtaG4sF+OADICAA+OwzPSNxXh07dsS1a9fwwQcfIDIyEjVr1sTq1avvTQQ/f/48PDxS8t5bt26hT58+iIyMRL58+RAaGort27ejimxtb1jBwUD79lw9unQpu3z//bfeUYnsunGDc5iKFgW+/Ra4etU5m/QJ4SzGjgUWLgQmTtS35YZJ2d7kwO7dEN59F5gwAfjyS+Ctt+x9dJFdMTExCAoKQnR0tMxnchKxsUC+fEDBgtw2Y+FClqmF8Rw6xOZ8q1cDnTpxftrkyXpHJYRze/NN4IcfgFOn+DmogUyXW+g2JJecDHzzDVCggCRLQmQmd26gb1/OfTl7lvuQOVNDN2E76x5yK1ZwIqs0qRQicyNH8jPv88/1i0G3hGnCBH5YjB2rVwRCGMs337Ch28aNfO9IhcmYjh3jIpdFi5gE27jgVQi3VqgQ0L8/MGNG+tsKOYJuCdOkSZzk3b+/XhEIYSxeXmzsat1GMCJC33hE9kREAMWLcy7TgAF6RyOEcfTvD0RHc0qCHnRJmHbs4NDCCy9w9Y8QwjZffcVrX182OhTGs3cvN1Bu0QKoUEHvaIQwjjJleNIYHq7PlARd0pV33uH1l1/q8egiI+Hh4ahSpQrq1q2rdygiA2XLAuXLc0huzx69oxFZdfs2cPw4+2gNHKh3NEIYz8CBwL59+nz+OXyVXGwsEBjI3goHD9rjiMLeZJWcc5s7lz2ZrB2/ZSsN49i6FWjcmE0rL16UTZSFyCqzmRuTN2kCzJ5t10M73yo560z3ceMc/chCuIbu3ZksJSSk3WNOOD9r0/4BAyRZEiI7PD05l2nhQs4DdCSHJ0xz5gB58wJhYY5+ZCFcR6dOvJ44UdcwRBbNncuKYL9+ekcihHH16sXCy6xZjn1chyZMixcDMTHAf7t6CCGyyTr5e8ECfeMQtrt6la0g6tSRVgJC5EShQsBLLwFTp3LHEEdxaMI0ejRXxclwnBA5kz8/8MgjXG0l26QYg3WRy3vv6RuHEK5g0CDg9Glg7VrHPabDEqZLl4AjR4BGjQB/f0c9qhCu64sveN2nj75xCNvMmwfkyiXTEYSwh8cfB2rVAqZMcdxjOixhsm6u+/HHjnpEIVxb+/acALlzp96RiMxcucKTxqZNZVWjEPZgMnEu0++/s5mlIzgsYVq+nGdXDRo46hGFcH1Vq7In0759ekciHmbMGF4PG6ZrGEK4lLAw7ku7erVjHs8hCdOdO8D584D0QxTCvqzzAQcP1jcOkTGluG9cnjzswSSEsI+SJYGaNbmRtSM4JGH67jt+aMhcCyHsq00bLqTYtQtIStI7GpGebduAW7eAZ5/VOxIhXE9YGPDbb475/HNIwjRvHj/Uu3RxxKOJ7JKtUYypenUurf3lF70jEen56CNev/66vnEI4YratAGiooA//9T+sTTfGsViYVfi8uWBo0ezcwThaLI1irHMnMnqbdmywD//yKRiZxIZCRQrxpXB0dGy2bgQ9maxACVKAB07pvSnyyb9t0b59VdOyurYUetHEsI99ezJJOnUKcf2JBGZs1aX2rWTZEkILXh4sMq0YgWn/mj6WNoePqVHwhtvaP1IQrgnLy9uRmkyAcOHO7bzrcjYqVPAtGl8Pjp00DsaIVxXWBjfb1rvral5wrRtGxASwv3jhBDa6NCBZ1cHDnBFltDfqFFAQADg5wc0b653NEK4rqee4rC31qvlNE2YDh8GYmOBp5/W8lGEEEOG8LpIEeD999mbSehn/37gp5+4Z1zz5kychBDa8PMDWrY0eMI0aRKvhw7V8lGEEIUL83L3Ljd4nTxZ74jcl1LA229zmPTkSc6vEEJoq00b7npw44Z2j6FpwrRzJ+Djw27EQght1aoFxMQAAweyyvTPP3pH5J6mTwc2bgT69uX8pYYN9Y5ICNfXsCFPViIitHsMTROmM2c4f0kIob0nn+R1ixZ83/XqJRPAHe3sWVaX+vThJPyAAKBiRb2jEsL1lSvHbvqGTJji4zl/qUYNrR5BCJFa+/a8Xr8e+OEHNnL79lt9Y3InSgG9ewP58wNffMEP7po1uUGyEEJbHh5A7doGTZhWreK1rA4RwnZKATdvcnf7O3eydt/y5dliYMcOoEkT4LXXgBEjZGjOUaZNAzZsYCPRwEB+cIeG6h2VEO6jTh2DJ0zWs14hRPoiI4GPPwaaNgXy5QMKFODqqoAAJkGdO3P1h9mc+bGKFOFEYwAYP55dpjt0YLVXaCciAnjzTc5batGCXb3/+UcSJiEcKTSUw+JaTfzWLGHau5cTvosW1eoRhL3JXnKOdesW57qULMmEKV8+4N13uSfc8uXA7NnAc88BJ04AbdsCZcoA//vfw49ZtSq/rJOTmXAtXQqcPs1u4DKfSRuRkXx+Hn0U+OYb/mzfPl5LwiSE41jfb1pVmTRLmGTCt/EMGjQIR48exZ49e/QOxeVt3AhUq8YE6NNPgX//BZYsAUaOZFU2LAzo0QOYOJFfvnv28MOgY0fgxRe5Gi49TZvyesMGXlerxs2vFy9O2aZD2E9CAvD880xGly5lPxiAH9j+/kClSvrGJ4Q70XrityYJk0z4FiJjK1eymWvlysCRI8Bbb7G69DB16jChWrSI+8U1b85K0v2sQ+C//prys7ZtgXHjgDFjmDgJ+1AK6N+fTSqXLUtbTY+IYJsHmfAthONoPfFbk4RJJnwLkb79+1khatMG+P13oHjxrN3/pZeAzZs5P+allx7cbDL1xO/URo7k7bt35+o5kXPjxnHYdOZM4LHH0v5OJnwLoQ8tJ35rkjD9/juvZcK3ECkSEzmXqFIlYMECwNs7e8epVSul0jR9+oO/L1KEG1GmZjIBs2bxi/3ZZ4Hdu7P32IImTAA++IBJ08svp/2dTPgWQj9aTvzWJGHat49fBjLhW4gUU6cCR4+yKuHrm7NjtWzJCePvvPPg0FzlyukP1/n7cziwWjWgVStJmrLryy85Of+993i536FDvK5Vy7FxCSFS3ncHD9r/2JokTNeuceKVEIIsFiA8nMNxNWva55hjxrBX048/pv15uXIcqrt+/cH75M7NCnDlyhwyl+G5rBk3jp28hw/PeBL9v//yulQpx8UlhKCSJXltfR/akyYJU3Q0kDevFkcWwpi2b+cwTf/+9jtm0aJAu3YcakutTBleZ3SGFRQErFnDyZEtW7KNgXi4pCQ2Ah01Chg7FvjkEw5zpufSJVbz5KRRCMfz9+dn3KVL9j+2JgnT3btAoUJaHFkIY9q5k2/kBg3se9ynngIOHOB7zqpyZV4fP57x/fLkYaWpbVtWvUaPlj5NGblxg6sav/uOw6qjRmWcLAHA5ctMZh92GyGEdooW5fvQ3rzsfUCLhU3zZP6SECn27cv5MnOlgG3b2Bpg714mStYO3oULc3VInTpA48b8WWZbouTKxcnn1asD77/PuTdz53LYTtCRI+yJFR3NPfqaNMn8PpcuSQ86IfQUEmKQCtOFC7yW8XshUty4we1OskMpYOFCzn1q1IgJU7FirHR89RVv88IL3FJl/nx+wQOsamXGZGLLgWXLgHXrgCeeAP7+O3txuprFi4F69dgxfc8e25IlgGe2kjAJoZ+QEG0qTHZPmKwrRMqXt/eRhdZkaxTtmEwP9kyyRWQkO0l37gyUKMFWAmfPMoF6910OpwFAp06ci3ThArdVMZmYMHXrxs18MxMWxtvHx7Ph7Fdf2bZ3nSu6cQPo0oV78LVqxflnpUvbfn/rkJwQQh9FixqkwnTsGK9lSwDjka1RtBMSApw/n7X7HD/OIbYdO9jl+9dfubGrR6p37blzvLZWrzw9mfwEBHAe4a+/skpiy2NXrcphvv79uRKscWP3qzYtXQpUqQKsXs0tZX7+OetDlDIkJ4S+rBWm7JykPozdE6bTp3n96KP2PrIQxhUaylVriYm23f7MGU7ozpuXSczzz6d/u4gI9nSqWjXtz/PmZYUoIoIrvJ56CrhyJfPH9fcHvv4a2LIFuHqV1aYJE7hnmiu7dIlVpRdeYIJ55AjQtWvWJ27fucP5TlJhEkI/RYvyvXj7tn2Pa/eEyXrGW7CgvY8shHE1bMhkac2azG9rNvPL28+Pm+g+rFqxciXw+OOAj0/anxcoAMTFscXAxo2cHN6rl+1nXI0aAX/9xWrT8OGsGP/4o+sN00VFASNGsHfVmjWsKi1blv0KkXXehFSYhNCP9f1n72E5uydMkZEPfngL4e5q1uTw2pQpmd/266+BXbuYoBQpkvHtTpzgyq0+fR78XUhISjWrdGlgxgzgt9+AOXNsj9labTp0iPF3786VfqtW2b/U7Wh377JyVqYM8M03wJtvcjuZ7FSVUrMmTFJhEkI/1vefvSd+2z1hunGDH7RCiLSGDOHcmLVrM75NbCwbIw4a9PCeTUoBQ4cyMbJO/E6teHHeJiaG/9+mDSeOv/ceh+iyokoVzu3ZsQPInx9o3Zqr6RYsMN5Q3ZUrbDpZvjxXB3bqxETp44/t02w3MpLXD0t0hRDass7ptL4f7cXuCVNSUvY3FRXClXXpwu1IevfOeOXa/PkcSnvnnYcfa9YsVnqmTUt/XzrrRGVrnyaAq+ouXeIwXnbUqwds2sRKlb8/qzElSjDxsA7FOyOluAVMly6M96OPuPrt2DFW/Ow5fGZNIOWkUQj95MrFa3uf0Nk9YTKbc9acTwhXZTIBM2dyOOiZZ4Bbtx68zezZwHPPpeyHlJ6VKzm3qHdvVo7SYz1pSd0BvGZNoH59PkZ2mUyMfcMGbiTcuTP3yCtTBnj2WQ79adH/JKuU4iT7ceM4cb1RIzb7/OwzJo3ff895S/Zmrd552b0lsBDCVtaVxFmtpmd6XPsejp2+PTTZcEUI4ytVipOLT55k8rJrV8rvEhPZEbx58/Tvm5TE4brnn2ei9LD5UNaqU3x82p83b87HtMccpMqVOf/n0iVuGxIby0SuaFFORB83jhPHHTVR/O5dNt987TXgkUeYKH3+OeNcu5ZtGt58E8iXT7sYkpN5LSeNQujHZOJJi/X9aC+abI0iHxZCZKx2bW5x0q0b5wK9+ioweDATi8REtiBILSGBTSk//5zL3YcPB8aMeXgVw1phur8kHRrKdgH//st5TvYQEMCJ5336cA7jb78BK1awmjNqFIcHa9XiY1sv5cvnrApz5w4rSBERrBxFRLDiZTazOhcWxqSySZP0hyy1kpzMf5fsIyeEvnRLmEwmkyk6OtqmA5rNuWEyATExsZnfWDiFhIQEJKT6Zr39X/OKGOuMYWF3RYuy0jRlig/Cw30wY4YHgoPNADzx/vtJyJWLb/aoKBP++ccTUVEmNGmSjPXr41G7tgV37jz8+BaLN4BcuHkzFjExKbvqFi7sASA3/v47DoGB9i/9eHtzQ9+2bZms7dzpiX37PHHggCeWLfPExIksP3t4KBQqpFCggEJgoIK/P+DtrWAyPdgV3WLhse7cMSEqyoRr1zwQFWX67/EUKle2oFYtM3r1MuOxx8yoWtVyL2FJSHDsxPTYWG94efkhJsbODWCEEFni5ZUHsbEJiImxrfldUFBQIIDbSj2k/q6UyvQCIBCAsu1ySwGnbbytXOTi7pfcChikgLMKUP9dLKkuSgFJClipgIZZOO7Q/+57/30q/ffzJxz4b/RQQBUFdFPAVAUcVEBUqn9f6oslncv9t0lSwFUFbFXAOAWEKaCoEzyXUMBrCohzgjjkIhd3v0Qpfg5m6X6BD8uFTA9Lpqz+qzBZMr0hgJIlAxAdfREXLngiMDDQlrtkSd26dTXZukOr4xrh2PdXmC5fvozHHnsMR48eRbFixXIa4gPkOaRVq7wwZIgfrl61ICxMoVIlMz7/3A8REbEoV86C+HjgyBEPbNvmhTlzvHHypCdat07CV1/Fo0iRh79vx42zYMKEvFi58goaN8517+e7d3uiRYsA/PlnLKpXt+ktna6H/T2UAo4d88Bvv3lh/Xov/PWXJ+7cYcmnbFkzata0oEIFM4KDFUJCFIoUsSAkRKFgQQVPzwePHR8PXLliwpUrHrh82YQrV0z4918PHDnigf37PXH9OqtWhQtb8NhjZjz9dDJatUpG4cJp/0aOeH3MmOGNkSP9cO2afSpMMTExKFGiBC5cuGD3z1OjvV+0PrZWx9XyOQSM9/dw1LFDQvLggw8SMGCAzRWmIGRSYbJpSO6hJar7eHpaAHggMDBQkxeHp6c2iZhWxzXysfPkyWOov7VR/s5mM/D665y03bo1cPx4CyxfvgGRkd74/HPg9OncqF0bCAwEChcGmjbl0v3Fi4HBg71Rr543li9n9/CM4+Vs7/z5AxAYmLIZ2j//cGw/NDQ3/Pyy/2+4/++RnAz88QfnLq1Ywa1dAgKAli253UidOpzHFBTkCcATQMa9R+4/tvXvUL36g7dVCrh4kXOYIiI8sHmzB157jcd+/HHOZWrblr2kHPH6yJOHfwt7P44Wn6dGeb846thaxgxo8xwCxvx7OOLYyclA7tx+CAy07YNOKZXpHBS7T/rmCjntZn0PGjTIUMc18rG14s7PoVJAv37so/Tdd0DfvsCUKS8AYLO1YsWA3buBDh3S3s/Dgw0qmzbldatWXPmVUXPLhARWdO6f8LxnD1CtGnKULAEpf4/ISLZKmDaNiUvx4ikTrp98MnuPk5W/tcnE3kolSgDt2vFn166lTDz/+GMmm/XrA3XrfoP4+Jz/2x8Ws5cX51wZYbWwEd4vjjy2ET9LAWP+PRxxbOsCDLuyZQ7TfxebFC5sVsBlFR0dbetdhJO5cOGCAqAuXLigdygu57vvOA9n7tz0f//qq0oVL65UUlLGx7hzR6lGjZQKDlbqxo30bzN4cLwClDp4MOV9GBurVFCQUsOHZz9+pZSyWJTavFmpl15SystLqVy5lOrdW6ndu/k7Z3L3rlKLFyvVvDn/7gULKjVsmFKnT2vzeHPm8HESEuxzvOjoaAVAPk8NTJ5DxzOb+T6cOTNLd8s0D7L7OZCXlwm5cuWBryPX8gq7sj538hza17lzwNtvc/l9t27p32bgQFZqHtaNO1cuYOFCzu154430b5OczCpvYGDKc7hgAbdK6dcvu/8CbuT7+OOsHv31F/Dll+zDNGMGULeu8y2n9/PjkOC6dezD9PLLrOyVLQv06GH/DuXWdg72apjn6+uL0aNHy3vRwOQ5dDxrOwGnrzBVqKBUQECWsjrhZOSMSBsDBihVpIhSmf1ZmzRRqlw5peLiHn67GTN4FnX06IO/e/FF/s5s5v/fuMGKVPv22QpdRUQo1bIlj/nYY0qtWeN81SRbxcUpNXkynwsfH6XeeEOpq1ftc+zffuPf6Px5+xxPCJF1ly/zfbh8eZbu5vgKU5EiD3YXFsLdxcQAP/7I6k5mcx2t84FGjnz47bp1AwoVYsXkfpcusYGshwfnTb3+Ot+XkyZlLe7z57lBbWgoqzGLFwM7d3Iyt7NVk2zl78/NjU+eZGPNWbNYcfroo5x/dmm1S7oQwnbW95/1/Wgvdk+YihXjKiB7d9gUwsjWrePWIb16ZX7bihWBTz/ltiOTJ2d8O19fDjEtXvzg765e5XCUUtxOZf584Ntvbf8AUQqYPp0TxLdu5ZDb4cMc3jJqonS/3LmB998HTp3iMOm4cezCvnt39o9p3cj30iX7xCiEyDrr+8+eG2sDGiRMjzzC6xMn7H1kIYwrIoLJSqlStt3+9deBt97ivmijRmU8J6Z+fW5zcuVK2p/fugXkycM5TmPGAJ98wuTKFufOsYLUrx/QsSO3HOnd23U3lC1YkHOxIiLYCqF+fW4/k51qU8GC/DtJhUkI/Vy+zBO7IkXse1y7J0wVKvD68GF7H1kI4zpyBHj0UdtvbzIBX3zBqsf48ZxonV6ftxo1Uo5vpRQQHQ3cvMnhvcmTgREjMn9MpVhJql6dE6RXr+b/BwXZHreRVasG7NjBv/nXX7PalNXeeh4ebA8hFSYh9HPpEnu32fskz+4JU5UqvD5xAvj444/xxBNPwN/fH3nz5rX3QwlhGHfusOKTFSYT8N57nDOUmAg89hirHzNnpmw0az3m7dvcjHbKFKBChVgkJQFJSTFITKyKYsWWZfpYCQkcLuzbF3jpJZ7wtGqV9X+n0Xl5MbmMiOBcp4YNgTlzsnaMkJCcVZjGjx+PunXrIk+ePChcuDDatWuHE1KyN5ypU6fi0Ucfvdewsn79+vj999/1DsstXL5s/+E4QIOEydqN98wZIDExES+++CIGDBhg74cRGggPD0eVKlVQt25dvUNxOb6+wN272btvnTrAgQPAkiWcd9OnD1C1KiePh4byNi+9xGrT4MFA3rzcmbdJk1sAjmZ6/MhINsT86SdOTJ85032qShmpVg3Yto0T63v2BIYOtX1eZtGiOUuYtmzZgkGDBmHnzp1Yt24dkpKS0LJlS8TFxWX/oMLhihcvjk8//RQRERHYu3cvnnrqKbRt2xZHUpeDhSYuX7b/hG8A9m8roJRSJpNSzZql/P+sWbNUUFBQltb3Cf1IWwH7e/NNpcqUsc+xbt1SasMGpb74gs0jAaVGjVJqyxalYmKU2rmTP/vwQ6UAqKVLl2Z4rD17lCpWTKmQEKV27bJPfK7EYlFq0iSlPD3ZVuHmzczv07+/UjVr2i+Gq1evKgBqy5Yt9juo0EW+fPnUzCx2UxRZV7cumwBnkePbCgA8m75/EqoQ7qxOHeD0aW7bkVN58wJPPcWqR4UKnGj84YdA48Ycojv6X1GpUqWHH2f5cqBRI56J7d3LIT+RlsnEiferV3M+U716mTe7zGmF6X7R0dEAgPz589vvoMKhzGYzFi5ciLi4ONSvX1/vcFyeVhUmTRKm3LmBGze0OLIQxtS8ObtAz5tnv2OazWwX8OyzaZf6//03r6tWzfi+//sf96t77jlgyxaNytcupHlzthtISmJieupUxrcNCWFbB3u0VrFYLBgyZAgaNGiAatWq5fyAwqEOHTqE3Llzw9fXF/3798fSpUtRxTrRV2jCYuE0A6eYwzR8+HCYTKaHXgICEhCT6b6/QriPwoW5aW54OCdw28PKlZwreP8+lmfP8rpixfTv99NPQOfObEi5cCG3WhGZK1cO+OMP9rdq3JiNL9MTEsIVh5GROX/MQYMG4fDhw1i4cGHODyYcrmLFijhw4AB27dqFAQMGoEePHjh6NPN5hSL7rl/nyYpTJExDhw7FsWPHHnopV84LcXHSvFKI1IYN43DOxx/n/FjR0ezV1KIF93BL7cgRVrPSW1K7eDEnMnfrBsye7bq9lbRSvDgrcnnyAM2apT88V64cr48dy9ljDR48GL/++is2bdqE4sWL5+xgQhc+Pj4oV64cQkNDMX78eNSoUQPffPON3mG5NGs+an0f2lOWPy4LFSqEQoUKPfQ2zZoBGzYAmzbxA10IwT5M773HhOnJJ7kyLTssFjaVvHWLfZLu77x95kz6Q2zr1rGq9OKLwPffc+sUkXXBwfx8a9yYc8l27uQWNVblyzOhiojI3uefUgqvvfYali5dis2bN6N06dL2C17oymKxICEhQe8wXJq1JUhmczizQ5M5TB068HrBghgcOHAA58+fh9lsxoEDB3DgwAHExsZq8bBCOL2RI5kotWkDrF+f9fsnJrJf0s8/cw+0+zuHx8cDsbEKjzwSjQMHDgAAzpw5g+XLj6JDBwuaNwfmzpVkKaeKFWPSdPs2P+9SD7N6eLDpZURE9o49aNAgzJs3DwsWLECePHkQGRmJyMhI3M1uXwqhixEjRuCPP/7A2bNncejQIYwYMQKbN29G165d9Q7NpUVEADVralQ9t2UpncpiWwGluAy3QIFTCsADl02bNmV5vZ9wHGkroK3YWKVatWL7jTffVCouzrb7RUQoVb26Ul5eSi1YkP5tlixhSwHgtVTvuUAFHFWBgZfUrVt2+2cIpdTWrUp5e7OVQGpvvaVU6dLZO2Z6n5kA1KxZs3Icr3CcXr16qVKlSikfHx9VqFAh1axZM7V27Vq9w3J5FSoo9dpr2bprpnmQSSllc26VlUSsWDF2N751Kyv3Es4gJiYGQUFBiI6ORmBgoN7huCSzmdtvvP8+h2969+a+bVWrcv6R1c2b3Px2+nTg9985rDd7Ns+g0tOnDxtPXrjA+TZmMxAWxiaMu3enbF0k7Of77/n8TZkCWHv0LlgAdO3K1cLSDUAIx4iOZtuV2bOBHj2yfPdMtxXXZEgO4Ad/VJRM/BYiPZ6ewNtvc4J2167A1KlArVrs3l21Kv+7VCmgQAGgXTuuuJo5k0lPRskSwN97ezNZAjhnavVqYNEiSZa08uqrnID/+uucEA6kdGDft0+/uIRwN/v389r6/rM3zdbINGnCSaYy8ds4wsPDER4eDrPZrHcobqNsWWDiROCTTzj2HhHBBpfJyaw81ajBppflyz84uTs9qSd8b9kCfPYZL+64L5wjffkl9/Lr1g04dCjtxO/mzfWOTgj3EBHBNilaTPgGoN2Q3IkTDPq114BJk7IemNCPDMkZU3w8Pyxat2Z/pUcf5dD45s2ciCy0de4c96Dr3JlDqE2asP/Wzz/rHZkQ7qFLF/ah2749W3fXb0iuYkUOO+zYodUjCCFSs26E3qwZMHw4twf44QdJlhylVCngiy/Y6mHtWg4LZHelnBAi6yIitBuOAzRMmACgSJGMu+EKIezrt994Xbw4MHkyMH68Ns3bRMb69uUQXO/eQJUqHCK9eVPvqIRwfdHR3BaqTh3tHkPThMk68dteW0EIITK2axcnfA8fzk11X3tN74jcj8nEyfm3bgEbN/JnUmUSQnvWBRaGrTC1b8/rWbO0fBQhhMUCHD8O5MvHqsZ338lQnF5KlQLGjOHKxEKFgDVr9I5ICNe3Zg3fb5Ura/cYmn6kvvIKz7gkYRJCW6tXA0lJ7DzdsyeHg4R+Bg3iasU8eYDly9lKVAihnRUrgOee03YXA00TJh8fLpu29kYQQmhj8mReJyezuiH05ecHfPghW0ScPMlVw0IIbZw8yc2uw8K0fRzNi/YvvMA5TBs2aP1IQrivrVt5/frrQIkS+sYiqHt3rhb28ODZrxBCGytXAr6+2vd81DxhevNNXn/7rdaPJIR7OnYMiI3lZpMjRugdjbDy8uJKRYsFmDdP72iEcF0rVrCdSu7c2j6O5glTcDAnYlm3DBBC2Nfnn/P65Ze5lYpwHu3a8TPw8GHg2jW9oxHC9Vj329R6OA5wQMIEsEwWFQWcOuWIRxPZFR4ejipVqqBu3bp6hyKyYNkyXr//vq5hiHSYTJwArhQ35BVC2Nfvv3OT8dattX8szbZGSW3vXqBuXe7kPWVKdo8iHEW2RjGOqCi2EihYUCoYzioqipW/SpW42bIQwn46duTiij17cnwo/bZGSa1OHcDfH/j1V0c8mhDuY9gwXvfrp28cImN587KZnnWumRDCPhIT2VKlTRvHPJ7DWts9/jhw4YJ8YAhhT//7H69Hj9Y3DvFww4ZxWM4630wIkXN//AHExDhm/hLgwISpf39ef/aZox5RCNcWHc3hniJFuCWKcF4vvMC+dDKPSQj7WbQIKFkSqFHDMY/nsITppZc4LDdtmqMeUQjX1rs3r996S984ROZMJuCJJ7jw5epVvaMRwviiongC0rs331+O4NDdpjp25MRU66aUQojsW7GCHxRvv613JMIW77zD67Fj9Y1DCFcwdy7nMFlPHB3BIavkrK5fBwoX5iTw3btzejShFVkl5/yWLuUwT9Wq7PEjnJ/Fwip7rlzAjRuyObIQ2aUUN9mtUYPDcnbiHKvkrAoWBGrXZpuBmzcd+chCuBbr6rj33tM3DmE7Dw+geXMOJaxbp3c0QhjXpk3cn3HgQMc+rsPPccaPZ3Y4fLijH1kI13DpEvDPP/xvRy2nFfbx8su8/vprfeMQwsimTAGqVAEaN3bs4zo8YWrRgk3cZLWIENkzdCivH3lE+72ThH1Zm+ivWwecO6dvLEIY0b//cneDgQMdN9nbSpdR9L59gbg42ZBSiKyyWDh/ydsbaNRI72hEVpUpAwQF8fmbPl3vaIQwnhkzAD8/oFs3xz+2LgnTBx8Anp7ARx/p8egiI7KXnPObMgVISGDiFBqqdzQiq0wmPm/FigEzZ/K5FELYJimJJxrdugF6rEfSJWHy8wOaNgX+/hs4c0aPCER6Bg0ahKNHj2KPHTblEdqYMAHw8uJmk5IwGVNoKCvsV68CS5boHY0QxrF8OXD5Mvel1YNuC1utkx4dPctdCKPasQM4fx6oXp3/X7OmruGIbKpdG7hyBWjYkJ+Dtnd2EcJ9KcX3S8OGwKOP6hODbglTtWrsIbNmjUx+FMIW3btzSKdlSy6ckAnfxlSyJK+7d+cO68uW6RqOEIawahWwfbu+rVR0bZ02bx6zxs6d9YxCCOe3fDlw8iTQvj03sA4J0TsikV3W565UKa4afu89IDlZ35iEcGZmMzBiBPDkk0CrVvrFoWvCVLMm91fasQPYt0/PSIRwbv37c+7SrFkcwy9aVO+IRHZZE6bLl9mX7tgxbvMghEjfggXc0WD8eMe3EkhN9+b8P/3EP4AeSwSFMILp04HISKBPHw7DXb4sFSYj8/MD8uXj8xgayo3JR48G4uP1jkwI55OQwJX1zz8P1Kunbyy6J0wlSwLPPQccPSrbBdzv5s2b6Nq1KwIDA5E3b168+uqriI2Nfeh9nnzySZhMpjSX/v37OyhiYW8WC/Duu/ySnTSJP7t0SSpMRle0KJ9HABg3jslTeLi+MQnhjKZN42KXjz/WOxInSJgA4Mcf2Zfp1Vf1jsS5dO3aFUeOHMG6devw66+/4o8//kDfvn0zvV+fPn1w+fLle5fPP//cAdEKLXz0ERAdza2EvLz4s+hoIG9eXcMSOZQvH59HAChfnjuuf/JJys+EEMDt2zyh6NmTm+3qzSkSprx5OSR34YJ0/7Y6duwYVq9ejZkzZ+Lxxx9Hw4YN8e2332LhwoW4ZD01zYC/vz+Cg4PvXQL16PAlciw5Gfj0U3aGHjUq7c+9vfWLS+Scl1faid4ffADcvcs+W0II+uorICYGGDNG70jIKRImAJg6FfDxAd54Q+9InMOOHTuQN29e1KlT597PmjdvDg8PD+zateuh950/fz4KFiyIatWqYcSIEbhz585Db5+QkICYmJg0F6G/11/nvJbPP+dO91bJySnVJmFM3t7sWmxVtCg/+77+msNzQri7q1eBL74ABg8GSpTQOxpymoTJzw8YMgS4eVPOsgAgMjIShQsXTvMzLy8v5M+fH5GRkRner0uXLpg3bx42bdqEESNG4Mcff8TL1i3SMzB+/HgEBQXdu5RwllenG4uN5Z5JwcHcezE1PVeJCPtQKm0SDADDhgG+vhyCEMLdffIJ3yMjRugdSQqnSZgALhnMnZvDD65a5Bg+fPgDk7Lvvxw/fjzbx+/bty9atWqF6tWro2vXrpg7dy6WLl2KU6dOZXifESNGIDo6+t7lwoUL2X58YR/PP89K0nffPfg7L6+01QlhPElJD1YJ8+YFRo7kJFdpsyLc2cGD3Ddz2DA26XUWTlXY9/AAZs8GOnQAnnkG2LZN74jsb+jQoejZs+dDb1OmTBkEBwfj6tWraX6enJyMmzdvIjg42ObHe/zxxwEAJ0+eRNmyZdO9ja+vL3x9fW0+ptDW4sXA+vVA/fpA27YP/t7bG0hMdHxcwn4SE9MfVn3jDWD+fE5y3buX0xSEcCdJSXz9V6wIDB2qdzRpOVXCBLCTcbNmwIYNbObWvbveEdlXoUKFUKhQoUxvV79+fURFRSEiIgKh/+2yunHjRlgslntJkC0OHDgAAAiRxj2GcOcOF0D4+gKrV6d/m4IFgWvXHBuXsK9r14AGDR78ubc3Txrr1OEKyY8+cnhoQuhq/HhWmHbt4uegM3GqITmrFSuAXLk4d8NVh+YyU7lyZTz99NPo06cPdu/ejW3btmHw4MHo1KkTiv7XhOfff/9FpUqVsHv3bgDAqVOn8NFHHyEiIgJnz57FihUr0L17dzRu3BiP6rVbociStm25WmraNCCjxY0hITIx2MiUenjz0Ro1OC1h/HgZmhPu5a+/eJIwfDibujobp0yY/P3ZmykhgUNz7mr+/PmoVKkSmjVrhmeffRYNGzbE9OnT7/0+KSkJJ06cuLcKzsfHB+vXr0fLli1RqVIlDB06FO3bt8fKlSv1+ieILEg9FNejR8a3CwlJaXoojOf2bSAu7uHd2keMAKpX59CEDL8Kd5CUBLzyClCpUto2Ks7EpJSy9bY239Bemjfn0NycOa43NOfMYmJiEBQUhOjoaOnh5CB37nCozWLhctqH/dmHDgV+/RU4ccJx8Qn7OXGCXwqbNwNNmmR8u7/+4tDc8OEyNCdc39ixvOzapVt1KdP1x05ZYbKSoTnhLmwZirOSCpOxWZ+7zKYVytCccBfOPhRn5dQJkwzNCXdg61CcValS7NN03yJKYRCnTrGXVrFimd9WhuaEqzPCUJyVUydMQMqque3buWpOCFdiy6q4+9WqxeuICO3iEtqJiOCXQ0BA5re1rpo7dkyG5YRrsq6Kmz3b+VbF3c/pEyYgZWiuTx/uWiyEq2jWzPahOKuyZbm/nCRMxhQRkbVhh9RDc3/8oV1cQjja9u3sbO/sQ3FWhkiY/P2BpUtZkq5bV0rTWgkPD0eVKlVQt25dvUNxC6+9BuzcCTz7rG1DcVYmEz9cJGEynqQknk2n2iLSJiNGAI0bs+J+7pw2sQnhSBcuAC+8ADz+ODefNgJDJEwA0KoVM9GrV/nBIexv0KBBOHr0KPbs2aN3KC5v3jxg8mSgdGkgO10fJGEypiNHOCczq2fT3t7A//4H5MkDhIVxDpsQRnXnDtCuHYfgFi82Tkd7wyRMAPDee/yw2LULGDhQ72iEyJ6DB1lRCgjg6qf7N2G1RZ06PEOTBpbGsns3n++aNbN+34IFOT3h9GlOArdY7B2dENpTCujVCzh+HFi+HLhvj3mnZqiECeDQXNmywNSp7M8khJHExKRsifHHH9xwNTuaNuXQ3O+/2y004QC//cbVkLlzZ+/+1aqxOrl4sUwCF8Y0fjywaBEXcWXnxEFPhkuYPDy4KWXu3MxS/9sqTQinZ7FwKCY2Fpg5E6hdO/vHKlQIeOIJVhyEMdy9C6xdyyp5TrRty+kJY8YAS5bYJTQhHGLFCuD994HRozkfz2gMlzABPCvfupX/3bAhEBWlZzRC2KZ9e+DkSaB/f/Ydyak2bfgFfPduzo8ltLdhA5+rNm1yfqyRI4GXXmJLir/+yvnxhNDakSNA166cu2SUSd73M2TCBLCU98MP3JMpNFTG84VzGz8eWLaMqzynTrXPMcPC+AW8caN9jie0tXIlUK4cezDllMnEz7+KFVlxunYt58cUQis3bvDzqnRpDsVlZ96mMzBo2NSjByd/nz6d8zK3EFpZs4YLFgoVAv78037HrVSJX8DLltnvmEIbZjMTprAwJjv2EBDA5/7uXaBDB2m3IpxTUhKroTExHJLL7vw9Z2DohAkAwsM5iXLVKs5pEsKZbNsGPPccl4Xv3Wvf5bMmEz+I/vc/WWbu7Nas4YrGl16y73FLluQ8pp07gU6d+OUkhLNITuYw3NatwC+/AI88ondEOWP4hAngWXvlysCsWWwGKIQz2LsXePJJJjbbtvHLzd769mWytGCB/Y8t7GfKFE7yf+wx+x+7QQOumlu5EujendUsIfRmsbCIsWQJT+qaNNE7opxziYTJw4O9bcqVYzPAd9/VOyLh7g4f5heZUsDmzVnv7GyrUqWA1q1ZaVVKm8cQOXPmDNsJDBxov+G4+7VuDSxcCPz8M9C7t8zpFPqyWIB+/YD583lp107viOzDJRImAPDy4iz8kiWBCROMOwtfGN+JE0yQkpM5FGPtu6SVgQN5wrB9u7aPI7Lnu++491/nzto+Tvv2nFA7Zw5fE5JACz0oBbzxBvD99xz16dhR74jsx2USJoDzQ44dA0JC2NRt5Ei9IzIW2Usu5w4f5mapiYnsYtusmfaP2aIFq6vh4do/lsiau3f5xfHKK9wTU2tduvDxpk1j+wqpNAlHsliAwYM50vPddxwidiUmZftpiGHOV+7c4XLbixeBIUOAr7/WOyJjiYmJQVBQEKKjoxEYGKh3OIaxdy+rScnJTJZat3bcY0+ezLO6w4c5n084hy+/BIYN4zYQ5co57nFnzwZefRV4+WW2H/D0dNxjC/dkNnMY7ocfgBkz+PozmEwHzF0yYQKA+HigShXOH+jXj9musI0kTFm3bRsneCvFYThHVJZSS0jgSULt2tL92VlERwNlygAvvqjP589PP7GxZYcOwI8/cqWmEFpITub+hj/9xCHhl1/WO6JsyTRhcqkhudT8/HhWV7Eiy9Ndu+odkXBV69alrADZvNnxyRLAXb/HjuVeizt3Ov7xxYMmTOCQnF7zKTt35uqkJUvYziA+Xp84hGtLSOBrbdEiLjwwaLJkE5etMFklJ7MT+MGDvP7zTyZTImNSYbLd558Dw4fz7P3PP9nJWy9mMzvgFygAbNqk3YoskbnISG4S/vrr7PKup1WrWGV69FEm1EWL6huPcB2RkcALLwD79jFhattW74hyxH0rTFZeXsD+/VxBEhEBFC/OYTohcsJi4Wtq2DAgf35WM/WeK+/pyS/nLVu4jF3oZ8wYVv2GDdM7EjZO3boV+Pdfrt7cvVvviIQr2LuXr6ezZ4E//jB8smQTl0+YAPZp+uUXVgNu3uQwnXyhiOyKieHE6iVLWLW8eJF7JDmD557jkOCAAYxTON6WLZwGMGYMNwp3BnXq8AuuVCmgcWPOaRIiu376CWjUCChWjK8rLRqyOiO3SJis3nkHWL+eCdRzz3HOhxBZcfgwq5R//80GgXv3OtcQr8kEzJwJ3LoFvP223tG4n7g4djdu2JDLq51JcDDn2HXuzOXe774rXcFF1lgsbNfTpQsXM2zZ4l5DvG6VMAHAU08BJ08CRYoAo0dz6bf0KhG2+OknzhGKi+Oy2Rkz9I4ofY88wgnHM2YAa9fqHY17GTGCe8bNmuWcO7L7+nLZ99dfs+VBmzZczSdEZmJi2LH7s8/4+TJnjnOdLDqCy0/6zkhyMlc2bd/O4ZR9+5ynfK43mfT9oDffBCZOZPPBzZv1n6+UGaXY0PLvv1kVk6dRe1u2sLXExInsieXs1q5lF+YiRbiLfIUKekcknNXJk5yjdPEiV8I984zeEWlCJn1nxMuLvXMGD+Yk8OLFgT179I5KOJvERA6vTJzI+R8XLjh/sgRwaO777zk017evbJOhtatXOczVqJFxNgBv2ZITwE0mzkH5/Xe9IxLOaN06vj6SkoBdu1w2WbKJ2yZMVt9+y9JifDxfFL16ue8QnWyNktavv3IF3LZtQKtWwOnT/H+jKFWKwy+LFgGffqp3NK4rMZHL9uPjudGoMw7FZaR8efbtatgQePZZbqdy+7beUQlnEBvLgkLLlvxu3LULqFRJ76j05bZDcvc7cwZo3pxfioUKcRWdVjvMOzt3H5JLTORY/e+/s7/S5Mms0hjVBx8A48Zxu5Y2bfSOxrUoxSRj1iz2vtJ6o2WtWCxc2ffOO+zj9f33/DwU7mnzZhYPrlzhydagQcY6EcgmGZKzVenSwKlTwKhRwI0bHHZx52qTu7JWlX7/nQlzZKSxkyWAy9vbtuXKliNH9I7GtUydCkyfzq1PjJosAfwyHDAAOHSIDTdbtJBqkzuyVpWaNgVKlGDD59dec4tkySZSYUqHu1eb3LHC5GpVpfvFxgJPPMEVftu2cYm5yJk1a9ieZPBgznFzFVJtck9uWlVKTSpM2SHVJvfiilWl++XOzZVQ8fFsbHntmt4RGdvmzUywn3kG+OILvaOxL6k2uRepKtlO/iQPMXYsl1OWKcM5CsHBnPgmXEN8PCe6tmnDCtO0aVwpaaSJ3VnxyCPAhg3A9ev8Erx+Xe+IjGnrVvZva9QI+Plnrrh1RaVLs9HvlCnAvHlAtWpcMSVcx6ZN3GNw1ixg0iT+f9myekflvCRhysT91aZ69YD69WU/OiOzWIChQ9mbyJWrSumpVIlJ06VLPKO8ckXviIxlwwbg6ae5amjZMtdv3Hd/tallS55gHDqkd2QiJ44cYYX0qaekqpQV8uex0dixTJIaNOAy3LJlWZ24eVPvyERWTJgA5MkDfPUVEBTEJmyuXFVKT7VqbLJ44wb3FTt5Uu+IjGHxYs5ZatwYWLWKTUzdhbXa9NNPwLFjQI0a7Dt19qzekYmsOH8eeOUVVpUOHmTlUKpKtpOEKQtKlgT+/BPYvx+oUoXViUKFgB49OLwjnNe8eZzA+u67bNT39decx9Oxo96R6aNyZQ4tAayWrF+vbzzOzGLhSsMOHXhWvmwZkCuXzkHpwMMD6NSJCVN4ODuFV6wIDBkic+Kc3fXrwFtvse/Wb78B33wDHD8OdO0qVaUsUUrZehH3Wb9eqZIllQKU8vFR6u23lTKb9Y4q56KjoxUAFR0drXcoObZ2rVIlSqQ8R+++6xrPkb3cuqXU008r5emp1DffKGWx6B2Rc7l9W6kXXlDKZFLq44/l75NabCz/JoGBSuXOrdSYMUrFxOgdlUjt9m2lPvqIz1GePEqNHcufiXRlmgdJwmQHCxYoVaAA/5oBAUp9+aXeEeWMKyRM+/crVaUKnxNPT6V69lTq7l29o3JOyclM9gGlevWSv5PV6dNKPfook4Hly/WOxnldv67U0KFK+foqVaiQUpMmKRUfr3dU7i0hQanJk5UqUoQnim++qdTVq3pH5fQkYXKkL79kwgTwQ3bAAKWMlHNMnjxZVa5cWVWoUMGwCdOcOUqVL8/nwGRSqnVrpW7c0DsqY5gzh196VaootXu33tHox2JRato0npGXKaPU4cN6R2QM588z4fbwUKpYMaU+/FCpS5f0jsq9XL6s1LhxrKqbTEr16KHU2bN6R2UYkjA5mtms1MiRSuXPn/KlXa+eUlu36h2Z7YxWYbpxQ6lXX01JVj09lWrZUj4osuPQIaVCQ/mlN3y4+1UKzp1TqkULvo5691YqKkrviIzn6FGl+vZVyt9fKS8vpV58UalNm2Q4UysWi1JbtijVsSP/3rly8bUriX6WScKkp99+U6pWLf6VAaUKF+Y4f0KC3pE9nFESpvXr+eVuMvHvW7CgUu+/7/x/X2eXlMS5Kd7erDbt2qV3RNozm5WaPp1VpeLFlVq9Wu+IjC8qisNzlSrx/Vm5slLffitJqL1ERysVHq5U1ar8+1aooNTEiZyXKLJFEiZncOUKS6O5cvEv7uWl1DPPKHX8uN6Rpc+ZE6aEBFbwrHPGTCal6tRRasMGvSNzPdZqk8mkVNeuSp06pXdE9mexKLVmjVK1a0tVSSsWi1IbNyrVoQOrvwEBSvXrp9SBA3pHZkwHD3K6R+7c/Hu+8AJPHqWCl2OSMDkTs1mpmTM5L8JadSpWTKn+/ZU6eVLv6FI4W8J09y7PnEJD+QFhnVzfu7fMT9JaUhLn8xQtyorT4MFKRUbqHZV97N6t1FNP8fX0xBNK/fGH3hG5vosXlRo9WqmQEP7d69VT6vPPnffk0VmcOKHUhAl8nQJKBQcr9cEHSl24oHdkLiXTPEg239XJiRPAm2+ygeCdO/xZUBDQpAn39WnRQr/YnGHz3YsX2Vxy2TI2x1OK/UIqVABGjGDTPOE4d+4A337LTTmTktgVeMAA9iYzEqWAHTuAL78EliwBqlYFxo/nViemTLfeFPaSlAQsXw7MmcMeYPHxfG+3aQOEhXGjaFfdcsYWZjNfpytW8HLiBLvKN2vGz77nn+cm4cKuMv0EkITJCezbx+Rg3Trg6lX+zMcHqFkT6NkTePVV/r+j6JUwbdvGXd83bkzpoO7nx61LevUCunVz7w9RZ3DrFpOmqVOBuDh+wQ0cyN3snbkBXmwssGAB90X76y92Nh41Cnj5ZcDTU+/o3NudO9xyZsUKYOVKbteTPz+7qoeFAa1asTu/q4uNZTPQFSu4IfiNG0DhwnyPtWnD91hAgN5RujRJmIzm6lUmDb/8wi0rlOKZb5ky3Mfu6afZbTh3bu1icETCZLHwDGrZMnac/uuvlG7p+fNzn7MhQ4CGDTV5eJFDsbHA/Pns+HzoEFCuHNCnD898y5fXOzqynqUvWgTMncuYW7dmgteihXMneO7KYuFWRStXMnE4dIgni02bct+z0FCgdm0gXz69I825qCieLEdE8CRx40ZuAl6tWkql7bHH5HXqQJIwGdmHH47H3LkmnDnzJJSqASBlPwZ/fw6HhIamTaKUUhg9ejRmzJiBqKgoNGjQAFOnTkX5LHyL2TthSp0cbdsG/P03K0jWl56HB/8t7dpxmNJowzzuTClg+3ZWbhYvBhISuMFvWBgv9eo5toKT3ll6cDArlH37AqVKOS4WkXNnz6YkTzt38vkFeAIZGsrqsxGSqNTJUUQEsHcvN3UHWDWqVy+lklSmjK6hujNJmIxs9OjRyJs3Ly5evIjvv/8eJ09GYckSDt0dOMB5Pqn3sPP3BwICruPWrU1o3746atbMhdWrv8TZs6tx/PhB+Nm4tXp2EiaLBbh0iRs6Hj/O6tiBA+knR4UKcS+zRo2AF17g0KMwvrg4zkexDq1cu8b9++rV45ea9VK0qH3mC5nNnNuxd2/KF9GePSln6WFh/AKSs3TXYLHw8yR10rF//4NJVGgoK54hIXythYQAvr7ax5eQAFy+zMulS8Dp0+knR7VqpSR6oaGcuyXDwk5BEiZXMHv2bAwZMgRRUVEP/O76daRKohROnkwA4Iu0z72CyaSQK5cH8uThmVhwMFCsGM+4c+dm2dt6sVjuYODAfpg0aTo8PXMhMZFfQvHx/DA4dw6IjORjR0dzDkJy8oNxS3LkvsxmYPdublC9Zw+/OKwbtBYpwtdESEjaL7UiRThnzdubr53kZF5u3075IrJ+Gf37L3DkCJM0gF+QdeoA9etz2E3O0t1D6iTKmjinTqKs8udP+1pL/d+BgZwbmfri6cnXsPU1mPq1eOlS2tei9do679LKmhylroRJcuTUJGFyBQ9LmFI7ffo0ypYtix079sPLqyaOHuUZ+IwZa+HtXQaBgeVw8ybf9AkJ/LDJLh8ffiAEBTEpKlqUyVe5ckCVKkD16pywKATACuPFi2nPuFN/6dz/BZeeggXTftFVrcovoVq1gLx5Nf8nCINQislL6oTm/uTGep2QkPXj+/qmfR3ef2397/z5ZeWlwWT6bMmaIxcSGRkJAChVqghCQnhWAwD//DMTJpMJixYtSnP7xESenX399Uz88MOPYGXKF4AfAE/Mm/cD8uXzh68vz/x9fYFHHuEXlxBZYTIBJUrw0q7dg7+/fZuroxITeSZvNqec7QcEsCLqyJWiwrhMJg4FFyjAodmMKMW5RXFxbHNwfzXp/qqTtzdfi3nzSiLkriRhcrDhw4fjs88+e+htjh07hkqVKmkei48PP1CmTOmGr79+6d7PY2JiUKJECbRpMx06tWESbiZPHvdYOi6ch8nE6QnOPFlcOBdJmBxs6NCh6Nmz50NvUyabEzCCg4MBAFeuXEFISMi9n1+5cgU1HzJ5yNfXF76OmBUphBBCGJQkTA5WqFAhFCpUSJNjly5dGsHBwdiwYcO9BCkmJga7du3CgAEDNHlMIYQQwh1kZdK3cDCTyVQSQH4AYQDeAdDov1+dVErF/neb4wBGKKWW/vf/wwAMB9ADwBkAHwF4FEAVpVQ8bGAymQIBRAMIUkrF2O9fJIQQQhiTVJic21gw8bHa/991UwCb//vvigCCUt3mcwABAKYDyAvgTwBP25os/ef2f8e8neWIhRBCCBckFSYhhBBCiExI/1shhBBCiExIwiSEEEIIkQlJmIQQQgghMiEJkxBCCCFEJiRhEkIIIYTIhCRMQgghhBCZkIRJCCGEECITkjAJIYQQQmRCEiYhhBBCiExIwiSEEEIIkQlJmIQQQgghMvF/s7qsY0r/UvUAAAAASUVORK5CYII=\n", "text/plain": [ - "Graphics object consisting of 8 graphics primitives" + "(-1, 4, 4, 6)" ] }, "execution_count": 15, @@ -398,7 +373,79 @@ } ], "source": [ - "draw_circles_from_coords(W) + draw_circles_from_coords(S4 * W)" + "S2 * root" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(-1, 2, 2, 6)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "S3 * root" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For some reason, I have no idea why, $s_4$ doesn't seem to work. Instead this weird matrix does appear to work." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(-1, 2, 4, 8)" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "S4 * root" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(-1, 2, 4, 14)" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "matrix([\n", + " [1, 0, 0, 0],\n", + " [0, 1, 0, 0],\n", + " [0, 0, 1, 0],\n", + " [4, 4, 4, -1]\n", + "]) * root" ] }, { @@ -412,7 +459,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 19, "metadata": {}, "outputs": [], "source": [ @@ -436,7 +483,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 20, "metadata": {}, "outputs": [ { @@ -494,14 +541,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Futile attempt to derive the quadratic form for the cubic packing\n", + "# The quadratic form for the cubical packing\n", "\n", - "We attempted to derive the quadratic form for the cubic packing based on the data collected above, but the results don't match the results from the Stange at all and are significantly uglier. There's probably a mistake in here somewhere, but I'm not sure where." + "We first made a matrix, $W_c$, whose rows are the abbc coordinates of the circles in the root octuple. Then we found the row echelon form of that matrix, resulting in a system of linear relations the coordinates must satisfy. Then, from there, we could derive the rest of the coordinates from the first four (we chose the first four to be the \"cubicle\" from the Stange), allowing us to derive the quadratic form." ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 21, "metadata": {}, "outputs": [], "source": [ @@ -509,76 +556,98 @@ " return [(h1^2 + h2^2 - 1) / b, b, h1, h2]" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The idea here is that by multiplying $W_c$ by an arbitrary vector in $\\mathbf{R}^8$ and setting that equal to $\\vec{0}$, we can find a basis for the null space of $W_c$, which will end up giving us a bunch of linear relations the curvatures must satisfy. The thing to notice is that $W_c\\vec{v} = 0$ for some $\\vec{v}\\in\\mathbf{R}^8$ gives the system of equations with coefficient matrix $W^T$. So finding the row echelon form of $W^T$ will give us the coefficients in the simplified system of linear relations, with each row equal to 0." + ] + }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[ 1 -3 -3 -5/4 -1 -1 -5 3/4]\n", - "[ -3 1 -3 1/2 -1 -1 -1 -3/2]\n", - "[ -3 -3 1 -5/4 -5 -1 -1 3/4]\n", - "[-5/4 1/2 -5/4 1 3/4 -3/2 3/4 -1]\n", - "[ -1 -1 -5 3/4 1 -3 -3 -5/4]\n", - "[ -1 -1 -1 -3/2 -3 1 -3 1/2]\n", - "[ -5 -1 -1 3/4 -3 -3 1 -5/4]\n", - "[ 3/4 -3/2 3/4 -1 -5/4 1/2 -5/4 1]" + "[ 1 0 0 0 1/2 1/2 1/2 -1/2]\n", + "[ 0 1 0 0 1/2 1/2 -1/2 1/2]\n", + "[ 0 0 1 0 1/2 -1/2 1/2 1/2]\n", + "[ 0 0 0 1 -1/2 1/2 1/2 1/2]" ] }, - "execution_count": 19, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Wc = matrix([\n", - " abbc_coords(sqrt(2), -1, sqrt(2)),\n", - " abbc_coords(2*sqrt(2), 1, 0),\n", - " abbc_coords(sqrt(2), -1, -sqrt(2)),\n", - " [sqrt(2) / 4, 0, 1, 0],\n", - " abbc_coords(sqrt(2), 1, sqrt(2)),\n", - " abbc_coords(2*sqrt(2), -1, 0),\n", - " abbc_coords(sqrt(2), 1, -sqrt(2)),\n", - " [sqrt(2) / 4, 0, -1, 0]\n", - "\n", + " [4, 0, 0, 1],\n", + " [0, 2, 0, 1],\n", + " [2, 1, -sqrt(2), -1],\n", + " [2, 1, sqrt(2), -1],\n", + " [2, 1, -sqrt(2), 1],\n", + " [2, 1, sqrt(2), 1],\n", + " [4, 0, 0, -1],\n", + " [0, 2, 0, -1],\n", "])\n", "\n", - "Wc * P * Wc.transpose()" + "Wc.transpose().echelon_form()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "From here we can conclude that $$\n", + "\\begin{align*}\n", + " 2b_1 &= -b_5 - b_6 - b_7 + b_8\\\\\n", + " 2b_2 &= -b_5 - b_6 + b_7 - b_8\\\\\n", + " 2b_3 &= -b_5 + b_6 - b_7 - b_8\\\\\n", + " 2b_4 &= b_5 - b_6 - b_7 - b_8\n", + "\\end{align*}\n", + "$$\n", + "This differs slightly from Stange's system of equations because we put the circles in a slightly different order." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In any case, this gives us the tools to derive the full octuple from just four coordinates, letting us use those four coordinates to represent the entire octuple, and finally making $WPW^T$ nonsingular, allowing us to find the quadratic form the same way we did for the Descartes quadratic form and the octahedral quadratic form." ] }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "W = matrix([\n", - " abbc_coords(sqrt(2), -1, sqrt(2)),\n", - " abbc_coords(2*sqrt(2), 1, 0),\n", - " abbc_coords(sqrt(2), -1, -sqrt(2)),\n", - " [sqrt(2) / 4, 0, 1, 0]\n", + " [4, 0, 0, 1],\n", + " [0, 2, 0, 1],\n", + " [2, 1, -sqrt(2), -1],\n", + " [2, 1, sqrt(2), -1],\n", "])" ] }, { "cell_type": "code", - "execution_count": 21, - "metadata": { - "scrolled": true - }, + "execution_count": 24, + "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[ 1 -3 -3 -5/4]\n", - "[ -3 1 -3 1/2]\n", - "[ -3 -3 1 -5/4]\n", - "[-5/4 1/2 -5/4 1]" + "[ 1 -3 -3 -3]\n", + "[-3 1 -3 -3]\n", + "[-3 -3 1 -3]\n", + "[-3 -3 -3 1]" ] }, - "execution_count": 21, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" } @@ -590,23 +659,346 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "97*b1^2 - 304*b1*b2 + 328*b2^2 - 290*b1*b3 - 304*b2*b3 + 97*b3^2 + 32*b1*b_avg - 1088*b2*b_avg + 32*b3*b_avg + 1280*b_avg^2" + "[ 5 -3 -3 -3]\n", + "[-3 5 -3 -3]\n", + "[-3 -3 5 -3]\n", + "[-3 -3 -3 5]" ] }, - "execution_count": 22, + "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "D = W2.transpose() * m.inverse() * W2\n", - "968*factor(simplify(D[1][1]))" + "32 * m.inverse()" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[bt1 b1 h11 h21]\n", + "[bt2 b2 h12 h22]\n", + "[bt3 b3 h13 h23]\n", + "[bt4 b4 h14 h24]" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "W2 = matrix([\n", + " [\n", + " var('bt' + str(i)),\n", + " var('b' + str(i)),\n", + " var('h1' + str(i)),\n", + " var('h2' + str(i)),\n", + " ] for i in range(1, 5)\n", + "])\n", + "W2" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "D = W2.transpose() * m.inverse() * W2" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "5*b1^2 - 6*b1*b2 + 5*b2^2 - 6*b1*b3 - 6*b2*b3 + 5*b3^2 - 6*b1*b4 - 6*b2*b4 - 6*b3*b4 + 5*b4^2" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "32 * simplify(expand(D[1][1]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This spits out the lovely quadratic form $$\n", + " 8(b_1^2 + b_2^2 + b_3^2 + b_4^2) = 3(b_1 + b_2 + b_3 + b_4)^2\n", + ".$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's the generators for the Weyl group. They might or might not be generators for the packing itself, however." + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[\n", + "[-1 6 6 6] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 6 -1 6 6] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 6 6 -1 6] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 6 6 6 -1]\n", + "]" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "S_i = weyl_generators(m, standard_basis(4))\n", + "S_i" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[ 1 -3 -3 -3 -1 -1 -1 -5]\n", + "[-3 1 -3 -3 -1 -1 -5 -1]\n", + "[-3 -3 1 -3 -1 -5 -1 -1]\n", + "[-3 -3 -3 1 -5 -1 -1 -1]\n", + "[-1 -1 -1 -5 1 -3 -3 -3]\n", + "[-1 -1 -5 -1 -3 1 -3 -3]\n", + "[-1 -5 -1 -1 -3 -3 1 -3]\n", + "[-5 -1 -1 -1 -3 -3 -3 1]" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Wc * P * Wc.transpose()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Generalization of this method\n", + "This method seems remarkably general. Given a matrix representing a root unit of a packing, we can find the linear relation between the coordinates, and thus represent the packing with only four coordinates, allowing us to find the quadratic form." + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "def quadform_from_root(root_matrix):\n", + " P = matrix([\n", + " [0, -1/2, 0, 0],\n", + " [-1/2, 0, 0, 0],\n", + " [0, 0, 1, 0],\n", + " [0, 0, 0, 1],\n", + " ])\n", + " \n", + " # step 1: find linear relation between coords\n", + " relation = root_matrix.transpose().rref() * vector([ var('b' + str(i)) for i in range(1, root_matrix.dimensions()[0] + 1)])\n", + " \n", + " # step 2: find matrix of quadratic form\n", + " W = root_matrix[0:4]\n", + " M = W * P * W.transpose()\n", + " \n", + " # step 3: repeat with arbitrary matrix\n", + " W2 = matrix([\n", + " [\n", + " var('bt' + str(i)),\n", + " var('b' + str(i)),\n", + " var('h1' + str(i)),\n", + " var('h2' + str(i)),\n", + " ] for i in range(1, 5)\n", + " ])\n", + " D = factor(simplify(expand(W2.transpose() * M.inverse() * W2)))\n", + " \n", + " return relation, M.inverse(), D[1][1]" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(2*b1 + b5 + b6 + b7 - b8, 2*b2 + b5 + b6 - b7 + b8, 2*b3 + b5 - b6 + b7 + b8, 2*b4 - b5 + b6 + b7 + b8)\n", + "[ 5 -3 -3 -3]\n", + "[-3 5 -3 -3]\n", + "[-3 -3 5 -3]\n", + "[-3 -3 -3 5]\n", + "5*b1^2 - 6*b1*b2 + 5*b2^2 - 6*b1*b3 - 6*b2*b3 + 5*b3^2 - 6*b1*b4 - 6*b2*b4 - 6*b3*b4 + 5*b4^2\n" + ] + }, + { + "data": { + "text/plain": [ + "[\n", + "[-9 6 6 6] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 6 -9 6 6] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 6 6 -9 6] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 6 6 6 -9]\n", + "]" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# cubical\n", + "relation, mat, equation = quadform_from_root(Wc)\n", + "print(2 * relation)\n", + "print(32 * mat)\n", + "print(32 * equation)\n", + "weyl_generators(32 * mat, standard_basis(4))" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(b1 + b5 + b6, b2 - b5, b3 - b6, b4 + b5 + b6)\n", + "[ 1 -2 -2 -1]\n", + "[-2 4 0 -2]\n", + "[-2 0 4 -2]\n", + "[-1 -2 -2 1]\n", + "b1^2 - 4*b1*b2 + 4*b2^2 - 4*b1*b3 + 4*b3^2 - 2*b1*b4 - 4*b2*b4 - 4*b3*b4 + b4^2\n" + ] + }, + { + "data": { + "text/plain": [ + "[\n", + "[-1 4 4 2] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 4 -7 0 4] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 4 0 -7 4] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 2 4 4 -1]\n", + "]" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# octahedral\n", + "relation, mat, equation = quadform_from_root(matrix([\n", + " [2, 0, 0, 1],\n", + " [2, 0, 0, -1],\n", + " [-1, 1, 0, 0],\n", + " [4, 2, 2*sqrt(2), -1],\n", + " [4, 2, 2*sqrt(2), 1],\n", + " [7, 1, 2*sqrt(2), 0],\n", + "]))\n", + "print(relation)\n", + "print(8 * mat)\n", + "print(8 * equation)\n", + "weyl_generators(8 * mat, standard_basis(4))" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(b1, b2, b3, b4)\n", + "[ 1 -1 -1 -1]\n", + "[-1 1 -1 -1]\n", + "[-1 -1 1 -1]\n", + "[-1 -1 -1 1]\n", + "b1^2 - 2*b1*b2 + b2^2 - 2*b1*b3 - 2*b2*b3 + b3^2 - 2*b1*b4 - 2*b2*b4 - 2*b3*b4 + b4^2\n" + ] + }, + { + "data": { + "text/plain": [ + "[\n", + "[-1 2 2 2] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 2 -1 2 2] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 2 2 -1 2] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 2 2 2 -1]\n", + "]" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# tetrahedral\n", + "relation, mat, equation = quadform_from_root(matrix([\n", + " [2, 0, 0, 1],\n", + " [2, 0, 0, -1],\n", + " [-1, 1, 0, 0],\n", + " [3, 1, 2, 0]\n", + "]))\n", + "print(relation)\n", + "print(4 * mat)\n", + "print(4 * equation)\n", + "weyl_generators(4 * mat, standard_basis(4))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I'm not sure if this works or not. I suspect that it is highly dependent on the order of the circles in the root. Interestingly, it looks like the only relation it was able to deduce for the tetrahedral packing is $b_1, b_2, b_3, b_4 = 0$, meaning there is no null space, as we'd expect. The octahedral quadratic form we get is very different, which isn't surprising, I guess, since this is a very different coordinate system, but I'm not sure if it's right. Likewise, for some reason, we get the cubical quadratic form from the root tetrahedral packing, which doesn't make sense. For whatever reason, though, we get the right matrix." ] } ], diff --git a/.ipynb_checkpoints/Homework-checkpoint.ipynb b/.ipynb_checkpoints/Homework-checkpoint.ipynb new file mode 100644 index 0000000..56de4a4 --- /dev/null +++ b/.ipynb_checkpoints/Homework-checkpoint.ipynb @@ -0,0 +1,557 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n" + ], + "text/plain": [ + "Graphics3d Object" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "var('x,y,z')\n", + "implicit_plot3d(x^2 + y^2 - z^2 == 0, (x,-3,3), (y,-3,3), (z,-3,3))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.2", + "language": "sage", + "name": "sagemath" + }, + "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.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/Apollonian Circle Packings.ipynb b/Apollonian Circle Packings.ipynb index 791e551..40bfa46 100644 --- a/Apollonian Circle Packings.ipynb +++ b/Apollonian Circle Packings.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Calculations to find the quadratic form of the octahedron packing" + "# Calculations to find the quadratic form of the octahedral packing" ] }, { @@ -55,7 +55,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "For an Apollonian packing, you can make a matrix, $W$, where the rows are the coordinates of each of the circles in a quadruple. In an octahedral packing, this is impossible, as there are six circles in a unit. You can try creating a $6\\times 4$ matrix, but, later, we end up needing to invert $WPW^T$, the product of a $6\\times 4$, $4\\times 4$, and $4\\times 6$ matrix, which is singular. Fortunately, the sextuples come in three pairs of circles. Each pair consists of circles that aren't tangent to each other. For whatever reason (we aren't sure yet why, but the Guettler and Mallows says this is the case), the average of the coordinates of the circles in a pair is the same across the sextuple. So, we can make a matrix where the first three rows are the coordinates of circles from different pairs, i.e. three mutually tangent circles, and the fourth is the average of the coordinates in a pair. From this, you can recover the coordinates for all the circles in a sextuple.\n", + "For an Apollonian packing, you can make a matrix, $W$, where the rows are the coordinates of each of the circles in a quadruple. In an octahedral packing, this is impossible, as there are six circles in a unit. You can try creating a $6\\times 4$ matrix, but, later, we end up needing to invert $WPW^T$, the product of a $6\\times 4$, $4\\times 4$, and $4\\times 6$ matrix, which is singular. Fortunately, the sextuples come in three pairs of circles. Each pair consists of circles that aren't tangent to each other. It turns out that the average of the coordinates of the circles in a pair is the same across the sextuple. So, we can make a matrix where the first three rows are the coordinates of circles from different pairs, i.e. three mutually tangent circles, and the fourth is the average of the coordinates in a pair. From this, you can recover the coordinates for all the circles in a sextuple.\n", "\n", "Here $W$ is such a matrix computed for the $(0, 0, 1, 1, 2, 2)$ root sextuple." ] @@ -70,7 +70,7 @@ " [2, 0, 0, 1],\n", " [2, 0, 0, -1],\n", " [-1, 1, 0, 0],\n", - " [3, 1, sqrt(2), 0]\n", + " [6, 2, 2*sqrt(2), 0]\n", "])" ] }, @@ -107,10 +107,10 @@ { "data": { "text/plain": [ - "[ 1 -1 -1 -1]\n", - "[-1 1 -1 -1]\n", - "[-1 -1 1 -1]\n", - "[-1 -1 -1 -1]" + "[ 1 -1 -1 -2]\n", + "[-1 1 -1 -2]\n", + "[-1 -1 1 -2]\n", + "[-2 -2 -2 -4]" ] }, "execution_count": 4, @@ -123,6 +123,21 @@ "M" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$M^{-1}$ is the matrix of the quadratic form, although it is nice to normalize it to have ones along the diagonal, which is fine since it is equal to zero.\n", + "$$\n", + "\\left(\\begin{matrix}\n", + " 1 & 0 & 0 & -1/2\\\\\n", + " 0 & 1 & 0 & -1/2\\\\\n", + " 0 & 0 & 1 & -1/2\\\\\n", + " -1/2 & -1/2 & -1/2 & 1/4\n", + "\\end{matrix}\\right)\n", + "$$" + ] + }, { "cell_type": "code", "execution_count": 5, @@ -133,10 +148,10 @@ { "data": { "text/plain": [ - "[ 1/2 0 0 -1/2]\n", - "[ 0 1/2 0 -1/2]\n", - "[ 0 0 1/2 -1/2]\n", - "[-1/2 -1/2 -1/2 1/2]" + "[ 1/2 0 0 -1/4]\n", + "[ 0 1/2 0 -1/4]\n", + "[ 0 0 1/2 -1/4]\n", + "[-1/4 -1/4 -1/4 1/8]" ] }, "execution_count": 5, @@ -220,7 +235,7 @@ { "data": { "text/plain": [ - "[b_avg == b1 + b2 + b3 - sqrt(2*b1*b2 + 2*(b1 + b2)*b3), b_avg == b1 + b2 + b3 + sqrt(2*b1*b2 + 2*(b1 + b2)*b3)]" + "b1^2 + b2^2 + b3^2 - b1*b_avg - b2*b_avg - b3*b_avg + 1/4*b_avg^2" ] }, "execution_count": 9, @@ -229,28 +244,7 @@ } ], "source": [ - "quad = 2 * factor(simplify(D[1][1]))\n", - "solve(quad == 0, b_avg)" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "b1^2 + b2^2 + b3^2 - 2*b1*b_avg - 2*b2*b_avg - 2*b3*b_avg + b_avg^2" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "quad" + "2 * factor(simplify(D[1][1]))" ] }, { @@ -258,32 +252,12 @@ "metadata": {}, "source": [ "So, we end up deriving the quadratic form $$\n", - " b_1^2 + b_2^2 + b_3^2 + b_{\\text{avg}}^2 - 2b_{\\text{avg}}(b_1 + b_2 + b_3) = 0\n", + " b_1^2 + b_2^2 + b_3^2 + b_{\\text{sum}}^2/4 - b_{\\text{sum}}(b_1 + b_2 + b_3) = 0\n", ".$$\n", "\n", "This means that, given three mutually tangent circles with curvatures $b_1,b_2,b_3$, there are two solutions for $b_{\\text{avg}}$, allowing us to derive two new sets of three mutually tangent circles with curvatures $b_1' = 2b_{\\text{avg}} - b_1$ etc." ] }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "1/2*b1*h11 - 1/2*b_avg*h11 - 1/2*b1*h1_avg - 1/2*b2*h1_avg - 1/2*b3*h1_avg + 1/2*b_avg*h1_avg + 1/2*b2*h21 - 1/2*b_avg*h21 + 1/2*b3*h31 - 1/2*b_avg*h31" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "factor(simplify(D[2][1]))" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -295,38 +269,51 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ - "S1 = matrix([\n", - " [-1, 2, 2, 2],\n", - " [0, 1, 0, 0],\n", - " [0, 0, 1, 0],\n", - " [0, 0, 0, 1],\n", - "])\n", - "\n", - "S2 = matrix([\n", - " [1, 0, 0, 0],\n", - " [2, -1, 2, 2],\n", - " [0, 0, 1, 0],\n", - " [0, 0, 0, 1]\n", - "])\n", - "\n", - "S3 = matrix([\n", - " [1, 0, 0, 0],\n", - " [0, 1, 0, 0],\n", - " [2, 2, -1, 2],\n", - " [0, 0, 0, 1]\n", - "])\n", - "\n", - "# this is the only one that differs from the generators of the Apollonian group\n", - "S4 = matrix([\n", - " [1, 0, 0, 0],\n", - " [0, 1, 0, 0],\n", - " [0, 0, 1, 0],\n", - " [2, 2, 2, 3]\n", - "])" + "def weyl_generators(matrix, alphas):\n", + " retval = []\n", + " for alpha in alphas:\n", + " retval.append(identity_matrix(len(alphas)) - 2 * alpha * alpha.transpose() * matrix)\n", + " return retval" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "def standard_basis(dim):\n", + " return [ matrix(dim, 1, [0] * i + [1] + [0] * (dim - i - 1)) for i in range(dim) ]" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[\n", + "[-1 0 0 1] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 0 -1 0 1] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 0 0 -1 1] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 1 1 1 1/2]\n", + "]" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "S_i = weyl_generators(2 * M.inverse(), standard_basis(4))\n", + "S_i" ] }, { @@ -335,49 +322,38 @@ "metadata": {}, "outputs": [], "source": [ - "W = matrix([\n", - " [2, 0, 0, 1],\n", - " [2, 0, 0, -1],\n", - " [-1, 1, 0, 0],\n", - " [3, 1, sqrt(2), 0]\n", - "])" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [], - "source": [ - "def draw_circles_from_coords(coords, color=(0,0,1)):\n", - " drawing = None\n", - " averages = coords[-1]\n", - " avgb = averages[1]\n", - " for row in coords[:-1]:\n", - " b = row[1]\n", - " if not b == 0:\n", - " x = row[2] / b\n", - " y = row[3] / b\n", - " if drawing is not None:\n", - " drawing += circle((x, y), 1 / b, rgbcolor=color)\n", - " else:\n", - " drawing = circle((x, y), 1 / b, rgbcolor=color)\n", - " newb = 2 * avgb - b\n", - " if not newb == 0:\n", - " newx = (2 * averages[2] - row[2]) / newb\n", - " newy = (2 * averages[3] - row[3]) / newb\n", - " if drawing is not None:\n", - " drawing += circle((newx, newy), 1 / newb, rgbcolor=color)\n", - " else:\n", - " drawing = circle((newx, newy), 1 / newb, rgbcolor=color)\n", - " return drawing" + "S1 = S_i[0]\n", + "S2 = S_i[1]\n", + "S3 = S_i[2]\n", + "S4 = S_i[3]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Hmm. Doesn't look like it worked. And my crude code to draw the circles doesn't seem to be the problem either. I calculated a few of these circles by hand using the matrices and it comes out exactly where it's drawn here, so my code isn't the problem. Unless, of course, I computed the generators of the Weyl group wrong." + "We can test this out on the packing on page 3 of the Guettler and Mallows. The root sextuple is (-1, 2, 2, 4, 4 7), so the coordinates would be (-1, 2, 4, 6). After applying $s_1$, it should swap out $s_1$ for its pair, resulting in (7, 2, 4, 6). Likewise for $s_2$ and $s_3$. Then, for $s_4$, it should give the average between the triple (-1, 2, 4) and the other tangent triple, i.e. (-1, 2, 4, 14)." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(7, 2, 4, 6)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "root = vector([-1, 2, 4, 6])\n", + "S1 * root" ] }, { @@ -387,9 +363,8 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAEECAYAAADah2tfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbr0lEQVR4nO3dd3iT1RcH8G+6aaFl0zJlbxkFBVkiywEFBWXJENmgoqgMRRBRVByIFGQoQ0D4KVuUPUQ2BWSjbBDKbksLXcn9/fE1tIWWpm3evHmT83mePNE2eXNo1nnPvfdck1IKQgghhBAiYx56ByCEEEII4ewkYRJCCCGEyIQkTEIIIYQQmZCESQghhBAiE5IwCSGEEEJkQhImIYQQQohMSMIkhBBCCJEJSZiEEEIIITIhCZN4gIkCTSaTSe9YhBBCCGfglYXbSktwNxEdHY2goCBER0frHYoQQgjhCJkWCKTC5OT++OMPtGnTBkWLFoXJZMKyZcsyvc/mzZtRu3Zt+Pr6oly5cpg9e7bmcQohhBCuTBImJxcXF4caNWogPDzcptufOXMGzz33HJo2bYoDBw5gyJAh6N27N9asWaNxpEIIIYTrMmVh810ZktOZyWTC0qVL0a5duwxvM2zYMKxatQqHDx++97NOnTohKioKq1evtulxYmJi7g3JBQYG5jRsIYQQwtllOiSXlTlMwgB27NiB5s2bp/lZq1at8Prr72LvXuDoUeDECeDsWeDiReDaNSAx0YLkZAWz2QSzGUhODgBwGRUr5oa3N+DhAXh6Al5eQFAQULQoUKoUUK4cUKkSUKMGULiwLv9cYUB37/L1d/lyyuXSJeDqVSAhAUhKAsxmvt68vYGAACAkhJeiRVOuS5Tg61KIzCgF3LyZ8lpLfR0bCyQn85KUxGuzOeUzz/o69PICcudO/7WYPz8gS2Rcn1SYDORhFSaLBdi+HWjTZiby5m0JH5+SuHWLHwbx8Qp8mtO+oz08AF9fwGxOQGLiXQBmAJb/LgoFChSGUh6wWHh8sxlITOSHSnp8fAB/fyZVBQsymXr6aSAsDJBClXu6cwf46y8gIiLlcvQoX0tWgYH84ilSBPDzw70k3foldvt2SmKVmJhyv4AAoFYtIDQUqFOH1xUqSBLlrm7dAvbtA/bvT0nIrUnR/a8dgElOSAhff/cnRp6e+O/kMW0iFRPDY928mfZYPj5AcHDaJKpUqZTXZ758DvsziOzLNOWVhMlArAlTWFg7/PknsHw5k6S//+aHhfWpNJks8Pf3QJ48/FDw8rqOgwdXYeTILqhWzRtVqzKZ8fHh7RMSEpCQkHDvcWJiYlCiRIkMh+QsFn4QHTwIHD8OnDwJnD8PREayYhUdzS/K1IlVrlz8AKldG2jVCmjXTpIoV3X2LLByJbBiBbB5M79ovL2B6tVTEptKlVK+XAICbDuuUnydX7oE/PsvX3979zIJO3WKt8mXD3juOSbprVrJa8xVWZMjaxK+dy9w+jR/FxAAlCnzYBUo9XVwMJPz7EpI4Ofd/dWq1Ana6dM8YQWA0qVTXvuhofwczJ8/538HYVeSMLmC8+eByZOBCRO2I3fuOoiL87mXHHl4cDiscmWgcWNgxYqeaNw4LyZOnHjv/rNmzcKQIUNsbhNgrzlMN28CS5YAa9cCBw4AFy4A8fEpv0+dRPXsCbRoke2HEjo7fhyYP59J0sGDTJCefBJo3Rpo0ACoVo3VTK1Yv0A3b34whrZtgS5d5CzfqJKTgZ07eXJoTZCtyVHu3ClVHGulsXx556gyWizAP/+kxBwRwddo6iTKGvcTT/DiJZNk9CQJk1H9+ScwcSKwaVPq8m8y8uZNRq1afmjcGHjhBeDRR9Peb9iwYfjtt99w6NChez/r0qULbt686RSTvh+WRPn4cD5Uz55A794pFTDhnJKSmJxMmQJs3Ohc1R1rlWvlSr6HvL2Bzp2BgQP5BSWc2+3bwJo1fP5WrQJu3GByVLt2SpJhHYL1MNBab2sSlXqIet8+/nud6f3jpjKfhaaUsvUiNJSUpNTs2Uo1bKiUn59SHIBQKn9+s2rR4qaaNeuEAqC++uortX//fnXu3DmllFLDhw9X3bp1u3ec06dPK39/f/XOO++oY8eOqfDwcOXp6alWr15tcyzR0dEKgIqOjrb7vzM9V64oNXKkUhUqKGUy8d9tMilVtqxS77yj1OXLDglD2OjKFaXGjFGqaFE+Vw0aKDV/vlLx8XpHlr7ISKU+/lipEiUY72OP8b2WmKh3ZCK18+eVCg9XqlUrpXx8+FxVr67Ue+8ptWuXUmaz3hFqw2xWavdupUaNUqpGDf67vb2VatlSqW+/VersWb0jdBuZ5kGSMOnoxg2+SSpVUsrDIyVRKF1aqTffVOrCBaU2bdqkwOpemkuPHj2UUkr16NFDNWnSJM1xN23apGrWrKl8fHxUmTJl1KxZs7IUl6MTptSSkpSaMUOpJ55Qytc3JXEsWFCpTp34wSL0ER3N12tAAC/9+il14IDeUdkuKUmp5cv5hQwwIf/pJ9f9InZ2FotSe/cq9cEHStWqxefEy0upZs2U+uYbpU6f1jtCfZw9y0SpRQsmTgATqVGjlNqzh383oQlJmJzR2rVKhYamVFO8vfn/kyYpdfeu3tHpmzDdb9MmpZ5/Xql8+VKSp0KFlBo9WqmEBL2jcw/x8Up99ZVSBQqw+vnuu0z2jezgQaXatOHrqVYtpVavli8iR4mOVmryZKWqVOHfP29epbp0UWrhQqVu3dI7OucSFaXUokVKde2a8hlYuTK/K6Ki9I7O5UjC5Czu3lVqxAh+6VgrSXXrKrVypfOd4TpTwpTa6dNK9eihVK5cKWejTz+t1NGjekfmupYuVapkSaU8PZXq00epixf1jsi+tm7lkCKg1FNPKXX8uN4Rua6DB5Xq35/VSU9Ppdq3V2rNGhkatVVSklLr1in14ov87AsIUKpvX2NVeZ2cJEx6O3RIqebN+QEB8EXep49zn6E7a8JkZTYrNXOmUmXKpFSdypThUJ6zJZ9Gdf26Up0782/bpo1rJxIWC09cypdnBe2LL5RKTtY7KteQkKDUggWcmwkoFRLC6rCrJd6O9u+/Sn34oXHmERqEJEx6MJs5Bl2yZMoXeoUKSv34o96R2cbZE6bUjh9X6tlnecYFsPrUvTsnJovsWbJEqcKFOQQwf777DFXFxXHuoMnEOXSunCRq7dw5LuQoXJjvy6ZNlfr5Z6km2VtiolKLF3Pel3W6wvDhSp05o3dkhiQJkyOZzTx7sg4ZeXsr1a6d8SYvGilhskpK4hmX9QPa+iEtK+xsFxPDuSSAUm3buu/fbutWpcqVY7Vp0iT3SRjt4Z9/uDjDw0OpPHmUGjxYqSNH9I7KPRw7ptTrrysVGMik/8UXlTpxQu+oDEUSJkf57ju+UAGlgoKUGj+eX+JGMnnyZFW5cmVVoUIFwyVMqW3dmrLqxmRSqkMHpW7f1jsq53b6tFLVqvFLbt48SRLi4pR67TW+hnr2lKGOzFy6pNSAAaz0FivG9gDyntNHbCy/j0qU4FSQvn05hCcyJQmT1pYsUSo4mH9JPz9WOYw+j8aIFab0bN/OSoF1gviAAcZLYh1h40YuRihbVqoB9/vxR7a3qF/ffStuDxMVxaE3f38O4X7+uVJ37ugdlVCKC42+/FKp/Pk56jFsmFI3b+odlVOThEkrf/7JLxjrl/HAga7zZewqCZPVsmVpk9oxY4yf1NpLeDjPQps1c+6FCHratYuTlYsVYx8cwS/jCRNSvoxHjJCWAM4qKorNP/392cLhs88kqc2AJEz2dvx42uGeF190vdKzqyVMVtOmcbjUOmw6bZreEenHYmEXdYDzHlwl2dfKv/+yDUiuXOzZ5K6SkpT6/nulihdnot2/vwz3GMXlyzyx9/Li6rrp0+V9fx9JmOzl7l2lnnkmZULxU0+5boneVRMmpVhZ+vDDlO1nQkLcr3u42czJuIBSEyfqHY1x3Lmj1HPPcduO5cv1jsbxNmxg00RAqY4dlfr7b70jEtlx8mRKy5AKFdhIWSilbMiDDLRtoX5WrAAKFgR+/52bwx4/DmzYAAQH6x2ZyCoPD+CDD7jZ5YABwJUrwGOPAa+8wo0xXZ1SwGuvAZMnA999B7zxht4RGUeuXNw4unVroH17bgzrDmJjuWlxs2ZAoULA3r3AwoVA+fJ6Ryayo2xZYMECYP9+oFgxoGVLoF8/ICZG78icnyRMDxEfDzz9NNC2LXdmnzEDOHAAqFhR78hETnl5AVOmACdPAuXKAbNnA0WKAHv26B2ZdpQC3n6b/+7p0/khKbLGxwdYtIg7ynfoAKxZo3dE2tq4EaheHZgzh0n2pk1AaKjeUQl7qFmTJ/5TpzKBql4dWLdO76icmyRMGbBWldasAerWBS5fBnr31jsqYW+lSwP//AOMHg3cvMlqU69erlltmjAB+Oor4NtvgT599I7GuLy8gJ9+Alq0AJ5/HoiI0Dsi+4uNBQYNYlWpVCng0CH+v4d8Y7gUkwno35/Pb7lyUm3KjEkpZettbb6hkcXH80Nw9WqeTYaHu1+iFBMTg6CgIERHRyMwMFDvcBzmzBl+CZ46xaGHVauYLLuCVauANm2A4cOBTz7ROxrXcPcu0LgxEBnJyqSrDNFv2sSThqtXgc8+43CcJEquz2IBpk0D3nkHKFAA+P57oHlzvaNyKFNmN5C3QSrWqtLq1VJVckelS3OIbvRo4MYN16k2HTsGdOnCuTfjxukdjevIlQtYtgxITuacpoQEvSPKGWtV6amnWFU6eBAYPFiSJXfh4cF5nYcPs9rUooVUm+4nFSbwCzEsjGfh7lpVSs1dK0yp3V9t2rgRqFZN76iy7tYtJn6+vsD27YCbPp2a2rULaNKESen333OYw2j+/BPo1k2qSoLurzbNncvXuIuTClNmrl7l2dSqVUCdOlJVEnR/talmTc5bMRKl+CV+8yawfLkkS1p5/HFOop81ixPqjWb6dKBpU6B4cakqCUpdbSpdmnPZpkzhZ4o7c+u3xY4dwCOPABcvMpPeswfIn1/vqPQTHh6OKlWqoK6rTNyxgzFjgJ07WaHp0gV48029I7LdtGkcXl6wgEuJhXa6d+dw1jvvcBGBESQlMeZ+/XjZuFFeJyKtRx4B1q9nEj1oECeIJybqHZV+3HZIbsqUlDOphQu5RFiQDMk96OZNoHZt4Nw5oEEDYPNmrpZyVmfPcplwly5MnIT24uLYpy0kBNiyxbmrNNevAy++yKG48HCgb1+9IxLO7ocfmDDVqwcsXsypCi5GhuTS06MHs+XAQC6nlGRJZCZ/fuD0afbl2rYNKFECuHRJ76jSpxTw6quMecIEvaNxHwEB/FL5809g0iS9o8nYoUNc1HL4MPvwSLIkbNGrF1dQnjjB189ff+kdkeO5VcIUH88zwLlzgSpV+IVXubLeUQmj8PBgt/dRo7iUvEwZVpqczbRpHF6ZOVPmLTla48bA668DI0c659DcsmVA/fp8Xezdy3iFsFWDBnzdFCgAPPEEK03uxG0Spn/+AYoW5aTGTp2AI0cAf3+9oxJGNHYsJ1FbLFyC/eWXekeU4upV4N132ZiyRQu9o3FPn3zCz5pBg/SOJIVSwEcfscfc009zxWSpUnpHJYyoRAlg61b2devQgfM8jd56xVZukTD9+isrSlFRwMSJxlvtJJxPWBj7G+XLx+1GOnXSOyL6+GMuax8/Xu9I3FdAAPDFF9xmYsMGvaNhg82OHbmH4ocfAv/7H2MUIrv8/fk9+vHHPIHs0AG4c0fvqLTn8pO+Fy0COndmf6W1a6UEbQuZ9G27+HhOgvzrL3bF1XMvpjNnuM/hmDEcEhL6UYpDFsnJwO7d+vVmiotjJWDnTmDePOCFF/SJQ7iuFSu4uCQ0lMWJPHn0jijb3HvS95w5TJZy5eIQnCRLwt78/IB9+zg0t3490KiRfuXp0aM5t+CNN/R5fJHCZAI+/ZTzPfSa5xETA7RqxXYpa9ZIsiS0ERbGYsSBA3y9RUfrHZF2XDZhmjoV6NmTpedjx6S/iNCOhweHXp55hiukHnvM8UnTwYOsIHzwgQy3OIsmTThf6L33WGlypFu3OIft8GFWPRs1cuzjC/fyxBP8DDx+nE0ub9zQOyJtuGTCNHEiW/sHBXEJZMmSekck3MFvv/EsPiKCncEdmTR98gk78kqXeucyfjzw99+cN+QoN26w4nnyJFdL1qvnuMcW7qtOHbYdOHeOneOvXdM7IvtzuYRp6lR2Y86fnx8YRYvqHZFwJ4sXA127stdNnTqOSZoiI/m4r78OeHtr/3jCdjVrMnlx1JYpUVGsLP37L1te1K7tmMcVAmDbni1buFq3RQs2/HUlLpUwzZrFylLevKwsFSyod0TGIluj2Me8eeyivH8/+5ZonTTNnMlEqUcPbR9HZM/AgWx2qnWjv9u3OQR49iyH4apX1/bxhEhPlSqcz3nxIl+PrjSnyWVWyVlXw+XOzRJ4cLDeERmXrJKzj7ZtuYKkaVMOjWghOZlDcU8/DcyYoc1jiJxJSuKeXG3aAN99p81jxMVxDt1ff3EuSZ062jyOELbav5/V1SpVuOggd269I8qUe6ySW7EiZTXc4cOSLAnnsHw5y9KbNgGtW2vzGL/+yjO5AQO0Ob7IOW9vbj8yb542Z9sJCUC7dlytuXq1JEvCOdSqxUTp0CGeLNy9q3dEOWf4hOnECaB9e+4mf/CgTPAWzmXtWg7LrVrFDtz29v33XJUnc1WcW58+7Nm1aJF9j6sUk+WtW/kaq1/fvscXIicee4zbSe3axZMG2we0nJOhE6Y7d7gCxGzmmKm0DhDO6I8/gOLFuRHuL7/Y77hxcZyr0rGj/Y4ptFG0KNCwIavh9jRpEuduzpjBNgZCOJsGDbgp9bx5zrWNVHYYOmGqV4+rQiZN4pMihDPy8GCrgVy5uIXKsWP2Oe66dRyOadPGPscT2goL44ldXJx9jrd2LfDWW8A77wDdutnnmEJooVMn7j7w7rtsv2JUhk2YrEu3u3UDBg/WOxohHq5wYU7GtVg4bBIbm/NjrlgBVKoElC+f82MJ7bVpwwTXHtvn/PMPK4utWsm+gcIYPvqIczk7d2aDSyMyZML07bfAggVcNjt3rt7RCGGb+vWB8HBO/H388Zwdy2zmhO+wMPvEJrRXvjwT3JUrc3ac6Gg+70WKcANUT0/7xCeEljw8OCxXvDhfv7du6R1R1hkuYdq2jXtl5c3LDSWFMJIBA9gv6ehRnmll15497KQrw3HGEhbGRDe7k1/NZlbXL19mhTEoyL7xCaGlwEC+bq9f5+efo7cMyilDJUyRkdwR3tOTO4D7++sdkRBZN3s2O+IuXMhtfLJjxw6uDM1ppUo4VpMm7IJ8+nT27v/ee1x1tGgRUKGCfWMTwhHKlgV+/pnz+YYN0zuarDFMwpSczP4i8fHcl0nmbQgj27mT2/e89RZX0WWVdb862QrFWEJDeR0RkfX7zp8PfPYZ8MUXnLskhFE1a8aTxa++4gmkURgmYWrdmvsjjRgBPP+83tG4JtkaxXH8/Fgl9fTkl19MTNbuHxGR8uUrjKNIEc7hyGrCdPAg8OqrHM4dMkST0IRwqEGDuFl4v35sumoEhtga5ZdfuDdX/frA9u16ReE+ZGsUx8nOa/v2bc5dmTkT6NVL2/iE/bVrx1WS69fbdvukJDYANJuZZPv5aRqeEA6TmMj2QElJwN69nGagI+NvjXLnDtC9O/+Qq1frHY0Q9tWhA+fl7dgBzJlj230OHOCkYakwGVNoKCtMtp6rfvIJW6jMni3JknAtPj58XZ84wbYDzs7pE6awMO5BM3MmZ9gL4WpWruQChn79bBuaO3qUS3SrVNE+NmF/jz7KhruXL2d+27/+AsaNY9M/2f5GuKJHHwVGjQI+/TR7c/scyakTpl9+YbO/Bg2Al1/WOxohtOHnx/4kCQnA009nfvtLlzgXRiZ8G1OxYrzOLGFKSgJ69gQqVwbef1/zsITQzfDhTJx69uTnoLNy2oQp9VCckVupC2GL55+3fWju8mXuTSaMyfrcZZYwpR6K8/HRPCwhdOPtbYyhOadNmNq0kaE44V5SD81FRWV8u8uXgZAQh4Ul7KxwYcBkenjCdOCADMUJ95J6aG7vXr2jSZ9TJky//AJs3ChDccK9pB6ae/bZjG936ZJUmIzMy4tDqpcupf97GYoT7so6NPfKK845NOd0CZMMxT0oPDwcjzzyCPz8/PD4449j9+7dGd529uzZMJlMaS5+srTGMGwZmouMlAqT0YWEZFxh+uQT4PBhGYoT7sfZh+acLmFq21aG4lJbtGgR3nrrLYwePRr79u1DjRo10KpVK1y9ejXD+wQGBuLy5cv3LufOnXNgxCKnUg/NxcY++Pu4OCB3bsfHJewnd26eHN7v4EEZihPuLfXQnLM1tHSqhOnAATZzq1tXhuKsvvrqK/Tp0wevvPIKqlSpgu+++w7+/v744YcfMryPyWRCcHDwvUuRIkUcGLHIKT8/VpcSEliavl9yMod1hHF5eaW/8ejbbwPlyslQnHBvw4dzSHro0OxvVK0Fp0qYXn6ZkyF//lnvSJxDYmIiIiIi0Lx583s/8/DwQPPmzbFjx44M7xcbG4tSpUqhRIkSaNu2LY4cOfLQx0lISEBMTEyai9BXhw784ly8mENwqUnCZHzpJUwbNgDr1nFITobihDvz9gbGjwc2bwbWrNE7mhROkzBt2AAcOQI88wxQqpTe0TiH69evw2w2P1AhKlKkCCLv/xb9T8WKFfHDDz9g+fLlmDdvHiwWC5544glcvHgxw8cZP348goKC7l1KlChh13+HyJ65c3l21bVr2p97egIWiz4xCfuwWPg8WinFs+p69bh1ihDu7rnngIYNuX+ss3zeOU3C1KsXuxf/+KPekRhb/fr10b17d9SsWRNNmjTBkiVLUKhQIUybNi3D+4wYMQLR0dH3LhcuXHBgxCIj9etzG42NG4Fjx1J+7uXFlVTCuJKS0lYJFy/mUupPP2WVXQh3ZzLx/XDgALBokd7RkFMkTPPnA+fPc0guf369o3EeBQsWhKenJ65cuZLm51euXEFwcLBNx/D29katWrVw8uTJDG/j6+uLwMDANBfhHObP53XqKpOvLxdGCOOKj0/ZaDQ5GXjvPXZ5b9JE37iEcCYNGrAn4/vvc6NevTlFwvTGGxyzf0gRxC35+PggNDQUGzZsuPczi8WCDRs2oH79+jYdw2w249ChQwiRdeiGVLEi0KwZsH8/sG0bf1akCHBfDi0MJjKSzyMAzJoF/P0352wIIdL65BPgzBmunNeb7gnTV18BN24Ar78uO3Gn56233sKMGTMwZ84cHDt2DAMGDEBcXBxe+W/5VPfu3TFixIh7tx87dizWrl2L06dPY9++fXj55Zdx7tw59O7dW69/gsihefM4XN2jB///YT18hPNTKqVb+507wJgxQOfOQM2aekcmhPOpVg3o1g0YOzb9NiuOpOtaG4sF+OADICAA+OwzPSNxXh07dsS1a9fwwQcfIDIyEjVr1sTq1avvTQQ/f/48PDxS8t5bt26hT58+iIyMRL58+RAaGort27ejimxtb1jBwUD79lw9unQpu3z//bfeUYnsunGDc5iKFgW+/Ra4etU5m/QJ4SzGjgUWLgQmTtS35YZJ2d7kwO7dEN59F5gwAfjyS+Ctt+x9dJFdMTExCAoKQnR0tMxnchKxsUC+fEDBgtw2Y+FClqmF8Rw6xOZ8q1cDnTpxftrkyXpHJYRze/NN4IcfgFOn+DmogUyXW+g2JJecDHzzDVCggCRLQmQmd26gb1/OfTl7lvuQOVNDN2E76x5yK1ZwIqs0qRQicyNH8jPv88/1i0G3hGnCBH5YjB2rVwRCGMs337Ch28aNfO9IhcmYjh3jIpdFi5gE27jgVQi3VqgQ0L8/MGNG+tsKOYJuCdOkSZzk3b+/XhEIYSxeXmzsat1GMCJC33hE9kREAMWLcy7TgAF6RyOEcfTvD0RHc0qCHnRJmHbs4NDCCy9w9Y8QwjZffcVrX182OhTGs3cvN1Bu0QKoUEHvaIQwjjJleNIYHq7PlARd0pV33uH1l1/q8egiI+Hh4ahSpQrq1q2rdygiA2XLAuXLc0huzx69oxFZdfs2cPw4+2gNHKh3NEIYz8CBwL59+nz+OXyVXGwsEBjI3goHD9rjiMLeZJWcc5s7lz2ZrB2/ZSsN49i6FWjcmE0rL16UTZSFyCqzmRuTN2kCzJ5t10M73yo560z3ceMc/chCuIbu3ZksJSSk3WNOOD9r0/4BAyRZEiI7PD05l2nhQs4DdCSHJ0xz5gB58wJhYY5+ZCFcR6dOvJ44UdcwRBbNncuKYL9+ekcihHH16sXCy6xZjn1chyZMixcDMTHAf7t6CCGyyTr5e8ECfeMQtrt6la0g6tSRVgJC5EShQsBLLwFTp3LHEEdxaMI0ejRXxclwnBA5kz8/8MgjXG0l26QYg3WRy3vv6RuHEK5g0CDg9Glg7VrHPabDEqZLl4AjR4BGjQB/f0c9qhCu64sveN2nj75xCNvMmwfkyiXTEYSwh8cfB2rVAqZMcdxjOixhsm6u+/HHjnpEIVxb+/acALlzp96RiMxcucKTxqZNZVWjEPZgMnEu0++/s5mlIzgsYVq+nGdXDRo46hGFcH1Vq7In0759ekciHmbMGF4PG6ZrGEK4lLAw7ku7erVjHs8hCdOdO8D584D0QxTCvqzzAQcP1jcOkTGluG9cnjzswSSEsI+SJYGaNbmRtSM4JGH67jt+aMhcCyHsq00bLqTYtQtIStI7GpGebduAW7eAZ5/VOxIhXE9YGPDbb475/HNIwjRvHj/Uu3RxxKOJ7JKtUYypenUurf3lF70jEen56CNev/66vnEI4YratAGiooA//9T+sTTfGsViYVfi8uWBo0ezcwThaLI1irHMnMnqbdmywD//yKRiZxIZCRQrxpXB0dGy2bgQ9maxACVKAB07pvSnyyb9t0b59VdOyurYUetHEsI99ezJJOnUKcf2JBGZs1aX2rWTZEkILXh4sMq0YgWn/mj6WNoePqVHwhtvaP1IQrgnLy9uRmkyAcOHO7bzrcjYqVPAtGl8Pjp00DsaIVxXWBjfb1rvral5wrRtGxASwv3jhBDa6NCBZ1cHDnBFltDfqFFAQADg5wc0b653NEK4rqee4rC31qvlNE2YDh8GYmOBp5/W8lGEEEOG8LpIEeD999mbSehn/37gp5+4Z1zz5kychBDa8PMDWrY0eMI0aRKvhw7V8lGEEIUL83L3Ljd4nTxZ74jcl1LA229zmPTkSc6vEEJoq00b7npw44Z2j6FpwrRzJ+Djw27EQght1aoFxMQAAweyyvTPP3pH5J6mTwc2bgT69uX8pYYN9Y5ICNfXsCFPViIitHsMTROmM2c4f0kIob0nn+R1ixZ83/XqJRPAHe3sWVaX+vThJPyAAKBiRb2jEsL1lSvHbvqGTJji4zl/qUYNrR5BCJFa+/a8Xr8e+OEHNnL79lt9Y3InSgG9ewP58wNffMEP7po1uUGyEEJbHh5A7doGTZhWreK1rA4RwnZKATdvcnf7O3eydt/y5dliYMcOoEkT4LXXgBEjZGjOUaZNAzZsYCPRwEB+cIeG6h2VEO6jTh2DJ0zWs14hRPoiI4GPPwaaNgXy5QMKFODqqoAAJkGdO3P1h9mc+bGKFOFEYwAYP55dpjt0YLVXaCciAnjzTc5batGCXb3/+UcSJiEcKTSUw+JaTfzWLGHau5cTvosW1eoRhL3JXnKOdesW57qULMmEKV8+4N13uSfc8uXA7NnAc88BJ04AbdsCZcoA//vfw49ZtSq/rJOTmXAtXQqcPs1u4DKfSRuRkXx+Hn0U+OYb/mzfPl5LwiSE41jfb1pVmTRLmGTCt/EMGjQIR48exZ49e/QOxeVt3AhUq8YE6NNPgX//BZYsAUaOZFU2LAzo0QOYOJFfvnv28MOgY0fgxRe5Gi49TZvyesMGXlerxs2vFy9O2aZD2E9CAvD880xGly5lPxiAH9j+/kClSvrGJ4Q70XrityYJk0z4FiJjK1eymWvlysCRI8Bbb7G69DB16jChWrSI+8U1b85K0v2sQ+C//prys7ZtgXHjgDFjmDgJ+1AK6N+fTSqXLUtbTY+IYJsHmfAthONoPfFbk4RJJnwLkb79+1khatMG+P13oHjxrN3/pZeAzZs5P+allx7cbDL1xO/URo7k7bt35+o5kXPjxnHYdOZM4LHH0v5OJnwLoQ8tJ35rkjD9/juvZcK3ECkSEzmXqFIlYMECwNs7e8epVSul0jR9+oO/L1KEG1GmZjIBs2bxi/3ZZ4Hdu7P32IImTAA++IBJ08svp/2dTPgWQj9aTvzWJGHat49fBjLhW4gUU6cCR4+yKuHrm7NjtWzJCePvvPPg0FzlyukP1/n7cziwWjWgVStJmrLryy85Of+993i536FDvK5Vy7FxCSFS3ncHD9r/2JokTNeuceKVEIIsFiA8nMNxNWva55hjxrBX048/pv15uXIcqrt+/cH75M7NCnDlyhwyl+G5rBk3jp28hw/PeBL9v//yulQpx8UlhKCSJXltfR/akyYJU3Q0kDevFkcWwpi2b+cwTf/+9jtm0aJAu3YcakutTBleZ3SGFRQErFnDyZEtW7KNgXi4pCQ2Ah01Chg7FvjkEw5zpufSJVbz5KRRCMfz9+dn3KVL9j+2JgnT3btAoUJaHFkIY9q5k2/kBg3se9ynngIOHOB7zqpyZV4fP57x/fLkYaWpbVtWvUaPlj5NGblxg6sav/uOw6qjRmWcLAHA5ctMZh92GyGEdooW5fvQ3rzsfUCLhU3zZP6SECn27cv5MnOlgG3b2Bpg714mStYO3oULc3VInTpA48b8WWZbouTKxcnn1asD77/PuTdz53LYTtCRI+yJFR3NPfqaNMn8PpcuSQ86IfQUEmKQCtOFC7yW8XshUty4we1OskMpYOFCzn1q1IgJU7FirHR89RVv88IL3FJl/nx+wQOsamXGZGLLgWXLgHXrgCeeAP7+O3txuprFi4F69dgxfc8e25IlgGe2kjAJoZ+QEG0qTHZPmKwrRMqXt/eRhdZkaxTtmEwP9kyyRWQkO0l37gyUKMFWAmfPMoF6910OpwFAp06ci3ThArdVMZmYMHXrxs18MxMWxtvHx7Ph7Fdf2bZ3nSu6cQPo0oV78LVqxflnpUvbfn/rkJwQQh9FixqkwnTsGK9lSwDjka1RtBMSApw/n7X7HD/OIbYdO9jl+9dfubGrR6p37blzvLZWrzw9mfwEBHAe4a+/skpiy2NXrcphvv79uRKscWP3qzYtXQpUqQKsXs0tZX7+OetDlDIkJ4S+rBWm7JykPozdE6bTp3n96KP2PrIQxhUaylVriYm23f7MGU7ozpuXSczzz6d/u4gI9nSqWjXtz/PmZYUoIoIrvJ56CrhyJfPH9fcHvv4a2LIFuHqV1aYJE7hnmiu7dIlVpRdeYIJ55AjQtWvWJ27fucP5TlJhEkI/RYvyvXj7tn2Pa/eEyXrGW7CgvY8shHE1bMhkac2azG9rNvPL28+Pm+g+rFqxciXw+OOAj0/anxcoAMTFscXAxo2cHN6rl+1nXI0aAX/9xWrT8OGsGP/4o+sN00VFASNGsHfVmjWsKi1blv0KkXXehFSYhNCP9f1n72E5uydMkZEPfngL4e5q1uTw2pQpmd/266+BXbuYoBQpkvHtTpzgyq0+fR78XUhISjWrdGlgxgzgt9+AOXNsj9labTp0iPF3786VfqtW2b/U7Wh377JyVqYM8M03wJtvcjuZ7FSVUrMmTFJhEkI/1vefvSd+2z1hunGDH7RCiLSGDOHcmLVrM75NbCwbIw4a9PCeTUoBQ4cyMbJO/E6teHHeJiaG/9+mDSeOv/ceh+iyokoVzu3ZsQPInx9o3Zqr6RYsMN5Q3ZUrbDpZvjxXB3bqxETp44/t02w3MpLXD0t0hRDass7ptL4f7cXuCVNSUvY3FRXClXXpwu1IevfOeOXa/PkcSnvnnYcfa9YsVnqmTUt/XzrrRGVrnyaAq+ouXeIwXnbUqwds2sRKlb8/qzElSjDxsA7FOyOluAVMly6M96OPuPrt2DFW/Ow5fGZNIOWkUQj95MrFa3uf0Nk9YTKbc9acTwhXZTIBM2dyOOiZZ4Bbtx68zezZwHPPpeyHlJ6VKzm3qHdvVo7SYz1pSd0BvGZNoH59PkZ2mUyMfcMGbiTcuTP3yCtTBnj2WQ79adH/JKuU4iT7ceM4cb1RIzb7/OwzJo3ff895S/Zmrd552b0lsBDCVtaVxFmtpmd6XPsejp2+PTTZcEUI4ytVipOLT55k8rJrV8rvEhPZEbx58/Tvm5TE4brnn2ei9LD5UNaqU3x82p83b87HtMccpMqVOf/n0iVuGxIby0SuaFFORB83jhPHHTVR/O5dNt987TXgkUeYKH3+OeNcu5ZtGt58E8iXT7sYkpN5LSeNQujHZOJJi/X9aC+abI0iHxZCZKx2bW5x0q0b5wK9+ioweDATi8REtiBILSGBTSk//5zL3YcPB8aMeXgVw1phur8kHRrKdgH//st5TvYQEMCJ5336cA7jb78BK1awmjNqFIcHa9XiY1sv5cvnrApz5w4rSBERrBxFRLDiZTazOhcWxqSySZP0hyy1kpzMf5fsIyeEvnRLmEwmkyk6OtqmA5rNuWEyATExsZnfWDiFhIQEJKT6Zr39X/OKGOuMYWF3RYuy0jRlig/Cw30wY4YHgoPNADzx/vtJyJWLb/aoKBP++ccTUVEmNGmSjPXr41G7tgV37jz8+BaLN4BcuHkzFjExKbvqFi7sASA3/v47DoGB9i/9eHtzQ9+2bZms7dzpiX37PHHggCeWLfPExIksP3t4KBQqpFCggEJgoIK/P+DtrWAyPdgV3WLhse7cMSEqyoRr1zwQFWX67/EUKle2oFYtM3r1MuOxx8yoWtVyL2FJSHDsxPTYWG94efkhJsbODWCEEFni5ZUHsbEJiImxrfldUFBQIIDbSj2k/q6UyvQCIBCAsu1ySwGnbbytXOTi7pfcChikgLMKUP9dLKkuSgFJClipgIZZOO7Q/+57/30q/ffzJxz4b/RQQBUFdFPAVAUcVEBUqn9f6oslncv9t0lSwFUFbFXAOAWEKaCoEzyXUMBrCohzgjjkIhd3v0Qpfg5m6X6BD8uFTA9Lpqz+qzBZMr0hgJIlAxAdfREXLngiMDDQlrtkSd26dTXZukOr4xrh2PdXmC5fvozHHnsMR48eRbFixXIa4gPkOaRVq7wwZIgfrl61ICxMoVIlMz7/3A8REbEoV86C+HjgyBEPbNvmhTlzvHHypCdat07CV1/Fo0iRh79vx42zYMKEvFi58goaN8517+e7d3uiRYsA/PlnLKpXt+ktna6H/T2UAo4d88Bvv3lh/Xov/PWXJ+7cYcmnbFkzata0oEIFM4KDFUJCFIoUsSAkRKFgQQVPzwePHR8PXLliwpUrHrh82YQrV0z4918PHDnigf37PXH9OqtWhQtb8NhjZjz9dDJatUpG4cJp/0aOeH3MmOGNkSP9cO2afSpMMTExKFGiBC5cuGD3z1OjvV+0PrZWx9XyOQSM9/dw1LFDQvLggw8SMGCAzRWmIGRSYbJpSO6hJar7eHpaAHggMDBQkxeHp6c2iZhWxzXysfPkyWOov7VR/s5mM/D665y03bo1cPx4CyxfvgGRkd74/HPg9OncqF0bCAwEChcGmjbl0v3Fi4HBg71Rr543li9n9/CM4+Vs7/z5AxAYmLIZ2j//cGw/NDQ3/Pyy/2+4/++RnAz88QfnLq1Ywa1dAgKAli253UidOpzHFBTkCcATQMa9R+4/tvXvUL36g7dVCrh4kXOYIiI8sHmzB157jcd+/HHOZWrblr2kHPH6yJOHfwt7P44Wn6dGeb846thaxgxo8xwCxvx7OOLYyclA7tx+CAy07YNOKZXpHBS7T/rmCjntZn0PGjTIUMc18rG14s7PoVJAv37so/Tdd0DfvsCUKS8AYLO1YsWA3buBDh3S3s/Dgw0qmzbldatWXPmVUXPLhARWdO6f8LxnD1CtGnKULAEpf4/ISLZKmDaNiUvx4ikTrp98MnuPk5W/tcnE3kolSgDt2vFn166lTDz/+GMmm/XrA3XrfoP4+Jz/2x8Ws5cX51wZYbWwEd4vjjy2ET9LAWP+PRxxbOsCDLuyZQ7TfxebFC5sVsBlFR0dbetdhJO5cOGCAqAuXLigdygu57vvOA9n7tz0f//qq0oVL65UUlLGx7hzR6lGjZQKDlbqxo30bzN4cLwClDp4MOV9GBurVFCQUsOHZz9+pZSyWJTavFmpl15SystLqVy5lOrdW6ndu/k7Z3L3rlKLFyvVvDn/7gULKjVsmFKnT2vzeHPm8HESEuxzvOjoaAVAPk8NTJ5DxzOb+T6cOTNLd8s0D7L7OZCXlwm5cuWBryPX8gq7sj538hza17lzwNtvc/l9t27p32bgQFZqHtaNO1cuYOFCzu154430b5OczCpvYGDKc7hgAbdK6dcvu/8CbuT7+OOsHv31F/Dll+zDNGMGULeu8y2n9/PjkOC6dezD9PLLrOyVLQv06GH/DuXWdg72apjn6+uL0aNHy3vRwOQ5dDxrOwGnrzBVqKBUQECWsjrhZOSMSBsDBihVpIhSmf1ZmzRRqlw5peLiHn67GTN4FnX06IO/e/FF/s5s5v/fuMGKVPv22QpdRUQo1bIlj/nYY0qtWeN81SRbxcUpNXkynwsfH6XeeEOpq1ftc+zffuPf6Px5+xxPCJF1ly/zfbh8eZbu5vgKU5EiD3YXFsLdxcQAP/7I6k5mcx2t84FGjnz47bp1AwoVYsXkfpcusYGshwfnTb3+Ot+XkyZlLe7z57lBbWgoqzGLFwM7d3Iyt7NVk2zl78/NjU+eZGPNWbNYcfroo5x/dmm1S7oQwnbW95/1/Wgvdk+YihXjKiB7d9gUwsjWrePWIb16ZX7bihWBTz/ltiOTJ2d8O19fDjEtXvzg765e5XCUUtxOZf584Ntvbf8AUQqYPp0TxLdu5ZDb4cMc3jJqonS/3LmB998HTp3iMOm4cezCvnt39o9p3cj30iX7xCiEyDrr+8+eG2sDGiRMjzzC6xMn7H1kIYwrIoLJSqlStt3+9deBt97ivmijRmU8J6Z+fW5zcuVK2p/fugXkycM5TmPGAJ98wuTKFufOsYLUrx/QsSO3HOnd23U3lC1YkHOxIiLYCqF+fW4/k51qU8GC/DtJhUkI/Vy+zBO7IkXse1y7J0wVKvD68GF7H1kI4zpyBHj0UdtvbzIBX3zBqsf48ZxonV6ftxo1Uo5vpRQQHQ3cvMnhvcmTgREjMn9MpVhJql6dE6RXr+b/BwXZHreRVasG7NjBv/nXX7PalNXeeh4ebA8hFSYh9HPpEnu32fskz+4JU5UqvD5xAvj444/xxBNPwN/fH3nz5rX3QwlhGHfusOKTFSYT8N57nDOUmAg89hirHzNnpmw0az3m7dvcjHbKFKBChVgkJQFJSTFITKyKYsWWZfpYCQkcLuzbF3jpJZ7wtGqV9X+n0Xl5MbmMiOBcp4YNgTlzsnaMkJCcVZjGjx+PunXrIk+ePChcuDDatWuHE1KyN5ypU6fi0Ucfvdewsn79+vj999/1DsstXL5s/+E4QIOEydqN98wZIDExES+++CIGDBhg74cRGggPD0eVKlVQt25dvUNxOb6+wN272btvnTrAgQPAkiWcd9OnD1C1KiePh4byNi+9xGrT4MFA3rzcmbdJk1sAjmZ6/MhINsT86SdOTJ85032qShmpVg3Yto0T63v2BIYOtX1eZtGiOUuYtmzZgkGDBmHnzp1Yt24dkpKS0LJlS8TFxWX/oMLhihcvjk8//RQRERHYu3cvnnrqKbRt2xZHUpeDhSYuX7b/hG8A9m8roJRSJpNSzZql/P+sWbNUUFBQltb3Cf1IWwH7e/NNpcqUsc+xbt1SasMGpb74gs0jAaVGjVJqyxalYmKU2rmTP/vwQ6UAqKVLl2Z4rD17lCpWTKmQEKV27bJPfK7EYlFq0iSlPD3ZVuHmzczv07+/UjVr2i+Gq1evKgBqy5Yt9juo0EW+fPnUzCx2UxRZV7cumwBnkePbCgA8m75/EqoQ7qxOHeD0aW7bkVN58wJPPcWqR4UKnGj84YdA48Ycojv6X1GpUqWHH2f5cqBRI56J7d3LIT+RlsnEiferV3M+U716mTe7zGmF6X7R0dEAgPz589vvoMKhzGYzFi5ciLi4ONSvX1/vcFyeVhUmTRKm3LmBGze0OLIQxtS8ObtAz5tnv2OazWwX8OyzaZf6//03r6tWzfi+//sf96t77jlgyxaNytcupHlzthtISmJieupUxrcNCWFbB3u0VrFYLBgyZAgaNGiAatWq5fyAwqEOHTqE3Llzw9fXF/3798fSpUtRxTrRV2jCYuE0A6eYwzR8+HCYTKaHXgICEhCT6b6/QriPwoW5aW54OCdw28PKlZwreP8+lmfP8rpixfTv99NPQOfObEi5cCG3WhGZK1cO+OMP9rdq3JiNL9MTEsIVh5GROX/MQYMG4fDhw1i4cGHODyYcrmLFijhw4AB27dqFAQMGoEePHjh6NPN5hSL7rl/nyYpTJExDhw7FsWPHHnopV84LcXHSvFKI1IYN43DOxx/n/FjR0ezV1KIF93BL7cgRVrPSW1K7eDEnMnfrBsye7bq9lbRSvDgrcnnyAM2apT88V64cr48dy9ljDR48GL/++is2bdqE4sWL5+xgQhc+Pj4oV64cQkNDMX78eNSoUQPffPON3mG5NGs+an0f2lOWPy4LFSqEQoUKPfQ2zZoBGzYAmzbxA10IwT5M773HhOnJJ7kyLTssFjaVvHWLfZLu77x95kz6Q2zr1rGq9OKLwPffc+sUkXXBwfx8a9yYc8l27uQWNVblyzOhiojI3uefUgqvvfYali5dis2bN6N06dL2C17oymKxICEhQe8wXJq1JUhmczizQ5M5TB068HrBghgcOHAA58+fh9lsxoEDB3DgwAHExsZq8bBCOL2RI5kotWkDrF+f9fsnJrJf0s8/cw+0+zuHx8cDsbEKjzwSjQMHDgAAzpw5g+XLj6JDBwuaNwfmzpVkKaeKFWPSdPs2P+9SD7N6eLDpZURE9o49aNAgzJs3DwsWLECePHkQGRmJyMhI3M1uXwqhixEjRuCPP/7A2bNncejQIYwYMQKbN29G165d9Q7NpUVEADVralQ9t2UpncpiWwGluAy3QIFTCsADl02bNmV5vZ9wHGkroK3YWKVatWL7jTffVCouzrb7RUQoVb26Ul5eSi1YkP5tlixhSwHgtVTvuUAFHFWBgZfUrVt2+2cIpdTWrUp5e7OVQGpvvaVU6dLZO2Z6n5kA1KxZs3Icr3CcXr16qVKlSikfHx9VqFAh1axZM7V27Vq9w3J5FSoo9dpr2bprpnmQSSllc26VlUSsWDF2N751Kyv3Es4gJiYGQUFBiI6ORmBgoN7huCSzmdtvvP8+h2969+a+bVWrcv6R1c2b3Px2+nTg9985rDd7Ns+g0tOnDxtPXrjA+TZmMxAWxiaMu3enbF0k7Of77/n8TZkCWHv0LlgAdO3K1cLSDUAIx4iOZtuV2bOBHj2yfPdMtxXXZEgO4Ad/VJRM/BYiPZ6ewNtvc4J2167A1KlArVrs3l21Kv+7VCmgQAGgXTuuuJo5k0lPRskSwN97ezNZAjhnavVqYNEiSZa08uqrnID/+uucEA6kdGDft0+/uIRwN/v389r6/rM3zdbINGnCSaYy8ds4wsPDER4eDrPZrHcobqNsWWDiROCTTzj2HhHBBpfJyaw81ajBppflyz84uTs9qSd8b9kCfPYZL+64L5wjffkl9/Lr1g04dCjtxO/mzfWOTgj3EBHBNilaTPgGoN2Q3IkTDPq114BJk7IemNCPDMkZU3w8Pyxat2Z/pUcf5dD45s2ciCy0de4c96Dr3JlDqE2asP/Wzz/rHZkQ7qFLF/ah2749W3fXb0iuYkUOO+zYodUjCCFSs26E3qwZMHw4twf44QdJlhylVCngiy/Y6mHtWg4LZHelnBAi6yIitBuOAzRMmACgSJGMu+EKIezrt994Xbw4MHkyMH68Ns3bRMb69uUQXO/eQJUqHCK9eVPvqIRwfdHR3BaqTh3tHkPThMk68dteW0EIITK2axcnfA8fzk11X3tN74jcj8nEyfm3bgEbN/JnUmUSQnvWBRaGrTC1b8/rWbO0fBQhhMUCHD8O5MvHqsZ338lQnF5KlQLGjOHKxEKFgDVr9I5ICNe3Zg3fb5Ura/cYmn6kvvIKz7gkYRJCW6tXA0lJ7DzdsyeHg4R+Bg3iasU8eYDly9lKVAihnRUrgOee03YXA00TJh8fLpu29kYQQmhj8mReJyezuiH05ecHfPghW0ScPMlVw0IIbZw8yc2uw8K0fRzNi/YvvMA5TBs2aP1IQrivrVt5/frrQIkS+sYiqHt3rhb28ODZrxBCGytXAr6+2vd81DxhevNNXn/7rdaPJIR7OnYMiI3lZpMjRugdjbDy8uJKRYsFmDdP72iEcF0rVrCdSu7c2j6O5glTcDAnYlm3DBBC2Nfnn/P65Ze5lYpwHu3a8TPw8GHg2jW9oxHC9Vj329R6OA5wQMIEsEwWFQWcOuWIRxPZFR4ejipVqqBu3bp6hyKyYNkyXr//vq5hiHSYTJwArhQ35BVC2Nfvv3OT8dattX8szbZGSW3vXqBuXe7kPWVKdo8iHEW2RjGOqCi2EihYUCoYzioqipW/SpW42bIQwn46duTiij17cnwo/bZGSa1OHcDfH/j1V0c8mhDuY9gwXvfrp28cImN587KZnnWumRDCPhIT2VKlTRvHPJ7DWts9/jhw4YJ8YAhhT//7H69Hj9Y3DvFww4ZxWM4630wIkXN//AHExDhm/hLgwISpf39ef/aZox5RCNcWHc3hniJFuCWKcF4vvMC+dDKPSQj7WbQIKFkSqFHDMY/nsITppZc4LDdtmqMeUQjX1rs3r996S984ROZMJuCJJ7jw5epVvaMRwviiongC0rs331+O4NDdpjp25MRU66aUQojsW7GCHxRvv613JMIW77zD67Fj9Y1DCFcwdy7nMFlPHB3BIavkrK5fBwoX5iTw3btzejShFVkl5/yWLuUwT9Wq7PEjnJ/Fwip7rlzAjRuyObIQ2aUUN9mtUYPDcnbiHKvkrAoWBGrXZpuBmzcd+chCuBbr6rj33tM3DmE7Dw+geXMOJaxbp3c0QhjXpk3cn3HgQMc+rsPPccaPZ3Y4fLijH1kI13DpEvDPP/xvRy2nFfbx8su8/vprfeMQwsimTAGqVAEaN3bs4zo8YWrRgk3cZLWIENkzdCivH3lE+72ThH1Zm+ivWwecO6dvLEIY0b//cneDgQMdN9nbSpdR9L59gbg42ZBSiKyyWDh/ydsbaNRI72hEVpUpAwQF8fmbPl3vaIQwnhkzAD8/oFs3xz+2LgnTBx8Anp7ARx/p8egiI7KXnPObMgVISGDiFBqqdzQiq0wmPm/FigEzZ/K5FELYJimJJxrdugF6rEfSJWHy8wOaNgX+/hs4c0aPCER6Bg0ahKNHj2KPHTblEdqYMAHw8uJmk5IwGVNoKCvsV68CS5boHY0QxrF8OXD5Mvel1YNuC1utkx4dPctdCKPasQM4fx6oXp3/X7OmruGIbKpdG7hyBWjYkJ+Dtnd2EcJ9KcX3S8OGwKOP6hODbglTtWrsIbNmjUx+FMIW3btzSKdlSy6ckAnfxlSyJK+7d+cO68uW6RqOEIawahWwfbu+rVR0bZ02bx6zxs6d9YxCCOe3fDlw8iTQvj03sA4J0TsikV3W565UKa4afu89IDlZ35iEcGZmMzBiBPDkk0CrVvrFoWvCVLMm91fasQPYt0/PSIRwbv37c+7SrFkcwy9aVO+IRHZZE6bLl9mX7tgxbvMghEjfggXc0WD8eMe3EkhN9+b8P/3EP4AeSwSFMILp04HISKBPHw7DXb4sFSYj8/MD8uXj8xgayo3JR48G4uP1jkwI55OQwJX1zz8P1Kunbyy6J0wlSwLPPQccPSrbBdzv5s2b6Nq1KwIDA5E3b168+uqriI2Nfeh9nnzySZhMpjSX/v37OyhiYW8WC/Duu/ySnTSJP7t0SSpMRle0KJ9HABg3jslTeLi+MQnhjKZN42KXjz/WOxInSJgA4Mcf2Zfp1Vf1jsS5dO3aFUeOHMG6devw66+/4o8//kDfvn0zvV+fPn1w+fLle5fPP//cAdEKLXz0ERAdza2EvLz4s+hoIG9eXcMSOZQvH59HAChfnjuuf/JJys+EEMDt2zyh6NmTm+3qzSkSprx5OSR34YJ0/7Y6duwYVq9ejZkzZ+Lxxx9Hw4YN8e2332LhwoW4ZD01zYC/vz+Cg4PvXQL16PAlciw5Gfj0U3aGHjUq7c+9vfWLS+Scl1faid4ffADcvcs+W0II+uorICYGGDNG70jIKRImAJg6FfDxAd54Q+9InMOOHTuQN29e1KlT597PmjdvDg8PD+zateuh950/fz4KFiyIatWqYcSIEbhz585Db5+QkICYmJg0F6G/11/nvJbPP+dO91bJySnVJmFM3t7sWmxVtCg/+77+msNzQri7q1eBL74ABg8GSpTQOxpymoTJzw8YMgS4eVPOsgAgMjIShQsXTvMzLy8v5M+fH5GRkRner0uXLpg3bx42bdqEESNG4Mcff8TL1i3SMzB+/HgEBQXdu5RwllenG4uN5Z5JwcHcezE1PVeJCPtQKm0SDADDhgG+vhyCEMLdffIJ3yMjRugdSQqnSZgALhnMnZvDD65a5Bg+fPgDk7Lvvxw/fjzbx+/bty9atWqF6tWro2vXrpg7dy6WLl2KU6dOZXifESNGIDo6+t7lwoUL2X58YR/PP89K0nffPfg7L6+01QlhPElJD1YJ8+YFRo7kJFdpsyLc2cGD3Ddz2DA26XUWTlXY9/AAZs8GOnQAnnkG2LZN74jsb+jQoejZs+dDb1OmTBkEBwfj6tWraX6enJyMmzdvIjg42ObHe/zxxwEAJ0+eRNmyZdO9ja+vL3x9fW0+ptDW4sXA+vVA/fpA27YP/t7bG0hMdHxcwn4SE9MfVn3jDWD+fE5y3buX0xSEcCdJSXz9V6wIDB2qdzRpOVXCBLCTcbNmwIYNbObWvbveEdlXoUKFUKhQoUxvV79+fURFRSEiIgKh/+2yunHjRlgslntJkC0OHDgAAAiRxj2GcOcOF0D4+gKrV6d/m4IFgWvXHBuXsK9r14AGDR78ubc3Txrr1OEKyY8+cnhoQuhq/HhWmHbt4uegM3GqITmrFSuAXLk4d8NVh+YyU7lyZTz99NPo06cPdu/ejW3btmHw4MHo1KkTiv7XhOfff/9FpUqVsHv3bgDAqVOn8NFHHyEiIgJnz57FihUr0L17dzRu3BiP6rVbociStm25WmraNCCjxY0hITIx2MiUenjz0Ro1OC1h/HgZmhPu5a+/eJIwfDibujobp0yY/P3ZmykhgUNz7mr+/PmoVKkSmjVrhmeffRYNGzbE9OnT7/0+KSkJJ06cuLcKzsfHB+vXr0fLli1RqVIlDB06FO3bt8fKlSv1+ieILEg9FNejR8a3CwlJaXoojOf2bSAu7uHd2keMAKpX59CEDL8Kd5CUBLzyClCpUto2Ks7EpJSy9bY239Bemjfn0NycOa43NOfMYmJiEBQUhOjoaOnh5CB37nCozWLhctqH/dmHDgV+/RU4ccJx8Qn7OXGCXwqbNwNNmmR8u7/+4tDc8OEyNCdc39ixvOzapVt1KdP1x05ZYbKSoTnhLmwZirOSCpOxWZ+7zKYVytCccBfOPhRn5dQJkwzNCXdg61CcValS7NN03yJKYRCnTrGXVrFimd9WhuaEqzPCUJyVUydMQMqque3buWpOCFdiy6q4+9WqxeuICO3iEtqJiOCXQ0BA5re1rpo7dkyG5YRrsq6Kmz3b+VbF3c/pEyYgZWiuTx/uWiyEq2jWzPahOKuyZbm/nCRMxhQRkbVhh9RDc3/8oV1cQjja9u3sbO/sQ3FWhkiY/P2BpUtZkq5bV0rTWgkPD0eVKlVQt25dvUNxC6+9BuzcCTz7rG1DcVYmEz9cJGEynqQknk2n2iLSJiNGAI0bs+J+7pw2sQnhSBcuAC+8ADz+ODefNgJDJEwA0KoVM9GrV/nBIexv0KBBOHr0KPbs2aN3KC5v3jxg8mSgdGkgO10fJGEypiNHOCczq2fT3t7A//4H5MkDhIVxDpsQRnXnDtCuHYfgFi82Tkd7wyRMAPDee/yw2LULGDhQ72iEyJ6DB1lRCgjg6qf7N2G1RZ06PEOTBpbGsns3n++aNbN+34IFOT3h9GlOArdY7B2dENpTCujVCzh+HFi+HLhvj3mnZqiECeDQXNmywNSp7M8khJHExKRsifHHH9xwNTuaNuXQ3O+/2y004QC//cbVkLlzZ+/+1aqxOrl4sUwCF8Y0fjywaBEXcWXnxEFPhkuYPDy4KWXu3MxS/9sqTQinZ7FwKCY2Fpg5E6hdO/vHKlQIeOIJVhyEMdy9C6xdyyp5TrRty+kJY8YAS5bYJTQhHGLFCuD994HRozkfz2gMlzABPCvfupX/3bAhEBWlZzRC2KZ9e+DkSaB/f/Ydyak2bfgFfPduzo8ltLdhA5+rNm1yfqyRI4GXXmJLir/+yvnxhNDakSNA166cu2SUSd73M2TCBLCU98MP3JMpNFTG84VzGz8eWLaMqzynTrXPMcPC+AW8caN9jie0tXIlUK4cezDllMnEz7+KFVlxunYt58cUQis3bvDzqnRpDsVlZ96mMzBo2NSjByd/nz6d8zK3EFpZs4YLFgoVAv78037HrVSJX8DLltnvmEIbZjMTprAwJjv2EBDA5/7uXaBDB2m3IpxTUhKroTExHJLL7vw9Z2DohAkAwsM5iXLVKs5pEsKZbNsGPPccl4Xv3Wvf5bMmEz+I/vc/WWbu7Nas4YrGl16y73FLluQ8pp07gU6d+OUkhLNITuYw3NatwC+/AI88ondEOWP4hAngWXvlysCsWWwGKIQz2LsXePJJJjbbtvHLzd769mWytGCB/Y8t7GfKFE7yf+wx+x+7QQOumlu5EujendUsIfRmsbCIsWQJT+qaNNE7opxziYTJw4O9bcqVYzPAd9/VOyLh7g4f5heZUsDmzVnv7GyrUqWA1q1ZaVVKm8cQOXPmDNsJDBxov+G4+7VuDSxcCPz8M9C7t8zpFPqyWIB+/YD583lp107viOzDJRImAPDy4iz8kiWBCROMOwtfGN+JE0yQkpM5FGPtu6SVgQN5wrB9u7aPI7Lnu++491/nzto+Tvv2nFA7Zw5fE5JACz0oBbzxBvD99xz16dhR74jsx2USJoDzQ44dA0JC2NRt5Ei9IzIW2Usu5w4f5mapiYnsYtusmfaP2aIFq6vh4do/lsiau3f5xfHKK9wTU2tduvDxpk1j+wqpNAlHsliAwYM50vPddxwidiUmZftpiGHOV+7c4XLbixeBIUOAr7/WOyJjiYmJQVBQEKKjoxEYGKh3OIaxdy+rScnJTJZat3bcY0+ezLO6w4c5n084hy+/BIYN4zYQ5co57nFnzwZefRV4+WW2H/D0dNxjC/dkNnMY7ocfgBkz+PozmEwHzF0yYQKA+HigShXOH+jXj9musI0kTFm3bRsneCvFYThHVJZSS0jgSULt2tL92VlERwNlygAvvqjP589PP7GxZYcOwI8/cqWmEFpITub+hj/9xCHhl1/WO6JsyTRhcqkhudT8/HhWV7Eiy9Ndu+odkXBV69alrADZvNnxyRLAXb/HjuVeizt3Ov7xxYMmTOCQnF7zKTt35uqkJUvYziA+Xp84hGtLSOBrbdEiLjwwaLJkE5etMFklJ7MT+MGDvP7zTyZTImNSYbLd558Dw4fz7P3PP9nJWy9mMzvgFygAbNqk3YoskbnISG4S/vrr7PKup1WrWGV69FEm1EWL6huPcB2RkcALLwD79jFhattW74hyxH0rTFZeXsD+/VxBEhEBFC/OYTohcsJi4Wtq2DAgf35WM/WeK+/pyS/nLVu4jF3oZ8wYVv2GDdM7EjZO3boV+Pdfrt7cvVvviIQr2LuXr6ezZ4E//jB8smQTl0+YAPZp+uUXVgNu3uQwnXyhiOyKieHE6iVLWLW8eJF7JDmD557jkOCAAYxTON6WLZwGMGYMNwp3BnXq8AuuVCmgcWPOaRIiu376CWjUCChWjK8rLRqyOiO3SJis3nkHWL+eCdRzz3HOhxBZcfgwq5R//80GgXv3OtcQr8kEzJwJ3LoFvP223tG4n7g4djdu2JDLq51JcDDn2HXuzOXe774rXcFF1lgsbNfTpQsXM2zZ4l5DvG6VMAHAU08BJ08CRYoAo0dz6bf0KhG2+OknzhGKi+Oy2Rkz9I4ofY88wgnHM2YAa9fqHY17GTGCe8bNmuWcO7L7+nLZ99dfs+VBmzZczSdEZmJi2LH7s8/4+TJnjnOdLDqCy0/6zkhyMlc2bd/O4ZR9+5ynfK43mfT9oDffBCZOZPPBzZv1n6+UGaXY0PLvv1kVk6dRe1u2sLXExInsieXs1q5lF+YiRbiLfIUKekcknNXJk5yjdPEiV8I984zeEWlCJn1nxMuLvXMGD+Yk8OLFgT179I5KOJvERA6vTJzI+R8XLjh/sgRwaO777zk017evbJOhtatXOczVqJFxNgBv2ZITwE0mzkH5/Xe9IxLOaN06vj6SkoBdu1w2WbKJ2yZMVt9+y9JifDxfFL16ue8QnWyNktavv3IF3LZtQKtWwOnT/H+jKFWKwy+LFgGffqp3NK4rMZHL9uPjudGoMw7FZaR8efbtatgQePZZbqdy+7beUQlnEBvLgkLLlvxu3LULqFRJ76j05bZDcvc7cwZo3pxfioUKcRWdVjvMOzt3H5JLTORY/e+/s7/S5Mms0hjVBx8A48Zxu5Y2bfSOxrUoxSRj1iz2vtJ6o2WtWCxc2ffOO+zj9f33/DwU7mnzZhYPrlzhydagQcY6EcgmGZKzVenSwKlTwKhRwI0bHHZx52qTu7JWlX7/nQlzZKSxkyWAy9vbtuXKliNH9I7GtUydCkyfzq1PjJosAfwyHDAAOHSIDTdbtJBqkzuyVpWaNgVKlGDD59dec4tkySZSYUqHu1eb3LHC5GpVpfvFxgJPPMEVftu2cYm5yJk1a9ieZPBgznFzFVJtck9uWlVKTSpM2SHVJvfiilWl++XOzZVQ8fFsbHntmt4RGdvmzUywn3kG+OILvaOxL6k2uRepKtlO/iQPMXYsl1OWKcM5CsHBnPgmXEN8PCe6tmnDCtO0aVwpaaSJ3VnxyCPAhg3A9ev8Erx+Xe+IjGnrVvZva9QI+Plnrrh1RaVLs9HvlCnAvHlAtWpcMSVcx6ZN3GNw1ixg0iT+f9myekflvCRhysT91aZ69YD69WU/OiOzWIChQ9mbyJWrSumpVIlJ06VLPKO8ckXviIxlwwbg6ae5amjZMtdv3Hd/tallS55gHDqkd2QiJ44cYYX0qaekqpQV8uex0dixTJIaNOAy3LJlWZ24eVPvyERWTJgA5MkDfPUVEBTEJmyuXFVKT7VqbLJ44wb3FTt5Uu+IjGHxYs5ZatwYWLWKTUzdhbXa9NNPwLFjQI0a7Dt19qzekYmsOH8eeOUVVpUOHmTlUKpKtpOEKQtKlgT+/BPYvx+oUoXViUKFgB49OLwjnNe8eZzA+u67bNT39decx9Oxo96R6aNyZQ4tAayWrF+vbzzOzGLhSsMOHXhWvmwZkCuXzkHpwMMD6NSJCVN4ODuFV6wIDBkic+Kc3fXrwFtvse/Wb78B33wDHD8OdO0qVaUsUUrZehH3Wb9eqZIllQKU8vFR6u23lTKb9Y4q56KjoxUAFR0drXcoObZ2rVIlSqQ8R+++6xrPkb3cuqXU008r5emp1DffKGWx6B2Rc7l9W6kXXlDKZFLq44/l75NabCz/JoGBSuXOrdSYMUrFxOgdlUjt9m2lPvqIz1GePEqNHcufiXRlmgdJwmQHCxYoVaAA/5oBAUp9+aXeEeWMKyRM+/crVaUKnxNPT6V69lTq7l29o3JOyclM9gGlevWSv5PV6dNKPfook4Hly/WOxnldv67U0KFK+foqVaiQUpMmKRUfr3dU7i0hQanJk5UqUoQnim++qdTVq3pH5fQkYXKkL79kwgTwQ3bAAKWMlHNMnjxZVa5cWVWoUMGwCdOcOUqVL8/nwGRSqnVrpW7c0DsqY5gzh196VaootXu33tHox2JRato0npGXKaPU4cN6R2QM588z4fbwUKpYMaU+/FCpS5f0jsq9XL6s1LhxrKqbTEr16KHU2bN6R2UYkjA5mtms1MiRSuXPn/KlXa+eUlu36h2Z7YxWYbpxQ6lXX01JVj09lWrZUj4osuPQIaVCQ/mlN3y4+1UKzp1TqkULvo5691YqKkrviIzn6FGl+vZVyt9fKS8vpV58UalNm2Q4UysWi1JbtijVsSP/3rly8bUriX6WScKkp99+U6pWLf6VAaUKF+Y4f0KC3pE9nFESpvXr+eVuMvHvW7CgUu+/7/x/X2eXlMS5Kd7erDbt2qV3RNozm5WaPp1VpeLFlVq9Wu+IjC8qisNzlSrx/Vm5slLffitJqL1ERysVHq5U1ar8+1aooNTEiZyXKLJFEiZncOUKS6O5cvEv7uWl1DPPKHX8uN6Rpc+ZE6aEBFbwrHPGTCal6tRRasMGvSNzPdZqk8mkVNeuSp06pXdE9mexKLVmjVK1a0tVSSsWi1IbNyrVoQOrvwEBSvXrp9SBA3pHZkwHD3K6R+7c/Hu+8AJPHqWCl2OSMDkTs1mpmTM5L8JadSpWTKn+/ZU6eVLv6FI4W8J09y7PnEJD+QFhnVzfu7fMT9JaUhLn8xQtyorT4MFKRUbqHZV97N6t1FNP8fX0xBNK/fGH3hG5vosXlRo9WqmQEP7d69VT6vPPnffk0VmcOKHUhAl8nQJKBQcr9cEHSl24oHdkLiXTPEg239XJiRPAm2+ygeCdO/xZUBDQpAn39WnRQr/YnGHz3YsX2Vxy2TI2x1OK/UIqVABGjGDTPOE4d+4A337LTTmTktgVeMAA9iYzEqWAHTuAL78EliwBqlYFxo/nViemTLfeFPaSlAQsXw7MmcMeYPHxfG+3aQOEhXGjaFfdcsYWZjNfpytW8HLiBLvKN2vGz77nn+cm4cKuMv0EkITJCezbx+Rg3Trg6lX+zMcHqFkT6NkTePVV/r+j6JUwbdvGXd83bkzpoO7nx61LevUCunVz7w9RZ3DrFpOmqVOBuDh+wQ0cyN3snbkBXmwssGAB90X76y92Nh41Cnj5ZcDTU+/o3NudO9xyZsUKYOVKbteTPz+7qoeFAa1asTu/q4uNZTPQFSu4IfiNG0DhwnyPtWnD91hAgN5RujRJmIzm6lUmDb/8wi0rlOKZb5ky3Mfu6afZbTh3bu1icETCZLHwDGrZMnac/uuvlG7p+fNzn7MhQ4CGDTV5eJFDsbHA/Pns+HzoEFCuHNCnD898y5fXOzqynqUvWgTMncuYW7dmgteihXMneO7KYuFWRStXMnE4dIgni02bct+z0FCgdm0gXz69I825qCieLEdE8CRx40ZuAl6tWkql7bHH5HXqQJIwGdmHH47H3LkmnDnzJJSqASBlPwZ/fw6HhIamTaKUUhg9ejRmzJiBqKgoNGjQAFOnTkX5LHyL2TthSp0cbdsG/P03K0jWl56HB/8t7dpxmNJowzzuTClg+3ZWbhYvBhISuMFvWBgv9eo5toKT3ll6cDArlH37AqVKOS4WkXNnz6YkTzt38vkFeAIZGsrqsxGSqNTJUUQEsHcvN3UHWDWqVy+lklSmjK6hujNJmIxs9OjRyJs3Ly5evIjvv/8eJ09GYckSDt0dOMB5Pqn3sPP3BwICruPWrU1o3746atbMhdWrv8TZs6tx/PhB+Nm4tXp2EiaLBbh0iRs6Hj/O6tiBA+knR4UKcS+zRo2AF17g0KMwvrg4zkexDq1cu8b9++rV45ea9VK0qH3mC5nNnNuxd2/KF9GePSln6WFh/AKSs3TXYLHw8yR10rF//4NJVGgoK54hIXythYQAvr7ax5eQAFy+zMulS8Dp0+knR7VqpSR6oaGcuyXDwk5BEiZXMHv2bAwZMgRRUVEP/O76daRKohROnkwA4Iu0z72CyaSQK5cH8uThmVhwMFCsGM+4c+dm2dt6sVjuYODAfpg0aTo8PXMhMZFfQvHx/DA4dw6IjORjR0dzDkJy8oNxS3LkvsxmYPdublC9Zw+/OKwbtBYpwtdESEjaL7UiRThnzdubr53kZF5u3075IrJ+Gf37L3DkCJM0gF+QdeoA9etz2E3O0t1D6iTKmjinTqKs8udP+1pL/d+BgZwbmfri6cnXsPU1mPq1eOlS2tei9do679LKmhylroRJcuTUJGFyBQ9LmFI7ffo0ypYtix079sPLqyaOHuUZ+IwZa+HtXQaBgeVw8ybf9AkJ/LDJLh8ffiAEBTEpKlqUyVe5ckCVKkD16pywKATACuPFi2nPuFN/6dz/BZeeggXTftFVrcovoVq1gLx5Nf8nCINQislL6oTm/uTGep2QkPXj+/qmfR3ef2397/z5ZeWlwWT6bMmaIxcSGRkJAChVqghCQnhWAwD//DMTJpMJixYtSnP7xESenX399Uz88MOPYGXKF4AfAE/Mm/cD8uXzh68vz/x9fYFHHuEXlxBZYTIBJUrw0q7dg7+/fZuroxITeSZvNqec7QcEsCLqyJWiwrhMJg4FFyjAodmMKMW5RXFxbHNwfzXp/qqTtzdfi3nzSiLkriRhcrDhw4fjs88+e+htjh07hkqVKmkei48PP1CmTOmGr79+6d7PY2JiUKJECbRpMx06tWESbiZPHvdYOi6ch8nE6QnOPFlcOBdJmBxs6NCh6Nmz50NvUyabEzCCg4MBAFeuXEFISMi9n1+5cgU1HzJ5yNfXF76OmBUphBBCGJQkTA5WqFAhFCpUSJNjly5dGsHBwdiwYcO9BCkmJga7du3CgAEDNHlMIYQQwh1kZdK3cDCTyVQSQH4AYQDeAdDov1+dVErF/neb4wBGKKWW/vf/wwAMB9ADwBkAHwF4FEAVpVQ8bGAymQIBRAMIUkrF2O9fJIQQQhiTVJic21gw8bHa/991UwCb//vvigCCUt3mcwABAKYDyAvgTwBP25os/ef2f8e8neWIhRBCCBckFSYhhBBCiExI/1shhBBCiExIwiSEEEIIkQlJmIQQQgghMiEJkxBCCCFEJiRhEkIIIYTIhCRMQgghhBCZkIRJCCGEECITkjAJIYQQQmRCEiYhhBBCiExIwiSEEEIIkQlJmIQQQgghMvF/s7qsY0r/UvUAAAAASUVORK5CYII=\n", "text/plain": [ - "Graphics object consisting of 8 graphics primitives" + "(-1, 4, 4, 6)" ] }, "execution_count": 15, @@ -398,7 +373,79 @@ } ], "source": [ - "draw_circles_from_coords(W) + draw_circles_from_coords(S4 * W)" + "S2 * root" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(-1, 2, 2, 6)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "S3 * root" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For some reason, I have no idea why, $s_4$ doesn't seem to work. Instead this weird matrix does appear to work." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(-1, 2, 4, 8)" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "S4 * root" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(-1, 2, 4, 14)" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "matrix([\n", + " [1, 0, 0, 0],\n", + " [0, 1, 0, 0],\n", + " [0, 0, 1, 0],\n", + " [4, 4, 4, -1]\n", + "]) * root" ] }, { @@ -412,7 +459,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 19, "metadata": {}, "outputs": [], "source": [ @@ -436,7 +483,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 20, "metadata": {}, "outputs": [ { @@ -494,14 +541,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Futile attempt to derive the quadratic form for the cubic packing\n", + "# The quadratic form for the cubical packing\n", "\n", - "We attempted to derive the quadratic form for the cubic packing based on the data collected above, but the results don't match the results from the Stange at all and are significantly uglier. There's probably a mistake in here somewhere, but I'm not sure where." + "We first made a matrix, $W_c$, whose rows are the abbc coordinates of the circles in the root octuple. Then we found the row echelon form of that matrix, resulting in a system of linear relations the coordinates must satisfy. Then, from there, we could derive the rest of the coordinates from the first four (we chose the first four to be the \"cubicle\" from the Stange), allowing us to derive the quadratic form." ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 21, "metadata": {}, "outputs": [], "source": [ @@ -509,76 +556,98 @@ " return [(h1^2 + h2^2 - 1) / b, b, h1, h2]" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The idea here is that by multiplying $W_c$ by an arbitrary vector in $\\mathbf{R}^8$ and setting that equal to $\\vec{0}$, we can find a basis for the null space of $W_c$, which will end up giving us a bunch of linear relations the curvatures must satisfy. The thing to notice is that $W_c\\vec{v} = 0$ for some $\\vec{v}\\in\\mathbf{R}^8$ gives the system of equations with coefficient matrix $W^T$. So finding the row echelon form of $W^T$ will give us the coefficients in the simplified system of linear relations, with each row equal to 0." + ] + }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[ 1 -3 -3 -5/4 -1 -1 -5 3/4]\n", - "[ -3 1 -3 1/2 -1 -1 -1 -3/2]\n", - "[ -3 -3 1 -5/4 -5 -1 -1 3/4]\n", - "[-5/4 1/2 -5/4 1 3/4 -3/2 3/4 -1]\n", - "[ -1 -1 -5 3/4 1 -3 -3 -5/4]\n", - "[ -1 -1 -1 -3/2 -3 1 -3 1/2]\n", - "[ -5 -1 -1 3/4 -3 -3 1 -5/4]\n", - "[ 3/4 -3/2 3/4 -1 -5/4 1/2 -5/4 1]" + "[ 1 0 0 0 1/2 1/2 1/2 -1/2]\n", + "[ 0 1 0 0 1/2 1/2 -1/2 1/2]\n", + "[ 0 0 1 0 1/2 -1/2 1/2 1/2]\n", + "[ 0 0 0 1 -1/2 1/2 1/2 1/2]" ] }, - "execution_count": 19, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Wc = matrix([\n", - " abbc_coords(sqrt(2), -1, sqrt(2)),\n", - " abbc_coords(2*sqrt(2), 1, 0),\n", - " abbc_coords(sqrt(2), -1, -sqrt(2)),\n", - " [sqrt(2) / 4, 0, 1, 0],\n", - " abbc_coords(sqrt(2), 1, sqrt(2)),\n", - " abbc_coords(2*sqrt(2), -1, 0),\n", - " abbc_coords(sqrt(2), 1, -sqrt(2)),\n", - " [sqrt(2) / 4, 0, -1, 0]\n", - "\n", + " [4, 0, 0, 1],\n", + " [0, 2, 0, 1],\n", + " [2, 1, -sqrt(2), -1],\n", + " [2, 1, sqrt(2), -1],\n", + " [2, 1, -sqrt(2), 1],\n", + " [2, 1, sqrt(2), 1],\n", + " [4, 0, 0, -1],\n", + " [0, 2, 0, -1],\n", "])\n", "\n", - "Wc * P * Wc.transpose()" + "Wc.transpose().echelon_form()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "From here we can conclude that $$\n", + "\\begin{align*}\n", + " 2b_1 &= -b_5 - b_6 - b_7 + b_8\\\\\n", + " 2b_2 &= -b_5 - b_6 + b_7 - b_8\\\\\n", + " 2b_3 &= -b_5 + b_6 - b_7 - b_8\\\\\n", + " 2b_4 &= b_5 - b_6 - b_7 - b_8\n", + "\\end{align*}\n", + "$$\n", + "This differs slightly from Stange's system of equations because we put the circles in a slightly different order." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In any case, this gives us the tools to derive the full octuple from just four coordinates, letting us use those four coordinates to represent the entire octuple, and finally making $WPW^T$ nonsingular, allowing us to find the quadratic form the same way we did for the Descartes quadratic form and the octahedral quadratic form." ] }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "W = matrix([\n", - " abbc_coords(sqrt(2), -1, sqrt(2)),\n", - " abbc_coords(2*sqrt(2), 1, 0),\n", - " abbc_coords(sqrt(2), -1, -sqrt(2)),\n", - " [sqrt(2) / 4, 0, 1, 0]\n", + " [4, 0, 0, 1],\n", + " [0, 2, 0, 1],\n", + " [2, 1, -sqrt(2), -1],\n", + " [2, 1, sqrt(2), -1],\n", "])" ] }, { "cell_type": "code", - "execution_count": 21, - "metadata": { - "scrolled": true - }, + "execution_count": 24, + "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[ 1 -3 -3 -5/4]\n", - "[ -3 1 -3 1/2]\n", - "[ -3 -3 1 -5/4]\n", - "[-5/4 1/2 -5/4 1]" + "[ 1 -3 -3 -3]\n", + "[-3 1 -3 -3]\n", + "[-3 -3 1 -3]\n", + "[-3 -3 -3 1]" ] }, - "execution_count": 21, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" } @@ -590,23 +659,346 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "97*b1^2 - 304*b1*b2 + 328*b2^2 - 290*b1*b3 - 304*b2*b3 + 97*b3^2 + 32*b1*b_avg - 1088*b2*b_avg + 32*b3*b_avg + 1280*b_avg^2" + "[ 5 -3 -3 -3]\n", + "[-3 5 -3 -3]\n", + "[-3 -3 5 -3]\n", + "[-3 -3 -3 5]" ] }, - "execution_count": 22, + "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "D = W2.transpose() * m.inverse() * W2\n", - "968*factor(simplify(D[1][1]))" + "32 * m.inverse()" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[bt1 b1 h11 h21]\n", + "[bt2 b2 h12 h22]\n", + "[bt3 b3 h13 h23]\n", + "[bt4 b4 h14 h24]" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "W2 = matrix([\n", + " [\n", + " var('bt' + str(i)),\n", + " var('b' + str(i)),\n", + " var('h1' + str(i)),\n", + " var('h2' + str(i)),\n", + " ] for i in range(1, 5)\n", + "])\n", + "W2" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "D = W2.transpose() * m.inverse() * W2" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "5*b1^2 - 6*b1*b2 + 5*b2^2 - 6*b1*b3 - 6*b2*b3 + 5*b3^2 - 6*b1*b4 - 6*b2*b4 - 6*b3*b4 + 5*b4^2" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "32 * simplify(expand(D[1][1]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This spits out the lovely quadratic form $$\n", + " 8(b_1^2 + b_2^2 + b_3^2 + b_4^2) = 3(b_1 + b_2 + b_3 + b_4)^2\n", + ".$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's the generators for the Weyl group. They might or might not be generators for the packing itself, however." + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[\n", + "[-1 6 6 6] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 6 -1 6 6] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 6 6 -1 6] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 6 6 6 -1]\n", + "]" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "S_i = weyl_generators(m, standard_basis(4))\n", + "S_i" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[ 1 -3 -3 -3 -1 -1 -1 -5]\n", + "[-3 1 -3 -3 -1 -1 -5 -1]\n", + "[-3 -3 1 -3 -1 -5 -1 -1]\n", + "[-3 -3 -3 1 -5 -1 -1 -1]\n", + "[-1 -1 -1 -5 1 -3 -3 -3]\n", + "[-1 -1 -5 -1 -3 1 -3 -3]\n", + "[-1 -5 -1 -1 -3 -3 1 -3]\n", + "[-5 -1 -1 -1 -3 -3 -3 1]" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Wc * P * Wc.transpose()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Generalization of this method\n", + "This method seems remarkably general. Given a matrix representing a root unit of a packing, we can find the linear relation between the coordinates, and thus represent the packing with only four coordinates, allowing us to find the quadratic form." + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "def quadform_from_root(root_matrix):\n", + " P = matrix([\n", + " [0, -1/2, 0, 0],\n", + " [-1/2, 0, 0, 0],\n", + " [0, 0, 1, 0],\n", + " [0, 0, 0, 1],\n", + " ])\n", + " \n", + " # step 1: find linear relation between coords\n", + " relation = root_matrix.transpose().rref() * vector([ var('b' + str(i)) for i in range(1, root_matrix.dimensions()[0] + 1)])\n", + " \n", + " # step 2: find matrix of quadratic form\n", + " W = root_matrix[0:4]\n", + " M = W * P * W.transpose()\n", + " \n", + " # step 3: repeat with arbitrary matrix\n", + " W2 = matrix([\n", + " [\n", + " var('bt' + str(i)),\n", + " var('b' + str(i)),\n", + " var('h1' + str(i)),\n", + " var('h2' + str(i)),\n", + " ] for i in range(1, 5)\n", + " ])\n", + " D = factor(simplify(expand(W2.transpose() * M.inverse() * W2)))\n", + " \n", + " return relation, M.inverse(), D[1][1]" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(2*b1 + b5 + b6 + b7 - b8, 2*b2 + b5 + b6 - b7 + b8, 2*b3 + b5 - b6 + b7 + b8, 2*b4 - b5 + b6 + b7 + b8)\n", + "[ 5 -3 -3 -3]\n", + "[-3 5 -3 -3]\n", + "[-3 -3 5 -3]\n", + "[-3 -3 -3 5]\n", + "5*b1^2 - 6*b1*b2 + 5*b2^2 - 6*b1*b3 - 6*b2*b3 + 5*b3^2 - 6*b1*b4 - 6*b2*b4 - 6*b3*b4 + 5*b4^2\n" + ] + }, + { + "data": { + "text/plain": [ + "[\n", + "[-9 6 6 6] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 6 -9 6 6] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 6 6 -9 6] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 6 6 6 -9]\n", + "]" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# cubical\n", + "relation, mat, equation = quadform_from_root(Wc)\n", + "print(2 * relation)\n", + "print(32 * mat)\n", + "print(32 * equation)\n", + "weyl_generators(32 * mat, standard_basis(4))" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(b1 + b5 + b6, b2 - b5, b3 - b6, b4 + b5 + b6)\n", + "[ 1 -2 -2 -1]\n", + "[-2 4 0 -2]\n", + "[-2 0 4 -2]\n", + "[-1 -2 -2 1]\n", + "b1^2 - 4*b1*b2 + 4*b2^2 - 4*b1*b3 + 4*b3^2 - 2*b1*b4 - 4*b2*b4 - 4*b3*b4 + b4^2\n" + ] + }, + { + "data": { + "text/plain": [ + "[\n", + "[-1 4 4 2] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 4 -7 0 4] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 4 0 -7 4] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 2 4 4 -1]\n", + "]" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# octahedral\n", + "relation, mat, equation = quadform_from_root(matrix([\n", + " [2, 0, 0, 1],\n", + " [2, 0, 0, -1],\n", + " [-1, 1, 0, 0],\n", + " [4, 2, 2*sqrt(2), -1],\n", + " [4, 2, 2*sqrt(2), 1],\n", + " [7, 1, 2*sqrt(2), 0],\n", + "]))\n", + "print(relation)\n", + "print(8 * mat)\n", + "print(8 * equation)\n", + "weyl_generators(8 * mat, standard_basis(4))" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(b1, b2, b3, b4)\n", + "[ 1 -1 -1 -1]\n", + "[-1 1 -1 -1]\n", + "[-1 -1 1 -1]\n", + "[-1 -1 -1 1]\n", + "b1^2 - 2*b1*b2 + b2^2 - 2*b1*b3 - 2*b2*b3 + b3^2 - 2*b1*b4 - 2*b2*b4 - 2*b3*b4 + b4^2\n" + ] + }, + { + "data": { + "text/plain": [ + "[\n", + "[-1 2 2 2] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n", + "[ 0 1 0 0] [ 2 -1 2 2] [ 0 1 0 0] [ 0 1 0 0]\n", + "[ 0 0 1 0] [ 0 0 1 0] [ 2 2 -1 2] [ 0 0 1 0]\n", + "[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 2 2 2 -1]\n", + "]" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# tetrahedral\n", + "relation, mat, equation = quadform_from_root(matrix([\n", + " [2, 0, 0, 1],\n", + " [2, 0, 0, -1],\n", + " [-1, 1, 0, 0],\n", + " [3, 1, 2, 0]\n", + "]))\n", + "print(relation)\n", + "print(4 * mat)\n", + "print(4 * equation)\n", + "weyl_generators(4 * mat, standard_basis(4))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I'm not sure if this works or not. I suspect that it is highly dependent on the order of the circles in the root. Interestingly, it looks like the only relation it was able to deduce for the tetrahedral packing is $b_1, b_2, b_3, b_4 = 0$, meaning there is no null space, as we'd expect. The octahedral quadratic form we get is very different, which isn't surprising, I guess, since this is a very different coordinate system, but I'm not sure if it's right. Likewise, for some reason, we get the cubical quadratic form from the root tetrahedral packing, which doesn't make sense. For whatever reason, though, we get the right matrix." ] } ],