DeRhamComputation/crystalline_p2.ipynb

32 KiB
Raw Blame History

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 fraction_pol(g):
        RxX.<x, X> = PolynomialRing(QQ)
        g = RxX(g)
        result = 0
        for a in g.monomials():
            c = g.monomial_coefficient(a)
            c = c%(p^2)
            dX = a.degree(X)
            dx = a.degree(x)
            if dX%p == 0:
                result += c*x^(dX//p + dx)
            else:
                result += c*X^(dX + dx*p)
        return result
    
def teichmuller(f):
    Rx.<x> = PolynomialRing(QQ)
    RxX.<x, X> = PolynomialRing(QQ, 2)
    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])
    f1 = sum(coeffs[i]*RxX(x^(i)) for i in range(0, len(coeffs))) + p*RxX(ff.coordinates[1](x = X))
    f1 = fraction_pol(f1)
    #RXp.<Xp> = PolynomialRing(Integers(p^2))
    #f1 = RXp(f1)
    return f1
def basis_de_rham_degrees(f, m, p):
    r = f.degree()
    delta = GCD(r, m)
    Rx.<x> = PolynomialRing(QQ)
    Rxy.<x, y> = PolynomialRing(QQ, 2)
    Fxy = FractionField(Rxy)
    basis_holo = holomorphic_differentials_basis(f, m, p)
    basis = []
    for k in range(0, len(basis_holo)):
        basis += [(basis_holo[k], Rx(0))]

    ## non-holomorphic elts of H^1_dR
    t = len(basis)
    degrees0 = {}
    degrees1 = {}
    for j in range(1, m):
        for i in range(1, r):
            if (r*(m-j) - m*i >= delta): 
                s = Rx(m-j)*Rx(x)*Rx(f.derivative()) - Rx(m)*Rx(i)*f
                psi = Rx(cut(s, i))
                basis += [(Fxy(psi/y^j), Fxy(m*y^(m-j)/x^i))]
                degrees0[t] = (psi.degree(), j)
                degrees1[t] = (-i, m-j)
                t += 1
    return basis, degrees0, degrees1

def de_rham_basis(f, m, p):
    basis, degrees0, degrees1 = basis_de_rham_degrees(f, m, p)
    return basis

def degrees_de_rham0(f, m, p):
    basis, degrees0, degrees1 = basis_de_rham_degrees(f, m, p)
    return degrees0

def degrees_de_rham1(f, m, p):
    basis, degrees0, degrees1 = basis_de_rham_degrees(f, m, p)
    return degrees1


class superelliptic:
    
    def __init__(self, f, m, p):
        Rx.<x> = PolynomialRing(QQ)
        Rxy.<x, y> = PolynomialRing(QQ, 2)
        Fxy = FractionField(Rxy)
        self.polynomial = Rx(f)
        self.exponent = m
        self.characteristic = p
        
        r = Rx(f).degree()
        delta = GCD(r, m)
        self.degree_holo = degrees_holomorphic_differentials(f, m, p)
        self.degree_de_rham0 = degrees_de_rham0(f, m, p)
        self.degree_de_rham1 = degrees_de_rham1(f, m, p)
        
        holo_basis = holomorphic_differentials_basis(f, m, p)
        holo_basis_converted = []
        for a in holo_basis:
            holo_basis_converted += [superelliptic_form(self, a)]
        
        self.basis_holomorphic_differentials = holo_basis_converted
        

        dr_basis = de_rham_basis(f, m, p)
        dr_basis_converted = []
        for (a, b) in dr_basis:
            dr_basis_converted += [superelliptic_cech(self, superelliptic_form(self, a), superelliptic_function(self, b))]
        
        self.basis_de_rham = dr_basis_converted
        
    def __repr__(self):
        f = self.polynomial
        m = self.exponent
        p = self.characteristic
        return 'Superelliptic curve with the equation y^' + str(m) + ' = ' + str(f)+' over finite field with ' + str(p) + ' elements.'
    
    def is_smooth(self):
        f = self.polynomial
        if f.discriminant() == 0:
            return 0
        return 1
    
    def genus(self):
        r = self.polynomial.degree()
        m = self.exponent
        delta = GCD(r, m)
        return 1/2*((r-1)*(m-1) - delta + 1)
        
class superelliptic_function:
    def __init__(self, C, g):
        p = C.characteristic
        RxXy.<x, X, y> = PolynomialRing(QQ, 3)
        FxXy = FractionField(RXy)
        f = C.polynomial
        r = f.degree()
        m = C.exponent
        
        self.curve = C
        self.function = g
        
    def __repr__(self):
        return str(self.function)
    
    def jth_component(self, j):
        g = self.function
        C = self.curve
        p = C.characteristic
        RxX.<x,X> = PolynomialRing(QQ, 2)
        FxX.<x> = FractionField(RxX)
        FxXRy.<y> = PolynomialRing(FxX)
        g = FxXRy(g)
        return coff(g, j)
    
    def __add__(self, other):
        C = self.curve
        g1 = self.function
        g2 = other.function
        g = reduction(C, g1 + g2)
        return superelliptic_function(C, g)
    
    def __sub__(self, other):
        C = self.curve
        g1 = self.function
        g2 = other.function
        g = reduction(C, g1 - g2)
        return superelliptic_function(C, g)
    
    def __mul__(self, other):
        C = self.curve
        g1 = self.function
        g2 = other.function
        g = reduction(C, g1 * g2)
        return superelliptic_function(C, g)
    
    def __truediv__(self, other):
        C = self.curve
        g1 = self.function
        g2 = other.function
        return superelliptic_function(C, g1 / g2)
    
def diffn(self):
    C = self.curve
    f = C.polynomial
    m = C.exponent
    p = C.characteristic
    g = self.function
    RXy.<X, y> = PolynomialRing(QQ, 2)
    FXy = FractionField(RXy)
    g = RXy(g)
    A = g.derivative(X)*X^(-(p-1))/p
    t = teichmuller(f)
    B = g.derivative(y)*t.derivative()/(m*y^(m-1))*X^(-(p-1))/p
    return superelliptic_form(C, A+B)
        
class superelliptic_form:
    def __init__(self, C, g):
        p = C.characteristic
        Rxy.<x, y> = PolynomialRing(QQ, 2)
        Fxy = FractionField(Rxy)
        g = Fxy(reduction_form(C, g))
        self.form = g
        self.curve = C      
        
    def __add__(self, other):
        C = self.curve
        g1 = self.form
        g2 = other.form
        g = reduction(C, g1 + g2)
        return superelliptic_form(C, g)
    
    def __sub__(self, other):
        C = self.curve
        g1 = self.form
        g2 = other.form
        g = reduction(C, g1 - g2)
        return superelliptic_form(C, g)
    
    def __repr__(self):
        g = self.form
        if len(str(g)) == 1:
            return str(g) + ' dx'
        return '('+str(g) + ') dx'

    def __rmul__(self, constant):
        C = self.curve
        omega = self.form
        return superelliptic_form(C, constant*omega)        

    
    def coordinates(self):
        C = self.curve
        p = C.characteristic
        m = C.exponent
        Rx.<x> = PolynomialRing(QQ)
        Fx = FractionField(Rx)
        FxRy.<y> = PolynomialRing(Fx)
        g = C.genus()
        degrees_holo = C.degree_holo
        degrees_holo_inv = {b:a for a, b in degrees_holo.items()}
        basis = C.basis_holomorphic_differentials
        
        for j in range(1, m):
            omega_j = Fx(self.jth_component(j))
            if omega_j != Fx(0):
                d = degree_of_rational_fctn(omega_j, p)
                index = degrees_holo_inv[(d, j)]
                a = coeff_of_rational_fctn(omega_j, p)
                a1 = coeff_of_rational_fctn(basis[index].jth_component(j), p)
                elt = self - (a/a1)*basis[index]
                return elt.coordinates() + a/a1*vector([QQ(i == index) for i in range(0, g)])
            
        return vector(g*[0])
    
    def jth_component(self, j):
        g = self.form
        C = self.curve
        p = C.characteristic
        Rx.<x> = PolynomialRing(QQ)
        Fx = FractionField(Rx)
        FxRy.<y> = PolynomialRing(Fx)
        Fxy = FractionField(FxRy)
        Ryinv.<y_inv> = PolynomialRing(Fx)
        g = Fxy(g)
        g = g(y = 1/y_inv)
        g = Ryinv(g)
        return coff(g, j)
    
    def is_regular_on_U0(self):
        C = self.curve
        p = C.characteristic
        m = C.exponent
        Rx.<x> = PolynomialRing(QQ)
        for j in range(1, m):
            if self.jth_component(j) not in Rx:
                return 0
        return 1
    
    def is_regular_on_Uinfty(self):
        C = self.curve
        p = C.characteristic
        m = C.exponent
        f = C.polynomial
        r = f.degree()
        delta = GCD(m, r)
        M = m/delta
        R = r/delta
        
        for j in range(1, m):
            A = self.jth_component(j)
            d = degree_of_rational_fctn(A, p)
            if(-d*M + j*R -(M+1)<0):
                return 0
        return 1
    
    
class superelliptic_cech:
    def __init__(self, C, omega, fct):
        self.omega0 = omega
        self.omega8 = omega - diffn(fct)
        self.f = fct
        self.curve = C
        
    def __add__(self, other):
        C = self.curve
        return superelliptic_cech(C, self.omega0 + other.omega0, self.f + other.f)
    
    def __sub__(self, other):
        C = self.curve
        return superelliptic_cech(C, self.omega0 - other.omega0, self.f - other.f)

    def __rmul__(self, constant):
        C = self.curve
        w1 = self.omega0.form
        f1 = self.f.function
        w2 = superelliptic_form(C, constant*w1)
        f2 = superelliptic_function(C, constant*f1)
        return superelliptic_cech(C, w2, f2)    
    
    def __repr__(self):
        return "(" + str(self.omega0) + ", " + str(self.f) + ", " + str(self.omega8) + ")" 

    def coordinates(self):
        C = self.curve
        p = C.characteristic
        m = C.exponent
        Rx.<x> = PolynomialRing(QQ)
        Fx = FractionField(Rx)
        FxRy.<y> = PolynomialRing(Fx)
        g = C.genus()
        degrees_holo = C.degree_holo
        degrees_holo_inv = {b:a for a, b in degrees_holo.items()}
        degrees0 = C.degree_de_rham0
        degrees0_inv = {b:a for a, b in degrees0.items()}
        degrees1 = C.degree_de_rham1
        degrees1_inv = {b:a for a, b in degrees1.items()}
        basis = C.basis_de_rham
        
        omega = self.omega0
        fct = self.f
        
        if fct.function == Rx(0) and omega.form != Rx(0):
            for j in range(1, m):
                omega_j = Fx(omega.jth_component(j))
                if omega_j != Fx(0):
                    d = degree_of_rational_fctn(omega_j, p)
                    index = degrees_holo_inv[(d, j)]
                    a = coeff_of_rational_fctn(omega_j, p)
                    a1 = coeff_of_rational_fctn(basis[index].omega0.jth_component(j), p)
                    elt = self - (a/a1)*basis[index]
                    return elt.coordinates() + a/a1*vector([QQ(i == index) for i in range(0, 2*g)])
                    
        for j in range(1, m):
            fct_j = Fx(fct.jth_component(j))
            if (fct_j != Rx(0)):
                d = degree_of_rational_fctn(fct_j, p)
            
                if (d, j) in degrees1.values():
                    index = degrees1_inv[(d, j)]
                    a = coeff_of_rational_fctn(fct_j, p)
                    elt = self - (a/m)*basis[index]
                    return elt.coordinates() + a/m*vector([QQ(i == index) for i in range(0, 2*g)])
                
                if d<0:
                    a = coeff_of_rational_fctn(fct_j, p)
                    h = superelliptic_function(C, FxRy(a*y^j*x^d))
                    elt = superelliptic_cech(C, self.omega0, self.f - h)
                    return elt.coordinates()
            
                if (fct_j != Rx(0)):
                    G = superelliptic_function(C, y^j*x^d)
                    a = coeff_of_rational_fctn(fct_j, p)
                    elt =self - a*superelliptic_cech(C, diffn(G), G)
                    return elt.coordinates()

        return vector(2*g*[0])
    
    def is_cocycle(self):
        w0 = self.omega0
        w8 = self.omega8
        fct = self.f
        if not w0.is_regular_on_U0() and not w8.is_regular_on_Uinfty():
            return('w0 & w8')
        if not w0.is_regular_on_U0():
            return('w0')
        if not w8.is_regular_on_Uinfty():
            return('w8')
        if w0.is_regular_on_U0() and w8.is_regular_on_Uinfty():
            return 1
        return 0
        
def degree_of_rational_fctn(f, p):
    Rx.<x> = PolynomialRing(QQ)
    Fx = FractionField(Rx)
    f = Fx(f)
    f1 = f.numerator()
    f2 = f.denominator()
    d1 = f1.degree()
    d2 = f2.degree()
    return(d1 - d2)

def coeff_of_rational_fctn(f, p):
    Rx.<x> = PolynomialRing(QQ)
    Fx = FractionField(Rx)
    f = Fx(f)
    if f == Rx(0):
        return 0
    f1 = f.numerator()
    f2 = f.denominator()
    d1 = f1.degree()
    d2 = f2.degree()
    a1 = f1.coefficients(sparse = false)[d1]
    a2 = f2.coefficients(sparse = false)[d2]
    return(a1/a2)

def coff(f, d):
    lista = f.coefficients(sparse = false)
    if len(lista) <= d:
        return 0
    return lista[d]

def cut(f, i):
    R = f.parent()
    coeff = f.coefficients(sparse = false)
    return sum(R(x^(j-i-1)) * coeff[j] for j in range(i+1, f.degree() + 1))

def polynomial_part(p, h):
    Rx.<x> = PolynomialRing(QQ)
    h = Rx(h)
    result = Rx(0)
    for i in range(0, h.degree()+1):
        if (i%p) == p-1:
            power = Integer((i-(p-1))/p)
            result += Integer(h[i]) * x^(power)    
    return result
Rx.<x> = PolynomialRing(QQ)
f = Rx(x^3 - x)
teichmuller(f)
-720*X^9 - 15312*X^7 - 107160*X^5 + x^3 - 249984*X^3 - x
6*X^7 + 3*X^5 + x^3 + 8*x
Rx.<x> = PolynomialRing(QQ)
C = superelliptic(x^3 - x, 2, 3)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/tmp/ipykernel_1111/3447231159.py in <module>
      1 Rx = PolynomialRing(QQ, names=('x',)); (x,) = Rx._first_ngens(1)
----> 2 C = superelliptic(x**Integer(3) - x, Integer(2), Integer(3))

/tmp/ipykernel_1111/3436947063.py in __init__(self, f, m, p)
     68         dr_basis_converted = []
     69         for (a, b) in dr_basis:
---> 70             dr_basis_converted += [superelliptic_cech(self, superelliptic_form(self, a), superelliptic_function(self, b))]
     71 
     72         self.basis_de_rham = dr_basis_converted

/tmp/ipykernel_1111/3436947063.py in __init__(self, C, omega, fct)
    260     def __init__(self, C, omega, fct):
    261         self.omega0 = omega
--> 262         self.omega8 = omega - diffn(fct)
    263         self.f = fct
    264         self.curve = C

/tmp/ipykernel_1111/3436947063.py in diffn(self)
    153     A = g.derivative(Xp)*Xp**(-(p-Integer(1)))/p
    154     t = teichmuller(f)
--> 155     B = g.derivative(y)*t.derivative()/(m*y**(m-Integer(1)))*Xp**(-(p-Integer(1)))/p
    156     return superelliptic_form(C, A+B)
    157 

/ext/sage/9.5/local/var/lib/sage/venv-python3.9.9/lib/python3.9/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__mul__ (build/cythonized/sage/structure/element.c:12253)()
   1514             return (<Element>left)._mul_(right)
   1515         if BOTH_ARE_ELEMENT(cl):
-> 1516             return coercion_model.bin_op(left, right, mul)
   1517 
   1518         cdef long value

/ext/sage/9.5/local/var/lib/sage/venv-python3.9.9/lib/python3.9/site-packages/sage/structure/coerce.pyx in sage.structure.coerce.CoercionModel.bin_op (build/cythonized/sage/structure/coerce.c:11751)()
   1246         # We should really include the underlying error.
   1247         # This causes so much headache.
-> 1248         raise bin_op_exception(op, x, y)
   1249 
   1250     cpdef canonical_coercion(self, x, y):

TypeError: unsupported operand parent(s) for *: 'Multivariate Polynomial Ring in Xp, y over Rational Field' and 'Univariate Polynomial Ring in Xp over Ring of integers modulo 9'
C.basis_de_rham
[((1/y) dx, 0, (1/y) dx), ((x/y) dx, 2/x*y, ((-1)/(x*y)) dx)]
g = basis_de_rham_degrees(x^3 - x, 2, 3)
RXpy.<Xp, y> = PolynomialRing(QQ, 2)
FXpy = FractionField(RXpy)
g(x = Xp^p, y = y)
2*y/Xp^3
t = teichmuller(x^3 - x)
t.derivative()
6*Xp^6 + 6*Xp^4 + 6*Xp^2
RxX.<x, X> = PolynomialRing(QQ)
g = RxX(x*X + 2*x*X^3)
g.monomials()[0].degree(X)
3
type(g.monomial_coefficient(X))
<class 'sage.rings.rational.Rational'>