{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:black; font-weight: bold; font-size:26px\">\n",
    "    Asymptotics of graphs and digraphs\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 Sections A.7-A.8, Tables 5-11"
   ]
  },
  {
   "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": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 15 # Our accuracy\n",
    "disp = 11 # How many terms we want to display\n",
    "disps = 5 # How many terms we want to display (short version)\n",
    "P.<a> = PolynomialRing(QQ)\n",
    "PP.<t> = PolynomialRing(P)\n",
    "PR.<w> = PolynomialRing(PP)\n",
    "R.<z> = PowerSeriesRing(PR,N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, a, a^3, a^6, a^10, a^15, a^21, a^28, a^36, a^45]"
      ]
     },
     "execution_count": 3,
     "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": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, 1/a, 1/a^3, 1/a^6, 1/a^10, 1/a^15, 1/a^21, 1/a^28, 1/a^36, 1/a^45]"
      ]
     },
     "execution_count": 29,
     "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(11)]\n",
    "PTTO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, 1/a^2, 1/a^6, 1/a^12, 1/a^20, 1/a^30, 1/a^42, 1/a^56, 1/a^72, 1/a^90]"
      ]
     },
     "execution_count": 24,
     "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(11)]\n",
    "HATSET"
   ]
  },
  {
   "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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, a^2, a^6, a^12, a^20, a^30, a^42, a^56, a^72, a^90]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of labeled graphs/tournaments\n",
    "# we take the 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": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0,\n",
       " 1,\n",
       " a^2 - 2,\n",
       " a^6 - 6*a^2 + 6,\n",
       " a^12 - 8*a^6 - 6*a^4 + 36*a^2 - 24,\n",
       " a^20 - 10*a^12 - 20*a^8 + 60*a^6 + 90*a^4 - 240*a^2 + 120,\n",
       " a^30 - 12*a^20 - 30*a^14 + 70*a^12 + 360*a^8 - 390*a^6 - 1080*a^4 + 1800*a^2 - 720,\n",
       " a^42 - 14*a^30 - 42*a^22 + 126*a^20 - 70*a^18 + 630*a^14 - 420*a^12 + 630*a^10 - 5040*a^8 + 1680*a^6 + 12600*a^4 - 15120*a^2 + 5040,\n",
       " a^56 - 16*a^42 - 56*a^32 + 168*a^30 - 112*a^26 - 70*a^24 + 1008*a^22 - 1344*a^20 + 1680*a^18 + 1260*a^16 - 8400*a^14 + 1680*a^12 - 20160*a^10 + 64680*a^8 + 10080*a^6 - 151200*a^4 + 141120*a^2 - 40320,\n",
       " a^72 - 18*a^56 - 72*a^44 + 216*a^42 - 168*a^36 + 1260*a^32 - 2016*a^30 + 3024*a^26 + 4158*a^24 - 18144*a^22 + 22680*a^20 - 28560*a^18 - 45360*a^16 + 90720*a^14 - 20160*a^12 + 453600*a^10 - 793800*a^8 - 483840*a^6 + 1905120*a^4 - 1451520*a^2 + 362880,\n",
       " a^90 - 20*a^72 - 90*a^58 + 270*a^56 - 240*a^48 + 2160*a^44 - 3300*a^42 - 252*a^40 + 5040*a^36 + 3780*a^34 - 22680*a^32 + 25200*a^30 + 15120*a^28 - 51030*a^26 - 115920*a^24 + 302400*a^22 - 483840*a^20 + 361200*a^18 + 982800*a^16 - 756000*a^14 + 1058400*a^12 - 8958600*a^10 + 9298800*a^8 + 11037600*a^6 - 25401600*a^4 + 16329600*a^2 - 3628800]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of irreducible labeled tournaments \n",
    "# we take the parameter a = sqrt{2}\n",
    "it = 1 - 1/g\n",
    "IT = [it[i] * i.factorial() for i in srange(disp)]\n",
    "IT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 1, 0, 2, 24, 544, 22320, 1677488, 236522496, 64026088576, 33832910196480]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The numbers of irreducible labeled tournaments\n",
    "ITsub = [it[i][0].subs(a=sqrt(2)) * i.factorial() for i in srange(disp)]\n",
    "ITsub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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,\n",
       " a^40 - 10*a^32 - 20*a^28 + 60*a^26 + 95*a^24 - 240*a^22 + 120*a^20 - 40*a^18 - 20*a^16 + 160*a^14 - 160*a^12 + 180*a^10 - 210*a^8 + 180*a^4 - 120*a^2 + 24,\n",
       " a^60 - 12*a^50 - 30*a^44 + 70*a^42 + 6*a^40 + 360*a^38 - 390*a^36 - 1080*a^34 + 1740*a^32 - 720*a^30 - 105*a^28 + 330*a^26 + 580*a^24 - 1560*a^22 + 750*a^20 + 600*a^18 - 1140*a^16 + 840*a^14 - 930*a^12 + 1980*a^10 - 1530*a^8 - 600*a^6 + 1440*a^4 - 720*a^2 + 120,\n",
       " a^84 - 14*a^72 - 42*a^64 + 126*a^62 - 63*a^60 + 630*a^56 - 420*a^54 + 630*a^52 - 5124*a^50 + 1680*a^48 + 12600*a^46 - 15309*a^44 + 5488*a^42 + 42*a^40 + 2520*a^38 - 2905*a^36 - 7140*a^34 + 11550*a^32 - 3010*a^30 - 1470*a^28 - 3780*a^26 + 15470*a^24 - 25620*a^22 + 18270*a^20 + 4200*a^18 - 15540*a^16 + 11340*a^14 - 16170*a^12 + 23940*a^10 - 10080*a^8 - 10080*a^6 + 12600*a^4 - 5040*a^2 + 720,\n",
       " a^112 - 16*a^98 - 56*a^88 + 168*a^86 + 8*a^84 - 112*a^82 - 70*a^80 + 1008*a^78 - 1344*a^76 + 1680*a^74 + 1148*a^72 - 8400*a^70 + 1680*a^68 - 20160*a^66 + 64372*a^64 + 11032*a^62 - 151704*a^60 + 141120*a^58 - 35280*a^56 - 3696*a^54 + 5768*a^52 - 40992*a^50 + 12299*a^48 + 104776*a^46 - 126784*a^44 + 53088*a^42 - 28924*a^40 + 14280*a^38 + 98560*a^36 - 211680*a^34 + 105840*a^32 + 99680*a^30 - 80920*a^28 - 103600*a^26 + 239400*a^24 - 290640*a^22 + 178080*a^20 + 47040*a^18 - 127050*a^16 + 139440*a^14 - 270480*a^12 + 282240*a^10 - 45360*a^8 - 141120*a^6 + 120960*a^4 - 40320*a^2 + 5040,\n",
       " a^144 - 18*a^128 - 72*a^116 + 216*a^114 + 9*a^112 - 168*a^108 + 1260*a^104 - 2016*a^102 + 2880*a^98 + 4158*a^96 - 18144*a^94 + 22680*a^92 - 28560*a^90 - 45828*a^88 + 92160*a^86 - 20088*a^84 + 452592*a^82 - 794430*a^80 - 474768*a^78 + 1892520*a^76 - 1435392*a^74 + 373296*a^72 - 75600*a^70 + 13104*a^68 - 173376*a^66 + 568134*a^64 + 102816*a^62 - 1342656*a^60 + 1214640*a^58 - 275184*a^56 - 252672*a^54 + 491400*a^52 - 38808*a^50 - 1328922*a^48 + 1918224*a^46 - 1010520*a^44 + 620928*a^42 - 1473696*a^40 + 1658160*a^38 + 108920*a^36 - 1895040*a^34 + 525420*a^32 + 1975680*a^30 - 695520*a^28 - 3296160*a^26 + 5720400*a^24 - 4959360*a^22 + 1648080*a^20 + 1043280*a^18 - 1459080*a^16 + 2600640*a^14 - 4505760*a^12 + 3265920*a^10 + 272160*a^8 - 1935360*a^6 + 1270080*a^4 - 362880*a^2 + 40320,\n",
       " a^180 - 20*a^162 - 90*a^148 + 270*a^146 + 10*a^144 - 240*a^138 + 2160*a^134 - 3300*a^132 - 252*a^130 - 180*a^128 + 5040*a^126 + 3780*a^124 - 22680*a^122 + 25200*a^120 + 15120*a^118 - 51705*a^116 - 113850*a^114 + 302490*a^112 - 483840*a^110 + 359520*a^108 + 982800*a^106 - 743400*a^104 + 1037520*a^102 - 8957160*a^100 + 9327600*a^98 + 11079300*a^96 - 25583040*a^94 + 16553160*a^92 - 3901080*a^90 - 473040*a^88 + 915840*a^86 - 195420*a^84 + 4577580*a^82 - 8085294*a^80 - 4562880*a^78 + 18819360*a^76 - 14815080*a^74 + 4579680*a^72 - 1935360*a^70 + 5062680*a^68 - 7081620*a^66 - 2570400*a^64 + 21299040*a^62 - 28350420*a^60 + 20107080*a^58 - 8638560*a^56 - 12048960*a^54 + 33927390*a^52 - 15133860*a^50 - 35130270*a^48 + 46090800*a^46 - 15397200*a^44 + 6347040*a^42 - 18579960*a^40 + 11415600*a^38 + 10027500*a^36 - 10773000*a^34 - 21621600*a^32 + 44553600*a^30 - 9563400*a^28 - 62748000*a^26 + 107604000*a^24 - 83764800*a^22 + 21795480*a^20 + 11264400*a^18 - 22415400*a^16 + 55339200*a^14 - 74239200*a^12 + 36469440*a^10 + 14061600*a^8 - 27216000*a^6 + 14515200*a^4 - 3628800*a^2 + 362880]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of strongly connected labeled digraphs\n",
    "# we take the 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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0,\n",
       " 1,\n",
       " 1,\n",
       " 18,\n",
       " 1606,\n",
       " 565080,\n",
       " 734774776,\n",
       " 3523091615568,\n",
       " 63519209389664176,\n",
       " 4400410978376102609280,\n",
       " 1190433705317814685295399296]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The numbers of strongly connected labeled digraphs\n",
    "SCDsub = [scd[i][0][0].subs(a=sqrt(2)) * i.factorial() for i in srange(disp)]\n",
    "SCDsub"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:black; font-weight: bold; font-size:18px\">\n",
    "    Semistrong digraphs, counting (strongly) connected components\n",
    "</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of semi-strong labeled digraphs\n",
    "# with the marking variable t for strongly connected components\n",
    "# we take the parameter a = sqrt{2}\n",
    "ssdt = exp(t*scd)\n",
    "SSDT = [ssdt[i] * i.factorial() for i in srange(disps)]\n",
    "SSDT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " t,\n",
       " t^2 + t,\n",
       " t^3 + 3*t^2 + 18*t,\n",
       " t^4 + 6*t^3 + 75*t^2 + 1606*t,\n",
       " t^5 + 10*t^4 + 195*t^3 + 8210*t^2 + 565080*t,\n",
       " t^6 + 15*t^5 + 405*t^4 + 25185*t^3 + 3417810*t^2 + 734774776*t,\n",
       " t^7 + 21*t^6 + 735*t^5 + 60095*t^4 + 12059880*t^3 + 5156301892*t^2 + 3523091615568*t,\n",
       " t^8 + 28*t^7 + 1218*t^6 + 122920*t^5 + 32424945*t^4 + 20677149388*t^3 + 28205966492172*t^2 + 63519209389664176*t,\n",
       " t^9 + 36*t^8 + 1890*t^7 + 226296*t^6 + 73564785*t^5 + 62188580244*t^4 + 127022651926380*t^3 + 571800941131927824*t^2 + 4400410978376102609280*t,\n",
       " t^10 + 45*t^9 + 2790*t^8 + 385770*t^7 + 148370985*t^6 + 155867605005*t^5 + 423729028084320*t^4 + 2859645482479547100*t^3 + 44006976046105332767760*t^2 + 1190433705317814685295399296*t]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of semi-strong labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "SSDTsub = [sum(ssdt[i][0][j].subs(a=sqrt(2)) * i.factorial() * t^j for j in srange(N)) for i in srange(disp)]\n",
    "SSDTsub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[t,\n",
       " -2*a^2*t*w,\n",
       " ((-2*a^7 + 6*a^5)*t)*w^2 + (2*a^3*t^2 + 2*a^3*t)*w,\n",
       " ((-2*a^15 + 18*a^11 - 24*a^9)*t)*w^3 + (-12*a^5*t^2 - 12*a^5*t)*w^2,\n",
       " ((-2*a^26 + 24*a^20 + 18*a^18 - 144*a^16 + 120*a^14)*t)*w^4 + ((-24*a^10 + 72*a^8)*t^2 + (-24*a^10 + 72*a^8)*t)*w^3 + (12*a^6*t^3 + (12*a^10 - 24*a^8 + 36*a^6)*t^2 + (12*a^10 - 24*a^8 + 24*a^6)*t)*w^2]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CGF of semi-strong labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "# we take the parameter a = sqrt{2}\n",
    "qssdt = t*exp((t+1)*scd.subs(z=a^3*z^2*w))*exp_had_prod(PhiTwoOneTwo,(1-it.subs(z=a^2*z*w))*(1-it.subs(z=a^2*z*w)),N)\n",
    "QSSDT = [qssdt[i] * i.factorial() for i in srange(disps)]\n",
    "QSSDT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[t,\n",
       " -4*t*w,\n",
       " 4*sqrt(2)*t^2*w + 8*sqrt(2)*t*w^2 + 4*sqrt(2)*t*w,\n",
       " -48*sqrt(2)*t^2*w^2 - 64*sqrt(2)*t*w^3 - 48*sqrt(2)*t*w^2,\n",
       " 96*t^3*w^2 + 384*t^2*w^3 - 4096*t*w^4 + 288*t^2*w^2 + 384*t*w^3 + 192*t*w^2]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CGF of semi-strong labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "QSSDTsub = [sum(sum(QSSDT[i][k][j].subs(a=sqrt(2)) * t^j * w^k for j in srange(N)) for k in srange(N)) for i in srange(disps)]\n",
    "QSSDTsub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[t,\n",
       " -2*a^2*t*w,\n",
       " ((-a^8 + 3*a^6)*t)*w^2 + (a^4*t^2 + a^4*t)*w,\n",
       " ((-1/3*a^18 + 3*a^14 - 4*a^12)*t)*w^3 + (-2*a^8*t^2 - 2*a^8*t)*w^2,\n",
       " ((-1/12*a^32 + a^26 + 3/4*a^24 - 6*a^22 + 5*a^20)*t)*w^4 + ((-a^16 + 3*a^14)*t^2 + (-a^16 + 3*a^14)*t)*w^3 + (1/2*a^12*t^3 + (1/2*a^16 - a^14 + 3/2*a^12)*t^2 + (1/2*a^16 - a^14 + a^12)*t)*w^2,\n",
       " ((-1/60*a^50 + 1/4*a^42 + 1/2*a^38 - 2*a^36 - 3*a^34 + 10*a^32 - 6*a^30)*t)*w^5 + ((-1/3*a^28 + 3*a^24 - 4*a^22)*t^2 + (-1/3*a^28 + 3*a^24 - 4*a^22)*t)*w^4 + (-a^18*t^3 + (-a^22 + 2*a^20 - 3*a^18)*t^2 + (-a^22 + 2*a^20 - 2*a^18)*t)*w^3,\n",
       " ((-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)*t)*w^6 + ((-1/12*a^44 + a^38 + 3/4*a^36 - 6*a^34 + 5*a^32)*t^2 + (-1/12*a^44 + a^38 + 3/4*a^36 - 6*a^34 + 5*a^32)*t)*w^5 + ((-1/2*a^28 + 3/2*a^26)*t^3 + (-1/2*a^32 + 5/2*a^30 - 9/2*a^28 + 9/2*a^26)*t^2 + (-1/2*a^32 + 5/2*a^30 - 4*a^28 + 3*a^26)*t)*w^4 + (1/6*a^24*t^4 + (1/2*a^28 - a^26 + a^24)*t^3 + (1/6*a^36 - a^32 + a^30 + 3/2*a^28 - 3*a^26 + 11/6*a^24)*t^2 + (1/6*a^36 - a^32 + a^30 + a^28 - 2*a^26 + a^24)*t)*w^3,\n",
       " ((-1/2520*a^98 + 1/120*a^86 + 1/40*a^78 - 1/10*a^76 + 1/24*a^74 - 1/2*a^70 + 1/2*a^68 - 1/2*a^66 + 5*a^64 - 5/2*a^62 - 15*a^60 + 21*a^58 - 8*a^56)*t)*w^7 + ((-1/60*a^64 + 1/4*a^56 + 1/2*a^52 - 2*a^50 - 3*a^48 + 10*a^46 - 6*a^44)*t^2 + (-1/60*a^64 + 1/4*a^56 + 1/2*a^52 - 2*a^50 - 3*a^48 + 10*a^46 - 6*a^44)*t)*w^6 + ((-1/6*a^42 + 3/2*a^38 - 2*a^36)*t^3 + (-1/6*a^46 + 1/3*a^44 + a^42 - 5*a^40 + 17/2*a^38 - 6*a^36)*t^2 + (-1/6*a^46 + 1/3*a^44 + 7/6*a^42 - 5*a^40 + 7*a^38 - 4*a^36)*t)*w^5 + (-1/3*a^32*t^4 + (-a^36 + 2*a^34 - 2*a^32)*t^3 + (-1/3*a^44 + 2*a^40 - 2*a^38 - 3*a^36 + 6*a^34 - 11/3*a^32)*t^2 + (-1/3*a^44 + 2*a^40 - 2*a^38 - 2*a^36 + 4*a^34 - 2*a^32)*t)*w^4,\n",
       " ((-1/20160*a^128 + 1/840*a^114 + 1/240*a^104 - 1/60*a^102 + 1/120*a^98 + 1/192*a^96 - 1/10*a^94 + 1/6*a^92 - 1/6*a^90 - 1/8*a^88 + 13/12*a^86 - 5/12*a^84 + 5/2*a^82 - 155/16*a^80 - 1/2*a^78 + 105/4*a^76 - 28*a^74 + 9*a^72)*t)*w^8 + ((-1/360*a^88 + 1/20*a^78 + 1/8*a^72 - 5/12*a^70 - 2*a^66 + 17/6*a^64 + 15/2*a^62 - 15*a^60 + 7*a^58)*t^2 + (-1/360*a^88 + 1/20*a^78 + 1/8*a^72 - 5/12*a^70 - 2*a^66 + 17/6*a^64 + 15/2*a^62 - 15*a^60 + 7*a^58)*t)*w^7 + ((-1/24*a^60 + 1/2*a^54 + 3/8*a^52 - 3*a^50 + 5/2*a^48)*t^3 + (-1/24*a^64 + 1/12*a^62 - 1/8*a^60 + 1/2*a^58 - 5/8*a^56 - 9/4*a^54 + 77/8*a^52 - 14*a^50 + 15/2*a^48)*t^2 + (-1/24*a^64 + 1/12*a^62 - 1/12*a^60 + 1/2*a^58 - 5/8*a^56 - 11/4*a^54 + 37/4*a^52 - 11*a^50 + 5*a^48)*t)*w^6 + ((-1/6*a^44 + 1/2*a^42)*t^4 + (-1/2*a^48 + 5/2*a^46 - 4*a^44 + 3*a^42)*t^3 + (-1/6*a^56 + 1/2*a^54 + a^52 - 4*a^50 + 3/2*a^48 + 15/2*a^46 - 65/6*a^44 + 11/2*a^42)*t^2 + (-1/6*a^56 + 1/2*a^54 + a^52 - 4*a^50 + 2*a^48 + 5*a^46 - 7*a^44 + 3*a^42)*t)*w^5 + (1/24*a^40*t^5 + (1/4*a^44 - 1/2*a^42 + 5/12*a^40)*t^4 + (1/6*a^52 - 7/8*a^48 + 1/2*a^46 + 2*a^44 - 3*a^42 + 35/24*a^40)*t^3 + (1/24*a^64 - 1/3*a^58 - 1/4*a^56 + 3/2*a^54 - 1/2*a^52 - 21/8*a^48 + 3/2*a^46 + 17/4*a^44 - 11/2*a^42 + 25/12*a^40)*t^2 + (1/24*a^64 - 1/3*a^58 - 1/4*a^56 + 3/2*a^54 - 2/3*a^52 - 7/4*a^48 + a^46 + 5/2*a^44 - 3*a^42 + a^40)*t)*w^4,\n",
       " ((-1/181440*a^162 + 1/6720*a^146 + 1/1680*a^134 - 1/420*a^132 + 1/720*a^126 - 7/480*a^122 + 1/36*a^120 - 1/30*a^116 - 11/240*a^114 + 1/4*a^112 - 1/3*a^110 + 43/108*a^108 + 5/8*a^106 - 5/3*a^104 + 1/2*a^102 - 15/2*a^100 + 125/8*a^98 + 49/6*a^96 - 42*a^94 + 36*a^92 - 10*a^90)*t)*w^9 + ((-1/2520*a^116 + 1/120*a^104 + 1/40*a^96 - 1/10*a^94 + 1/24*a^92 - 1/2*a^88 + 1/2*a^86 - 1/2*a^84 + 5*a^82 - 5/2*a^80 - 15*a^78 + 21*a^76 - 8*a^74)*t^2 + (-1/2520*a^116 + 1/120*a^104 + 1/40*a^96 - 1/10*a^94 + 1/24*a^92 - 1/2*a^88 + 1/2*a^86 - 1/2*a^84 + 5*a^82 - 5/2*a^80 - 15*a^78 + 21*a^76 - 8*a^74)*t)*w^8 + ((-1/120*a^82 + 1/8*a^74 + 1/4*a^70 - a^68 - 3/2*a^66 + 5*a^64 - 3*a^62)*t^3 + (-1/120*a^86 + 1/60*a^84 - 1/40*a^82 + 1/8*a^78 - 1/4*a^76 + 5/8*a^74 - 3/2*a^72 + 5/4*a^70 + 5*a^68 - 35/2*a^66 + 21*a^64 - 9*a^62)*t^2 + (-1/120*a^86 + 1/60*a^84 - 1/60*a^82 + 1/8*a^78 - 1/4*a^76 + 1/2*a^74 - 3/2*a^72 + a^70 + 6*a^68 - 16*a^66 + 16*a^64 - 6*a^62)*t)*w^7 + ((-1/18*a^60 + 1/2*a^56 - 2/3*a^54)*t^4 + (-1/6*a^64 + 1/3*a^62 + 7/6*a^60 - 5*a^58 + 7*a^56 - 4*a^54)*t^3 + (-1/18*a^72 + 5/6*a^68 - a^66 - 7/2*a^64 + 8*a^62 - 1/9*a^60 - 15*a^58 + 35/2*a^56 - 22/3*a^54)*t^2 + (-1/18*a^72 + 5/6*a^68 - a^66 - 10/3*a^64 + 23/3*a^62 - 4/3*a^60 - 10*a^58 + 11*a^56 - 4*a^54)*t)*w^6 + (-1/12*a^50*t^5 + (-1/2*a^54 + a^52 - 5/6*a^50)*t^4 + (-1/3*a^62 + 7/4*a^58 - a^56 - 4*a^54 + 6*a^52 - 35/12*a^50)*t^3 + (-1/12*a^74 + 2/3*a^68 + 1/2*a^66 - 3*a^64 + a^62 + 21/4*a^58 - 3*a^56 - 17/2*a^54 + 11*a^52 - 25/6*a^50)*t^2 + (-1/12*a^74 + 2/3*a^68 + 1/2*a^66 - 3*a^64 + 4/3*a^62 + 7/2*a^58 - 2*a^56 - 5*a^54 + 6*a^52 - 2*a^50)*t)*w^5,\n",
       " ((-1/1814400*a^200 + 1/60480*a^182 + 1/13440*a^168 - 1/3360*a^166 + 1/5040*a^158 - 1/420*a^154 + 29/6720*a^152 + 1/4800*a^150 - 1/180*a^146 - 1/240*a^144 + 1/30*a^142 - 1/24*a^140 - 1/60*a^138 + 7/96*a^136 + 47/288*a^134 - 1/2*a^132 + 23/30*a^130 - 275/432*a^128 - 5/3*a^126 + 15/8*a^124 - 23/12*a^122 + 277/16*a^120 - 343/16*a^118 - 23*a^116 + 63*a^114 - 45*a^112 + 11*a^110)*t)*w^10 + ((-1/20160*a^148 + 1/840*a^134 + 1/240*a^124 - 1/60*a^122 + 1/120*a^118 + 1/192*a^116 - 1/10*a^114 + 1/6*a^112 - 1/6*a^110 - 1/8*a^108 + 13/12*a^106 - 5/12*a^104 + 5/2*a^102 - 155/16*a^100 - 1/2*a^98 + 105/4*a^96 - 28*a^94 + 9*a^92)*t^2 + (-1/20160*a^148 + 1/840*a^134 + 1/240*a^124 - 1/60*a^122 + 1/120*a^118 + 1/192*a^116 - 1/10*a^114 + 1/6*a^112 - 1/6*a^110 - 1/8*a^108 + 13/12*a^106 - 5/12*a^104 + 5/2*a^102 - 155/16*a^100 - 1/2*a^98 + 105/4*a^96 - 28*a^94 + 9*a^92)*t)*w^9 + ((-1/720*a^108 + 1/40*a^98 + 1/16*a^92 - 5/24*a^90 - a^86 + 17/12*a^84 + 15/4*a^82 - 15/2*a^80 + 7/2*a^78)*t^3 + (-1/720*a^112 + 1/360*a^110 - 1/240*a^108 + 1/40*a^102 - 1/20*a^100 + 3/40*a^98 + 1/16*a^96 - 1/3*a^94 + 29/48*a^92 - 13/8*a^90 + 41/12*a^88 - 25/12*a^86 - 43/4*a^84 + 119/4*a^82 - 59/2*a^80 + 21/2*a^78)*t^2 + (-1/720*a^112 + 1/360*a^110 - 1/360*a^108 + 1/40*a^102 - 1/20*a^100 + 1/20*a^98 + 1/16*a^96 - 1/3*a^94 + 13/24*a^92 - 17/12*a^90 + 41/12*a^88 - 13/12*a^86 - 73/6*a^84 + 26*a^82 - 22*a^80 + 7*a^78)*t)*w^8 + ((-1/72*a^80 + 1/6*a^74 + 1/8*a^72 - a^70 + 5/6*a^68)*t^4 + (-1/24*a^84 + 1/12*a^82 - 1/12*a^80 + 1/2*a^78 - 5/8*a^76 - 11/4*a^74 + 37/4*a^72 - 11*a^70 + 5*a^68)*t^3 + (-1/72*a^92 + 1/12*a^88 + 1/12*a^86 - 7/4*a^82 + 67/72*a^80 + 33/4*a^78 - 103/8*a^76 - 53/12*a^74 + 215/8*a^72 - 26*a^70 + 55/6*a^68)*t^2 + (-1/72*a^92 + 1/12*a^88 + 1/12*a^86 + 1/24*a^84 - 11/6*a^82 + a^80 + 31/4*a^78 - 49/4*a^76 - 3/2*a^74 + 71/4*a^72 - 16*a^70 + 5*a^68)*t)*w^7 + ((-1/24*a^64 + 1/8*a^62)*t^5 + (-1/4*a^68 + 5/4*a^66 - 23/12*a^64 + 5/4*a^62)*t^4 + (-1/6*a^76 + 1/2*a^74 + 7/8*a^72 - 25/8*a^70 - 1/2*a^68 + 9*a^66 - 251/24*a^64 + 35/8*a^62)*t^3 + (-1/24*a^88 + 1/8*a^86 + 1/3*a^82 - 3/4*a^80 - 9/4*a^78 + 5*a^76 - 3/2*a^74 + 21/8*a^72 - 75/8*a^70 + 1/4*a^68 + 73/4*a^66 - 223/12*a^64 + 25/4*a^62)*t^2 + (-1/24*a^88 + 1/8*a^86 + 1/3*a^82 - 3/4*a^80 - 9/4*a^78 + 31/6*a^76 - 2*a^74 + 7/4*a^72 - 25/4*a^70 + 1/2*a^68 + 21/2*a^66 - 10*a^64 + 3*a^62)*t)*w^6 + (1/120*a^60*t^6 + (1/12*a^64 - 1/6*a^62 + 1/8*a^60)*t^5 + (1/12*a^72 - 3/8*a^68 + 4/3*a^64 - 5/3*a^62 + 17/24*a^60)*t^4 + (1/24*a^84 - 1/3*a^78 - 1/6*a^76 + 4/3*a^74 - a^72 + 3/2*a^70 - 13/4*a^68 + 71/12*a^64 - 35/6*a^62 + 15/8*a^60)*t^3 + (1/120*a^100 - 1/12*a^92 - 1/6*a^88 + 1/2*a^86 + 7/8*a^84 - 2*a^82 + a^80 - a^78 - 1/2*a^76 + 4*a^74 - 43/12*a^72 + 9/2*a^70 - 57/8*a^68 + 29/3*a^64 - 25/3*a^62 + 137/60*a^60)*t^2 + (1/120*a^100 - 1/12*a^92 - 1/6*a^88 + 1/2*a^86 + 5/6*a^84 - 2*a^82 + a^80 - 2/3*a^78 - 1/3*a^76 + 8/3*a^74 - 5/2*a^72 + 3*a^70 - 17/4*a^68 + 5*a^64 - 4*a^62 + a^60)*t)*w^5]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of semi-strong labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "# we take the parameter a = sqrt{2}\n",
    "QSSDTzero = [qssdt[i] * a^(i*(i-1)/2) for i in srange(disp)]\n",
    "QSSDTzero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[t,\n",
       " -4*t*w,\n",
       " 4*t^2*w + 8*t*w^2 + 4*t*w,\n",
       " -32*t^2*w^2 - 128/3*t*w^3 - 32*t*w^2,\n",
       " 32*t^3*w^2 + 128*t^2*w^3 - 4096/3*t*w^4 + 96*t^2*w^2 + 128*t*w^3 + 64*t*w^2,\n",
       " -512*t^3*w^3 - 4096/3*t^2*w^4 - 3473408/15*t*w^5 - 1536*t^2*w^3 - 4096/3*t*w^4 - 1024*t*w^3,\n",
       " 2048/3*t^4*w^3 + 4096*t^3*w^4 - 262144/3*t^2*w^5 - 4984930304/45*t*w^6 + 4096*t^3*w^3 + 12288*t^2*w^4 - 262144/3*t*w^5 + 18432*t^2*w^3 + 8192*t*w^4 + 45056/3*t*w^3,\n",
       " -65536/3*t^4*w^4 - 262144/3*t^3*w^5 - 444596224/15*t^2*w^6 - 50988241125376/315*t*w^7 - 131072*t^3*w^4 - 262144*t^2*w^5 - 444596224/15*t*w^6 - 589824*t^2*w^4 - 524288/3*t*w^5 - 1441792/3*t*w^4,\n",
       " 131072/3*t^5*w^4 + 1048576/3*t^4*w^5 - 33554432/3*t^3*w^6 - 1276142157824/45*t^2*w^7 - 239467516496183296/315*t*w^8 + 1310720/3*t^4*w^4 + 2097152*t^3*w^5 - 33554432*t^2*w^6 - 1276142157824/45*t*w^7 + 4325376*t^3*w^4 + 9437184*t^2*w^5 - 67108864/3*t*w^6 + 233046016/3*t^2*w^4 + 23068672/3*t*w^5 + 221249536/3*t*w^4,\n",
       " -8388608/3*t^5*w^5 - 134217728/9*t^4*w^6 - 113816633344/15*t^3*w^7 - 26105979456192512/315*t^2*w^8 - 34025482048081491918848/2835*t*w^9 - 83886080/3*t^4*w^5 - 268435456/3*t^3*w^6 - 113816633344/5*t^2*w^7 - 26105979456192512/315*t*w^8 - 276824064*t^3*w^5 - 402653184*t^2*w^6 - 227633266688/15*t*w^7 - 14914945024/3*t^2*w^5 - 2952790016/9*t*w^6 - 14159970304/3*t*w^5,\n",
       " 134217728/15*t^6*w^5 + 268435456/3*t^5*w^6 - 34359738368/9*t^4*w^7 - 653384784805888/45*t^3*w^8 - 245214736892091695104/315*t^2*w^9 - 9342744711155730849185398784/14175*t*w^10 + 134217728*t^5*w^5 + 2684354560/3*t^4*w^6 - 68719476736/3*t^3*w^7 - 653384784805888/15*t^2*w^8 - 245214736892091695104/315*t*w^9 + 6576668672/3*t^4*w^5 + 8858370048*t^3*w^6 - 103079215104*t^2*w^7 - 1306769569611776/45*t*w^8 + 79322677248*t^3*w^5 + 477278240768/3*t^2*w^6 - 755914244096/9*t*w^7 + 15626433200128/3*t^2*w^5 + 453119049728/3*t*w^6 + 76973330137088/15*t*w^5]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of semi-strong labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "QSSDTzerosub = [sum(sum(QSSDTzero[i][k][j].subs(a=sqrt(2)) * t^j * w^k for j in srange(N)) for k in srange(N)) for i in srange(disp)]\n",
    "QSSDTzerosub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[t, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, -4*t, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 4*t^2 + 4*t, 8*t, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, -32*t^2 - 32*t, -128/3*t, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 32*t^3 + 96*t^2 + 64*t, 128*t^2 + 128*t, -4096/3*t, 0, 0, 0, 0, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  -512*t^3 - 1536*t^2 - 1024*t,\n",
       "  -4096/3*t^2 - 4096/3*t,\n",
       "  -3473408/15*t,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  2048/3*t^4 + 4096*t^3 + 18432*t^2 + 45056/3*t,\n",
       "  4096*t^3 + 12288*t^2 + 8192*t,\n",
       "  -262144/3*t^2 - 262144/3*t,\n",
       "  -4984930304/45*t,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  -65536/3*t^4 - 131072*t^3 - 589824*t^2 - 1441792/3*t,\n",
       "  -262144/3*t^3 - 262144*t^2 - 524288/3*t,\n",
       "  -444596224/15*t^2 - 444596224/15*t,\n",
       "  -50988241125376/315*t,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  131072/3*t^5 + 1310720/3*t^4 + 4325376*t^3 + 233046016/3*t^2 + 221249536/3*t,\n",
       "  1048576/3*t^4 + 2097152*t^3 + 9437184*t^2 + 23068672/3*t,\n",
       "  -33554432/3*t^3 - 33554432*t^2 - 67108864/3*t,\n",
       "  -1276142157824/45*t^2 - 1276142157824/45*t,\n",
       "  -239467516496183296/315*t,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  -8388608/3*t^5 - 83886080/3*t^4 - 276824064*t^3 - 14914945024/3*t^2 - 14159970304/3*t,\n",
       "  -134217728/9*t^4 - 268435456/3*t^3 - 402653184*t^2 - 2952790016/9*t,\n",
       "  -113816633344/15*t^3 - 113816633344/5*t^2 - 227633266688/15*t,\n",
       "  -26105979456192512/315*t^2 - 26105979456192512/315*t,\n",
       "  -34025482048081491918848/2835*t,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  134217728/15*t^6 + 134217728*t^5 + 6576668672/3*t^4 + 79322677248*t^3 + 15626433200128/3*t^2 + 76973330137088/15*t,\n",
       "  268435456/3*t^5 + 2684354560/3*t^4 + 8858370048*t^3 + 477278240768/3*t^2 + 453119049728/3*t,\n",
       "  -34359738368/9*t^4 - 68719476736/3*t^3 - 103079215104*t^2 - 755914244096/9*t,\n",
       "  -653384784805888/45*t^3 - 653384784805888/15*t^2 - 1306769569611776/45*t,\n",
       "  -245214736892091695104/315*t^2 - 245214736892091695104/315*t,\n",
       "  -9342744711155730849185398784/14175*t]]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of semi-strong labeled digraphs in matrix form (Table 5)\n",
    "# the marking variable t is for strongly connected components\n",
    "QSSDTzeroMatrix = [[QSSDTzero[i][j].subs(a=sqrt(2)) for j in srange(disp)] for i in srange(disp)]\n",
    "QSSDTzeroMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, -4, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, -64, -128/3, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 192, 256, -4096/3, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, -3072, -8192/3, -3473408/15, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 114688/3, 24576, -524288/3, -4984930304/45, 0, 0, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  -3670016/3,\n",
       "  -524288,\n",
       "  -889192448/15,\n",
       "  -50988241125376/315,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  156237824,\n",
       "  58720256/3,\n",
       "  -67108864,\n",
       "  -2552284315648/45,\n",
       "  -239467516496183296/315,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  -9999220736,\n",
       "  -7516192768/9,\n",
       "  -227633266688/5,\n",
       "  -52211958912385024/315,\n",
       "  -34025482048081491918848/2835,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  156330367123456/15,\n",
       "  319975063552,\n",
       "  -1924145348608/9,\n",
       "  -1306769569611776/15,\n",
       "  -490429473784183390208/315,\n",
       "  -9342744711155730849185398784/14175]]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of semi-strong labeled digraphs in matrix form (Table 6)\n",
    "QSSDzeroMatrix = [[QSSDTzeroMatrix[i][j].subs(t=1) for j in srange(disp)] for i in srange(disp)]\n",
    "QSSDzeroMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, -4, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 4, 8, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, -32, -128/3, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 64, 128, -4096/3, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, -1024, -4096/3, -3473408/15, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 45056/3, 8192, -262144/3, -4984930304/45, 0, 0, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  -1441792/3,\n",
       "  -524288/3,\n",
       "  -444596224/15,\n",
       "  -50988241125376/315,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  221249536/3,\n",
       "  23068672/3,\n",
       "  -67108864/3,\n",
       "  -1276142157824/45,\n",
       "  -239467516496183296/315,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  -14159970304/3,\n",
       "  -2952790016/9,\n",
       "  -227633266688/15,\n",
       "  -26105979456192512/315,\n",
       "  -34025482048081491918848/2835,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  76973330137088/15,\n",
       "  453119049728/3,\n",
       "  -755914244096/9,\n",
       "  -1306769569611776/45,\n",
       "  -245214736892091695104/315,\n",
       "  -9342744711155730849185398784/14175]]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of semi-strong labeled digraphs with one strongly connected component (Table 4)\n",
    "QSSDTzeroMatrixPrim = [[QSSDTzeroMatrix[i][j].diff(t).subs(t=0) for j in srange(disp)] for i in srange(disp)]\n",
    "QSSDTzeroMatrixPrim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, -32, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 96, 128, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, -1536, -4096/3, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 18432, 12288, -262144/3, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, -589824, -262144, -444596224/15, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 233046016/3, 9437184, -33554432, -1276142157824/45, 0, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  -14914945024/3,\n",
       "  -402653184,\n",
       "  -113816633344/5,\n",
       "  -26105979456192512/315,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  15626433200128/3,\n",
       "  477278240768/3,\n",
       "  -103079215104,\n",
       "  -653384784805888/15,\n",
       "  -245214736892091695104/315,\n",
       "  0]]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of semi-strong labeled digraphs with two strongly connected components (Table 7)\n",
    "QSSDTzeroMatrixSec = [[QSSDTzeroMatrix[i][j].diff(t,2).subs(t=0) / 2 for j in srange(disp)] for i in srange(disp)]\n",
    "QSSDTzeroMatrixSec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 32, 0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, -512, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 4096, 4096, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, -131072, -262144/3, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 4325376, 2097152, -33554432/3, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, -276824064, -268435456/3, -113816633344/15, 0, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  79322677248,\n",
       "  8858370048,\n",
       "  -68719476736/3,\n",
       "  -653384784805888/45,\n",
       "  0,\n",
       "  0]]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of semi-strong labeled digraphs with three strongly connected components (Table 8)\n",
    "QSSDTzeroMatrixThir = [[QSSDTzeroMatrix[i][j].diff(t,3).subs(t=0) / 6 for j in srange(disp)] for i in srange(disp)]\n",
    "QSSDTzeroMatrixThir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:black; font-weight: bold; font-size:18px\">\n",
    "    Labeled digraphs, counting (strongly) connected components\n",
    "</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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,\n",
       " ((120*a^20 - 240*a^18 + 90*a^16 + 60*a^14 - 20*a^12 - 10*a^8 + 1)/a^20)*t^5 + ((240*a^22 - 660*a^20 + 420*a^18 + 260*a^16 - 340*a^14 + 140*a^12 - 120*a^10 + 60*a^8 - 10*a^4 + 20*a^2 - 10)/a^20)*t^4 + ((60*a^26 + 70*a^24 - 720*a^22 + 920*a^20 - 60*a^18 - 480*a^16 + 360*a^14 - 390*a^12 + 360*a^10 - 155*a^8 + 120*a^4 - 120*a^2 + 35)/a^20)*t^3 + ((10*a^32 + 20*a^28 - 120*a^26 - 165*a^24 + 720*a^22 - 500*a^20 - 80*a^18 + 150*a^16 - 240*a^14 + 430*a^12 - 420*a^10 + 315*a^8 - 290*a^4 + 220*a^2 - 50)/a^20)*t^2 + ((a^40 - 10*a^32 - 20*a^28 + 60*a^26 + 95*a^24 - 240*a^22 + 120*a^20 - 40*a^18 - 20*a^16 + 160*a^14 - 160*a^12 + 180*a^10 - 210*a^8 + 180*a^4 - 120*a^2 + 24)/a^20)*t,\n",
       " ((720*a^30 - 1800*a^28 + 1080*a^26 + 390*a^24 - 360*a^22 - 70*a^18 + 30*a^16 + 12*a^10 - 1)/a^30)*t^6 + ((1800*a^32 - 5760*a^30 + 4950*a^28 + 1620*a^26 - 4050*a^24 + 1860*a^22 - 1050*a^20 + 840*a^18 - 210*a^16 - 120*a^14 + 240*a^12 - 120*a^10 + 15*a^4 - 30*a^2 + 15)/a^30)*t^5 + ((480*a^36 + 720*a^34 - 7470*a^32 + 11200*a^30 - 2220*a^28 - 6450*a^26 + 5160*a^24 - 3540*a^22 + 4260*a^20 - 3070*a^18 + 510*a^16 + 1440*a^14 - 1460*a^12 + 420*a^10 + 75*a^8 + 60*a^6 - 330*a^4 + 300*a^2 - 85)/a^30)*t^4 + ((90*a^42 - 30*a^40 + 360*a^38 - 1350*a^36 - 2820*a^34 + 11370*a^32 - 9240*a^30 - 1380*a^28 + 4110*a^26 - 2535*a^24 + 4800*a^22 - 7230*a^20 + 5700*a^18 - 540*a^16 - 3060*a^14 + 2055*a^12 + 390*a^10 - 990*a^8 - 480*a^6 + 1605*a^4 - 1050*a^2 + 225)/a^30)*t^3 + ((12*a^50 + 30*a^44 - 160*a^42 + 24*a^40 - 720*a^38 + 1260*a^36 + 3180*a^34 - 7440*a^32 + 3800*a^30 + 555*a^28 - 690*a^26 + 455*a^24 - 1200*a^22 + 3270*a^20 - 4000*a^18 + 1350*a^16 + 900*a^14 + 95*a^12 - 2682*a^10 + 2445*a^8 + 1020*a^6 - 2730*a^4 + 1500*a^2 - 274)/a^30)*t^2 + ((a^60 - 12*a^50 - 30*a^44 + 70*a^42 + 6*a^40 + 360*a^38 - 390*a^36 - 1080*a^34 + 1740*a^32 - 720*a^30 - 105*a^28 + 330*a^26 + 580*a^24 - 1560*a^22 + 750*a^20 + 600*a^18 - 1140*a^16 + 840*a^14 - 930*a^12 + 1980*a^10 - 1530*a^8 - 600*a^6 + 1440*a^4 - 720*a^2 + 120)/a^30)*t,\n",
       " ((5040*a^42 - 15120*a^40 + 12600*a^38 + 1680*a^36 - 5040*a^34 + 630*a^32 - 420*a^30 + 630*a^28 - 70*a^24 + 126*a^22 - 42*a^20 - 14*a^12 + 1)/a^42)*t^7 + ((15120*a^44 - 55440*a^42 + 60480*a^40 + 5040*a^38 - 48510*a^36 + 28980*a^34 - 12600*a^32 + 11340*a^30 - 3780*a^28 - 2520*a^26 + 3612*a^24 - 2184*a^22 + 462*a^20 + 210*a^16 - 420*a^14 + 210*a^12 - 21*a^4 + 42*a^2 - 21)/a^42)*t^6 + ((4200*a^48 + 7560*a^46 - 82530*a^44 + 143640*a^42 - 49770*a^40 - 85680*a^38 + 89950*a^36 - 52920*a^34 + 57540*a^32 - 43470*a^30 + 2310*a^28 + 24780*a^26 - 25130*a^24 + 11130*a^22 - 840*a^20 + 840*a^18 - 4620*a^16 + 4200*a^14 - 1155*a^12 - 105*a^8 - 210*a^6 + 735*a^4 - 630*a^2 + 175)/a^42)*t^5 + ((840*a^54 - 630*a^52 + 5040*a^50 - 15470*a^48 - 40950*a^46 + 169680*a^44 - 159740*a^42 - 18060*a^40 + 102900*a^38 - 77210*a^36 + 100380*a^34 - 127050*a^32 + 74130*a^30 + 13860*a^28 - 63840*a^26 + 51835*a^24 - 9660*a^22 - 10290*a^20 - 6440*a^18 + 22470*a^16 - 15540*a^14 + 4795*a^12 - 3150*a^10 + 2310*a^8 + 2940*a^6 - 5985*a^4 + 3570*a^2 - 735)/a^42)*t^4 + ((126*a^62 - 42*a^60 + 630*a^56 - 2100*a^54 + 1386*a^52 - 14700*a^50 + 19040*a^48 + 70770*a^46 - 164850*a^44 + 91980*a^42 + 26061*a^40 - 44100*a^38 + 58800*a^36 - 108780*a^34 + 116970*a^32 - 64260*a^30 - 10815*a^28 + 56910*a^26 - 34580*a^24 - 25956*a^22 + 33957*a^20 + 18060*a^18 - 45990*a^16 + 29190*a^14 - 17101*a^12 + 21420*a^10 - 11235*a^8 - 12390*a^6 + 19425*a^4 - 9450*a^2 + 1624)/a^42)*t^3 + ((14*a^72 + 42*a^64 - 252*a^62 + 105*a^60 - 1260*a^56 + 1680*a^54 - 1386*a^52 + 14784*a^50 - 9450*a^48 - 49980*a^46 + 77889*a^44 - 30968*a^42 - 3633*a^40 + 6720*a^38 - 21805*a^36 + 44520*a^34 - 47040*a^32 + 25690*a^30 - 735*a^28 - 11550*a^26 - 11137*a^24 + 52164*a^22 - 41517*a^20 - 16660*a^18 + 43470*a^16 - 28770*a^14 + 29435*a^12 - 42210*a^10 + 19110*a^8 + 19740*a^6 - 26754*a^4 + 11508*a^2 - 1764)/a^42)*t^2 + ((a^84 - 14*a^72 - 42*a^64 + 126*a^62 - 63*a^60 + 630*a^56 - 420*a^54 + 630*a^52 - 5124*a^50 + 1680*a^48 + 12600*a^46 - 15309*a^44 + 5488*a^42 + 42*a^40 + 2520*a^38 - 2905*a^36 - 7140*a^34 + 11550*a^32 - 3010*a^30 - 1470*a^28 - 3780*a^26 + 15470*a^24 - 25620*a^22 + 18270*a^20 + 4200*a^18 - 15540*a^16 + 11340*a^14 - 16170*a^12 + 23940*a^10 - 10080*a^8 - 10080*a^6 + 12600*a^4 - 5040*a^2 + 720)/a^42)*t,\n",
       " ((40320*a^56 - 141120*a^54 + 151200*a^52 - 10080*a^50 - 64680*a^48 + 20160*a^46 - 1680*a^44 + 8400*a^42 - 1260*a^40 - 1680*a^38 + 1344*a^36 - 1008*a^34 + 70*a^32 + 112*a^30 - 168*a^26 + 56*a^24 + 16*a^14 - 1)/a^56)*t^8 + ((141120*a^58 - 584640*a^56 + 776160*a^54 - 104160*a^52 - 588000*a^50 + 470400*a^48 - 179760*a^46 + 137760*a^44 - 63840*a^42 - 33600*a^40 + 53088*a^38 - 36456*a^36 + 11312*a^34 + 2072*a^32 + 1064*a^30 - 5936*a^28 + 4312*a^26 - 896*a^24 - 336*a^18 + 672*a^16 - 336*a^14 + 28*a^4 - 56*a^2 + 28)/a^56)*t^7 + ((40320*a^62 + 84000*a^60 - 977760*a^58 + 1950480*a^56 - 991200*a^54 - 1098720*a^52 + 1559040*a^50 - 871920*a^48 + 720720*a^46 - 594160*a^44 + 49952*a^42 + 362460*a^40 - 408912*a^38 + 215992*a^36 - 13944*a^34 - 16772*a^32 - 54656*a^30 + 73920*a^28 - 33880*a^26 + 5600*a^24 - 1680*a^22 - 3360*a^20 + 11760*a^18 - 10080*a^16 + 2800*a^14 - 56*a^12 + 126*a^8 + 504*a^6 - 1428*a^4 + 1176*a^2 - 322)/a^56)*t^6 + ((8400*a^68 - 10080*a^66 + 68460*a^64 - 189840*a^62 - 584640*a^60 + 2578800*a^58 - 2816660*a^56 - 54320*a^54 + 2101680*a^52 - 1723400*a^50 + 1595860*a^48 - 1904000*a^46 + 1138760*a^44 + 239120*a^42 - 1127140*a^40 + 1040760*a^38 - 248920*a^36 - 273280*a^34 + 47880*a^32 + 334040*a^30 - 336560*a^28 + 149240*a^26 - 67690*a^24 + 36960*a^22 + 47040*a^20 - 96320*a^18 + 57260*a^16 - 10360*a^14 - 3780*a^12 + 7560*a^10 - 3570*a^8 - 10920*a^6 + 17500*a^4 - 9800*a^2 + 1960)/a^56)*t^5 + ((1344*a^76 - 1008*a^74 + 10192*a^70 - 29400*a^68 + 37632*a^66 - 258384*a^64 + 274400*a^62 + 1388520*a^60 - 3363920*a^58 + 2121700*a^56 + 720216*a^54 - 1504440*a^52 + 1512840*a^50 - 2233490*a^48 + 2265760*a^46 - 1128400*a^44 - 435120*a^42 + 1486744*a^40 - 957320*a^38 - 566104*a^36 + 1003128*a^34 + 25830*a^32 - 878192*a^30 + 756700*a^28 - 428792*a^26 + 364644*a^24 - 162960*a^22 - 206640*a^20 + 296240*a^18 - 120785*a^16 - 12376*a^14 + 69160*a^12 - 92400*a^10 + 25830*a^8 + 76440*a^6 - 91980*a^4 + 41160*a^2 - 6769)/a^56)*t^4 + ((168*a^86 - 56*a^84 + 1008*a^78 - 4032*a^76 + 2912*a^74 + 1932*a^72 - 28896*a^70 + 33544*a^68 - 68320*a^66 + 383712*a^64 - 140000*a^62 - 1499932*a^60 + 2358944*a^58 - 958328*a^56 - 349048*a^54 + 547120*a^52 - 909776*a^50 + 1301440*a^48 - 1148560*a^46 + 612248*a^44 + 253064*a^42 - 781844*a^40 - 71568*a^38 + 1475936*a^36 - 1402912*a^34 - 105182*a^32 + 1071336*a^30 - 846664*a^28 + 585088*a^26 - 570234*a^24 + 110040*a^22 + 430080*a^20 - 356944*a^18 + 28378*a^16 + 146776*a^14 - 315420*a^12 + 365400*a^10 - 77070*a^8 - 231000*a^6 + 234472*a^4 - 90944*a^2 + 13132)/a^56)*t^3 + ((16*a^98 + 56*a^88 - 336*a^86 + 48*a^84 + 112*a^82 + 70*a^80 - 2016*a^78 + 4032*a^76 - 3584*a^74 - 3080*a^72 + 27104*a^70 - 14224*a^68 + 60928*a^66 - 258160*a^64 + 4088*a^62 + 763756*a^60 - 878304*a^58 + 282408*a^56 + 43008*a^54 - 98448*a^52 + 200368*a^50 - 209909*a^48 + 120904*a^46 - 37744*a^44 - 104664*a^42 + 123564*a^40 + 331352*a^38 - 940352*a^36 + 888384*a^34 - 59738*a^32 - 573384*a^30 + 439460*a^28 - 172200*a^26 + 29120*a^24 + 308280*a^22 - 445200*a^20 + 98560*a^18 + 171605*a^16 - 265960*a^14 + 520576*a^12 - 562800*a^10 + 100044*a^8 + 306096*a^6 - 279552*a^4 + 98784*a^2 - 13068)/a^56)*t^2 + ((a^112 - 16*a^98 - 56*a^88 + 168*a^86 + 8*a^84 - 112*a^82 - 70*a^80 + 1008*a^78 - 1344*a^76 + 1680*a^74 + 1148*a^72 - 8400*a^70 + 1680*a^68 - 20160*a^66 + 64372*a^64 + 11032*a^62 - 151704*a^60 + 141120*a^58 - 35280*a^56 - 3696*a^54 + 5768*a^52 - 40992*a^50 + 12299*a^48 + 104776*a^46 - 126784*a^44 + 53088*a^42 - 28924*a^40 + 14280*a^38 + 98560*a^36 - 211680*a^34 + 105840*a^32 + 99680*a^30 - 80920*a^28 - 103600*a^26 + 239400*a^24 - 290640*a^22 + 178080*a^20 + 47040*a^18 - 127050*a^16 + 139440*a^14 - 270480*a^12 + 282240*a^10 - 45360*a^8 - 141120*a^6 + 120960*a^4 - 40320*a^2 + 5040)/a^56)*t,\n",
       " ((362880*a^72 - 1451520*a^70 + 1905120*a^68 - 483840*a^66 - 793800*a^64 + 453600*a^62 - 20160*a^60 + 90720*a^58 - 45360*a^56 - 28560*a^54 + 22680*a^52 - 18144*a^50 + 4158*a^48 + 3024*a^46 - 2016*a^42 + 1260*a^40 - 168*a^36 + 216*a^30 - 72*a^28 - 18*a^16 + 1)/a^72)*t^9 + ((1451520*a^74 - 6713280*a^72 + 10523520*a^70 - 3538080*a^68 - 7166880*a^66 + 7832160*a^64 - 3144960*a^62 + 1753920*a^60 - 1103760*a^58 - 378000*a^56 + 910224*a^54 - 675864*a^52 + 260064*a^50 + 28728*a^48 - 9072*a^46 - 80640*a^44 + 70560*a^42 - 17136*a^40 - 6048*a^38 + 3024*a^36 - 4536*a^34 + 10656*a^32 - 7704*a^30 + 1584*a^28 + 504*a^20 - 1008*a^18 + 504*a^16 - 36*a^4 + 72*a^2 - 36)/a^72)*t^8 + ((423360*a^78 + 997920*a^76 - 12428640*a^74 + 28067760*a^72 - 19020960*a^70 - 13456800*a^68 + 27311760*a^66 - 16579080*a^64 + 10236240*a^62 - 8888040*a^60 + 1424304*a^58 + 6095628*a^56 - 7617456*a^54 + 4496184*a^52 - 619920*a^50 - 479052*a^48 - 586656*a^46 + 1117368*a^44 - 569016*a^42 + 20664*a^40 + 58968*a^38 - 61488*a^36 + 179928*a^34 - 198072*a^32 + 89208*a^30 - 15120*a^28 + 2268*a^24 + 9072*a^22 - 25704*a^20 + 21168*a^18 - 5796*a^16 + 84*a^12 - 126*a^8 - 1008*a^6 + 2520*a^4 - 2016*a^2 + 546)/a^72)*t^7 + ((90720*a^84 - 151200*a^82 + 952560*a^80 - 2509920*a^78 - 8550360*a^76 + 40642560*a^74 - 51030000*a^72 + 4974480*a^70 + 40786200*a^68 - 38858400*a^66 + 28963872*a^64 - 31948560*a^62 + 20853000*a^60 + 5174064*a^58 - 23848272*a^56 + 23180472*a^54 - 7975800*a^52 - 3795120*a^50 + 2014488*a^48 + 4458384*a^46 - 5311152*a^44 + 2057832*a^42 - 352548*a^40 + 287280*a^38 + 630000*a^36 - 1680840*a^34 + 1360800*a^32 - 461160*a^30 - 2520*a^28 + 136080*a^26 - 64386*a^24 - 196560*a^22 + 315000*a^20 - 175392*a^18 + 34776*a^16 - 2016*a^14 + 7560*a^12 - 15120*a^10 + 2646*a^8 + 31248*a^6 - 43344*a^4 + 23184*a^2 - 4536)/a^72)*t^6 + ((15120*a^92 - 18144*a^90 + 2268*a^88 + 154224*a^86 - 443520*a^84 + 804384*a^82 - 4406220*a^80 + 4147920*a^78 + 26111232*a^76 - 67208400*a^74 + 48271860*a^72 + 16273656*a^70 - 43285032*a^68 + 39589200*a^66 - 49826070*a^64 + 52552080*a^62 - 24635520*a^60 - 16619400*a^58 + 41954472*a^56 - 30663360*a^54 - 1716120*a^52 + 15054984*a^50 - 300258*a^48 - 14649264*a^46 + 12370680*a^44 - 4879224*a^42 + 3877146*a^40 - 2819880*a^38 - 2894472*a^36 + 6156360*a^34 - 3891510*a^32 + 703584*a^30 + 1073772*a^28 - 1658160*a^26 + 479010*a^24 + 1335600*a^22 - 1634220*a^20 + 781200*a^18 - 211617*a^16 + 136080*a^14 - 253680*a^12 + 287280*a^10 - 6930*a^8 - 327600*a^6 + 335160*a^4 - 141120*a^2 + 22449)/a^72)*t^5 + ((2016*a^102 - 1512*a^100 + 168*a^96 + 18144*a^94 - 65016*a^92 + 65016*a^90 + 51912*a^88 - 555408*a^86 + 715680*a^84 - 2029104*a^82 + 8481816*a^80 - 1542576*a^78 - 38264184*a^76 + 62637624*a^74 - 26911584*a^72 - 15769656*a^70 + 24697008*a^68 - 32813928*a^66 + 48401640*a^64 - 47118960*a^62 + 22057140*a^60 + 15539832*a^58 - 35120232*a^56 + 14938056*a^54 + 17504424*a^52 - 19161576*a^50 - 7819308*a^48 + 24211152*a^46 - 15733116*a^44 + 7275912*a^42 - 9064944*a^40 + 4148928*a^38 + 7590996*a^36 - 10214064*a^34 + 4367034*a^32 + 1021104*a^30 - 5400864*a^28 + 6758640*a^26 - 1948170*a^24 - 3391920*a^22 + 3877776*a^20 - 2057832*a^18 + 1030806*a^16 - 1199520*a^14 + 2076480*a^12 - 1890000*a^10 - 58590*a^8 + 1597680*a^6 - 1355004*a^4 + 487368*a^2 - 67284)/a^72)*t^4 + ((216*a^114 - 72*a^112 + 1512*a^104 - 6048*a^102 + 1872*a^100 + 4032*a^98 + 3654*a^96 - 54432*a^94 + 104328*a^92 - 107688*a^90 - 147420*a^88 + 734832*a^86 - 503964*a^84 + 2727648*a^82 - 8302644*a^80 - 1993824*a^78 + 29496096*a^76 - 34321896*a^74 + 10571736*a^72 + 5004720*a^70 - 8072064*a^68 + 14283864*a^66 - 22345092*a^64 + 25863264*a^62 - 16403688*a^60 - 3900960*a^58 + 12075084*a^56 + 2124192*a^54 - 17926272*a^52 + 12195792*a^50 + 9939321*a^48 - 20431152*a^46 + 11304720*a^44 - 4447296*a^42 + 4994766*a^40 + 1586592*a^38 - 9263072*a^36 + 5885712*a^34 - 359478*a^32 - 2390688*a^30 + 8754300*a^28 - 12358080*a^26 + 6247122*a^24 + 1014048*a^22 - 3329676*a^20 + 3264912*a^18 - 2712069*a^16 + 3976560*a^14 - 6822564*a^12 + 5518800*a^10 + 324576*a^8 - 3904992*a^6 + 2928240*a^4 - 945504*a^2 + 118124)/a^72)*t^3 + ((18*a^128 + 72*a^116 - 432*a^114 + 63*a^112 + 168*a^108 - 2772*a^104 + 6048*a^102 - 360*a^100 - 6912*a^98 - 7980*a^96 + 54432*a^94 - 77112*a^92 + 89376*a^90 + 139068*a^88 - 425808*a^86 + 161172*a^84 - 1804320*a^82 + 4068918*a^80 + 1949808*a^78 - 11683224*a^76 + 10662624*a^74 - 2992668*a^72 - 458640*a^70 + 950544*a^68 - 1688400*a^66 + 3778236*a^64 - 6995520*a^62 + 6626004*a^60 - 1819440*a^58 - 458136*a^56 - 2590896*a^54 + 5779368*a^52 - 3877272*a^50 - 2059155*a^48 + 5085360*a^46 - 2657340*a^44 - 127680*a^42 + 2014488*a^40 - 4914000*a^38 + 3886260*a^36 + 1572480*a^34 - 1814850*a^32 - 930240*a^30 - 3715560*a^28 + 10417680*a^26 - 10436244*a^24 + 6189120*a^22 - 851760*a^20 - 2876328*a^18 + 3322494*a^16 - 5511744*a^14 + 9497880*a^12 - 7166880*a^10 - 533736*a^8 + 4540032*a^6 - 3137616*a^4 + 940896*a^2 - 109584)/a^72)*t^2 + ((a^144 - 18*a^128 - 72*a^116 + 216*a^114 + 9*a^112 - 168*a^108 + 1260*a^104 - 2016*a^102 + 2880*a^98 + 4158*a^96 - 18144*a^94 + 22680*a^92 - 28560*a^90 - 45828*a^88 + 92160*a^86 - 20088*a^84 + 452592*a^82 - 794430*a^80 - 474768*a^78 + 1892520*a^76 - 1435392*a^74 + 373296*a^72 - 75600*a^70 + 13104*a^68 - 173376*a^66 + 568134*a^64 + 102816*a^62 - 1342656*a^60 + 1214640*a^58 - 275184*a^56 - 252672*a^54 + 491400*a^52 - 38808*a^50 - 1328922*a^48 + 1918224*a^46 - 1010520*a^44 + 620928*a^42 - 1473696*a^40 + 1658160*a^38 + 108920*a^36 - 1895040*a^34 + 525420*a^32 + 1975680*a^30 - 695520*a^28 - 3296160*a^26 + 5720400*a^24 - 4959360*a^22 + 1648080*a^20 + 1043280*a^18 - 1459080*a^16 + 2600640*a^14 - 4505760*a^12 + 3265920*a^10 + 272160*a^8 - 1935360*a^6 + 1270080*a^4 - 362880*a^2 + 40320)/a^72)*t,\n",
       " ((3628800*a^90 - 16329600*a^88 + 25401600*a^86 - 11037600*a^84 - 9298800*a^82 + 8958600*a^80 - 1058400*a^78 + 756000*a^76 - 982800*a^74 - 361200*a^72 + 483840*a^70 - 302400*a^68 + 115920*a^66 + 51030*a^64 - 15120*a^62 - 25200*a^60 + 22680*a^58 - 3780*a^56 - 5040*a^54 + 252*a^50 + 3300*a^48 - 2160*a^46 + 240*a^42 - 270*a^34 + 90*a^32 + 20*a^18 - 1)/a^90)*t^10 + ((16329600*a^92 - 83462400*a^90 + 151048800*a^88 - 79833600*a^86 - 86070600*a^84 + 133131600*a^82 - 62785800*a^80 + 24796800*a^78 - 17766000*a^76 - 3477600*a^74 + 16254000*a^72 - 12882240*a^70 + 5445090*a^68 + 147420*a^66 - 708750*a^64 - 907200*a^62 + 1168020*a^60 - 400680*a^58 - 117180*a^56 + 85680*a^54 - 59220*a^52 + 181080*a^50 - 164340*a^48 + 41760*a^46 + 11520*a^44 - 5760*a^42 + 7560*a^38 - 17730*a^36 + 12780*a^34 - 2610*a^32 - 720*a^22 + 1440*a^20 - 720*a^18 + 45*a^4 - 90*a^2 + 45)/a^90)*t^9 + ((4838400*a^96 + 12700800*a^94 - 169192800*a^92 + 427744800*a^90 - 364316400*a^88 - 148024800*a^86 + 484394400*a^84 - 342165600*a^82 + 168210000*a^80 - 133282800*a^78 + 30678480*a^76 + 106747200*a^74 - 144573660*a^72 + 91037520*a^70 - 18045720*a^68 - 11990160*a^66 - 2158380*a^64 + 16077600*a^62 - 10694880*a^60 + 713160*a^58 + 1605240*a^56 - 1039680*a^54 + 2968560*a^52 - 3935160*a^50 + 2014920*a^48 - 170640*a^46 - 208080*a^44 + 91620*a^42 + 127260*a^40 - 441000*a^38 + 461160*a^36 - 202860*a^34 + 31500*a^32 + 1680*a^30 - 2520*a^26 - 20160*a^24 + 50400*a^22 - 40320*a^20 + 10920*a^18 - 120*a^12 + 90*a^8 + 1800*a^6 - 4140*a^4 + 3240*a^2 - 870)/a^90)*t^8 + ((1058400*a^102 - 2268000*a^100 + 13834800*a^98 - 35708400*a^96 - 130032000*a^94 + 669286800*a^92 - 956503800*a^90 + 219183300*a^88 + 770742000*a^86 - 878371200*a^84 + 579670560*a^82 - 540974700*a^80 + 373605120*a^78 + 111173580*a^76 - 491674680*a^74 + 489984600*a^72 - 201625200*a^70 - 57432060*a^68 + 73788120*a^66 + 42348600*a^64 - 86441040*a^62 + 39506040*a^60 - 3693060*a^58 + 3093300*a^56 + 9945600*a^54 - 32460120*a^52 + 29657880*a^50 - 9565920*a^48 - 2255400*a^46 + 4000500*a^44 - 1936620*a^42 - 2835000*a^40 + 6116040*a^38 - 4570020*a^36 + 1564920*a^34 - 245700*a^32 + 151200*a^30 - 302400*a^28 + 52920*a^26 + 625170*a^24 - 866880*a^22 + 463680*a^20 - 92400*a^18 + 1260*a^16 + 2520*a^14 - 13650*a^12 + 26460*a^10 + 5040*a^8 - 75600*a^6 + 95130*a^4 - 49140*a^2 + 9450)/a^90)*t^7 + ((181440*a^110 - 302400*a^108 + 90720*a^106 + 2328480*a^104 - 7121520*a^102 + 15951600*a^100 - 76408920*a^98 + 67677120*a^96 + 490991760*a^94 - 1357454700*a^92 + 1117141704*a^90 + 305203500*a^88 - 1129962960*a^86 + 1018476900*a^84 - 1058853600*a^82 + 1097998650*a^80 - 470594880*a^78 - 483349860*a^76 + 1024274160*a^74 - 785014020*a^72 + 84929040*a^70 + 311647140*a^68 - 113053500*a^66 - 214680690*a^64 + 235869480*a^62 - 96897780*a^60 + 64474200*a^58 - 50062320*a^56 - 65071440*a^54 + 145028520*a^52 - 92493954*a^50 + 5773950*a^48 + 34518960*a^46 - 39545100*a^44 + 16914660*a^42 + 20726748*a^40 - 36300600*a^38 + 23713200*a^36 - 8562330*a^34 + 3509730*a^32 - 5073600*a^30 + 5747490*a^28 - 147420*a^26 - 6581190*a^24 + 6788880*a^22 - 2871540*a^20 + 358260*a^18 + 220815*a^16 - 393120*a^14 + 734370*a^12 - 706860*a^10 - 205380*a^8 + 1108800*a^6 - 1015560*a^4 + 408240*a^2 - 63273)/a^90)*t^6 + ((25200*a^120 - 30240*a^118 + 3780*a^116 + 5040*a^114 + 302400*a^112 - 1088640*a^110 + 1330140*a^108 + 1043280*a^106 - 10417680*a^104 + 15975120*a^102 - 50712480*a^100 + 179177040*a^98 - 10973340*a^96 - 916753320*a^94 + 1590815520*a^92 - 753738720*a^90 - 526027320*a^88 + 894055680*a^86 - 1008878850*a^84 + 1353258900*a^82 - 1226751750*a^80 + 371038080*a^78 + 666750420*a^76 - 1103540760*a^74 + 569337720*a^72 + 300086640*a^70 - 497530530*a^68 - 55778520*a^66 + 493431750*a^64 - 377435520*a^62 + 195380430*a^60 - 207997020*a^58 + 57365280*a^56 + 260902320*a^54 - 323858430*a^52 + 110111400*a^50 + 77649390*a^48 - 160884360*a^46 + 170528400*a^44 - 72934680*a^42 - 59821020*a^40 + 104459040*a^38 - 70473270*a^36 + 32365620*a^34 - 25507440*a^32 + 41584200*a^30 - 37878750*a^28 - 1845900*a^26 + 33829950*a^24 - 29405880*a^22 + 10409805*a^20 + 518070*a^18 - 3450825*a^16 + 6073200*a^14 - 9507750*a^12 + 6868260*a^10 + 2280600*a^8 - 7938000*a^6 + 5949405*a^4 - 2020410*a^2 + 269325)/a^90)*t^5 + ((2880*a^132 - 2160*a^130 + 240*a^126 + 30240*a^122 - 108360*a^120 + 58320*a^118 + 83880*a^116 + 94920*a^114 - 1216530*a^112 + 2474640*a^110 - 2546040*a^108 - 4480560*a^106 + 17745840*a^104 - 17926560*a^102 + 89705700*a^100 - 225077580*a^98 - 102812220*a^96 + 950949720*a^94 - 1146612600*a^92 + 344632680*a^90 + 297105480*a^88 - 450339120*a^86 + 697075560*a^84 - 969506370*a^82 + 847885500*a^80 - 245377440*a^78 - 455038920*a^76 + 587190240*a^74 - 12825960*a^72 - 553301280*a^70 + 387986760*a^68 + 274208340*a^66 - 573092100*a^64 + 371946960*a^62 - 251049960*a^60 + 204301440*a^58 + 159769260*a^56 - 506035600*a^54 + 348683580*a^52 + 15924510*a^50 - 205876905*a^48 + 319384440*a^46 - 361328940*a^44 + 193651680*a^42 + 40795020*a^40 - 140779800*a^38 + 116003090*a^36 - 64270440*a^34 + 79073100*a^32 - 141734880*a^30 + 110731950*a^28 + 22726620*a^26 - 107185890*a^24 + 84218400*a^22 - 25288830*a^20 - 8278220*a^18 + 18807075*a^16 - 36061200*a^14 + 52325070*a^12 - 31997700*a^10 - 11373390*a^8 + 30706200*a^6 - 20085660*a^4 + 6055560*a^2 - 723680)/a^90)*t^4 + ((270*a^146 - 90*a^144 + 2160*a^134 - 8640*a^132 + 2700*a^130 + 1440*a^128 + 4320*a^126 + 3780*a^124 - 83160*a^122 + 161400*a^120 - 16560*a^118 - 208350*a^116 - 336840*a^114 + 1815075*a^112 - 2918160*a^110 + 2732100*a^108 + 6480720*a^106 - 14385240*a^104 + 12169440*a^102 - 88373880*a^100 + 158487120*a^98 + 127531320*a^96 - 565301520*a^94 + 517699980*a^92 - 129051720*a^90 - 70975440*a^88 + 131737320*a^86 - 255179940*a^84 + 381764250*a^82 - 339612840*a^80 + 57830640*a^78 + 201981780*a^76 - 115814160*a^74 - 225781500*a^72 + 402161760*a^70 - 163426410*a^68 - 208797540*a^66 + 296050860*a^64 - 186288480*a^62 + 137868150*a^60 + 6237000*a^58 - 341570880*a^56 + 449344560*a^54 - 148483755*a^52 - 116808570*a^50 + 162964740*a^48 - 258951600*a^46 + 365664420*a^44 - 264210780*a^42 + 55131300*a^40 + 78568560*a^38 - 90511380*a^36 + 52013880*a^34 - 116183250*a^32 + 227703000*a^30 - 148664250*a^28 - 85903020*a^26 + 219157680*a^24 - 166740120*a^22 + 44951895*a^20 + 23729250*a^18 - 46824435*a^16 + 99159480*a^14 - 138125400*a^12 + 75932640*a^10 + 28105560*a^8 - 64562400*a^6 + 38408220*a^4 - 10631160*a^2 + 1172700)/a^90)*t^3 + ((20*a^162 + 90*a^148 - 540*a^146 + 80*a^144 + 240*a^138 - 4320*a^134 + 9060*a^132 - 288*a^130 - 1260*a^128 - 9600*a^126 - 7560*a^124 + 75600*a^122 - 103440*a^120 - 26640*a^118 + 172395*a^116 + 350730*a^114 - 1203435*a^112 + 1834560*a^110 - 1573320*a^108 - 4116960*a^106 + 5472000*a^104 - 5192400*a^102 + 44654220*a^100 - 59340060*a^98 - 61632180*a^96 + 183027600*a^94 - 137424960*a^92 + 33509736*a^90 + 5580720*a^88 - 14691960*a^86 + 39786750*a^84 - 72578520*a^82 + 55157634*a^80 + 27605760*a^78 - 74004840*a^76 + 12093480*a^74 + 88400340*a^72 - 108954720*a^70 + 26595450*a^68 + 48441540*a^66 - 38671920*a^64 + 5894280*a^62 + 13095600*a^60 - 83764800*a^58 + 178559640*a^56 - 136077440*a^54 - 25746525*a^52 + 72496422*a^50 + 2331135*a^48 + 22228200*a^46 - 123725520*a^44 + 122082600*a^42 - 35544348*a^40 - 23045400*a^38 + 15367450*a^36 - 2148300*a^34 + 80946180*a^32 - 167185200*a^30 + 79929360*a^28 + 127867320*a^26 - 247429560*a^24 + 189720720*a^22 - 49421610*a^20 - 27509580*a^18 + 53661510*a^16 - 124120080*a^14 + 168826680*a^12 - 86592240*a^10 - 32874120*a^8 + 67975200*a^6 - 37862640*a^4 + 9862560*a^2 - 1026576)/a^90)*t^2 + ((a^180 - 20*a^162 - 90*a^148 + 270*a^146 + 10*a^144 - 240*a^138 + 2160*a^134 - 3300*a^132 - 252*a^130 - 180*a^128 + 5040*a^126 + 3780*a^124 - 22680*a^122 + 25200*a^120 + 15120*a^118 - 51705*a^116 - 113850*a^114 + 302490*a^112 - 483840*a^110 + 359520*a^108 + 982800*a^106 - 743400*a^104 + 1037520*a^102 - 8957160*a^100 + 9327600*a^98 + 11079300*a^96 - 25583040*a^94 + 16553160*a^92 - 3901080*a^90 - 473040*a^88 + 915840*a^86 - 195420*a^84 + 4577580*a^82 - 8085294*a^80 - 4562880*a^78 + 18819360*a^76 - 14815080*a^74 + 4579680*a^72 - 1935360*a^70 + 5062680*a^68 - 7081620*a^66 - 2570400*a^64 + 21299040*a^62 - 28350420*a^60 + 20107080*a^58 - 8638560*a^56 - 12048960*a^54 + 33927390*a^52 - 15133860*a^50 - 35130270*a^48 + 46090800*a^46 - 15397200*a^44 + 6347040*a^42 - 18579960*a^40 + 11415600*a^38 + 10027500*a^36 - 10773000*a^34 - 21621600*a^32 + 44553600*a^30 - 9563400*a^28 - 62748000*a^26 + 107604000*a^24 - 83764800*a^22 + 21795480*a^20 + 11264400*a^18 - 22415400*a^16 + 55339200*a^14 - 74239200*a^12 + 36469440*a^10 + 14061600*a^8 - 27216000*a^6 + 14515200*a^4 - 3628800*a^2 + 362880)/a^90)*t]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EGF of labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "# we take the 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": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " t,\n",
       " 3/2*t^2 + 1/2*t,\n",
       " 25/8*t^3 + 21/8*t^2 + 9/4*t,\n",
       " 543/64*t^4 + 387/32*t^3 + 1173/64*t^2 + 803/32*t,\n",
       " 29281/1024*t^5 + 29915/512*t^4 + 122595/1024*t^3 + 135895/512*t^2 + 70635/128*t,\n",
       " 3781503/32768*t^6 + 10110735/32768*t^5 + 25349355/32768*t^4 + 70500705/32768*t^3 + 114612375/16384*t^2 + 91846847/4096*t,\n",
       " 1138779265/2097152*t^7 + 3767987307/2097152*t^6 + 11002159455/2097152*t^5 + 35130437825/2097152*t^4 + 17265271485/262144*t^3 + 171448339951/524288*t^2 + 220193225973/131072*t,\n",
       " 783702329343/268435456*t^8 + 773027042823/67108864*t^7 + 5109415147935/134217728*t^6 + 4499199458565/33554432*t^5 + 156618999659535/268435456*t^4 + 231722601662643/67108864*t^3 + 1851196751372157/67108864*t^2 + 3969950586854011/16777216*t,\n",
       " 1213442454842881/68719476736*t^9 + 1390244741807607/17179869184*t^8 + 10210204164239793/34359738368*t^7 + 9713471222508417/8589934592*t^6 + 360041417744945265/68719476736*t^5 + 578921222060374683/17179869184*t^4 + 5517968688758637531/17179869184*t^3 + 18568937552061936879/4294967296*t^2 + 34378210768563301635/536870912*t,\n",
       " 4175098976430598143/35184372088832*t^10 + 21795313802880245805/35184372088832*t^9 + 43893098221016031885/17592186044416*t^8 + 178285533097741095285/17592186044416*t^7 + 1726485549525658743255/35184372088832*t^6 + 11581243315062995535885/35184372088832*t^5 + 3740600636087000673645/1099511627776*t^4 + 484767599023193364298935/8796093022208*t^3 + 2839611643411126665349695/2199023255552*t^2 + 9300263322795427228870307/274877906944*t]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DTsub = [sum(dt[i][0][k].subs(a=sqrt(2)) * i.factorial() * t^k for k in srange(N)) for i in srange(disp)]\n",
    "DTsub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, 2, 8, 64, 1024, 32768, 2097152, 268435456, 68719476736, 35184372088832]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The numbers of labeled digraphs\n",
    "DTsubsub = [DTsub[i].subs(t=1) for i in srange(disp)]\n",
    "DTsubsub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "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": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CGF of labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "# we take the parameter a = sqrt{2}\n",
    "num = t*exp((1-t)*scd.subs(z=a^3*z^2*w))*exp_had_prod(PhiTwoOneTwo,(1-it.subs(z=a^2*z*w))*(1-it.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(disps)]\n",
    "QDT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[t,\n",
       " 4*t^2*w - 4*t*w,\n",
       " 20*t^3*w^2 - 28*t^2*w^2 - 4*t^2*w + 8*t*w^2 + 4*t*w,\n",
       " 122*t^4*w^3 - 174*t^3*w^3 - 48*t^3*w^2 + 84*t^2*w^3 + 72*t^2*w^2 - 32*t*w^3 - 24*t*w^2,\n",
       " 1775/2*t^5*w^4 - 1085*t^4*w^4 - 480*t^4*w^3 + 1605/2*t^3*w^4 + 768*t^3*w^3 - 93*t^2*w^4 + 12*t^3*w^2 - 336*t^2*w^3 - 512*t*w^4 - 36*t^2*w^2 + 48*t*w^3 + 24*t*w^2]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CGF of labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "QDTsub = [sum(sum(QDT[i][k][j].subs(a=sqrt(2)) * t^j * w^k for j in srange(N)) for k in srange(N)) for i in srange(disps)]\n",
    "QDTsub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many terms we want to display\n",
    "disp = 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[t,\n",
       " (2*a^2*t^2 - 2*a^2*t)*w,\n",
       " ((3*a^6 - a^4)*t^3 + (a^8 - 6*a^6 + a^4)*t^2 + (-a^8 + 3*a^6)*t)*w^2 + (-a^4*t^2 + a^4*t)*w,\n",
       " ((4*a^12 - 3*a^10 + 1/3*a^6)*t^4 + (3*a^14 - 12*a^12 + 4*a^10 + 2*a^8 - a^6)*t^3 + (1/3*a^18 - 6*a^14 + 12*a^12 - a^10 - 2*a^8 + 2/3*a^6)*t^2 + (-1/3*a^18 + 3*a^14 - 4*a^12)*t)*w^3 + (-2*a^10*t^3 + (2*a^10 + 2*a^8)*t^2 - 2*a^8*t)*w^2,\n",
       " ((5*a^20 - 6*a^18 + 3/4*a^16 + a^14 - 1/12*a^8)*t^5 + (6*a^22 - 43/2*a^20 + 12*a^18 + 9/2*a^16 - 11/3*a^14 + 1/2*a^12 - a^10 + 1/2*a^8)*t^4 + (a^26 + 3/4*a^24 - 18*a^22 + 193/6*a^20 - 7*a^18 - 15/2*a^16 + 3*a^14 - 5/2*a^12 + 3*a^10 - 11/12*a^8)*t^3 + (1/12*a^32 - 2*a^26 - 3/2*a^24 + 18*a^22 - 62/3*a^20 + a^18 + 9/4*a^16 - 1/3*a^14 + 2*a^12 - 2*a^10 + 1/2*a^8)*t^2 + (-1/12*a^32 + a^26 + 3/4*a^24 - 6*a^22 + 5*a^20)*t)*w^4 + ((-3*a^18 + a^16)*t^4 + (-a^20 + 5*a^18 + 2*a^16)*t^3 + (a^20 - 2*a^18 - 2*a^16 - 3*a^14)*t^2 + (-a^16 + 3*a^14)*t)*w^3 + (1/2*a^12*t^3 + (-1/2*a^16 + a^14 - 3/2*a^12)*t^2 + (1/2*a^16 - a^14 + a^12)*t)*w^2,\n",
       " ((6*a^30 - 10*a^28 + 3*a^26 + 2*a^24 - 1/2*a^22 - 1/4*a^18 + 1/60*a^10)*t^6 + (10*a^32 - 36*a^30 + 28*a^28 + 13/2*a^26 - 12*a^24 + 7/2*a^22 - 3*a^20 + 5/3*a^18 - 1/6*a^14 + 1/3*a^12 - 1/6*a^10)*t^5 + (2*a^36 + 5/2*a^34 - 40*a^32 + 145/2*a^30 - 24*a^28 - 283/12*a^26 + 17*a^24 - 65/6*a^22 + 11*a^20 - 9/2*a^18 + 2*a^14 - 2*a^12 + 7/12*a^10)*t^4 + (1/4*a^42 + 1/2*a^38 - 6*a^36 - 33/4*a^34 + 60*a^32 - 407/6*a^30 + 17/3*a^28 + 211/12*a^26 - 29/3*a^24 + 29/2*a^22 - 15*a^20 + 91/12*a^18 - 29/6*a^14 + 11/3*a^12 - 5/6*a^10)*t^3 + (1/60*a^50 - 1/2*a^42 - a^38 + 6*a^36 + 35/4*a^34 - 40*a^32 + 94/3*a^30 + 1/3*a^28 - 7/2*a^26 + 8/3*a^24 - 20/3*a^22 + 7*a^20 - 9/2*a^18 + 3*a^14 - 2*a^12 + 2/5*a^10)*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)*t)*w^5 + ((-4*a^28 + 3*a^26 - 1/3*a^22)*t^5 + (-3*a^30 + 10*a^28 + a^26 - 4*a^24 + 4/3*a^22)*t^4 + (-1/3*a^34 + 5*a^30 - 6*a^28 - 7*a^26 + 2*a^24 - 5/3*a^22)*t^3 + (1/3*a^34 - 2*a^30 + 1/3*a^28 + 3*a^26 - a^24 + 14/3*a^22)*t^2 + (-1/3*a^28 + 3*a^24 - 4*a^22)*t)*w^4 + (a^22*t^4 + (-a^26 + 2*a^24 - 3*a^22 - a^18)*t^3 + (a^26 - 2*a^24 + 3*a^22 - 2*a^20 + 3*a^18)*t^2 + (-a^22 + 2*a^20 - 2*a^18)*t)*w^3,\n",
       " ((7*a^42 - 15*a^40 + 15/2*a^38 + 17/6*a^36 - 2*a^34 - 5/12*a^30 + 1/8*a^28 + 1/20*a^22 - 1/360*a^12)*t^7 + (15*a^44 - 57*a^42 + 113/2*a^40 + 4*a^38 - 57/2*a^36 + 23/2*a^34 - 47/8*a^32 + 19/4*a^30 - 7/8*a^28 - 1/2*a^26 + a^24 - 8/15*a^22 + 1/24*a^16 - 1/12*a^14 + 1/24*a^12)*t^6 + (10/3*a^48 + 11/2*a^46 - 153/2*a^44 + 889/6*a^42 - 269/4*a^40 - 53*a^38 + 1303/24*a^36 - 28*a^34 + 85/3*a^32 - 209/12*a^30 + 17/8*a^28 + 19/3*a^26 - 121/18*a^24 + 25/12*a^22 + 5/24*a^20 + 1/6*a^18 - 11/12*a^16 + 5/6*a^14 - 17/72*a^12)*t^5 + (1/2*a^54 - 1/8*a^52 + 2*a^50 - 77/6*a^48 - 103/4*a^46 + 617/4*a^44 - 3305/18*a^42 + 197/8*a^40 + 273/4*a^38 - 1075/24*a^36 + 277/6*a^34 - 437/8*a^32 + 199/6*a^30 - 53/24*a^28 - 52/3*a^26 + 107/8*a^24 - 11/12*a^22 - 11/4*a^20 - 4/3*a^18 + 107/24*a^16 - 35/12*a^14 + 5/8*a^12)*t^4 + (1/20*a^62 + 1/8*a^56 - 17/12*a^54 + 23/120*a^52 - 6*a^50 + 37/2*a^48 + 169/4*a^46 - 461/3*a^44 + 368/3*a^42 + 10/3*a^40 - 385/12*a^38 + 194/9*a^36 - 110/3*a^34 + 511/12*a^32 - 121/4*a^30 + 4*a^28 + 91/6*a^26 - 653/72*a^24 - 323/60*a^22 + 163/24*a^20 + 17/6*a^18 - 91/12*a^16 + 25/6*a^14 - 137/180*a^12)*t^3 + (1/360*a^72 - 1/10*a^62 - 1/4*a^56 + 4/3*a^54 - 1/15*a^52 + 6*a^50 - 71/6*a^48 - 59/2*a^46 + 911/12*a^44 - 398/9*a^42 - 53/24*a^40 + 16/3*a^38 - 97/18*a^36 + 9*a^34 - 125/12*a^32 + 61/6*a^30 - 19/6*a^28 - 11/3*a^26 + 17/12*a^24 + 47/10*a^22 - 17/4*a^20 - 5/3*a^18 + 4*a^16 - 2*a^14 + 1/3*a^12)*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)*t)*w^6 + ((-5*a^40 + 6*a^38 - 3/4*a^36 - a^34 + 1/12*a^28)*t^6 + (-6*a^42 + 37/2*a^40 - 4*a^38 - 39/4*a^36 + 4*a^34 + 1/6*a^32 + a^30 - 7/12*a^28)*t^5 + (-a^46 - 3/4*a^44 + 15*a^42 - 53/3*a^40 - 11*a^38 + 8*a^36 + 3*a^34 + 1/3*a^32 - 4*a^30 + 17/12*a^28)*t^4 + (-1/12*a^52 + 5/3*a^46 + 23/12*a^44 - 12*a^42 + 7/2*a^40 + 8*a^38 + 29/4*a^36 - 6*a^34 - 7/6*a^32 + 5*a^30 - 17/12*a^28)*t^3 + (1/12*a^52 - 2/3*a^46 - 13/12*a^44 + 3*a^42 + 2/3*a^40 - 11/2*a^36 + 6*a^34 - 13/3*a^32 - 2*a^30 + 1/2*a^28)*t^2 + (-1/12*a^44 + a^38 + 3/4*a^36 - 6*a^34 + 5*a^32)*t)*w^5 + ((3/2*a^34 - 1/2*a^32)*t^5 + (-3/2*a^38 + 4*a^36 - 13/2*a^34 + 2*a^32 - 2*a^30)*t^4 + (-1/2*a^40 + 7/2*a^38 - 15/2*a^36 + 10*a^34 - 13/2*a^32 + 6*a^30 - 1/2*a^28 + 3/2*a^26)*t^3 + (1/2*a^40 - 2*a^38 + 7/2*a^36 - 5*a^34 + 11/2*a^32 - 13/2*a^30 + 9/2*a^28 - 9/2*a^26)*t^2 + (-1/2*a^32 + 5/2*a^30 - 4*a^28 + 3*a^26)*t)*w^4 + (-1/6*a^24*t^4 + (1/2*a^28 - a^26 + a^24)*t^3 + (-1/6*a^36 + a^32 - a^30 - 3/2*a^28 + 3*a^26 - 11/6*a^24)*t^2 + (1/6*a^36 - a^32 + a^30 + a^28 - 2*a^26 + a^24)*t)*w^3,\n",
       " ((8*a^56 - 21*a^54 + 15*a^52 + 5/2*a^50 - 5*a^48 + 1/2*a^46 - 1/2*a^44 + 1/2*a^42 - 1/24*a^38 + 1/10*a^36 - 1/40*a^34 - 1/120*a^26 + 1/2520*a^14)*t^8 + (21*a^58 - 86*a^56 + 207/2*a^54 - 10*a^52 - 167/3*a^50 + 32*a^48 - 12*a^46 + 65/6*a^44 - 27/8*a^42 - 7/4*a^40 + 53/20*a^38 - 33/20*a^36 + 11/40*a^34 + 1/8*a^30 - 1/4*a^28 + 47/360*a^26 - 1/120*a^18 + 1/60*a^16 - 1/120*a^14)*t^7 + (5*a^62 + 10*a^60 - 133*a^58 + 841/3*a^56 - 327/2*a^54 - 195/2*a^52 + 3475/24*a^50 - 149/2*a^48 + 509/8*a^46 - 139/3*a^44 + 21/4*a^42 + 75/4*a^40 - 799/40*a^38 + 171/20*a^36 - 1/2*a^34 + 1/2*a^32 - 17/6*a^30 + 8/3*a^28 - 7/9*a^26 - 1/24*a^22 - 1/12*a^20 + 7/24*a^18 - 1/4*a^16 + 5/72*a^14)*t^6 + (5/6*a^68 - 1/2*a^66 + 5*a^64 - 563/24*a^62 - 121/2*a^60 + 2719/8*a^58 - 1313/3*a^56 + 345/4*a^54 + 203*a^52 - 3013/18*a^50 + 545/4*a^48 - 1251/8*a^46 + 1139/12*a^44 + 49/24*a^42 - 697/12*a^40 + 1787/36*a^38 - 11*a^36 - 157/24*a^34 - 38/9*a^32 + 365/24*a^30 - 43/4*a^28 + 3*a^26 - 5/4*a^24 + 11/12*a^22 + 7/6*a^20 - 19/8*a^18 + 17/12*a^16 - 7/24*a^14)*t^5 + (1/10*a^76 - 1/40*a^74 + 1/2*a^70 - 3*a^68 + 63/40*a^66 - 79/4*a^64 + 757/18*a^62 + 541/4*a^60 - 3605/8*a^58 + 1142/3*a^56 + 823/60*a^54 - 1913/12*a^52 + 123*a^50 - 331/2*a^48 + 2213/12*a^46 - 219/2*a^44 - 77/24*a^42 + 437/6*a^40 - 3553/72*a^38 - 373/30*a^36 + 503/20*a^34 + 38/3*a^32 - 139/4*a^30 + 259/12*a^28 - 3167/360*a^26 + 17/2*a^24 - 107/24*a^22 - 59/12*a^20 + 185/24*a^18 - 15/4*a^16 + 29/45*a^14)*t^4 + (1/120*a^86 + 1/40*a^78 - 3/10*a^76 + 29/360*a^74 - 3/2*a^70 + 4*a^68 - 53/24*a^66 + 1777/60*a^64 - 2617/72*a^62 - 144*a^60 + 39619/120*a^58 - 17611/90*a^56 - 1409/60*a^54 + 227/4*a^52 - 2113/36*a^50 + 1151/12*a^48 - 1231/12*a^46 + 1217/18*a^44 - 95/24*a^42 - 439/12*a^40 + 4381/360*a^38 + 182/5*a^36 - 4093/120*a^34 - 251/18*a^32 + 437/12*a^30 - 87/4*a^28 + 4873/360*a^26 - 67/4*a^24 + 91/12*a^22 + 47/6*a^20 - 637/60*a^18 + 137/30*a^16 - 7/10*a^14)*t^3 + (1/2520*a^98 - 1/60*a^86 - 1/20*a^78 + 3/10*a^76 - 7/72*a^74 + 3/2*a^70 - 7/3*a^68 + 49/30*a^66 - 298/15*a^64 + 61/4*a^62 + 297/4*a^60 - 15409/120*a^58 + 5251/90*a^56 + 271/60*a^54 - 47/6*a^52 + 275/24*a^50 - 115/6*a^48 + 269/12*a^46 - 613/36*a^44 + 11/4*a^42 + 29/6*a^40 + 883/180*a^38 - 599/30*a^36 + 63/4*a^34 + 5*a^32 - 85/6*a^30 + 17/2*a^28 - 85/12*a^26 + 19/2*a^24 - 4*a^22 - 4*a^20 + 5*a^18 - 2*a^16 + 2/7*a^14)*t^2 + (-1/2520*a^98 + 1/120*a^86 + 1/40*a^78 - 1/10*a^76 + 1/24*a^74 - 1/2*a^70 + 1/2*a^68 - 1/2*a^66 + 5*a^64 - 5/2*a^62 - 15*a^60 + 21*a^58 - 8*a^56)*t)*w^7 + ((-6*a^54 + 10*a^52 - 3*a^50 - 2*a^48 + 1/2*a^46 + 1/4*a^42 - 1/60*a^34)*t^7 + (-10*a^56 + 32*a^54 - 16*a^52 - 17*a^50 + 27/2*a^48 - 2*a^46 + 3*a^44 - 7/4*a^42 - 1/6*a^40 + 1/6*a^38 - 1/3*a^36 + 11/60*a^34)*t^6 + (-2*a^60 - 5/2*a^58 + 34*a^56 - 89/2*a^54 - 14*a^52 + 125/4*a^50 - 7/2*a^48 + 17/3*a^46 - 12*a^44 + 3*a^42 + 7/6*a^40 - 13/6*a^38 + 7/3*a^36 - 3/4*a^34)*t^5 + (-1/4*a^66 - 1/2*a^62 + 5*a^60 + 103/12*a^58 - 81/2*a^56 + 17*a^54 + 82/3*a^52 - 23/6*a^50 + 2/3*a^48 - 68/3*a^46 + 16*a^44 - 5/4*a^42 - 17/6*a^40 + 41/6*a^38 - 17/3*a^36 + 17/12*a^34)*t^4 + (-1/60*a^74 + 5/12*a^66 + 1/6*a^64 + 5/6*a^62 - 4*a^60 - 17/2*a^58 + 113/6*a^56 + 17/6*a^54 - 20/3*a^52 - 43/12*a^50 - 113/6*a^48 + 143/6*a^46 - 8*a^44 - 3/4*a^42 + 17/6*a^40 - 47/6*a^38 + 17/3*a^36 - 37/30*a^34)*t^3 + (1/60*a^74 - 1/6*a^66 - 3/20*a^64 - 1/3*a^62 + a^60 + 29/12*a^58 - 31/12*a^56 - 4/3*a^54 - 7/6*a^52 - 11/6*a^50 + 79/6*a^48 - 46/3*a^46 + 7*a^44 + 1/2*a^42 - a^40 + 3*a^38 - 2*a^36 + 2/5*a^34)*t^2 + (-1/60*a^64 + 1/4*a^56 + 1/2*a^52 - 2*a^50 - 3*a^48 + 10*a^46 - 6*a^44)*t)*w^6 + ((2*a^48 - 3/2*a^46 + 1/6*a^42)*t^6 + (-2*a^52 + 7*a^50 - 12*a^48 + 16/3*a^46 - 5/3*a^44)*t^5 + (-4/3*a^54 + 8*a^52 - 18*a^50 + 21*a^48 - 65/6*a^46 + 23/3*a^44 - 17/6*a^42 + 3*a^40)*t^4 + (-1/6*a^58 + 1/3*a^56 + 2*a^54 - 9*a^52 + 31/2*a^50 - 15*a^48 + 35/3*a^46 - 43/3*a^44 + 71/6*a^42 - 9*a^40 + 3/2*a^38 - 2*a^36)*t^3 + (1/6*a^58 - 1/3*a^56 - 2/3*a^54 + 3*a^52 - 9/2*a^50 + 4*a^48 - 9/2*a^46 + 8*a^44 - 31/3*a^42 + 11*a^40 - 17/2*a^38 + 6*a^36)*t^2 + (-1/6*a^46 + 1/3*a^44 + 7/6*a^42 - 5*a^40 + 7*a^38 - 4*a^36)*t)*w^5 + (-1/3*a^38*t^5 + (a^42 - 2*a^40 + 2*a^38 + 1/3*a^32)*t^4 + (-1/3*a^50 + 2*a^46 - 2*a^44 - 3*a^42 + 6*a^40 - 11/3*a^38 - a^36 + 2*a^34 - 2*a^32)*t^3 + (1/3*a^50 - 2*a^46 + 7/3*a^44 + 2*a^42 - 6*a^40 + 4*a^38 + 3*a^36 - 6*a^34 + 11/3*a^32)*t^2 + (-1/3*a^44 + 2*a^40 - 2*a^38 - 2*a^36 + 4*a^34 - 2*a^32)*t)*w^4]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "# we take the parameter a = sqrt{2}\n",
    "QDTzero = [qdt[i] * a^(i*(i-1)) for i in srange(disp)]\n",
    "QDTzero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[t,\n",
       " 4*t^2*w - 4*t*w,\n",
       " 20*t^3*w^2 - 28*t^2*w^2 - 4*t^2*w + 8*t*w^2 + 4*t*w,\n",
       " 488/3*t^4*w^3 - 232*t^3*w^3 - 64*t^3*w^2 + 112*t^2*w^3 + 96*t^2*w^2 - 128/3*t*w^3 - 32*t*w^2,\n",
       " 7100/3*t^5*w^4 - 8680/3*t^4*w^4 - 1280*t^4*w^3 + 2140*t^3*w^4 + 2048*t^3*w^3 - 248*t^2*w^4 + 32*t^3*w^2 - 896*t^2*w^3 - 4096/3*t*w^4 - 96*t^2*w^2 + 128*t*w^3 + 64*t*w^2,\n",
       " 965768/15*t^6*w^5 - 56496*t^5*w^5 - 124928/3*t^5*w^4 + 229048/3*t^4*w^5 + 180224/3*t^4*w^4 + 56720*t^3*w^5 + 2048*t^4*w^3 - 30720*t^3*w^4 + 271808/3*t^2*w^5 - 6656*t^3*w^3 + 40960/3*t^2*w^4 - 3473408/15*t*w^5 + 5632*t^2*w^3 - 4096/3*t*w^4 - 1024*t*w^3,\n",
       " 30174104/9*t^7*w^6 - 7449928/5*t^6*w^6 - 7270400/3*t^6*w^5 + 5533112*t^5*w^6 + 8163328/3*t^5*w^5 + 87461224/9*t^4*w^6 + 163840*t^5*w^4 - 5894144/3*t^4*w^5 + 87938960/3*t^3*w^6 - 524288*t^4*w^4 + 913408/3*t^3*w^5 + 2895728576/45*t^2*w^6 - 2048/3*t^4*w^3 + 430080*t^3*w^4 + 1449984*t^2*w^5 - 4984930304/45*t*w^6 + 4096*t^3*w^3 - 77824*t^2*w^4 - 262144/3*t*w^5 - 18432*t^2*w^3 + 8192*t*w^4 + 45056/3*t*w^3,\n",
       " 106687687696/315*t^8*w^7 + 681737168/45*t^7*w^7 - 3955785728/15*t^7*w^6 + 12074075696/15*t^6*w^7 + 924614656/5*t^6*w^6 + 18568018096/9*t^5*w^7 + 63963136/3*t^6*w^5 - 281116672*t^5*w^6 + 70331744128/9*t^4*w^7 - 62914560*t^5*w^5 - 731348992/3*t^4*w^6 + 1476072351296/45*t^3*w^7 - 524288/3*t^5*w^4 + 134217728/3*t^4*w^5 - 347275264*t^3*w^6 + 5311318221568/45*t^2*w^7 + 3211264/3*t^4*w^4 - 41156608/3*t^3*w^5 + 14709161984/15*t^2*w^6 - 50988241125376/315*t*w^7 - 4849664*t^3*w^4 + 10747904*t^2*w^5 - 444596224/15*t*w^6 + 13303808/3*t^2*w^4 - 524288/3*t*w^5 - 1441792/3*t*w^4]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs,\n",
    "# with the marking variable t for strongly connected components\n",
    "QDTzerosub = [sum(sum(QDTzero[i][k][j].subs(a=sqrt(2)) * t^j * w^k for j in srange(N)) for k in srange(N)) for i in srange(disp)]\n",
    "QDTzerosub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[t, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 4*t^2 - 4*t, 0, 0, 0, 0, 0, 0],\n",
       " [0, -4*t^2 + 4*t, 20*t^3 - 28*t^2 + 8*t, 0, 0, 0, 0, 0],\n",
       " [0,\n",
       "  0,\n",
       "  -64*t^3 + 96*t^2 - 32*t,\n",
       "  488/3*t^4 - 232*t^3 + 112*t^2 - 128/3*t,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  32*t^3 - 96*t^2 + 64*t,\n",
       "  -1280*t^4 + 2048*t^3 - 896*t^2 + 128*t,\n",
       "  7100/3*t^5 - 8680/3*t^4 + 2140*t^3 - 248*t^2 - 4096/3*t,\n",
       "  0,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  2048*t^4 - 6656*t^3 + 5632*t^2 - 1024*t,\n",
       "  -124928/3*t^5 + 180224/3*t^4 - 30720*t^3 + 40960/3*t^2 - 4096/3*t,\n",
       "  965768/15*t^6 - 56496*t^5 + 229048/3*t^4 + 56720*t^3 + 271808/3*t^2 - 3473408/15*t,\n",
       "  0,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  -2048/3*t^4 + 4096*t^3 - 18432*t^2 + 45056/3*t,\n",
       "  163840*t^5 - 524288*t^4 + 430080*t^3 - 77824*t^2 + 8192*t,\n",
       "  -7270400/3*t^6 + 8163328/3*t^5 - 5894144/3*t^4 + 913408/3*t^3 + 1449984*t^2 - 262144/3*t,\n",
       "  30174104/9*t^7 - 7449928/5*t^6 + 5533112*t^5 + 87461224/9*t^4 + 87938960/3*t^3 + 2895728576/45*t^2 - 4984930304/45*t,\n",
       "  0],\n",
       " [0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  -524288/3*t^5 + 3211264/3*t^4 - 4849664*t^3 + 13303808/3*t^2 - 1441792/3*t,\n",
       "  63963136/3*t^6 - 62914560*t^5 + 134217728/3*t^4 - 41156608/3*t^3 + 10747904*t^2 - 524288/3*t,\n",
       "  -3955785728/15*t^7 + 924614656/5*t^6 - 281116672*t^5 - 731348992/3*t^4 - 347275264*t^3 + 14709161984/15*t^2 - 444596224/15*t,\n",
       "  106687687696/315*t^8 + 681737168/45*t^7 + 12074075696/15*t^6 + 18568018096/9*t^5 + 70331744128/9*t^4 + 1476072351296/45*t^3 + 5311318221568/45*t^2 - 50988241125376/315*t]]"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs in matrix form, (Table 9)\n",
    "# the marking variable t is for strongly connected components\n",
    "QDTzeroMatrix = [[QDTzero[i][j].subs(a=sqrt(2)) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTzeroMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0]]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs in matrix form\n",
    "QDzeroMatrix = [[QDTzeroMatrix[i][j].subs(t=1) for j in srange(disp)] for i in srange(disp)]\n",
    "QDzeroMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, -4, 0, 0, 0, 0, 0, 0],\n",
       " [0, 4, 8, 0, 0, 0, 0, 0],\n",
       " [0, 0, -32, -128/3, 0, 0, 0, 0],\n",
       " [0, 0, 64, 128, -4096/3, 0, 0, 0],\n",
       " [0, 0, 0, -1024, -4096/3, -3473408/15, 0, 0],\n",
       " [0, 0, 0, 45056/3, 8192, -262144/3, -4984930304/45, 0],\n",
       " [0, 0, 0, 0, -1441792/3, -524288/3, -444596224/15, -50988241125376/315]]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with one strongly connected component (Table 4)\n",
    "QDTzeroMatrixPrim = [[QDTzeroMatrix[i][j].diff(t).subs(t=0) for j in srange(disp)] for i in srange(disp)]\n",
    "QDTzeroMatrixPrim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 4, 0, 0, 0, 0, 0, 0],\n",
       " [0, -4, -28, 0, 0, 0, 0, 0],\n",
       " [0, 0, 96, 112, 0, 0, 0, 0],\n",
       " [0, 0, -96, -896, -248, 0, 0, 0],\n",
       " [0, 0, 0, 5632, 40960/3, 271808/3, 0, 0],\n",
       " [0, 0, 0, -18432, -77824, 1449984, 2895728576/45, 0],\n",
       " [0, 0, 0, 0, 13303808/3, 10747904, 14709161984/15, 5311318221568/45]]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with two strongly connected components (Table 10)\n",
    "QDTzeroMatrixSec = [[QDTzeroMatrix[i][j].diff(t,2).subs(t=0) / 2 for j in srange(disp)] for i in srange(disp)]\n",
    "QDTzeroMatrixSec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 20, 0, 0, 0, 0, 0],\n",
       " [0, 0, -64, -232, 0, 0, 0, 0],\n",
       " [0, 0, 32, 2048, 2140, 0, 0, 0],\n",
       " [0, 0, 0, -6656, -30720, 56720, 0, 0],\n",
       " [0, 0, 0, 4096, 430080, 913408/3, 87938960/3, 0],\n",
       " [0, 0, 0, 0, -4849664, -41156608/3, -347275264, 1476072351296/45]]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients of the CGF of labeled digraphs with three strongly connected components (Table 11)\n",
    "QDTzeroMatrixThir = [[QDTzeroMatrix[i][j].diff(t,3).subs(t=0) / 6 for j in srange(disp)] for i in srange(disp)]\n",
    "QDTzeroMatrixThir"
   ]
  },
  {
   "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
}
