Circle-Packings/Apollonian Circle Packings.ipynb

{
 "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": 3,
   "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": 4,
   "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": 4,
     "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": 5,
   "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": 6,
   "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": 6,
     "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": 7,
   "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": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.inverse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = W2.transpose() * M.inverse() * W2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.inverse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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": 11,
     "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": 12,
   "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": 13,
   "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": 14,
   "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": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "4 * M.inverse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_i = weyl_generators(4 * M.inverse(), standard_basis(4))\n",
    "S_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7, 2, 4, 6)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "root = vector([-1, 2, 4, 6])\n",
    "S1 * root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1, 4, 4, 6)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S2 * root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1, 2, 2, 6)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S3 * root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1, 2, 4, 14)"
      ]
     },
     "execution_count": 20,
     "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": 21,
   "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": 22,
   "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/
      "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": 23,
   "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": 24,
   "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": 24,
     "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": 25,
   "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": 26,
   "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": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = W * P * W.transpose()\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "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": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "32 * m.inverse() * 2/5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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": 28,
     "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": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = W2.transpose() * m.inverse() * W2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "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": 30,
     "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": 31,
   "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": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_i = weyl_generators(m.inverse(), standard_basis(4))\n",
    "S_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "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": 32,
     "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": 33,
   "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[-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": 34,
   "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",
       "[ -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": 34,
     "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": 35,
   "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]  [ 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",
       "]"
      ]
     },
     "execution_count": 35,
     "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": 36,
   "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": 36,
     "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 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. It feels very strange, since it's only really two dimensional, rather than four, and it doesn't quite look like the equations hold up. But it does absolutely work for the tetrahedral packing, which is awesome."
   ]
  }
 ],
 "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
}