DeRhamComputation/crystalline_p2.ipynb

N = 2
p = 3
RQ = PolynomialRing(QQ, 'X', 2*N)
X = RQ.gens()[:N]
Y = RQ.gens()[N:]
Rpx.<x> = PolynomialRing(GF(p), 1)
#RQx.<x> = PolynomialRing(QQ, 1)

def witt_pol(lista):
    n = len(lista)
    return sum(p^i*lista[i]^(p^(n-i-1)) for i in range(0, n))

def witt_sum(n):
    if n == 0:
        return X[0] + Y[0]
    return 1/p^n*(witt_pol(X[:n+1]) + witt_pol(Y[:n+1]) - sum(p^k*witt_sum(k)^(p^(n-k)) for k in range(0, n)))

class witt:
    def __init__(self, coordinates):
        self.coordinates = coordinates
    def __repr__(self):
        lista = [Rpx(a) for a in self.coordinates]
        return str(lista)
    def __add__(self, other):
        lista = []
        for i in range(0, N):
            lista+= [witt_sum(i)(self.coordinates + other.coordinates)]
        return witt(lista)
    def __rmul__(self, constant):
        if constant<0:
            m = (-constant)*(p^N-1)
            return m*self
        if constant == 0:
            return witt(N*[0])
        return self + (constant - 1)*self
    def __sub__(self, other):
        return self + (-1)*other

def teichmuller(f):
    Rx.<x> = PolynomialRing(QQ)
    RXp.<Xp> = PolynomialRing(QQ)
    f = Rx(f)
    ff = witt([f, 0])
    coeffs = f.coefficients(sparse=false)
    for i in range(0, len(coeffs)):
        ff -= coeffs[i]*witt([Rx(x^i), 0])
    print(ff)
    f1 = sum(coeffs[i]*RXp(Xp^(3*i)) for i in range(0, len(coeffs))) + p*RXp(ff.coordinates[1](x = Xp))
    RXp.<Xp> = PolynomialRing(Integers(p^2))
    f1 = RXp(f1)
    return f1

Rx.<x> = PolynomialRing(QQ)
f = Rx(x^3 - x)
teichmuller(f)

[0, -x^7 + x^5]

Xp^9 + 6*Xp^7 + 3*Xp^5 + 8*Xp^3

f = Rx(x^3 - x)

f.coefficients(sparse=false)

[0, -1, 0, 1]