{{{id=1|
def M(E,p,b):
    r"""
    INPUT:
    - 'E' - Hyperelliptic Curve of the form y^2 = f(x)
    - 'p' - a prime number
    - 'b' - q=p^b
    OUTPUT:
    - 'M' The matrix M = (c_(pi-j)), f(x)^((p-1)/2) = \sum c_i x^i
  
  Next step  - '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
    if not p.is_prime():
        return 'p must be prime'
    Fq=GF(p^b,'a');
    C=E.change_ring(Fq)
    f,h = C.hyperelliptic_polynomials()    
    if h != 0:
        return 'E must be of the form y^2 = f(x)'
    d = f.degree()
    if d%2 == 0:
        return 'the degree of f is even'
    df=f.derivative()
    R=df.resultant(f)
    if R == 0:
        return 'so curve is not smooth'    
    F = f^((p-1)/2)
    #coefficients returns a_0, ... , a_n where f(x) = a_n x^n + ... + a_0
    C = F.list()
    g = E.genus()
    M=[];
    for j in range(1,g+1):
        H=Sequence([(C[i]) for i in range((p*j-1-g), (p*j-1))])
        H.reverse()
        M.append(H);
    return matrix(Fq,M), g;
///
}}}

{{{id=9|
def a(g, M):
    a=g-rank(M)
    return a
///
}}}

{{{id=4|
P.<x>=QQ[]
///
}}}

{{{id=2|
E=HyperellipticCurve(x^9+1,0)
///
}}}

{{{id=3|
M(E,5,2)[1]
///
4
}}}

{{{id=5|
a(M(E,5,2)[1],M(E,5,2)[0])
///
2
}}}

{{{id=8|

///
}}}