Weierstrass_Points.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the labelling of these sections is based on the version of the paper completed on 15/08/22. This cell should be updated in the event of a change to the labelling in the paper.\n",
    "\n",
    "Currently running on Sagemath version 9.7.beta8, Release Date: 2022-08-07.\n",
    "Using Python 3.10.5.\n",
    "\n",
    "This notebook is intended to be ran sequentially from the beginning. Changing the order in which cells are ran may lead to errors. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remark 3.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x^12 - 32*x^11 - 114*x^10 - 200*x^9 + 100*x^8 + 48*x^7 - 936*x^6 + 1728*x^5 - 2000*x^4 + 3200*x^3 - 2624*x^2 + 768*x - 64 \n",
      "\n",
      "x^24 - 24*x^23 + 1306*x^22 - 2864*x^21 + 10096*x^20 - 32704*x^19 - 5704*x^18 - 41824*x^17 + 43056*x^16 + 831616*x^15 + 837856*x^14 + 992256*x^13 + 2603136*x^12 + 1238016*x^11 + 1560576*x^10 + 5584896*x^9 + 3357696*x^8 + 3838976*x^7 + 5856256*x^6 + 2543616*x^5 + 2200576*x^4 + 1355776*x^3 + 454656*x^2 + 65536*x + 4096 \n",
      "\n",
      "x^24 + 56*x^23 + 1176*x^22 - 1784*x^21 - 3904*x^20 + 36096*x^19 + 12776*x^18 - 211904*x^17 + 304736*x^16 + 431616*x^15 + 339456*x^14 - 1985664*x^13 - 625344*x^12 + 1034496*x^11 + 3512576*x^10 - 584704*x^9 - 3572224*x^8 - 2018304*x^7 + 3303936*x^6 + 3055616*x^5 + 1099776*x^4 + 45056*x^3 + 229376*x^2 - 16384*x + 4096 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# We will find the polynomials in the HC model which give the base coordinate of the Weierstrass points.\n",
    "R.<x, y> = QQ[]\n",
    "f = x*(y^5+1) + (x*y)^2 - x^4*y - 2*y^3\n",
    "g = 4\n",
    "\n",
    "# We initialise some values which will be useful in the calculation of the Wronskian determinant. \n",
    "FF = R.fraction_field()\n",
    "dyf = FF(f.derivative(y))\n",
    "dxf = FF(f.derivative(x))\n",
    "vs = [FF(vi)/dyf for vi in Curve(f).riemann_surface().cohomology_basis()]\n",
    "rat = dxf/dyf\n",
    " \n",
    "# We build the matrix we need, iteratively taking the derivatives \n",
    "# Away from the branch locus we can take x to be a local coordinate. \n",
    "WM = [vs]\n",
    "for i in range(g-1):\n",
    "    WM.append([w.derivative(x)-rat*w.derivative(y) for w in WM[-1]])\n",
    "WM = Matrix(WM)\n",
    "\n",
    "# We do this because of our knowledge that we expect a factor of (dyf)^(g*(g+1)/2) in the denominator, and we \n",
    "# want to make the calculation quicker\n",
    "for i in range(g):\n",
    "    WM[i,:] *= dyf^(i+1)\n",
    "\n",
    "# We can now calculate the Wronskian determinant, and we will print it's factors that don't \n",
    "# come from factors of the discriminant. \n",
    "W = R(WM.det())\n",
    "\n",
    "W_factors = [fac[0] for fac in f.resultant(W, y).factor()]\n",
    "bl_factors = set([fac[0] for fac in f.resultant(f.derivative(y), y).factor()])\n",
    "\n",
    "W_factors = [wf for wf in W_factors if wf not in bl_factors]\n",
    "\n",
    "for fi in W_factors:\n",
    "    print(fi,\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C4 x S3\n",
      "C4 x S3\n",
      "C4 x S3\n"
     ]
    }
   ],
   "source": [
    "# We check the Galois group for all these polynomials. \n",
    "for fac in W_factors:\n",
    "    SFi = QQ['x'](fac).splitting_field(names='g')\n",
    "    G = SFi.galois_group()\n",
    "    print(G.structure_description())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "# We check that the claimed field extension does split these factors. \n",
    "K.<z> = CyclotomicField(5)\n",
    "s = polygen(K)\n",
    "L.<ir2> = K.extension(s^2+2)\n",
    "s = polygen(L)\n",
    "M.<cr> = L.extension(s^3+7*s^2+8*s+4)\n",
    "\n",
    "print(all([len(p.change_ring(M).factor())==p.degree() for p in W_factors]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proposition 3.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "verbose 0 (3893: multi_polynomial_ideal.py, groebner_basis) Warning: falling back to very slow toy implementation.\n",
      "verbose 0 (1082: multi_polynomial_ideal.py, dimension) Warning: falling back to very slow toy implementation.\n",
      "Defining polynomial for delta: y1^2 + ((-11/14*a - 6/7)*b - 6/7*a - 4/7)*y1 + (3565/6272*a + 2231/1568)*b + 2231/1568*a + 989/392 \n",
      "\n",
      "Intersection points of osculating plane at W_{345}\n",
      "[((-d + (11/14*a + 6/7)*b + 6/7*a + 4/7 : d : (-43/112*a^2 - 113/112*a - 23/28)*b - 13/28*a^2 - 27/28*a - 5/7 : (43/112*a^2 + 25/112*a - 1/28)*b + 13/28*a^2 + 3/28*a - 6/7 : 1), 1), ((d : -d + (11/14*a + 6/7)*b + 6/7*a + 4/7 : (-43/112*a^2 - 113/112*a - 23/28)*b - 13/28*a^2 - 27/28*a - 5/7 : (43/112*a^2 + 25/112*a - 1/28)*b + 13/28*a^2 + 3/28*a - 6/7 : 1), 1), ((-1/4*a*b : -1/4*a*b : -1/4*a^2*b : (1/4*a^2 + 1/2*a)*b - 1 : 1), 4)]\n"
     ]
    }
   ],
   "source": [
    "# We first verify the points P, P' using computer algebra\n",
    "\n",
    "# We build the extension of QQ which contains the Weierstrass points in the P4 model. \n",
    "s = polygen(QQ)\n",
    "Ka.<a> = QQ.extension(s^3+2*s^2+3*s+4)\n",
    "# We can then add a second root\n",
    "s = polygen(Ka)\n",
    "Kab.<b> = Ka.extension(s^2+(a+2)*s+a^2+2*a+3)\n",
    "# The final root is given by the condition -(a+b+c)==2\n",
    "c = -2-a-b\n",
    "\n",
    "# We now want to find the points where the osculating plane at a WP intersects the canonical embedding\n",
    "# we know there will be an intersection of order 4 at the WP, and then two more points. \n",
    "\n",
    "# First let's create the ideal and get a groebner basis.\n",
    "Q5 = PolynomialRing(Kab, 5, names='x')\n",
    "gens = Q5.gens()\n",
    "Q5.inject_variables(verbose=False)\n",
    "Hs = [sum([xi^k for xi in gens]) for k in range(1,6)]\n",
    "wp = [1,1,a,b,c]\n",
    "# The equation for the osculating plane comes from Edge 1978. \n",
    "OP = Hs[3](*wp)*sum([wp[i]^2*gens[i] for i in range(5)])+2*Hs[4](*wp)*sum([wp[i]*gens[i] for i in range(5)])\n",
    "II = Q5.ideal(Hs[:3] + [OP])\n",
    "bs = II.groebner_basis()\n",
    "\n",
    "# As this is an a subscheme of projective space really, to deal with this as an ideal we need it to\n",
    "# have dimension 1, this degree of freedom coming from scaling. \n",
    "if II.dimension()!=1:\n",
    "    raise ValueError(\"Dimension of the ideal is inccorect, mistake somewhere.\")\n",
    "    \n",
    "if not len(bs)==4:\n",
    "    raise ValueError(\"Incorrect number of generating polynomials for this ideal.\")\n",
    "    \n",
    "# Now we know we have 4 polynomials in the basis to deal with we can do this manually and exactly.\n",
    "# Note that becuase of the ordering implicitly taken x4 is the homogenising variable for each. \n",
    "# We create a polynomial ring over Kab on the affine patch where x4 != 0\n",
    "# Because we know the only points where we can have x4 == 0, this will not be a restriction. \n",
    "i = 4\n",
    "Qi = PolynomialRing(Kab, i, names='y')\n",
    "Qi.inject_variables(verbose=False)\n",
    "# substitute into the groebner basis these affine coordinates. \n",
    "sgb = [bi.subs({x0:y0, x1:y1, x2:y2, x3:y3, x4:1}) for bi in bs]\n",
    "# Create a projective space to store the results\n",
    "P4_Kab = ProjectiveSpace(4, names='x', R=Kab)\n",
    "\n",
    "\n",
    "# One can check that for y3 == y3s_and_mults_exact[0][0], there were no solutions \n",
    "# for y1 over Kab. For y3 == y3s_and_mults_exact[0][0] the solutions gave only the order-4 intersection\n",
    "# at the Weierstrass point. \n",
    "# Hence we enlarge our field. \n",
    "y3s_and_mults_exact = Kab['y3'](sgb[0](1,1,1,y3)).roots()\n",
    "y3i = y3s_and_mults_exact[0][0]\n",
    "pol = Kab['y1'](sgb[1](1,y1,1,y3i))\n",
    "# Look at the polynomial again\n",
    "print(\"Defining polynomial for delta:\", pol, \"\\n\")\n",
    "Kabd.<d> = Kab.extension(pol)\n",
    "\n",
    "# Now repeat the process, to get all the solutions. \n",
    "P4_Kabd = ProjectiveSpace(4, names='x', R=Kabd)\n",
    "sols_with_mults_exact = []\n",
    "for y3i, mi in y3s_and_mults_exact:\n",
    "    pol = Kabd['y1'](sgb[1](1,y1,1,y3i))\n",
    "    y1s_and_mults = pol.roots()\n",
    "    for y1j, mj in y1s_and_mults:\n",
    "        y0ij = Kabd['y0'](sgb[2](y0,y1j,1,y3i)).roots(multiplicities=False)[0]\n",
    "        y2ij = Kabd['y2'](sgb[3](1,y1j,y2,y3i)).roots(multiplicities=False)[0]\n",
    "        sols_with_mults_exact.append((P4_Kabd([y0ij, y1j, y2ij, y3i, 1]), mi*mj))\n",
    "print(\"Intersection points of osculating plane at W_{345}\")\n",
    "print(sols_with_mults_exact)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "# We next check using the Abel-Jacobi map that these points do give a canonical divisor\n",
    "\n",
    "# We initialise the Riemann surface\n",
    "load(\"riemann_surface.py\")\n",
    "R.<x, y> = QQ[]\n",
    "f = x*(y^5+1) + (x*y)^2 - x^4*y - 2*y^3\n",
    "integration_method = 'rigorous'\n",
    "prec = 200\n",
    "S = RiemannSurface(f, prec=prec, integration_method=integration_method)\n",
    "\n",
    "# Convert the points to the HC model in a float field. \n",
    "embd = Kabd.embeddings(S._CC)[0]\n",
    "z = exp(2*S._CC.pi()*S._CC(sqrt(-1))/5)\n",
    "Ainv = matrix(4, 4, [1/5*z^2 - 1/5*z,\n",
    " -1/5*z + 1/5,\n",
    " 1/5*z^3 - 1/5*z,\n",
    " -1/5*z^3 - 1/5*z^2 - 2/5*z - 1/5,\n",
    " 1/5*z - 1/5,\n",
    " -1/5*z^2 + 1/5*z,\n",
    " 1/5*z^3 + 1/5*z^2 + 2/5*z + 1/5,\n",
    " -1/5*z^3 + 1/5*z,\n",
    " -1/5*z^3 + 1/5*z,\n",
    " 1/5*z^3 + 1/5*z^2 + 2/5*z + 1/5,\n",
    " 1/5*z - 1/5,\n",
    " -1/5*z^2 + 1/5*z,\n",
    " -1/5*z^3 - 1/5*z^2 - 2/5*z - 1/5,\n",
    " 1/5*z^3 - 1/5*z,\n",
    " 1/5*z^2 - 1/5*z,\n",
    " -1/5*z + 1/5])\n",
    "\n",
    "points_x4 = [sol for sol, _ in sols_with_mults_exact]\n",
    "points_L3 = [Ainv*vector([embd(co) for co in sol[:4]]) \n",
    "              for sol, _ in sols_with_mults_exact]\n",
    "\n",
    "points_XYZ = [(L[1]^2-L[0]*L[2], L[0]*L[3], L[1]*L[3]) for L in points_L3]\n",
    "points_xy = [(X/Z, Y/Z) for X, Y, Z in points_XYZ]\n",
    "p1, p2, W = points_xy\n",
    "\n",
    "# Calculate the AJ map to the points\n",
    "AJ_bp_p1 = S._aj_based(p1)\n",
    "AJ_bp_p2 = S._aj_based(p2)\n",
    "AJ_bp_W = S._aj_based(W)\n",
    "\n",
    "# Now compute the canonical divisor (using a known theta characteristic)\n",
    "F = Curve(f).function_field()\n",
    "divx = F(x).divisor()\n",
    "d, c, b, a = divx.support()\n",
    "avoid = [a, b, c, d]\n",
    "a_avoid, _ = S.strong_approximation(1*a, avoid)\n",
    "a_list = S.divisor_to_divisor_list(a_avoid)\n",
    "AJ_bp_a = S.abel_jacobi(a_list)\n",
    "Delta = a+3*b-d\n",
    "Delta_avoid, _ = S.strong_approximation(Delta, avoid)\n",
    "Delta_list = S.divisor_to_divisor_list(Delta_avoid)\n",
    "AJ_bp_D = S.abel_jacobi(Delta_list)\n",
    "\n",
    "# Check that the difference between the images of the two divisors in the Jacobian\n",
    "# is small wrt to the lattice. \n",
    "print(S.reduce_over_period_lattice(4*AJ_bp_W+AJ_bp_p1+AJ_bp_p2-2*AJ_bp_D).norm() < 1e-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.7.beta8",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}