### CoCalc provides the best real-time collaborative environment for Jupyter Notebooks, LaTeX documents, and SageMath, scalable from individual use to large groups and classes! Also, H100 GPUs starting at \$2/hour.

Views: 521
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)
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);

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
# 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()

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
```
```