{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1a1128ed-285f-4a20-a3b0-8b48542e698a",
   "metadata": {},
   "source": [
    "Let $r\\geq 1$ be any positive integer, $R=\\mathbb{Z}[s]/(\\Phi_r(s))$ denote the $r$th cyclotomic integers and $Q$ its field of fractions. Then the Laurent polynomial $I_r\\in R[x^{\\pm 1},y^{\\pm 1},t]$, defined by the trace\n",
    "\\begin{equation}\n",
    "    I_r=\\operatorname{Tr}[A(s^{r-1})\\cdot A(s^{r-2})\\cdot\\ldots\\cdot A(s^2)\\cdot A(s)\\cdot A(1)]-(t^{r}+1),\n",
    "\\end{equation}\n",
    "where\n",
    "\\begin{equation}\n",
    " A(z)=A_0+z\\,A_1+z^2 A_2,\\qquad A_2=\\begin{bmatrix}\n",
    " 1 & 0\\\\\n",
    " 0 & 0\n",
    " \\end{bmatrix},\n",
    "\\end{equation}\n",
    "with\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "A_0&=\\begin{bmatrix}\n",
    "    t+x-xy & -w\\,x\\\\\n",
    "    w^{-1}(t+x-ty-2xy+xy^2) & x(y-1)\n",
    "\\end{bmatrix},\\\\   \n",
    "A_1&=\\begin{bmatrix}\n",
    "    y-x+x/y-1-t/x & w\\\\\n",
    "    w^{-1}(y-2x-1+xy+x/y-t/x) & 1\n",
    "\\end{bmatrix},\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "defines an integral of motion of $qP_I$.\n",
    "\n",
    "We implement this integral as an element of the field of rational functions $Q(t)(x, y)$ in $(x,y)$ over the field of rational functions $Q(t)$. The following function returns the integral of motion for any input value $r$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9a73c347-cfd6-4745-b42b-8e0368d101fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "//Input: a positive integer r.\n",
    "//Output: integral of motion with s an rth root of unity\n",
    "integral := function(r)\n",
    "    Q<s> := CyclotomicField(r); //Construct cyclotomic field of order r\n",
    "    Rt<t> := PolynomialRing(Q,1); \n",
    "    Qt<t> := FieldOfFractions(Rt); //Construct field Qt of rational functions in t over cyclotomic field\n",
    "    K<x,y> := PolynomialRing(Qt, 2); //Construct polynomial ring K in x,y over Qt\n",
    "    gl2 := MatrixRing(K, 2); //Construct ring of 2x2 matrices over K\n",
    "    //Define matrix coefficients of polynomial A(z) multiplied by x*y\n",
    "    A2 := Matrix(K, 2, 2, [x*y, 0, 0, 0]);\n",
    "    A1 := Matrix(K, 2, 2, [x^2 - t*y - x*y - x^2*y + x*y^2, x*y, x^2 - t*y - x*y - 2*x^2*y + x*y^2 + x^2*y^2, x*y]);\n",
    "    A0 := Matrix(K, 2, 2, [t*x*y + x^2*y - x^2*y^2, -x^2*y, t*x*y + x^2*y - t*x*y^2 - 2*x^2*y^2 + x^2*y^3,-x^2*y+x^2*y^2]);\n",
    "    linprobpol<z> := PolynomialRing(gl2); //Construct corresponding polynomial ring for A(z)\n",
    "    A := A0+A1*z+A2*z^2;  //Define A(z)\n",
    "    //Construct product of matrices\n",
    "    Aprod := Evaluate(A,1);\n",
    "    for i in [1..r-1] do\n",
    "        Aprod := Evaluate(A,Q!s^i)*Aprod;\n",
    "    end for;\n",
    "    //Take trace to obtain numerator of integral\n",
    "    integralnum := Trace(Aprod) - x^r*y^r*(1+t^r);\n",
    "    //Denominator is given by\n",
    "    integralden := x^r*y^r;\n",
    "    //ratio\n",
    "    Kfrac<x,y> := FieldOfFractions(K);\n",
    "    integralrat := (Kfrac!integralnum)/(Kfrac!integralden);\n",
    "    return integralrat;\n",
    "end function;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61c3e583-d3a1-4080-a75f-431944bbeb4b",
   "metadata": {},
   "source": [
    "Let's compute the first six integrals of motion $I_r$, $1\\leq r \\leq 6$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "70fce7bc-0377-4fdb-ab4f-7acf2b7a0a56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral for r=1:\n",
      "(-1*x^2*y + x^2 + x*y^2 + -t*y)/(x*y)\n",
      "Integral for r=2:\n",
      "(-1*x^4*y^2 + 2*x^4*y + -1*x^4 + -2*x^3*y^3 + 2*x^3*y^2 + -1*x^2*y^4 + 2*t*x^2*y + 2*t*x*y^3 + -t^2*y^2)/(x^2*y^2)\n",
      "Integral for r=3:\n",
      "(-1*x^6*y^3 + 3*x^6*y^2 + -3*x^6*y + x^6 + (-3*s - 3)*x^5*y^4 + (6*s + 6)*x^5*y^3 + (-3*s - 3)*x^5*y^2 + -3*s*x^4*y^5 + 3*s*x^4*y^4 + 3*t*x^4*y^2 + -3*t*x^4*y + x^3*y^6 + 3*s*t*x^3*y^4 + -3*t*x^2*y^5 + 3*t^2*x^2*y^2 + 3*t^2*x*y^4 + -t^3*y^3)/(x^3*y^3)\n",
      "Integral for r=4:\n",
      "(-1*x^8*y^4 + 4*x^8*y^3 + -6*x^8*y^2 + 4*x^8*y + -1*x^8 + -4*s*x^7*y^5 + 12*s*x^7*y^4 + -12*s*x^7*y^3 + 4*s*x^7*y^2 + 6*x^6*y^6 + -12*x^6*y^5 + 6*x^6*y^4 + 4*t*x^6*y^3 + -8*t*x^6*y^2 + 4*t*x^6*y + 4*s*x^5*y^7 + -4*s*x^5*y^6 + -4*t*x^5*y^5 + (4*s + 4)*t*x^5*y^4 + -4*s*t*x^5*y^3 + -1*x^4*y^8 + -8*s*t*x^4*y^6 + 4*s*t*x^4*y^5 + 4*t^2*x^4*y^3 + -6*t^2*x^4*y^2 + 4*t*x^3*y^7 + 4*s*t^2*x^3*y^5 + -4*t^2*x^3*y^4 + -6*t^2*x^2*y^6 + 4*t^3*x^2*y^3 + 4*t^3*x*y^5 + -t^4*y^4)/(x^4*y^4)\n",
      "Integral for r=5:\n",
      "(-1*x^10*y^5 + 5*x^10*y^4 + -10*x^10*y^3 + 10*x^10*y^2 + -5*x^10*y + x^10 + (-5*s^3 - 5*s^2 - 5*s - 5)*x^9*y^6 + (20*s^3 + 20*s^2 + 20*s + 20)*x^9*y^5 + (-30*s^3 - 30*s^2 - 30*s - 30)*x^9*y^4 + (20*s^3 + 20*s^2 + 20*s + 20)*x^9*y^3 + (-5*s^3 - 5*s^2 - 5*s - 5)*x^9*y^2 + -10*s^3*x^8*y^7 + 30*s^3*x^8*y^6 + -30*s^3*x^8*y^5 + (5*t + 10*s^3)*x^8*y^4 + -15*t*x^8*y^3 + 15*t*x^8*y^2 + -5*t*x^8*y + 10*s^2*x^7*y^8 + -20*s^2*x^7*y^7 + (5*s^3*t + 10*s^2)*x^7*y^6 + (10*s^2 + 10*s + 10)*t*x^7*y^5 + (-15*s^3 - 20*s^2 - 20*s - 20)*t*x^7*y^4 + (10*s^3 + 10*s^2 + 10*s + 10)*t*x^7*y^3 + -5*s*x^6*y^9 + 5*s*x^6*y^8 + -15*s^2*t*x^6*y^7 + (5*s^3 + 20*s^2)*t*x^6*y^6 + (-5*s^3 - 5*s^2)*t*x^6*y^5 + 5*t^2*x^6*y^4 + -15*t^2*x^6*y^3 + 10*t^2*x^6*y^2 + x^5*y^10 + 15*s*t*x^5*y^8 + -10*s*t*x^5*y^7 + 5*s^2*t^2*x^5*y^6 + (-5*s^3 - 5*s^2 - 5*s)*t^2*x^5*y^4 + -5*t*x^4*y^9 + -15*s*t^2*x^4*y^7 + (5*s + 5)*t^2*x^4*y^6 + 5*t^3*x^4*y^4 + -10*t^3*x^4*y^3 + 10*t^2*x^3*y^8 + 5*s*t^3*x^3*y^6 + -10*t^3*x^3*y^5 + -10*t^3*x^2*y^7 + 5*t^4*x^2*y^4 + 5*t^4*x*y^6 + -t^5*y^5)/(x^5*y^5)\n",
      "Integral for r=6:\n",
      "(-1*x^12*y^6 + 6*x^12*y^5 + -15*x^12*y^4 + 20*x^12*y^3 + -15*x^12*y^2 + 6*x^12*y + -1*x^12 + (-6*s + 6)*x^11*y^7 + (30*s - 30)*x^11*y^6 + (-60*s + 60)*x^11*y^5 + (60*s - 60)*x^11*y^4 + (-30*s + 30)*x^11*y^3 + (6*s - 6)*x^11*y^2 + 15*s*x^10*y^8 + -60*s*x^10*y^7 + 90*s*x^10*y^6 + (6*t - 60*s)*x^10*y^5 + (-24*t + 15*s)*x^10*y^4 + 36*t*x^10*y^3 + -24*t*x^10*y^2 + 6*t*x^10*y + -20*x^9*y^9 + 60*x^9*y^8 + (-6*s*t - 60)*x^9*y^7 + ((36*s - 18)*t + 20)*x^9*y^6 + (-72*s + 54)*t*x^9*y^5 + (60*s - 54)*t*x^9*y^4 + (-18*s + 18)*t*x^9*y^3 + (-15*s + 15)*x^8*y^10 + (30*s - 30)*x^8*y^9 + (24*t - 15*s + 15)*x^8*y^8 + (-18*s - 54)*t*x^8*y^7 + (36*s + 36)*t*x^8*y^6 + (6*t^2 + (-18*s - 6)*t)*x^8*y^5 + -27*t^2*x^8*y^4 + 36*t^2*x^8*y^3 + -15*t^2*x^8*y^2 + 6*s*x^7*y^11 + -6*s*x^7*y^10 + (36*s - 36)*t*x^7*y^9 + (-54*s + 60)*t*x^7*y^8 + (-6*t^2 + (18*s - 24)*t)*x^7*y^7 + 12*s*t^2*x^7*y^6 + (-30*s + 30)*t^2*x^7*y^5 + (18*s - 24)*t^2*x^7*y^4 + -1*x^6*y^12 + -24*s*t*x^6*y^10 + 18*s*t*x^6*y^9 + (-27*s + 27)*t^2*x^6*y^8 + (30*s - 30)*t^2*x^6*y^7 + 6*t^3*x^6*y^5 + -24*t^3*x^6*y^4 + 20*t^3*x^6*y^3 + 6*t*x^5*y^11 + 36*s*t^2*x^5*y^9 + (-18*s - 6)*t^2*x^5*y^8 + (6*s - 6)*t^3*x^5*y^7 + (-6*s - 6)*t^3*x^5*y^6 + (-6*s + 24)*t^3*x^5*y^5 + -15*t^2*x^4*y^10 + -24*s*t^3*x^4*y^8 + (6*s + 18)*t^3*x^4*y^7 + 6*t^4*x^4*y^5 + -15*t^4*x^4*y^4 + 20*t^3*x^3*y^9 + 6*s*t^4*x^3*y^7 + -18*t^4*x^3*y^6 + -15*t^4*x^2*y^8 + 6*t^5*x^2*y^5 + 6*t^5*x*y^7 + -t^6*y^6)/(x^6*y^6)\n"
     ]
    }
   ],
   "source": [
    "    for r in [1..6] do\n",
    "        print(\"Integral for r=\" cat Sprint(r) cat \":\");\n",
    "        integral(r);\n",
    "    end for;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13659fa4-1bc9-4a0a-8791-72d2c7c4268f",
   "metadata": {},
   "source": [
    "We want to check the genus of generic fibres of the integrals. To do this, we consider\n",
    "\\begin{equation}\n",
    "x^r y^r I_r(x,y)-c x^r y^r=0\n",
    "\\end{equation}\n",
    "as defining a curve in $(x,y)\\in \\mathbb{P}^1\\times\\mathbb{P}^1$ over $Q(t,c)$, the field of rational functions in $(t,c)$ over $Q$. (In the code below we take the canonical projective completion of the curve in $\\mathbb{P}^2$ instead, but this is irrelevant when it comes to the geometric genus of the curve).\n",
    "\n",
    "The following function returns this genus for any input value of $r$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e3a64379-1e1a-427f-8a01-39c5d381d8dc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "//Input: a positive integer r.\n",
    "//Output: genus of generic fibres of integral of motion when s an rth root of unity.\n",
    "fibre_genus := function(r)\n",
    "    Q<s> := CyclotomicField(r); //Construct cyclotomic field of order r\n",
    "    Rct<t,c> := PolynomialRing(Q,2); \n",
    "    Qct<t,c> := FieldOfFractions(Rct); //Construct field Qct of rational functions in t,c over cyclotomic field\n",
    "    K<x,y> := PolynomialRing(Qct, 2); //Construct polynomial ring K in x,y over Qct\n",
    "    gl2 := MatrixRing(K, 2); //Construct ring of 2x2 matrices over K\n",
    "    //Define matrix coefficients of polynomial A(z) multiplied by x*y\n",
    "    A2 := Matrix(K, 2, 2, [x*y, 0, 0, 0]);\n",
    "    A1 := Matrix(K, 2, 2, [x^2 - t*y - x*y - x^2*y + x*y^2, x*y, x^2 - t*y - x*y - 2*x^2*y + x*y^2 + x^2*y^2, x*y]);\n",
    "    A0 := Matrix(K, 2, 2, [t*x*y + x^2*y - x^2*y^2, -x^2*y, t*x*y + x^2*y - t*x*y^2 - 2*x^2*y^2 + x^2*y^3,-x^2*y+x^2*y^2]);\n",
    "    linprobpol<z> := PolynomialRing(gl2); //Construct corresponding polynomial ring for A(z)\n",
    "    A := A0+A1*z+A2*z^2;  //Define A(z)\n",
    "    //Construct product of matrices\n",
    "    Aprod := Evaluate(A,1);\n",
    "    for i in [1..r-1] do\n",
    "        Aprod := Evaluate(A,Q!s^i)*Aprod;\n",
    "    end for;\n",
    "    //Take trace to obtain invariant pencil\n",
    "    pencilinhom := Trace(Aprod) - x^r*y^r*(1+t^r) - c*x^r*y^r;\n",
    "    //Change to homogeneous coordinates\n",
    "    POLXYZ<X,Y,Z> := PolynomialRing(Qct, 3);\n",
    "    phi := hom< K -> POLXYZ | [X, Y] >;\n",
    "    pencilinhomXYZ := phi(pencilinhom); //Coerce inhomogeneous pencil\n",
    "    pencilhom := Homogenization(pencilinhomXYZ,Z); //compute homogeneous form of pencil\n",
    "    //Construct projective curve over Qct\n",
    "    P2<X, Y, Z> := ProjectiveSpace(Qct, 2); \n",
    "    C := Curve(P2, pencilhom);  \n",
    "    g := Genus(C); //Compute genus\n",
    "    return g;\n",
    "end function;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c74e5366-0a58-4022-bb7d-2e00ca54fa18",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "1\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "fibre_genus(1);\n",
    "fibre_genus(2);\n",
    "fibre_genus(3);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18587457-ea11-43bb-8156-94ec59904f72",
   "metadata": {},
   "source": [
    "Since fibre_genus(r) takes forever to compute the genus when $r\\geq 4$, we implement an analogous function that computes the genus for specific values of $c,t$. To be precise, it compute the genus for $c$ ranging from $0$ to $10$ and $t$ ranging from $1$ to $10$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b2b054c8-2eaf-41d4-83b1-3dd5a22eaa3f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "//Input: a positive integer r.\n",
    "//Ouput: list with genera of fibres of integral of motion for c ranging from 0 to 10 and t from 1 to 10.\n",
    "fibre_genus_num := function(r)\n",
    "    Q<s> := CyclotomicField(r); //Construct cyclotomic field of order r\n",
    "    K<x,y,t,c> := PolynomialRing(Q, 4); //Construct polynomial ring K in x,y,t,c over Q\n",
    "    gl2 := MatrixRing(K, 2); //Construct ring of 2x2 matrices over K\n",
    "    //Define matrix coefficients of polynomial A(z) multiplied by x*y\n",
    "    A2 := Matrix(K, 2, 2, [x*y, 0, 0, 0]);\n",
    "    A1 := Matrix(K, 2, 2, [x^2 - t*y - x*y - x^2*y + x*y^2, x*y, x^2 - t*y - x*y - 2*x^2*y + x*y^2 + x^2*y^2, x*y]);\n",
    "    A0 := Matrix(K, 2, 2, [t*x*y + x^2*y - x^2*y^2, -x^2*y, t*x*y + x^2*y - t*x*y^2 - 2*x^2*y^2 + x^2*y^3,-x^2*y+x^2*y^2]);\n",
    "    linprobpol<z> := PolynomialRing(gl2); //Construct corresponding polynomial ring for A(z)\n",
    "    A := A0+A1*z+A2*z^2;  //Define A(z)\n",
    "    //Construct product of matrices\n",
    "    Aprod := Evaluate(A,1);\n",
    "    for i in [1..r-1] do\n",
    "        Aprod := Evaluate(A,Q!s^i)*Aprod;\n",
    "    end for;\n",
    "    //Take trace to obtain invariant pencil\n",
    "    pencilinhom := Trace(Aprod) - x^r*y^r*(1+t^r) - c*x^r*y^r;\n",
    "    //Change to homogeneous coordinates\n",
    "    POLXYZ<X,Y,Z> := PolynomialRing(Q, 3);\n",
    "    P2<X, Y, Z> := ProjectiveSpace(Q, 2); \n",
    "    genusdata := [];\n",
    "    for cval in [0..10] do\n",
    "        for tval in [1..10] do\n",
    "            poleval := Evaluate(pencilinhom, [X, Y, tval, cval]);  \n",
    "            polXYZ := POLXYZ!poleval;\n",
    "            polXYZhom := Homogenization(polXYZ,POLXYZ!Z); //compute homogeneous form of pencil\n",
    "            //Construct projective curve over Q\n",
    "            C := Curve(P2, polXYZhom);  \n",
    "            g := Genus(C);\n",
    "            Append(~genusdata, g);\n",
    "        end for;\n",
    "    end for;\n",
    "    //Compute genus\n",
    "    return genusdata;\n",
    "end function;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "353fbf71-2384-45ed-889d-265e864a3d08",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]\n",
      "[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]\n",
      "[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]\n",
      "[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]\n"
     ]
    }
   ],
   "source": [
    "fibre_genus_num(1);\n",
    "fibre_genus_num(2);\n",
    "fibre_genus_num(3);\n",
    "fibre_genus_num(4);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Magma",
   "language": "magma",
   "name": "magma"
  },
  "language_info": {
   "codemirror_mode": "pascal",
   "file_extension": ".m",
   "mimetype": "text/x-pascal",
   "name": "magma",
   "version": "2.28-5\r"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}