{{{id=1|
for i in srange(1,4):
    p=next_prime(10^i)
    for j in srange(2,5):
        b=j
        Fq=GF(p^b,'a');
        R.<x>=Fq[];
        f = x^p - x
        te = walltime()
        d=exponentiating_f(p,b,f);
        t = walltime()
        d=f^((p-1)/2)
        print walltime(te),walltime(t)
///
0.0772469043732 0.000319957733154
0.0801980495453 0.000555992126465
0.0959780216217 0.000894784927368
7.04978108406 0.018355846405
10.6509320736 0.0304780006409
12.040555954 0.0518951416016
}}}

{{{id=2|
def exponentiating_f(p,b,f):

    """
    INPUT:
    - f - Function over F_p^b in x
    - p - Size of the base field
    - b - Degree of the extension of the base field
    OUTPUT:
    F = f^((p-1)/2)
    """
    Fq=GF(p^b,'a');
    R.<x>=Fq[];

    # calculate the frobenius
    frob = R.hom([x^p])
    fp = frob(f)
    # to compute f^(p-1) use the double slash divisor which does not include a remainder, as there is not one
    ftmp = fp//f
    # compute the square root
    prec = (f.degree() * (p-1)/2) +1
    K.<x> = PowerSeriesRing(Fq, 'x',prec)
    F = K(ftmp).square_root()
    # convert from a power series to a polynomial
    F = F.truncate()
    return F
///
}}}

{{{id=5|
f^((p-1)/2)
///
x^3 + 3*x^2 + 3*x + 1
}}}

{{{id=6|
p=next_prime(10^2)
b=3
Fq=GF(p^b,'a')
R.<x>=Fq[]
f = x^p - x
%prun d=exponentiating_f(p,b,f)
///
}}}

{{{id=7|

///
}}}