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

%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)

%time foo = X.cartier_manin_matrices(2^11)
%time foo = X.cartier_manin_matrices(2^11, compressed=False)

%time foo = X.cartier_manin_matrices(2^12)

%time foo = X.cartier_manin_matrices(2^13)

%time foo = X.cartier_manin_matrices(2^14)

%time foo = X.cartier_manin_matrices(2^14,  compressed=False)