{{{id=1|
def N(E,p,b):
    r"""
    INPUT:
    - 'E' - Hyperelliptic Curve of the form y^2 = f(x)
    - 'p' - a prime number
    - 'b' - q=p^b
    
    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
    """
    
    #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 = GF(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()
    
    #solving weird indexing issue (inserting zeros when necessary-- that is, when p < 2g-1)
    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 the Mk matrices together to get 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=4|
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=9|
N(E2,5,1)
///
[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 2 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]
}}}

{{{id=12|
EBig = HyperellipticCurve(x^101-x,0)
///
}}}

{{{id=13|
N(EBig,7,1)
///
50 x 50 dense matrix over Finite Field of size 7
}}}