platonic-solids-and-prime-numbers-10-02-02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 『正多面体と素数』10.2.2の計算ノート"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(x, y, X)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var('x y X')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}X^{11} + 11 \\, X^{6} - X</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}X^{11} + 11 \\, X^{6} - X\n",
       "\\end{math}"
      ],
      "text/plain": [
       "X^11 + 11*X^6 - X"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "F_12(x,y)=x^11*y+11*x^6*y^6-x*y^11\n",
    "show(F_12(X,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ヘッシアンの行列式を計算\n",
    "def det_hessian(f):\n",
    "    return det(f.hessian()).expand().simplify_full()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_20_0(x,y)=det_hessian(F_12(x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}-121 \\, x^{20} + 27588 \\, x^{15} y^{5} - 59774 \\, x^{10} y^{10} - 27588 \\, x^{5} y^{15} - 121 \\, y^{20}</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}-121 \\, x^{20} + 27588 \\, x^{15} y^{5} - 59774 \\, x^{10} y^{10} - 27588 \\, x^{5} y^{15} - 121 \\, y^{20}\n",
       "\\end{math}"
      ],
      "text/plain": [
       "-121*x^20 + 27588*x^15*y^5 - 59774*x^10*y^10 - 27588*x^5*y^15 - 121*y^20"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(F_20_0(x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x^{20} - 228 \\, x^{15} y^{5} + 494 \\, x^{10} y^{10} + 228 \\, x^{5} y^{15} + y^{20}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#最高次の係数を1へ\n",
    "F_20(x,y) = (F_20_0(x,y)/F_20_0(x,y).coefficient(x^20)).expand().simplify_full()\n",
    "latex(F_20(x,y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_{20}(x,y)=x^{20} - 228 \\, x^{15} y^{5} + 494 \\, x^{10} y^{10} + 228 \\, x^{5} y^{15} + y^{20}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "X^{20} - 228 \\, X^{15} + 494 \\, X^{10} + 228 \\, X^{5} + 1"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latex(F_20(X,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "『正多面体と素数』p.128\n",
    "\n",
    "$$\n",
    "H(X,1)=F_R^{\\{2\\}}(X)=X^{20} - 228 \\, X^{15} + 494 \\, X^{10} + 228 \\, X^{5} + 1\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#二つの関数からヤコビアンを計算\n",
    "def det_jacobian(f,g):\n",
    "    return det(jacobian([f,g],[x,y])).expand().simplify_full()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_30_0(x,y)=det_jacobian(F_12(x,y), F_20(x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}-20 \\, x^{30} - 10440 \\, x^{25} y^{5} + 200100 \\, x^{20} y^{10} + 200100 \\, x^{10} y^{20} + 10440 \\, x^{5} y^{25} - 20 \\, y^{30}</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}-20 \\, x^{30} - 10440 \\, x^{25} y^{5} + 200100 \\, x^{20} y^{10} + 200100 \\, x^{10} y^{20} + 10440 \\, x^{5} y^{25} - 20 \\, y^{30}\n",
       "\\end{math}"
      ],
      "text/plain": [
       "-20*x^30 - 10440*x^25*y^5 + 200100*x^20*y^10 + 200100*x^10*y^20 + 10440*x^5*y^25 - 20*y^30"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(F_30_0(x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x^{30} + 522 \\, x^{25} y^{5} - 10005 \\, x^{20} y^{10} - 10005 \\, x^{10} y^{20} - 522 \\, x^{5} y^{25} + y^{30}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#最高次の係数を1へ\n",
    "F_30(x,y) = (F_30_0(x,y)/F_30_0(x,y).coefficient(x^30)).expand().simplify_full()\n",
    "latex(F_30(x,y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_{30}(x,y)=x^{30} + 522 \\, x^{25} y^{5} - 10005 \\, x^{20} y^{10} - 10005 \\, x^{10} y^{20} - 522 \\, x^{5} y^{25} + y^{30}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "X^{30} + 522 \\, X^{25} - 10005 \\, X^{20} - 10005 \\, X^{10} - 522 \\, X^{5} + 1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latex(F_30(X,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "『正多面体と素数』p.128\n",
    "\n",
    "$$\n",
    "T(X,1)=F_R^{\\{1\\}}(X)=X^{30} + 522 \\, X^{25} - 10005 \\, X^{20} - 10005 \\, X^{10} - 522 \\, X^{5} + 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-1728, 0]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#多項式として除算\n",
    "(F_20(X,1)^3-F_30(X,1)^2).maxima_methods().divide((F_12(X,1)^5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "すなわち、\n",
    "$$\n",
    "1728F_{12}^5 + F_{20}^3 - F_{30}^2=0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#有理数の分子と分母を素因数分解してlatex形式で返す関数\n",
    "def get_factorized_rational(numerator, denominator):\n",
    "    factors_a = list(factor(numerator))\n",
    "    factors_b = list(factor(denominator))\n",
    "    num = \"\"\n",
    "    if(len(factors_a)==0):\n",
    "        num = \"1\"\n",
    "    else:\n",
    "        for i, t in enumerate(factors_a):\n",
    "            p, k = t\n",
    "            if(k==1):\n",
    "                num += f\"{p}\"\n",
    "            else:\n",
    "                num += f\"{p}^{{{k}}}\"\n",
    "            if(i<len(factors_a)-1):\n",
    "                num += r\" \\!\\cdot\\! \"\n",
    "                \n",
    "    den = \"\"\n",
    "    if(len(factors_b)==0):\n",
    "        den = \"1\"\n",
    "    else:\n",
    "        for i, t in enumerate(factors_b):\n",
    "            p, k = t\n",
    "            if(k==1):\n",
    "                den += f\"{p}\"\n",
    "            else:\n",
    "                den += f\"{p}^{{{k}}}\"\n",
    "            if(i<len(factors_b)-1):\n",
    "                den += r\" \\!\\cdot\\! \"\n",
    "    return \"\\\\frac{{{0}}}{{{1}}}\".format(num, den)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 整数の平方根をいい感じにlatex形式で返す関数\n",
    "def format_factors(squared_factors):\n",
    "    if(len(squared_factors)==0):\n",
    "        return '1'\n",
    "    integer_parts = []\n",
    "    sqrt_parts = []\n",
    "    for i, t in enumerate(squared_factors):\n",
    "        p,k = t\n",
    "        if(k%2==0):\n",
    "            if(k//2==1):\n",
    "                integer_parts.append(f'{p}')\n",
    "            else:\n",
    "                integer_parts.append(f'{p}^{k//2}')\n",
    "        else:\n",
    "            if(k//2==1):\n",
    "                integer_parts.append(f'{p}')\n",
    "            elif(k//2>1):\n",
    "                integer_parts.append(f'{p}^{k//2}')\n",
    "            sqrt_parts.append(f'{p}')\n",
    "    integer_part = ' \\\\!\\\\cdot\\\\! '.join(integer_parts)\n",
    "    if(len(sqrt_parts)==0):\n",
    "        sqrt_part = ''\n",
    "    else:\n",
    "        sqrt_part_inner = ' \\\\!\\\\cdot\\\\! '.join(sqrt_parts)\n",
    "        sqrt_part = f'\\\\sqrt{{{sqrt_part_inner}}}'\n",
    "    return f'{integer_part}{sqrt_part}'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#x,yの同次（deg次）多項式fをいい感じに表示する\n",
    "def show_polynomial(f, deg):\n",
    "    for i in reversed(range(deg+1)):\n",
    "        j = deg-i\n",
    "        if(i<j):break\n",
    "        c=f.coefficient(x^i*y^j)\n",
    "        p = numerator(c*c).pyobject()\n",
    "        q = denominator(c*c).pyobject()\n",
    "        if(p>0):\n",
    "            sgn1 = '+' if f.coefficient(x^i*y^j)/f.coefficient(x^j*y^i) == 1 else '-'\n",
    "            if sign(c)>0:\n",
    "                sgn2=\"+\"\n",
    "            else:\n",
    "                sgn2=\"-\"\n",
    "            a = format_factors(p.factor())\n",
    "            b = format_factors(q.factor())\n",
    "            print( sgn2 + (f\"{a}\" if b==\"1\" else \"\\\\frac{{{0}}}{{{1}}} \".format(a,b)) \\\n",
    "                    + (\"\\\\left(x^{{{0}}}y^{{{1}}} {2} x^{{{1}}}y^{{{0}}}\\\\right)\".format(i,j,sgn1) if i != j\n",
    "                          else \"x^{{{0}}}y^{{{1}}}\".format(i, deg-i)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+1\\left(x^{11}y^{1} - x^{1}y^{11}\\right)\n",
      "+11x^{6}y^{6}\n"
     ]
    }
   ],
   "source": [
    "show_polynomial(F_12(x,y), 12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+1\\left(x^{20}y^{0} + x^{0}y^{20}\\right)\n",
      "-2^2 \\!\\cdot\\! 3 \\!\\cdot\\! 19\\left(x^{15}y^{5} - x^{5}y^{15}\\right)\n",
      "+2 \\!\\cdot\\! 13 \\!\\cdot\\! 19x^{10}y^{10}\n"
     ]
    }
   ],
   "source": [
    "show_polynomial(F_20(x,y),20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+1\\left(x^{30}y^{0} + x^{0}y^{30}\\right)\n",
      "+2 \\!\\cdot\\! 3^2 \\!\\cdot\\! 29\\left(x^{25}y^{5} - x^{5}y^{25}\\right)\n",
      "-3 \\!\\cdot\\! 5 \\!\\cdot\\! 23 \\!\\cdot\\! 29\\left(x^{20}y^{10} + x^{10}y^{20}\\right)\n"
     ]
    }
   ],
   "source": [
    "show_polynomial(F_30(x,y),30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_32(x,y)=(F_12(x,y)*F_20(x,y)).expand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+1\\left(x^{31}y^{1} - x^{1}y^{31}\\right)\n",
      "-7 \\!\\cdot\\! 31\\left(x^{26}y^{6} + x^{6}y^{26}\\right)\n",
      "-5 \\!\\cdot\\! 13 \\!\\cdot\\! 31\\left(x^{21}y^{11} - x^{11}y^{21}\\right)\n",
      "+2 \\!\\cdot\\! 5 \\!\\cdot\\! 19 \\!\\cdot\\! 31x^{16}y^{16}\n"
     ]
    }
   ],
   "source": [
    "show_polynomial(F_32(x,y),32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_42(x,y)=(F_12(x,y)*F_30(x,y)).expand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+1\\left(x^{41}y^{1} - x^{1}y^{41}\\right)\n",
      "+13 \\!\\cdot\\! 41\\left(x^{36}y^{6} + x^{6}y^{36}\\right)\n",
      "-2^3 \\!\\cdot\\! 13 \\!\\cdot\\! 41\\left(x^{31}y^{11} - x^{11}y^{31}\\right)\n",
      "-3 \\!\\cdot\\! 29 \\!\\cdot\\! 31 \\!\\cdot\\! 41\\left(x^{26}y^{16} + x^{16}y^{26}\\right)\n"
     ]
    }
   ],
   "source": [
    "show_polynomial(F_42(x,y),42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_50(x,y)=(F_20(x,y)*F_30(x,y)).expand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+1\\left(x^{50}y^{0} + x^{0}y^{50}\\right)\n",
      "+2 \\!\\cdot\\! 3 \\!\\cdot\\! 7^2\\left(x^{45}y^{5} - x^{5}y^{45}\\right)\n",
      "-7^2 \\!\\cdot\\! 43 \\!\\cdot\\! 61\\left(x^{40}y^{10} + x^{10}y^{40}\\right)\n",
      "+2^2 \\!\\cdot\\! 3 \\!\\cdot\\! 7 \\!\\cdot\\! 19 \\!\\cdot\\! 37 \\!\\cdot\\! 43\\left(x^{35}y^{15} - x^{15}y^{35}\\right)\n",
      "-2 \\!\\cdot\\! 7^2 \\!\\cdot\\! 31 \\!\\cdot\\! 37 \\!\\cdot\\! 43\\left(x^{30}y^{20} + x^{20}y^{30}\\right)\n"
     ]
    }
   ],
   "source": [
    "show_polynomial(F_50(x,y),50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_62(x,y)=(F_12(x,y)*F_20(x,y)*F_30(x,y)).expand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+1\\left(x^{61}y^{1} - x^{1}y^{61}\\right)\n",
      "+5 \\!\\cdot\\! 61\\left(x^{56}y^{6} + x^{6}y^{56}\\right)\n",
      "-2 \\!\\cdot\\! 13 \\!\\cdot\\! 61 \\!\\cdot\\! 79\\left(x^{51}y^{11} - x^{11}y^{51}\\right)\n",
      "+5 \\!\\cdot\\! 7 \\!\\cdot\\! 17 \\!\\cdot\\! 31 \\!\\cdot\\! 61\\left(x^{46}y^{16} + x^{16}y^{46}\\right)\n",
      "+5 \\!\\cdot\\! 7 \\!\\cdot\\! 11 \\!\\cdot\\! 23 \\!\\cdot\\! 43 \\!\\cdot\\! 61\\left(x^{41}y^{21} - x^{21}y^{41}\\right)\n",
      "-2 \\!\\cdot\\! 7 \\!\\cdot\\! 37 \\!\\cdot\\! 41 \\!\\cdot\\! 43 \\!\\cdot\\! 61\\left(x^{36}y^{26} + x^{26}y^{36}\\right)\n"
     ]
    }
   ],
   "source": [
    "show_polynomial(F_62(x,y),62)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 同時多項式の係数の正規化\n",
    "def get_normalized_coefficients(f, deg):\n",
    "    v = []\n",
    "    for i in reversed(range(deg+1)):\n",
    "        j = deg-i\n",
    "        c = f.coefficient(x^i*y^j) * sqrt(factorial(i)*factorial(j))\n",
    "        v.append(c)\n",
    "    return vector(v).normalized()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(0,\\,\\frac{1}{55} \\, \\sqrt{77} \\sqrt{11},\\,0,\\,0,\\,0,\\,0,\\,\\frac{1}{5} \\, \\sqrt{11},\\,0,\\,0,\\,0,\\,0,\\,-\\frac{1}{55} \\, \\sqrt{77} \\sqrt{11},\\,0\\right)</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(0,\\,\\frac{1}{55} \\, \\sqrt{77} \\sqrt{11},\\,0,\\,0,\\,0,\\,0,\\,\\frac{1}{5} \\, \\sqrt{11},\\,0,\\,0,\\,0,\\,0,\\,-\\frac{1}{55} \\, \\sqrt{77} \\sqrt{11},\\,0\\right)\n",
       "\\end{math}"
      ],
      "text/plain": [
       "(0, 1/55*sqrt(77)*sqrt(11), 0, 0, 0, 0, 1/5*sqrt(11), 0, 0, 0, 0, -1/55*sqrt(77)*sqrt(11), 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "u = get_normalized_coefficients(F_12(x,y), 12)\n",
    "show(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 正規化された係数の列を球面調和関数の重ね合わせとして表示\n",
    "def show_as_spherical_harmonics_superposition(v, deg):\n",
    "    out = \"\"\n",
    "    for i in range(deg+1):\n",
    "        j = deg-i\n",
    "        if(i>j):break\n",
    "        c = v[i]\n",
    "        if(c!=0):\n",
    "            sgn1 = \"+\" if sign(c)>0 else \"-\"\n",
    "            sgn2 = \"+\" if v[i]==v[j] else \"-\"\n",
    "            p = numerator(c*c).pyobject()\n",
    "            q = denominator(c*c).pyobject()\n",
    "            r = get_factorized_rational(p,q)\n",
    "            if(i==j):\n",
    "                out += sgn1 + \"\\\\sqrt{{{0}}}\".format(r) + \"Y^{{{0}}}_{{{1}}}\".format(deg/2, -(i-j)/2)\n",
    "            else:\n",
    "                out += (sgn1 + \"\\\\sqrt{{{0}}}\".format(r)\n",
    "                      + \"\\\\left(Y^{{{0}}}_{{{1}}} {2} Y^{{{0}}}_{{{3}}}\\\\right)\".format(deg/2,-(i-j)/2 , sgn2, (i-j)/2))\n",
    "            out += \"\\n\"\n",
    "    print(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+\\sqrt{\\frac{7}{5^{2}}}\\left(Y^{6}_{5} - Y^{6}_{-5}\\right)\n",
      "+\\sqrt{\\frac{11}{5^{2}}}Y^{6}_{0}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "show_as_spherical_harmonics_superposition(u, 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\mathcal{F}_{12}(\\theta,\\varphi)\n",
    "=\n",
    "+\\sqrt{\\frac{7}{5^{2}}}\\left(Y^{6}_{-5} - Y^{6}_{5}\\right)\n",
    "+\\sqrt{\\frac{11}{5^{2}}}Y^{6}_{0}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\frac{1}{18525} \\, \\sqrt{46189} \\sqrt{741},\\,0,\\,0,\\,0,\\,0,\\,-\\frac{1}{975} \\, \\sqrt{741} \\sqrt{429},\\,0,\\,0,\\,0,\\,0,\\,\\frac{1}{75} \\, \\sqrt{741},\\,0,\\,0,\\,0,\\,0,\\,\\frac{1}{975} \\, \\sqrt{741} \\sqrt{429},\\,0,\\,0,\\,0,\\,0,\\,\\frac{1}{18525} \\, \\sqrt{46189} \\sqrt{741}\\right)</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\frac{1}{18525} \\, \\sqrt{46189} \\sqrt{741},\\,0,\\,0,\\,0,\\,0,\\,-\\frac{1}{975} \\, \\sqrt{741} \\sqrt{429},\\,0,\\,0,\\,0,\\,0,\\,\\frac{1}{75} \\, \\sqrt{741},\\,0,\\,0,\\,0,\\,0,\\,\\frac{1}{975} \\, \\sqrt{741} \\sqrt{429},\\,0,\\,0,\\,0,\\,0,\\,\\frac{1}{18525} \\, \\sqrt{46189} \\sqrt{741}\\right)\n",
       "\\end{math}"
      ],
      "text/plain": [
       "(1/18525*sqrt(46189)*sqrt(741), 0, 0, 0, 0, -1/975*sqrt(741)*sqrt(429), 0, 0, 0, 0, 1/75*sqrt(741), 0, 0, 0, 0, 1/975*sqrt(741)*sqrt(429), 0, 0, 0, 0, 1/18525*sqrt(46189)*sqrt(741))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "v = get_normalized_coefficients(F_20(x,y), 20)\n",
    "show(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+\\sqrt{\\frac{11 \\!\\cdot\\! 17}{3 \\!\\cdot\\! 5^{4}}}\\left(Y^{10}_{10} + Y^{10}_{-10}\\right)\n",
      "-\\sqrt{\\frac{11 \\!\\cdot\\! 19}{5^{4}}}\\left(Y^{10}_{5} - Y^{10}_{-5}\\right)\n",
      "+\\sqrt{\\frac{13 \\!\\cdot\\! 19}{3 \\!\\cdot\\! 5^{4}}}Y^{10}_{0}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "show_as_spherical_harmonics_superposition(v, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "w=get_normalized_coefficients(F_30(x,y),30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+\\sqrt{\\frac{7 \\!\\cdot\\! 11 \\!\\cdot\\! 13}{2^{2} \\!\\cdot\\! 5^{5}}}\\left(Y^{15}_{15} + Y^{15}_{-15}\\right)\n",
      "+\\sqrt{\\frac{3 \\!\\cdot\\! 11 \\!\\cdot\\! 29}{2 \\!\\cdot\\! 5^{5}}}\\left(Y^{15}_{10} - Y^{15}_{-10}\\right)\n",
      "-\\sqrt{\\frac{23 \\!\\cdot\\! 29}{2^{2} \\!\\cdot\\! 5^{4}}}\\left(Y^{15}_{5} + Y^{15}_{-5}\\right)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "show_as_spherical_harmonics_superposition(w, 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+\\sqrt{\\frac{3 \\!\\cdot\\! 17 \\!\\cdot\\! 23 \\!\\cdot\\! 29}{2^{2} \\!\\cdot\\! 5^{7}}}\\left(Y^{16}_{15} - Y^{16}_{-15}\\right)\n",
      "-\\sqrt{\\frac{7 \\!\\cdot\\! 17 \\!\\cdot\\! 23 \\!\\cdot\\! 31}{2 \\!\\cdot\\! 3 \\!\\cdot\\! 5^{7}}}\\left(Y^{16}_{10} + Y^{16}_{-10}\\right)\n",
      "-\\sqrt{\\frac{13 \\!\\cdot\\! 17 \\!\\cdot\\! 31}{2^{2} \\!\\cdot\\! 5^{6}}}\\left(Y^{16}_{5} - Y^{16}_{-5}\\right)\n",
      "+\\sqrt{\\frac{2^{4} \\!\\cdot\\! 19 \\!\\cdot\\! 31}{3 \\!\\cdot\\! 5^{6}}}Y^{16}_{0}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "v32=get_normalized_coefficients(F_32(x,y),32)\n",
    "show_as_spherical_harmonics_superposition(v32, 32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+\\sqrt{\\frac{17 \\!\\cdot\\! 19 \\!\\cdot\\! 37}{2 \\!\\cdot\\! 5^{7}}}\\left(Y^{21}_{20} - Y^{21}_{-20}\\right)\n",
      "+\\sqrt{\\frac{3 \\!\\cdot\\! 13 \\!\\cdot\\! 17 \\!\\cdot\\! 41}{2 \\!\\cdot\\! 5^{7}}}\\left(Y^{21}_{15} + Y^{21}_{-15}\\right)\n",
      "-\\sqrt{\\frac{2 \\!\\cdot\\! 13 \\!\\cdot\\! 41}{5^{7}}}\\left(Y^{21}_{10} - Y^{21}_{-10}\\right)\n",
      "-\\sqrt{\\frac{29 \\!\\cdot\\! 31 \\!\\cdot\\! 41}{2 \\!\\cdot\\! 5^{7}}}\\left(Y^{21}_{5} + Y^{21}_{-5}\\right)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "v42=get_normalized_coefficients(F_42(x,y),42)\n",
    "show_as_spherical_harmonics_superposition(v42, 42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+\\sqrt{\\frac{11 \\!\\cdot\\! 23 \\!\\cdot\\! 41 \\!\\cdot\\! 47}{3 \\!\\cdot\\! 5^{9}}}\\left(Y^{25}_{25} + Y^{25}_{-25}\\right)\n",
      "+\\sqrt{\\frac{3 \\!\\cdot\\! 7^{2} \\!\\cdot\\! 11 \\!\\cdot\\! 41}{2 \\!\\cdot\\! 5^{10}}}\\left(Y^{25}_{20} - Y^{25}_{-20}\\right)\n",
      "-\\sqrt{\\frac{7^{2} \\!\\cdot\\! 43 \\!\\cdot\\! 61^{2}}{2 \\!\\cdot\\! 3 \\!\\cdot\\! 5^{10}}}\\left(Y^{25}_{15} + Y^{25}_{-15}\\right)\n",
      "+\\sqrt{\\frac{7 \\!\\cdot\\! 11 \\!\\cdot\\! 19 \\!\\cdot\\! 37 \\!\\cdot\\! 43}{5^{10}}}\\left(Y^{25}_{10} - Y^{25}_{-10}\\right)\n",
      "-\\sqrt{\\frac{7^{2} \\!\\cdot\\! 31 \\!\\cdot\\! 37 \\!\\cdot\\! 43}{2 \\!\\cdot\\! 3 \\!\\cdot\\! 5^{10}}}\\left(Y^{25}_{5} + Y^{25}_{-5}\\right)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "v50=get_normalized_coefficients(F_50(x,y),50)\n",
    "show_as_spherical_harmonics_superposition(v50, 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+\\sqrt{\\frac{19 \\!\\cdot\\! 29 \\!\\cdot\\! 47 \\!\\cdot\\! 53 \\!\\cdot\\! 59}{2^{2} \\!\\cdot\\! 5^{12}}}\\left(Y^{31}_{30} - Y^{31}_{-30}\\right)\n",
      "+\\sqrt{\\frac{47 \\!\\cdot\\! 53 \\!\\cdot\\! 61}{2 \\!\\cdot\\! 5^{10}}}\\left(Y^{31}_{25} + Y^{31}_{-25}\\right)\n",
      "-\\sqrt{\\frac{13 \\!\\cdot\\! 47 \\!\\cdot\\! 61 \\!\\cdot\\! 79^{2}}{2 \\!\\cdot\\! 3 \\!\\cdot\\! 5^{12}}}\\left(Y^{31}_{20} - Y^{31}_{-20}\\right)\n",
      "+\\sqrt{\\frac{7 \\!\\cdot\\! 17 \\!\\cdot\\! 31^{2} \\!\\cdot\\! 61}{2 \\!\\cdot\\! 3 \\!\\cdot\\! 5^{11}}}\\left(Y^{31}_{15} + Y^{31}_{-15}\\right)\n",
      "+\\sqrt{\\frac{7 \\!\\cdot\\! 11 \\!\\cdot\\! 19 \\!\\cdot\\! 23 \\!\\cdot\\! 43 \\!\\cdot\\! 61}{2^{2} \\!\\cdot\\! 3 \\!\\cdot\\! 5^{11}}}\\left(Y^{31}_{10} - Y^{31}_{-10}\\right)\n",
      "-\\sqrt{\\frac{2 \\!\\cdot\\! 7 \\!\\cdot\\! 37 \\!\\cdot\\! 41 \\!\\cdot\\! 43 \\!\\cdot\\! 61}{3 \\!\\cdot\\! 5^{12}}}\\left(Y^{31}_{5} + Y^{31}_{-5}\\right)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "v62=get_normalized_coefficients(F_62(x,y),62)\n",
    "show_as_spherical_harmonics_superposition(v62, 62)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAEuCAYAAAAwQP9DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAD9ElEQVR4nO3UQQ0AIBDAMMC/50MDL7KkVbDX9swsgJLzOwDglXEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAjnEBOcYF5BgXkGNcQI5xATnGBeQYF5BjXECOcQE5xgXkGBeQY1xAzgVxXwVZtTeflwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#https://scipython.com/book/chapter-8-scipy/examples/visualizing-the-spherical-harmonics/\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "from matplotlib import cm, colors\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import numpy as np\n",
    "from numpy.core.records import array\n",
    "from scipy.special import sph_harm\n",
    "\n",
    "def superpose(coefficients, theta, phi):\n",
    "    val = 0\n",
    "    for t in coefficients:\n",
    "        l,m,c = t\n",
    "        val += c*sph_harm(m,l,theta,phi)\n",
    "    return val\n",
    "\n",
    "fig = plt.figure(figsize=plt.figaspect(1.), facecolor='none')\n",
    "fig.subplots_adjust(left=0, right=1, bottom=0, top=1)\n",
    "\n",
    "ax = fig.add_subplot(111, projection='3d', facecolor='none')\n",
    "ax.set_axis_off()\n",
    "ax.set_box_aspect((1,1,1)) #https://sabopy.com/py/matplotlib-3d-50/\n",
    "\n",
    "fcolors = np.array([])\n",
    "#resolution=dict(phi=100, theta=100)\n",
    "\n",
    "def plot_spherical_harmonics(ax, coefficients, deg, elev=0.0, azim=0.0, is_first=True):\n",
    "    global fcolors\n",
    "    #global resolution\n",
    "    phi = np.linspace(0, np.pi, 360)\n",
    "    theta = np.linspace(0, 2*np.pi, 180)\n",
    "    phi, theta = np.meshgrid(phi, theta)\n",
    "\n",
    "    x = np.sin(phi) * np.cos(theta)\n",
    "    y = np.sin(phi) * np.sin(theta)\n",
    "    z = np.cos(phi)\n",
    "\n",
    "    if(is_first):\n",
    "        if(deg%4==0):\n",
    "            fcolors = superpose(coefficients, theta, phi).real\n",
    "        elif(deg%4==2):\n",
    "            fcolors = superpose(coefficients, theta, phi).imag\n",
    "        \n",
    "    fmax, fmin = fcolors.max(), fcolors.min()\n",
    "    #fcolors = (fcolors - fmin)/(fmax - fmin)\n",
    "    norm = colors.Normalize(vmin=fmin, vmax=fmax)\n",
    "\n",
    "    ax.view_init(elev=elev, azim=azim)\n",
    "    ax.plot_surface(x, y, z,  rstride=1, cstride=1, facecolors=plt.cm.jet(norm(fcolors)), shade=False)\n",
    "    plt.plot(1,1)\n",
    "\n",
    "interval = 150\n",
    "omega = 12\n",
    "frames = 6\n",
    "    \n",
    "def plot(i, deg, coefficients):\n",
    "    print(f\"{i+1}/{frames}...\")\n",
    "    plot_spherical_harmonics(ax, coefficients, deg, elev=30, azim=omega*i, is_first=True if i==0 else False)\n",
    "\n",
    "def animate(f, deg, file=\"out.gif\"):\n",
    "    v = get_normalized_coefficients(f, deg)\n",
    "    coefficients = []\n",
    "    for i, c in enumerate(v):\n",
    "        coefficients.append((deg//2, deg//2-i, c.n()))\n",
    "    ani = animation.FuncAnimation(fig, lambda i:plot(i,deg,coefficients), interval=interval, frames=frames)\n",
    "    #ani.repeat()\n",
    "    ani.save(file, writer=\"imagemagick\", savefig_kwargs=dict(facecolor=\"black\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "#F = [(12,F_12(x,y)), (20,F_20(x,y)),(30, F_30(x,y)), (32, F_32(x,y)), (42, F_42(x,y)), (50, F_50(x,y)), (62, F_62(x,y))]\n",
    "#for a in F:\n",
    "#    deg, f = a\n",
    "#    animate(f,deg,file=f\"rotating-icosahedron-deg-{deg}-e.gif\")"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}