Final Python Code for N and M
system:sage


{{{id=1|
%python
def Cartier_matrix(E,p,b):
    r"""
    INPUT:
    - 'E' - Hyperelliptic Curve of the form y^2 = f(x), defined over Q[x]
    E will be redefined over the finite field of order q
    - 'p' - a prime number, q = p^b
    - 'b' - q = p^b
    OUTPUT:
    - 'Cartier-Manin matrix' - The matrix M = (c_(pi-j)), f(x)^((p-1)/2) = \sum c_i x^i
    """
    
    #check primality and form of elliptic curve
    if not p.is_prime():
        return 'p must be prime'
    if p == 2:
        return 'p must be odd'  
          
    #change elliptic curve to be defined over Fq

    Fq = FiniteField(p**b,name='a')
    C = E.change_ring(Fq)
    g = C.genus()
        
    f,h = C.hyperelliptic_polynomials()
    
    if h != 0:
        return 'E must be of the form y^2 = f(x)'
    
    #check degree
    d = f.degree()
    if d%2 == 0:
        return 'error: the degree of f is even'
        
    #check smoothness
    df = f.derivative()
    R = df.resultant(f)
    if R == 0:
        return 'error: f(x) is not smooth'
        
    F = f**((p-1)/2)
    
    
    #coeff: creates a list of the coefficients of F: a_0, ... , a_n 
    #where F = a_n x^n + ... + a_0
    coeff = F.list()
    
    #inserting zeros when necessary-- that is, when deg(F) < p*g-1, 
    #which is the highest powered coefficient needed for our matrix
    zeros = [0 for i in range(p*g-len(coeff))]
    coeff = coeff + zeros
     
    M = [];
    for i in range(1,g+1):
        H = [(coeff[j]) for j in range( p*i-1, p*i-g-1, -1)]
        M.append(H);
         
    return M
///
}}}

{{{id=17|
%python
def Hasse_Witt(E,p,b):
    r"""
    INPUT:
    - 'E' - Hyperelliptic Curve of the form y^2 = f(x), defined over Q[x]
    - 'p' - a prime number, q = p^b; E will be redefined over the finite field of order q
    - 'b' - q = p^b
    OUTPUT:
    - 'Hasse-Witt matrix' -
    """
    
    #check primality and form of elliptic curve
    if not p.is_prime():
        return 'p must be prime'
    if p == 2:
        return 'p must be odd'  
          
    #change elliptic curve to be defined over Fq

    Fq = FiniteField(p**b,name='a')
    C = E.change_ring(Fq)
    g = C.genus()
        
    f,h = C.hyperelliptic_polynomials()
    
    if h != 0:
        return 'E must be of the form y^2 = f(x)'
    
    #check degree
    d = f.degree()
    if d%2 == 0:
        return 'error: the degree of f is even'
        
    #check smoothness
    df = f.derivative()
    R = df.resultant(f)
    if R == 0:
        return 'error: f(x) is not smooth'
        
    F = f**((p-1)/2)
    
    
    #coeff: creates a list of the coefficients of F: a_0, ... , a_n 
    #where F = a_n x^n + ... + a_0
    coeff = F.list()
    
    #inserting zeros when necessary-- that is, when deg(F) < p*g-1, 
    #which is the highest powered coefficient needed for our matrix
    zeros = [0 for i in range(p*g-len(coeff))]
    coeff = coeff + zeros
    
    Mall = []
    for k in range(g):  
        Mk = [];
        a = p**k
        for i in range(1,g+1):
            H = [(coeff[j])**a for j in range( p*i-1, p*i-g-1, -1)]
            Mk.append(H);
        Mall.append(matrix(Fq,Mk));
         
    #multiply all of the Mk matrices together to get the matrix N
    N = identity_matrix(Fq,g)
    for l in Mall:
         N = N*l;
    return N
///
}}}

{{{id=2|
P.<x>=QQ[];
E = HyperellipticCurve(x^11+x+1,0);
E2 = HyperellipticCurve(x^23-x,0);
E3 = HyperellipticCurve(x^51-x,0);
///
}}}

{{{id=24|
Cartier_matrix(E,13,1)
///
([1, 6, 2, 7, 2, 6, 1, 0, 0, 0, 0, 6, 4, 8, 8, 4, 6, 0, 0, 0, 0, 0, 2, 8, 12, 8, 2, 0, 0, 0, 0, 0, 0, 7, 8, 8, 7, 0, 0, 0, 0, 0, 0, 0, 2, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], [[4, 6, 0, 0, 0], [8, 12, 8, 2, 0], [0, 0, 7, 8, 8], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]])
}}}

{{{id=25|
N(E,13,1)
///
(
[ 4  6  0  0  0]  [12  5  7  6  7]
[ 8 12  8  2  0]  [11 10  8  7  9]
[ 0  0  7  8  8]  [ 0  0 11  7  7]
[ 0  0  0  0  0]  [ 0  0  0  0  0]
[ 0  0  0  0  0], [ 0  0  0  0  0]
)
}}}

{{{id=16|
%hide
%python
def Hasse_Vitt(E,p,b):
    r"""
    INPUT:
    - 'M'- Cartier matrix.
    - 'g' -genus
    - 'p' -over the field of char p
    - 'b' - power of p, Fq,  
    
    OUTPUT:
    - 'N' - The matrix N = M M^(p)...M^(p^(g-1)) where M = (c_(pi-j)), f(x)^((p-1)/2) = \sum c_i x^i
    """
    
    #Call the Cartier_matrix command to compute M
    coeff,M = Cartier_matrix(E,p,b)
    print coeff
    print M
    
    g = E.genus();
    
    from sage.matrix.constructor import matrix
    
    Fq = FiniteField(p**b,name='c');
    
    def frob_mat(Coeffs, k):
        a = p**k
        mat = []
        Coeffs_pow = [c**a for c in Coeffs]
        for i in range(1,g+1):
            H = [(Coeffs_pow[j]) for j in range(p*i-1, p*i-g-1, -1)]
            mat.append(H);
        return matrix(Fq,mat)
        
    Mall = [M] + [frob_mat(coeff,k) for k in range(1,g)]
    
    #multiply to create the identity matrix    
    N = identity_matrix(Fq,g)
    for l in Mall:
         N = N*l;
    return N
///
}}}

{{{id=19|
N(E,7,1)
///
[0 6 0 0 6]
[4 1 4 3 4]
[0 0 0 0 0]
[2 3 2 5 0]
[3 6 3 3 0]
}}}

{{{id=20|
timeit('N(E,103,1)')
///
125 loops, best of 3: 5.11 ms per loop
}}}

{{{id=21|
timeit('prime_range(10^6)')
///
25 loops, best of 3: 10.3 ms per loop
}}}

{{{id=22|
N(E,103,1)
///
[52 36 16 12 31]
[48 66 47 33 73]
[77 99 89 14  0]
[97 52 44 74 29]
[62 73 59 58 98]
}}}

{{{id=23|

///
}}}