CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

| Download
Views: 558
License: MIT
Image: ubuntu2004
Kernel: SageMath 9.4
# This notebook is associated to the paper "Counting points on smooth plane quartics" # by Edgar Costa, David Harvey, and Andrew V. Sutherland # It includes a demonstration implementation of Algorithms 5.8,5.11,5.14 # For algorithm 5.16 we give only the simple O(p log p loglog p) implementation # we also implement the uncompressed variation, see remark 5.17. # For the average polynomial time algorithms we implement a simple remainder tree, not a remainder forest (so kappa=0), # We implement: # - the compressed version algorithm 5.22 with remark 5.26 incorporated # - the uncompressed version, see remark 5.24 # - optimised Harvey's algorithm with remark 5.29 incorporated. # We only address primes p <= N that do not divide D from sage.all import ( GCD, GF, Matrix, PolynomialRing, ZZ, cached_method, identity_matrix, inverse_mod, is_prime, polygen, prime_range, prod, vector, ) def decrease(w, i): return tuple(c - int(i==j) for j, c in enumerate(w)) @cached_method def monomials(d): return [ (e0, e1, e2) for e0 in range(d + 1) for e1 in range(d + 1) for e2 in range(d + 1) if e0 + e1 + e2 == d ] def invfactorialproduct(R, y): if any(elt not in ZZ or elt < 0 for elt in y): return 0 y = [elt for elt in y if y != 0] if not y: return 1 return 1/prod(prod( R(i) for i in range(1, elt + 1) ) for elt in y) def random_int(bits): return ZZ.random_element(2 ** (bits + 1) - 1) - 2**bits + 1 def random_quartic(bits): x0, x1, x2 = PolynomialRing(ZZ, "x", 3).gens() return sum( [ random_int(bits) * x0**e0 * x1**e1 * x2**e2 for (e0, e1, e2) in monomials(4) ] ) def art(M, m, c=1): """Computes [prod([M[i] for i in range(k+1)])/c^k % m[k] for k in range(len(M))] assuming M[i]*M[i+1] is divisible by c""" N = len(M) if N <= 1: return [M[0] % m[0]] if N == 1 else [] n = (N + (N % 2)) // 2 M.append(c) m.append(1) # pad C = art( [(M[2 * k - 1] if k else c) * M[2 * k] / c for k in range(n)], [m[2 * k] * m[2 * k + 1] for k in range(n)], c=c, ) C = sum( [ [ C[k] % m[2 * k], C[k] * M[2 * k + 1] * inverse_mod(c, m[2 * k + 1]) % m[2 * k + 1], ] for k in range(n) ], [], )[:N] return C class Quartic: # F = homogeneous quartic in Z[x0,x1,x2] def __init__(self, F): self.F = F self.d = d = F.degree() assert d == 4 self.ZZV = PolynomialRing(ZZ, 3, "V") (V0, V1, V2) = V = self.ZZV.gens() self.D = [monomials(i) for i in range(2*d)] self.e = [vector(elt) for elt in reversed(sorted(self.D[1]))] # for any v with |v| = 4m+6, # let b(v) = ((F^(n+1))_{v-u})_{u in T2}, ((m + 1) (F^n)_{v-u})_{u in B6} (vector in Z^16) # where B6 is a certain 10-element subset of T6. # note that ell=2d-2=6 in the paper # compute: # - the set B6 # - the nonzero integer lambda6, # - a 16*28 integer matrix PI # - a 28*16 matrix matrix PSI # PSI over Z[V0,V1,V2] (with entries degree at most 1) such that # for any |v| = 4m+6 we have # PI a(v) = b(v) # PSI(v) b(v) = d (m + 1) lambda6 a(v). # In particular, # PI PSI(v) = d (m + 1) lambda6 I_{16} (identity matrix) # PSI(v) PI a(v) = d (m + 1) lambda_6 a(v) # (the idea is that PI compresses a length 28 vector into a length 16 vector, # and PSI(v) performs the reverse operation) # compute the set B6 T = Matrix(ZZ, 18, 28) for i in range(3): for u_index, u in enumerate(self.D[2]): #u = self.D[2][u_index] for t in self.D[4]: j = self.D[6].index((u[0] + t[0], u[1] + t[1], u[2] + t[2])) T[6 * i + u_index, j] = t[i] * F[t] B6 = T.nonpivots() J = Matrix(ZZ, 10, 28) for i, v in enumerate(B6): J[i, v] = 1 S = T.stack(J); Sadj = S.adjugate() Sadj /= GCD(Sadj.list()) lambda6I = Sadj*S lambda6 = ZZ(lambda6I[0,0]) assert lambda6I == lambda6*identity_matrix(28) PI = Matrix(ZZ, 16, 28) for u_index in range(6): u = self.D[2][u_index] for t in self.D[4]: PI[u_index, self.D[6].index((u[0] + t[0], u[1] + t[1], u[2] + t[2]))] = F[ t ] for j in range(10): PI[j + 6, B6[j]] = 1 self.PI = PI m = -2 # taking advantage that m = p - 2 = -2 mod p PSI = Matrix(self.ZZV, 28, 16) V = self.ZZV.gens() for u_index, u in enumerate(self.D[6]): # how to write x^u row = Sadj[u_index] # the term in D_2 for t_index, t in enumerate(self.D[2]): PSI[u_index, t_index] = sum((V[i] - t[i]) * row[6*i + t_index] for i in range(3)) # the term in B_6 for b_index, b in enumerate(B6): PSI[u_index, 6 + b_index] = (m + 1) * row[6*3 + b_index] self.PSI = PSI # self.PSIeval = PSI(0, 2*p-1, 2*p - 1) mod p self.PSIeval = PSI.substitute({V0:0, V1: -1, V2: - 1}).change_ring(ZZ) self.c = lambda6 # this is the c in Remark 5.15 self.theta0 = self.phi_generic(0, 1)[0] self.theta1 = self.phi_generic(1, 0)[0] self.denominator = 2 * lambda6 * self.theta0 * self.theta1 @cached_method def phi_generic(self, i, j): d = self.d ell = 2*d-2 A = Matrix(2*d, 2*d) Mconst, Mi, Mk = [Matrix(2*d, (3*d - 1)*d//2) for _ in range(3)] k = (-i - j) % 3 m = -2 # taking advantage that m = p - 2 = -2 mod p U = sorted([tuple(a * self.e[k] - ((a + 1) - 4) * self.e[j]) for a in range(d)]) S = sorted([tuple(c * self.e[k] - (c + 1 - 8) * self.e[j]) for c in range(2*d)]) # T = list({tuple(vector(u) - vector(t)) for u in U for t in self.D[4]}.difference(S)) T = sorted({tuple(vector(u) + vector(t)) for u in U for t in self.D[d] if t[i] != 0}) assert len(T) == (3*d - 1)*d//2 for a, u in enumerate(U): for t, ft in self.F.dict().items(): uplust = tuple(vector(u) + vector(t)) if t[i] == 0: Gcoordinate = S.index(uplust) A[2 * a, Gcoordinate] = ft A[2 * a + 1, Gcoordinate] = t[k] * ft else: Gcoordinate = T.index(uplust) # M[2*l][Gcoordinate] = ((m + 1)*t[i] - w[i])*f[t] # = ((m + 1)*t[i] - (v[i] + 1))*f[t] Mconst[2 * a, Gcoordinate] = ((m + 1) * t[i] - 1) * ft Mi[2 * a, Gcoordinate] = -ft # Mk[2*l][Gcoordinate] = 0 # M[2*l + 1][Gcoordinate] = # = (w[k]*t[i] - w[i]*t[k])*f[t] # = ((v[k] - u[k])*t[i] - (v[i] + 1)*t[k])*f[t] Mconst[2 * a + 1, Gcoordinate] = (-u[k] * t[i] -t[k]) * ft Mi[2 * a + 1, Gcoordinate] = -t[k] * ft Mk[2 * a + 1, Gcoordinate] = t[i] * ft # Bdegree[2*l+1] = 0 thetai = A.determinant() adjA = A.adjugate() Bconst, Bi, Bk = [adjA * elt for elt in [Mconst, Mi, Mk]] PHIconst, PHIi, PHIk = [Matrix(len(self.D[ell + 1]), len(self.D[ell])) for _ in range(3)] for s_coordinate, s in enumerate(self.D[ell + 1]): if s[i] > 0: # i.e., s - ei in D(v, 6) si = decrease(s, i) si_coordinate = self.D[ell].index(si) PHIconst[s_coordinate, si_coordinate] = thetai PHIi[s_coordinate, si_coordinate] = thetai if s[i] == 0: assert s in S B_row = S.index(s) for j, t in enumerate(T): ti = decrease(t, i) ti_coordinate = self.D[ell].index(ti) PHIconst[s_coordinate, ti_coordinate] = Bconst[B_row, j] PHIi[s_coordinate, ti_coordinate] = Bi[B_row, j] PHIk[s_coordinate, ti_coordinate] = Bk[B_row, j] return thetai, PHIconst, PHIi, PHIk @cached_method def tau_generic(self, i, j): thetai, PHIconst, PHIi, PHIk = self.phi_generic(i,j) projT = Matrix([[ int(t==decrease(w,j)) for w in self.D[7]] for t in self.D[6]]) tauconst, taui, tauk = [projT*elt for elt in [PHIconst, PHIi, PHIk]] return thetai, tauconst, taui, tauk def tau(self, v, i, j): k = (-i - j) % 3 thetai, tauconst, taui, tauk = self.tau_generic(i, j) return thetai, tauconst + taui*v[i] + tauk*v[k], taui def T(self, v, i, j): thetai, tauconst, taui, tauk = self.tau_generic(i, j) V = self.ZZV.gens() k = (-i - j) % 3 tau = tauconst + taui*V[i] + tauk*V[k] calT = self.PI.change_ring(self.ZZV)*tau*self.PSI x = polygen(ZZ) v = list(v) v[i] += x v[j] -= x calTv = calT(*v) T = [Matrix(16,16) for _ in range(3)] for r in range(16): for c in range(16): for d in range(3): T[d][r,c] = calTv[r,c][d] return T @cached_method def Tk_gen(self): return self.T((0,-1,-1), 0, 1) def Tk(self, k): T0, T1, T2 = self.Tk_gen() # we realign the k variable so that # newk = p - 2 - k, so that Tk(0) ... Tk(p-2) # matches with T_{w + (p-2)(e0 -e1), 0, 1} ... T_{w, 0, 1} newk = -2 -k return T0 + newk*T1 + newk**2*T2 def tauk(self, k): _, T0, T1 = self.tau((0,-1,-1), 0, 1) # we realign the k variable so that # newk = p - 2 - k, so that tauk(0) ... tauk(p-2) # matches with tau_{w + (p-2)(e0 -e1), 0, 1} ... tau_{w, 0, 1} newk = -2 -k return T0 + newk*T1 @cached_method def equation_matrix(self, v): ell = 2*self.d-2 m = -2 M = Matrix(2*len(self.D[ell - self.d]), len(self.D[ell])) for s_coordinate, s in enumerate(self.D[ell - self.d]): for t, ft in self.F.dict().items(): spt = tuple(vector(s) + vector(t)) j = self.D[ell].index(spt) for i in range(2): M[2*s_coordinate + i, j] = ((v[i] - s[i]) - (m + 1) * t[i])*ft return M @cached_method def fpower_naive(self, p): return self.F.change_ring(GF(p))**(p-2) def naive_f_at_v(self, p, v): fm = self.fpower_naive(p) ell = sum(v) - self.d*(p-2) coordinates = [tuple(vector(v) - vector(elt)) for elt in self.D[ell]] return vector([fm[elt] if all(m >= 0 for m in elt) else 0 for elt in coordinates]) def cm_naive(self, p): """Computes the Cartier-Manin matrix diretly from the definition""" Fpow = self.F.change_ring(GF(p)) ** (p - 1) # Apply equation (2.6) in the paper return Matrix( [ [Fpow[p-1, p-1, 2*p-2], Fpow[2*p-1, p-1, p-2], Fpow[p-1, 2*p-1, p-2]], [Fpow[p-2, p-1, 2*p-1], Fpow[2*p-2, p-1, p-1], Fpow[p-2, 2*p-1, p-1]], [Fpow[p-1, p-2, 2*p-1], Fpow[2*p-1, p-2, p-1], Fpow[p-1, 2*p-2, p-1]] ] ) # This is Algorithm 5.8 in the paper with steps (1) and (2) interleaved def cm_final(self, p, Cp): """Given Cp for a prime p not dividing D, computes the Cartier-Manin matrix Cp""" K = GF(p) A = Matrix(K, 3, 3) R = PolynomialRing(K, "x") F = self.F assert Cp.ncols() == Cp.nrows() and Cp.ncols() in [16, 28] compressed = Cp.ncols() == 16 # (a) compute 1st column of A # step 1a # compute a1 = a(w_1) = a(0, 2p-1, 2p-1) g = R([F[0, 4, 0], F[0, 3, 1], F[0, 2, 2], F[0, 1, 3], F[0, 0, 4]]) ** 2 # compute starting vector F^(p-2)|D(w1,6) via Section 4.2 using g=F(0,1,x)^2 # using coefficients of series (J/J0)^(-1) = 1 + G1*x + G2*x^2 + ... Q = Matrix(K, 8, 8) for j in range(8): Q[0, j] = -g[j + 1] / g[0] for j in range(7): Q[j + 1, j] = 1 Qpow = Q ** (p - 1) Qpow2 = Qpow * Qpow * Q # M^(2p-1) # deduce (F^(p-2))_{0, 2p-1-(6-j), 2p-1-j} for 0 <= j <= 6 a1 = vector(K, 28) for j in range(7): u = (0, 6 - j, j) a1[self.D[6].index(u)] = ( Qpow[j, 0] * F[0, 3, 1] / F[0, 4, 0] ** 2 + Qpow2[j, 0] / F[0, 4, 0] ) # step 2a # temp = b(p-1, p, 2p-1) if compressed: # compress it by using PI temp = Cp * self.PI * a1 if not compressed: # go from F^(p-2) to F^(p-1) by applying PI temp = self.PI * Cp * a1 # first col of CM matrix A[0, 0] = temp[self.D[2].index((0, 1, 1))] A[1, 0] = temp[self.D[2].index((1, 1, 0))] A[2, 0] = temp[self.D[2].index((0, 2, 0))] # (b) compute 2nd column of A # step 1b # compute a2 = a(w_2) = a((3p-1, 0, p-1)) # compute starting vector F^(p-2)|D(w2,6) via Section 4.2 using g=F(1,0,x)^2 g = R([F[4, 0, 0], F[3, 0, 1], F[2, 0, 2], F[1, 0, 3], F[0, 0, 4]]) ** 2 Q = Matrix(K, 8, 8) for j in range(8): Q[0, j] = -g[j + 1] / g[0] for j in range(7): Q[j + 1, j] = 1 Qpow_2 = Q ** (p - 1) # deduce (F^(p-2))_{3p-1-(6-j), 0, p-1-j} for 0 <= j <= 6 a2 = vector(K, 28) for j in range(7): u = (6 - j, 0, j) a2[self.D[6].index(u)] = Qpow_2[j, 0] / F[4, 0, 0] # step 2b # temp = b(v_2) = b((2p, p-1, p-1)) # (NOTE: in this next line, we don't need the matrix inverse; just the solution to the system is enough) if compressed: # compress it by using PI temp = Cp.solve_right(self.PI * a2) if not compressed: # go from F^(p-2) to F^(p-1) by applying PI # note: Cp is not invertible, but a2 is in its row space # so if we add the equations that temp must also solve, we get a unique solution M = Cp.stack(self.equation_matrix((0, -1, -1))) assert M.rank() == 28 temp = self.PI*M.solve_right(vector(GF(p), a2.list() + [0]*12)) # second col of CM matrix A[0, 1] = temp[self.D[2].index((1, 0, 1))] A[1, 1] = temp[self.D[2].index((2, 0, 0))] A[2, 1] = temp[self.D[2].index((1, 1, 0))] # (c) compute 3rd column of A # step 1c # compute a(w_3) = a((0, 3p-1, p-1)) # we already computed starting vector F^(p-2)|D(w3,6) via Section 4.2 using g=F(0,1,x)^2 in step 1a # deduce (F^(p-2))_{0, 3p-1-(6-j), p-1-j} for 0 <= j <= 6 a3 = vector(K, 28) for j in range(7): u = (0, 6 - j, j) a3[self.D[6].index(u)] = Qpow[j, 0] / F[0, 4, 0] # step 2c # temp = b(v_1) = b(p-1, 2p, p-1) = (p-1, 2p, p-1) if compressed: # compress it by using PI temp = Cp * self.PI * a3 if not compressed: # go from F^(p-2) to F^(p-1) by applying PI temp = self.PI * Cp * a3 # third col of CM matrix A[0, 2] = temp[self.D[2].index((0, 1, 1))] A[1, 2] = temp[self.D[2].index((1, 1, 0))] A[2, 2] = temp[self.D[2].index((0, 2, 0))] return A def good_primes(self, N): return [p for p in prime_range(N + 1) if self.denominator % p != 0] def mn(self, n): p = n + 2 return p if is_prime(p) and self.denominator % p != 0 else 1 def Cp(self, p, compressed=True): r""" return: let v = (0,-1,-1) = (0, 2p-1, 2p-1) mod p = (0, 3p-1, p-1) mod p if compressed: C[p] = - T_{v + (p-2)(e0 - e1), 0, 1} ... T_{v, 0, 1} else: C[p] = - tau_{v + (p-2)(e0 - e1), 0, 1} ... tau_{v, 0, 1} i.e., the transition matrix from b(0, 2p-1, 2p-1) to b(p-1, p, 2p-1) modulo p (the -1 comes from the denominator (p-1)! mod p) """ Mn = self.Tk if compressed else self.tauk return -prod([Mn(n).change_ring(GF(p)) for n in range(p - 1)]) def cartier_manin_matrices(self, N, algorithm="CHS", compressed=True): """returns dictionary p -> Cartier Manin matrix at p, for p < N not dividing D""" if algorithm == "naive": # This is the O(p log p loglog p) implementation of Algorithm 5.11 return {p: self.cm_final(p, self.Cp(p)) for p in self.good_primes(N - 1)} elif algorithm == "CHS": # This is a remainder tree implementation of Algorithm 5.14 that uses Remark 5.16 if compressed: M = [self.Tk(n) for n in range(N)] else: M = [self.tauk(n) for n in range(N)] m = [self.mn(n) for n in range(N - 1)] c = self.c if compressed else 1 C = art(M, m, c=c) return { p: self.cm_final(p, -C[p - 2].change_ring(GF(p)) / c) for p in self.good_primes(N) } elif algorithm == "Harvey": d = self.d h = 3*d-2 z = vector([0,0,h-d]) Bh = [vector(elt) for elt in monomials(h)] # interior of Newton polygon Bds = [vector([1,1,2]), vector([2,1,1]), vector([1,2,1])] V0 = Matrix(len(Bds), len(Bh)) # pick the 3 columns we want for i, elt in enumerate(Bds): V0[i, Bh.index(z + elt)] = 1 M = [self.Qk(n) for n in range(N)] M[0] = V0*M[0] m = [self.mn(n) for n in range(N - 1)] Q = art(M, m, c=1) return { p: self.cm_harvey(p, Q[p - 2].change_ring(GF(p))) for p in self.good_primes(N) } else: raise ValueError("not a valid algorithm option") @cached_method def Qmatrix(self): d = self.d h = 3*d-2 Bh = [vector(elt) for elt in monomials(h)] n = len(Bh) Q0 = Matrix(n, n) Q1 = Matrix(n, n) o = [vector(elt) for elt in [(d,0,0), (0,d,0), (0,0,d)]]; for t, bt in enumerate(Bh): i = min([i for i in range(3) if bt[i] >= d]) for u, bu in enumerate(Bh): a = vector(bu) - vector(bt) + o[i] if any(elt < 0 for elt in a): continue Q0[t, u] = self.F[a]*((h if i == 2 else d) - bt[i] - a[i]) Q1[t, u] = - self.F[a]*a[i] return Q0, Q1 def Qk(self, k): # we realign the k variable so that # newk = p - 2 - k, so that Qk(0) ... Qk(p-2) # matches with Q(p-2) ... Q(0) newk = -2 -k Q0, Q1 = self.Qmatrix() return Q0 + newk*Q1 def cm_harvey(self, p, Qp): d = self.d h = 3*d-2 z = vector([0,0,h-d]) Bh = [vector(elt) for elt in monomials(h)] # interior of Newton polygon Bds = [vector([1,1,2]), vector([2,1,1]), vector([1,2,1])] V = Matrix(len(Bh), len(Bds)) for i, bi in enumerate(Bh): for j, bj in enumerate(Bds): V[i, j] = invfactorialproduct(GF(p), (p*bj + z -bi)/d) return Qp*V
#F = random_quartic(10) R.<x0,x1,x2> = ZZ[] F = -60*x0^4 + 982*x0^3*x1 - 1003*x0^2*x1^2 + 918*x0*x1^3 - 265*x1^4 - 889*x0^3*x2 - 995*x0^2*x1*x2 + 116*x0*x1^2*x2 + 457*x1^3*x2 - 938*x0^2*x2^2 + 522*x0*x1*x2^2 - 912*x1^2*x2^2 + 471*x0*x2^3 + 828*x1*x2^3 - 760*x2^4; X = Quartic(F) expected = {7: 6, 13: 10, 17: 7, 23: 3, 29: 23, 31: 4, 37: 32, 41: 5, 43: 0, 47: 5, 59: 53, 61: 57, 67: 4, 71: 12, 73: 66, 79: 11, 89: 78, 97: 2, 101: 9, 103: 99, 107: 1, 109: 32, 113: 24, 127: 7, 131: 121, 137: 0, 139: 13, 149: 23, 151: 14, 163: 16, 167: 7, 173: 9, 179: 10, 181: 0, 191: 6, 193: 6, 197: 17, 199: 184} %time CM_naive = X.cartier_manin_matrices(200, algorithm="naive") %time CM_CHS = X.cartier_manin_matrices(200, algorithm="CHS") %time CM_H = X.cartier_manin_matrices(200, algorithm="Harvey") print(CM_naive == CM_CHS, CM_naive == CM_H) print({p: A.trace() for p, A in CM_CHS.items()} == {p:expected[p] for p in expected if p <= 200})
CPU times: user 682 ms, sys: 48 µs, total: 682 ms Wall time: 696 ms CPU times: user 481 ms, sys: 0 ns, total: 481 ms Wall time: 504 ms CPU times: user 1.49 s, sys: 4.99 ms, total: 1.49 s Wall time: 1.52 s True True True
%time CM_compressed = X.cartier_manin_matrices(2^10) %time CM_uncompressed = X.cartier_manin_matrices(2^10, compressed=False) %time CM_H = X.cartier_manin_matrices(2^10, algorithm="naive") print(CM_compressed == CM_uncompressed, CM_compressed == CM_H)
CPU times: user 3.57 s, sys: 26 ms, total: 3.59 s Wall time: 3.64 s CPU times: user 7.76 s, sys: 34.5 ms, total: 7.79 s Wall time: 8.17 s CPU times: user 14.7 s, sys: 17.1 ms, total: 14.7 s Wall time: 15.2 s True True
%time foo = X.cartier_manin_matrices(2^11) %time foo = X.cartier_manin_matrices(2^11, compressed=False)
CPU times: user 7.36 s, sys: 54.1 ms, total: 7.41 s Wall time: 7.52 s CPU times: user 19.6 s, sys: 113 ms, total: 19.7 s Wall time: 20.6 s
%time foo = X.cartier_manin_matrices(2^12)
CPU times: user 16.9 s, sys: 78.1 ms, total: 17 s Wall time: 17.6 s
%time foo = X.cartier_manin_matrices(2^13)
CPU times: user 38.8 s, sys: 445 ms, total: 39.2 s Wall time: 40.9 s
%time foo = X.cartier_manin_matrices(2^14)
CPU times: user 1min 26s, sys: 841 ms, total: 1min 27s Wall time: 1min 30s
%time foo = X.cartier_manin_matrices(2^14, compressed=False)
CPU times: user 4min 18s, sys: 1.41 s, total: 4min 20s Wall time: 4min 30s