{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8565ca20-c37d-4b6c-8aa0-a32074c8e6e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# g = Graph(17)\n",
    "# a = graph_editor(g)\n",
    "# a.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "118f3510-30b1-4a0c-81e3-146325f5c8d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A =g.adjacency_matrix()\n",
    "# A[0, 0] = -13\n",
    "# for i in (1, 3, 5, 6, 7):\n",
    "#     A[i, i] = -2\n",
    "#     A[i + 8, i + 8] = -2\n",
    "# for i in (2,):\n",
    "#     A[i, i] = -4\n",
    "#     A[i + 8, i + 8] = -4\n",
    "# for i in (4, 8):\n",
    "#     A[i, i] = -1\n",
    "#     A[i + 8, i + 8] = -1\n",
    "# # A = A.change_ring(QQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1c1485e6-9433-4b5b-93c4-20dc6d95926c",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = Matrix(QQ, 17, [-13, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -4, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -4, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -2, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1])\n",
    "Vf = vector(4 * [0] + [1] + 7 * [0] + [1] + 4 * [0])\n",
    "Vg = vector(8 * [0] + [1] + 7 * [0] + [1])\n",
    "Vd = vector([13] + 7 * [0] + [1] + 7 * [0] + [1])\n",
    "R.<x, y, z> = QQ[]\n",
    "S.<s> = QQ[]\n",
    "T.<t> = QQ[]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2d89bcf0-311e-4799-872c-72b29b23ee83",
   "metadata": {},
   "outputs": [],
   "source": [
    "def trinomiales(A, Vf, Vg, Vd, total=True):\n",
    "    n = A.nrows()\n",
    "    A1 = -A^-1\n",
    "    Nf, Ng, Nd = [A1 * N0 for N0 in (Vf, Vg, Vd)]\n",
    "    diag = vector(A.diagonal())\n",
    "    unos = vector(n * [1])\n",
    "    if total:\n",
    "        VV = Vf + Vg\n",
    "    else:\n",
    "        VV = Vf.parent()(0)\n",
    "    euler = (-A + 2 * identity_matrix(n)) * unos + diag - VV - Vd\n",
    "    nus = A1 * (diag + 2 * unos) + unos\n",
    "    expos = {}\n",
    "    for j in range(n):\n",
    "        a, b = (Nf[j], Ng[j])\n",
    "        eq = 3 * Nf[j] == 2 * Ng[j]\n",
    "        nnd = min(3 * Nf[j], 2 * Ng[j]) == Nd[j]\n",
    "        if nnd:\n",
    "            expos[j] = ((a, b), (a, b), nus[j], eq, nnd)\n",
    "        else:\n",
    "            expos[j] = ((a, b), (Nf[j], Nd[j] - Ng[j]), nus[j], eq, nnd)\n",
    "    cadenas = ('d',)\n",
    "    if total:\n",
    "        expos['f'] = ((1, 0), (0, 0), 1, False, True)\n",
    "        expos['g'] = ((0, 1), (0, 0), 1, False, True)\n",
    "        cadenas = ('f', 'g') + cadenas\n",
    "    expos['d'] = ((0, 0), (0, 1), 1, True, False)\n",
    "    pols = {}\n",
    "    for j in range(n):\n",
    "        if expos[j][-1]:\n",
    "            a, b = expos[j][0]\n",
    "            pols[j] = (z^3 + x^a * z + x^b, (expos[j][2], 1, 1), euler[j])\n",
    "        else:\n",
    "            a, b = expos[j][1]\n",
    "            pols[j] = (z^3 + x^a * z^2 + x^b, (expos[j][2], 1, 1), euler[j])\n",
    "    interdiv = [(i, j) for i, j in Combinations(n, 2) if A[i, j] != 0]\n",
    "    interstrict = [(i, j) for i in range(n) for j in cadenas if eval('V' + j)[i] != 0]\n",
    "    for i,j in interdiv + interstrict:\n",
    "        if j in ZZ:\n",
    "            ep = 1\n",
    "        else:\n",
    "            ep = eval('V' + j)[i]\n",
    "        eqi, nndi = expos[i][-2:]\n",
    "        eqj, nndj = expos[j][-2:]\n",
    "        nd = eqi and eqj and not (nndi and nndj)\n",
    "        if not nd:\n",
    "            ai, bi = expos[i][0]\n",
    "            aj, bj = expos[j][0]\n",
    "            pols[i, j] = (z^3 + x^ai * y^aj * z + x^bi * y^bj,\n",
    "                          (expos[i][2], expos[j][2], 1), ep)\n",
    "        else:\n",
    "            ai, bi = expos[i][1]\n",
    "            aj, bj = expos[j][1]\n",
    "            pols[i, j] = (z^3 + x^ai * y^aj * z^2 + x^bi * y^bj, \n",
    "                          (expos[i][2], expos[j][2], 1), ep)\n",
    "    return n, Nf, Ng, Nd, euler, nus, expos, pols\n",
    "\n",
    "def estrato(pol, pesos, euler):\n",
    "    Z1 = ZetaFunctions(pol)\n",
    "    zetatop = euler * Z1.topological_zeta(local=True, weights=pesos)\n",
    "    zetamon = Z1.monodromy_zeta()^euler\n",
    "    return zetatop, zetamon\n",
    "\n",
    "def zetas(n, pols):\n",
    "    zetatop=0\n",
    "    zetamon=1\n",
    "    for k in pols:\n",
    "        zt, zm = estrato(*pols[k])\n",
    "        zetatop += zt\n",
    "        zetamon *= zm\n",
    "    return zetatop, zetamon / (1 - t^3)^n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "779f6c34-62bc-4484-8caa-f5ef789693a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "n, Nf, Ng, Nd, euler, nus, expos, pols = trinomiales(A, Vf, Vg, Vd)\n",
    "zt0, zm0 = zetas(n, pols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5d85cf37-2286-4a0a-a13a-3e5cd781d450",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-7/39) * (s + 19/39)^-1 * (s + 1/2)^-1 * (s + 1)^-1 * (s^2 - 1207/546*s - 19/14)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zt0.factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ae05ef5a-5687-43a8-bac6-7cc9361149d3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 4\n",
      "1 5\n",
      "6 4\n",
      "3 16\n",
      "26 4\n",
      "13 4\n",
      "78 4\n",
      "39 15\n"
     ]
    }
   ],
   "source": [
    "for p, j in zm0.factor():\n",
    "    print(p.is_cyclotomic(certificate=True), j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "016e3e91-9506-4223-9010-1b72c474175d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, [76/65/(s + 19/39), -2/39/(s + 1/2), -253/195/(s + 1)])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zt0.partial_fraction_decomposition()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "cb10bd5c-a177-4d82-8dec-ce05bd9491bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "n1, Nf1, Ng1, Nd1, euler1, nus1, expos1, pols1 = trinomiales(A, Vf, Vg, Vd, total=False)\n",
    "zt1, zm1 = zetas(n1, pols1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d03c12a4-58ef-414a-ba00-e5566f8b2201",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zt0 - zt1, zm0 / zm1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9ee974ba-af17-4998-85a9-d806a92a7242",
   "metadata": {},
   "outputs": [],
   "source": [
    "B = A.delete_columns([3, 4, 11, 12]).delete_rows([3, 4, 11, 12])\n",
    "B[2, 2] = -2\n",
    "B[8, 8] = -2\n",
    "Wf = vector(Vf.list()[:2] + [2] + Vf.list()[5:10] + [2] + Vf.list()[13:])\n",
    "Wg = vector(Vg.list()[:3] + Vg.list()[5:11] + Vg.list()[13:])\n",
    "Wd = vector(Vd.list()[:3] + Vd.list()[5:11] + Vd.list()[13:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e0e38fde-16ba-4e5b-826e-088da1bd1200",
   "metadata": {},
   "outputs": [],
   "source": [
    "n2, Nf2, Ng2, Nd2, euler2, nus2, expos2, pols2 = trinomiales(B, Wf, Wg, Wd, total=False)\n",
    "zt2, zm2 = zetas(n2, pols2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "de6c6fe8-c3b6-4e9e-b65f-10a41dc2ed05",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zt0 - zt2, zm1 / zm2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f09d5d7e-c788-4e6d-96b4-8cf999675ff8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.3.rc3",
   "language": "sage",
   "name": "sagemath-310"
  },
  "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}