Source Code

Below is the implementation in SAGE [SJ05]:

R.<y,z,x> = PolynomialRing(QQ,3)
S.<y,z,x> = R.quotient([x^3 -x+1/4 - y^2, y*z - 1])
Q = x^3 - x + S(1/4)
Q5 = (x^5)^3 - (x^5) + S(1/4)
Qprime = 3*x^2 - S(1)
G = S(1/5) * (Q5 - Q^5)

M = 3
p = 5

def F(i):
    a = 0
    for k in range(M):
        b = binomial(-1/2, k) * p^(k+1) * G^k * x^(p*(i+1)-1) * z^((2*k+1)*p-1)
        a += b
    return a

V = VectorSpace(QQ,30)

def coefficient_vector(f):
    v = V(0)
    for e, c in f.dict().iteritems():
        v[3*e[1] + e[2]] = c
    return v

def d(i,j):
   i = i
   j = j
   if i == 0:
       return S(1/2 * -j) * Qprime * z^(1+j)
   return S(1/2) * (2*i*x^(i-1)*z^(j-1)-  j*x^i*Qprime*z^(j+1))

   def Wspace():
    dA = [coefficient_vector(reduce_mod(d(i, j),5)) \
                              for i in range(3) for j in range(1,5 +1) if j%2!=0]
    H1 = [coefficient_vector(reduce_mod(S(1),5)), coefficient_vector(reduce_mod(x,5))]
    B = V.subspace(dA).basis()
    print len(dA), len(B)
    print B
    print H1
    return V.subspace_with_basis(B + H1)

def reduce_mod(f, n):
    """
    Return the reduction of f modulo n.
    """
    R.<y,z,x> = PolynomialRing(IntegerModRing(n),3)
    return R(sage_eval(str(f),  {'x':x, 'y':y, 'z':z}))

def rd(i,j,n):
    return reduce_mod(d(i,j),n)

def coeff_list_of_pow(f, n, var_z_number=1, var_x_number=2, d=2):
    """
    if f is a polynomial in n variables and z is the var_number-th
    generator of the parent of f, this returns the degree-d list of
    coefficients of z^n up to degree d.
    """
    z = f.parent().gen(var_z_number)
    c = f.coefficient(z^n)
    x = f.parent().gen(var_x_number)
    return [c.coefficient(x^i) for i in range(d+1)]

def lowest_d_matrix(j, n, p=5):
    """
    Return the 3x3 matrix mod p of coefficients of the coefficients
    of z^n in the d(i,j) for i=0,1,2.  Here j > 0.
    """
    MS = MatrixSpace(GF(p), 3)
    v = []
    for i in range(3):
        v.append(coeff_list_of_pow(rd(i,j,p), n))
    return MS(v).transpose()

def find_lincomb_needed_to_subtract_off(f, n, p=5):
    j = n-1
    A = lowest_d_matrix(j,n,p)
    v = coeff_list_of_pow(f, n)
    w = Vector(GF(p), v)
    return A^(-1)*w

f0 = reduce_mod(F(0),25)
find_lincomb_needed_to_subtract_off(reduce_mod((1/5)*F(0),5), 4)
f01 = f0 - 5*rd(0,3,25)-10*rd(1,3,25)-20*rd(2,3,25)
find_lincomb_needed_to_subtract_off(reduce_mod((1/5)*(10*z^2+5*z^2*x),5),2)
f02 = f01 - 5*(2*rd(0,1,25) + 1*rd(1,1,25) + 2*rd(2,1,25))

f1 = reduce_mod(F(1),25)
f11 = f1-Mod(1/3,25)*reduce_mod(d(4,1),25)
find_lincomb_needed_to_subtract_off(reduce_mod((1/5)*(10*z^4+20*z^4*x+15*z^4*x^2),5),4)
f12 = f11 - 10*rd(0,3,25)-15*rd(1,3,25)-5*rd(2,3,25)
f13 = f12 - 20*rd(0,1,25)-23*rd(1,1,25)-13*rd(2,1,25)