Even_Characteristics.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.4.rc0, Release Date: 2021-07-27.\n",
    "Using Python 3.9.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": [
    "## Theorem 5.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# odd spin structures: 120\n",
      "# even spin structures: 136 \n",
      "\n",
      "Partitions:\n",
      "120 = [20, 20, 60, 20]\n",
      "136 = [30, 10, 30, 30, 5, 5, 10, 10, 5, 1]\n",
      "\n",
      "Space of invariant spin structures: Affine space p + W where:\n",
      "  p = (0, 0, 0, 0, 1, 0, 0, 1)\n",
      "  W = Vector space of degree 8 and dimension 0 over Finite Field of size 2\n",
      "Basis matrix:\n",
      "[]\n"
     ]
    }
   ],
   "source": [
    "# We will first use the method of Kallel and Sjerve\n",
    "\n",
    "from sage.geometry.hyperplane_arrangement.affine_subspace import AffineSubspace\n",
    "from sage.groups.perm_gps.partn_ref.refinement_graphs import get_orbits\n",
    "\n",
    "# Initialise the Curve\n",
    "from sage.schemes.riemann_surfaces.riemann_surface import RiemannSurface\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 = 60\n",
    "S = RiemannSurface(f, prec=prec, integration_method=integration_method)\n",
    "Mlist = S.symplectic_isomorphisms()\n",
    "g = S.genus\n",
    "\n",
    "# Initialise some useful objects\n",
    "I2g = matrix.identity(2*g)\n",
    "Ig = matrix.identity(g)\n",
    "Jg = block_matrix([[matrix.zero(g), Ig], [-Ig, matrix.zero(g)]])\n",
    "F2 = FiniteField(2)\n",
    "Vbars = [vector([sum([A[j,i]*A[k,i]*Jg[j,k] for j in range(2*g) for k in range(j, 2*g)]) \n",
    "                 for i in range(2*g)], F2) for A in Mlist] \n",
    "Mbars = [matrix(m, ring=F2) for m in Mlist]\n",
    "# XS will be the vectors corresponding to theta characteristics\n",
    "XS = [vector(Integer(n).digits(2, padto=2*g), F2) for n in range(4^g)]\n",
    "# split by parity (as to include odd characteristics here)\n",
    "odds = [X for X in XS if sum([X[i]*X[i+g] for i in range(g)])]\n",
    "evens = [X for X in XS if not sum([X[i]*X[i+g] for i in range(g)])]\n",
    "print(\"# odd spin structures:\", 2^(g-1)*(2^g-1))\n",
    "print(\"# even spin structures:\", 2^(g-1)*(2^g+1),\"\\n\")\n",
    "\n",
    "# Use get orbits to have the orbit decomposition.\n",
    "print(\"Partitions:\")\n",
    "for parity in [odds, evens]:\n",
    "    SS_Perms_exact = [[parity.index(A.transpose()*xi+V) for xi in parity] for A, V in zip(Mbars, Vbars)]\n",
    "    SS_orbits_exact = get_orbits(SS_Perms_exact, len(parity))\n",
    "    lengths = [len(o) for o in SS_orbits_exact]\n",
    "    print(sum(lengths),\"=\",lengths)\n",
    "\n",
    "print()\n",
    "# We can also intersect the corresponding affine subspaces. \n",
    "MIs = [M-I2g for M in Mbars]\n",
    "X0s = [MI.solve_left(-V) for MI, V in zip(MIs, Vbars)]\n",
    "Kers = [MI.kernel() for MI in MIs]\n",
    "AS = [AffineSubspace(X0, K) for X0, K in zip(X0s, Kers)]\n",
    "Sols = AS[0]\n",
    "for AI in AS[1:]:\n",
    "    Sols = Sols.intersection(AI)\n",
    "print(\"Space of invariant spin structures:\", Sols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "removing redundancy\n",
      "warning, next part of cell will take approximately 2.9884307861328123 minutes to run\n",
      "256 = [1, 30, 10, 5, 30, 10, 20, 60, 5, 20, 20, 5, 10, 30]\n"
     ]
    }
   ],
   "source": [
    "# To do the calculation we need some legwork completed in the automorphism group notebook,\n",
    "# namely to calculate the automorphism group in the HC model as a group of substitutions. \n",
    "# See that notebook for more of a discussion of this first bit of code. \n",
    "\n",
    "# Because we will want this Riemann surface to have the AJ map, we will redefine. \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 = 60\n",
    "S = RiemannSurface(f, prec=prec, integration_method=integration_method)\n",
    "Mlist = S.symplectic_isomorphisms()\n",
    "\n",
    "#####\n",
    "K.<z> = CyclotomicField(5)\n",
    "Tnum = S.tangent_representation_numerical(Mlist)\n",
    "zc = S._CC(z)\n",
    "indeps = [zc^m for m in range(z.multiplicative_order()-1)]\n",
    "def lincomb(v):\n",
    "    ls = [v]+indeps\n",
    "    w = vector(K, pari.lindep(ls).sage()) \n",
    "    if w[0] != 0 :\n",
    "        return K(list(-w[1:]/w[0]))    \n",
    "Talg = [matrix(K, g, g, [lincomb(a) for a in tnum.list()]) for tnum in Tnum]\n",
    "L = PolynomialRing(K, names=['x','y'])\n",
    "vB = vector([L(b) for b in S.cohomology_basis()]) \n",
    "def mat_to_subs(M):\n",
    "    w = M*vB\n",
    "    yp = w[3]/w[0]\n",
    "    xp = (w[1]+yp*w[3])/w[2]\n",
    "    return (xp, yp)\n",
    "auts_as_subs = [mat_to_subs(talg) for talg in Talg]\n",
    "\n",
    "half_len = len(auts_as_subs)/2\n",
    "if all([auts_as_subs.index(auts_as_subs[-i])==half_len-i for i in range(1, half_len)]):\n",
    "    print(\"removing redundancy\")\n",
    "    auts_as_subs = auts_as_subs[:half_len]   \n",
    "####\n",
    "# We calculate the AJ map of a certain theta characteristic as done in the Weierstass notebook\n",
    "# Now compute the canonical divisor (using a known theta characteristic)\n",
    "F = Curve(f).function_field()\n",
    "divx = F(x).divisor()\n",
    "c, d, b, a = sorted(divx.support(), key = lambda si: divx.valuation(si))\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",
    "\n",
    "# We store the point that is the basepoint of the Abel-Jacobi map \n",
    "bpCC = (S._vertices[S._basepoint[0]], S._wvalues[S._basepoint[0]][S._basepoint[1]])\n",
    "# We will calculate those elements of Tnum that are indeed the pullback of an automorphism\n",
    "# on the differentials, and not this pullback multiplied -1 in the Jacobian. These are \n",
    "# equivalently the orientation preserving operations. \n",
    "Tnum_op = []\n",
    "\n",
    "# We expect the output of the symplectic isomorphisms to be s.t. T[i]+T[120+i]==0, \n",
    "# and we rely on this for the next calculation, so we check it. \n",
    "if not all([Talg[i]+Talg[120+i]==0 for i in range(120)]):\n",
    "    raise ValueError(\"Matrices are not distributed as expected\")\n",
    "\n",
    "# We pick out the correct matrices in Tnum as described below\n",
    "for sigma, L1, L2 in zip(auts_as_subs, Tnum[:120], Tnum[120:]):\n",
    "    sigma_bp = tuple(sigmaj(*bpCC) for sigmaj in sigma)\n",
    "    tA1 = (L1-Ig)*2*AJ_bp_D+2*(g-1)*S._aj_based(sigma_bp)\n",
    "    tA2 = (L2-Ig)*2*AJ_bp_D+2*(g-1)*S._aj_based(sigma_bp)\n",
    "    if S.reduce_over_period_lattice(tA1).norm()<1e-10:\n",
    "        Tnum_op.append(L1)\n",
    "    elif S.reduce_over_period_lattice(tA2).norm()<1e-10:\n",
    "        Tnum_op.append(L2)\n",
    "    else:\n",
    "        raise ValueError(\"neither transform is orientation preserving\")\n",
    "\n",
    "if (len(Tnum_op)!=120):\n",
    "    raise ValueError(\"Incorrect number of automorphisms found.\")\n",
    "\n",
    "    \n",
    "# As shown below we need to only calculate the transform of the half-lattice vectors. \n",
    "V = VectorSpace(S._CC, 2*g)\n",
    "XS = [V(vector(Integer(n).digits(2, padto=2*g), F2)) for n in range(4^g)]\n",
    "PM  = S.period_matrix()\n",
    "half_lattices = [PM*xi/2 for xi in XS]\n",
    "\n",
    "# The following function will give a way to find which half lattice one is transformed to \n",
    "def index_of_transformed_hl(hl, L):\n",
    "    hl2 = L*hl\n",
    "    ds = [S.reduce_over_period_lattice(hl2-hls).norm() for hls in half_lattices]\n",
    "    # If using checks, insert a check clause here, for now we assume these have all worked\n",
    "    return ds.index(min(ds))\n",
    "\n",
    "# To generate the whole group we will just take the action of R, S, and U\n",
    "# We need the matrices corresponding to R and S, and then their associated transform\n",
    "K.<z> = CyclotomicField(5)\n",
    "rs = (polygen(K)^2-5).roots(multiplicities=False)\n",
    "s5 = CC(sqrt(5))\n",
    "remainders = [(K.embeddings(CC)[0](r)-s5).abs() for r in rs]\n",
    "s5zs = rs[remainders.index(min(remainders))]\n",
    "# To avoid overlap of notation we will call them Rp and Sp\n",
    "Sp = Matrix([[1,0,0],[0,z^(1),0],[0,0,z^(-1)]])\n",
    "Rp = Matrix([[1,2,2],[1,z^2+z^(-2),z^1+z^(-1)],[1,z^1+z^(-1),z^2+z^(-2)]])/s5zs\n",
    "# Now we have the matrices get the action. \n",
    "L = PolynomialRing(K, names=['x','y'])\n",
    "RP_action = Rp*vector([L.gen(0), L.gen(1), 1])\n",
    "RP_as_subs = (RP_action[0]/RP_action[2], RP_action[1]/RP_action[2])\n",
    "SP_action = Sp*vector([L.gen(0), L.gen(1), 1])\n",
    "SP_as_subs = (SP_action[0]/SP_action[2], SP_action[1]/SP_action[2])\n",
    "# We get U as our simplest non-linear as before\n",
    "non_lin = []\n",
    "for i in range(120):\n",
    "    xp, yp = auts_as_subs[i]\n",
    "    xnum = xp.numerator()\n",
    "    xden = xp.denominator()\n",
    "    ynum = yp.numerator()\n",
    "    yden = yp.denominator()\n",
    "    nds = [xnum, xden, ynum, yden]\n",
    "    if not all([ndi.degree()<=1 for ndi in nds]):\n",
    "        non_lin.append(auts_as_subs[i])\n",
    "if len(non_lin) != 60:\n",
    "    raise ValueError    \n",
    "lengths = [len(aut[0]._repr_())+len(aut[1]._repr_()) for aut in non_lin]\n",
    "nl = non_lin[lengths.index(min(lengths))]\n",
    "\n",
    "# With this we can collect the information\n",
    "RSU_as_subs = [RP_as_subs, SP_as_subs, nl]\n",
    "RSU_as_matrices = [Tnum_op[auts_as_subs.index(geni)] for geni in RSU_as_subs]\n",
    "\n",
    "# We take note of the time we expect the cell will take to run to be responsible. \n",
    "import time\n",
    "ct = time.time()\n",
    "index_of_transformed_hl(half_lattices[0], RSU_as_matrices[0])\n",
    "ct2 = time.time()\n",
    "print('warning, next part of cell will\\\n",
    " take approximately {} minutes to run'.format((3*len(half_lattices)*(ct2-ct))/60))\n",
    "\n",
    "HL_Perms_num = [[index_of_transformed_hl(hli, Li) for hli in half_lattices] \n",
    "                for Li in RSU_as_matrices]\n",
    "# This next line takes a while \n",
    "nSS = 4^g\n",
    "HL_orbits_num = get_orbits(HL_Perms_num, nSS)\n",
    "lengths = [len(o) for o in HL_orbits_num]\n",
    "print(sum(lengths),\"=\",lengths)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's worth writing down what is going on in the above calculation. We have \n",
    "$$\n",
    "A_\\ast(K) = \\int_\\ast^K {\\omega}\n",
    "$$\n",
    "and so when we act on the left with $L$ the matrix corresponding to automorphism $\\sigma$\n",
    "$$\n",
    "L A_\\ast(K) = \\int_\\ast^K L \\omega = \\int_\\ast^K \\sigma^\\ast \\omega = \\int_{\\sigma(\\ast)}^{\\sigma(K)} \\omega = \\int_\\ast^{\\sigma(K)} \\omega - (\\deg K) \\int_\\ast^{\\sigma(\\ast)} \\omega \n",
    "$$\n",
    "and hence \n",
    "$$\n",
    "(L-I)A_\\ast(K)+2(g-1)A_\\ast(\\sigma(\\ast)) = \\int_K^{\\sigma(K)} \\omega\n",
    "$$\n",
    "Then $\\sigma(K) \\sim K$ iff the RHS is 0 in the Jacobian. As an automorphism must preserve the canonical divisor, this let's us pick out the correct matrices. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now given a spin structure $D$ s.t. $2D=K$ we should have $2A_\\ast(D) = A_\\ast(2D) = A_\\ast(K)$ and to the AJ images of the spin structures can differ only by a half period $h = \\Omega n$, that is $A_\\ast(D) = A_\\ast(\\Delta)+h$. Now as spin structures are mapped to each other we expect $\\sigma(D) = D^\\prime$, and so \n",
    "$$\n",
    "L A_\\ast(D) = \\int_{\\sigma(\\ast)}^{\\sigma(D)} \\omega = A_\\ast(D^\\prime) - (g-1) A_\\ast(\\sigma(\\ast))\n",
    "$$\n",
    "\n",
    "\n",
    "As in the previous case we see \n",
    "$$\n",
    "L [A_\\ast(\\Delta)+h] = [LA_\\ast(\\Delta)] + Lh = [A_\\ast(\\Delta) - (g-1) A_\\ast(\\sigma(\\ast))] + Lh = A_\\ast(\\Delta) + [Lh - (g-1) A_\\ast(\\sigma(\\ast))]\n",
    "$$\n",
    "Thus in order to find the transform of the spin structures we need to just look at the transforms of the half-lattice vectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theorem 5.11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# The above numerical calculation was certainly much more involved than that of KS06, and slower. \n",
    "# The main benefit is that it will let us identify whether a candidate theta characteristic is invariant\n",
    "# directly, which KS06 will not. We demonstrate this by proving theorem 5.11\n",
    "for sigma, Li in zip(RSU_as_subs, RSU_as_matrices):\n",
    "    sigma_bp = tuple(sigmaj(*bpCC) for sigmaj in sigma)\n",
    "    tA = (Li-Ig)*AJ_bp_D+(g-1)*S._aj_based(sigma_bp)\n",
    "    print((S.reduce_over_period_lattice(tA).norm()<1e-10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proposition 5.13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.6684134041982483463310551740431753521118841261082271773679e-59\n",
      "3.5763230795065398350081479773261584288238596070347201420422e-58\n",
      "5.2958839222882894254427149485302839233349623611678864662229e-24\n"
     ]
    }
   ],
   "source": [
    "# We now want to verify that, in the RR homology basis, we have that the unique invariant theta \n",
    "# characteristic is the divisor of the Szego kernel. \n",
    "# This will require high-precision calculation, so we will repeat lots of the calculation above. \n",
    "from sage.schemes.riemann_surfaces.riemann_surface import numerical_inverse\n",
    "\n",
    "# Initialise the curve \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",
    "g = S.genus\n",
    "\n",
    "# Define our canonical divisor. \n",
    "K_can = F(y^3-x).divisor() + F(x).differential().divisor() - F(f.derivative(y)).divisor()\n",
    "# We make some check to ensure a, b, c, d, have been assigned correctly\n",
    "if K_can-3*a-2*c-1*d:\n",
    "    raise ValueError(\"a,b,c,d have been assigned incorrectly\")\n",
    "    \n",
    "# We recalculate the period matrix and normalise, as this is required for the RCV\n",
    "PM = S.period_matrix()\n",
    "AM = PM[:,0:g]\n",
    "AInv = numerical_inverse(AM)\n",
    "PM_normalised = AInv*PM\n",
    "\n",
    "# We can calculate all possible half-periods in the normalised lattice\n",
    "V = VectorSpace(S._CC, 2*g)\n",
    "XS = [V(vector(Integer(n).digits(2, padto=2*g), F2)) for n in range(4^g)]\n",
    "half_lattices_normalised = [PM_normalised*xi/2 for xi in XS]\n",
    "K_can_avoid, _ = S.strong_approximation(K_can, avoid)\n",
    "K_can_list = S.divisor_to_divisor_list(K_can_avoid)\n",
    "AJ_bp_can = S.abel_jacobi(K_can_list)\n",
    "\n",
    "# Let's set the basepoin to be a, and normalise\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",
    "AJ_a_can = AJ_bp_can - 2*(g-1)*AJ_bp_a\n",
    "AJ_a_can_normalised = AInv*AJ_a_can\n",
    "\n",
    "# Use these to create a list of candidate RCVs using the fact that 2K == A(can). \n",
    "RCV_a0 = AJ_a_can_normalised/2\n",
    "RCV_a_list = [RCV_a0+hl for hl in half_lattices_normalised]\n",
    "\n",
    "# Create the theta function\n",
    "from riemann_theta.riemann_theta import RiemannTheta\n",
    "Tau = AInv*PM[:,g:]\n",
    "RT = RiemannTheta(Tau)\n",
    "\n",
    "# We evaluate the theta function at the RCV (with the correct sign for convention) for all candidates.\n",
    "# We then record the index of the value which gives the smallest value, as this will be the correct one. \n",
    "# Note we need to set normaliesd=True as this is the condition required on the homology basis. \n",
    "Thetas = [RT(S.reduce_over_period_lattice(-Ki, normalised=True)).norm() for Ki in RCV_a_list]\n",
    "mT = min(Thetas)\n",
    "# make a check of accuracy:\n",
    "if mT>1e-10:\n",
    "    raise ValueError(\"A root of theta unlikely to be found\")\n",
    "ind = Thetas.index(mT)\n",
    "RCV_a = RCV_a_list[ind]\n",
    "\n",
    "# We now want to loop over homology transforms\n",
    "# we copy the code from the period matrix notebook, with fewer comments\n",
    "####\n",
    "from sage.geometry.hyperbolic_space.hyperbolic_model import moebius_transform\n",
    "j0 = -29^3*5*2^(-5)\n",
    "alp = S._CCz.gen()\n",
    "roots = (4*j0*alp*(1-alp)-1728).roots(multiplicities=False)\n",
    "t1, t2 = (I*hypergeometric([1/6,5/6],[1],1-al).hypergeometric_simplify()/hypergeometric([1/6,5/6],[1],al).hypergeometric_simplify() \n",
    "          for al in roots)\n",
    "if (t1.real_part()-1/2).abs()>1e-10:\n",
    "    print(\"Swapping roots\")\n",
    "    t1 = t2\n",
    "if (t1.real_part()-1/2).abs()>1e-10:\n",
    "    raise ValueError(\"Invalid taus found, check calculation.\")\n",
    "t0 = moebius_transform(matrix(2,2,[-1,1,-2,1]), S._CC(t1)) - 1\n",
    "if ((elliptic_j(t0)+29^3*5/2^5).norm()+(elliptic_j(5*t0)+25/2).norm())>1e-10:\n",
    "    raise ValueError(\"Incorrect tau0 found\")\n",
    "RM_Braden = t0*matrix([[4,1,-1,1],[1,4,1,-1],[-1,1,4,1],[1,-1,1,4]])\n",
    "Rs = integer_matrix_relations(Tau, RM_Braden)\n",
    "r = len(Rs)\n",
    "A = PolynomialRing(QQ, r, 'x')\n",
    "gensA = A.gens()\n",
    "R = sum(gensA[i] * Rs[i].change_ring(A) for i in range(r))\n",
    "tr = (R * S.rosati_involution(R)).trace()\n",
    "M = Matrix(ZZ, r, r, [tr.derivative(gen1).derivative(gen2)\n",
    "                        for gen1 in gensA for gen2 in gensA])\n",
    "vs = M.__pari__().qfminim(4*g)[2].sage().transpose()\n",
    "vs = [v for v in vs if v * M * v == 4*g]\n",
    "vs += [-v for v in vs]\n",
    "RsIso = []\n",
    "for v in vs:\n",
    "    R = sum(v[i] * Rs[i] for i in range(r))\n",
    "    if R * S.rosati_involution(R) == 1:\n",
    "        RsIso.append(R)\n",
    "if len(RsIso)!=240:\n",
    "    raise ValueError(\"Not enough transforms found, error somewhere.\")\n",
    "####\n",
    "\n",
    "# We recreate the function that reduces over the period lattice, now for an arbitrary complex torus\n",
    "# of dimension g\n",
    "def reduce_over_lattice(vector, Lambda):\n",
    "    VR = VectorSpace(S._RR, 2*S.genus)\n",
    "    VC = VectorSpace(S._CC, S.genus)\n",
    "    I = S._CC(0,1)\n",
    "    PM = Lambda        \n",
    "    def C2R(v):\n",
    "        return VR([z.real_part() for z in v]+[z.imag_part() for z in v])        \n",
    "    u = C2R(vector)\n",
    "    basis_vecs = [C2R(c) for c in PM.columns()]\n",
    "    M = Matrix([[ei.dot_product(ej) for ei in basis_vecs] for ej in basis_vecs])\n",
    "    v_dot_e = VR([u.dot_product(e) for e in basis_vecs])\n",
    "    coeffs = M.solve_right(v_dot_e)\n",
    "    u -= sum([t.round()*e for t, e in zip(coeffs, basis_vecs)])\n",
    "    reduced = VC([S._CC(u[i]+I*u[i+S.genus]) for i in range(S.genus)])    \n",
    "    return reduced\n",
    "\n",
    "# Calculate the AJ image of the invariant theta characteristic and normalise \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",
    "AJ_a_D = AJ_bp_D - Delta.degree()*AJ_bp_a\n",
    "AJ_a_D_normalised = AInv*AJ_a_D\n",
    "\n",
    "\n",
    "# We now loop over the homology changes, checking the valid ones against Braden's result, \n",
    "# and then comparing AJ images\n",
    "vers = []\n",
    "for R0 in RsIso:\n",
    "    D0 = R0[0:g, 0:g]\n",
    "    B0 = R0[0:g, g:]\n",
    "    C0 = R0[g:, 0:g]\n",
    "    A0 = R0[g:, g:]\n",
    "    if (max([el.norm() for el in (B0+Tau*A0-(D0+Tau*C0)*RM_Braden).list()])>1e-10):\n",
    "        raise ValueError(\"Homoloy transform is invalid\")\n",
    "    mat1 = numerical_inverse(C0.transpose()*Tau+D0.transpose())\n",
    "    mat2 = A0.transpose()*Tau+B0.transpose()\n",
    "    char1 = vector([(C0.transpose()*D0)[j,j] for j in range(g)])/2\n",
    "    char2 = vector([(A0.transpose()*B0)[j,j] for j in range(g)])/2\n",
    "    RCV_a_RR = RCV_a*mat1 - char1*mat2*mat1-char2\n",
    "\n",
    "    # We check that this still satisfies the theta condition.\n",
    "    RTRR = RiemannTheta(RM_Braden)\n",
    "    if RTRR(-RCV_a_RR).norm()>1e-10:\n",
    "        raise ValueError(\"Error in transform of RCV\")\n",
    "    vers.append(reduce_over_lattice(vector([-12,-3,3,-3])/5+t0*vector([-6,-6,3,0])-2*RCV_a_RR, \n",
    "                          block_matrix(1,2,[Ig, RM_Braden])).norm())\n",
    "    \n",
    "for i in range(len(RsIso)):\n",
    "    if vers[i]>1e-10:\n",
    "        continue\n",
    "    R0 = RsIso[i]\n",
    "    D0 = R0[0:g, 0:g]\n",
    "    B0 = R0[0:g, g:]\n",
    "    C0 = R0[g:, 0:g]\n",
    "    A0 = R0[g:, g:]\n",
    "    mat1 = numerical_inverse(C0.transpose()*Tau+D0.transpose())\n",
    "    mat2 = A0.transpose()*Tau+B0.transpose()\n",
    "    char1 = vector([(C0.transpose()*D0)[j,j] for j in range(g)])/2\n",
    "    char2 = vector([(A0.transpose()*B0)[j,j] for j in range(g)])/2\n",
    "    RCV_a_RR = RCV_a*mat1 - char1*mat2*mat1-char2\n",
    "\n",
    "    # We check thath this still satisfies the theta condition.\n",
    "    RTRR = RiemannTheta(RM_Braden)\n",
    "    if RTRR(-RCV_a_RR).norm()>1e-10:\n",
    "        raise ValueError(\"Error in transform of RCV\")\n",
    "        \n",
    "    t3 = (reduce_over_lattice(vector([3,2,-2,-3])/10+t0.imag_part()*I*vector([1,-2,-2,1])-RCV_a_RR, \n",
    "                          block_matrix(1,2,[Ig, RM_Braden])).norm())\n",
    "    tt = (reduce_over_lattice(10*RCV_a_RR, block_matrix(1,2,[Ig, RM_Braden])).norm())\n",
    "    \n",
    "    tl = (reduce_over_lattice(AJ_a_D_normalised*mat1-RCV_a_RR, block_matrix(1,2,[Ig, RM_Braden])).norm())\n",
    "    \n",
    "    print(t3)\n",
    "    print(tt)\n",
    "    print(tl)\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.4.rc0",
   "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}