{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:black; font-weight: bold; font-size:26px\">\n",
    "    Labeled digraphs with (strongly) connected components\n",
    "</span>\n",
    "\n",
    "This is a supplementary file for the paper\n",
    "    \n",
    "    \"Asymptotics for graphically divergent series: dense digraphs and 2-SAT formulae\"\n",
    "    by Sergey Dovgal and Khaydar Nurligareev.\n",
    "    \n",
    "Here, you can find the code for obtaining asymptotic coefficients from Section A.10, Tables 18-22"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We mark\n",
    "<ul>\n",
    "    <li>\n",
    "        all strongly connected components by variable t,\n",
    "    </li>\n",
    "    <li>\n",
    "        purely source-like strongly connected components by variable u,\n",
    "    </li>\n",
    "    <li>\n",
    "        purely sink-like strongly connected components by variable v,\n",
    "    </li>\n",
    "    <li>\n",
    "        isolated strongly connected components by variable y.\n",
    "    </li>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:black; font-weight: bold; font-size:18px\">\n",
    "    Preliminary section\n",
    "</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 7 # Our accuracy\n",
    "disp = 7 # How many terms we want to display\n",
    "P.<a> = PolynomialRing(QQ)\n",
    "PU.<u> = PolynomialRing(P)\n",
    "PV.<v> = PolynomialRing(PU)\n",
    "PY.<y> = PolynomialRing(PV)\n",
    "PT.<t> = PolynomialRing(PY)\n",
    "PW.<w> = PolynomialRing(PT)\n",
    "R.<z> = PowerSeriesRing(PW,N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exponential Hadamard product of two series\n",
    "# order is the accuracy of this operation\n",
    "def exp_had_prod(f, g, order):\n",
    "    return sum(f[n] * g[n] * factorial(n) * z^n for n in range(order))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, 1/a^2, 1/a^6, 1/a^12, 1/a^20, 1/a^30]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Transformation function hat(set)\n",
    "hatset = sum(((1 / a^(i*(i-1))) / i.factorial()) * z^i for i in srange(N))\n",
    "HATSET = [hatset[i] * i.factorial() for i in srange(disp)]\n",
    "HATSET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, 1/a, 1/a^3, 1/a^6, 1/a^10, 1/a^15]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Transformation function Phi_{a^2}^{2,1}\n",
    "PhiTwoTwoOne = sum((1 / a^(i*(i-1)/2) / i.factorial()) * z^i for i in srange(N))\n",
    "PTTO = [PhiTwoTwoOne[i] * i.factorial() for i in srange(disp)]\n",
    "PTTO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, a, a^3, a^6, a^10, a^15]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Transformation function Phi_{a^2}^{1,2}\n",
    "PhiTwoOneTwo = sum((a^(i*(i-1)/2) / i.factorial()) * z^i for i in srange(N))\n",
    "PTOT = [PhiTwoOneTwo[i] * i.factorial() for i in srange(disp)]\n",
    "PTOT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:black; font-weight: bold; font-size:18px\">\n",
    "    Defining supplementary EGFs\n",
    "</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "disp = 5 # How many terms we want to display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, a^2, a^6, a^12]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of labeled graphs/tournaments or GGF of labeled digraphs (we take a parameter a = sqrt{2})\n",
    "g = sum((a^(i*(i-1)) / i.factorial()) * z^i for i in srange(N))\n",
    "G = [g[i] * i.factorial() for i in srange(disp)]\n",
    "G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0,\n",
       " 1,\n",
       " a^4 - 2*a^2 + 1,\n",
       " a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2,\n",
       " a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of strongly connected labeled digraphs (we take a parameter a = sqrt{2})\n",
    "scd = -log(exp_had_prod(g,1/g,N))\n",
    "SCD = [scd[i] * i.factorial() for i in srange(disp)]\n",
    "SCD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " -2*a^2*w,\n",
       " (-2*a^7 + 6*a^5)*w^2 + 2*a^3*w,\n",
       " (-2*a^15 + 18*a^11 - 24*a^9)*w^3 - 12*a^5*w^2,\n",
       " (-2*a^26 + 24*a^20 + 18*a^18 - 144*a^16 + 120*a^14)*w^4 + (-24*a^10 + 72*a^8)*w^3 + (12*a^10 - 24*a^8 + 24*a^6)*w^2]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CGF of strongly connected labeled digraphs (we take a parameter a = sqrt{2})\n",
    "qscd = exp(scd.subs(z=a^3*z^2*w))*exp_had_prod(PhiTwoOneTwo,(1/g.subs(z=a^2*z*w))*(1/g.subs(z=a^2*z*w)),N)\n",
    "QSCD = [qscd[i] * i.factorial() for i in srange(disp)]\n",
    "QSCD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " -4*w,\n",
       " 8*sqrt(2)*w^2 + 4*sqrt(2)*w,\n",
       " -64*sqrt(2)*w^3 - 48*sqrt(2)*w^2,\n",
       " -4096*w^4 + 384*w^3 + 192*w^2]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QSCDsub = [QSCD[i].subs(a=sqrt(2)) for i in srange(disp)]\n",
    "QSCDsub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1, 0, 0, 0, 0],\n",
       " [0, -4, 0, 0, 0],\n",
       " [0, 4*sqrt(2), 8*sqrt(2), 0, 0],\n",
       " [0, 0, -48*sqrt(2), -64*sqrt(2), 0],\n",
       " [0, 0, 192, 384, -4096]]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QSCDMatrix = [[QSCD[i][j].subs(a=sqrt(2)) for j in srange(disp)] for i in srange(disp)]\n",
    "QSCDMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " t,\n",
       " ((2*a^2 - 1)/a^2)*t^2 + ((a^4 - 2*a^2 + 1)/a^2)*t,\n",
       " ((6*a^6 - 6*a^4 + 1)/a^6)*t^3 + ((6*a^8 - 12*a^6 + 3*a^4 + 6*a^2 - 3)/a^6)*t^2 + ((a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2)/a^6)*t,\n",
       " ((24*a^12 - 36*a^10 + 6*a^8 + 8*a^6 - 1)/a^12)*t^4 + ((36*a^14 - 84*a^12 + 36*a^10 + 36*a^8 - 24*a^6 + 6*a^4 - 12*a^2 + 6)/a^12)*t^3 + ((8*a^18 + 6*a^16 - 72*a^14 + 80*a^12 - 21*a^8 + 4*a^6 - 30*a^4 + 36*a^2 - 11)/a^12)*t^2 + ((a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6)/a^12)*t]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of labeled digraphs,\n",
    "# with the marking variable t for strongly connected components (we take a parameter a = sqrt{2})\n",
    "dt = 1 / (exp_had_prod(exp(-t*scd),hatset,N))\n",
    "DT = [dt[i] * i.factorial() for i in srange(disp)]\n",
    "DT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[t,\n",
       " (2*a^2*t^2 - 2*a^2*t)*w,\n",
       " ((6*a^4 - 2*a^2)*t^3 + (2*a^6 - 12*a^4 + 2*a^2)*t^2 + (-2*a^6 + 6*a^4)*t)*w^2 + (-2*a^2*t^2 + 2*a^2*t)*w,\n",
       " ((24*a^6 - 18*a^4 + 2)*t^4 + (18*a^8 - 72*a^6 + 24*a^4 + 12*a^2 - 6)*t^3 + (2*a^12 - 36*a^8 + 72*a^6 - 6*a^4 - 12*a^2 + 4)*t^2 + (-2*a^12 + 18*a^8 - 24*a^6)*t)*w^3 + (-12*a^4*t^3 + (12*a^4 + 12*a^2)*t^2 - 12*a^2*t)*w^2,\n",
       " (((120*a^12 - 144*a^10 + 18*a^8 + 24*a^6 - 2)/a^4)*t^5 + ((144*a^14 - 516*a^12 + 288*a^10 + 108*a^8 - 88*a^6 + 12*a^4 - 24*a^2 + 12)/a^4)*t^4 + ((24*a^18 + 18*a^16 - 432*a^14 + 772*a^12 - 168*a^10 - 180*a^8 + 72*a^6 - 60*a^4 + 72*a^2 - 22)/a^4)*t^3 + ((2*a^24 - 48*a^18 - 36*a^16 + 432*a^14 - 496*a^12 + 24*a^10 + 54*a^8 - 8*a^6 + 48*a^4 - 48*a^2 + 12)/a^4)*t^2 + (-2*a^20 + 24*a^14 + 18*a^12 - 144*a^10 + 120*a^8)*t)*w^4 + ((-72*a^6 + 24*a^4)*t^4 + (-24*a^8 + 120*a^6 + 48*a^4)*t^3 + (24*a^8 - 48*a^6 - 48*a^4 - 72*a^2)*t^2 + (-24*a^4 + 72*a^2)*t)*w^3 + (12*t^3 + (-12*a^4 + 24*a^2 - 36)*t^2 + (12*a^4 - 24*a^2 + 24)*t)*w^2]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CGF of labeled digraphs,\n",
    "# with the marking variable t for strongly connected components (we take a parameter a = sqrt{2})\n",
    "num = t*exp((1-t)*scd.subs(z=a^3*z^2*w))*exp_had_prod(PhiTwoOneTwo,(1/g.subs(z=a^2*z*w))*(1/g.subs(z=a^2*z*w)),N)\n",
    "qdt = dt.subs(z=a^2*z*w)*dt.subs(z=a^2*z*w)*exp_had_prod(PhiTwoTwoOne,num,N)\n",
    "QDT = [qdt[i] * i.factorial() for i in srange(disp)]\n",
    "QDT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " t,\n",
       " t^2 + (a^4 - 2*a^2 + 1)*t,\n",
       " t^3 + (3*a^4 - 6*a^2 + 3)*t^2 + (a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2)*t,\n",
       " t^4 + (6*a^4 - 12*a^2 + 6)*t^3 + (4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*t^2 + (a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6)*t]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of semi-strong labeled digraphs,\n",
    "# with the marking variable t for strongly connected components (we take a parameter a = sqrt{2})\n",
    "ssdt = exp(t*scd)\n",
    "SSDT = [ssdt[i] * i.factorial() for i in srange(disp)]\n",
    "SSDT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " (u - 1)*t,\n",
       " (u^2 - 2*u + 1)*t^2 + ((a^4 - 2*a^2 + 1)*u - a^4 + 2*a^2 - 1)*t,\n",
       " (u^3 - 3*u^2 + 3*u - 1)*t^3 + ((3*a^4 - 6*a^2 + 3)*u^2 + (-6*a^4 + 12*a^2 - 6)*u + 3*a^4 - 6*a^2 + 3)*t^2 + ((a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2)*u - a^12 + 6*a^8 - 6*a^6 - 3*a^4 + 6*a^2 - 2)*t,\n",
       " (u^4 - 4*u^3 + 6*u^2 - 4*u + 1)*t^4 + ((6*a^4 - 12*a^2 + 6)*u^3 + (-18*a^4 + 36*a^2 - 18)*u^2 + (18*a^4 - 36*a^2 + 18)*u - 6*a^4 + 12*a^2 - 6)*t^3 + ((4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*u^2 + (-8*a^12 + 42*a^8 - 24*a^6 - 60*a^4 + 72*a^2 - 22)*u + 4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*t^2 + ((a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6)*u - a^24 + 8*a^18 + 6*a^16 - 36*a^14 + 20*a^12 + 21*a^8 - 12*a^6 - 24*a^4 + 24*a^2 - 6)*t]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Supplementary function based on the EGF of semi-strong labeled digraphs,\n",
    "# with the marking variable t for strongly connected components (we take a parameter a = sqrt{2})\n",
    "ssdtu = exp((u-1)*t*scd)\n",
    "partu = exp_had_prod(ssdtu,hatset,N)\n",
    "SSDTU = [ssdtu[i] * i.factorial() for i in srange(disp)]\n",
    "SSDTU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " (v - 1)*t,\n",
       " (v^2 - 2*v + 1)*t^2 + ((a^4 - 2*a^2 + 1)*v - a^4 + 2*a^2 - 1)*t,\n",
       " (v^3 - 3*v^2 + 3*v - 1)*t^3 + ((3*a^4 - 6*a^2 + 3)*v^2 + (-6*a^4 + 12*a^2 - 6)*v + 3*a^4 - 6*a^2 + 3)*t^2 + ((a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2)*v - a^12 + 6*a^8 - 6*a^6 - 3*a^4 + 6*a^2 - 2)*t,\n",
       " (v^4 - 4*v^3 + 6*v^2 - 4*v + 1)*t^4 + ((6*a^4 - 12*a^2 + 6)*v^3 + (-18*a^4 + 36*a^2 - 18)*v^2 + (18*a^4 - 36*a^2 + 18)*v - 6*a^4 + 12*a^2 - 6)*t^3 + ((4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*v^2 + (-8*a^12 + 42*a^8 - 24*a^6 - 60*a^4 + 72*a^2 - 22)*v + 4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*t^2 + ((a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6)*v - a^24 + 8*a^18 + 6*a^16 - 36*a^14 + 20*a^12 + 21*a^8 - 12*a^6 - 24*a^4 + 24*a^2 - 6)*t]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Supplementary function based on the EGF of semi-strong labeled digraphs,\n",
    "# with the marking variable t for strongly connected components (we take a parameter a = sqrt{2})\n",
    "ssdtv = exp((v-1)*t*scd)\n",
    "partv = exp_had_prod(ssdtv,hatset,N)\n",
    "SSDTV = [ssdtv[i] * i.factorial() for i in srange(disp)]\n",
    "SSDTV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " (y - v - u + 1)*t,\n",
       " (y^2 + (-2*v - 2*u + 2)*y + v^2 + (2*u - 2)*v + u^2 - 2*u + 1)*t^2 + ((a^4 - 2*a^2 + 1)*y + (-a^4 + 2*a^2 - 1)*v + (-a^4 + 2*a^2 - 1)*u + a^4 - 2*a^2 + 1)*t,\n",
       " (y^3 + (-3*v - 3*u + 3)*y^2 + (3*v^2 + (6*u - 6)*v + 3*u^2 - 6*u + 3)*y - v^3 + (-3*u + 3)*v^2 + (-3*u^2 + 6*u - 3)*v - u^3 + 3*u^2 - 3*u + 1)*t^3 + ((3*a^4 - 6*a^2 + 3)*y^2 + ((-6*a^4 + 12*a^2 - 6)*v + (-6*a^4 + 12*a^2 - 6)*u + 6*a^4 - 12*a^2 + 6)*y + (3*a^4 - 6*a^2 + 3)*v^2 + ((6*a^4 - 12*a^2 + 6)*u - 6*a^4 + 12*a^2 - 6)*v + (3*a^4 - 6*a^2 + 3)*u^2 + (-6*a^4 + 12*a^2 - 6)*u + 3*a^4 - 6*a^2 + 3)*t^2 + ((a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2)*y + (-a^12 + 6*a^8 - 6*a^6 - 3*a^4 + 6*a^2 - 2)*v + (-a^12 + 6*a^8 - 6*a^6 - 3*a^4 + 6*a^2 - 2)*u + a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2)*t,\n",
       " (y^4 + (-4*v - 4*u + 4)*y^3 + (6*v^2 + (12*u - 12)*v + 6*u^2 - 12*u + 6)*y^2 + (-4*v^3 + (-12*u + 12)*v^2 + (-12*u^2 + 24*u - 12)*v - 4*u^3 + 12*u^2 - 12*u + 4)*y + v^4 + (4*u - 4)*v^3 + (6*u^2 - 12*u + 6)*v^2 + (4*u^3 - 12*u^2 + 12*u - 4)*v + u^4 - 4*u^3 + 6*u^2 - 4*u + 1)*t^4 + ((6*a^4 - 12*a^2 + 6)*y^3 + ((-18*a^4 + 36*a^2 - 18)*v + (-18*a^4 + 36*a^2 - 18)*u + 18*a^4 - 36*a^2 + 18)*y^2 + ((18*a^4 - 36*a^2 + 18)*v^2 + ((36*a^4 - 72*a^2 + 36)*u - 36*a^4 + 72*a^2 - 36)*v + (18*a^4 - 36*a^2 + 18)*u^2 + (-36*a^4 + 72*a^2 - 36)*u + 18*a^4 - 36*a^2 + 18)*y + (-6*a^4 + 12*a^2 - 6)*v^3 + ((-18*a^4 + 36*a^2 - 18)*u + 18*a^4 - 36*a^2 + 18)*v^2 + ((-18*a^4 + 36*a^2 - 18)*u^2 + (36*a^4 - 72*a^2 + 36)*u - 18*a^4 + 36*a^2 - 18)*v + (-6*a^4 + 12*a^2 - 6)*u^3 + (18*a^4 - 36*a^2 + 18)*u^2 + (-18*a^4 + 36*a^2 - 18)*u + 6*a^4 - 12*a^2 + 6)*t^3 + ((4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*y^2 + ((-8*a^12 + 42*a^8 - 24*a^6 - 60*a^4 + 72*a^2 - 22)*v + (-8*a^12 + 42*a^8 - 24*a^6 - 60*a^4 + 72*a^2 - 22)*u + 8*a^12 - 42*a^8 + 24*a^6 + 60*a^4 - 72*a^2 + 22)*y + (4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*v^2 + ((8*a^12 - 42*a^8 + 24*a^6 + 60*a^4 - 72*a^2 + 22)*u - 8*a^12 + 42*a^8 - 24*a^6 - 60*a^4 + 72*a^2 - 22)*v + (4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*u^2 + (-8*a^12 + 42*a^8 - 24*a^6 - 60*a^4 + 72*a^2 - 22)*u + 4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*t^2 + ((a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6)*y + (-a^24 + 8*a^18 + 6*a^16 - 36*a^14 + 20*a^12 + 21*a^8 - 12*a^6 - 24*a^4 + 24*a^2 - 6)*v + (-a^24 + 8*a^18 + 6*a^16 - 36*a^14 + 20*a^12 + 21*a^8 - 12*a^6 - 24*a^4 + 24*a^2 - 6)*u + a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6)*t]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Supplementary function based on the EGF of semi-strong labeled digraphs,\n",
    "# with the marking variable t for strongly connected components (we take a parameter a = sqrt{2})\n",
    "ssdtall = exp((y-u-v+1)*t*scd)\n",
    "SSDTALL = [ssdtall[i] * i.factorial() for i in srange(disp)]\n",
    "SSDTALL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:black; font-weight: bold; font-size:18px\">\n",
    "    Principal section\n",
    "</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " (v + u - 1)*t,\n",
       " (1/a^2*v^2 + (2*u - 2/a^2)*v + 1/a^2*u^2 - 2/a^2*u + 1/a^2)*t^2 + (((a^4 - 2*a^2 + 1)/a^2)*v + ((a^4 - 2*a^2 + 1)/a^2)*u + (-a^4 + 2*a^2 - 1)/a^2)*t,\n",
       " (1/a^6*v^3 + (3/a^2*u - 3/a^6)*v^2 + (3/a^2*u^2 + ((6*a^2 - 12)/a^2)*u + 3/a^6)*v + 1/a^6*u^3 - 3/a^6*u^2 + 3/a^6*u - 1/a^6)*t^3 + (((3*a^4 - 6*a^2 + 3)/a^6)*v^2 + (((6*a^4 - 12*a^2 + 6)/a^2)*u + (-6*a^4 + 12*a^2 - 6)/a^6)*v + ((3*a^4 - 6*a^2 + 3)/a^6)*u^2 + ((-6*a^4 + 12*a^2 - 6)/a^6)*u + (3*a^4 - 6*a^2 + 3)/a^6)*t^2 + (((a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2)/a^6)*v + ((a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2)/a^6)*u + (-a^12 + 6*a^8 - 6*a^6 - 3*a^4 + 6*a^2 - 2)/a^6)*t,\n",
       " (1/a^12*v^4 + (4/a^6*u - 4/a^12)*v^3 + (6/a^4*u^2 + ((12*a^4 - 12*a^2 - 12)/a^6)*u + 6/a^12)*v^2 + (4/a^6*u^3 + ((12*a^4 - 12*a^2 - 12)/a^6)*u^2 + ((24*a^6 - 60*a^4 + 24*a^2 + 24)/a^6)*u - 4/a^12)*v + 1/a^12*u^4 - 4/a^12*u^3 + 6/a^12*u^2 - 4/a^12*u + 1/a^12)*t^4 + (((6*a^4 - 12*a^2 + 6)/a^12)*v^3 + (((6*a^6 - 18*a^2 + 12)/a^6)*u + (-18*a^4 + 36*a^2 - 18)/a^12)*v^2 + (((6*a^6 - 18*a^2 + 12)/a^6)*u^2 + ((36*a^8 - 96*a^6 + 36*a^4 + 72*a^2 - 48)/a^6)*u + (18*a^4 - 36*a^2 + 18)/a^12)*v + ((6*a^4 - 12*a^2 + 6)/a^12)*u^3 + ((-18*a^4 + 36*a^2 - 18)/a^12)*u^2 + ((18*a^4 - 36*a^2 + 18)/a^12)*u + (-6*a^4 + 12*a^2 - 6)/a^12)*t^3 + (((4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)/a^12)*v^2 + (((8*a^12 + 6*a^10 - 72*a^8 + 84*a^6 - 42*a^2 + 16)/a^6)*u + (-8*a^12 + 42*a^8 - 24*a^6 - 60*a^4 + 72*a^2 - 22)/a^12)*v + ((4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)/a^12)*u^2 + ((-8*a^12 + 42*a^8 - 24*a^6 - 60*a^4 + 72*a^2 - 22)/a^12)*u + (4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)/a^12)*t^2 + (((a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6)/a^12)*v + ((a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6)/a^12)*u + (-a^24 + 8*a^18 + 6*a^16 - 36*a^14 + 20*a^12 + 21*a^8 - 12*a^6 - 24*a^4 + 24*a^2 - 6)/a^12)*t]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inner = partu*partv*dt\n",
    "INNER = [inner[i] * i.factorial() for i in srange(disp)]\n",
    "INNER"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " y*t,\n",
       " (y^2 + ((2*a^2 - 2)*u)*v)*t^2 + ((a^4 - 2*a^2 + 1)*y)*t,\n",
       " (y^3 + (((6*a^2 - 6)*u)*v)*y + ((3*a^4 - 6*a^2 + 3)*u)*v^2 + ((3*a^4 - 6*a^2 + 3)*u^2 + (6*a^6 - 12*a^4 + 6*a^2)*u)*v)*t^3 + ((3*a^4 - 6*a^2 + 3)*y^2 + ((6*a^8 - 12*a^6 + 12*a^2 - 6)*u)*v)*t^2 + ((a^12 - 6*a^8 + 6*a^6 + 3*a^4 - 6*a^2 + 2)*y)*t,\n",
       " (y^4 + (((12*a^2 - 12)*u)*v)*y^2 + (((12*a^4 - 24*a^2 + 12)*u)*v^2 + ((12*a^4 - 24*a^2 + 12)*u^2 + (24*a^6 - 48*a^4 + 24*a^2)*u)*v)*y + ((4*a^6 - 12*a^4 + 12*a^2 - 4)*u)*v^3 + ((6*a^8 - 24*a^4 + 24*a^2 - 6)*u^2 + (12*a^10 - 12*a^8 - 36*a^6 + 60*a^4 - 24*a^2)*u)*v^2 + ((4*a^6 - 12*a^4 + 12*a^2 - 4)*u^3 + (12*a^10 - 12*a^8 - 36*a^6 + 60*a^4 - 24*a^2)*u^2 + (24*a^12 - 60*a^10 + 24*a^8 + 48*a^6 - 48*a^4 + 12*a^2)*u)*v)*t^4 + ((6*a^4 - 12*a^2 + 6)*y^3 + (((24*a^8 - 36*a^6 - 36*a^4 + 84*a^2 - 36)*u)*v)*y + ((6*a^12 - 42*a^8 + 48*a^6 + 18*a^4 - 48*a^2 + 18)*u)*v^2 + ((6*a^12 - 42*a^8 + 48*a^6 + 18*a^4 - 48*a^2 + 18)*u^2 + (36*a^14 - 96*a^12 + 36*a^10 + 96*a^8 - 84*a^6 + 12*a^2)*u)*v)*t^3 + ((4*a^12 - 21*a^8 + 12*a^6 + 30*a^4 - 36*a^2 + 11)*y^2 + ((8*a^18 + 6*a^16 - 72*a^14 + 76*a^12 - 8*a^6 - 60*a^4 + 72*a^2 - 22)*u)*v)*t^2 + ((a^24 - 8*a^18 - 6*a^16 + 36*a^14 - 20*a^12 - 21*a^8 + 12*a^6 + 24*a^4 - 24*a^2 + 6)*y)*t]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of labeled digraphs (we take a parameter a = sqrt{2}),\n",
    "# with\n",
    "# the marking variable t for strongly connected components,\n",
    "# the marking variable u for purely source-like components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "# the marking variable y for isolated components\n",
    "dtall = ssdtall * exp_had_prod(inner,g,N)\n",
    "DTALL = [dtall[i] * i.factorial() for i in srange(disp)]\n",
    "DTALL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " t*y,\n",
       " 2*t^2*u*v + t^2*y^2 + t*y,\n",
       " 3*t^3*u^2*v + 3*t^3*u*v^2 + 6*t^3*u*v*y + t^3*y^3 + 12*t^3*u*v + 18*t^2*u*v + 3*t^2*y^2 + 18*t*y,\n",
       " 4*t^4*u^3*v + 42*t^4*u^2*v^2 + 4*t^4*u*v^3 + 12*t^4*u^2*v*y + 12*t^4*u*v^2*y + 12*t^4*u*v*y^2 + t^4*y^4 + 96*t^4*u^2*v + 96*t^4*u*v^2 + 48*t^4*u*v*y + 216*t^4*u*v + 90*t^3*u^2*v + 90*t^3*u*v^2 + 84*t^3*u*v*y + 6*t^3*y^3 + 504*t^3*u*v + 1098*t^2*u*v + 75*t^2*y^2 + 1606*t*y]"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DTALLsub = [sum(sum(sum(sum(dtall[i][0][k][l][m][n].subs(a=sqrt(2)) * i.factorial() * u^n * v^m * y^l * t^k for n in srange(N)) for m in srange(N)) for l in srange(N)) for k in srange(N)) for i in srange(disp)]\n",
    "DTALLsub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1, 0, 0, 0, 0],\n",
       " [t*y, 0, 0, 0, 0],\n",
       " [(2*u*v + y^2)*t^2 + t*y, 0, 0, 0, 0],\n",
       " [(3*u*v^2 + 6*u*v*y + y^3 + 3*(u^2 + 4*u)*v)*t^3 + 3*(6*u*v + y^2)*t^2 + 18*t*y,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [(4*u*v^3 + 12*u*v*y^2 + y^4 + 6*(7*u^2 + 16*u)*v^2 + 4*(u^3 + 24*u^2 + 54*u)*v + 12*(u*v^2 + (u^2 + 4*u)*v)*y)*t^4 + 6*(15*u*v^2 + 14*u*v*y + y^3 + 3*(5*u^2 + 28*u)*v)*t^3 + 3*(366*u*v + 25*y^2)*t^2 + 1606*t*y,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0]]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DTALLMatrix = [[DTALL[i][j].subs(a=sqrt(2)) for j in srange(disp)] for i in srange(disp)]\n",
    "DTALLMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "disp = 3 # How many terms we want to display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(y - v - u + 1)*t,\n",
       " ((-2*a^2*y + 2*a^2*v + 2*a^2*u - 2*a^2)*t)*w,\n",
       " (((-2*a^7 + 6*a^5)*y + (2*a^7 - 6*a^5)*v + (2*a^7 - 6*a^5)*u - 2*a^7 + 6*a^5)*t)*w^2 + ((2*a^3*y^2 + (-2*a^3*v - 2*a^3*u + 2*a^3)*y)*t^2 + (2*a^3*y - 2*a^3*v - 2*a^3*u + 2*a^3)*t)*w]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#dOne = (y-u-v+1) * t * dtall.subs(a^3*z^2*w) * qscd\n",
    "dOne = (y-u-v+1) * t * dtall.subs(z=a^3*z^2*w) * qscd\n",
    "DONE = [dOne[i] * i.factorial() for i in srange(disp)]\n",
    "DONE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[-t*(u + v - y - 1), 0, 0],\n",
       " [0, 4*t*(u + v - y - 1), 0],\n",
       " [0,\n",
       "  4*(sqrt(2)*y^2 - (sqrt(2)*u + sqrt(2)*v - sqrt(2))*y)*t^2 - 4*(sqrt(2)*u + sqrt(2)*v - sqrt(2)*y - sqrt(2))*t,\n",
       "  -8*(sqrt(2)*u + sqrt(2)*v - sqrt(2)*y - sqrt(2))*t]]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DONEMatrix = [[DONE[i][j].subs(a=sqrt(2)) for j in srange(disp)] for i in srange(disp)]\n",
    "DONEMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(y - v - u + 1)*t,\n",
       " ((-2*a^2*y + 2*a^2*v + 2*a^2*u - 2*a^2)*t)*w,\n",
       " (((-a^8 + 3*a^6)*y + (a^8 - 3*a^6)*v + (a^8 - 3*a^6)*u - a^8 + 3*a^6)*t)*w^2 + ((a^4*y^2 + (-a^4*v - a^4*u + a^4)*y)*t^2 + (a^4*y - a^4*v - a^4*u + a^4)*t)*w]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DONEzero = [dOne[i] * a^(i*(i-1)/2) for i in srange(disp)]\n",
    "DONEzero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[-t*(u + v - y - 1), 0, 0],\n",
       " [0, 4*t*(u + v - y - 1), 0],\n",
       " [0, -4*((u + v - 1)*y - y^2)*t^2 - 4*t*(u + v - y - 1), -8*t*(u + v - y - 1)]]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DONEzeroMatrix = [[DONEzero[i][j].subs(a=sqrt(2)) for j in srange(disp)] for i in srange(disp)]\n",
    "DONEzeroMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 0, ((2*a^3*y - 2*a^3*v - 2*a^3*u + 2*a^3)*t)*w]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dTwoZero = ssdtall.subs(z=a^3*z^2*w)\n",
    "DTWOZERO = [dTwoZero[i] * i.factorial() for i in srange(disp)]\n",
    "DTWOZERO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(u - 1)*t,\n",
       " (((2*u - 2)*v)*t^2 + (-2*a^2*u + 2*a^2)*t)*w,\n",
       " (((4/a^2*u - 4/a^2)*v^2 + (((8*a^2 - 8)/a^2)*u + (-8*a^2 + 8)/a^2)*v)*t^3 + ((((-4*a^4 - 8*a^2 + 4)/a^2)*u + (4*a^4 + 8*a^2 - 4)/a^2)*v)*t^2 + ((-2*a^6 + 6*a^4)*u + 2*a^6 - 6*a^4)*t)*w^2 + ((2*a^2*u^2 - 4*a^2*u + 2*a^2)*t^2 + (2*a^2*u - 2*a^2)*t)*w]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qpartu = (u-1) * t * exp_had_prod(PhiTwoTwoOne,ssdtu.subs(z=a^3*z^2*w)*qscd,N)\n",
    "dTwoOne = qpartu * partv.subs(z=2*z*w) * dt.subs(z=2*z*w)\n",
    "DTWOONE = [dTwoOne[i] * i.factorial() for i in srange(disp)]\n",
    "DTWOONE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(v - 1)*t,\n",
       " ((2*u*v - 2*u)*t^2 + (-2*a^2*v + 2*a^2)*t)*w,\n",
       " (((4/a^2*u^2 + ((8*a^2 - 8)/a^2)*u)*v - 4/a^2*u^2 + ((-8*a^2 + 8)/a^2)*u)*t^3 + ((((-4*a^4 - 8*a^2 + 4)/a^2)*u)*v + ((4*a^4 + 8*a^2 - 4)/a^2)*u)*t^2 + ((-2*a^6 + 6*a^4)*v + 2*a^6 - 6*a^4)*t)*w^2 + ((2*a^2*v^2 - 4*a^2*v + 2*a^2)*t^2 + (2*a^2*v - 2*a^2)*t)*w]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qpartv = (v-1) * t * exp_had_prod(PhiTwoTwoOne,ssdtv.subs(z=a^3*z^2*w)*qscd,N)\n",
    "dTwoTwo = qpartv * partu.subs(z=2*z*w) * dt.subs(z=2*z*w)\n",
    "DTWOTWO = [dTwoTwo[i] * i.factorial() for i in srange(disp)]\n",
    "DTWOTWO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[t,\n",
       " ((2*v + 2*u + 2*a^2 - 4)*t^2 - 2*a^2*t)*w,\n",
       " ((4/a^2*v^2 + (8*u + (8*a^4 - 8*a^2 - 8)/a^2)*v + 4/a^2*u^2 + ((8*a^4 - 8*a^2 - 8)/a^2)*u + (6*a^6 - 18*a^4 + 8*a^2 + 8)/a^2)*t^3 + (((-4*a^4 - 8*a^2 + 4)/a^2)*v + ((-4*a^4 - 8*a^2 + 4)/a^2)*u + (2*a^8 - 12*a^6 + 10*a^4 + 16*a^2 - 8)/a^2)*t^2 + (-2*a^6 + 6*a^4)*t)*w^2 + (-2*a^2*t^2 + 2*a^2*t)*w]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dTwoThr = partu.subs(z=2*z*w) * partv.subs(z=2*z*w) * qdt\n",
    "DTWOTHR = [dTwoThr[i] * i.factorial() for i in srange(disp)]\n",
    "DTWOTHR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(v + u - 1)*t,\n",
       " ((4*u*v + 2*a^2 - 4)*t^2 + (-2*a^2*v - 2*a^2*u + 2*a^2)*t)*w,\n",
       " ((4/a*u*v^2 + (4/a*u^2 + ((24*a^2 - 16)/a)*u + 8*a^3 - 16*a)*v + (8*a^3 - 16*a)*u + (6*a^6 - 18*a^4 + 8*a^2 + 8)/a)*t^3 + ((((-8*a^4 - 16*a^2 + 8)/a)*u)*v + (2*a^8 - 12*a^6 + 10*a^4 + 16*a^2 - 8)/a)*t^2 + ((-2*a^7 + 6*a^5)*v + (-2*a^7 + 6*a^5)*u + 2*a^7 - 6*a^5)*t)*w^2 + ((2*a^3*v^2 - 4*a^3*v + 2*a^3*u^2 - 4*a^3*u + 2*a^3)*t^2 + (2*a^3*v + 2*a^3*u - 2*a^3)*t)*w]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dTwo = exp_had_prod(PhiTwoOneTwo,dTwoOne+dTwoTwo+dTwoThr,N)\n",
    "DTWO = [dTwo[i] * i.factorial() for i in srange(disp)]\n",
    "DTWO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[y*t,\n",
       " ((4*u*v + 2*a^2 - 4)*t^2 - 2*a^2*y*t)*w,\n",
       " ((4/a*u*v^2 + (4/a*u^2 + ((24*a^2 - 16)/a)*u + 8*a^3 - 16*a)*v + (8*a^3 - 16*a)*u + (6*a^6 - 18*a^4 + 8*a^2 + 8)/a)*t^3 + ((((-8*a^4 - 16*a^2 + 8)/a)*u)*v + (2*a^8 - 12*a^6 + 10*a^4 + 16*a^2 - 8)/a)*t^2 + ((-2*a^7 + 6*a^5)*y)*t)*w^2 + ((2*a^3*y^2 - 4*a^3*u*v)*t^2 + 2*a^3*y*t)*w]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CGF of labeled digraphs (we take a parameter a = sqrt{2}), with\n",
    "# the marking variable t for strongly connected components,\n",
    "# the marking variable u for purely source-like components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "# the marking variable y for isolated components\n",
    "qdtall = dOne + dTwoZero*dTwo\n",
    "QDTALL = [qdtall[i] * i.factorial() for i in srange(disp)]\n",
    "QDTALL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "disp = 6 # How many terms we want to display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[y*t,\n",
       " ((4*u*v + 2*a^2 - 4)*t^2 - 2*a^2*y*t)*w,\n",
       " ((2*u*v^2 + (2*u^2 + (12*a^2 - 8)*u + 4*a^4 - 8*a^2)*v + (4*a^4 - 8*a^2)*u + 3*a^6 - 9*a^4 + 4*a^2 + 4)*t^3 + (((-4*a^4 - 8*a^2 + 4)*u)*v + a^8 - 6*a^6 + 5*a^4 + 8*a^2 - 4)*t^2 + ((-a^8 + 3*a^6)*y)*t)*w^2 + ((a^4*y^2 - 2*a^4*u*v)*t^2 + a^4*y*t)*w,\n",
       " ((4/3*u*v^3 + ((8*a^4 - 4)*u + 4*a^6 - 8*a^4)*v^2 + (4/3*u^3 + (8*a^4 - 4)*u^2 + (8*a^8 + 16*a^6 - 40*a^4 + 8)*u + 6*a^10 - 10*a^8 - 16*a^6 + 24*a^4)*v + (4*a^6 - 8*a^4)*u^2 + (6*a^10 - 10*a^8 - 16*a^6 + 24*a^4)*u + 4*a^12 - 15*a^10 + 12*a^8 + 25/3*a^6 - 8*a^4 - 8/3)*t^4 + (((-4*a^6 + 4*a^4 - 8*a^2 + 4)*u)*v^2 + ((-4*a^6 + 4*a^4 - 8*a^2 + 4)*u^2 + (-32*a^6 + 8*a^4 + 32*a^2 - 16)*u + 2*a^12 - 8*a^10 - 6*a^8 + 36*a^6 - 16*a^4)*v + (2*a^12 - 8*a^10 - 6*a^8 + 36*a^6 - 16*a^4)*u + 3*a^14 - 16*a^12 + 20*a^10 + 14*a^8 - 33*a^6 + 16*a^4 - 16*a^2 + 8)*t^3 + (((-4/3*a^12 + 4*a^10 + 8*a^6 + 8*a^4 - 16*a^2 + 16/3)*u)*v + 1/3*a^18 - 6*a^14 + 40/3*a^12 - 5*a^10 - 2*a^8 - 22/3*a^6 - 8*a^4 + 16*a^2 - 16/3)*t^2 + ((-1/3*a^18 + 3*a^14 - 4*a^12)*y)*t)*w^3 + (((4*a^6*u*v + 2*a^8 - 4*a^6)*y + ((2*a^8 - 4*a^6)*u)*v^2 + ((2*a^8 - 4*a^6)*u^2 + (-8*a^8 + 4*a^6)*u - 2*a^8 + 4*a^6)*v + (-2*a^8 + 4*a^6)*u - 2*a^10 + 6*a^8 - 4*a^6)*t^3 + (-2*a^8*y^2 + 8*a^8*u*v + 2*a^10 - 4*a^8)*t^2 - 2*a^8*y*t)*w^2,\n",
       " ((2/3*u*v^4 + ((16/3*a^6 - 8/3)*u + 8/3*a^8 - 16/3*a^6)*v^3 + (4*a^8*u^2 + (8*a^12 + 8*a^10 - 12*a^8 - 16*a^6 + 4)*u + 6*a^14 - 10*a^12 - 8*a^10 + 16*a^6)*v^2 + (2/3*u^4 + (16/3*a^6 - 8/3)*u^3 + (8*a^12 + 8*a^10 - 12*a^8 - 16*a^6 + 4)*u^2 + (12*a^16 - 4*a^14 - 64*a^10 + 32*a^8 + 112/3*a^6 - 16/3)*u + 8*a^18 - 18*a^16 - 8*a^14 + 38/3*a^12 + 32*a^10 - 8*a^8 - 64/3*a^6)*v + (8/3*a^8 - 16/3*a^6)*u^3 + (6*a^14 - 10*a^12 - 8*a^10 + 16*a^6)*u^2 + (8*a^18 - 18*a^16 - 8*a^14 + 38/3*a^12 + 32*a^10 - 8*a^8 - 64/3*a^6)*u + 5*a^20 - 22*a^18 + 99/4*a^16 + 9*a^14 - 64/3*a^12 - 17/12*a^8 + 16/3*a^6 + 4/3)*t^5 + (((-8/3*a^8 + 4*a^4 - 8*a^2 + 4)*u)*v^3 + ((-8*a^12 - 16*a^8 + 16*a^6 - 12*a^4 + 24*a^2 - 12)*u + 2*a^16 - 12*a^14 + 18*a^12 - 16*a^10 + 32*a^8 - 16*a^6)*v^2 + ((-8/3*a^8 + 4*a^4 - 8*a^2 + 4)*u^3 + (-8*a^12 - 16*a^8 + 16*a^6 - 12*a^4 + 24*a^2 - 12)*u^2 + (4*a^18 - 8*a^16 - 12*a^14 - 40*a^12 + 48*a^10 + 104*a^8 - 80*a^6 + 24*a^4 - 48*a^2 + 24)*u + 6*a^20 - 22*a^18 - 2*a^16 + 58*a^14 - 16*a^12 + 8*a^10 - 88*a^8 + 48*a^6)*v + (2*a^16 - 12*a^14 + 18*a^12 - 16*a^10 + 32*a^8 - 16*a^6)*u^2 + (6*a^20 - 22*a^18 - 2*a^16 + 58*a^14 - 16*a^12 + 8*a^10 - 88*a^8 + 48*a^6)*u + 6*a^22 - 67/2*a^20 + 52*a^18 + 25/2*a^16 - 251/3*a^14 + 105/2*a^12 - 33*a^10 + 275/6*a^8 - 16*a^6 - 8*a^4 + 16*a^2 - 8)*t^4 + (((-2*a^16 + 6*a^14 - 16/3*a^12 + 16*a^10 - 22*a^8 + 8*a^6 + 20*a^4 - 24*a^2 + 22/3)*u)*v^2 + ((-2*a^16 + 6*a^14 - 16/3*a^12 + 16*a^10 - 22*a^8 + 8*a^6 + 20*a^4 - 24*a^2 + 22/3)*u^2 + (4*a^18 + 8*a^16 - 72*a^14 + 424/3*a^12 - 64*a^10 + 4*a^8 - 80*a^4 + 96*a^2 - 88/3)*u + 2/3*a^24 - 22/3*a^20 + 16/3*a^18 + 10*a^16 + 28*a^14 - 206/3*a^12 - 16*a^10 + 184/3*a^8 - 64/3*a^6)*v + (2/3*a^24 - 22/3*a^20 + 16/3*a^18 + 10*a^16 + 28*a^14 - 206/3*a^12 - 16*a^10 + 184/3*a^8 - 64/3*a^6)*u + a^26 - 7/12*a^24 - 18*a^22 + 281/6*a^20 - 65/3*a^18 - 63/2*a^16 + 7*a^14 + 25/6*a^12 + 67*a^10 - 1003/12*a^8 + 80/3*a^6 + 40*a^4 - 48*a^2 + 44/3)*t^3 + (((8/3*a^20 - 20/3*a^18 - 4*a^16 + 28*a^14 - 128/3*a^12 + 32*a^10 - 116/3*a^8 + 16*a^6 + 32*a^4 - 32*a^2 + 8)*u)*v + 1/12*a^32 - 2*a^26 - 3/2*a^24 + 18*a^22 - 70/3*a^20 + 23/3*a^18 + 25/4*a^16 - 85/3*a^14 + 134/3*a^12 - 34*a^10 + 235/6*a^8 - 16*a^6 - 32*a^4 + 32*a^2 - 8)*t^2 + ((-1/12*a^32 + a^26 + 3/4*a^24 - 6*a^22 + 5*a^20)*y)*t)*w^4 + (((2*a^8*u*v^2 + (2*a^8*u^2 + (12*a^10 - 8*a^8)*u + 4*a^12 - 8*a^10)*v + (4*a^12 - 8*a^10)*u + 3*a^14 - 9*a^12 + 4*a^10 + 4*a^8)*y - 2*a^8*u*v^3 + ((4*a^12 - 4*a^8)*u^2 + (4*a^14 - 8*a^12 - 12*a^10 + 10*a^8)*u - 4*a^12 + 8*a^10)*v^2 + (-2*a^8*u^3 + (4*a^14 - 8*a^12 - 12*a^10 + 10*a^8)*u^2 + (-20*a^14 + 8*a^12 + 28*a^10 - 8*a^8)*u - 4*a^16 + 5*a^14 + 13*a^12 - 12*a^10 - 4*a^8)*v + (-4*a^12 + 8*a^10)*u^2 + (-4*a^16 + 5*a^14 + 13*a^12 - 12*a^10 - 4*a^8)*u - 3*a^18 + 9*a^16 - a^14 - 13*a^12 + 4*a^10 + 4*a^8)*t^4 + ((((-4*a^12 - 8*a^10 + 4*a^8)*u)*v + a^16 - 6*a^14 + 5*a^12 + 8*a^10 - 4*a^8)*y + ((2*a^16 - 8*a^14 + 8*a^12 + 8*a^10 - 4*a^8)*u)*v^2 + ((2*a^16 - 8*a^14 + 8*a^12 + 8*a^10 - 4*a^8)*u^2 + (-8*a^16 + 44*a^14 - 20*a^12 - 8*a^10 + 4*a^8)*u + 3*a^16 - 2*a^14 - 5*a^12 - 8*a^10 + 4*a^8)*v + (3*a^16 - 2*a^14 - 5*a^12 - 8*a^10 + 4*a^8)*u - a^20 + 5*a^18 - a^16 - 18*a^14 + 13*a^12 + 8*a^10 - 4*a^8)*t^3 + ((-a^16 + 3*a^14)*y^2 + ((6*a^16 - 22*a^14 + 4*a^12)*u)*v + a^20 - 2*a^18 - 7*a^16 + 16*a^14 - 4*a^12)*t^2 + ((-a^16 + 3*a^14)*y)*t)*w^3 + ((1/2*a^12*y^3 + (((a^14 - 3*a^12)*u)*v)*y + ((-a^14 + 3/2*a^12)*u)*v^2 + ((-a^14 + 3/2*a^12)*u^2 + a^14*u)*v)*t^3 + ((1/2*a^16 - a^14 + 3/2*a^12)*y^2 + ((-a^16 + 2*a^14 - 3*a^12)*u)*v)*t^2 + ((1/2*a^16 - a^14 + a^12)*y)*t)*w^2,\n",
       " ((4/15*u*v^5 + ((8/3*a^8 - 4/3)*u + 4/3*a^10 - 8/3*a^8)*v^4 + (8/3*a^12*u^2 + (16/3*a^16 + 16/3*a^14 - 8*a^12 - 32/3*a^8 + 8/3)*u + 4*a^18 - 20/3*a^16 - 16/3*a^14 + 16/3*a^12 - 16/3*a^10 + 32/3*a^8)*v^3 + (8/3*a^12*u^3 + (8*a^18 - 16*a^12)*u^2 + (12*a^22 - 4*a^20 - 32*a^16 - 16*a^14 + 104/3*a^12 + 16*a^8 - 8/3)*u + 8*a^24 - 18*a^22 + 4*a^20 - 58/3*a^18 + 36*a^16 + 16*a^14 - 64/3*a^12 + 8*a^10 - 16*a^8)*v^2 + (4/15*u^5 + (8/3*a^8 - 4/3)*u^4 + (16/3*a^16 + 16/3*a^14 - 8*a^12 - 32/3*a^8 + 8/3)*u^3 + (12*a^22 - 4*a^20 - 32*a^16 - 16*a^14 + 104/3*a^12 + 16*a^8 - 8/3)*u^2 + (16*a^26 - 12*a^24 - 48*a^22 + 244/3*a^20 - 128*a^18 + 112*a^16 + 160/3*a^14 - 176/3*a^12 - 24*a^8 + 8/3)*u + 10*a^28 - 28*a^26 - 5/2*a^24 + 50*a^22 - 136/3*a^20 + 224/3*a^18 - 131/2*a^16 - 80/3*a^14 + 80/3*a^12 - 16/3*a^10 + 40/3*a^8)*v + (4/3*a^10 - 8/3*a^8)*u^4 + (4*a^18 - 20/3*a^16 - 16/3*a^14 + 16/3*a^12 - 16/3*a^10 + 32/3*a^8)*u^3 + (8*a^24 - 18*a^22 + 4*a^20 - 58/3*a^18 + 36*a^16 + 16*a^14 - 64/3*a^12 + 8*a^10 - 16*a^8)*u^2 + (10*a^28 - 28*a^26 - 5/2*a^24 + 50*a^22 - 136/3*a^20 + 224/3*a^18 - 131/2*a^16 - 80/3*a^14 + 80/3*a^12 - 16/3*a^10 + 40/3*a^8)*u + 6*a^30 - 30*a^28 + 43*a^26 + 3*a^24 - 81/2*a^22 + 28/3*a^20 + 13/12*a^18 + 41/3*a^16 - 16/3*a^12 + 161/60*a^10 - 8/3*a^8 - 8/15)*t^6 + (((-4/3*a^10 + 8/3*a^4 - 16/3*a^2 + 8/3)*u)*v^4 + ((-16/3*a^16 - 32/3*a^14 + 64/3*a^12 - 80/3*a^10 + 16*a^8 - 32/3*a^4 + 64/3*a^2 - 32/3)*u + 4/3*a^20 - 8*a^18 + 20/3*a^16 + 56/3*a^14 - 112/3*a^12 + 40*a^10 - 16*a^8)*v^3 + ((-8*a^18 + 16*a^16 - 32*a^14 + 16*a^12)*u^2 + (4*a^24 - 16*a^22 + 4*a^20 - 16*a^18 - 32*a^16 + 144*a^14 - 112*a^12 + 88*a^10 - 48*a^8 + 16*a^4 - 32*a^2 + 16)*u + 6*a^26 - 28*a^24 + 36*a^22 - 32*a^20 + 66*a^18 + 8*a^16 - 136*a^14 + 144*a^12 - 120*a^10 + 48*a^8)*v^2 + ((-4/3*a^10 + 8/3*a^4 - 16/3*a^2 + 8/3)*u^4 + (-16/3*a^16 - 32/3*a^14 + 64/3*a^12 - 80/3*a^10 + 16*a^8 - 32/3*a^4 + 64/3*a^2 - 32/3)*u^3 + (4*a^24 - 16*a^22 + 4*a^20 - 16*a^18 - 32*a^16 + 144*a^14 - 112*a^12 + 88*a^10 - 48*a^8 + 16*a^4 - 32*a^2 + 16)*u^2 + (12*a^28 - 24*a^26 - 56*a^24 + 176*a^22 - 284*a^20 + 288*a^18 + 464/3*a^16 - 1376/3*a^14 + 832/3*a^12 - 640/3*a^10 + 112*a^8 - 64/3*a^4 + 128/3*a^2 - 64/3)*u + 12*a^30 - 47*a^28 + 14*a^26 + 129*a^24 - 542/3*a^22 + 611/3*a^20 - 688/3*a^18 - 137/3*a^16 + 776/3*a^14 - 592/3*a^12 + 152*a^10 - 64*a^8)*v + (4/3*a^20 - 8*a^18 + 20/3*a^16 + 56/3*a^14 - 112/3*a^12 + 40*a^10 - 16*a^8)*u^3 + (6*a^26 - 28*a^24 + 36*a^22 - 32*a^20 + 66*a^18 + 8*a^16 - 136*a^14 + 144*a^12 - 120*a^10 + 48*a^8)*u^2 + (12*a^30 - 47*a^28 + 14*a^26 + 129*a^24 - 542/3*a^22 + 611/3*a^20 - 688/3*a^18 - 137/3*a^16 + 776/3*a^14 - 592/3*a^12 + 152*a^10 - 64*a^8)*u + 10*a^32 - 60*a^30 + 110*a^28 - 19/2*a^26 - 166*a^24 + 893/6*a^22 - 73*a^20 + 289/3*a^18 - 34*a^16 - 353/6*a^14 + 209/3*a^12 - 305/6*a^10 + 16*a^8 + 16/3*a^4 - 32/3*a^2 + 16/3)*t^5 + (((-4/3*a^20 + 4*a^18 - 8*a^14 + 56/3*a^12 - 8*a^10 - 12*a^8 + 32*a^4 - 32*a^2 + 28/3)*u)*v^3 + ((-8/3*a^24 + 8*a^22 - 4/3*a^20 + 4*a^18 + 24*a^16 - 40*a^14 + 152/3*a^12 - 72*a^10 + 196/3*a^8 - 96*a^4 + 96*a^2 - 28)*u + 2/3*a^30 - 12*a^26 + 92/3*a^24 - 110/3*a^22 + 172/3*a^20 - 224/3*a^18 - 4*a^16 + 104*a^14 - 464/3*a^12 + 332/3*a^10 - 88/3*a^8)*v^2 + ((-4/3*a^20 + 4*a^18 - 8*a^14 + 56/3*a^12 - 8*a^10 - 12*a^8 + 32*a^4 - 32*a^2 + 28/3)*u^3 + (-8/3*a^24 + 8*a^22 - 4/3*a^20 + 4*a^18 + 24*a^16 - 40*a^14 + 152/3*a^12 - 72*a^10 + 196/3*a^8 - 96*a^4 + 96*a^2 - 28)*u^2 + (4/3*a^32 - 16/3*a^28 + 8*a^26 + 76/3*a^24 - 192*a^22 + 1240/3*a^20 - 208*a^18 - 728/3*a^16 + 336*a^14 - 1264/3*a^12 + 432*a^10 - 656/3*a^8 + 192*a^4 - 192*a^2 + 56)*u + 2*a^34 + 1/6*a^32 - 82/3*a^30 + 97/3*a^28 + 48*a^26 - 95*a^24 + 436/3*a^22 - 1043/3*a^20 + 868/3*a^18 + 469/6*a^16 - 304*a^14 + 416*a^12 - 952/3*a^10 + 88*a^8)*v + (2/3*a^30 - 12*a^26 + 92/3*a^24 - 110/3*a^22 + 172/3*a^20 - 224/3*a^18 - 4*a^16 + 104*a^14 - 464/3*a^12 + 332/3*a^10 - 88/3*a^8)*u^2 + (2*a^34 + 1/6*a^32 - 82/3*a^30 + 97/3*a^28 + 48*a^26 - 95*a^24 + 436/3*a^22 - 1043/3*a^20 + 868/3*a^18 + 469/6*a^16 - 304*a^14 + 416*a^12 - 952/3*a^10 + 88*a^8)*u + 2*a^36 - 3/2*a^34 - 125/3*a^32 + 755/6*a^30 - 250/3*a^28 - 1243/12*a^26 + 377/3*a^24 - 313/6*a^22 + 551/3*a^20 - 1451/6*a^18 + 139/3*a^16 + 162*a^14 - 242*a^12 + 1703/12*a^10 - 16/3*a^8 - 64*a^4 + 64*a^2 - 56/3)*t^4 + (((-2/3*a^30 + 6*a^26 - 32/3*a^24 + 44/3*a^22 - 28*a^20 + 88/3*a^18 - 64/3*a^16 + 8/3*a^14 + 8*a^12 + 100/3*a^10 - 68*a^8 + 232/3*a^4 - 176/3*a^2 + 40/3)*u)*v^2 + ((-2/3*a^30 + 6*a^26 - 32/3*a^24 + 44/3*a^22 - 28*a^20 + 88/3*a^18 - 64/3*a^16 + 8/3*a^14 + 8*a^12 + 100/3*a^10 - 68*a^8 + 232/3*a^4 - 176/3*a^2 + 40/3)*u^2 + (4*a^32 + 8/3*a^30 - 4*a^28 - 72*a^26 + 104/3*a^24 + 760/3*a^22 - 288*a^20 - 136/3*a^18 + 376/3*a^16 - 224/3*a^14 + 192*a^12 - 976/3*a^10 + 320*a^8 - 928/3*a^4 + 704/3*a^2 - 160/3)*u + 1/6*a^40 - 2*a^34 - 17/3*a^32 + 50/3*a^30 - 8*a^28 + 20*a^26 + 179/6*a^24 - 252*a^22 + 956/3*a^20 - 308/3*a^18 - 7*a^16 + 224/3*a^14 - 544/3*a^12 + 136*a^10 - 32*a^8)*v + (1/6*a^40 - 2*a^34 - 17/3*a^32 + 50/3*a^30 - 8*a^28 + 20*a^26 + 179/6*a^24 - 252*a^22 + 956/3*a^20 - 308/3*a^18 - 7*a^16 + 224/3*a^14 - 544/3*a^12 + 136*a^10 - 32*a^8)*u + 1/4*a^42 - 1/3*a^40 + 1/2*a^38 - 6*a^36 - 17/4*a^34 + 202/3*a^32 - 205/2*a^30 + 77/3*a^28 + 451/12*a^26 - 248/3*a^24 + 1415/6*a^22 - 925/3*a^20 + 2395/12*a^18 - 206/3*a^16 - 509/6*a^14 + 475/3*a^12 - 85/6*a^10 - 120*a^8 + 464/3*a^4 - 352/3*a^2 + 80/3)*t^3 + (((1/5*a^40 - 7/3*a^32 - 16/3*a^30 + 20/3*a^28 + 28*a^26 - 64/3*a^24 - 104/3*a^22 + 32/3*a^20 + 152/3*a^18 - 128/3*a^16 + 64/3*a^14 - 64/3*a^12 + 80*a^10 - 112*a^8 + 96*a^4 - 64*a^2 + 64/5)*u)*v + 1/60*a^50 - 1/2*a^42 - 1/5*a^40 - a^38 + 6*a^36 + 35/4*a^34 - 113/3*a^32 + 110/3*a^30 - 19/3*a^28 - 63/2*a^26 + 24*a^24 + 28*a^22 - 11/3*a^20 - 331/6*a^18 + 128/3*a^16 - 55/3*a^14 + 58/3*a^12 - 398/5*a^10 + 112*a^8 - 96*a^4 + 64*a^2 - 64/5)*t^2 + ((-1/60*a^50 + 1/4*a^42 + 1/2*a^38 - 2*a^36 - 3*a^34 + 10*a^32 - 6*a^30)*y)*t)*w^5 + (((4/3*a^10*u*v^3 + ((8*a^14 - 4*a^10)*u + 4*a^16 - 8*a^14)*v^2 + (4/3*a^10*u^3 + (8*a^14 - 4*a^10)*u^2 + (8*a^18 + 16*a^16 - 40*a^14 + 8*a^10)*u + 6*a^20 - 10*a^18 - 16*a^16 + 24*a^14)*v + (4*a^16 - 8*a^14)*u^2 + (6*a^20 - 10*a^18 - 16*a^16 + 24*a^14)*u + 4*a^22 - 15*a^20 + 12*a^18 + 25/3*a^16 - 8*a^14 - 8/3*a^10)*y - 4/3*a^10*u*v^4 + ((4/3*a^16 - 4/3*a^10)*u^2 + (-8/3*a^16 - 8*a^14 + 16/3*a^10)*u - 4*a^16 + 8*a^14)*v^3 + ((4/3*a^16 - 4/3*a^10)*u^3 + (8*a^20 - 8*a^16 - 16*a^14 + 8*a^10)*u^2 + (8*a^22 - 24*a^20 - 8*a^18 - 8*a^16 + 56*a^14 - 12*a^10)*u - 4*a^22 + 2*a^20 + 10*a^18 + 20*a^16 - 32*a^14)*v^2 + (-4/3*a^10*u^4 + (-8/3*a^16 - 8*a^14 + 16/3*a^10)*u^3 + (8*a^22 - 24*a^20 - 8*a^18 - 8*a^16 + 56*a^14 - 12*a^10)*u^2 + (-8*a^24 - 32*a^22 + 52*a^20 + 28*a^18 + 32*a^16 - 88*a^14 + 8*a^10)*u - 6*a^26 + 10*a^24 + 12*a^22 - 3*a^20 - 22*a^18 - 73/3*a^16 + 32*a^14 + 8/3*a^10)*v + (-4*a^16 + 8*a^14)*u^3 + (-4*a^22 + 2*a^20 + 10*a^18 + 20*a^16 - 32*a^14)*u^2 + (-6*a^26 + 10*a^24 + 12*a^22 - 3*a^20 - 22*a^18 - 73/3*a^16 + 32*a^14 + 8/3*a^10)*u - 4*a^28 + 15*a^26 - 12*a^24 - 13/3*a^22 - 7*a^20 + 12*a^18 + 11*a^16 - 8*a^14 - 8/3*a^10)*t^5 + ((((-4*a^16 + 4*a^14 - 8*a^12 + 4*a^10)*u)*v^2 + ((-4*a^16 + 4*a^14 - 8*a^12 + 4*a^10)*u^2 + (-32*a^16 + 8*a^14 + 32*a^12 - 16*a^10)*u + 2*a^22 - 8*a^20 - 6*a^18 + 36*a^16 - 16*a^14)*v + (2*a^22 - 8*a^20 - 6*a^18 + 36*a^16 - 16*a^14)*u + 3*a^24 - 16*a^22 + 20*a^20 + 14*a^18 - 33*a^16 + 16*a^14 - 16*a^12 + 8*a^10)*y + ((16/3*a^16 - 4*a^14 + 8*a^12 - 4*a^10)*u)*v^3 + ((-16*a^18 + 16*a^16 - 8*a^14 + 16*a^12 - 8*a^10)*u^2 + (8*a^24 - 24*a^22 + 16*a^20 + 32*a^18 + 8*a^16 - 4*a^14 - 40*a^12 + 20*a^10)*u + 2*a^22 + 6*a^18 - 36*a^16 + 16*a^14)*v^2 + ((16/3*a^16 - 4*a^14 + 8*a^12 - 4*a^10)*u^3 + (8*a^24 - 24*a^22 + 16*a^20 + 32*a^18 + 8*a^16 - 4*a^14 - 40*a^12 + 20*a^10)*u^2 + (-32*a^24 + 132*a^22 - 64*a^20 - 52*a^18 - 64*a^16 + 40*a^14 + 32*a^12 - 16*a^10)*u - 2*a^28 + 6*a^26 + 17*a^24 - 50*a^22 + 12*a^20 - 20*a^18 + 69*a^16 - 32*a^14 + 16*a^12 - 8*a^10)*v + (2*a^22 + 6*a^18 - 36*a^16 + 16*a^14)*u^2 + (-2*a^28 + 6*a^26 + 17*a^24 - 50*a^22 + 12*a^20 - 20*a^18 + 69*a^16 - 32*a^14 + 16*a^12 - 8*a^10)*u - 3*a^30 + 14*a^28 - 11*a^26 - 25*a^24 + 76/3*a^22 + 4*a^20 + 30*a^18 - 131/3*a^16 + 16*a^14 - 16*a^12 + 8*a^10)*t^4 + ((((-4/3*a^22 + 4*a^20 + 8*a^16 + 8*a^14 - 16*a^12 + 16/3*a^10)*u)*v + 1/3*a^28 - 6*a^24 + 40/3*a^22 - 5*a^20 - 2*a^18 - 22/3*a^16 - 8*a^14 + 16*a^12 - 16/3*a^10)*y + ((4/3*a^28 - 14*a^24 + 70/3*a^22 - 4*a^20 - 16*a^18 - 4/3*a^16 - 8*a^14 + 16*a^12 - 16/3*a^10)*u)*v^2 + ((4/3*a^28 - 14*a^24 + 70/3*a^22 - 4*a^20 - 16*a^18 - 4/3*a^16 - 8*a^14 + 16*a^12 - 16/3*a^10)*u^2 + (-16/3*a^28 + 80*a^24 - 484/3*a^22 + 28*a^20 + 64*a^18 - 56/3*a^16 + 8*a^14 - 16*a^12 + 16/3*a^10)*u + 5/3*a^28 - 24*a^24 + 116/3*a^22 - 11*a^20 + 2*a^18 + 22/3*a^16 + 8*a^14 - 16*a^12 + 16/3*a^10)*v + (5/3*a^28 - 24*a^24 + 116/3*a^22 - 11*a^20 + 2*a^18 + 22/3*a^16 + 8*a^14 - 16*a^12 + 16/3*a^10)*u - 1/3*a^34 + 5*a^30 - 7*a^28 - 7*a^26 + 4*a^24 + 71/3*a^22 + 3*a^20 - 34*a^18 + 6*a^16 - 8*a^14 + 16*a^12 - 16/3*a^10)*t^3 + ((-1/3*a^28 + 3*a^24 - 4*a^22)*y^2 + ((10/3*a^28 - 34*a^24 + 52*a^22 - 16*a^18 + 16/3*a^16)*u)*v + 1/3*a^34 - 2*a^30 - 8/3*a^28 + 3*a^26 + 30*a^24 - 130/3*a^22 + 16*a^18 - 16/3*a^16)*t^2 + ((-1/3*a^28 + 3*a^24 - 4*a^22)*y)*t)*w^4 + (((2*a^16*u*v + a^18 - 2*a^16)*y^2 + (((2*a^18 - 4*a^16)*u)*v^2 + ((2*a^18 - 4*a^16)*u^2 + (-8*a^18 + 4*a^16)*u - 2*a^18 + 4*a^16)*v + (-2*a^18 + 4*a^16)*u - 2*a^20 + 6*a^18 - 4*a^16)*y + ((a^20 - 2*a^18 + 2*a^16)*u)*v^3 + ((-4*a^18 + 4*a^16)*u^2 + (-3*a^20 + 10*a^18 - 4*a^16)*u + a^18 - 2*a^16)*v^2 + ((a^20 - 2*a^18 + 2*a^16)*u^3 + (-3*a^20 + 10*a^18 - 4*a^16)*u^2 + (6*a^20 - 6*a^18 - 2*a^16)*u + 2*a^20 - 6*a^18 + 4*a^16)*v + (a^18 - 2*a^16)*u^2 + (2*a^20 - 6*a^18 + 4*a^16)*u + a^22 - 4*a^20 + 5*a^18 - 2*a^16)*t^4 + (-a^18*y^3 + (((6*a^18 + 2*a^16)*u)*v + a^22 - 2*a^20 + a^18 - 2*a^16)*y + ((a^24 - 2*a^22 + 3*a^20 - 3*a^18 - 2*a^16)*u)*v^2 + ((a^24 - 2*a^22 + 3*a^20 - 3*a^18 - 2*a^16)*u^2 + (-4*a^24 + 8*a^22 - 12*a^20 + 2*a^16)*u - a^22 + 2*a^20 - a^18 + 2*a^16)*v + (-a^22 + 2*a^20 - a^18 + 2*a^16)*u - a^26 + 4*a^24 - 6*a^22 + 4*a^20 + a^18 - 2*a^16)*t^3 + ((-a^22 + 2*a^20 - 3*a^18)*y^2 + ((2*a^24 - 2*a^22 + 6*a^18)*u)*v + a^26 - 4*a^24 + 6*a^22 - 4*a^20)*t^2 + ((-a^22 + 2*a^20 - 2*a^18)*y)*t)*w^3]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs (we take a parameter a = sqrt{2}),\n",
    "# with\n",
    "# the marking variable t for strongly connected components,\n",
    "# the marking variable u for purely source-like components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "# the marking variable y for isolated components\n",
    "QDTALLzero = [qdtall[i] * a^(i*(i-1)/2) for i in srange(disp)]\n",
    "QDTALLzero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[t*y, 0, 0, 0, 0, 0],\n",
       " [0, 4*t^2*u*v - 4*t*y, 0, 0, 0, 0],\n",
       " [0,\n",
       "  -4*(2*u*v - y^2)*t^2 + 4*t*y,\n",
       "  2*(u*v^2 + (u^2 + 8*u)*v)*t^3 - 28*t^2*u*v + 8*t*y,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  32*(u*v*y - 3*u*v)*t^3 + 32*(4*u*v - y^2)*t^2 - 32*t*y,\n",
       "  4/3*(u*v^3 + 21*u*v^2 + (u^3 + 21*u^2 + 78*u)*v)*t^4 - 4*(7*u*v^2 + (7*u^2 + 44*u)*v)*t^3 + 112*t^2*u*v - 128/3*t*y,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  -32*(u*v^2 + 2*u*v*y - y^3 + (u^2 - 4*u)*v)*t^3 - 96*(2*u*v - y^2)*t^2 + 64*t*y,\n",
       "  -32*(u*v^3 - (6*u^2 - 7*u)*v^2 + (u^3 + 7*u^2 + 40*u)*v - (u*v^2 + (u^2 + 8*u)*v)*y)*t^4 + 64*(3*u*v^2 - 7*u*v*y + 3*(u^2 + 11*u)*v)*t^3 - 128*(8*u*v - y^2)*t^2 + 128*t*y,\n",
       "  2/3*(u*v^4 + 60*u*v^3 + 6*(16*u^2 + 113*u)*v^2 + (u^4 + 60*u^3 + 678*u^2 + 1976*u)*v)*t^5 - 4/3*(29*u*v^3 + 489*u*v^2 + (29*u^3 + 489*u^2 + 1134*u)*v)*t^4 + 2*(89*u*v^2 + (89*u^2 + 892*u)*v)*t^3 - 248*t^2*u*v - 4096/3*t*y,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  512*(u*v^3 + u*v*y^2 - 2*(u^2 - u)*v^2 - 6*u*v*y + (u^3 + 2*u^2 + 5*u)*v)*t^4 + 512*(2*u*v^2 + 7*u*v*y - y^3 + (2*u^2 - 23*u)*v)*t^3 + 512*(14*u*v - 3*y^2)*t^2 - 1024*t*y,\n",
       "  -128/3*(u*v^4 - (7*u^2 - 36*u)*v^3 - (7*u^3 + 102*u^2 - 177*u)*v^2 + (u^4 + 36*u^3 + 177*u^2 + 786*u)*v - (u*v^3 + 21*u*v^2 + (u^3 + 21*u^2 + 78*u)*v)*y)*t^5 + 128/3*(29*u*v^3 - 3*(34*u^2 - 125*u)*v^2 + (29*u^3 + 375*u^2 + 876*u)*v - 3*(7*u*v^2 + (7*u^2 + 44*u)*v)*y)*t^4 - 512*(u*v^2 - 7*u*v*y + (u^2 + 65*u)*v)*t^3 + 4096/3*(11*u*v - y^2)*t^2 - 4096/3*t*y,\n",
       "  4/15*(u*v^5 + 155*u*v^4 + 10*(64*u^2 + 513*u)*v^3 + 10*(64*u^3 + 1152*u^2 + 4767*u)*v^2 + (u^5 + 155*u^4 + 5130*u^3 + 47670*u^2 + 122730*u)*v)*t^6 - 8/3*(15*u*v^4 + 740*u*v^3 + 6*(192*u^2 + 959*u)*v^2 + (15*u^4 + 740*u^3 + 5754*u^2 + 7016*u)*v)*t^5 + 4/3*(359*u*v^3 + 6699*u*v^2 + (359*u^3 + 6699*u^2 + 43146*u)*v)*t^4 - 8*(553*u*v^2 + (553*u^2 - 8196*u)*v)*t^3 + 271808/3*t^2*u*v - 3473408/15*t*y]]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs (Section A.10),\n",
    "# with\n",
    "# the marking variable t for strongly connected components,\n",
    "# the marking variable u for purely source-like components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "# the marking variable y for isolated components\n",
    "QDTALLzeroMatrix = [[QDTALLzero[i][j].subs(a=sqrt(2)) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTALLzeroMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[y, 0, 0, 0, 0, 0],\n",
       " [0, -4*y, 0, 0, 0, 0],\n",
       " [0, 4*y, 8*y, 0, 0, 0],\n",
       " [0, 0, -32*y, -128/3*y, 0, 0],\n",
       " [0, 0, 64*y, 128*y, -4096/3*y, 0],\n",
       " [0, 0, 0, -1024*y, -4096/3*y, -3473408/15*y]]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with one strongly connected component (Table 4)\n",
    "# the marking variable u for purely source-like components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "# the marking variable y for isolated components\n",
    "QDTALLzeroMatrixPrimT = [[QDTALLzeroMatrix[i][j].diff(t).subs(t=0) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTALLzeroMatrixPrimT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 0, 0, 0, 0],\n",
       " [0, 4*u*v, 0, 0, 0, 0],\n",
       " [0, -8*u*v + 4*y^2, -28*u*v, 0, 0, 0],\n",
       " [0, 0, 128*u*v - 32*y^2, 112*u*v, 0, 0],\n",
       " [0, 0, -192*u*v + 96*y^2, -1024*u*v + 128*y^2, -248*u*v, 0],\n",
       " [0, 0, 0, 7168*u*v - 1536*y^2, 45056/3*u*v - 4096/3*y^2, 271808/3*u*v]]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with two strongly connected components (Table 19)\n",
    "# the marking variable u for purely source-like components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "# the marking variable y for isolated components\n",
    "QDTALLzeroMatrixSecT = [[QDTALLzeroMatrix[i][j].diff(t,2).subs(t=0) / 2 for j in srange(disp)] for i in srange(disp)]\n",
    "QDTALLzeroMatrixSecT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 2*u*v^2 + 2*(u^2 + 8*u)*v, 0, 0, 0],\n",
       " [0, 0, 32*u*v*y - 96*u*v, -28*u*v^2 - 4*(7*u^2 + 44*u)*v, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  -32*u*v^2 - 64*u*v*y + 32*y^3 - 32*(u^2 - 4*u)*v,\n",
       "  192*u*v^2 - 448*u*v*y + 192*(u^2 + 11*u)*v,\n",
       "  178*u*v^2 + 2*(89*u^2 + 892*u)*v,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  1024*u*v^2 + 3584*u*v*y - 512*y^3 + 512*(2*u^2 - 23*u)*v,\n",
       "  -512*u*v^2 + 3584*u*v*y - 512*(u^2 + 65*u)*v,\n",
       "  -4424*u*v^2 - 8*(553*u^2 - 8196*u)*v]]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with three strongly connected components (Table 20)\n",
    "# the marking variable u for purely source-like components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "# the marking variable y for isolated components\n",
    "QDTALLzeroMatrixThirT = [[QDTALLzeroMatrix[i][j].diff(t,3).subs(t=0) / 6 for j in srange(disp)] for i in srange(disp)]\n",
    "QDTALLzeroMatrixThirT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[t*y, 0, 0, 0, 0, 0],\n",
       " [0, 4*t^2*y - 4*t*y, 0, 0, 0, 0],\n",
       " [0, 4*(y^2 - 2*y)*t^2 + 4*t*y, 2*(y^2 + 9*y)*t^3 - 28*t^2*y + 8*t*y, 0, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  32*(y^2 - 3*y)*t^3 - 32*(y^2 - 4*y)*t^2 - 32*t*y,\n",
       "  4/3*(y^3 + 21*y^2 + 100*y)*t^4 - 4*(7*y^2 + 51*y)*t^3 + 112*t^2*y - 128/3*t*y,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  32*(y^3 - 3*y^2 + 3*y)*t^3 + 96*(y^2 - 2*y)*t^2 + 64*t*y,\n",
       "  -32*(y^3 - (y^2 + 9*y)*y + y^2 + 48*y)*t^4 - 256*(y^2 - 9*y)*t^3 + 128*(y^2 - 8*y)*t^2 + 128*t*y,\n",
       "  2/3*(y^4 + 60*y^3 + 774*y^2 + 2715*y)*t^5 - 4/3*(29*y^3 + 489*y^2 + 1652*y)*t^4 + 2*(89*y^2 + 981*y)*t^3 - 248*t^2*y - 4096/3*t*y,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  1024*(y^3 - 3*y^2 + 4*y)*t^4 - 512*(y^3 - 9*y^2 + 21*y)*t^3 - 512*(3*y^2 - 14*y)*t^2 - 1024*t*y,\n",
       "  -128/3*(y^4 + 29*y^3 - (y^3 + 21*y^2 + 100*y)*y + 68*y^2 + 1000*y)*t^5 + 128/3*(29*y^3 - 3*(7*y^2 + 51*y)*y + 273*y^2 + 1280*y)*t^4 + 3072*(y^2 - 11*y)*t^3 - 4096/3*(y^2 - 11*y)*t^2 - 4096/3*t*y,\n",
       "  4/15*(y^5 + 155*y^4 + 5770*y^3 + 59830*y^2 + 175686*y)*t^6 - 8/3*(15*y^4 + 740*y^3 + 6906*y^2 + 13525*y)*t^5 + 4/3*(359*y^3 + 6699*y^2 + 50204*y)*t^4 - 8*(553*y^2 - 7643*y)*t^3 + 271808/3*t^2*y - 3473408/15*t*y]]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs (Section A.9 with s=y)\n",
    "# the marking variable t for all strongly connected components,\n",
    "# the marking variable y for source-like components,\n",
    "QDTSMatrix = [[QDTALLzeroMatrix[i][j].subs(u=y,v=1) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTSMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[t*y, 0, 0, 0, 0, 0],\n",
       " [0, 4*t^2*y^2 - 4*t*y, 0, 0, 0, 0],\n",
       " [0,\n",
       "  -4*t^2*y^2 + 4*t*y,\n",
       "  2*(y^3 + (y^2 + 8*y)*y)*t^3 - 28*t^2*y^2 + 8*t*y,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  32*(y^3 - 3*y^2)*t^3 + 96*t^2*y^2 - 32*t*y,\n",
       "  4/3*(y^4 + 21*y^3 + (y^3 + 21*y^2 + 78*y)*y)*t^4 - 4*(7*y^3 + (7*y^2 + 44*y)*y)*t^3 + 112*t^2*y^2 - 128/3*t*y,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  -32*(2*y^3 + (y^2 - 4*y)*y)*t^3 - 96*t^2*y^2 + 64*t*y,\n",
       "  -32*(y^4 - (6*y^2 - 7*y)*y^2 - (y^3 + (y^2 + 8*y)*y)*y + (y^3 + 7*y^2 + 40*y)*y)*t^4 - 64*(4*y^3 - 3*(y^2 + 11*y)*y)*t^3 - 896*t^2*y^2 + 128*t*y,\n",
       "  2/3*(y^5 + 60*y^4 + 6*(16*y^2 + 113*y)*y^2 + (y^4 + 60*y^3 + 678*y^2 + 1976*y)*y)*t^5 - 4/3*(29*y^4 + 489*y^3 + (29*y^3 + 489*y^2 + 1134*y)*y)*t^4 + 2*(89*y^3 + (89*y^2 + 892*y)*y)*t^3 - 248*t^2*y^2 - 4096/3*t*y,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  512*(2*y^4 - 2*(y^2 - y)*y^2 - 6*y^3 + (y^3 + 2*y^2 + 5*y)*y)*t^4 + 512*(8*y^3 + (2*y^2 - 23*y)*y)*t^3 + 5632*t^2*y^2 - 1024*t*y,\n",
       "  -128/3*(y^5 - (7*y^2 - 36*y)*y^3 - (7*y^3 + 102*y^2 - 177*y)*y^2 + (y^4 + 36*y^3 + 177*y^2 + 786*y)*y - (y^4 + 21*y^3 + (y^3 + 21*y^2 + 78*y)*y)*y)*t^5 + 128/3*(29*y^4 - 3*(34*y^2 - 125*y)*y^2 + (29*y^3 + 375*y^2 + 876*y)*y - 3*(7*y^3 + (7*y^2 + 44*y)*y)*y)*t^4 + 512*(6*y^3 - (y^2 + 65*y)*y)*t^3 + 40960/3*t^2*y^2 - 4096/3*t*y,\n",
       "  4/15*(y^6 + 155*y^5 + 10*(64*y^2 + 513*y)*y^3 + 10*(64*y^3 + 1152*y^2 + 4767*y)*y^2 + (y^5 + 155*y^4 + 5130*y^3 + 47670*y^2 + 122730*y)*y)*t^6 - 8/3*(15*y^5 + 740*y^4 + 6*(192*y^2 + 959*y)*y^2 + (15*y^4 + 740*y^3 + 5754*y^2 + 7016*y)*y)*t^5 + 4/3*(359*y^4 + 6699*y^3 + (359*y^3 + 6699*y^2 + 43146*y)*y)*t^4 - 8*(553*y^3 + (553*y^2 - 8196*y)*y)*t^3 + 271808/3*t^2*y^2 - 3473408/15*t*y]]"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs\n",
    "# the marking variable t for all strongly connected components,\n",
    "# the marking variable y for components that are either source-like or sink-like\n",
    "QDTOMatrix = [[QDTALLzeroMatrix[i][j].subs(u=y,v=y) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTOMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 0, 0, 0, 0],\n",
       " [0, 4*t^2*v, 0, 0, 0, 0],\n",
       " [0, -8*t^2*v, 2*(v^2 + 8*v)*t^3 - 28*t^2*v, 0, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  32*(v*y - 3*v)*t^3 + 128*t^2*v,\n",
       "  4/3*(v^3 + 21*v^2 + 78*v)*t^4 - 4*(7*v^2 + 44*v)*t^3 + 112*t^2*v,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  -32*(v^2 + 2*v*y - 4*v)*t^3 - 192*t^2*v,\n",
       "  -32*(v^3 + 7*v^2 - (v^2 + 8*v)*y + 40*v)*t^4 + 64*(3*v^2 - 7*v*y + 33*v)*t^3 - 1024*t^2*v,\n",
       "  2/3*(v^4 + 60*v^3 + 678*v^2 + 1976*v)*t^5 - 4/3*(29*v^3 + 489*v^2 + 1134*v)*t^4 + 2*(89*v^2 + 892*v)*t^3 - 248*t^2*v,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  512*(v^3 + v*y^2 + 2*v^2 - 6*v*y + 5*v)*t^4 + 512*(2*v^2 + 7*v*y - 23*v)*t^3 + 7168*t^2*v,\n",
       "  -128/3*(v^4 + 36*v^3 + 177*v^2 - (v^3 + 21*v^2 + 78*v)*y + 786*v)*t^5 + 128/3*(29*v^3 + 375*v^2 - 3*(7*v^2 + 44*v)*y + 876*v)*t^4 - 512*(v^2 - 7*v*y + 65*v)*t^3 + 45056/3*t^2*v,\n",
       "  4/15*(v^5 + 155*v^4 + 5130*v^3 + 47670*v^2 + 122730*v)*t^6 - 8/3*(15*v^4 + 740*v^3 + 5754*v^2 + 7016*v)*t^5 + 4/3*(359*v^3 + 6699*v^2 + 43146*v)*t^4 - 8*(553*v^2 - 8196*v)*t^3 + 271808/3*t^2*v]]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with one purely source-like strongly connected component\n",
    "# the marking variable t for all strongly connected components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "# the marking variable y for isolated components\n",
    "QDTALLzeroMatrixPrimU = [[QDTALLzeroMatrix[i][j].diff(u).subs(u=0) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTALLzeroMatrixPrimU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 0, 0, 0, 0],\n",
       " [0, 4*t^2, 0, 0, 0, 0],\n",
       " [0, -8*t^2, 18*t^3 - 28*t^2, 0, 0, 0],\n",
       " [0, 0, -64*t^3 + 128*t^2, 400/3*t^4 - 204*t^3 + 112*t^2, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  32*t^3 - 192*t^2,\n",
       "  -1248*t^4 + 1856*t^3 - 1024*t^2,\n",
       "  1810*t^5 - 6608/3*t^4 + 1962*t^3 - 248*t^2,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  1536*t^4 - 7168*t^3 + 7168*t^2,\n",
       "  -38400*t^5 + 144256/3*t^4 - 30208*t^3 + 45056/3*t^2,\n",
       "  234248/5*t^6 - 108200/3*t^5 + 200816/3*t^4 + 61144*t^3 + 271808/3*t^2]]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with one purely source-like strongly connected component (Table 21)\n",
    "# the marking variable t for all strongly connected components,\n",
    "QDTALLzeroMatrixPrimU = [[QDTALLzeroMatrix[i][j].diff(u).subs(u=0,v=1,y=1) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTALLzeroMatrixPrimU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[t, 0, 0, 0, 0, 0],\n",
       " [0, -4*t, 0, 0, 0, 0],\n",
       " [0, 4*t, 8*t, 0, 0, 0],\n",
       " [0, 0, 32*t^3*u*v - 32*t, -128/3*t, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  -64*t^3*u*v + 64*t,\n",
       "  32*(u*v^2 + (u^2 + 8*u)*v)*t^4 - 448*t^3*u*v + 128*t,\n",
       "  -4096/3*t,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  -3072*t^4*u*v + 3584*t^3*u*v - 1024*t,\n",
       "  128/3*(u*v^3 + 21*u*v^2 + (u^3 + 21*u^2 + 78*u)*v)*t^5 - 128*(7*u*v^2 + (7*u^2 + 44*u)*v)*t^4 + 3584*t^3*u*v - 4096/3*t,\n",
       "  -3473408/15*t]]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with one isolated strongly connected component\n",
    "# the marking variable t for all strongly connected components,\n",
    "# the marking variable v for purely source-like components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "QDTALLzeroMatrixPrimY = [[QDTALLzeroMatrix[i][j].diff(y).subs(y=0) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTALLzeroMatrixPrimY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[t, 0, 0, 0, 0, 0],\n",
       " [0, -4*t, 0, 0, 0, 0],\n",
       " [0, 4*t, 8*t, 0, 0, 0],\n",
       " [0, 0, 32*t^3 - 32*t, -128/3*t, 0, 0],\n",
       " [0, 0, -64*t^3 + 64*t, 320*t^4 - 448*t^3 + 128*t, -4096/3*t, 0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  -3072*t^4 + 3584*t^3 - 1024*t,\n",
       "  15616/3*t^5 - 7424*t^4 + 3584*t^3 - 4096/3*t,\n",
       "  -3473408/15*t]]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with one isolated strongly connected component (Table 22)\n",
    "# the marking variable t for all strongly connected components,\n",
    "QDTALLzeroMatrixPrimY = [[QDTALLzeroMatrix[i][j].diff(y).subs(y=0,u=1,v=1) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTALLzeroMatrixPrimY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "256"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd([15616, 7424, 3584, 4096])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[61, 29, 14, 16]"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[15616/256, 7424/256, 3584/256, 4096/256]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[y*t,\n",
       " ((4*u*v + 2*a^2 - 4)*t^2 - 2*a^2*y*t)*w,\n",
       " ((2*u*v^2 + (2*u^2 + (12*a^2 - 8)*u + 4*a^4 - 8*a^2)*v + (4*a^4 - 8*a^2)*u + 3*a^6 - 9*a^4 + 4*a^2 + 4)*t^3 + (((-4*a^4 - 8*a^2 + 4)*u)*v + a^8 - 6*a^6 + 5*a^4 + 8*a^2 - 4)*t^2 + ((-a^8 + 3*a^6)*y)*t)*w^2 + ((a^4*y^2 - 2*a^4*u*v)*t^2 + a^4*y*t)*w,\n",
       " ((4/3*u*v^3 + ((8*a^4 - 4)*u + 4*a^6 - 8*a^4)*v^2 + (4/3*u^3 + (8*a^4 - 4)*u^2 + (8*a^8 + 16*a^6 - 40*a^4 + 8)*u + 6*a^10 - 10*a^8 - 16*a^6 + 24*a^4)*v + (4*a^6 - 8*a^4)*u^2 + (6*a^10 - 10*a^8 - 16*a^6 + 24*a^4)*u + 4*a^12 - 15*a^10 + 12*a^8 + 25/3*a^6 - 8*a^4 - 8/3)*t^4 + (((-4*a^6 + 4*a^4 - 8*a^2 + 4)*u)*v^2 + ((-4*a^6 + 4*a^4 - 8*a^2 + 4)*u^2 + (-32*a^6 + 8*a^4 + 32*a^2 - 16)*u + 2*a^12 - 8*a^10 - 6*a^8 + 36*a^6 - 16*a^4)*v + (2*a^12 - 8*a^10 - 6*a^8 + 36*a^6 - 16*a^4)*u + 3*a^14 - 16*a^12 + 20*a^10 + 14*a^8 - 33*a^6 + 16*a^4 - 16*a^2 + 8)*t^3 + (((-4/3*a^12 + 4*a^10 + 8*a^6 + 8*a^4 - 16*a^2 + 16/3)*u)*v + 1/3*a^18 - 6*a^14 + 40/3*a^12 - 5*a^10 - 2*a^8 - 22/3*a^6 - 8*a^4 + 16*a^2 - 16/3)*t^2 + ((-1/3*a^18 + 3*a^14 - 4*a^12)*y)*t)*w^3 + (((4*a^6*u*v + 2*a^8 - 4*a^6)*y + ((2*a^8 - 4*a^6)*u)*v^2 + ((2*a^8 - 4*a^6)*u^2 + (-8*a^8 + 4*a^6)*u - 2*a^8 + 4*a^6)*v + (-2*a^8 + 4*a^6)*u - 2*a^10 + 6*a^8 - 4*a^6)*t^3 + (-2*a^8*y^2 + 8*a^8*u*v + 2*a^10 - 4*a^8)*t^2 - 2*a^8*y*t)*w^2,\n",
       " ((2/3*u*v^4 + ((16/3*a^6 - 8/3)*u + 8/3*a^8 - 16/3*a^6)*v^3 + (4*a^8*u^2 + (8*a^12 + 8*a^10 - 12*a^8 - 16*a^6 + 4)*u + 6*a^14 - 10*a^12 - 8*a^10 + 16*a^6)*v^2 + (2/3*u^4 + (16/3*a^6 - 8/3)*u^3 + (8*a^12 + 8*a^10 - 12*a^8 - 16*a^6 + 4)*u^2 + (12*a^16 - 4*a^14 - 64*a^10 + 32*a^8 + 112/3*a^6 - 16/3)*u + 8*a^18 - 18*a^16 - 8*a^14 + 38/3*a^12 + 32*a^10 - 8*a^8 - 64/3*a^6)*v + (8/3*a^8 - 16/3*a^6)*u^3 + (6*a^14 - 10*a^12 - 8*a^10 + 16*a^6)*u^2 + (8*a^18 - 18*a^16 - 8*a^14 + 38/3*a^12 + 32*a^10 - 8*a^8 - 64/3*a^6)*u + 5*a^20 - 22*a^18 + 99/4*a^16 + 9*a^14 - 64/3*a^12 - 17/12*a^8 + 16/3*a^6 + 4/3)*t^5 + (((-8/3*a^8 + 4*a^4 - 8*a^2 + 4)*u)*v^3 + ((-8*a^12 - 16*a^8 + 16*a^6 - 12*a^4 + 24*a^2 - 12)*u + 2*a^16 - 12*a^14 + 18*a^12 - 16*a^10 + 32*a^8 - 16*a^6)*v^2 + ((-8/3*a^8 + 4*a^4 - 8*a^2 + 4)*u^3 + (-8*a^12 - 16*a^8 + 16*a^6 - 12*a^4 + 24*a^2 - 12)*u^2 + (4*a^18 - 8*a^16 - 12*a^14 - 40*a^12 + 48*a^10 + 104*a^8 - 80*a^6 + 24*a^4 - 48*a^2 + 24)*u + 6*a^20 - 22*a^18 - 2*a^16 + 58*a^14 - 16*a^12 + 8*a^10 - 88*a^8 + 48*a^6)*v + (2*a^16 - 12*a^14 + 18*a^12 - 16*a^10 + 32*a^8 - 16*a^6)*u^2 + (6*a^20 - 22*a^18 - 2*a^16 + 58*a^14 - 16*a^12 + 8*a^10 - 88*a^8 + 48*a^6)*u + 6*a^22 - 67/2*a^20 + 52*a^18 + 25/2*a^16 - 251/3*a^14 + 105/2*a^12 - 33*a^10 + 275/6*a^8 - 16*a^6 - 8*a^4 + 16*a^2 - 8)*t^4 + (((-2*a^16 + 6*a^14 - 16/3*a^12 + 16*a^10 - 22*a^8 + 8*a^6 + 20*a^4 - 24*a^2 + 22/3)*u)*v^2 + ((-2*a^16 + 6*a^14 - 16/3*a^12 + 16*a^10 - 22*a^8 + 8*a^6 + 20*a^4 - 24*a^2 + 22/3)*u^2 + (4*a^18 + 8*a^16 - 72*a^14 + 424/3*a^12 - 64*a^10 + 4*a^8 - 80*a^4 + 96*a^2 - 88/3)*u + 2/3*a^24 - 22/3*a^20 + 16/3*a^18 + 10*a^16 + 28*a^14 - 206/3*a^12 - 16*a^10 + 184/3*a^8 - 64/3*a^6)*v + (2/3*a^24 - 22/3*a^20 + 16/3*a^18 + 10*a^16 + 28*a^14 - 206/3*a^12 - 16*a^10 + 184/3*a^8 - 64/3*a^6)*u + a^26 - 7/12*a^24 - 18*a^22 + 281/6*a^20 - 65/3*a^18 - 63/2*a^16 + 7*a^14 + 25/6*a^12 + 67*a^10 - 1003/12*a^8 + 80/3*a^6 + 40*a^4 - 48*a^2 + 44/3)*t^3 + (((8/3*a^20 - 20/3*a^18 - 4*a^16 + 28*a^14 - 128/3*a^12 + 32*a^10 - 116/3*a^8 + 16*a^6 + 32*a^4 - 32*a^2 + 8)*u)*v + 1/12*a^32 - 2*a^26 - 3/2*a^24 + 18*a^22 - 70/3*a^20 + 23/3*a^18 + 25/4*a^16 - 85/3*a^14 + 134/3*a^12 - 34*a^10 + 235/6*a^8 - 16*a^6 - 32*a^4 + 32*a^2 - 8)*t^2 + ((-1/12*a^32 + a^26 + 3/4*a^24 - 6*a^22 + 5*a^20)*y)*t)*w^4 + (((2*a^8*u*v^2 + (2*a^8*u^2 + (12*a^10 - 8*a^8)*u + 4*a^12 - 8*a^10)*v + (4*a^12 - 8*a^10)*u + 3*a^14 - 9*a^12 + 4*a^10 + 4*a^8)*y - 2*a^8*u*v^3 + ((4*a^12 - 4*a^8)*u^2 + (4*a^14 - 8*a^12 - 12*a^10 + 10*a^8)*u - 4*a^12 + 8*a^10)*v^2 + (-2*a^8*u^3 + (4*a^14 - 8*a^12 - 12*a^10 + 10*a^8)*u^2 + (-20*a^14 + 8*a^12 + 28*a^10 - 8*a^8)*u - 4*a^16 + 5*a^14 + 13*a^12 - 12*a^10 - 4*a^8)*v + (-4*a^12 + 8*a^10)*u^2 + (-4*a^16 + 5*a^14 + 13*a^12 - 12*a^10 - 4*a^8)*u - 3*a^18 + 9*a^16 - a^14 - 13*a^12 + 4*a^10 + 4*a^8)*t^4 + ((((-4*a^12 - 8*a^10 + 4*a^8)*u)*v + a^16 - 6*a^14 + 5*a^12 + 8*a^10 - 4*a^8)*y + ((2*a^16 - 8*a^14 + 8*a^12 + 8*a^10 - 4*a^8)*u)*v^2 + ((2*a^16 - 8*a^14 + 8*a^12 + 8*a^10 - 4*a^8)*u^2 + (-8*a^16 + 44*a^14 - 20*a^12 - 8*a^10 + 4*a^8)*u + 3*a^16 - 2*a^14 - 5*a^12 - 8*a^10 + 4*a^8)*v + (3*a^16 - 2*a^14 - 5*a^12 - 8*a^10 + 4*a^8)*u - a^20 + 5*a^18 - a^16 - 18*a^14 + 13*a^12 + 8*a^10 - 4*a^8)*t^3 + ((-a^16 + 3*a^14)*y^2 + ((6*a^16 - 22*a^14 + 4*a^12)*u)*v + a^20 - 2*a^18 - 7*a^16 + 16*a^14 - 4*a^12)*t^2 + ((-a^16 + 3*a^14)*y)*t)*w^3 + ((1/2*a^12*y^3 + (((a^14 - 3*a^12)*u)*v)*y + ((-a^14 + 3/2*a^12)*u)*v^2 + ((-a^14 + 3/2*a^12)*u^2 + a^14*u)*v)*t^3 + ((1/2*a^16 - a^14 + 3/2*a^12)*y^2 + ((-a^16 + 2*a^14 - 3*a^12)*u)*v)*t^2 + ((1/2*a^16 - a^14 + a^12)*y)*t)*w^2,\n",
       " ((4/15*u*v^5 + ((8/3*a^8 - 4/3)*u + 4/3*a^10 - 8/3*a^8)*v^4 + (8/3*a^12*u^2 + (16/3*a^16 + 16/3*a^14 - 8*a^12 - 32/3*a^8 + 8/3)*u + 4*a^18 - 20/3*a^16 - 16/3*a^14 + 16/3*a^12 - 16/3*a^10 + 32/3*a^8)*v^3 + (8/3*a^12*u^3 + (8*a^18 - 16*a^12)*u^2 + (12*a^22 - 4*a^20 - 32*a^16 - 16*a^14 + 104/3*a^12 + 16*a^8 - 8/3)*u + 8*a^24 - 18*a^22 + 4*a^20 - 58/3*a^18 + 36*a^16 + 16*a^14 - 64/3*a^12 + 8*a^10 - 16*a^8)*v^2 + (4/15*u^5 + (8/3*a^8 - 4/3)*u^4 + (16/3*a^16 + 16/3*a^14 - 8*a^12 - 32/3*a^8 + 8/3)*u^3 + (12*a^22 - 4*a^20 - 32*a^16 - 16*a^14 + 104/3*a^12 + 16*a^8 - 8/3)*u^2 + (16*a^26 - 12*a^24 - 48*a^22 + 244/3*a^20 - 128*a^18 + 112*a^16 + 160/3*a^14 - 176/3*a^12 - 24*a^8 + 8/3)*u + 10*a^28 - 28*a^26 - 5/2*a^24 + 50*a^22 - 136/3*a^20 + 224/3*a^18 - 131/2*a^16 - 80/3*a^14 + 80/3*a^12 - 16/3*a^10 + 40/3*a^8)*v + (4/3*a^10 - 8/3*a^8)*u^4 + (4*a^18 - 20/3*a^16 - 16/3*a^14 + 16/3*a^12 - 16/3*a^10 + 32/3*a^8)*u^3 + (8*a^24 - 18*a^22 + 4*a^20 - 58/3*a^18 + 36*a^16 + 16*a^14 - 64/3*a^12 + 8*a^10 - 16*a^8)*u^2 + (10*a^28 - 28*a^26 - 5/2*a^24 + 50*a^22 - 136/3*a^20 + 224/3*a^18 - 131/2*a^16 - 80/3*a^14 + 80/3*a^12 - 16/3*a^10 + 40/3*a^8)*u + 6*a^30 - 30*a^28 + 43*a^26 + 3*a^24 - 81/2*a^22 + 28/3*a^20 + 13/12*a^18 + 41/3*a^16 - 16/3*a^12 + 161/60*a^10 - 8/3*a^8 - 8/15)*t^6 + (((-4/3*a^10 + 8/3*a^4 - 16/3*a^2 + 8/3)*u)*v^4 + ((-16/3*a^16 - 32/3*a^14 + 64/3*a^12 - 80/3*a^10 + 16*a^8 - 32/3*a^4 + 64/3*a^2 - 32/3)*u + 4/3*a^20 - 8*a^18 + 20/3*a^16 + 56/3*a^14 - 112/3*a^12 + 40*a^10 - 16*a^8)*v^3 + ((-8*a^18 + 16*a^16 - 32*a^14 + 16*a^12)*u^2 + (4*a^24 - 16*a^22 + 4*a^20 - 16*a^18 - 32*a^16 + 144*a^14 - 112*a^12 + 88*a^10 - 48*a^8 + 16*a^4 - 32*a^2 + 16)*u + 6*a^26 - 28*a^24 + 36*a^22 - 32*a^20 + 66*a^18 + 8*a^16 - 136*a^14 + 144*a^12 - 120*a^10 + 48*a^8)*v^2 + ((-4/3*a^10 + 8/3*a^4 - 16/3*a^2 + 8/3)*u^4 + (-16/3*a^16 - 32/3*a^14 + 64/3*a^12 - 80/3*a^10 + 16*a^8 - 32/3*a^4 + 64/3*a^2 - 32/3)*u^3 + (4*a^24 - 16*a^22 + 4*a^20 - 16*a^18 - 32*a^16 + 144*a^14 - 112*a^12 + 88*a^10 - 48*a^8 + 16*a^4 - 32*a^2 + 16)*u^2 + (12*a^28 - 24*a^26 - 56*a^24 + 176*a^22 - 284*a^20 + 288*a^18 + 464/3*a^16 - 1376/3*a^14 + 832/3*a^12 - 640/3*a^10 + 112*a^8 - 64/3*a^4 + 128/3*a^2 - 64/3)*u + 12*a^30 - 47*a^28 + 14*a^26 + 129*a^24 - 542/3*a^22 + 611/3*a^20 - 688/3*a^18 - 137/3*a^16 + 776/3*a^14 - 592/3*a^12 + 152*a^10 - 64*a^8)*v + (4/3*a^20 - 8*a^18 + 20/3*a^16 + 56/3*a^14 - 112/3*a^12 + 40*a^10 - 16*a^8)*u^3 + (6*a^26 - 28*a^24 + 36*a^22 - 32*a^20 + 66*a^18 + 8*a^16 - 136*a^14 + 144*a^12 - 120*a^10 + 48*a^8)*u^2 + (12*a^30 - 47*a^28 + 14*a^26 + 129*a^24 - 542/3*a^22 + 611/3*a^20 - 688/3*a^18 - 137/3*a^16 + 776/3*a^14 - 592/3*a^12 + 152*a^10 - 64*a^8)*u + 10*a^32 - 60*a^30 + 110*a^28 - 19/2*a^26 - 166*a^24 + 893/6*a^22 - 73*a^20 + 289/3*a^18 - 34*a^16 - 353/6*a^14 + 209/3*a^12 - 305/6*a^10 + 16*a^8 + 16/3*a^4 - 32/3*a^2 + 16/3)*t^5 + (((-4/3*a^20 + 4*a^18 - 8*a^14 + 56/3*a^12 - 8*a^10 - 12*a^8 + 32*a^4 - 32*a^2 + 28/3)*u)*v^3 + ((-8/3*a^24 + 8*a^22 - 4/3*a^20 + 4*a^18 + 24*a^16 - 40*a^14 + 152/3*a^12 - 72*a^10 + 196/3*a^8 - 96*a^4 + 96*a^2 - 28)*u + 2/3*a^30 - 12*a^26 + 92/3*a^24 - 110/3*a^22 + 172/3*a^20 - 224/3*a^18 - 4*a^16 + 104*a^14 - 464/3*a^12 + 332/3*a^10 - 88/3*a^8)*v^2 + ((-4/3*a^20 + 4*a^18 - 8*a^14 + 56/3*a^12 - 8*a^10 - 12*a^8 + 32*a^4 - 32*a^2 + 28/3)*u^3 + (-8/3*a^24 + 8*a^22 - 4/3*a^20 + 4*a^18 + 24*a^16 - 40*a^14 + 152/3*a^12 - 72*a^10 + 196/3*a^8 - 96*a^4 + 96*a^2 - 28)*u^2 + (4/3*a^32 - 16/3*a^28 + 8*a^26 + 76/3*a^24 - 192*a^22 + 1240/3*a^20 - 208*a^18 - 728/3*a^16 + 336*a^14 - 1264/3*a^12 + 432*a^10 - 656/3*a^8 + 192*a^4 - 192*a^2 + 56)*u + 2*a^34 + 1/6*a^32 - 82/3*a^30 + 97/3*a^28 + 48*a^26 - 95*a^24 + 436/3*a^22 - 1043/3*a^20 + 868/3*a^18 + 469/6*a^16 - 304*a^14 + 416*a^12 - 952/3*a^10 + 88*a^8)*v + (2/3*a^30 - 12*a^26 + 92/3*a^24 - 110/3*a^22 + 172/3*a^20 - 224/3*a^18 - 4*a^16 + 104*a^14 - 464/3*a^12 + 332/3*a^10 - 88/3*a^8)*u^2 + (2*a^34 + 1/6*a^32 - 82/3*a^30 + 97/3*a^28 + 48*a^26 - 95*a^24 + 436/3*a^22 - 1043/3*a^20 + 868/3*a^18 + 469/6*a^16 - 304*a^14 + 416*a^12 - 952/3*a^10 + 88*a^8)*u + 2*a^36 - 3/2*a^34 - 125/3*a^32 + 755/6*a^30 - 250/3*a^28 - 1243/12*a^26 + 377/3*a^24 - 313/6*a^22 + 551/3*a^20 - 1451/6*a^18 + 139/3*a^16 + 162*a^14 - 242*a^12 + 1703/12*a^10 - 16/3*a^8 - 64*a^4 + 64*a^2 - 56/3)*t^4 + (((-2/3*a^30 + 6*a^26 - 32/3*a^24 + 44/3*a^22 - 28*a^20 + 88/3*a^18 - 64/3*a^16 + 8/3*a^14 + 8*a^12 + 100/3*a^10 - 68*a^8 + 232/3*a^4 - 176/3*a^2 + 40/3)*u)*v^2 + ((-2/3*a^30 + 6*a^26 - 32/3*a^24 + 44/3*a^22 - 28*a^20 + 88/3*a^18 - 64/3*a^16 + 8/3*a^14 + 8*a^12 + 100/3*a^10 - 68*a^8 + 232/3*a^4 - 176/3*a^2 + 40/3)*u^2 + (4*a^32 + 8/3*a^30 - 4*a^28 - 72*a^26 + 104/3*a^24 + 760/3*a^22 - 288*a^20 - 136/3*a^18 + 376/3*a^16 - 224/3*a^14 + 192*a^12 - 976/3*a^10 + 320*a^8 - 928/3*a^4 + 704/3*a^2 - 160/3)*u + 1/6*a^40 - 2*a^34 - 17/3*a^32 + 50/3*a^30 - 8*a^28 + 20*a^26 + 179/6*a^24 - 252*a^22 + 956/3*a^20 - 308/3*a^18 - 7*a^16 + 224/3*a^14 - 544/3*a^12 + 136*a^10 - 32*a^8)*v + (1/6*a^40 - 2*a^34 - 17/3*a^32 + 50/3*a^30 - 8*a^28 + 20*a^26 + 179/6*a^24 - 252*a^22 + 956/3*a^20 - 308/3*a^18 - 7*a^16 + 224/3*a^14 - 544/3*a^12 + 136*a^10 - 32*a^8)*u + 1/4*a^42 - 1/3*a^40 + 1/2*a^38 - 6*a^36 - 17/4*a^34 + 202/3*a^32 - 205/2*a^30 + 77/3*a^28 + 451/12*a^26 - 248/3*a^24 + 1415/6*a^22 - 925/3*a^20 + 2395/12*a^18 - 206/3*a^16 - 509/6*a^14 + 475/3*a^12 - 85/6*a^10 - 120*a^8 + 464/3*a^4 - 352/3*a^2 + 80/3)*t^3 + (((1/5*a^40 - 7/3*a^32 - 16/3*a^30 + 20/3*a^28 + 28*a^26 - 64/3*a^24 - 104/3*a^22 + 32/3*a^20 + 152/3*a^18 - 128/3*a^16 + 64/3*a^14 - 64/3*a^12 + 80*a^10 - 112*a^8 + 96*a^4 - 64*a^2 + 64/5)*u)*v + 1/60*a^50 - 1/2*a^42 - 1/5*a^40 - a^38 + 6*a^36 + 35/4*a^34 - 113/3*a^32 + 110/3*a^30 - 19/3*a^28 - 63/2*a^26 + 24*a^24 + 28*a^22 - 11/3*a^20 - 331/6*a^18 + 128/3*a^16 - 55/3*a^14 + 58/3*a^12 - 398/5*a^10 + 112*a^8 - 96*a^4 + 64*a^2 - 64/5)*t^2 + ((-1/60*a^50 + 1/4*a^42 + 1/2*a^38 - 2*a^36 - 3*a^34 + 10*a^32 - 6*a^30)*y)*t)*w^5 + (((4/3*a^10*u*v^3 + ((8*a^14 - 4*a^10)*u + 4*a^16 - 8*a^14)*v^2 + (4/3*a^10*u^3 + (8*a^14 - 4*a^10)*u^2 + (8*a^18 + 16*a^16 - 40*a^14 + 8*a^10)*u + 6*a^20 - 10*a^18 - 16*a^16 + 24*a^14)*v + (4*a^16 - 8*a^14)*u^2 + (6*a^20 - 10*a^18 - 16*a^16 + 24*a^14)*u + 4*a^22 - 15*a^20 + 12*a^18 + 25/3*a^16 - 8*a^14 - 8/3*a^10)*y - 4/3*a^10*u*v^4 + ((4/3*a^16 - 4/3*a^10)*u^2 + (-8/3*a^16 - 8*a^14 + 16/3*a^10)*u - 4*a^16 + 8*a^14)*v^3 + ((4/3*a^16 - 4/3*a^10)*u^3 + (8*a^20 - 8*a^16 - 16*a^14 + 8*a^10)*u^2 + (8*a^22 - 24*a^20 - 8*a^18 - 8*a^16 + 56*a^14 - 12*a^10)*u - 4*a^22 + 2*a^20 + 10*a^18 + 20*a^16 - 32*a^14)*v^2 + (-4/3*a^10*u^4 + (-8/3*a^16 - 8*a^14 + 16/3*a^10)*u^3 + (8*a^22 - 24*a^20 - 8*a^18 - 8*a^16 + 56*a^14 - 12*a^10)*u^2 + (-8*a^24 - 32*a^22 + 52*a^20 + 28*a^18 + 32*a^16 - 88*a^14 + 8*a^10)*u - 6*a^26 + 10*a^24 + 12*a^22 - 3*a^20 - 22*a^18 - 73/3*a^16 + 32*a^14 + 8/3*a^10)*v + (-4*a^16 + 8*a^14)*u^3 + (-4*a^22 + 2*a^20 + 10*a^18 + 20*a^16 - 32*a^14)*u^2 + (-6*a^26 + 10*a^24 + 12*a^22 - 3*a^20 - 22*a^18 - 73/3*a^16 + 32*a^14 + 8/3*a^10)*u - 4*a^28 + 15*a^26 - 12*a^24 - 13/3*a^22 - 7*a^20 + 12*a^18 + 11*a^16 - 8*a^14 - 8/3*a^10)*t^5 + ((((-4*a^16 + 4*a^14 - 8*a^12 + 4*a^10)*u)*v^2 + ((-4*a^16 + 4*a^14 - 8*a^12 + 4*a^10)*u^2 + (-32*a^16 + 8*a^14 + 32*a^12 - 16*a^10)*u + 2*a^22 - 8*a^20 - 6*a^18 + 36*a^16 - 16*a^14)*v + (2*a^22 - 8*a^20 - 6*a^18 + 36*a^16 - 16*a^14)*u + 3*a^24 - 16*a^22 + 20*a^20 + 14*a^18 - 33*a^16 + 16*a^14 - 16*a^12 + 8*a^10)*y + ((16/3*a^16 - 4*a^14 + 8*a^12 - 4*a^10)*u)*v^3 + ((-16*a^18 + 16*a^16 - 8*a^14 + 16*a^12 - 8*a^10)*u^2 + (8*a^24 - 24*a^22 + 16*a^20 + 32*a^18 + 8*a^16 - 4*a^14 - 40*a^12 + 20*a^10)*u + 2*a^22 + 6*a^18 - 36*a^16 + 16*a^14)*v^2 + ((16/3*a^16 - 4*a^14 + 8*a^12 - 4*a^10)*u^3 + (8*a^24 - 24*a^22 + 16*a^20 + 32*a^18 + 8*a^16 - 4*a^14 - 40*a^12 + 20*a^10)*u^2 + (-32*a^24 + 132*a^22 - 64*a^20 - 52*a^18 - 64*a^16 + 40*a^14 + 32*a^12 - 16*a^10)*u - 2*a^28 + 6*a^26 + 17*a^24 - 50*a^22 + 12*a^20 - 20*a^18 + 69*a^16 - 32*a^14 + 16*a^12 - 8*a^10)*v + (2*a^22 + 6*a^18 - 36*a^16 + 16*a^14)*u^2 + (-2*a^28 + 6*a^26 + 17*a^24 - 50*a^22 + 12*a^20 - 20*a^18 + 69*a^16 - 32*a^14 + 16*a^12 - 8*a^10)*u - 3*a^30 + 14*a^28 - 11*a^26 - 25*a^24 + 76/3*a^22 + 4*a^20 + 30*a^18 - 131/3*a^16 + 16*a^14 - 16*a^12 + 8*a^10)*t^4 + ((((-4/3*a^22 + 4*a^20 + 8*a^16 + 8*a^14 - 16*a^12 + 16/3*a^10)*u)*v + 1/3*a^28 - 6*a^24 + 40/3*a^22 - 5*a^20 - 2*a^18 - 22/3*a^16 - 8*a^14 + 16*a^12 - 16/3*a^10)*y + ((4/3*a^28 - 14*a^24 + 70/3*a^22 - 4*a^20 - 16*a^18 - 4/3*a^16 - 8*a^14 + 16*a^12 - 16/3*a^10)*u)*v^2 + ((4/3*a^28 - 14*a^24 + 70/3*a^22 - 4*a^20 - 16*a^18 - 4/3*a^16 - 8*a^14 + 16*a^12 - 16/3*a^10)*u^2 + (-16/3*a^28 + 80*a^24 - 484/3*a^22 + 28*a^20 + 64*a^18 - 56/3*a^16 + 8*a^14 - 16*a^12 + 16/3*a^10)*u + 5/3*a^28 - 24*a^24 + 116/3*a^22 - 11*a^20 + 2*a^18 + 22/3*a^16 + 8*a^14 - 16*a^12 + 16/3*a^10)*v + (5/3*a^28 - 24*a^24 + 116/3*a^22 - 11*a^20 + 2*a^18 + 22/3*a^16 + 8*a^14 - 16*a^12 + 16/3*a^10)*u - 1/3*a^34 + 5*a^30 - 7*a^28 - 7*a^26 + 4*a^24 + 71/3*a^22 + 3*a^20 - 34*a^18 + 6*a^16 - 8*a^14 + 16*a^12 - 16/3*a^10)*t^3 + ((-1/3*a^28 + 3*a^24 - 4*a^22)*y^2 + ((10/3*a^28 - 34*a^24 + 52*a^22 - 16*a^18 + 16/3*a^16)*u)*v + 1/3*a^34 - 2*a^30 - 8/3*a^28 + 3*a^26 + 30*a^24 - 130/3*a^22 + 16*a^18 - 16/3*a^16)*t^2 + ((-1/3*a^28 + 3*a^24 - 4*a^22)*y)*t)*w^4 + (((2*a^16*u*v + a^18 - 2*a^16)*y^2 + (((2*a^18 - 4*a^16)*u)*v^2 + ((2*a^18 - 4*a^16)*u^2 + (-8*a^18 + 4*a^16)*u - 2*a^18 + 4*a^16)*v + (-2*a^18 + 4*a^16)*u - 2*a^20 + 6*a^18 - 4*a^16)*y + ((a^20 - 2*a^18 + 2*a^16)*u)*v^3 + ((-4*a^18 + 4*a^16)*u^2 + (-3*a^20 + 10*a^18 - 4*a^16)*u + a^18 - 2*a^16)*v^2 + ((a^20 - 2*a^18 + 2*a^16)*u^3 + (-3*a^20 + 10*a^18 - 4*a^16)*u^2 + (6*a^20 - 6*a^18 - 2*a^16)*u + 2*a^20 - 6*a^18 + 4*a^16)*v + (a^18 - 2*a^16)*u^2 + (2*a^20 - 6*a^18 + 4*a^16)*u + a^22 - 4*a^20 + 5*a^18 - 2*a^16)*t^4 + (-a^18*y^3 + (((6*a^18 + 2*a^16)*u)*v + a^22 - 2*a^20 + a^18 - 2*a^16)*y + ((a^24 - 2*a^22 + 3*a^20 - 3*a^18 - 2*a^16)*u)*v^2 + ((a^24 - 2*a^22 + 3*a^20 - 3*a^18 - 2*a^16)*u^2 + (-4*a^24 + 8*a^22 - 12*a^20 + 2*a^16)*u - a^22 + 2*a^20 - a^18 + 2*a^16)*v + (-a^22 + 2*a^20 - a^18 + 2*a^16)*u - a^26 + 4*a^24 - 6*a^22 + 4*a^20 + a^18 - 2*a^16)*t^3 + ((-a^22 + 2*a^20 - 3*a^18)*y^2 + ((2*a^24 - 2*a^22 + 6*a^18)*u)*v + a^26 - 4*a^24 + 6*a^22 - 4*a^20)*t^2 + ((-a^22 + 2*a^20 - 2*a^18)*y)*t)*w^3,\n",
       " ((4/45*u*v^6 + ((16/15*a^10 - 8/15)*u + 8/15*a^12 - 16/15*a^10)*v^5 + (4/3*a^16*u^2 + (8/3*a^20 + 8/3*a^18 - 4*a^16 - 16/3*a^10 + 4/3)*u + 2*a^22 - 10/3*a^20 - 8/3*a^18 + 8/3*a^16 - 8/3*a^12 + 16/3*a^10)*v^4 + (16/9*a^18*u^3 + (16/3*a^24 - 16/3*a^18 - 16/3*a^16)*u^2 + (8*a^28 - 8/3*a^26 - 32/3*a^22 - 32/3*a^20 - 32/9*a^18 + 16*a^16 + 32/3*a^10 - 16/9)*u + 16/3*a^30 - 12*a^28 + 8/3*a^26 - 44/9*a^24 + 8/3*a^22 + 40/3*a^20 + 64/9*a^18 - 32/3*a^16 + 16/3*a^12 - 32/3*a^10)*v^3 + (4/3*a^16*u^4 + (16/3*a^24 - 16/3*a^18 - 16/3*a^16)*u^3 + (12*a^30 - 4*a^28 - 32*a^24 + 16*a^18 + 16*a^16)*u^2 + (16*a^34 - 12*a^32 - 24*a^30 + 52/3*a^28 - 40*a^26 + 24*a^24 + 128/3*a^22 + 16*a^20 - 16/3*a^18 - 92/3*a^16 - 32/3*a^10 + 4/3)*u + 10*a^36 - 28*a^34 + 27/2*a^32 - 2*a^30 - 4/3*a^28 + 40*a^26 + 7/6*a^24 - 92/3*a^22 - 20*a^20 - 16/3*a^18 + 56/3*a^16 - 16/3*a^12 + 32/3*a^10)*v^2 + (4/45*u^6 + (16/15*a^10 - 8/15)*u^5 + (8/3*a^20 + 8/3*a^18 - 4*a^16 - 16/3*a^10 + 4/3)*u^4 + (8*a^28 - 8/3*a^26 - 32/3*a^22 - 32/3*a^20 - 32/9*a^18 + 16*a^16 + 32/3*a^10 - 16/9)*u^3 + (16*a^34 - 12*a^32 - 24*a^30 + 52/3*a^28 - 40*a^26 + 24*a^24 + 128/3*a^22 + 16*a^20 - 16/3*a^18 - 92/3*a^16 - 32/3*a^10 + 4/3)*u^2 + (20*a^38 - 24*a^36 - 61*a^34 + 52*a^32 + 176*a^30 - 1072/3*a^28 + 815/3*a^26 + 32*a^24 - 128*a^22 - 64/3*a^20 - 16/3*a^18 + 112/3*a^16 + 176/15*a^10 - 16/15)*u + 12*a^40 - 40*a^38 + 10*a^36 + 73*a^34 - 44*a^32 - 76*a^30 + 319/2*a^28 - 133*a^26 - 59/3*a^24 + 56*a^22 + 481/30*a^20 + 16/3*a^18 - 16*a^16 + 8/3*a^12 - 32/5*a^10)*v + (8/15*a^12 - 16/15*a^10)*u^5 + (2*a^22 - 10/3*a^20 - 8/3*a^18 + 8/3*a^16 - 8/3*a^12 + 16/3*a^10)*u^4 + (16/3*a^30 - 12*a^28 + 8/3*a^26 - 44/9*a^24 + 8/3*a^22 + 40/3*a^20 + 64/9*a^18 - 32/3*a^16 + 16/3*a^12 - 32/3*a^10)*u^3 + (10*a^36 - 28*a^34 + 27/2*a^32 - 2*a^30 - 4/3*a^28 + 40*a^26 + 7/6*a^24 - 92/3*a^22 - 20*a^20 - 16/3*a^18 + 56/3*a^16 - 16/3*a^12 + 32/3*a^10)*u^2 + (12*a^40 - 40*a^38 + 10*a^36 + 73*a^34 - 44*a^32 - 76*a^30 + 319/2*a^28 - 133*a^26 - 59/3*a^24 + 56*a^22 + 481/30*a^20 + 16/3*a^18 - 16*a^16 + 8/3*a^12 - 32/5*a^10)*u + 7*a^42 - 39*a^40 + 135/2*a^38 - 79/6*a^36 - 63*a^34 + 33*a^32 + 59/12*a^30 + 443/24*a^28 - 17/3*a^26 - 107/9*a^24 + 81/20*a^22 - 101/15*a^20 + 16/9*a^18 + 8/3*a^16 - 77/72*a^12 + 16/15*a^10 + 8/45)*t^7 + (((-8/15*a^12 + 4/3*a^4 - 8/3*a^2 + 4/3)*u)*v^5 + ((-8/3*a^20 - 16/3*a^18 + 8/3*a^16 + 32/3*a^14 - 56/3*a^12 + 32/3*a^10 - 20/3*a^4 + 40/3*a^2 - 20/3)*u + 2/3*a^24 - 4*a^22 + 10/3*a^20 + 16/3*a^18 + 8/3*a^16 - 64/3*a^14 + 80/3*a^12 - 32/3*a^10)*v^4 + ((-16/3*a^24 + 16/3*a^22 - 8/3*a^20 - 32/3*a^18 + 8*a^16)*u^2 + (8/3*a^30 - 32/3*a^28 - 8*a^26 + 16*a^24 - 64/3*a^22 + 8/3*a^20 + 208/3*a^18 - 104/3*a^16 - 128/3*a^14 + 80*a^12 - 128/3*a^10 + 40/3*a^4 - 80/3*a^2 + 40/3)*u + 4*a^32 - 56/3*a^30 + 16*a^28 + 68/3*a^26 - 128/3*a^24 + 52*a^22 - 20/3*a^20 - 176/3*a^18 + 16/3*a^16 + 256/3*a^14 - 320/3*a^12 + 128/3*a^10)*v^3 + ((-16/3*a^24 + 16/3*a^22 - 8/3*a^20 - 32/3*a^18 + 8*a^16)*u^3 + (4*a^32 - 24*a^30 + 36*a^28 - 64*a^26 + 64*a^24 - 32*a^22 + 16*a^20 + 64*a^18 - 48*a^16)*u^2 + (12*a^36 - 36*a^34 + 4*a^32 + 28*a^30 - 56*a^28 + 96*a^26 + 32/3*a^24 + 64/3*a^22 - 8/3*a^20 - 656/3*a^18 + 120*a^16 + 64*a^14 - 368/3*a^12 + 64*a^10 - 40/3*a^4 + 80/3*a^2 - 40/3)*u + 12*a^38 - 55*a^36 + 76*a^34 - 55*a^32 + 230/3*a^30 - 17/3*a^28 - 242/3*a^26 + 25*a^24 - 236/3*a^22 + 16/3*a^20 + 496/3*a^18 - 48*a^16 - 128*a^14 + 160*a^12 - 64*a^10)*v^2 + ((-8/15*a^12 + 4/3*a^4 - 8/3*a^2 + 4/3)*u^5 + (-8/3*a^20 - 16/3*a^18 + 8/3*a^16 + 32/3*a^14 - 56/3*a^12 + 32/3*a^10 - 20/3*a^4 + 40/3*a^2 - 20/3)*u^4 + (8/3*a^30 - 32/3*a^28 - 8*a^26 + 16*a^24 - 64/3*a^22 + 8/3*a^20 + 208/3*a^18 - 104/3*a^16 - 128/3*a^14 + 80*a^12 - 128/3*a^10 + 40/3*a^4 - 80/3*a^2 + 40/3)*u^3 + (12*a^36 - 36*a^34 + 4*a^32 + 28*a^30 - 56*a^28 + 96*a^26 + 32/3*a^24 + 64/3*a^22 - 8/3*a^20 - 656/3*a^18 + 120*a^16 + 64*a^14 - 368/3*a^12 + 64*a^10 - 40/3*a^4 + 80/3*a^2 - 40/3)*u^2 + (24*a^40 - 54*a^38 - 88*a^36 + 242*a^34 + 124*a^32 - 2314/3*a^30 + 2492/3*a^28 + 50*a^26 - 720*a^24 + 832/3*a^22 - 424/3*a^20 + 1328/3*a^18 - 200*a^16 - 96*a^14 + 560/3*a^12 - 96*a^10 + 40/3*a^4 - 80/3*a^2 + 40/3)*u + 20*a^42 - 86*a^40 + 54*a^38 + 421/2*a^36 - 323*a^34 + 109/2*a^32 + 718/3*a^30 - 667/2*a^28 + 83/3*a^26 + 1489/6*a^24 - 106/3*a^22 + 41*a^20 - 640/3*a^18 + 72*a^16 + 96*a^14 - 128*a^12 + 160/3*a^10)*v + (2/3*a^24 - 4*a^22 + 10/3*a^20 + 16/3*a^18 + 8/3*a^16 - 64/3*a^14 + 80/3*a^12 - 32/3*a^10)*u^4 + (4*a^32 - 56/3*a^30 + 16*a^28 + 68/3*a^26 - 128/3*a^24 + 52*a^22 - 20/3*a^20 - 176/3*a^18 + 16/3*a^16 + 256/3*a^14 - 320/3*a^12 + 128/3*a^10)*u^3 + (12*a^38 - 55*a^36 + 76*a^34 - 55*a^32 + 230/3*a^30 - 17/3*a^28 - 242/3*a^26 + 25*a^24 - 236/3*a^22 + 16/3*a^20 + 496/3*a^18 - 48*a^16 - 128*a^14 + 160*a^12 - 64*a^10)*u^2 + (20*a^42 - 86*a^40 + 54*a^38 + 421/2*a^36 - 323*a^34 + 109/2*a^32 + 718/3*a^30 - 667/2*a^28 + 83/3*a^26 + 1489/6*a^24 - 106/3*a^22 + 41*a^20 - 640/3*a^18 + 72*a^16 + 96*a^14 - 128*a^12 + 160/3*a^10)*u + 15*a^44 - 97*a^42 + 409/2*a^40 - 74*a^38 - 551/2*a^36 + 671/2*a^34 - 1191/8*a^32 + 1729/12*a^30 - 703/8*a^28 - 611/6*a^26 + 152*a^24 - 1868/15*a^22 + 50*a^20 + 80/3*a^18 - 191/24*a^16 - 385/12*a^14 + 3973/120*a^12 - 32/3*a^10 - 8/3*a^4 + 16/3*a^2 - 8/3)*t^6 + (((-2/3*a^24 + 2*a^22 - 16/3*a^16 + 32/3*a^14 - 32/9*a^12 - 20/3*a^8 - 16/3*a^6 + 88/3*a^4 - 80/3*a^2 + 68/9)*u)*v^4 + ((-16/9*a^30 + 16/3*a^28 - 8/3*a^24 + 40/3*a^22 + 128/3*a^20 - 944/9*a^18 + 112/3*a^16 + 256/3*a^14 - 1024/9*a^12 + 112/3*a^10 + 80/3*a^8 + 64/3*a^6 - 352/3*a^4 + 320/3*a^2 - 272/9)*u + 4/9*a^36 - 8*a^32 + 160/9*a^30 - 8/3*a^28 - 104/3*a^26 + 608/9*a^24 - 160/3*a^22 - 188/3*a^20 + 944/9*a^18 + 48*a^16 - 192*a^14 + 440/3*a^12 - 112/3*a^10)*v^3 + ((-4*a^32 + 12*a^30 - 64/3*a^28 + 80*a^26 - 152*a^24 + 128*a^22 + 16*a^20 - 80*a^18 + 88/3*a^16)*u^2 + (4/3*a^40 - 44/3*a^36 + 80/3*a^34 - 16/3*a^32 - 12*a^30 + 188/3*a^28 - 192*a^26 + 1268/3*a^24 - 1196/3*a^22 - 728/3*a^20 + 584*a^18 - 168*a^16 - 320*a^14 + 1088/3*a^12 - 112*a^10 - 40*a^8 - 32*a^6 + 176*a^4 - 160*a^2 + 136/3)*u + 2*a^42 + 1/6*a^40 - 36*a^38 + 269/3*a^36 - 314/3*a^34 + 419/3*a^32 - 156*a^30 - 61/3*a^28 + 290*a^26 - 3355/6*a^24 + 1274/3*a^22 + 838/3*a^20 - 504*a^18 - 256/3*a^16 + 576*a^14 - 440*a^12 + 112*a^10)*v^2 + ((-2/3*a^24 + 2*a^22 - 16/3*a^16 + 32/3*a^14 - 32/9*a^12 - 20/3*a^8 - 16/3*a^6 + 88/3*a^4 - 80/3*a^2 + 68/9)*u^4 + (-16/9*a^30 + 16/3*a^28 - 8/3*a^24 + 40/3*a^22 + 128/3*a^20 - 944/9*a^18 + 112/3*a^16 + 256/3*a^14 - 1024/9*a^12 + 112/3*a^10 + 80/3*a^8 + 64/3*a^6 - 352/3*a^4 + 320/3*a^2 - 272/9)*u^3 + (4/3*a^40 - 44/3*a^36 + 80/3*a^34 - 16/3*a^32 - 12*a^30 + 188/3*a^28 - 192*a^26 + 1268/3*a^24 - 1196/3*a^22 - 728/3*a^20 + 584*a^18 - 168*a^16 - 320*a^14 + 1088/3*a^12 - 112*a^10 - 40*a^8 - 32*a^6 + 176*a^4 - 160*a^2 + 136/3)*u^2 + (4*a^44 + 3*a^42 - 112/3*a^40 - 26/3*a^38 + 184*a^36 - 130/3*a^34 - 2348/3*a^32 + 12962/9*a^30 - 704*a^28 - 907/3*a^26 - 288*a^24 + 1424/3*a^22 + 3232/3*a^20 - 14480/9*a^18 + 1168/3*a^16 + 2432/3*a^14 - 7808/9*a^12 + 784/3*a^10 + 160/3*a^8 + 128/3*a^6 - 704/3*a^4 + 640/3*a^2 - 544/9)*u + 4*a^46 + a^44 - 209/3*a^42 + 107*a^40 + 368/3*a^38 - 6601/18*a^36 + 616/3*a^34 + 359/3*a^32 - 3026/9*a^30 + 769/3*a^28 - 469/3*a^26 + 5798/9*a^24 - 544*a^22 - 3449/6*a^20 + 7712/9*a^18 + 32/3*a^16 - 704*a^14 + 568*a^12 - 448/3*a^10)*v + (4/9*a^36 - 8*a^32 + 160/9*a^30 - 8/3*a^28 - 104/3*a^26 + 608/9*a^24 - 160/3*a^22 - 188/3*a^20 + 944/9*a^18 + 48*a^16 - 192*a^14 + 440/3*a^12 - 112/3*a^10)*u^3 + (2*a^42 + 1/6*a^40 - 36*a^38 + 269/3*a^36 - 314/3*a^34 + 419/3*a^32 - 156*a^30 - 61/3*a^28 + 290*a^26 - 3355/6*a^24 + 1274/3*a^22 + 838/3*a^20 - 504*a^18 - 256/3*a^16 + 576*a^14 - 440*a^12 + 112*a^10)*u^2 + (4*a^46 + a^44 - 209/3*a^42 + 107*a^40 + 368/3*a^38 - 6601/18*a^36 + 616/3*a^34 + 359/3*a^32 - 3026/9*a^30 + 769/3*a^28 - 469/3*a^26 + 5798/9*a^24 - 544*a^22 - 3449/6*a^20 + 7712/9*a^18 + 32/3*a^16 - 704*a^14 + 568*a^12 - 448/3*a^10)*u + 10/3*a^48 - 5/2*a^46 - 165/2*a^44 + 561/2*a^42 - 2963/12*a^40 - 653/3*a^38 + 32605/72*a^36 - 718/3*a^34 + 323*a^32 - 17755/36*a^30 + 2995/24*a^28 + 1244/3*a^26 - 12791/18*a^24 + 6137/12*a^22 + 557/24*a^20 - 3325/18*a^18 - 377/4*a^16 + 1669/6*a^14 - 12433/72*a^12 + 112/3*a^10 - 40/3*a^8 - 32/3*a^6 + 176/3*a^4 - 160/3*a^2 + 136/9)*t^5 + (((-4/9*a^36 + 4*a^32 - 16/3*a^30 - 4*a^28 + 20*a^26 - 32*a^24 + 12*a^22 + 24*a^20 - 32/3*a^18 - 200/3*a^16 + 304/3*a^14 - 212/3*a^12 + 88*a^10 - 88*a^8 - 128/3*a^6 + 428/3*a^4 - 280/3*a^2 + 20)*u)*v^3 + ((20/3*a^36 - 40/3*a^34 - 4*a^32 + 40/3*a^30 - 4/3*a^28 + 92/3*a^26 - 320/3*a^24 - 44/3*a^22 + 320*a^20 - 1312/3*a^18 + 736/3*a^16 + 16/3*a^14 - 68/3*a^12 - 632/3*a^10 + 264*a^8 + 128*a^6 - 428*a^4 + 280*a^2 - 60)*u + 1/6*a^48 - 4*a^42 - 5/3*a^40 + 100/3*a^38 - 64*a^36 + 286/3*a^34 - 833/6*a^32 + 362/3*a^30 - 62/3*a^28 - 476/3*a^26 + 893/3*a^24 + 20/3*a^22 - 1562/3*a^20 + 1408/3*a^18 + 328/3*a^16 - 1280/3*a^14 + 784/3*a^12 - 160/3*a^10)*v^2 + ((-4/9*a^36 + 4*a^32 - 16/3*a^30 - 4*a^28 + 20*a^26 - 32*a^24 + 12*a^22 + 24*a^20 - 32/3*a^18 - 200/3*a^16 + 304/3*a^14 - 212/3*a^12 + 88*a^10 - 88*a^8 - 128/3*a^6 + 428/3*a^4 - 280/3*a^2 + 20)*u^3 + (20/3*a^36 - 40/3*a^34 - 4*a^32 + 40/3*a^30 - 4/3*a^28 + 92/3*a^26 - 320/3*a^24 - 44/3*a^22 + 320*a^20 - 1312/3*a^18 + 736/3*a^16 + 16/3*a^14 - 68/3*a^12 - 632/3*a^10 + 264*a^8 + 128*a^6 - 428*a^4 + 280*a^2 - 60)*u^2 + (1/3*a^50 + 14/3*a^42 + 28/3*a^40 - 64/3*a^38 - 512/3*a^36 + 275/3*a^34 + 2776/3*a^32 - 3992/3*a^30 + 484/3*a^28 + 494/3*a^26 + 1616/3*a^24 + 1048/3*a^22 - 2168*a^20 + 7136/3*a^18 - 1880/3*a^16 - 2816/3*a^14 + 2248/3*a^12 + 784/3*a^10 - 528*a^8 - 256*a^6 + 856*a^4 - 560*a^2 + 120)*u + 1/2*a^52 - 1/3*a^50 + 2/3*a^48 - 8*a^46 - 47/3*a^44 + 227/3*a^42 - 77/2*a^40 - 290/3*a^38 + 1595/6*a^36 - 1159/3*a^34 + 224/3*a^32 + 250*a^30 + 19/3*a^28 + 237*a^26 - 1619/2*a^24 - 254/3*a^22 + 4763/3*a^20 - 1456*a^18 - 520/3*a^16 + 3488/3*a^14 - 2272/3*a^12 + 160*a^10)*v + (1/6*a^48 - 4*a^42 - 5/3*a^40 + 100/3*a^38 - 64*a^36 + 286/3*a^34 - 833/6*a^32 + 362/3*a^30 - 62/3*a^28 - 476/3*a^26 + 893/3*a^24 + 20/3*a^22 - 1562/3*a^20 + 1408/3*a^18 + 328/3*a^16 - 1280/3*a^14 + 784/3*a^12 - 160/3*a^10)*u^2 + (1/2*a^52 - 1/3*a^50 + 2/3*a^48 - 8*a^46 - 47/3*a^44 + 227/3*a^42 - 77/2*a^40 - 290/3*a^38 + 1595/6*a^36 - 1159/3*a^34 + 224/3*a^32 + 250*a^30 + 19/3*a^28 + 237*a^26 - 1619/2*a^24 - 254/3*a^22 + 4763/3*a^20 - 1456*a^18 - 520/3*a^16 + 3488/3*a^14 - 2272/3*a^12 + 160*a^10)*u + 1/2*a^54 - 9/8*a^52 + 7/3*a^50 - 29/2*a^48 - 39/4*a^46 + 2227/12*a^44 - 5969/18*a^42 + 765/8*a^40 + 865/4*a^38 - 20897/72*a^36 + 3379/6*a^34 - 6813/8*a^32 + 1213/2*a^30 - 2981/24*a^28 - 440*a^26 + 18617/24*a^24 - 2267/12*a^22 - 2627/4*a^20 + 1468/3*a^18 + 9643/24*a^16 - 8995/12*a^14 + 10319/24*a^12 - 688/3*a^10 + 176*a^8 + 256/3*a^6 - 856/3*a^4 + 560/3*a^2 - 40)*t^4 + (((-1/6*a^48 + 2*a^42 + 7/10*a^40 - 28/3*a^38 + 46/3*a^36 - 32*a^34 + 158/3*a^32 - 134/3*a^30 + 16/3*a^28 + 28*a^26 + 482/9*a^24 - 638/3*a^22 + 608/3*a^20 + 128/3*a^18 - 776/3*a^16 + 688/3*a^14 - 1468/9*a^12 + 264*a^10 - 652/3*a^8 - 272/3*a^6 + 728/3*a^4 - 400/3*a^2 + 1096/45)*u)*v^2 + ((-1/6*a^48 + 2*a^42 + 7/10*a^40 - 28/3*a^38 + 46/3*a^36 - 32*a^34 + 158/3*a^32 - 134/3*a^30 + 16/3*a^28 + 28*a^26 + 482/9*a^24 - 638/3*a^22 + 608/3*a^20 + 128/3*a^18 - 776/3*a^16 + 688/3*a^14 - 1468/9*a^12 + 264*a^10 - 652/3*a^8 - 272/3*a^6 + 728/3*a^4 - 400/3*a^2 + 1096/45)*u^2 + (11/5*a^50 + 2/3*a^48 - 4*a^44 - 125/3*a^42 + 6/5*a^40 + 112/3*a^38 + 584/3*a^36 + 232*a^34 - 3980/3*a^32 + 1100*a^30 - 40/3*a^28 - 80*a^26 - 3440/9*a^24 + 1496/3*a^22 + 496/3*a^20 - 3328/3*a^18 + 3248/3*a^16 - 1024/3*a^14 + 2416/9*a^12 - 4896/5*a^10 + 2608/3*a^8 + 1088/3*a^6 - 2912/3*a^4 + 1600/3*a^2 - 4384/45)*u + 1/30*a^60 - 3/10*a^52 - 9/5*a^50 - 13/18*a^48 + 4*a^46 + 59/6*a^44 - 56/9*a^42 + 9*a^40 + 140/3*a^38 - 4453/18*a^36 + 25/3*a^34 + 3839/6*a^32 - 4544/9*a^30 - 42*a^28 - 170/3*a^26 + 1363/9*a^24 + 860/3*a^22 - 11368/15*a^20 + 5632/9*a^18 + 64*a^16 - 448*a^14 + 1344/5*a^12 - 256/5*a^10)*v + (1/30*a^60 - 3/10*a^52 - 9/5*a^50 - 13/18*a^48 + 4*a^46 + 59/6*a^44 - 56/9*a^42 + 9*a^40 + 140/3*a^38 - 4453/18*a^36 + 25/3*a^34 + 3839/6*a^32 - 4544/9*a^30 - 42*a^28 - 170/3*a^26 + 1363/9*a^24 + 860/3*a^22 - 11368/15*a^20 + 5632/9*a^18 + 64*a^16 - 448*a^14 + 1344/5*a^12 - 256/5*a^10)*u + 1/20*a^62 - 1/15*a^60 + 1/8*a^56 - 17/12*a^54 + 19/24*a^52 - 23/5*a^50 + 353/18*a^48 + 137/4*a^46 - 508/3*a^44 + 1555/9*a^42 - 259/15*a^40 - 1729/12*a^38 + 291*a^36 - 664/3*a^34 - 63/4*a^32 - 1121/36*a^30 + 272/3*a^28 + 305/2*a^26 - 2653/72*a^24 - 13041/20*a^22 + 114223/120*a^20 - 4045/18*a^18 - 8411/12*a^16 + 4697/6*a^14 - 5767/12*a^12 + 2768/5*a^10 - 1304/3*a^8 - 544/3*a^6 + 1456/3*a^4 - 800/3*a^2 + 2192/45)*t^3 + (((1/9*a^60 - 2/5*a^52 - 22/15*a^50 - 5/9*a^48 + 5/3*a^44 + 190/9*a^42 - 79/15*a^40 - 56/3*a^38 - 196/9*a^36 - 96*a^34 + 868/3*a^32 - 620/3*a^30 + 56/3*a^28 + 16*a^26 + 1528/9*a^24 - 1336/3*a^22 + 1072/3*a^20 + 320/3*a^18 - 1184/3*a^16 + 832/3*a^14 - 2864/15*a^12 + 352*a^10 - 272*a^8 - 320/3*a^6 + 256*a^4 - 128*a^2 + 64/3)*u)*v + 1/360*a^72 - 1/10*a^62 - 1/9*a^60 - 1/4*a^56 + 4/3*a^54 + 1/3*a^52 + 112/15*a^50 - 203/18*a^48 - 59/2*a^46 + 297/4*a^44 - 196/3*a^42 + 367/120*a^40 + 24*a^38 + 295/18*a^36 + 105*a^34 - 1199/4*a^32 + 1301/6*a^30 - 131/6*a^28 - 59/3*a^26 - 6061/36*a^24 + 13501/30*a^22 - 4339/12*a^20 - 325/3*a^18 + 1196/3*a^16 - 838/3*a^14 + 2869/15*a^12 - 352*a^10 + 272*a^8 + 320/3*a^6 - 256*a^4 + 128*a^2 - 64/3)*t^2 + ((-1/360*a^72 + 1/20*a^62 + 1/8*a^56 - 5/12*a^54 - 2*a^50 + 17/6*a^48 + 15/2*a^46 - 15*a^44 + 7*a^42)*y)*t)*w^6 + (((2/3*a^12*u*v^4 + ((16/3*a^18 - 8/3*a^12)*u + 8/3*a^20 - 16/3*a^18)*v^3 + (4*a^20*u^2 + (8*a^24 + 8*a^22 - 12*a^20 - 16*a^18 + 4*a^12)*u + 6*a^26 - 10*a^24 - 8*a^22 + 16*a^18)*v^2 + (2/3*a^12*u^4 + (16/3*a^18 - 8/3*a^12)*u^3 + (8*a^24 + 8*a^22 - 12*a^20 - 16*a^18 + 4*a^12)*u^2 + (12*a^28 - 4*a^26 - 64*a^22 + 32*a^20 + 112/3*a^18 - 16/3*a^12)*u + 8*a^30 - 18*a^28 - 8*a^26 + 38/3*a^24 + 32*a^22 - 8*a^20 - 64/3*a^18)*v + (8/3*a^20 - 16/3*a^18)*u^3 + (6*a^26 - 10*a^24 - 8*a^22 + 16*a^18)*u^2 + (8*a^30 - 18*a^28 - 8*a^26 + 38/3*a^24 + 32*a^22 - 8*a^20 - 64/3*a^18)*u + 5*a^32 - 22*a^30 + 99/4*a^28 + 9*a^26 - 64/3*a^24 - 17/12*a^20 + 16/3*a^18 + 4/3*a^12)*y - 2/3*a^12*u*v^5 + ((2/3*a^20 - 2/3*a^12)*u^2 + (-4/3*a^20 - 16/3*a^18 + 10/3*a^12)*u - 8/3*a^20 + 16/3*a^18)*v^4 + ((8/3*a^26 - 20/3*a^20 - 16/3*a^18 + 8/3*a^12)*u^2 + (-8*a^26 - 8*a^24 - 8*a^22 + 44/3*a^20 + 80/3*a^18 - 20/3*a^12)*u - 8/3*a^28 - 2/3*a^26 + 10*a^24 + 8*a^22 + 8/3*a^20 - 64/3*a^18)*v^3 + ((2/3*a^20 - 2/3*a^12)*u^4 + (8/3*a^26 - 20/3*a^20 - 16/3*a^18 + 8/3*a^12)*u^3 + (16*a^30 - 12*a^28 - 16*a^26 - 16*a^24 - 16*a^22 + 36*a^20 + 32*a^18 - 8*a^12)*u^2 + (8*a^32 - 48*a^30 + 12*a^28 + 98/3*a^26 + 18*a^24 + 80*a^22 - 164/3*a^20 - 208/3*a^18 + 28/3*a^12)*u - 6*a^34 + 10*a^32 + 18*a^28 - 2*a^26 - 68/3*a^24 - 40*a^22 + 8*a^20 + 112/3*a^18)*v^2 + (-2/3*a^12*u^5 + (-4/3*a^20 - 16/3*a^18 + 10/3*a^12)*u^4 + (-8*a^26 - 8*a^24 - 8*a^22 + 44/3*a^20 + 80/3*a^18 - 20/3*a^12)*u^3 + (8*a^32 - 48*a^30 + 12*a^28 + 98/3*a^26 + 18*a^24 + 80*a^22 - 164/3*a^20 - 208/3*a^18 + 28/3*a^12)*u^2 + (-12*a^36 + 4*a^34 - 32*a^32 + 112*a^30 - 140/3*a^26 - 76/3*a^24 - 128*a^22 + 176/3*a^20 + 80*a^18 - 16/3*a^12)*u - 8*a^38 + 18*a^36 + 8*a^34 - 53/3*a^32 - 2*a^30 - 139/4*a^28 + 13/3*a^26 + 34*a^24 + 32*a^22 - 79/12*a^20 - 80/3*a^18 - 4/3*a^12)*v + (-8/3*a^20 + 16/3*a^18)*u^4 + (-8/3*a^28 - 2/3*a^26 + 10*a^24 + 8*a^22 + 8/3*a^20 - 64/3*a^18)*u^3 + (-6*a^34 + 10*a^32 + 18*a^28 - 2*a^26 - 68/3*a^24 - 40*a^22 + 8*a^20 + 112/3*a^18)*u^2 + (-8*a^38 + 18*a^36 + 8*a^34 - 53/3*a^32 - 2*a^30 - 139/4*a^28 + 13/3*a^26 + 34*a^24 + 32*a^22 - 79/12*a^20 - 80/3*a^18 - 4/3*a^12)*u - 5*a^40 + 22*a^38 - 99/4*a^36 - 9*a^34 + 79/3*a^32 - 22*a^30 + 157/6*a^28 + 11/3*a^26 - 64/3*a^24 - 11/4*a^20 + 16/3*a^18 + 4/3*a^12)*t^6 + ((((-8/3*a^20 + 4*a^16 - 8*a^14 + 4*a^12)*u)*v^3 + ((-8*a^24 - 16*a^20 + 16*a^18 - 12*a^16 + 24*a^14 - 12*a^12)*u + 2*a^28 - 12*a^26 + 18*a^24 - 16*a^22 + 32*a^20 - 16*a^18)*v^2 + ((-8/3*a^20 + 4*a^16 - 8*a^14 + 4*a^12)*u^3 + (-8*a^24 - 16*a^20 + 16*a^18 - 12*a^16 + 24*a^14 - 12*a^12)*u^2 + (4*a^30 - 8*a^28 - 12*a^26 - 40*a^24 + 48*a^22 + 104*a^20 - 80*a^18 + 24*a^16 - 48*a^14 + 24*a^12)*u + 6*a^32 - 22*a^30 - 2*a^28 + 58*a^26 - 16*a^24 + 8*a^22 - 88*a^20 + 48*a^18)*v + (2*a^28 - 12*a^26 + 18*a^24 - 16*a^22 + 32*a^20 - 16*a^18)*u^2 + (6*a^32 - 22*a^30 - 2*a^28 + 58*a^26 - 16*a^24 + 8*a^22 - 88*a^20 + 48*a^18)*u + 6*a^34 - 67/2*a^32 + 52*a^30 + 25/2*a^28 - 251/3*a^26 + 105/2*a^24 - 33*a^22 + 275/6*a^20 - 16*a^18 - 8*a^16 + 16*a^14 - 8*a^12)*y + ((10/3*a^20 - 4*a^16 + 8*a^14 - 4*a^12)*u)*v^4 + ((-8/3*a^26 + 4*a^24 - 8*a^22 + 20/3*a^20 - 4*a^16 + 8*a^14 - 4*a^12)*u^2 + (32/3*a^26 + 16*a^22 + 8/3*a^20 - 16*a^18 + 16*a^16 - 32*a^14 + 16*a^12)*u + 2/3*a^28 + 20/3*a^26 - 18*a^24 + 16*a^22 - 32*a^20 + 16*a^18)*v^3 + ((-8/3*a^26 + 4*a^24 - 8*a^22 + 20/3*a^20 - 4*a^16 + 8*a^14 - 4*a^12)*u^3 + (8*a^32 - 16*a^30 - 20*a^28 + 32*a^26 - 8*a^24 + 48*a^22 + 8*a^20 - 32*a^18 + 24*a^16 - 48*a^14 + 24*a^12)*u^2 + (24*a^34 - 80*a^32 + 76*a^30 + 66*a^28 - 64*a^26 + 50*a^24 - 104*a^22 - 112*a^20 + 112*a^18 - 36*a^16 + 72*a^14 - 36*a^12)*u - 2*a^36 + 10*a^34 - 10*a^32 + 14*a^30 - 28*a^28 - 38*a^26 + 34*a^24 - 24*a^22 + 120*a^20 - 64*a^18)*v^2 + ((10/3*a^20 - 4*a^16 + 8*a^14 - 4*a^12)*u^4 + (32/3*a^26 + 16*a^22 + 8/3*a^20 - 16*a^18 + 16*a^16 - 32*a^14 + 16*a^12)*u^3 + (24*a^34 - 80*a^32 + 76*a^30 + 66*a^28 - 64*a^26 + 50*a^24 - 104*a^22 - 112*a^20 + 112*a^18 - 36*a^16 + 72*a^14 - 36*a^12)*u^2 + (-4*a^38 + 4*a^36 - 56*a^34 + 292*a^32 - 240*a^30 - 164*a^28 + 208/3*a^26 - 56*a^24 + 128*a^22 + 680/3*a^20 - 176*a^18 + 24*a^16 - 48*a^14 + 24*a^12)*u - 6*a^40 + 18*a^38 + 16*a^36 - 60*a^34 - 71/6*a^32 - 2*a^30 + 131/2*a^28 + 217/3*a^26 - 137/2*a^24 + 41*a^22 - 803/6*a^20 + 64*a^18 + 8*a^16 - 16*a^14 + 8*a^12)*v + (2/3*a^28 + 20/3*a^26 - 18*a^24 + 16*a^22 - 32*a^20 + 16*a^18)*u^3 + (-2*a^36 + 10*a^34 - 10*a^32 + 14*a^30 - 28*a^28 - 38*a^26 + 34*a^24 - 24*a^22 + 120*a^20 - 64*a^18)*u^2 + (-6*a^40 + 18*a^38 + 16*a^36 - 60*a^34 - 71/6*a^32 - 2*a^30 + 131/2*a^28 + 217/3*a^26 - 137/2*a^24 + 41*a^22 - 803/6*a^20 + 64*a^18 + 8*a^16 - 16*a^14 + 8*a^12)*u - 6*a^42 + 61/2*a^40 - 36*a^38 - 167/4*a^36 + 106*a^34 - 224/3*a^32 + 69*a^30 - 353/12*a^28 - 203/3*a^26 + 121/2*a^24 - 49*a^22 + 331/6*a^20 - 16*a^18 - 8*a^16 + 16*a^14 - 8*a^12)*t^5 + ((((-2*a^28 + 6*a^26 - 16/3*a^24 + 16*a^22 - 22*a^20 + 8*a^18 + 20*a^16 - 24*a^14 + 22/3*a^12)*u)*v^2 + ((-2*a^28 + 6*a^26 - 16/3*a^24 + 16*a^22 - 22*a^20 + 8*a^18 + 20*a^16 - 24*a^14 + 22/3*a^12)*u^2 + (4*a^30 + 8*a^28 - 72*a^26 + 424/3*a^24 - 64*a^22 + 4*a^20 - 80*a^16 + 96*a^14 - 88/3*a^12)*u + 2/3*a^36 - 22/3*a^32 + 16/3*a^30 + 10*a^28 + 28*a^26 - 206/3*a^24 - 16*a^22 + 184/3*a^20 - 64/3*a^18)*v + (2/3*a^36 - 22/3*a^32 + 16/3*a^30 + 10*a^28 + 28*a^26 - 206/3*a^24 - 16*a^22 + 184/3*a^20 - 64/3*a^18)*u + a^38 - 7/12*a^36 - 18*a^34 + 281/6*a^32 - 65/3*a^30 - 63/2*a^28 + 7*a^26 + 25/6*a^24 + 67*a^22 - 1003/12*a^20 + 80/3*a^18 + 40*a^16 - 48*a^14 + 44/3*a^12)*y + ((2*a^28 - 26/3*a^26 + 28/3*a^24 - 24*a^22 + 26*a^20 - 8*a^18 - 20*a^16 + 24*a^14 - 22/3*a^12)*u)*v^3 + ((4/3*a^32 - 4*a^30 + 8*a^28 - 12*a^26 + 152/3*a^24 - 80*a^22 + 176/3*a^20 - 16*a^18 - 40*a^16 + 48*a^14 - 44/3*a^12)*u^2 + (16/3*a^38 + 4*a^36 - 68*a^34 + 316/3*a^32 - 28*a^30 - 70*a^28 + 338/3*a^26 - 716/3*a^24 + 200*a^22 - 202/3*a^20 + 8*a^18 + 100*a^16 - 120*a^14 + 110/3*a^12)*u + 4/3*a^36 - 4*a^34 + 4/3*a^32 - 16/3*a^30 + 22*a^28 - 44*a^26 + 206/3*a^24 + 16*a^22 - 184/3*a^20 + 64/3*a^18)*v^2 + ((2*a^28 - 26/3*a^26 + 28/3*a^24 - 24*a^22 + 26*a^20 - 8*a^18 - 20*a^16 + 24*a^14 - 22/3*a^12)*u^3 + (16/3*a^38 + 4*a^36 - 68*a^34 + 316/3*a^32 - 28*a^30 - 70*a^28 + 338/3*a^26 - 716/3*a^24 + 200*a^22 - 202/3*a^20 + 8*a^18 + 100*a^16 - 120*a^14 + 110/3*a^12)*u^2 + (-68/3*a^38 - 40/3*a^36 + 344*a^34 - 632*a^32 + 556/3*a^30 + 264*a^28 - 832/3*a^26 + 1388/3*a^24 - 272*a^22 - 36*a^20 + 128/3*a^18 - 80*a^16 + 96*a^14 - 88/3*a^12)*u - 2/3*a^44 + 16/3*a^40 + 23/3*a^38 - 75/4*a^36 - 100*a^34 + 413/2*a^32 - 29*a^30 - 647/6*a^28 + 271/3*a^26 - 437/6*a^24 - 83*a^22 + 1739/12*a^20 - 48*a^18 - 40*a^16 + 48*a^14 - 44/3*a^12)*v + (4/3*a^36 - 4*a^34 + 4/3*a^32 - 16/3*a^30 + 22*a^28 - 44*a^26 + 206/3*a^24 + 16*a^22 - 184/3*a^20 + 64/3*a^18)*u^2 + (-2/3*a^44 + 16/3*a^40 + 23/3*a^38 - 75/4*a^36 - 100*a^34 + 413/2*a^32 - 29*a^30 - 647/6*a^28 + 271/3*a^26 - 437/6*a^24 - 83*a^22 + 1739/12*a^20 - 48*a^18 - 40*a^16 + 48*a^14 - 44/3*a^12)*u - a^46 + 7/12*a^44 + 15*a^42 - 85/3*a^40 - 46/3*a^38 + 569/12*a^36 + 21*a^34 - 55/2*a^32 - 161/3*a^30 + 535/12*a^28 - 43/3*a^26 - 263/6*a^24 + 131*a^22 - 425/4*a^20 + 80/3*a^18 + 40*a^16 - 48*a^14 + 44/3*a^12)*t^4 + ((((8/3*a^32 - 20/3*a^30 - 4*a^28 + 28*a^26 - 128/3*a^24 + 32*a^22 - 116/3*a^20 + 16*a^18 + 32*a^16 - 32*a^14 + 8*a^12)*u)*v + 1/12*a^44 - 2*a^38 - 3/2*a^36 + 18*a^34 - 70/3*a^32 + 23/3*a^30 + 25/4*a^28 - 85/3*a^26 + 134/3*a^24 - 34*a^22 + 235/6*a^20 - 16*a^18 - 32*a^16 + 32*a^14 - 8*a^12)*y + ((2/3*a^44 - 26/3*a^38 - 6*a^36 + 56*a^34 - 160/3*a^32 + 8/3*a^30 + 8*a^28 - 76/3*a^26 + 236/3*a^24 - 72*a^22 + 50*a^20 - 16*a^18 - 32*a^16 + 32*a^14 - 8*a^12)*u)*v^2 + ((2/3*a^44 - 26/3*a^38 - 6*a^36 + 56*a^34 - 160/3*a^32 + 8/3*a^30 + 8*a^28 - 76/3*a^26 + 236/3*a^24 - 72*a^22 + 50*a^20 - 16*a^18 - 32*a^16 + 32*a^14 - 8*a^12)*u^2 + (-8/3*a^44 + 152/3*a^38 + 36*a^36 - 428*a^34 + 1516/3*a^32 - 116/3*a^30 - 104*a^28 + 148/3*a^26 - 560/3*a^24 + 192*a^22 - 84*a^20 + 16*a^18 + 32*a^16 - 32*a^14 + 8*a^12)*u + 7/12*a^44 + 2/3*a^40 - 50/3*a^38 - 25/2*a^36 + 146*a^34 - 520/3*a^32 + 25/3*a^30 + 661/12*a^28 + 7*a^26 - 134/3*a^24 + 34*a^22 - 235/6*a^20 + 16*a^18 + 32*a^16 - 32*a^14 + 8*a^12)*v + (7/12*a^44 + 2/3*a^40 - 50/3*a^38 - 25/2*a^36 + 146*a^34 - 520/3*a^32 + 25/3*a^30 + 661/12*a^28 + 7*a^26 - 134/3*a^24 + 34*a^22 - 235/6*a^20 + 16*a^18 + 32*a^16 - 32*a^14 + 8*a^12)*u - 1/12*a^52 + 5/3*a^46 + 2*a^44 - 12*a^42 + 13/6*a^40 + 10*a^38 + 39/4*a^36 - 65/2*a^32 + 62/3*a^30 - 155/6*a^28 - 37/3*a^26 + 350/3*a^24 - 114*a^22 + 371/6*a^20 - 16*a^18 - 32*a^16 + 32*a^14 - 8*a^12)*t^3 + ((-1/12*a^44 + a^38 + 3/4*a^36 - 6*a^34 + 5*a^32)*y^2 + ((3/2*a^44 - 58/3*a^38 - 27/2*a^36 + 124*a^34 - 338/3*a^32 - 4*a^30 + 4*a^28 + 16/3*a^26 + 32*a^24 - 32*a^22 + 8*a^20)*u)*v + 1/12*a^52 - 2/3*a^46 - 5/2*a^44 + 3*a^42 + 2/3*a^40 + 55/3*a^38 + 29/4*a^36 - 112*a^34 + 310/3*a^32 + 2*a^30 - 7/2*a^28 - 16/3*a^26 - 32*a^24 + 32*a^22 - 8*a^20)*t^2 + ((-1/12*a^44 + a^38 + 3/4*a^36 - 6*a^34 + 5*a^32)*y)*t)*w^5 + (((a^20*u*v^2 + (a^20*u^2 + (6*a^22 - 4*a^20)*u + 2*a^24 - 4*a^22)*v + (2*a^24 - 4*a^22)*u + 3/2*a^26 - 9/2*a^24 + 2*a^22 + 2*a^20)*y^2 + (-2*a^20*u*v^3 + ((4*a^24 - 4*a^20)*u^2 + (4*a^26 - 8*a^24 - 12*a^22 + 10*a^20)*u - 4*a^24 + 8*a^22)*v^2 + (-2*a^20*u^3 + (4*a^26 - 8*a^24 - 12*a^22 + 10*a^20)*u^2 + (-20*a^26 + 8*a^24 + 28*a^22 - 8*a^20)*u - 4*a^28 + 5*a^26 + 13*a^24 - 12*a^22 - 4*a^20)*v + (-4*a^24 + 8*a^22)*u^2 + (-4*a^28 + 5*a^26 + 13*a^24 - 12*a^22 - 4*a^20)*u - 3*a^30 + 9*a^28 - a^26 - 13*a^24 + 4*a^22 + 4*a^20)*y + a^20*u*v^4 + ((a^28 - 4*a^24 + 3*a^20)*u^2 + (2*a^30 - 2*a^28 - 4*a^26 + 8*a^24 + 6*a^22 - 6*a^20)*u + 2*a^24 - 4*a^22)*v^3 + ((a^28 - 4*a^24 + 3*a^20)*u^3 + (-6*a^28 - 8*a^26 + 20*a^24 + 12*a^22 - 12*a^20)*u^2 + (-6*a^30 + 9*a^28 + 24*a^26 - 18*a^24 - 24*a^22 + 9*a^20)*u + 4*a^28 - 13/2*a^26 - 17/2*a^24 + 10*a^22 + 2*a^20)*v^2 + (a^20*u^4 + (2*a^30 - 2*a^28 - 4*a^26 + 8*a^24 + 6*a^22 - 6*a^20)*u^3 + (-6*a^30 + 9*a^28 + 24*a^26 - 18*a^24 - 24*a^22 + 9*a^20)*u^2 + (14*a^30 - 4*a^28 - 33*a^26 - a^24 + 26*a^22)*u + 2*a^32 - a^30 - 13*a^28 + 9*a^26 + 15*a^24 - 8*a^22 - 4*a^20)*v + (2*a^24 - 4*a^22)*u^3 + (4*a^28 - 13/2*a^26 - 17/2*a^24 + 10*a^22 + 2*a^20)*u^2 + (2*a^32 - a^30 - 13*a^28 + 9*a^26 + 15*a^24 - 8*a^22 - 4*a^20)*u + 3/2*a^34 - 9/2*a^32 - a^30 + 11*a^28 - 5/2*a^26 - 17/2*a^24 + 2*a^22 + 2*a^20)*t^5 + ((((-2*a^24 - 4*a^22 + 2*a^20)*u)*v + 1/2*a^28 - 3*a^26 + 5/2*a^24 + 4*a^22 - 2*a^20)*y^2 + (((2*a^28 - 8*a^26 + 9*a^24 + 6*a^22 - 3*a^20)*u)*v^2 + ((2*a^28 - 8*a^26 + 9*a^24 + 6*a^22 - 3*a^20)*u^2 + (-8*a^28 + 50*a^26 - 36*a^24 + 6*a^22)*u + 5*a^28 - 10*a^26 + 5*a^24 - 12*a^22 + 4*a^20)*v + (5*a^28 - 10*a^26 + 5*a^24 - 12*a^22 + 4*a^20)*u - a^32 + 13/2*a^30 - 17/2*a^28 - 11/2*a^26 + 13/2*a^24 + 6*a^22 - 2*a^20)*y + ((a^32 - 2*a^30 - 3*a^28 + 8*a^26 - 7*a^24 - 2*a^22 + a^20)*u)*v^3 + ((2*a^32 - 4*a^30 + 2*a^28 + 16*a^26 - 14*a^24 - 4*a^22 + 2*a^20)*u^2 + (2*a^34 - 11*a^32 + 20*a^30 + a^28 - 58*a^26 + 41*a^24 - 8*a^22 + a^20)*u - 11/2*a^28 + 13*a^26 - 15/2*a^24 + 8*a^22 - 2*a^20)*v^2 + ((a^32 - 2*a^30 - 3*a^28 + 8*a^26 - 7*a^24 - 2*a^22 + a^20)*u^3 + (2*a^34 - 11*a^32 + 20*a^30 + a^28 - 58*a^26 + 41*a^24 - 8*a^22 + a^20)*u^2 + (-10*a^34 + 34*a^32 - 58*a^30 - a^28 + 76*a^26 - 49*a^24 + 26*a^22 - 6*a^20)*u - 2*a^36 + 8*a^34 - 13*a^32 + 11/2*a^30 + 29/2*a^28 - 21/2*a^26 + 7/2*a^24 - 10*a^22 + 2*a^20)*v + (-11/2*a^28 + 13*a^26 - 15/2*a^24 + 8*a^22 - 2*a^20)*u^2 + (-2*a^36 + 8*a^34 - 13*a^32 + 11/2*a^30 + 29/2*a^28 - 21/2*a^26 + 7/2*a^24 - 10*a^22 + 2*a^20)*u - 3/2*a^38 + 8*a^36 - 33/2*a^34 + 13*a^32 + 13/2*a^30 - 13*a^28 - 5/2*a^26 + 4*a^24 + 2*a^22)*t^4 + ((-1/2*a^28 + 3/2*a^26)*y^3 + (((-a^30 + 8*a^28 - 25*a^26 + 12*a^24 - 8*a^22 + 2*a^20)*u)*v + 3/2*a^32 - 6*a^30 + 2*a^28 + 12*a^26 - 23/2*a^24 + 8*a^22 - 2*a^20)*y + ((a^36 - 4*a^34 + 7*a^32 - 5*a^30 - 15/2*a^28 + 41/2*a^26 - 12*a^24 + 8*a^22 - 2*a^20)*u)*v^2 + ((a^36 - 4*a^34 + 7*a^32 - 5*a^30 - 15/2*a^28 + 41/2*a^26 - 12*a^24 + 8*a^22 - 2*a^20)*u^2 + (-4*a^36 + 22*a^34 - 40*a^32 + 35*a^30 + 9*a^28 - 16*a^26 + 12*a^24 - 8*a^22 + 2*a^20)*u + 2*a^36 - 8*a^34 + 21/2*a^32 - 2*a^30 - 2*a^28 - 12*a^26 + 23/2*a^24 - 8*a^22 + 2*a^20)*v + (2*a^36 - 8*a^34 + 21/2*a^32 - 2*a^30 - 2*a^28 - 12*a^26 + 23/2*a^24 - 8*a^22 + 2*a^20)*u - 1/2*a^40 + 7/2*a^38 - 19/2*a^36 + 12*a^34 - 3*a^32 - 8*a^30 + 12*a^26 - 23/2*a^24 + 8*a^22 - 2*a^20)*t^3 + ((-1/2*a^32 + 5/2*a^30 - 9/2*a^28 + 9/2*a^26)*y^2 + ((2*a^36 - 8*a^34 + 11*a^32 - 9*a^30 + 5*a^28 - 9*a^26)*u)*v + 1/2*a^40 - 2*a^38 + 3/2*a^36 + 3*a^34 - 5*a^32 + 4*a^28)*t^2 + ((-1/2*a^32 + 5/2*a^30 - 4*a^28 + 3*a^26)*y)*t)*w^4 + ((1/6*a^24*y^4 + (((a^26 - 2*a^24)*u)*v)*y^2 + (((1/2*a^28 - 2*a^26 + 2*a^24)*u)*v^2 + ((1/2*a^28 - 2*a^26 + 2*a^24)*u^2 + (a^30 - 2*a^28 + 2*a^26)*u)*v)*y + ((-1/2*a^28 + a^26 - 2/3*a^24)*u)*v^3 + ((-a^28 + 2*a^26 - a^24)*u^2 + (-a^30 + 5/2*a^28 - 2*a^26)*u)*v^2 + ((-1/2*a^28 + a^26 - 2/3*a^24)*u^3 + (-a^30 + 5/2*a^28 - 2*a^26)*u^2 + (a^30 - 2*a^28 + a^26)*u)*v)*t^4 + ((1/2*a^28 - a^26 + a^24)*y^3 + (((a^32 - 2*a^30 - 2*a^28 + 7*a^26 - 6*a^24)*u)*v)*y + ((-a^32 + 2*a^30 + 1/2*a^28 - 4*a^26 + 3*a^24)*u)*v^2 + ((-a^32 + 2*a^30 + 1/2*a^28 - 4*a^26 + 3*a^24)*u^2 + (a^32 - 2*a^30 + a^28 + a^26)*u)*v)*t^3 + ((1/6*a^36 - a^32 + a^30 + 3/2*a^28 - 3*a^26 + 11/6*a^24)*y^2 + ((-1/3*a^36 + 2*a^32 - 2*a^30 - 3*a^28 + 6*a^26 - 11/3*a^24)*u)*v)*t^2 + ((1/6*a^36 - a^32 + a^30 + a^28 - 2*a^26 + a^24)*y)*t)*w^3]"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs (we take a parameter a = sqrt{2}), with\n",
    "# the marking variable t for strongly connected components,\n",
    "# the marking variable u for purely source-like components,\n",
    "# the marking variable v for purely sink-like components,\n",
    "# the marking variable y for isolated components\n",
    "QDTALLzero = [qdtall[i] * a^(i*(i-1)/2) for i in srange(N)]\n",
    "QDTALLzero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[t*y, 0, 0, 0, 0, 0, 0],\n",
       " [0, 4*t^2*u*v - 4*t*y, 0, 0, 0, 0, 0],\n",
       " [0,\n",
       "  -4*(2*u*v - y^2)*t^2 + 4*t*y,\n",
       "  2*(u*v^2 + (u^2 + 8*u)*v)*t^3 - 28*t^2*u*v + 8*t*y,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  32*(u*v*y - 3*u*v)*t^3 + 32*(4*u*v - y^2)*t^2 - 32*t*y,\n",
       "  4/3*(u*v^3 + 21*u*v^2 + (u^3 + 21*u^2 + 78*u)*v)*t^4 - 4*(7*u*v^2 + (7*u^2 + 44*u)*v)*t^3 + 112*t^2*u*v - 128/3*t*y,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  -32*(u*v^2 + 2*u*v*y - y^3 + (u^2 - 4*u)*v)*t^3 - 96*(2*u*v - y^2)*t^2 + 64*t*y,\n",
       "  -32*(u*v^3 - (6*u^2 - 7*u)*v^2 + (u^3 + 7*u^2 + 40*u)*v - (u*v^2 + (u^2 + 8*u)*v)*y)*t^4 + 64*(3*u*v^2 - 7*u*v*y + 3*(u^2 + 11*u)*v)*t^3 - 128*(8*u*v - y^2)*t^2 + 128*t*y,\n",
       "  2/3*(u*v^4 + 60*u*v^3 + 6*(16*u^2 + 113*u)*v^2 + (u^4 + 60*u^3 + 678*u^2 + 1976*u)*v)*t^5 - 4/3*(29*u*v^3 + 489*u*v^2 + (29*u^3 + 489*u^2 + 1134*u)*v)*t^4 + 2*(89*u*v^2 + (89*u^2 + 892*u)*v)*t^3 - 248*t^2*u*v - 4096/3*t*y,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  512*(u*v^3 + u*v*y^2 - 2*(u^2 - u)*v^2 - 6*u*v*y + (u^3 + 2*u^2 + 5*u)*v)*t^4 + 512*(2*u*v^2 + 7*u*v*y - y^3 + (2*u^2 - 23*u)*v)*t^3 + 512*(14*u*v - 3*y^2)*t^2 - 1024*t*y,\n",
       "  -128/3*(u*v^4 - (7*u^2 - 36*u)*v^3 - (7*u^3 + 102*u^2 - 177*u)*v^2 + (u^4 + 36*u^3 + 177*u^2 + 786*u)*v - (u*v^3 + 21*u*v^2 + (u^3 + 21*u^2 + 78*u)*v)*y)*t^5 + 128/3*(29*u*v^3 - 3*(34*u^2 - 125*u)*v^2 + (29*u^3 + 375*u^2 + 876*u)*v - 3*(7*u*v^2 + (7*u^2 + 44*u)*v)*y)*t^4 - 512*(u*v^2 - 7*u*v*y + (u^2 + 65*u)*v)*t^3 + 4096/3*(11*u*v - y^2)*t^2 - 4096/3*t*y,\n",
       "  4/15*(u*v^5 + 155*u*v^4 + 10*(64*u^2 + 513*u)*v^3 + 10*(64*u^3 + 1152*u^2 + 4767*u)*v^2 + (u^5 + 155*u^4 + 5130*u^3 + 47670*u^2 + 122730*u)*v)*t^6 - 8/3*(15*u*v^4 + 740*u*v^3 + 6*(192*u^2 + 959*u)*v^2 + (15*u^4 + 740*u^3 + 5754*u^2 + 7016*u)*v)*t^5 + 4/3*(359*u*v^3 + 6699*u*v^2 + (359*u^3 + 6699*u^2 + 43146*u)*v)*t^4 - 8*(553*u*v^2 + (553*u^2 - 8196*u)*v)*t^3 + 271808/3*t^2*u*v - 3473408/15*t*y,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  -2048/3*(4*u*v^3 - y^4 + 6*(u^2 + 2*u)*v^2 - 24*u*v*y + 4*(u^3 + 3*u^2 - 3*u)*v)*t^4 - 4096*(3*u*v^2 - y^3 + 3*(u^2 - 2*u)*v)*t^3 - 18432*(2*u*v - y^2)*t^2 + 45056/3*t*y,\n",
       "  1024*(u*v^4 + (3*u^2 + 38*u)*v^3 + (3*u^3 - 68*u^2 + 33*u)*v^2 + (u*v^2 + (u^2 + 8*u)*v)*y^2 + (u^4 + 38*u^3 + 33*u^2 + 168*u)*v - 2*(u*v^3 - (6*u^2 - 7*u)*v^2 + (u^3 + 7*u^2 + 40*u)*v)*y)*t^5 - 1024*(15*u*v^3 + 14*u*v*y^2 - (98*u^2 - 107*u)*v^2 + (15*u^3 + 107*u^2 + 518*u)*v - (13*u*v^2 + (13*u^2 + 140*u)*v)*y)*t^4 + 2048*(21*u*v^2 - 35*u*v*y + 2*y^3 + 3*(7*u^2 + 67*u)*v)*t^3 - 4096*(22*u*v - 3*y^2)*t^2 + 8192*t*y,\n",
       "  -128/3*(u*v^5 - (15*u^2 - 91*u)*v^4 - 2*(146*u^2 - 1013*u)*v^3 - (15*u^4 + 292*u^3 + 3540*u^2 - 10258*u)*v^2 + (u^5 + 91*u^4 + 2026*u^3 + 10258*u^2 + 39752*u)*v - (u*v^4 + 60*u*v^3 + 6*(16*u^2 + 113*u)*v^2 + (u^4 + 60*u^3 + 678*u^2 + 1976*u)*v)*y)*t^6 + 256/3*(37*u*v^4 - (179*u^2 - 1356*u)*v^3 - (179*u^3 + 78*u^2 - 10341*u)*v^2 + (37*u^4 + 1356*u^3 + 10341*u^2 + 11026*u)*v - (29*u*v^3 + 489*u*v^2 + (29*u^3 + 489*u^2 + 1134*u)*v)*y)*t^5 - 128/3*(683*u*v^3 - 3*(622*u^2 + 573*u)*v^2 + (683*u^3 - 1719*u^2 + 53196*u)*v - 3*(89*u*v^2 + (89*u^2 + 892*u)*v)*y)*t^4 + 512/3*(4385*u*v^2 - 93*u*v*y + (4385*u^2 - 6893*u)*v)*t^3 + 8192/3*(563*u*v - 32*y^2)*t^2 - 262144/3*t*y,\n",
       "  4/45*(u*v^6 + 378*u*v^5 + 15*(256*u^2 + 2177*u)*v^4 + 20*(512*u^3 + 9984*u^2 + 44479*u)*v^3 + 15*(256*u^4 + 13312*u^3 + 156672*u^2 + 503809*u)*v^2 + (u^6 + 378*u^5 + 32655*u^4 + 889580*u^3 + 7557135*u^2 + 17990772*u)*v)*t^7 - 4/15*(123*u*v^5 + 16025*u*v^4 + 10*(6400*u^2 + 38779*u)*v^3 + 10*(6400*u^3 + 84480*u^2 + 194437*u)*v^2 + (123*u^5 + 16025*u^4 + 387790*u^3 + 1944370*u^2 - 81970*u)*v)*t^6 + 4/3*(795*u*v^4 + 37652*u*v^3 + 6*(7424*u^2 + 82523*u)*v^2 + (795*u^4 + 37652*u^3 + 495138*u^2 + 3038120*u)*v)*t^5 - 4/9*(57043*u*v^3 - 225657*u*v^2 + (57043*u^3 - 225657*u^2 - 21528078*u)*v)*t^4 - 136/3*(14569*u*v^2 + (14569*u^2 - 675748*u)*v)*t^3 + 2895728576/45*t^2*u*v - 4984930304/45*t*y]]"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QDTALLzeroMatrix = [[QDTALLzero[i][j].subs(a=sqrt(2)) for j in srange(N)] for i in srange(N)]\n",
    "QDTALLzeroMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 256, 0, 0, 0],\n",
       " [0, 0, 0, -3072, -5632, 0, 0],\n",
       " [0, 0, 0, 16384, 143360, 114176, 0]]"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs, particular example (Table 18)\n",
    "# the marking variable t=4 for all strongly connected components,\n",
    "# the marking variable u=1 for purely source-like components,\n",
    "# the marking variable v=1 for purely sink-like components,\n",
    "# the marking variable y=1 for components that are either source-like or sink-like\n",
    "QDTOMatrix = [[QDTALLzeroMatrix[i][j].diff(u).diff(v).diff(y).diff(t,4).subs(u=0,v=0,y=0,t=0)/24 for j in srange(N)] for i in srange(N)]\n",
    "QDTOMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.0",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
