{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Calculations to find the quadratic form of the octahedral packing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$P$ is the matrix of the quadratic form corresponding to $h_1^2 + h_2^2 - \\tilde{b}b = 1$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"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",
"])"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[ 0 -2 0 0]\n",
"[-2 0 0 0]\n",
"[ 0 0 1 0]\n",
"[ 0 0 0 1]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P.inverse()"
]
},
{
"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. 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."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"W = matrix([\n",
" [2, 0, 0, 1],\n",
" [2, 0, 0, -1],\n",
" [-1, 1, 0, 0],\n",
" [6, 2, 2*sqrt(2), 0]\n",
"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we have $$\n",
" WPW^T = \\left(\\begin{matrix}\n",
" 1 & -1 & -1 & -1\\\\\n",
" -1 & 1 & -1 & -1\\\\\n",
" -1 & -1 & 1 & -1\\\\\n",
" -1 & -1 & -1 & -1\n",
" \\end{matrix}\\right) = M\n",
"$$\n",
"So, inverting both sides, we get $$\n",
" (WPW^T)^{-1} = M^{-1}\n",
"$$ $$\n",
" (W^T)^{-1}P^{-1}W^{-1} = M^{-1}\n",
"$$ $$\n",
" P^{-1} = W^TM^{-1}W\n",
".$$\n",
"\n",
"Like in the case with Apollonian packings, this is true for any sextuple $W$. So, we can substitute an arbitrary sextuple for $W$ and it must be equal to $P^{-1}$, letting us derive some useful quadratic forms."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"[ 1 -1 -1 -2]\n",
"[-1 1 -1 -2]\n",
"[-1 -1 1 -2]\n",
"[-2 -2 -2 -4]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = W * P * W.transpose()\n",
"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,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"[ 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,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M.inverse()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"bt1 = var('bt1')\n",
"b1 = var('b1')\n",
"h11 = var('h11')\n",
"h12 = var('h12')\n",
"bt2 = var('bt2')\n",
"b2 = var('b2')\n",
"h21 = var('h21')\n",
"h22 = var('h22')\n",
"bt3 = var('bt3')\n",
"b3 = var('b3')\n",
"h31 = var('h31')\n",
"h32 = var('h32')\n",
"b5_avg = var('b5_avg')\n",
"b_avg = var('b_avg')\n",
"h1_avg = var('h1_avg')\n",
"h2_avg = var('h2_avg')\n",
"\n",
"\n",
"W2 = matrix([\n",
" [bt1, b1, h11, h12],\n",
" [bt2, b2, h21, h22],\n",
" [bt3, b3, h31, h32],\n",
" [b5_avg, b_avg, h1_avg, h2_avg],\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"D = W2.transpose() * M.inverse() * W2"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[ 0 -2 0 0]\n",
"[-2 0 0 0]\n",
"[ 0 0 1 0]\n",
"[ 0 0 0 1]"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P.inverse()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"b1^2 + b2^2 + b3^2 - b1*b_avg - b2*b_avg - b3*b_avg + 1/4*b_avg^2"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"2 * factor(simplify(D[1][1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So, we end up deriving the quadratic form $$\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": "markdown",
"metadata": {},
"source": [
"# Calculations to find the generators for the \"Apollonian\" Group for the octahedral packing\n",
"\n",
"We can try the Weyl group, since that is one way to derive the generators for the Apollonian group and see if it works. My worry is that it won't since the coordinate system is so different."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def weyl_generators(matrix, alphas):\n",
" retval = []\n",
" for alpha in alphas:\n",
" scale_factor = (alpha.transpose() * matrix * alpha)[0][0]\n",
" retval.append(identity_matrix(len(alphas)) - 2 * alpha * alpha.transpose() * matrix / scale_factor)\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": [
"[ 2 0 0 -1]\n",
"[ 0 2 0 -1]\n",
"[ 0 0 2 -1]\n",
"[ -1 -1 -1 1/2]"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"4 * M.inverse()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"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], [ 4 4 4 -1]\n",
"]"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S_i = weyl_generators(4 * M.inverse(), standard_basis(4))\n",
"S_i"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"S1 = S_i[0]\n",
"S2 = S_i[1]\n",
"S3 = S_i[2]\n",
"S4 = S_i[3]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(7, 2, 4, 6)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"root = vector([-1, 2, 4, 6])\n",
"S1 * root"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-1, 4, 4, 6)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S2 * root"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-1, 2, 2, 6)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S3 * 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": [
"S4 * root"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Calculations to find the curvatures of the dual circles of the octahedral packing\n",
"\n",
"Here the basic idea was to work out the intersection points of the circles and let the computer algebraically find the radii of the circles. The end up coming out, rather anticlimactically, to $(0, 0, \\sqrt 2, \\sqrt 2, \\sqrt 2, \\sqrt 2, 2\\sqrt 2, 2\\sqrt 2)$. It's rather satisfying and a bit surprising that the dual of the root sextuple of the octahedral packing is the root of the cubic packing, since the dual polyhedron of the octahedron is the cube."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"def circle_from_points(pta, ptb, ptc):\n",
" a = var('a')\n",
" b = var('b')\n",
" r = var('r')\n",
" x = var('x')\n",
" y = var('y')\n",
" \n",
" circle_func = (x - a)^2 + (y - b)^2 == r^2\n",
" \n",
" eq1 = circle_func.subs(x == pta[0]).subs(y == pta[1])\n",
" eq2 = circle_func.subs(x == ptb[0]).subs(y == ptb[1])\n",
" eq3 = circle_func.subs(x == ptc[0]).subs(y == ptc[1])\n",
" \n",
" res = solve([eq1, eq2, eq3], a, b, r)[1]\n",
" \n",
" return (res[0].rhs(), res[1].rhs(), res[2].rhs())"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(1/2*sqrt(2), 1)\n",
"sqrt(2)\n",
"(1/4*sqrt(2), 0)\n",
"2*sqrt(2)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAGECAYAAAAySIfuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAABdvElEQVR4nO2dd3gUVRfG300lARIIPRAgEFrohChFEEQRRKpipwiCCCgqKqAi8inSQRGC2MAKKL2qoBBERELoJfSeACEJmxBSd8/3x2GTLCTZNruzuzm/55lnk93ZO2dmZ+57yznnaogIgiAIgmDAQ20DBEEQBOdChEEQBEEwQoRBEARBMEKEQRAEQTBChEEQBEEwQoRBEARBMEKEQRAEQTBChEEQCqBhAjQajUZtWwRBLbys/J5ExQluiVarRWBgILRardqmCIK9MNnokR6DIAiCYIQIgyAIgmCECIMgCIJghAiD4DYsXLgQzZo1Q0BAAAICAtC2bVts3rxZbbMEweUQYRDchho1amDatGmIjY3F3r178dBDD6F37944evSo2qYJgkuhsTLttnglORO3bwNxcUDDhoC/v9rWOBVBQUGYOXMmhg4dWujnWVlZyMrKyvs/NTUVISEh0Gq1CAgIcJSZ7oncl86KeCWVCOLigIgIfhUAADqdDsuWLUN6ejratm1b5H5Tp05FYGBg3hYSEuJAK90cuS9dFhEGwa04fPgwypQpA19fX4wYMQKrV69GeHh4kftPmDABWq02b7t06ZIDrRUE58TaADdBcEoaNGiAAwcOQKvVYsWKFRg0aBCio6OLFAdfX1/4+vo62EpBcG5EGAS3wsfHB2FhYQCAiIgIxMTE4LPPPsOiRYtUtkwQXAcZShLcGr1ebzS5LAiCaaTHILgNEyZMQPfu3VGzZk2kpaXh559/xvbt2/H777+rbZoguBQiDILbcP36dQwcOBAJCQkIDAxEs2bN8Pvvv+ORRx5R2zRBcClEGAS34ZtvvlHbBEFwC2SOQRAEQTBChEEQBEEwQoRBEARBMEKEQRAEQTBChEEQACxYsADh4eGIjIxU2xRBUB0RBkEAMGrUKBw7dgwxMTFqmyIIqiPCIAiCIBghwiAIgiAYIcIgCIIgGCHCIAiCIBghwiAIgiAYIcIgCIIgGCHCIAiCIBghwiAIgiAYIcIgCIIgGCHCIAiCIBghwiAIkFxJglAQEQZBgORKEoSCiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIAiQlhiAURIRBECApMQShICIMgiAIghEiDIIgCIIRIgyCIAiCESIMgiAIghEiDIIgCIIRIgyCIAiCESIMgiAIghEiDIIgCIIRIgyCIAiCESIMgiAIghEiDIIAyZUkCAURYRAESK4kQSiICIMgCIJghAiDIAiCYIQIgyAIgmCECIMgCIJghAiDIAiCYIQIgyAIgmCECIPgNkydOhWRkZEoW7YsKleujD59+uDEiRNqmyUILocIg+A2REdHY9SoUdi9eze2bNmCnJwcdO3aFenp6WqbJgguhZfaBgiCUvz2229G/y9ZsgSVK1dGbGwsOnbsqJJVguB6iDAIbotWqwUABAUFFblPVlYWsrKy8v5PTU21u12C4OzIUJLgluj1erz++uto3749mjRpUuR+U6dORWBgYN4WEhLiQCsFwTkRYRDcklGjRuHIkSNYtmxZsftNmDABWq02b7t06ZKDLBQE50WGkgS3Y/To0diwYQN27NiBGjVqFLuvr68vfH19HWSZILgGIgyC20BEePXVV7F69Wps374doaGhapskCC6JCIPgNowaNQo///wz1q5di7Jly+Lq1asAgMDAQPj5+alsnSC4DjLHILgNCxcuhFarRadOnVCtWrW8bfny5WqbJgguhfQYBLeBiNQ2QRDcAukxCIIgCEaIMAiCIAhGiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIAoAFCxYgPDwckZGRapsiCKojwiAI4OC4Y8eOISYmRm1TBEF1RBgEQRAEI0QYBEEQBCMk8llQntxc4No1ID4eSEgAUlKAnByACPDyAvz8gKpVgeBgoFo1oGxZtS12XoiApCS+jvHxwPXrfC11OsDTE/DxASpVyr+WFSoAGo3aVgsujgiDYBtpacCePUBsbP529ixXaOYSGAi0bAlERPB2331A3br2s9lZyc0FDhwA9u7Nv5ZHjwLZ2eaX4esLhIfnX8vWrYEWLViQBcFM5G4RLOfyZWD9emDdOuCvv7jiKlOGK/eePbliMrRgg4OBoCCumDQaQK8H0tO5BWzYzp0D9u0DVq4EZs/mY9SvD/TqxVvbtu5bsaWlAb//ztdz40buHXh6Ak2acMU+aBAQEpJ/LStX5srfw4OvZVYW9yIMvbNLl4BDh1isFy/mnkWFCkCPHnwtu3aVHppgEjd92gTFyc4GVq8GoqKAHTu4on7wQWDmTOCRR4AGDbiyMoWnJxAQwFuDBvd+npQE/PMPV5Q//ADMmsWV4UsvAS+/DNSsqfy5ORoiYPduvpa//sqVe9OmfH6PPQa0asXDbabw8OD9atXi7W4yMlhwN21iEf/+exaV/v2BkSOBNm1k2EkoHCKyZhOcidhYIoBflSYxkej994mqVOFjPPgg0Y8/EqWkKH+su9HpiHbvJho9migggMjDg6hXL6Lt2+12SK1WSwBIq9UqX3h2NtHXXxO1aMHXsm5dohkziM6eVf5YhXH2LNH06UR16vDxW7Rge7Kz7XM8e96Xgi2YrONFGNwBezyAaWlEkycTlS1LVKYM0ahRREeOKFe+NfYsWkTUtCmf66OPEu3bp/hh7CIMOh3RsmVEYWFEGg1Rz55Emzfz+2qg0/Hxe/Zke8LC2D6l7RFhcFZM1vHirioYo9cDCxfy5O+UKcCwYTwHMH8+0LixenaVKQMMHw4cPAisWAGcP89DLs8+y3MezsqOHUBkJPDMMzx0duAAD+t062be0Js98PDg469bB+zfz/M5zzzDdu7YoY5NglMhwiDkc/o00KkTjz937w6cOsWTwRUrqm1ZPhoN8MQTwJEjwFdfAdHRLFjffmuZJ5S9uXULePVVnofx9uYKd8MGoFkztS0zpnlznvSOjmY7H3wQeO01dhAQSiwiDAL3EubN40rr8mX2NFqyxLkner28eEL66FGgb19g6FCeuLWy96BorqToaL6W33wDfPopsGsX0KGD7eXak44d2c5PPwW+/prtl95DiUWEoaRz6xbw5JPAmDHAkCHs6ti5s9pWmU/58ixiGzaw7a1aAX//bXExiuRKIuIe1kMPAdWr87DXmDHqDRlZiocH23vwILvGdu4MzJnjXD0xwSG4yB0r2IXz54F27YAtW4A1a3geoUwZta2yjh49uEJr3Bjo0oWHmRxJZiYweDDw1lvA228D27cD9eo51galqFeP7X/7bWDsWODFF/n8hBKDCENJZdcunmxMT2ef+t691bbIdipWBP74gyfMhw/n1q9eb//j3rjBczO//AL89BMwbRrHa7gynp58Hj/+CCxfzr2HGzfUtkpwECIMJZHt2zkCNjycI2TV9DZSGm9vYMEC9qz6/HOee9Dp7He8q1dZFM6d4yGs556z37HU4Pnnea7h7Fk+z2vX1LZIcAAiDCWNHTt4krZdO2DzZk6X4I6MGMGt9x9+4LkTe/Qcrl/n+YSUFJ5wbt1a+WM4A5GRfH4pKdxzuH5dbYsEOyPCUJLYu5fH4tu3Zx92f3+1LbIvzz7LwvDDD+w6quQkamoqpwJJSQG2bQMaNlSubGekYUM+z+RkPu/UVLUtEuyICENJIT6e5xGaNOGJ5lKl1LbIMTz7LPDll5yXKCpKmTJ1Oh4yunAB2LqVA8RKAvXrA3/+yef9/PP2HaITVEWEoSSQmcm+/hoNJ8IrXVptixzLSy8Br7/Ok9F//WV7ee+/z8Nwy5a51/yMOTRuDCxdykFxEyeqbY1gJ0QYSgIvv8w+/mvW8AI5JZGZM3k+oH9/nki1lqVL2VtnxgxOK1ES6d6dz3/qVBZHwe0QYXB3fvmF0y1/+aX7To6ag5cXu10GBnK8gTWT0Zcu8aT2c88Bb76puIkuxdixPEz38svOnatKsAoRBnfm+nVg1CiObB4wQG1r1Kd8ec6p9PffHMxXAJMpMYg4PqJMGXaHLenrGGg0PGdTpgxfF4mOditEGNyZUaP4dcECde1wJjp1AkaPBsaP56SBdzCZEuPbb3mlta++AsqVc4ipTk+5ctwT/e03TksiuA0iDO7Kxo2cnnrBAl4BTchn6lSeazEIpymSknjoZPBgjgER8unRg5cffeMNvk6CWyDC4I7odMCECdw67t9fbWucjzJleMnQP/4wz0tp6lS+ptOn2982V2TGDCA3V66PGyHC4I4sXQocPswVWkkfCy+Kvn2B++7jIaXixscvXeL5iLFjpedVFJUr8/X5/HOZiHYTRBjcjexs9i/v25cXexcKR6Nht9OYGGDVqqL3mzwZCAjgik8omrFjuSf24YdqWyIogAiDu7FyJafT/ugjtS1xfjp35tiGWbMK//zqVXb1fecdoGxZx9rmagQEAOPG8fWSRHsujwiDuxEVxZVdSYvItZZXX+W04/v23fvZN99w/MPQoY63yxUZOpSv1zffqG2JYCMiDO7EqVPAzp28ZrNgHo8/DoSEcJruguTmAosWcTBb+fLq2OZqlC/PQW9ffCF5lFwcEQZ3YsUKXpKxVy+1LXEdvLw4evennwCtNv/tP/7giWcRWcsYOZKv26ZNalsi2IAIgzvx11/cwvX2VtsS12LgQCAjgxPj3cFr7VrORNuqlYqGuSARETyMuXKl2pYINiDC4E4kJ7vHEp2OJiQEaNmS16gA4AnA6/ff5VpaS69ewIYNMpzkwogwuBOBgUDbtmpb4ZLsqVoVqcuXo23r1mgHwCMlBejZU22zXJNevTgK+tAhtS0RrESEwZ144AHXX4ReJe77+GME6PXYM3cuHgegr1yZl7QULOe++4AqVThZoeCSiDC4A+np/FqS02rbSosWQNmy8IyJQRsAunbtAA95PKzCwwPo0AE4ckRtSwQrkTvfHTh5kl8bNVLXDlfGwwNo1QqeBw6gJQBdy5ZqW+TaREQAcXFqWyFYiQiDO3D8OL+Ghqprh6sTEQHPmBiUBaBr0UJta1ybiIj8nqzgcogwuAOGlpmXl7p2uDoREfCIjwcA6Jo3V9kYFyciQm0LBBsQYXAH7lRmgo2EhQEAUgCJdraVoCDOnyS4JCIM7sCNG2pb4B5UqwYAuKmuFe5DpUpqWyBYiQiDq0MEJCaqbYXTsGPHDvTs2RPBwcHQaDRYs2aN+V+uUgUAICPjClGxotoWCFYiwuDqpKUBmZlqW+E0pKeno3nz5lhgzTrXPj4gLy9kK29WyUSEwWWR2UpXRzw/jOjevTu6d+9ufQEeHtArZ07JpnRptS0QrERRYbh9W1yXlebmTSAhgUeLbtzIf01KArKygIq3c7AewDuYhn8GhcGjHFChAjfWDFulSkDVqvx3yVrpsyXOnAksdKmFnBzg7Fng8OFcHDtGSE7WIDsb2JEDAPXQrZsf/P2BGjU4PCQ8nBPXlqzrVzQ6Hfs8JCYa35uJiUBqKn/++skKuB+18dzAMGgCAH9/43vS8FqjBlCqlNpn5H40bMjX3BoUFYa4OPFSczQ1wDVVEFJw7IhOJk6N2Ie33jK1j/EjwKs/5+Lff+1jUUlCCx28kIu4o7lyX6pAbKz1yYE1VNxC6EVT6Jekx2AZej1w7BiwYwdvp07x+6VLcyvVsNWsya2roKB7UyF5JSagWbdgAMDxH2OR0aiVUflaLbfiLl/m3+b4cd5SUnif6tWBBx8EOnbkrBDulLE7IqIVPvroM9y82QErVgAXLnDrtHlzbk01agSEhmbDyyt/VuGBLkHYl5uD29FalCnD7pZJSfnX7fhx4MABvq733w88+SRfO3cLIUlLA3bt4vvyn3/4f4AT0Rruy/r183uipUvf25uqPe4pBG399Z77EuDeblIScP06cOZM/rU9fZrXSPL05AwvHTvyFhzsoBN3I4rpMZju9xKRNZtgJbm5RJs2EQ0bRlS1KhFAFBRENGAA0U8/EZ08SaTTWVDgrVtcCEAUG2vWV/R6oosXiVatIhoxgig4mL9erhzRs88SLVtGlJFh3fk5C0ePEgFfko9PDnl5ET3zDFF0NJ97cei9vGgfQFqttsh9MjKIfviBqG1bvm7VqxP9739EyckKn4SDuXqV6LPPiLp0IfLy4nNr2ZLogw+I/vyTKCXFwgK7d7foviQiyswk2ruXaN48okcfJfL25iKaNiV67z2igwcttEEoDJN1vAiDg7h2jeiTT4hq1uSrXq8e0VtvcWWVk2Nj4aVLW/wAFkSv569OmsQVAUBUsSLRuHFEZ8/aaJuDyc7mSporlIv03HNHKSHBzC9nZBABdMSEMBRk3z4WeX9/omrViNats952NdDriXbsYOH09iby8eEKecECbjzYROvWNt2XRERaLdGvv3LDKSiIi2vfnhtRmZk22ldyEWFQE72eaOdOouee44euVCmiF18k2rNH4QPVqmXzA1iQEyeI3niDexAaDVGPHkQbN1rYk1GBgweJmjfPJU9PPb30UgIBPjRnzhzav38/XbhwwXQB584RAXTeAmEwcPky0WOP8c8wYIDz9x7S0oiiooiaNMlvqMydq7DdCt+X2dlEK1dyjwYgqlSJaMIEovPnFSm+JCHCoBa7dxN16sRXOCyMaPZsoqQkOx1MgZZZYaSnE339dX4volkzFghTwzGOpmAvITQ0jYBWBJ4Hy9sGDRpkuqBdu4gA0lohDER8XZYsIQoMdN7eQ0YG0Zw5RBUqEHl4EPXtS7Rli51E38aebHEcP0702mtEAQE87DViBFF8vOKHcVdEGBzNsWP8sAHcGluzxgEt7d697fYAEnGF9/ffRB068GE6duQ61BlISyN6+GEiT0+i99+3cXhh8WLSazREAKWePm11MQV7D5MnO4eQ5uYSLV7MQ5menjz8ZU4nymoSEiye+7KGW7eIZs4kKl+eyM+PexAWz4WUPEQYHMXVq0RDh3IrrFYtou++44fRIYwbxz/lv//a9TB6PfcYmjXjw/XuzcNOapGURNSmDVHZskTbtytQ4OjRlFu7NhFA6StW2FSUXk80ZQpfpzfeUFccNm4kCg9nW558kiguzgEH3bDBIcJgICWFRcHPj0Vi5kzuSQqFIsJgb/R6op9/5omxChWIPv1UhUmxJUv4p/zhB4ccTqcj+vFHotq1ed5k9mwHiuAdbt4katWKr/nevQoV2rYtZffvT0kAZbz3niJFLljAP82rrzpeHJKSiF54gY/fubMd5raKY/JkHlNzkDAYiI8neuUV7hW1bCleTEUgwmBPrl7NHzZ66imi69dVMuSff9iId9916GHT04nGjOEJ6nbtHNd7SE/nYa1y5RR88HNyiPz8KGPKFNoCUHaPHgoVTLRoEf88H3ygWJEmWbeO3aEDA7nd4PAeS69eRPff73BhMBATQ9S4Mc87/e9/0nu4CxEGe7F0KbdWK1VidzpViY3ln/KRR1Q5/I4dRHXr5vce7D2n8vTT7B6q6DzHzp1EAN3aupU+BEgfGKhobTJ9Ov9ES5YoVmShJCezVxTA3mSXL9v3eIWSnc2qPXy4asJAxD33d9+V3kMhiDAoTXY2d1VV7yUUxCAMpUsTZWWpYkLB3sPjj7P/uT345Rc+1Z9/Vrjgd94hqlyZtCkpFGEYG9+2TdFDDBzILfhLlxQtNo/Dh4lCQ1XsJRj46y++fj/+qKowGDD0HkqV4vgHQYRBURIT2QXV25voyy/VtqYABmEAiP74Q1VTNm1iF8JGjYhOnVK27GvXOPCuXz/lK72kKlXo18BAql+/PmkA0lWrxrPGCpKczFHm3bsrb/+aNURlyrBjwLlzypZtMa+/zuHge/c6hTAQsZvuwIFszrhxjp8TczJEGJTi4EGebK1UiYdOnAqDMFStSjRypNrW0PHjHDBVvjz7yCvFk0/y8N3Vq8qVSUQ8OQIQrVlDWq2WAFDW4ME8PqZwDW5w1vn2W2XK0+uJPvqIy3ziCXbfVRW9nrstI0bk35dOIAxEbNrs2ew52KMHOzCUUEQYlGDtWh6lad7cSaMsDQ/gwIFcG9++rbZFlJxM1LUrj+/On297ecuX8ykuXWp7WfcwbhyPv6Sn5wnDrXXr+IDR0YofbtAg7lXZOqSUkcHDmQDRhx86SWT6tm1s0F9/OZ0wGNi8mX/uhg2JzpxR2xpVEGGwlaVLuXLr25eDaZwSwwO4Zg0P8i9erLZFRMSOPmPGsGmffGJ9OdnZPATTt68dxs0zMnh86vXXiYjyhEF78yZRgwY8060wyckcGT1woPVlpKez8JYqRWRjyIWy9O/PNa4hAZcTCgMRx3LUrcsjXmrG4qiECIMtfPcddzsHDFAg0Z09KfgAdu9OFBmptkV56PXcmrXFXdMw4WwXr5Lvv+fC79QOecKg1XKqUS8vu+RamD6dyNeX560s5dYtnuvy9+esp07DlSt8vebN4/+dWBiI+Gdt1IioShWiI0fUtsahiDBYy08/ceN76FAn6aIXR8EHcP16/ttZclbcYdo0Nuujjyz/bqdORA88oLxNpNdznqkCbr5GwpCSwrXvxImKHzoxkYVhxgzLvnf7NgerlSnDaUqcivff5+tlGLx3cmEgYq/CZs1YHEpQz0GEwRpWruTho8GDXUAUiIwfwNxcTl7fqZNzJOkpwMcfs5kzZ5r/HV5XwQ7uqUS8IAVAtHVr3ltGwkBENHYs18J28EseMIDnac29x7KyiLp147rX6Rwgrl3j6zR2bP57LiAMRPzThofzsJKrpZm3EhEGS9mzh1tyTz/tQi5tdz+Ahl7Db7+pa1chvPsum2buuPjo0USVK9shzUhODo+F3xUUeI8w3LjBM8VjxihsAKe2AtjF1xR6PfdefXyMdMx5MKQ6vXEj/z0XEQYiHlYKC+OhJXvF4DgRIgyWEB/Pk5xt2rjYIiB3P4B6PY+9tGjhdF0evT4/cvnAgeL3zcjgBHl2yfTx9dd8ze5KtHSPMBBxV8fHR/EAAb2eI3L79DG977x55JDIaas4e5aDe6ZMMX7fhYSBiCekAwM5QNNlGoXWIcJgLhkZnNqlenUXzOte2ANoyJ+khK+owqSncwK8mjV5BKIo/vuPT+G//xQ2IDGRA1IK8TgqVBhu3WI3oh49FB+emzSJYzOKK3brVh7afPNNRQ+tDHo9X5fg4Hvd9lxMGIi496bRcKZWN0aEwRz0eh7vLVWKw+ddjqIewBEjOADDCZ21L17kIaIOHYrO4hEVxU4uiq8//fTTnA63kDU/CxUGIg5mAdhVTUEMo35FxcecPs2hKV27OqlnnCGzb2GrErmgMBDxHJjd5rWcAxEGc/j8c74SLptHpagHMDWVF4fo1MnphpSIuFPj7c3D04UxdCgHFSrKihWF/tjz58+nRo0aUf369QsXBiKi55/n5HBXrihmTnw8m7Ny5b2fZWTwYk/16jnpUqGXL/PYy4ABhX/uosJQsKF4+LDa1tgFEQZTnDrFi3s4QSYJ6ynuAfzzT7LYFciBfPYZFZmvrkULoiFDFDzYxYs8hNSnT5FjN0X2GIh4gYOqVXnCWsHme7Vqhc+jjBvHUxuHDil2KOXIyeGl86pVK1q1XFQYiFiUw8PZm9kpe2q2IcJQHDodD2WEhjpBjhlbMPUAjhvHkXpO6KWk0/FSoXf/BhkZPIwUFaXQgQpObBTjelqsMBDlD/gXdMu0kccf56GiguzezT/Z3fO5TsObb/J1KM5FyoWFgYjntpz6N7AeEYbiKK616lKYegBzc3mCMDDQKaN4Tp/mXtvo0fnvHTzIp7RzpwIHsMAVyqQwEOXfOAq5CL33Hjs9GMjIYLdJp22tLl7M52+IcC4KFxcGIqLx43m4082GlEQYisIwhFSwMnJZzHkAtVquberXd5JFJIy5W6R37eL/FXkgJ00ic4MnzBKGgkEFCiTZmzGDpy4MGIaQnDJNw/btbNxLL5n20HIDYXDTISWTdbwHSiijRgHVqgHTpqltiYMICADWrwdSU4FHHgGSk9W2yIjRo4GOHYHhw4GcHCAzk9/387Ox4NmzgcmTgU8+AZ54wmY7AQAaDbBgAdChA9CjB7B7t03F+fnln++RI8DMmcCkSUDjxgrYqiT//gs8/jif94IFfB3cnFKlgMWLgX37gPnz1bbGcZRIYfjzT+CPP7jOKF1abWscSN26wNatwJUrQJcuQGKi2hbl4eEBfP45cPo08M03LA4A4OVlQ6EzZgBvvQW89x4wYYIidubh6wusXQu0aAF07Qr884/VRXl7A9nZ/Pe77wKhocDbbytjpmLs3Ak8+ijQsiWft4+P2hY5jPvuA4YOBaZM4XZVSaDECQMRMH480KYN0Lu32taoQOPGwLZtQEIC8OCDwNmzaluUR7NmwPPPcwNfr+f3DBWmRej1/COPGwd88AHw0UeK2plH6dLA5s1ARARXmuvXW1VMVha3THfu5CI++ojFwmlYv57PLyIC2LSphLWmmEmTgFu3gFmz1LbEMZQ4YVixAti7l4eQSkBPuHCaNAF27OBmeWQk8NdfaluUx//+ByQlcaMUADIyLCwgNZUVf+bM/GEke/7QZcoAGzfy8Fzv3sDUqdz6sICMDBaG8eO5Qf7003ay1VKIeAiud28Who0b+XxLINWrA6+9BsyZA1y7prY1DsCciYhCNpckO5uDhbp3V9sShbF2ki85mX3yPT3Zw8RJsrG+9hoHbN+V+NQ0J0/yBHtAgHmZ6QrBrMnnwtDpOD03QPTMM+weayZjx3J4hFPlPUxPZ08ugCfvrQmQdIPJ54IkJ7OTgBs4rIhXUkEM0fv796tticLY8gDm5PDqZQC7tCoY1Wst166xMHh7m7legU7HOaFKl2blP37c6mNbLQwGli9nd7d69cxeMKFzZ/YkdppM6X//zalG/f2Jfv3V+nLcTBiIiKZO5fvS1mVZVUa8kgoyfz7QvTvPFwp38PIC5s7lsZvYWJ6D+P57i4dDlKRyZZ7sIwJiYkzsfO4cT6SPHg0MGMDn0LChQ+wslKeeAvbvBypWZDerN94Abt8ucne9HtizB9BqeShJ1eHN27fZ3o4dgUqV2BXnySdVNMj5GDmS/Q6++kptS+yMOepRyOZy7NnDjZf169W2xA4o1TJLSiJ64QUuq1MnDr9ViePH2YwqVYrY4eZNjgzz9+d8UDYuUmBWriRLyM0lmjWLE+7UqsVLiBaSy/nkST7PatVUTGeVm8v21arF9s6apUzeaTfsMRARvfIKD/1lZ6ttidXIUJKBwYP5vnfLPOtKP4AbN3L2NoCoXz+bhmZsoVEjNiElpcCbGRlEs2dzrmo/P44GS01V7Jg2DyXdzYkTnJsJ4DUkN2wwGi/64gv+aPJkZQ5nEXo929O0KRvRt6+ykfFuKgyHDvFp/fKL2pZYjQwlAezlsmwZMGIE4OmptjUuwGOPAQcOAN99x0Mz4eFAz57smmnwI3UAw4fz6/LlAM6fB155hceZ3noLCAoCOncGTp3i9996C/j5Z+DECYfaeA+ZmTw2tHAh2/Xhhzw+9NBD7M7y+OMcTzJrFnDrFn78kb82apQDbUxLA774gsdUH3+cr+W//wKrVgH16zvQENekaVOO8YuKUtsSO2KOehSyuRSzZ3MUf3GLwrg09myZZWbyamctW/IxQkN5RbOjR+0+U5oRn0wf4z265nnHZcewBQdz5r1HH+XtwQfZLsPnZcvyRPr69RZ3Ea3qMWi1RAsWcN4ELy+2wcuLewidO/NCzQ8/TNS2La+LfMdOvUZDB9CM3qrynf3zLWRnc76RkSP5+nh4cE9m61b7/Y5u2mMgIlq6lE/t6FG1LbEKGUoi4qU6+/VT2wo74ogHUK/nOYcBA/J9SevWJXrjDXYNVSL/UlYWn8P06ZzT6U4Fmo5SlNm9N9Hq1cUfJymJaMsWdh2JiODv16rF/5u5oIFFwnDsGA84lynDLr99+nA62D17il5dSK/nrIELFtCVGvdRLjzyhaRDB6Iff+TPba2s9XpOCPbTT0TPPcd+lgZRff99TkFub9xYGLKy2CtalSFA2zFZx2vIOu8T9VxWLOTqVSA4GFiyBBg4UG1r7MS+fRyVGhsLtGpl/+NlZnJQ3Pr1wLp1QHw8vx8SwnaEh/NFr1aNX4OCOJTXwwPIzQXS0zny2rCdO8fncPhwfqizRgO0aoXYR8aj9bQn8OqrGsybZ6GdMTHc31+2DAgMBBYtMhnunpqaisDAQGi1WgQEBBS+U3Y2B35NmcLeO8OG8VajhkXm1akDXL+QAe0Hs+D59ZfA5cv5H5Yrx79l8+ZcruFaVq7MbjGenoBOx2HT16/zb5CQwGUcPMj3glbLZbVoAfTqxVurVo5zfXL0felgnnkGOHPGDM8558P0DWCOehSyuQxff8295sREtS2xI2q2zAwt4OXLid55h6hLF6KQEHb2Ljj8U9QWFMTLtD3/PFH79vxex468wDzxSJCvLy9vabXXzpUrRD17ctnPP889iyIw2WPYv5/t9fIi+uCDotclNUFcHJvTvn2BN3//nfNv+/kR9erFPZB69fJ7aKY2QxxH37483Ld5s7qZdN24x0DEnTuAF7JzMUzW8bakKHMJ1q0D2rVjt3LBDmg0PJlaty778BvQ63nWPyGBM7nm5vJ73t6c/6FaNaBqVf773DlORJeQwC38l1/m3gW4YdylC6fo2bKFMzNYTHAwx2n8+CPnNdi+nbMohodbVs4PPwBDhvD39uzh/BVWMmMGv44cWeDNrl2BY8d4Iv2rrzgb7OHD3ENIS+NewfXrnMokN5djULy9uRcRHAyULWu1PYLldO/O9+eGDXzLuhXmqEchm0tw+zY3vsyKnnVlXLllFhfH495hYdzzKIRVq/j0GjZUwHf88mWeFK5QodDrVWSPISqKjRgyxOpegoG4OO5QeXgU42m7ejV3lbp25RvZFXHl+9JMOnUieuwxta2wGJN1vFu7q/79Nyco69FDbUuEQrlwAXj4YaB8ef6x6tYtdLeuXblxfPIkD+3bRPXq3GOoW5e7H8ePm/7Od99x037MGODrr21KOa3TAS++yA39Bx8sppHfpw/w22+ccvWpp/LzkAtOxeOPcxp/q7IAOzFuLQx79/KcY6NGalsi3ENWFj9VPj48RlS1apG7li7Nc7CNGwMff8whFjZRvjzHZFStyuMBaWlF77trFw8fDRvGqUNsnLj97DNe18fHB3jgARM7d+rEsQW//87DS4LT0aYN38pHj6ptibK4vTA40glDsIDJkzkYbc0anm8wQevW3Npu2BAYPFiBBnRQEHtV3bhR9Ko4GRncvL/vPg5Ys/FGOnGC1wwaMgS4eZMddkzy6KM8ITFvHhAdbdPxBeVp0YKnw2Jj1bZEWdxaGGJjzXz4BMcSEwNMn86rnzRtatZXIiKAuDiun48cYV2xmdq1ed2GRYuw9tVXER4ejsjIyPzP33+fh7sWL7Y5ZD4rizUmJISXbgAsuDdfe427F0OGsKuv4DSULs0jEnv3qm2JsritMCQmAhcvcktTcCJyc7nJ37Ilr7BmJhER7NTk4cGiMGWKQhkuX34ZeOgh9F63Dsf27kXMHad0z717eejoo49sztaq0wEvvMBu/d99x8JWpQpPd5iFhweLU0ICi5XgVBhCNdwJtxUGww8lPQYnY/16dsn84guLFnRu0oTH5WNjeV3k0aO5Tl+2zEZ7PDw48O3SJWDp0ry3fT79lAXhzTdtKl6v55xPq1dzzqe2bbl1GRFh4chUWBivW/3FF0BKik02CcrSujVw6JB7TUC7rTCcOMHu33XqqG2JYERUFNeOFnblfHw4v9uJE1yhfvYZR7I//zw3pm0iLIwTBy5YABChOgCvTZtYfWwYQsrJAQYNYvuWLMkPuj5xgifSLebll7n7sWSJ1TYJytO4MYvChQtqW6IcbisMCQkc8+Phtmfogpw4AWzdeldUl/kEB/PvCvDv+u233BofMoTnZ21KqjpyJLB/PzxjYzEMAPz8ePzHSrRajk9bvpx7NYaiiPgczB5GKkjlykD//jzRomYGWcGI4GB+Ndyb7oDbVpsJCWY5uwiO5JtvgAoVrF4VrKAwACwOUVHs6TNuHIdEnDtnpW2PPgqEhsJ7yRK8BCDnqaeAonIlmeD333noa9s2DrguGBCu1XKqKavvzZEjOdX4jh1WFiAojQiDCxEfn/+DCU7Czp0crVaqlFVfr1YtP1+fAY2GYxu2bgXOnmUnp6goKxrUnp7AY4/B8++/UR1AbvfuFtun1QIvvQR068bTE0eOcJhEQQz2W31vtmvHcRg7d1pZgKA0ZcsC/v733puujNsKg/QYnIzcXI5Ms8FNzNBjKCwhcJcunFZowABe9KZLFx65sojWreFxZ6BYZ8HC4EScL6dJEx46WrSIUzHVqnXvvoZWpdX35p2ss27nBuPCaDT39mZdHREGwTHExXHAmA1uYtWq8SRfcnLhn5cty8PvW7fykFLDhtx6X7eO52xNEhEBDRGuA6DKlU3urtUC8+fz5GPPnvm9hOHDi/Y4slkY7tjpdo7zLk5hvVlXxm2F4dYtTochOAn79/OrDRlJDUP+pmK8unThFEhLlrBnZ+/e7J32ySe8PkeRNGoE8vDAtWJ2IeKOz4gRPIH8+ussDH/9VXQvoSC3brGXrr9/8fsVS0QEr7uQmGhDIYKSBAS4V+yh26bdzsmxyE1esDdJSRwmauWELpD/e5qTDsPPj11FBw3ixvXChRyr9t57QGgo160RETyy1aoVZ8iAlxfg6YmsOxMUROyCGBtrvCUl8dDB22/znIIlHkaK3JeGvFLJybxQkKA6Xl7ulefQ4ltUo9FotIaVoZwUvR4gCkB2dgZSU93o1yqKW7fyX1NT1bWlCHxSUuDj64tbNtiXleUJoDRSUm4hNdX82eX69TmIeeJEYOtWLxw44Ik9ezyxYYMnMjN5zMfDg6DRAOk6DYBWCAoqC72eQMSfBwXpERGhw/DhOkRE6NGpUy68vbl8S07p1i0feHr6IjW1mMR9JvDQ6VAGwK0bN6B35vFSF7gvlYLIDxkZGqSm3lbbFJMEBgYGAEgjKmb5TnNycxfcAASAl/Z08o0IGOIEdsgGgN4GKNnmcjrc+V3rWfHdUgQMJGAtAVfulEMEJBNwkIBdBPxDmfCiGDQg4B8CdhNwnIDbd/bVEXCMgMUEdLLyHMYQkGbTdWhzx/jGTvC7ymbYVhCwyQnsMHsLKK6et3jN5zs9hkKba6mpqQgJCcGlS5eKXi/3LiIjI/Py0yi5f1BQWcycmYkvvmhhdvmW2OJM54oDBzi5f3Q0p3tU2BYlztV70SKUmjgRadeuFToza44927d7onfv0ggN7YIDB1abZcfZsxo89NAvAF5ESooHHnggF/fdp0PLljq0aKFD5cqZyM7Oyts/ODQU+3Nz4X/sGKrfGSPKzQXi4jxw4IAnDhzwRHS0J06e9ET9+jq89FI2Fi16APv2bTPLnnr15uLmzQ+QmGhej6Gw6+IZHY3SvXohbd8+0F1rWDjTPfxC06b48eJFs+9Luz0fsPxcLbUlOHgvOnbsgGXLMszaX81zDQwMDISJHoPFQ0nFdj/uEBAQYPaN5unpafa+luzv5QV4eflZVL6ltgDOca4oUyb/1VnPtUEDICsLAampnGLUCnsM6+N4epLJfbdt42jo334DPDyexBtveGDECCAszAvGt32BlXLS0kC5uSgNoEzZskbHaNeONwAg4rouKsoT777rB53ub7z9tj8mTOCErcXh4UHIydGgbNkAs3IlFXpd4uMBT0+UbdiQJ1NM7W8Ce93D3h4eCADMvi/t9nwUwNxztbxsb5Qq5Y2AAG+z9lbzXInI5Lie6l5Jo0aNssv+FSpwqn1LyrfUFkux17lag8PP1eCmWoT/vTn23LjBryNHPl3kPqmpeQlTkZjIeYpmzVqKWbM4JVKx7N8PDYAQoNgIOY2G19D55ReenH7sscNYt47jGExlq+jRow2IzM+DV+h1iY1lV6i7RKHI/RXEkvKferro38nWsq3Z355lV6hQHxUq2K98e/+u92DpHAMVs+ZzkevlqkBEBNGwYfYr35nO1d5r6ypyrno9UZUqRO+9Z3URM2YQBQQU/fkffxDVrElUpgzRF1/wIS1i9mzS+/gQAZSwbZtFX9VqiYYP55/hoYeIzp0rfL+dO3mfI0cstK0grVoRvfiiDQUwdr+HnWjNZ3ufa/XqRBMn2qVoizHjXB275rOvry8mTZoEX19fJYu1iuBg+wacONO52htFzlWjYd/QPXusLqKoNCdpadxL6NoVqFePI6BfftmKBddiYpB7J+2p/5EjFn01IIAjnrdsAU6fzu893D3warDf6nszI4NPUIF88nIPK4NOx/ExzpKCR5FzNUc9CtmcnuHDuWFVInCillmxzJ9P5OlJdPmyVV9/+mmizp2N37t2jahFCxt6CQZSUoj8/SnjvfdoF0A5dx/IAgr2HoYPJ8rNzf8sI4Pf/+47KwtfsoQLOH3aavschqvclzZy9Sqf5po1altiNo7tMTgT7hai7hYMGMAJ9Kxceu3uNCcXLwIdOgDXrgH//mtlL8HA998D2dnIGTgQUQC8tm0DTp60qihD72HJEuDrr3nNCMMiLqVKAeXK2XBvRkVxno+7vJEE9VAkzYmT4bbCUL06VxiZmWpbIuQREMDi8OWXVoWJXryYH2WckAB07swV7t9/89CN1RBxhduvH6hqVfwKQB8UxGNBNjBoEPDrr7x62wsv5Odrql6dz8Vi9u7loTgr17MQ7IPht7RqjQ0nxW2FoXlzft4PHVLbEsGIkSO5Vv/+e4u+lpICnD8PNGvGKSkeeYRFYft2BRrPa9ZwKtZXXgEAZAHIGTCAVwIqNrmSafr144yrq1Zxj4aI701D6iiLmDoVqFmTV5sTnIb9+zmlirPMMSiB2wpDs2YcyyDZiZ2Mpk15Tc6xY4ErV8z+2r59/BoRATzzDPcGt2wxnbTOJMnJLFaPP85BgnfIHjOG14Z95ZV7Z5AtpE8f1phvvgHmzOFzOHiQA+fMZsUKVpfp021ablRQnthY9quwehjTCbGbMPTq1Qs1a9ZEqVKlUK1aNQwYMADxDhz0L1WKXb3tLQznz5/H0KFDERoaCj8/P9StWxeTJk1CtjutDF6AKVOmoF27dvD390e5cuWsK+TTTzm96PDhZle6sbEcJ7VtG6fV/uknTnNtM2PG8HjjokVGTzZVqMBDSWvW8NqcNjJwIPDGG8D773Peu4wMzgBrFomJLF79+gEWxgYUxo4dO9CzZ08EBwdDo9FgzZo1NpfpjEydOhWRkZEoW7YsKleujD59+uCExYt0mCY2VhEnMZtYuHAhmjVrlhfU1rZtW2zevNnq8uwmDJ07d8Yvv/yCEydOYOXKlThz5gyetHJJR2tp3dr+whAXFwe9Xo9Fixbh6NGjmDt3Lr744gu8++679j2wSmRnZ6N///545c6wi1WUL8/zDJs28eysGezdCzRqxEt4DhvGrqk2s3Il8OOPwGefYcHq1QgPD0dkZGT+5088wRXx6NHApUs2H+7jjznoe/58/t+se1Ov5xzfej3PgyjQLE1PT0fz5s2xYMECm8tyZqKjozFq1Cjs3r0bW7ZsQU5ODrp27Yp0BfNjJySwI4HawlCjRg1MmzYNsbGx2Lt3Lx566CH07t0bR48eta5Ac1yXCtksZu3ataTRaCg7O9sGLyvLiIoi8vJiF0FHMmPGDAoNDXXcAVVwC1y8eDEFBgbaVsiIEey+um6dyV1DQ4lCQnhTJEZp+3aiUqWI+vc38nG9JzgoMZGoVi2ihg2Jrl+3+bA7dxJpNESVKhGNHm1iZ72e6NVX+QurVtl87MIAQKtXr7ZL2c7mrnr9+nUCQNHR0YqVuX49n2JRAY1qUr58efr6668L+8g53FWTk5Px008/oV27dvD2Ni+XiBK0bcvjuP/847BDAgC0Wi2CgoIce1BX5PPPeRWd/v2BtWuL3O3iRV6R7dIl7mDYsKQDs20bzyk88ADwww/Ft8IrVuSxq5QUzrNh4/qN7dvz4j5JSZzHqUj0euDNN/kaLVwI9O1r03EFfi4BKPpsbt/Oy2PYPNelIDqdDsuWLUN6ejratm1rXSHmqEchm1m888475O/vTwCoTZs2dOPGDauUz1r0eg5VHzPGccc8deoUBQQE0Jdffum4g7pqj4GIKCuL6Mknuefw2WdEOt09u8yZw6f3wgs2HkuvJ1q8mMjXl6hrV6Jbt+7Zpch0AseO8c1Upw7R3r02mZGeTlS1Kp/T2bOF7HDzJtFzz3FPISrKpmOZAiWkx6DT6ahHjx7Uvn17RcutV4/opZcULdJqDh06RKVLlyZPT08KDAykjRs3FrWryTreouyqGo1mGoBxxe1z/PhxNLwzK/j2229j6NChuHDhAiZPnoyBAwdiw4YN0Dho+l6jAXr14jV/5861bHh2/PjxmD59erH7FDxXALhy5Qq6deuG/v37Y9iwYdaa7XCsOdciuXEDOHWKB17j47mFnZLCcQtE7Crm58fNrOBgjgqaPJn/HjOGnf6/+YbX4rzDt9/y66RJNpzklSvsL7pxI/Dii9wKtyRlQKNGwM6dPAF8//3A+PG88o8VaQf8/YG33uLthx+ADz4o8OFvv/EkilbLEXJNmgAbNuQPZl+/ztdSp2PvJB8fns02XMvgYJ6Vt2ntUPdj1KhROHLkCHbu3KlYmSdO8K0+e7ZiRdpEgwYNcODAAWi1WqxYsQKDBg1CdHQ0wsPDLS7LovUYNBpNJQAVjh8v2p+iTp068DHkRy7A5cuXERISgl27dlnfvbGCzZvZ7fvIEfZSMpfExEQkJSUVu0/Bc42Pj0enTp3Qpk0bLFmyBB4eDvQE3rePZ79iY3mdSgux9FwBYMmSJRg7ZgySNmwAduzIX/eyYOSWjw9XVEFByFvuTKfjxXGvXjVOMVq6NOesvnQJyMriNThffx2pVBblynFdd+yYxacG3L7NMRPjx3NluWgR0LNnkbunpqYiMDAQWq228DTHOTnAtGm8TmiDBtzi6NLF4knh1FS+LDVrAmfPgoM0xo5ll9Rq1YCyZYEzZ/Kj4jw8gMqVeStViv/X69mj6vp13gxpXT09WcgMa5c+8gjbWggajQarV69Gnz59LLLfLGy8L5Vi9OjRWLt2LXbs2IHQ0FDFyp05kxsrN244pw4//PDDqFu3LhYtWnT3RyZvVot6DESUCMCqFcj1d27arKwsE3sqS+fOXOesW2eZMFSqVAmVzFxP98qVK+jcuTMiIiKwePFix4qCAlhyrrh9G9i8GQ989RVOpqUBHTsCgYFcATz9NL+Gh+cLQnEVZkYGC8S5c1yJxMZyRXfmDDejJ09GavX70I3eQ783e1h2UqdOAV98wXm3b97kiOtPP2WPKFvw9uaeQu/ewJAh+ZXuK69wqLOZLrwBAUD7+3MRuWsudLW/gOeFs/yBpycLQ+vW3KVo3hyoUYMFobjFonNzWRwuX+Ygib17+XouXcqRgPXrc/e5Vy+e6HCxe9QaiAivvvoqVq9eje3btysqCgDXKY884pyiAHCda3V9a854UyFbsezevZs+//xz2r9/P50/f57+/PNPateuHdWtW5cyMzOtG0CzgX79OA23Pbh8+TKFhYVRly5d6PLly5SQkJC3OQxHjOXGxRGNGUO6gAAigK5WqULTvb3pxA8/0P7YWEpLS1PuWAkJnGO7fn3S8wAU6UuXJurUiWjuXKJdu4ji4znxXUoK7//ffzweP3QoZ9UDiIKCiN5+26KEcxalZ9briaKjObuflxeRnx9Rx45Eb75J9NNPPC9x7Rq7Ud24QXThAtHmzUTvvkvUtCnpPDyJAMrVeBF16ED0yy88AaEk6ens9fXSS3kTG7q6deny2LF0KDqaANCcOXNo//79dOHCBWWPrfIcwyuvvEKBgYG0fft2o+fy9u3bNpedkEDk4UH01VcKGKoA48ePp+joaDp37hwdOnSIxo8fTxqNhv7444/CdjdZx9tFGA4dOkSdO3emoKAg8vX1pdq1a9OIESPospVZNW1lzRo+UxvnDAtl8eLFRa6r6jDs+QBu307UpQuXX7EirW/ShOoWcq7bLFy/wBzOnSMKQAqtqfkqr+VwRySK3Ly8WBSGDiX64Qer/JStztufkEA0axZPpIeGmrZVoyFq1IheL7WQygfkmi5fCXQ6oh076GqXLpQF0G2AvgKo5p3fcNCgQcoeT2VhKOq5XLx4sc1lf/QRtwNSUmwuShGGDBlCtWrVIh8fH6pUqRJ16dKlKFEgUksYnI2cHPZ/HzpUbUvshD0ewP37ibp143JbtSL68UciB/f2Rozgw+dpzsGD/CN6eRGVK8crMS1dSrRhA9GePYoErCi2oEtSEhu+di03K3v0YLsrVCCaMoV7EkT07LN8jvv322y6ZVy9SvTJJyy4Pj7suqdAnIYRTuSVpCQ5OUQ1ajiPN5IViDAY+PhjVvjkZLUtsQNKPoA3bhANHMjl1a9P9OuvNixyYD25uTwS5OFBdE9M5NmzRAMGcKs7LIzo998VO66iK33l5hLNnElUtiwvPffxx/e4yG7YwJf6mWdsP5xVpKVx8zcggO2cOdN4AQlbcFNhWL3a5U9LhMFAQgKRtzcPUbsdSj2Aq1YRVa5MVL48r3qTk6OMfVZgqDDDw4vZ6eBBXkcT4N6DDZX5/PnzqVGjRlS/fn1lhOHYMaL772fxeu01FtxCSExk8319FYrotpbERLZToyFq04bo+HHby3RTYXjkEb5ELowIQ0GeeYYDUpRqEDkNtj6AN2/mj2n07s0qqjLdunHGipdfNrGjXs8iVqYMjxf++adNx7W5x6DTcavb15d7XLt2mfxKjRpcH8+bZ90hFeWff9huX1+eM7Glt+iGwhAXx6f0/fdqW2ITIgwF2b3bLX7Ue7HlATxxgqhBA6LAQJ5HUGHY6G527aK8+VmzvT7OneN1Pz08iD791OrzsEkYbt3iyWeAaOxYIjO9X/r146H+4GDlnZKs4vZtth/gXFKFRIibhRsKwzPP8O/k6PxrCmOyjnd/Z+YC3H8/58b/4AOOoSrx/PYbcN99HGsQE8NrUKqcVJ6IY9HCwvhvs1dmq12bF2h4801ORjR0qGN/5AsXOD5g82YOUps1iyO8zaBpU45Nu349P/Oqqvj5sf2rVnEG3Ace4PMr4ezfzxnYP/yQYwzdmRIlDADwySccnHtvMGAJY9EioEcPrsx27wbq1VPbIgCsVTt2cPYKgAOAzcbTk8NRv/8e+Pln4OGHObWEvdm3D4iM5GP9+6/FCe/KlOEYtOHDeZG2ggHhqtK3L7BrFwcIRkZaueyc+zBhAscyvvii2pbYnxInDI0aAYMHc278tDS1rVGJzz7jHP8jR3L4ZmCg2hYB4FbzhAlAhw5Ay5b8npmNbmMGDOC0l0eOsDjYs6bds4dTYtSuzb2upk0tLsLPjwO+J05kgZgxQ3kzraZZMz7H2rU5u+yePWpbpArbtgG//w5MmVJ8ALq7UOKEAeCuYGoq95ZLHAsW8FDLO+8A8+Y51TKRS5dyNodp0/JTBFltXps2/DSfOwc8+qh9eg6xsbxiUHg4p+auWNGqYry8OAVT1aq8yttnn3FmC6ehUiUepgsP5/M1rLNaQtDreXgzMpJzKJYESqQwhITwUPS0aVYmZnNVfv6ZVyMznLwTLVJ74wab9cQTQLt2+WO4Nk0TtGjBFdqpU5wjSMnlVk+d4kqyYUOeV7BhkYjMzPye0dtvcwduxAibl5pWlsBAPs8GDThB0OnTalvkMKKiuKM0a5ZTPTJ2pUQKA8AT0HXq8LCSRYuyuyoxMZz0beBAp7zDR4/m38Ew+WoQhowMGwtu2ZLTVu/ezQdRorZNTeUkehUr2iwKAJ+j4XwDA3n6Z+NGTsntVAQE8CRQhQostKmpaltkd86c4eVkR47kfJElhRIrDKVKceLN2NgSMKQUH8/uWC1b3rPovTOwciWwfDmLQtWq/J5hVMbGBdOY9u15/YWvvuLmny3odMBzz/E1XbfO9myt4HOsUCH//169gBde4OUp4uNtLl5Zypfn875yhb3YDGN+bohez85tlSsDJpYrcT/M8WktZHMb3n6bU8UcOaK2JTZQnL94djaHaVavzhlJnYzERA627tPHOPRAr+fQio8/VvBgr7/OK8UVs+avyTiG99/nWIlNmxQzq337e1NiJCVxMtQePZwitOReNm3iQJOJE4vex8XjGD7/nM3/6y+1LVEcCXAzRUYGr/MeEeHCQSvFPYAff8yV4e7djrfLBHo9B3cFBRUebN25M1HfvgoeMCeH01uHhnKOoEIoVhh272ZR+OgjxUzKzSXy9+dg6btZu5Z/VkeuEmsR//sfX489ewr/3IWF4fhx/l1GjlTbErsgwmAOe/dyBoCBA520dWaKoh7AQ4c4QdSECerYZYKPP2azV60q/PO33iKqWVPhg546xdkUR40yettkriRDC6J1a0VzSB05UnyrdNgw7tGakVnD8eTkcIsqPLzwVpWLCkNyMqfOCQ8vsv3g6ogwmMuPP/LVmD1bbUusoLAHMDub02U3buzwdNnmYFgj48MPi95n2TLeJzFR4YN/9lmRtXGRPYZx4+wy5vjdd2zKzZuFf56VRfTAA5wy49IlRQ+tDIcP83UZP/7ez1xQGHJyiLp25TySFqzv5GqIMFjCuHHcM968WW1LLKSwBzAqiseAi+rmq8jhw5zz7oknOOdcUZw6xae1YYPCBuh0XNuGh9+TUbFQYThxgofjFJ3wYEaN4tZpcVy7xvkBIyKcJJfS3RiGK0+eNH7fBYXhzTf5VLZuVdsSuyLCYAm5uTzZFxjIWRRdhrsfwFu3eOZy4EB17SqExEQe4m/WzHQ3Xa/n0ZunnrKDIXv28DW7azWvQoWhf3+umRWehMrI4HV7Xn/d9L779/OY9zPPOOFwZ0YGp4h9+mnj911MGBYvZnOdIsutfRFhsBStlqhRIx7bPn9ebWvM5O4HcMoU7t6fO6eqWXdz8yYP0VesaL5pn33GC5/ZxaHqySf5hy5Q4d8jDDExfG2//Vbxw3//PRd94oR5+//6K+8/ZowTisM337BxBdfPdSFhWLuW77OXXnLCa6s8IgzWcOkSt2rr1iVSaZlqyyj4AKak8GpcY8aobZURaWlE7drxipz79pn/vZQUbin/7392MCoujscNPv007617hOHRR7mlYIdFi9q04UVfLCEqin/q8eOdrALLyeHu3aOP5r/nIsLw22/cjnriCVXXpnIkIgzWcu4cjx6EhRFdvKi2NSYo+ADOmcOeSE6w2I4BrZZ99cuWJfrvP8u/P2wYh2HY5aF94QWiOnXyJjuMhOHoUb6uP/2k+GH37eOiV6+2/Ltz5jipOBg8OI4d4/9dQBg2bWKPxMcf54n+EoIIgy2cOUNUqxZR7dq8zLDTYngAY2J4JvPZZ9W2KI/kZKL77uN5G2tDKfbv59NbsUJJy+7w779c+J2ANSNhGD2ao+/sUGO89BIPy1srdgZxcKphpcxMokqVeIlQIqcXhjVruA3Vq5dTOu7ZExEGW7lwgYeUqlRxUl9yovwHcMECfv37b7UtIiJ2UmnYkAPYbK0bOnfmhr21i4kViV5P1LIlNxmpgDBcucJdnPfeU/iArN+enoUHtVmC4ed++mkn8laaMIGHMm/dclph0Ot57srTk/0KsrPVtsjhiDAowbVrPBTi43OPE4tzYHgAu3QhatrUKZqQf/zB8wkNGijj4XXyJMelvfqq7WXdw1dfsWvvpUt5wnD788/Zd/nCBUUPlZnJoSWtWilTIf36K8/BtGrlJEOe58/zdfv2W6cUhsxMoqFD2aw33ywxcwp3I8KgFFlZ3P0HiN54w8luKMMD6OdHNHWqqqbo9TyX6+FB1K0bTx4rxdy5fJrbtytXJhGxkV5eRFFRecKQ3a0bUceOCh+I6N13efji0CHlyty/n52rqlQh+ucf5cq1mg4deHzGyYTh6tX8Bt6SJWpboyoiDEqi13NiLU9Pjo68dk1ti+5geAABVbMB3rpF9OKLbMZbb90TO2Yzhrg0uwwpPfQQUbdupNVqyR8gfalSRLNmKXoIwxCSHeLk6No1ro99fIgWLVK50zhjBjdSdu50GmH47z92JqlalaeVSjgiDPZg61b2xa9YkWj5crWtoXxhqF5dtRph+3ausP38OM2DvTAMKQ0bpuyp7njiCcrSaKhlWBj1Mojs3ZG8NqDVcqC1UkNIhZGVRTRiBJvevbuKKTTi4tgIQxdPRWHIzGTvLQ8PdoJwCfdz+yPCYC+uXWO/Z4DjpFTtPRiE4e7czQ7g1i0e9we4Na9gXVokX39NJvMsWcyZM0QApS9dSgsByg0LU6zo27eJHnyQPbOOHlWs2CLZsIEoOJiP9+23KrUV6tXLf0BUEoY9e1iMvb25l1YCJ5mLQoTB3ixfzmkNKlYkWrpUpYdwyxb+KR08v/DXX/m9hE8/LT7vkdJ88gmf8pw5ChWo1xNVqUKZb79NewHKev55RYrNzGSHJ8PIiqNITiYaNCi/96DwHLppBg3iwEAVhCE9Pb+XEBHBubkEI0QYHEHB3kObNnaYHDXFvHl88DVrHHK4Q4e4snNkL+Fu9Hr2jAQ4KloRQX7sMcp56CHKAui2AvMLaWlEDz9MVKoUR9eqwfr1RNWqcRDX2LFEN2446MDz5nFT3YHCkJ3N8yvBwdJLMIEIgyP54w8eQza00g4ccNCBX3mF7slTYwfOnSMaMIA9O+vW5R6SI3sJd6PXc1oogGPRbA5SmjiRdOXKEQF0y8b0mhcvciOhTBmibdtstMtGUlOJJk1iWwICuMJUfPL+bv75h/IcIuwsDHo90S+/ENWvz4d79lm3TpmtBCIMjkan4+GlevW4An32WQdkvu7Sxa4P4IkTnB7a25tdIqOinCt9QFQU2xYebuO1Xr2aCCAdQNqrV60qQq/nOZCAAG65OlPW82vXOCjZ8DvOns1DTnYhPZ0fADvel9nZHA3fujUfpls3dt0VTCLCoBaGbm3NmnyVIyN5IvD2bTscrGlTxR/AnBxeWe3hh7noihUd1NK0ksOHeTzZw4PHl63qPdxJYJRY1NKeJrh4kSsngGjwYDtWujZy7hxPAXh789zH0KF26mxWqGAXYbhyhR0PgoO5+A4d1O+VuRgiDGqTm0u0bl1+hVG+PEdcxsYqOFFdtapiD2BcHNHkyez5ChC1bUv0ww+usR52djaLl6H38NdfFl7jq1eJADploTBkZBAtXJjfS1B8YSE7kZDA1yskhH/r++7j7NmKrZjXoIFi92VGBtHGjewB6OlJVLo00csvO3C41r0QYXAmTp3iwK+gIMoLO3jlFV4xzurxcb2eo3atfABzcoh27GC7DGO0/v4c5e0EcUkW8fHHH1Pbtm3J1zeCPD1jCeD0EwsWcByBSXQ60ms0dNBMYTh7luidd/IbxoMGOW8voThycthvoWtXPg8PD26Fz5xp/loRhdK+vU3CkJjIEcr9+rEQAJx7a968opdCFczCZB2vISJYgVVfEpicHODvv4F163g7dw4oUwbo3BmIjAQiInirUsWMwpKTgQoV+O/YWKBVq2J3v3kT2LePd927F/jzTyApiY/VsydvDz8M+PvbfJoOZ9KkSShXrhwuX76Mr7/+BqtW3URUFLB2LeDnBwwYAAwdCjRvDnh5FV4G+fhgf04OwrRaBAQE3PP5rVvA9u3AF18AmzYBgYHAiy8CI0YA9evb9/wcQUICsHEj35dbtgCZmXxeDz7I92Tr1kCTJoCvrxmF9enDF9+M+1KvB86c4V1jY4Fdu4Ddu3n2+v77+b7s1Qto3BjQaBQ51ZKMySsowqAyRMDRo8D69VxJ79sHpKTwZzVq8MNYrx4QHAxUq5b/WrUqV94eV+OhqVEdAKDbE4usxq1w7RoQH88PeUIC/332LD9wZ85w2aVLAy1bAh078gMXGQl4eKh0ERRmyZIleP3113Hz5k0AwOXLwFdfAV9+CVy9yiLRvHm+ALduDdSsye97l/VFbHY2Qq5q4eUVgOPH8yur2Fjg+HH+zVq2BEaNAp591jVF1Bxu3+Z7cv16rqSPHQN0OsDbG2jaFGjRAggJufferFiRhdfjuWeA5ctBe2Oha94KWm3+/Wi4N69cAY4c4fs+NZWPW6sW/ybduwM9evC9LiiKY4Xh9m0gLs6a4gQDRPywHD/OW1wcP0iJidx6u5sQXMRF1EI3bMLv6H7P5z4+/KBWrQo0aACEhwONGnFF6OnpgBNSgXXr1mH27NmIjo42ej8nBzh0KP/aHj8OXLxIIMp/TjLhi8NogkjE5r3n48Pi3KgRb40bA2FhJa/lmpkJnDqVf+1OneL7MjmZW/x38ynG4GFsRRMcLbS8wEC+N2vVyr+2DRsC5cvb+URKCA0bFtlocaww7NvHLTDBcVTHZVxGCDahG17BQlxEbbVNcmlYGMoiEpcA+KltjkvzLV5EL6xDS+zDJdRS25wSRzEjeCaFoYiRVuto2JCNERyH1w0v4FHgMfyG335MRkaj2mqbpCjz5s3Dd98tKXafFStWIjQ0NO//onoMhZGdnY3s7Oy8/727EDS5SYiOzkGZMiIMthD6ThrK/5mM339MQkYjEQZH07Ch9d9VVBj8/U3OMQlKk5Xf727UCICbXf+ZM5/F+PFdi92nTp3q8PHJ///QoWR4eh408170ubMxRHp4gsfPC5l7FixBowXgnvelu6OoMAgq4OvLg7VardqW2IVKlSqhUqVKjjlYWho0Oh3McbgRzCAxUW0LBCsRYXAHKlVyW2GwhIsXLyI5ORkXL16ETqfDgQMHAABhYWEoU6aM6QISEgAA0lFQiBs31LZAsBIRBnegYkXg9Gm1rVCdDz74AN99913e/y1btgQAbNu2DZ06dTJdQHw8ACAIQI51ThmCASIRBhfGTTzXSzi1ZGIP4PiFwqI4zRIFADh6FOTpCT8AmgsX7Gmq+3P2LJCVpbYVgpWIMLgDjRrx661b6trh6sTGQn/HlcNz/36VjXFxxD3RpRFhcAcMwnDihLp2uDqxsdC1bo3LADzvzE8IVhIba2ZOF8EZEWFwB2rX5tejhUeYCmZw+zZw9Ch0LVogFoDnvn1qW+TaxMTkN1gEl0OEwR0wZITbvVtdO1yZP/8EdDroHngA2wF47t4tQ3PWkpYG/POPpEFwYUQY3InYWHFbtZKj06bhnI8PIp5/HusAaLKzOb2oYDlbtgDZ2ZyhUXBJRBjcidxc4Pff1bbC9dDr0fjsWYS+9hpiYmJwFoCuYUPOPS1Yzrp1nJu7Rg21LRGsRITBnWjQAFi1Sm0rXI///uN83D175r2V2707sGEDp2QVzCc7m69bgWspuB4iDO5E9+7A6tWSisBSvvySY0Hat897K+eppzhAa80a9exyRVav5pWfnn9ebUsEGxBhcCd69uTVdr79Vm1LXIekJGDZMuCVV4wWqNCHh/MY+YIFKhrngkRF8XJvjRurbYlgAyIM7kS5csAzz/C6kzqd2ta4BosX8yozQ4bc+9nIkUB0tLgBm8uRI8COHXzdBJdGhMHdGDkSOH9e5hrMISsL+Pxz4KmnOBHh3fTty0vfzZ7teNtckdmz+Xr16aO2JYKNiDC4G5GRPNfw/vvspSQUzRdf8ILQ771X+Oc+PsCECcB33/FalkLRHDsGfP898O67MFocQ3BJRBjckU8+AU6e5GESoXBSU4GPPwZefLH4pa5efpkXyC5KPATmvff4Og0frrYlggKIMLgjLVoAzz0HfPghkJ6utjXOyaxZHKH74YfF7+frC/zvf+xt8++/DjHN5fj3X/be+ugjvl6CyyPC4K589BGQnCwt3cI4ehSYPh144w3zgrCeew5o2ZJbw5JK2pisLGDYMF7T97nn1LZGUAgRBnelTh0eUpo3D/j7b7WtcR5yc3n4qE4dYNIk877j6cnDcnFxLLhCPh99lD9s6SHVibsgv6Q789prQLt2XBHevq22Nc7BrFmcU2rJEqBUqby3FyxYgPDwcERGRhb+vebNgYkTgWnTZK0BA3v38vWYOBFo1kxtawQF0ZB1SxjKuofOxL59nMkyNpa79AU5eZIrtaee4spQo1HFRKdg1y6gUyceQpo+vdBdUlNTERgYCK1Wi4CAu1Z/zskB7r+fs67+9x9Qvrz9bXZWUlL4WpQpw9fC2/vefYq7LwU1MVkJSI/B3alfH/j6a3YlnDtXbWvU49Iljkto08b64SBvb+CXXzhVxrPPltwgwtxcDqS8cQP49dfCRUFwaUQYSgLPPw+MGwe8/XbJzL56+zYHXZUqBaxYYZuffVgYi8PWrXxNSyLjxvH6Fb/+CtStq7Y1gh0QYSgpTJnCgW9PP81d/JJCTg57y8TFAWvXApUr217mww8Dc+ZwpG9Jy6W0YAGf+9y5QJcualsj2AkRhpKCpyfw88+cmvuRR4BDh9S2yP7k5gIDBgCbNnErv0UL5cp+9VWeqxg9mofqSgJffcXnazhvwW0RYShJBAQAv/3GKaY7d3bvnkN2Ns8DrFzJ2VN79FC2fI2GewwjR7If/8KFypbvbERFcRzHyJF83iXZiaEEIMJQ0ihfnseHw8JYHDZtUtsi5UlOZiFYt46FoV8/+xxHowHmzwfGjOEKc+JEztTqTuj1nHdr1Cg+z/nzRRRKACIMJZHy5Xny9MEHgccfB2bMAKxzW3Y+jh0D7rsP2L+fe0e9etn3eBoNj7dPncrzOP36caoNdyAtjc/nk0/4/ObOFVEoIYgwlFTKluX8NhMmsJfJgAGuX6GtXs3uqH5+QEwM94gcgUYDjB/PPZS//uKgwpMnHXNse3HyJNC2LZ/PunV8fiIKJQYRhpKMhwe3cpctY5Fo2pQrAlcjOZmFrV8/nljftQsIDXW8HY8/DuzezfmDmjfnFrarxTrodOx11Lw5e3T99x+fl1CiEGEQ2IX10CGuTLt04fFyV+k9rFvHy0iuX8+R3StWcG9ILcLDgQMHOF332LE8XHfihHr2WMKJE7yc6VtvASNG8HBco0ZqWyWogAiDwNSpw5PS8+fzwjR163ICPmfNJrp3L8cT9O7N6RaOHgUGDbJ6uMNkriRL8PcHPv0U2L4duHoVaNKE3TuvXrW9bHtw9SpPLjdpAly/zsuZzp3L5yGUTIjImk1wJmJjiQB+VYKLF4mGDCHy8CCqXZvou++IsrKUKdtWjh4levJJPt/wcKLVq4n0esWK12q1BIC0Wq0yBaanE02bRlSuHJG/P9H77xPduKFM2bZy4wbRe++xXeXKEU2fzvYqhdL3paAUJut46TEI9xISAnzzDXD4MK9DMGgQxz588AEvhelocnJ4iOihh3jY6L//OM3zoUOc6sKZJ0X9/Xly/+xZDoqbNQuoXh0YPBjYs8fx9hDxcQcPZjvmzOEsvGfPAu+8I70EAYAMJQnFER4OrFrFAtGvHw8v1K4N9OzJwnHtmv2OnZMDbNvGUba1agH9+/N7P/8MnDrFFZunp/2OrzTly3OK6osXgcmTeZjp/vs5++jUqTwUZi+XYSIuf+pUoHVrPm50NNtx4QK/X5IzxQr3IGm33QFHpTdOSwN++AFYuhT45x9+7/772ROodWvegoOtKzsjAzh4kM9h506OQbh5k8vr25ejbh2Q87/YtNtKotPxOX77LSc2TE/neZ7HHssXjPr1rRM/nY7dTWNjuXe1cSNw7hxQujTw6KPAkCFAt272F1ZJu+2smOxiizC4A2o8gImJHDW9di2vEHfjBr9ftSr3NIKDgWrV+DUoiFMze3hw/qL0dCAhIX87d44D03Q63q95c64ge/Xi83HgUJHDhKEgmZncO1q3DvjjDx7WAXitg6ZNeflRw7WsXJnXVfb05OuVlcUTxvHxfC0vX+YhNsNa33XqAF278rXs3NlocSK7I8LgrIgwlAjUfgCJeL2D2FjeTp3Kr6ji47k3UBAPD67gDJVdjRo8lxERwRWhigvKqyIMd5OSwr9pbCxw5Ej+tUxI4M/upnx5vpaG69mkCV/LVq3UHSJS+74UisKkMHg5wgrBzdFogJo1eevb1/gzIk5ol5vLeXe8vXlzpfkBR1O+PMeTFJbWOieHr2VuLuDlxZsslCMojAiDYF80Gu4BqNgLcCsMwioIdkS8kgRBEAQjRBgEQRAEI0QYBEEQBCNEGAQBCudKEgQXR4RBEACMGjUKx44dQ0xMjNqmCILqiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIbsH58+cxdOhQhIaGws/PD3Xr1sWkSZOQnZ2ttmmC4HJIEj3BLYiLi4Ner8eiRYsQFhaGI0eOYNiwYUhPT8esWbPUNk8QXAoRBsEt6NatG7p165b3f506dXDixAksXLhQhEEQLESEQXBbtFotgoKCit0nKysLWVlZef+npqba2yxBcHpkjkFwS06fPo3PP/8cL7/8crH7TZ06FYGBgXlbSEiIgywUBOdFhEFwasaPHw+NRlPsFhcXZ/SdK1euoFu3bujfvz+GDRtWbPkTJkyAVqvN2y5dumTP0xEEl0CGkgSnZuzYsRg8eHCx+9SpUyfv7/j4eHTu3Bnt2rXDl19+abJ8X19f+MrqcoJghAiD4NRUqlQJlSpVMmvfK1euoHPnzoiIiMDixYvh4SEdYkGwBhEGwS24cuUKOnXqhFq1amHWrFlITEzM+6xq1aoqWiYIrocIg+AWbNmyBadPn8bp06dRo0YNo8+ISCWrBME1kb624BYMHjwYRFToJgiCZYgwCIIgCEaIMAiCIAhGiDAIgiAIRogwCIIgCEaIMAgCgAULFiA8PByRkZFqmyIIqiPCIAgARo0ahWPHjiEmJkZtUwRBdUQYBEEQBCNEGARBEAQjRBgEQRAEI0QYBEEQBCNEGARBEAQjRBgEQRAEI0QYBEEQBCNEGARBEAQjRBgEQRAEI0QYBEEQBCNEGAQBkitJEAoiwiAIkFxJglAQEQZBEATBCBEGQRAEwQgRBkEQBMEIEQZBEATBCBEGQRAEwQgRBkEQBMEIEQZBEATBCBEGQRAEwQgRBkEQBMEIEQZBEATBCBEGQYDkShKEgogwCAIkV5IgFESEQRAEQTBChEEQBEEwQoRBEARBMEKEQRAEQTBChEEQBEEwQoRBEARBMEKEQRAEQTBChEEQBEEwQoRBEARBMEKEQRAEQTBChEEQILmSBKEgIgyCAMmVJAgFEWEQBEEQjBBhEARBEIwQYRAEQRCMEGEQBEEQjBBhEARBEIwQYRAEQRCMEGEQ3IZevXqhZs2aKFWqFKpVq4YBAwYgPj5ebbMEweUQYRDchs6dO+OXX37BiRMnsHLlSpw5cwZPPvmk2mYJgsvhpbYBgqAUb7zxRt7ftWrVwvjx49GnTx/k5OTA29tbRcsEwbUQYRDckuTkZPz0009o165dsaKQlZWFrKysvP9TU1MdYZ4gODUylCS4FePGjUPp0qVRoUIFXLx4EWvXri12/6lTpyIwMDBvCwkJcZClguC8aIjImu9Z9SXBTty+DcTFAQ0bAv7+alujKOPHj8f06dOL3ef48eNo2LAhAODGjRtITk7GhQsXMHnyZAQGBmLDhg3QaDSFfrewHkNISAi0Wi0CAgKUO5GSiBvfly5O4Q9DwR1EGARnJjExEUlJScXuU6dOHfj4+Nzz/uXLlxESEoJdu3ahbdu2Zh0vNTUVgYGBIgyCO2NSGGSOQXBqKlWqhEqVKln1Xb1eDwBGPQJBEEwjwiC4Bf/99x9iYmLwwAMPoHz58jhz5gwmTpyIunXrmt1bEASBkclnwS3w9/fHqlWr0KVLFzRo0ABDhw5Fs2bNEB0dDV9fX7XNEwSXQuYYBKEAMscglADsNvksCG6JRqMJAKAFEEhEEtQglEhEGAShABr2ay0LII3k4RBKKCIMgiAIghEy+SwIgiAYIcIgCIIgGCHCIAiCIBghwiAIgiAYIcIgCIIgGCHCIAiCIBghwiAIgiAYIcIgCIIgGPF/xuUYL+/4+VUAAAAASUVORK5CYII=\n",
"text/plain": [
"Graphics object consisting of 14 graphics primitives"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"c1 = line([(-3, 1), (3, 1)])\n",
"c2 = line([(-3, -1), (3, -1)])\n",
"c3 = circle((-sqrt(2), 0), 1)\n",
"c4 = circle((sqrt(2), 0), 1)\n",
"c5 = circle((0, 1/2), 1/2)\n",
"c6 = circle((0, -1/2), 1/2)\n",
"\n",
"b1 = line([(-sqrt(2), -3), (-sqrt(2), 3)], rgbcolor=(1, 0, 0))\n",
"b2 = line([(sqrt(2), -3), (sqrt(2), 3)], rgbcolor=(1, 0, 0))\n",
"\n",
"x, y, r = circle_from_points((sqrt(2) / 3, 1 / 3), (0, 1), (sqrt(2), 1))\n",
"print('({}, {})'.format(x, y))\n",
"print(1 / r)\n",
"\n",
"b3 = circle(( x, y), r, rgbcolor=(1, 0, 0))\n",
"b4 = circle((-x, y), r, rgbcolor=(1, 0, 0))\n",
"b5 = circle(( x, -y), r, rgbcolor=(1, 0, 0))\n",
"b6 = circle((-x, -y), r, rgbcolor=(1, 0, 0))\n",
"\n",
"x, y, r = circle_from_points((sqrt(2) / 3, 1 / 3), (0, 0), (sqrt(2) / 3, -1 / 3))\n",
"print('({}, {})'.format(x, y))\n",
"print(1 / r)\n",
"\n",
"b7 = circle(( x, y), r, rgbcolor=(1, 0, 0))\n",
"b8 = circle((-x, y), r, rgbcolor=(1, 0, 0))\n",
"\n",
"show(c1 + c2 + c3 + c4 + c5 + c6 + b1 + b2 + b3 + b4 + b5 + b6 + b7 + b8)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The quadratic form for the cubical packing\n",
"\n",
"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": 21,
"metadata": {},
"outputs": [],
"source": [
"def abbc_coords(b, h1, h2):\n",
" 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": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[ 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": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Wc = matrix([\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.transpose().rref()"
]
},
{
"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": 23,
"metadata": {},
"outputs": [],
"source": [
"W = matrix([\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": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[ 1 -3 -3 -3]\n",
"[-3 1 -3 -3]\n",
"[-3 -3 1 -3]\n",
"[-3 -3 -3 1]"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = W * P * W.transpose()\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[ 2 -6/5 -6/5 -6/5]\n",
"[-6/5 2 -6/5 -6/5]\n",
"[-6/5 -6/5 2 -6/5]\n",
"[-6/5 -6/5 -6/5 2]"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"32 * m.inverse() * 2/5"
]
},
{
"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/5 6/5 6/5] [ 1 0 0 0] [ 1 0 0 0]\n",
"[ 0 1 0 0] [6/5 -1 6/5 6/5] [ 0 1 0 0]\n",
"[ 0 0 1 0] [ 0 0 1 0] [6/5 6/5 -1 6/5]\n",
"[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1],\n",
"\n",
"[ 1 0 0 0]\n",
"[ 0 1 0 0]\n",
"[ 0 0 1 0]\n",
"[6/5 6/5 6/5 -1]\n",
"]"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S_i = weyl_generators(m.inverse(), standard_basis(4))\n",
"S_i"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"[ 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]\n",
"[ 0 0 0 0 0 0 0 0]\n",
"[ 0 0 0 0 0 0 0 0]\n",
"[ 0 0 0 0 0 0 0 0]\n",
"[ 0 0 0 0 0 0 0 0]"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(Wc * P * Wc.transpose()).rref()"
]
},
{
"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",
" n = root_matrix.dimensions()[1]\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",
" relations_temp = vector([ var('b' + str(i)) for i in range(1, n + 1)]) * root_matrix.transpose().rref()\n",
" relations = []\n",
" for i, expr in enumerate(relations_temp):\n",
" relations.append(var('b' + str(i + 1)) == expr)\n",
" \n",
" # step 2: find matrix of quadratic form\n",
" W = root_matrix[-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 relations[4:], M.inverse(), D[1][1]"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[2 \\, b_{5} = b_{1} + b_{2} + b_{3} - b_{4}, 2 \\, b_{6} = b_{1} + b_{2} - b_{3} + b_{4}, 2 \\, b_{7} = b_{1} - b_{2} + b_{3} + b_{4}, 2 \\, b_{8} = -b_{1} + b_{2} + b_{3} + b_{4}\\right]\n",
"\\end{math}"
],
"text/plain": [
"[2*b5 == b1 + b2 + b3 - b4,\n",
" 2*b6 == b1 + b2 - b3 + b4,\n",
" 2*b7 == b1 - b2 + b3 + b4,\n",
" 2*b8 == -b1 + b2 + b3 + b4]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"5 & -3 & -3 & -3 \\\\\n",
"-3 & 5 & -3 & -3 \\\\\n",
"-3 & -3 & 5 & -3 \\\\\n",
"-3 & -3 & -3 & 5\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ 5 -3 -3 -3]\n",
"[-3 5 -3 -3]\n",
"[-3 -3 5 -3]\n",
"[-3 -3 -3 5]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}5 \\, b_{1}^{2} - 6 \\, b_{1} b_{2} + 5 \\, b_{2}^{2} - 6 \\, b_{1} b_{3} - 6 \\, b_{2} b_{3} + 5 \\, b_{3}^{2} - 6 \\, b_{1} b_{4} - 6 \\, b_{2} b_{4} - 6 \\, b_{3} b_{4} + 5 \\, b_{4}^{2}\n",
"\\end{math}"
],
"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"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(\\begin{array}{rrrr}\n",
"-1 & \\frac{6}{5} & \\frac{6}{5} & \\frac{6}{5} \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"\\frac{6}{5} & -1 & \\frac{6}{5} & \\frac{6}{5} \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"\\frac{6}{5} & \\frac{6}{5} & -1 & \\frac{6}{5} \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"\\frac{6}{5} & \\frac{6}{5} & \\frac{6}{5} & -1\n",
"\\end{array}\\right)\\right]\n",
"\\end{math}"
],
"text/plain": [
"[\n",
"[ -1 6/5 6/5 6/5] [ 1 0 0 0] [ 1 0 0 0]\n",
"[ 0 1 0 0] [6/5 -1 6/5 6/5] [ 0 1 0 0]\n",
"[ 0 0 1 0] [ 0 0 1 0] [6/5 6/5 -1 6/5]\n",
"[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1],\n",
"\n",
"[ 1 0 0 0]\n",
"[ 0 1 0 0]\n",
"[ 0 0 1 0]\n",
"[6/5 6/5 6/5 -1]\n",
"]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# cubical\n",
"relation, mat, equation = quadform_from_root(Wc)\n",
"show([2 * eq for eq in relation])\n",
"show(32 * mat)\n",
"show(32 * equation)\n",
"show(weyl_generators(32 * mat, standard_basis(4)))"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[b_{5} = b_{1} - b_{2} + b_{4}, b_{6} = b_{1} - b_{3} + b_{4}\\right]\n",
"\\end{math}"
],
"text/plain": [
"[b5 == b1 - b2 + b4, b6 == b1 - b3 + b4]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"1 & -2 & -2 & -1 \\\\\n",
"-2 & 4 & 0 & -2 \\\\\n",
"-2 & 0 & 4 & -2 \\\\\n",
"-1 & -2 & -2 & 1\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ 1 -2 -2 -1]\n",
"[-2 4 0 -2]\n",
"[-2 0 4 -2]\n",
"[-1 -2 -2 1]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}b_{1}^{2} - 4 \\, b_{1} b_{2} + 4 \\, b_{2}^{2} - 4 \\, b_{1} b_{3} + 4 \\, b_{3}^{2} - 2 \\, b_{1} b_{4} - 4 \\, b_{2} b_{4} - 4 \\, b_{3} b_{4} + b_{4}^{2}\n",
"\\end{math}"
],
"text/plain": [
"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"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(\\begin{array}{rrrr}\n",
"-1 & 4 & 4 & 2 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"1 & -1 & 0 & 1 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"1 & 0 & -1 & 1 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"2 & 4 & 4 & -1\n",
"\\end{array}\\right)\\right]\n",
"\\end{math}"
],
"text/plain": [
"[\n",
"[-1 4 4 2] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n",
"[ 0 1 0 0] [ 1 -1 0 1] [ 0 1 0 0] [ 0 1 0 0]\n",
"[ 0 0 1 0] [ 0 0 1 0] [ 1 0 -1 1] [ 0 0 1 0]\n",
"[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 2 4 4 -1]\n",
"]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"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",
"show(relation)\n",
"show(8 * mat)\n",
"show(8 * equation)\n",
"show(weyl_generators(8 * mat, standard_basis(4)))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\right]\n",
"\\end{math}"
],
"text/plain": [
"[]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"1 & -1 & -1 & -1 \\\\\n",
"-1 & 1 & -1 & -1 \\\\\n",
"-1 & -1 & 1 & -1 \\\\\n",
"-1 & -1 & -1 & 1\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ 1 -1 -1 -1]\n",
"[-1 1 -1 -1]\n",
"[-1 -1 1 -1]\n",
"[-1 -1 -1 1]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}b_{1}^{2} - 2 \\, b_{1} b_{2} + b_{2}^{2} - 2 \\, b_{1} b_{3} - 2 \\, b_{2} b_{3} + b_{3}^{2} - 2 \\, b_{1} b_{4} - 2 \\, b_{2} b_{4} - 2 \\, b_{3} b_{4} + b_{4}^{2}\n",
"\\end{math}"
],
"text/plain": [
"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"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(\\begin{array}{rrrr}\n",
"-1 & 2 & 2 & 2 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"2 & -1 & 2 & 2 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"2 & 2 & -1 & 2 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"2 & 2 & 2 & -1\n",
"\\end{array}\\right)\\right]\n",
"\\end{math}"
],
"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",
"]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"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",
"show(relation)\n",
"show(4 * mat)\n",
"show(4 * equation)\n",
"show(weyl_generators(4 * mat, standard_basis(4)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# $n$-Gon Base Pyramid\n",
"\n",
"The goal here is to find the quadratic form for an arbitrary $n$-gon base pyramid. The key to the whole process is a magic formula Dylan found for the bilinear form between two circles in an $n$-gon base pyramidal packing, namely $$\n",
" \\frac{1 - \\cos\\left(\\frac{2\\pi}{n}\\right) - 4\\sin^2\\left(\\frac{p\\pi}{n}\\right)}{1-\\cos\\left(\\frac{2\\pi}{n}\\right)}\n",
"$$\n",
"\n",
"where $p$ is how many circles are between the two circles in question. If we are finding the bilinear form between the central circle and any other circle it will always be $-1$ so we don't need to worry about that case."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"# function to compute W^T*P*W for n-gon pyramidal packing\n",
"def wmatrix(n):\n",
" vals = []\n",
" for i in range(n+1):\n",
" row = []\n",
" for j in range(n+1):\n",
" if i == j: # same vertex bilinear form'd with itself, so 1\n",
" row.append(1)\n",
" elif i ==0 or j == 0: # vertex bilinear form'd with special point, so tangent and therefore -1\n",
" row.append(-1)\n",
" else:\n",
" p = abs(i - j) # otherwise Dylan's crazy formula\n",
" row.append(\n",
" (1 - cos(2 * pi / n) - 4 * sin(pi * p / n)^2) / (1 - cos(2 * pi / n))\n",
" )\n",
" vals.append(row)\n",
" return matrix(vals)"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"# just convenience function for quadratic forms\n",
"def qform(matrix, vector):\n",
" return vector * matrix * vector"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"# function to compute the linear relations and quadratic formula for an ngon base pyramid\n",
"# the circles will be numbered 1 in the center, then 2 through n winding around the circle by step, so\n",
"# for n = 4 with step = 1, we have\n",
"# b3\n",
"# b4 b1 b2\n",
"# b5\n",
"# and the quadratic form is in terms of b1, b2, b3, and b4\n",
"# and for n = 6 with step = 2, we have\n",
"# b4 b3\n",
"# b5 b1 b2\n",
"# b6 b7\n",
"# and the quadratic form is in terms of b1, b3, b5, and b7\n",
"def ngon_linear_relations_and_quadratic_form(n, step=1):\n",
" mat = wmatrix(n)\n",
" \n",
" # work out initial relations\n",
" relations_temp = vector([ var('b' + str(i)) for i in range(1, n + 2) ]) * mat.transpose().rref()\n",
" relations = []\n",
" for i in range(4, n + 1):\n",
" relations.append(var('b' + str(i + 1)) == relations_temp[i])\n",
" \n",
" # rewrite the relations in terms of the variables we care about, depends on the step\n",
" shuffled = (list(range(1, n + 2)) * step)[::step]\n",
" targets = [ var('b' + str(i)) for i in shuffled[4:] ]\n",
" relations = solve(relations, *targets)[0]\n",
" \n",
" # find the matrix corresponding to the quadratic form, picking the appropriate rows from the matrix\n",
" Q = mat[:4*step:step,:4*step:step].inverse()\n",
" \n",
" # find the quadratic form in variables; proper units will satisfy this being equal to zero\n",
" nqform = qform(Q, vector([ var('b' + str(i)) for i in range(1, 5) ]))\n",
" \n",
" return relations, Q, expand(nqform)"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}b_{5} = b_{2} - b_{3} + b_{4}\n",
"\\end{math}"
],
"text/plain": [
"b5 == b2 - b3 + b4"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"4 & -2 & 0 & -2 \\\\\n",
"-2 & 1 & -2 & -1 \\\\\n",
"0 & -2 & 4 & -2 \\\\\n",
"-2 & -1 & -2 & 1\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ 4 -2 0 -2]\n",
"[-2 1 -2 -1]\n",
"[ 0 -2 4 -2]\n",
"[-2 -1 -2 1]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}4 \\, b_{1}^{2} - 4 \\, b_{1} b_{2} + b_{2}^{2} - 4 \\, b_{2} b_{3} + 4 \\, b_{3}^{2} - 4 \\, b_{1} b_{4} - 2 \\, b_{2} b_{4} - 4 \\, b_{3} b_{4} + b_{4}^{2}\n",
"\\end{math}"
],
"text/plain": [
"4*b1^2 - 4*b1*b2 + b2^2 - 4*b2*b3 + 4*b3^2 - 4*b1*b4 - 2*b2*b4 - 4*b3*b4 + b4^2"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(\\begin{array}{rrrr}\n",
"-1 & 1 & 0 & 1 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"4 & -1 & 4 & 2 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 1 & -1 & 1 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"4 & 2 & 4 & -1\n",
"\\end{array}\\right)\\right]\n",
"\\end{math}"
],
"text/plain": [
"[\n",
"[-1 1 0 1] [ 1 0 0 0] [ 1 0 0 0] [ 1 0 0 0]\n",
"[ 0 1 0 0] [ 4 -1 4 2] [ 0 1 0 0] [ 0 1 0 0]\n",
"[ 0 0 1 0] [ 0 0 1 0] [ 0 1 -1 1] [ 0 0 1 0]\n",
"[ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 1], [ 4 2 4 -1]\n",
"]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"relations, Q, nqform = ngon_linear_relations_and_quadratic_form(4)\n",
"\n",
"show(relations)\n",
"show(8 * Q)\n",
"show(8 * nqform)\n",
"show(weyl_generators(8 * Q, standard_basis(4)))"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[3 \\, b_{2} = 2 \\, b_{3} - b_{5} + 2 \\, b_{7}, 3 \\, b_{4} = 2 \\, b_{3} + 2 \\, b_{5} - b_{7}, 3 \\, b_{6} = -b_{3} + 2 \\, b_{5} + 2 \\, b_{7}\\right]\n",
"\\end{math}"
],
"text/plain": [
"[3*b2 == 2*b3 - b5 + 2*b7, 3*b4 == 2*b3 + 2*b5 - b7, 3*b6 == -b3 + 2*b5 + 2*b7]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"9 & -1 & -1 & -1 \\\\\n",
"-1 & 1 & -1 & -1 \\\\\n",
"-1 & -1 & 1 & -1 \\\\\n",
"-1 & -1 & -1 & 1\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ 9 -1 -1 -1]\n",
"[-1 1 -1 -1]\n",
"[-1 -1 1 -1]\n",
"[-1 -1 -1 1]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}9 \\, b_{1}^{2} - 2 \\, b_{1} b_{2} + b_{2}^{2} - 2 \\, b_{1} b_{3} - 2 \\, b_{2} b_{3} + b_{3}^{2} - 2 \\, b_{1} b_{4} - 2 \\, b_{2} b_{4} - 2 \\, b_{3} b_{4} + b_{4}^{2}\n",
"\\end{math}"
],
"text/plain": [
"9*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"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(\\begin{array}{rrrr}\n",
"-1 & \\frac{2}{9} & \\frac{2}{9} & \\frac{2}{9} \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"2 & -1 & 2 & 2 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"2 & 2 & -1 & 2 \\\\\n",
"0 & 0 & 0 & 1\n",
"\\end{array}\\right), \\left(\\begin{array}{rrrr}\n",
"1 & 0 & 0 & 0 \\\\\n",
"0 & 1 & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"2 & 2 & 2 & -1\n",
"\\end{array}\\right)\\right]\n",
"\\end{math}"
],
"text/plain": [
"[\n",
"[ -1 2/9 2/9 2/9] [ 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",
"]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"relations, Q, nqform = ngon_linear_relations_and_quadratic_form(6, 2)\n",
"\n",
"show([3 * eq for eq in relations])\n",
"show(12 * Q)\n",
"show(12 * nqform)\n",
"show(weyl_generators(12 * Q, standard_basis(4)))"
]
}
],
"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"
},
"name": "Apollonian Circle Packings.ipynb"
},
"nbformat": 4,
"nbformat_minor": 4
}