Circle-Packings/Apollonian Circle Packings.ipynb

371 lines
39 KiB
Text
Raw Normal View History

2021-03-05 18:47:48 -08:00
{
"cells": [
{
"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": "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",
" [3, 1, sqrt(2), 0]\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"[ 1 -1 -1 -1]\n",
"[-1 1 -1 -1]\n",
"[-1 -1 1 -1]\n",
"[-1 -1 -1 -1]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = W * P * W.transpose()\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"-8"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m.determinant()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"[ 1/2 0 0 -1/2]\n",
"[ 0 1/2 0 -1/2]\n",
"[ 0 0 1/2 -1/2]\n",
"[-1/2 -1/2 -1/2 1/2]"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m.inverse()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"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",
"bt4 = var('bt4')\n",
"b4 = var('b4')\n",
"h41 = var('h41')\n",
"h42 = var('h42')\n",
"\n",
"\n",
"W2 = matrix([\n",
" [bt1, b1, h11, h12],\n",
" [bt2, b2, h21, h22],\n",
" [bt3, b3, h31, h32],\n",
" [bt4, b4, h41, h42],\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"D = W2.transpose() * m.inverse() * W2"
]
},
{
"cell_type": "code",
"execution_count": 11,
"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": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P.inverse()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[b4 == b1 + b2 + b3 - sqrt(2*b1*b2 + 2*(b1 + b2)*b3), b4 == b1 + b2 + b3 + sqrt(2*b1*b2 + 2*(b1 + b2)*b3)]"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"quad = 2 * factor(simplify(D[1][1]))\n",
"solve(quad == 0, b4)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"b1^2 + b2^2 + b3^2 - 2*b1*b4 - 2*b2*b4 - 2*b3*b4 + b4^2"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"quad"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1/2*b1*h11 - 1/2*b4*h11 + 1/2*b2*h21 - 1/2*b4*h21 + 1/2*b3*h31 - 1/2*b4*h31 - 1/2*b1*h41 - 1/2*b2*h41 - 1/2*b3*h41 + 1/2*b4*h41"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"factor(simplify(D[2][1]))"
]
},
{
"cell_type": "code",
"execution_count": 20,
"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": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1/2, 1/2, 1/2*sqrt(2))"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"circle_from_points((0, 1), (1, 0), (1, 1))"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"sqrt(2)\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(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(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)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "SageMath 9.2",
"language": "sage",
"name": "sagemath"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}