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

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 『正多面体と素数』10.2.3の計算ノート"
   ]
  },
  {
   "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^{12} - 33 \\, X^{6} + 11 \\, \\sqrt{5} {\\left(X^{9} - X^{3}\\right)} + 1</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}X^{12} - 33 \\, X^{6} + 11 \\, \\sqrt{5} {\\left(X^{9} - X^{3}\\right)} + 1\n",
       "\\end{math}"
      ],
      "text/plain": [
       "X^12 - 33*X^6 + 11*sqrt(5)*(X^9 - X^3) + 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "F_12(x,y)=x^12+y^12+11*sqrt(5)*(x^9*y^3-x^3*y^9)-33*x^6*y^6\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}}8712 \\, \\sqrt{5} x^{19} y - 248292 \\, \\sqrt{5} x^{13} y^{7} + 248292 \\, \\sqrt{5} x^{7} y^{13} - 8712 \\, \\sqrt{5} x y^{19} - 310365 \\, x^{16} y^{4} - 2689830 \\, x^{10} y^{10} - 310365 \\, x^{4} y^{16}</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}8712 \\, \\sqrt{5} x^{19} y - 248292 \\, \\sqrt{5} x^{13} y^{7} + 248292 \\, \\sqrt{5} x^{7} y^{13} - 8712 \\, \\sqrt{5} x y^{19} - 310365 \\, x^{16} y^{4} - 2689830 \\, x^{10} y^{10} - 310365 \\, x^{4} y^{16}\n",
       "\\end{math}"
      ],
      "text/plain": [
       "8712*sqrt(5)*x^19*y - 248292*sqrt(5)*x^13*y^7 + 248292*sqrt(5)*x^7*y^13 - 8712*sqrt(5)*x*y^19 - 310365*x^16*y^4 - 2689830*x^10*y^10 - 310365*x^4*y^16"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(F_20_0(x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-\\frac{57}{8} \\, \\sqrt{5} x^{16} y^{4} - \\frac{247}{4} \\, \\sqrt{5} x^{10} y^{10} - \\frac{57}{8} \\, \\sqrt{5} x^{4} y^{16} + x^{19} y - \\frac{57}{2} \\, x^{13} y^{7} + \\frac{57}{2} \\, x^{7} y^{13} - x y^{19}"
      ]
     },
     "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^19*y)).expand().simplify_full()\n",
    "latex(F_20(x,y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_{20}(x,y)=-\\frac{57}{8} \\, \\sqrt{5} x^{16} y^{4} - \\frac{247}{4} \\, \\sqrt{5} x^{10} y^{10} - \\frac{57}{8} \\, \\sqrt{5} x^{4} y^{16} + x^{19} y - \\frac{57}{2} \\, x^{13} y^{7} + \\frac{57}{2} \\, x^{7} y^{13} - x y^{19}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8 \\, X^{19} - 57 \\, \\sqrt{5} X^{16} - 228 \\, X^{13} - 494 \\, \\sqrt{5} X^{10} + 228 \\, X^{7} - 57 \\, \\sqrt{5} X^{4} - 8 \\, X"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latex(8*F_20(X,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "『正多面体と素数』p.133\n",
    "\n",
    "$$\n",
    "H(X,1)=F_R^{\\{1\\}}(X)=8 \\, X^{19} - 57 \\, \\sqrt{5} X^{16} - 228 \\, X^{13} - 494 \\, \\sqrt{5} X^{10} + 228 \\, X^{7} - 57 \\, \\sqrt{5} X^{4} - 8 \\, X\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}}-870 \\, \\sqrt{5} x^{27} y^{3} - 30015 \\, \\sqrt{5} x^{21} y^{9} + 30015 \\, \\sqrt{5} x^{9} y^{21} + 870 \\, \\sqrt{5} x^{3} y^{27} + 12 \\, x^{30} + \\frac{11745}{2} \\, x^{24} y^{6} - \\frac{570285}{2} \\, x^{18} y^{12} - \\frac{570285}{2} \\, x^{12} y^{18} + \\frac{11745}{2} \\, x^{6} y^{24} + 12 \\, y^{30}</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}-870 \\, \\sqrt{5} x^{27} y^{3} - 30015 \\, \\sqrt{5} x^{21} y^{9} + 30015 \\, \\sqrt{5} x^{9} y^{21} + 870 \\, \\sqrt{5} x^{3} y^{27} + 12 \\, x^{30} + \\frac{11745}{2} \\, x^{24} y^{6} - \\frac{570285}{2} \\, x^{18} y^{12} - \\frac{570285}{2} \\, x^{12} y^{18} + \\frac{11745}{2} \\, x^{6} y^{24} + 12 \\, y^{30}\n",
       "\\end{math}"
      ],
      "text/plain": [
       "-870*sqrt(5)*x^27*y^3 - 30015*sqrt(5)*x^21*y^9 + 30015*sqrt(5)*x^9*y^21 + 870*sqrt(5)*x^3*y^27 + 12*x^30 + 11745/2*x^24*y^6 - 570285/2*x^18*y^12 - 570285/2*x^12*y^18 + 11745/2*x^6*y^24 + 12*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": [
       "-\\frac{145}{2} \\, \\sqrt{5} x^{27} y^{3} - \\frac{10005}{4} \\, \\sqrt{5} x^{21} y^{9} + \\frac{10005}{4} \\, \\sqrt{5} x^{9} y^{21} + \\frac{145}{2} \\, \\sqrt{5} x^{3} y^{27} + x^{30} + \\frac{3915}{8} \\, x^{24} y^{6} - \\frac{190095}{8} \\, x^{18} y^{12} - \\frac{190095}{8} \\, x^{12} y^{18} + \\frac{3915}{8} \\, x^{6} y^{24} + 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)=-\\frac{145}{2} \\, \\sqrt{5} x^{27} y^{3} - \\frac{10005}{4} \\, \\sqrt{5} x^{21} y^{9} + \\frac{10005}{4} \\, \\sqrt{5} x^{9} y^{21} + \\frac{145}{2} \\, \\sqrt{5} x^{3} y^{27} + x^{30} + \\frac{3915}{8} \\, x^{24} y^{6} - \\frac{190095}{8} \\, x^{18} y^{12} - \\frac{190095}{8} \\, x^{12} y^{18} + \\frac{3915}{8} \\, x^{6} y^{24} + y^{30}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8 \\, X^{30} - 580 \\, \\sqrt{5} X^{27} + 3915 \\, X^{24} - 20010 \\, \\sqrt{5} X^{21} - 190095 \\, X^{18} - 190095 \\, X^{12} + 20010 \\, \\sqrt{5} X^{9} + 3915 \\, X^{6} + 580 \\, \\sqrt{5} X^{3} + 8"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latex(8*F_30(X,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "『正多面体と素数』p.134\n",
    "\n",
    "$$\n",
    "T(X,1)=F_R^{\\{1\\}}(X)=8 \\, X^{30} - 580 \\, \\sqrt{5} X^{27} + 3915 \\, X^{24} - 20010 \\, \\sqrt{5} X^{21} - 190095 \\, X^{18} - 190095 \\, X^{12} + 20010 \\, \\sqrt{5} X^{9} + 3915 \\, X^{6} + 580 \\, \\sqrt{5} X^{3} + 8\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[200*sqrt(5), 0]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#多項式として除算\n",
    "(F_12(X,1)^5-F_30(X,1)^2).maxima_methods().divide((F_20(X,1)^3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "すなわち、\n",
    "$$\n",
    "F_{12}^5 - F_{30}^2+200\\sqrt{5}\\,F_{20}^3=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^{12}y^{0} + x^{0}y^{12}\\right)\n",
      "+11\\sqrt{5}\\left(x^{9}y^{3} - x^{3}y^{9}\\right)\n",
      "-3 \\!\\cdot\\! 11x^{6}y^{6}\n"
     ]
    }
   ],
   "source": [
    "show_polynomial(F_12(x,y), 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_{12}(x,y)=&+1\\left(x^{12}y^{0} + x^{0}y^{12}\\right)\n",
    "\\\\&+11\\sqrt{5}\\left(x^{9}y^{3} - x^{3}y^{9}\\right)\n",
    "\\\\&-3 \\!\\cdot\\! 11x^{6}y^{6}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+1\\left(x^{19}y^{1} - x^{1}y^{19}\\right)\n",
      "-\\frac{3 \\!\\cdot\\! 19\\sqrt{5}}{2^3} \\left(x^{16}y^{4} + x^{4}y^{16}\\right)\n",
      "-\\frac{3 \\!\\cdot\\! 19}{2} \\left(x^{13}y^{7} - x^{7}y^{13}\\right)\n",
      "-\\frac{13 \\!\\cdot\\! 19\\sqrt{5}}{2^2} x^{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",
      "-\\frac{5 \\!\\cdot\\! 29\\sqrt{5}}{2} \\left(x^{27}y^{3} - x^{3}y^{27}\\right)\n",
      "+\\frac{3^3 \\!\\cdot\\! 5 \\!\\cdot\\! 29}{2^3} \\left(x^{24}y^{6} + x^{6}y^{24}\\right)\n",
      "-\\frac{3 \\!\\cdot\\! 5 \\!\\cdot\\! 23 \\!\\cdot\\! 29\\sqrt{5}}{2^2} \\left(x^{21}y^{9} - x^{9}y^{21}\\right)\n",
      "-\\frac{3 \\!\\cdot\\! 5 \\!\\cdot\\! 19 \\!\\cdot\\! 23 \\!\\cdot\\! 29}{2^3} \\left(x^{18}y^{12} + x^{12}y^{18}\\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",
      "+\\frac{31\\sqrt{5}}{2^3} \\left(x^{28}y^{4} + x^{4}y^{28}\\right)\n",
      "-\\frac{3^2 \\!\\cdot\\! 13 \\!\\cdot\\! 31}{2^3} \\left(x^{25}y^{7} - x^{7}y^{25}\\right)\n",
      "-\\frac{3 \\!\\cdot\\! 13 \\!\\cdot\\! 31\\sqrt{5}}{2^3} \\left(x^{22}y^{10} + x^{10}y^{22}\\right)\n",
      "-\\frac{3 \\!\\cdot\\! 5^2 \\!\\cdot\\! 7 \\!\\cdot\\! 31}{2^3} \\left(x^{19}y^{13} - x^{13}y^{19}\\right)\n",
      "+\\frac{3^2 \\!\\cdot\\! 19 \\!\\cdot\\! 31\\sqrt{5}}{2} x^{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^{42}y^{0} + x^{0}y^{42}\\right)\n",
      "-\\frac{3 \\!\\cdot\\! 41\\sqrt{5}}{2} \\left(x^{39}y^{3} - x^{3}y^{39}\\right)\n",
      "-\\frac{13 \\!\\cdot\\! 41 \\!\\cdot\\! 53}{2^3} \\left(x^{36}y^{6} + x^{6}y^{36}\\right)\n",
      "+\\frac{13 \\!\\cdot\\! 41 \\!\\cdot\\! 79\\sqrt{5}}{2^3} \\left(x^{33}y^{9} - x^{9}y^{33}\\right)\n",
      "-\\frac{3 \\!\\cdot\\! 7 \\!\\cdot\\! 13 \\!\\cdot\\! 31 \\!\\cdot\\! 41}{2} \\left(x^{30}y^{12} + x^{12}y^{30}\\right)\n",
      "-5 \\!\\cdot\\! 29 \\!\\cdot\\! 31 \\!\\cdot\\! 41\\sqrt{5}\\left(x^{27}y^{15} - x^{15}y^{27}\\right)\n",
      "+\\frac{3 \\!\\cdot\\! 5 \\!\\cdot\\! 13 \\!\\cdot\\! 29 \\!\\cdot\\! 31 \\!\\cdot\\! 41}{2^3} \\left(x^{24}y^{18} + x^{18}y^{24}\\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^{49}y^{1} - x^{1}y^{49}\\right)\n",
      "-\\frac{7^2 \\!\\cdot\\! 13\\sqrt{5}}{2^3} \\left(x^{46}y^{4} + x^{4}y^{46}\\right)\n",
      "+\\frac{3^2 \\!\\cdot\\! 7 \\!\\cdot\\! 773}{2^4} \\left(x^{43}y^{7} - x^{7}y^{43}\\right)\n",
      "-\\frac{7^2 \\!\\cdot\\! 11^2 \\!\\cdot\\! 43\\sqrt{5}}{2^6} \\left(x^{40}y^{10} + x^{10}y^{40}\\right)\n",
      "+\\frac{7^2 \\!\\cdot\\! 19 \\!\\cdot\\! 43 \\!\\cdot\\! 59}{2^5} \\left(x^{37}y^{13} - x^{13}y^{37}\\right)\n",
      "+\\frac{3^2 \\!\\cdot\\! 7^2 \\!\\cdot\\! 19 \\!\\cdot\\! 37 \\!\\cdot\\! 43\\sqrt{5}}{2^6} \\left(x^{34}y^{16} + x^{16}y^{34}\\right)\n",
      "+\\frac{7^2 \\!\\cdot\\! 37^2 \\!\\cdot\\! 43}{2} \\left(x^{31}y^{19} - x^{19}y^{31}\\right)\n",
      "+\\frac{5 \\!\\cdot\\! 7 \\!\\cdot\\! 29 \\!\\cdot\\! 31 \\!\\cdot\\! 37 \\!\\cdot\\! 43\\sqrt{5}}{2^5} \\left(x^{28}y^{22} + x^{22}y^{28}\\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",
      "-\\frac{3^2 \\!\\cdot\\! 61\\sqrt{5}}{2^3} \\left(x^{58}y^{4} + x^{4}y^{58}\\right)\n",
      "-\\frac{61 \\!\\cdot\\! 359}{2^4} \\left(x^{55}y^{7} - x^{7}y^{55}\\right)\n",
      "+\\frac{3 \\!\\cdot\\! 11 \\!\\cdot\\! 61 \\!\\cdot\\! 1021\\sqrt{5}}{2^6} \\left(x^{52}y^{10} + x^{10}y^{52}\\right)\n",
      "-\\frac{3^2 \\!\\cdot\\! 5 \\!\\cdot\\! 17 \\!\\cdot\\! 61 \\!\\cdot\\! 331}{2^6} \\left(x^{49}y^{13} - x^{13}y^{49}\\right)\n",
      "+\\frac{3 \\!\\cdot\\! 7 \\!\\cdot\\! 17 \\!\\cdot\\! 31 \\!\\cdot\\! 53 \\!\\cdot\\! 61\\sqrt{5}}{2^5} \\left(x^{46}y^{16} + x^{16}y^{46}\\right)\n",
      "+\\frac{3 \\!\\cdot\\! 5 \\!\\cdot\\! 7 \\!\\cdot\\! 11 \\!\\cdot\\! 23 \\!\\cdot\\! 61 \\!\\cdot\\! 211}{2^5} \\left(x^{43}y^{19} - x^{19}y^{43}\\right)\n",
      "+\\frac{3^2 \\!\\cdot\\! 7 \\!\\cdot\\! 23 \\!\\cdot\\! 41 \\!\\cdot\\! 43 \\!\\cdot\\! 61\\sqrt{5}}{2^4} \\left(x^{40}y^{22} + x^{22}y^{40}\\right)\n",
      "+\\frac{3 \\!\\cdot\\! 7 \\!\\cdot\\! 13 \\!\\cdot\\! 41 \\!\\cdot\\! 43 \\!\\cdot\\! 59 \\!\\cdot\\! 61}{2^6} \\left(x^{37}y^{25} - x^{25}y^{37}\\right)\n",
      "-\\frac{7 \\!\\cdot\\! 37 \\!\\cdot\\! 41 \\!\\cdot\\! 43 \\!\\cdot\\! 61 \\!\\cdot\\! 151\\sqrt{5}}{2^6} \\left(x^{34}y^{28} + x^{28}y^{34}\\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(\\frac{2}{297} \\, \\sqrt{231} \\sqrt{11},\\,0,\\,0,\\,\\frac{1}{135} \\, \\sqrt{105} \\sqrt{11} \\sqrt{5},\\,0,\\,0,\\,-\\frac{1}{9} \\, \\sqrt{11},\\,0,\\,0,\\,-\\frac{1}{135} \\, \\sqrt{105} \\sqrt{11} \\sqrt{5},\\,0,\\,0,\\,\\frac{2}{297} \\, \\sqrt{231} \\sqrt{11}\\right)</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\frac{2}{297} \\, \\sqrt{231} \\sqrt{11},\\,0,\\,0,\\,\\frac{1}{135} \\, \\sqrt{105} \\sqrt{11} \\sqrt{5},\\,0,\\,0,\\,-\\frac{1}{9} \\, \\sqrt{11},\\,0,\\,0,\\,-\\frac{1}{135} \\, \\sqrt{105} \\sqrt{11} \\sqrt{5},\\,0,\\,0,\\,\\frac{2}{297} \\, \\sqrt{231} \\sqrt{11}\\right)\n",
       "\\end{math}"
      ],
      "text/plain": [
       "(2/297*sqrt(231)*sqrt(11), 0, 0, 1/135*sqrt(105)*sqrt(11)*sqrt(5), 0, 0, -1/9*sqrt(11), 0, 0, -1/135*sqrt(105)*sqrt(11)*sqrt(5), 0, 0, 2/297*sqrt(231)*sqrt(11))"
      ]
     },
     "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 reversed(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{2^{2} \\!\\cdot\\! 7}{3^{5}}}\\left(Y^{6}_{6} + Y^{6}_{-6}\\right)\n",
      "-\\sqrt{\\frac{7 \\!\\cdot\\! 11}{3^{5}}}\\left(Y^{6}_{3} - Y^{6}_{-3}\\right)\n",
      "-\\sqrt{\\frac{11}{3^{4}}}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{2^{2} \\!\\cdot\\! 7}{3^{5}}}\\left(Y^{6}_{6} + Y^{6}_{-6}\\right)\n",
    "-\\sqrt{\\frac{7 \\!\\cdot\\! 11}{3^{5}}}\\left(Y^{6}_{3} - Y^{6}_{-3}\\right)\n",
    "-\\sqrt{\\frac{11}{3^{4}}}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(0,\\,\\frac{4}{300105} \\, \\sqrt{230945} \\sqrt{3705},\\,0,\\,0,\\,-\\frac{1}{15795} \\, \\sqrt{3705} \\sqrt{2145} \\sqrt{5},\\,0,\\,0,\\,-\\frac{1}{15795} \\, \\sqrt{3705} \\sqrt{2145},\\,0,\\,0,\\,-\\frac{1}{243} \\, \\sqrt{3705} \\sqrt{5},\\,0,\\,0,\\,\\frac{1}{15795} \\, \\sqrt{3705} \\sqrt{2145},\\,0,\\,0,\\,-\\frac{1}{15795} \\, \\sqrt{3705} \\sqrt{2145} \\sqrt{5},\\,0,\\,0,\\,-\\frac{4}{300105} \\, \\sqrt{230945} \\sqrt{3705},\\,0\\right)</script></html>"
      ],
      "text/latex": [
       "\\begin{math}\n",
       "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(0,\\,\\frac{4}{300105} \\, \\sqrt{230945} \\sqrt{3705},\\,0,\\,0,\\,-\\frac{1}{15795} \\, \\sqrt{3705} \\sqrt{2145} \\sqrt{5},\\,0,\\,0,\\,-\\frac{1}{15795} \\, \\sqrt{3705} \\sqrt{2145},\\,0,\\,0,\\,-\\frac{1}{243} \\, \\sqrt{3705} \\sqrt{5},\\,0,\\,0,\\,\\frac{1}{15795} \\, \\sqrt{3705} \\sqrt{2145},\\,0,\\,0,\\,-\\frac{1}{15795} \\, \\sqrt{3705} \\sqrt{2145} \\sqrt{5},\\,0,\\,0,\\,-\\frac{4}{300105} \\, \\sqrt{230945} \\sqrt{3705},\\,0\\right)\n",
       "\\end{math}"
      ],
      "text/plain": [
       "(0, 4/300105*sqrt(230945)*sqrt(3705), 0, 0, -1/15795*sqrt(3705)*sqrt(2145)*sqrt(5), 0, 0, -1/15795*sqrt(3705)*sqrt(2145), 0, 0, -1/243*sqrt(3705)*sqrt(5), 0, 0, 1/15795*sqrt(3705)*sqrt(2145), 0, 0, -1/15795*sqrt(3705)*sqrt(2145)*sqrt(5), 0, 0, -4/300105*sqrt(230945)*sqrt(3705), 0)"
      ]
     },
     "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{2^{4} \\!\\cdot\\! 11 \\!\\cdot\\! 17}{3^{9}}}\\left(Y^{10}_{9} - Y^{10}_{-9}\\right)\n",
      "-\\sqrt{\\frac{5 \\!\\cdot\\! 11 \\!\\cdot\\! 19}{3^{8}}}\\left(Y^{10}_{6} + Y^{10}_{-6}\\right)\n",
      "+\\sqrt{\\frac{11 \\!\\cdot\\! 19}{3^{8}}}\\left(Y^{10}_{3} - Y^{10}_{-3}\\right)\n",
      "-\\sqrt{\\frac{5^{2} \\!\\cdot\\! 13 \\!\\cdot\\! 19}{3^{9}}}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{2^{4} \\!\\cdot\\! 7 \\!\\cdot\\! 11 \\!\\cdot\\! 13}{3^{12}}}\\left(Y^{15}_{15} + Y^{15}_{-15}\\right)\n",
      "+\\sqrt{\\frac{5^{2} \\!\\cdot\\! 11 \\!\\cdot\\! 13 \\!\\cdot\\! 29}{3^{12}}}\\left(Y^{15}_{12} - Y^{15}_{-12}\\right)\n",
      "+\\sqrt{\\frac{11 \\!\\cdot\\! 29}{2^{2} \\!\\cdot\\! 3^{8}}}\\left(Y^{15}_{9} + Y^{15}_{-9}\\right)\n",
      "+\\sqrt{\\frac{5 \\!\\cdot\\! 7 \\!\\cdot\\! 23 \\!\\cdot\\! 29}{2 \\!\\cdot\\! 3^{11}}}\\left(Y^{15}_{6} - Y^{15}_{-6}\\right)\n",
      "-\\sqrt{\\frac{11 \\!\\cdot\\! 19 \\!\\cdot\\! 23 \\!\\cdot\\! 29}{2^{2} \\!\\cdot\\! 3^{11}}}\\left(Y^{15}_{3} + Y^{15}_{-3}\\right)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "show_as_spherical_harmonics_superposition(w, 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "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 = 10\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": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "#animate(F_12(x,y),12,file='rotating-icosahedron-deg-12-f.gif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "#animate(F_20(x,y),20,file='rotating-icosahedron-deg-20-f.gif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "#animate(F_30(x,y),30,file='rotating-icosahedron-deg-30-f.gif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "#animate(F_62(x,y),62,file='rotating-icosahedron-deg-62-f.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
}