"""
Broumas (Compositio 107, 1997) gave an explicit formula for the
p-descent map.  I wanted to implement it in Sage to answer a question
described below.

Context:  
K is a function field of characteristic p, E an elliptic curve over K.
Assume j(E) is not a pth-power in K.  (Equivalently, E does not have a
point of order p defined over a separable extension of K.)  General
theory gives an injective homomorphism

	E(K)/pE(K) -> K

(the target is the additive group of K).  It actually lands in a
certain subspace, but that will not be needed here.  Broumas gives
an explicit formula taking the (x,y)-coordinates of a point (on a
special model of E) to K.  It involves field operations and a
derivative operator.  When p>3, "a special model" just means the usual
Weierstrass form y^2=x^3+ax+b.

My problem:
p>2.  Consider y^2=x(x+1)(x+t) over Fp(t).
Fix f>0 and let q=p^f, d=q+1.
Let K=F_{q^2}(u) and make it an extension of Fp(t) by setting t=u^d.
I have a bunch of points on E(K):  Let P(u)=(u,u(u+1)^((q+1)/2)), let
zeta be a primitive d-th root of unity in K, and set

     P_i=P(zeta^i u)  i=0,...,d-1

These points are almost independent (sum of P_i over even i is
torsion, simialrly for summing over odd i) and generate a subgroup of
E(K) of rank d-2=q-1.  There are also a few explicit torsion points
around, which we will ignore below because they are of order prime to
p.  Let V be the subgroup of E(K) generated by the P_i (and the
torsion points).  Then it turns out that the index [E(K):V] is finite
of p-power order.  My problem is to compute in specific cases whether
the index is 1 or not.  This amounts to computing the dimension as a
vector space over Fp of the image of V under the descent map:

      V/pV -> E(K)/pE(K) -> K

If the image has dimension q-1 then the index is 1 and we have the
full Mordell-Weil group explicitly.  If it is less, we need to find
more points.

Based on data produced by Sage and other considerations, I guess that
the index is >1 exactly when f>2.

The code below is divided into two cases, p=3 and p>3.
"""
#################################
##p=3:
p=3; f=3; q=3^f
k=FiniteField(q^2,name='a')
zeta=k.multiplicative_generator()^(q-1)  # a primitive (q+1)st root of unity
R, u = PolynomialRing(k, 'u').objgen()
K=R.fraction_field()
t=u^(q+1)
EL=EllipticCurve(K,[0,t+1,0,t,0])
##
## special weierstrass form suitable for p=3 descent map 
## (a1=a3=a4=0)
E=EL.change_weierstrass_model([1/(t+1),t/(t+1),0,0])
## so (x,y)=(x'/(t+1)^2+t/(t+1),y'/(t+1)^3)
## and (x',y')=((t+1)^2*x-t*(t+1),(t+1)^3*y)
##
a2=E.a2(); a6=E.a6(); A=a2 # Hasse invariant
def wp(z):
    return z^p-A*z
##
## Broumas' derivative operator (script D in the paper):
def DD(z):
    return (t+1)^2*(t-1)*u*derivative(z,u)
## DD characterized by DD(a6)/a6==a2
##
## Broumas' descent map:
## The formula corrects a typo in Broumas 
##   (missing x in last term of formula for mu)
def mu(P):
    if 2*P==E(0):
        return 0
    else:
        x=P[0]/P[2]; y=P[1]/P[2]
        return wp(DD(x)/(2*y))+y+wp((x^2+(-DD(a2)/a2+DD(a6)/a6)*x)/(2*y))
##
#### Torsion points
P0=E([-t*(t+1),0]); P1=E([t^2-1,0]); Pt=E([-t^3+t,0])  
Q4=E([u^((q+1)/2)*(t+1)^2-t^2-t,(t+1)^3*u^((q+1)/2)*(u^((q+1)/2)+1)]) # 2*Q4==P0
##
## Generators
def point(i):
    return E([(t+1)^2*zeta^i*u-t*(t+1),(t+1)^3*zeta^i*u*(zeta^i*u+1)^((q+1)/2)])
P=[]
for i in range(q+1):
    P.append(point(i))
print 'points created'    
## checks when q=3:    
## P[0]+P[2]==Q4-P1
## P[1]+P[3]==Q4
##
## Descents:
L=[]
for i in range(q+1):
    L.append(mu(P[i]))
print 'descents done'    
##
## the L[i] (miraculously) turn out to be polynomials
## strip out the coefficients and turn them into vectors with Fp entries:
N=L[0].numerator().degree()+1
m=matrix(GF(p),q+1,2*f*N)
for pt in range(q+1):
    x=L[pt]; print 'row', pt, 'created'
    y=x.numerator().coeffs()
    for i in range(N):
        z=y[i]._vector_()
        for j in range(2*f):
            m[pt,2*f*i+j]=z[j]
##
## rank of m is the dimension of the image of the points under the descent map
(q-1)-rank(m)


#################################################
##p>3:
p=5; f=3; q=p^f
k=FiniteField(q^2,name='a')
zeta=k.multiplicative_generator()^(q-1)  # a primitive (q+1)st root of unity
R, u = PolynomialRing(k, 'u').objgen()
K=R.fraction_field()
t=u^(q+1)
EL=EllipticCurve(K,[0,t+1,0,t,0])
E=EL.change_weierstrass_model([1,-(t+1)/3,0,0])
##
## Basic invariants of E:
a4=E.a4(); a6=E.a6(); delta=E.discriminant(); j=E.j_invariant(); j_prime=derivative(j,u)
##
## Hasse invariant, Weierstrass P, and M function:
S, zz = PolynomialRing(K, 'zz').objgen()
g=(zz^3+E.a4()*zz+E.a6())^((p-1)/2)
A=g.coeffs()[p-1]
M=g.quo_rem(zz^p)[0]
def wp(z):
    return z^p-A*z
##
## Broumas' derivative operator (script D in the paper):
def DD(z):
    return 18*(a6/a4)*(j/j_prime)*derivative(z,u)
##
## Broumas' descent map
def mu(P):
    if 2*P==E(0):
        return 0
    else:
        x=P[0]/P[2]; y=P[1]/P[2]
        return wp(DD(x)/(2*y))+y*M(x)+wp((-2*x^2-(1/6)*(DD(delta)/delta)*x-(4/3)*a4)/(2*y))
##
## Torsion points:
P0=E([(t+1)/3,0]); P1=E([-1+(t+1)/3,0]); Pt=E([-t+(t+1)/3,0])  
Q4=E([u^((q+1)/2)+(t+1)/3,u^((q+1)/2)*(u^((q+1)/2)+1)]) # a point of order 4
##
## Generators:
def point(i):
    return E([zeta^i*u+(t+1)/3,zeta^i*u*(zeta^i*u+1)^((q+1)/2)])
P=[]
for i in range(q+1):
    P.append(point(i))
print 'points created'    
##
## checks when q=7:
## P[0]+P[2]+P[4]+P[6] == Q4+P1
## P[1]+P[3]+P[5]+P[7]==Q4
##
## Descents:
L=[]
for i in range(q+1):
    L.append(mu(P[i]))
print 'descents done'    
##
## the L[i] (miraculously) turn out to be polynomials
## strip out the coefficients and turn them into vectors with Fp entries:
N=L[0].numerator().degree()+1
m=matrix(GF(p),q+1,2*f*N)
for pt in range(q+1):
    x=L[pt]; print 'row', pt, 'created'
    y=x.numerator().coeffs()
    for i in range(N):
        z=y[i]._vector_()
        for j in range(2*f):
            m[pt,2*f*i+j]=z[j]
##
## rank of m is the dimension of the image of the points under the descent map
(q-1)-rank(m)