{{{id=1|
def weird_matrix_thing(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');
    E.change_ring(Fq)
    f,h = E.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'
   
    g=derivative(f,x)
    R=g.resultant(f)
    if R == 0:
        return 'f(x) has repeated roots, so it 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(QQ,M)
///
}}}

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

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

{{{id=22|
timeit('weird_matrix_thing(E,13,2)')
///
125 loops, best of 3: 4.84 ms per loop
}}}

{{{id=26|
p=11;b=2;
///
}}}

{{{id=25|
Fq=GF(p^b,'a');
E.change_ring(Fq)
///
Hyperelliptic Curve over Finite Field in a of size 11^2 defined by y^2 = x^9 + x + 1
}}}

{{{id=27|

///
}}}