# load auxiliary functions
load mchinclude2.sage
#
#
#
# THE MAIN ROUTINE FOLLOWS
#
#E is an elliptic curve over the rationals, N is a positive integer with conductor(E) dividing N,
#gIndex is the nonnegative integer corresponding to the index of g in the cuspidal space of modular symbols,
#and n is the number of fourier coefficients taken in the computations.

def chowheegner(E,N,gIndex,n):

#
#
#We start by introducing the basic ingredients that can be computed by means of the above functions:
#the rational function u=etagen(N), the list of generators of Gamma_0(N), the lattice of E, etc.
#
#
    R = LaurentSeriesRing(CC, 'q')
    q=R.0
    Q = LaurentSeriesRing(QQ, 't')
    t=Q.0
    L=E.period_lattice()
    f=E.newform().q_expansion(n)
    cond=E.conductor()
    if cond not in divisors(N):
        return 'The conductor of E must divide N' #really, we want this to raise an error, but I don't know how to program that in
    div=len(divisors(N/cond))
#
#
#
    print 'Computing etaproduct...'
    u=etagen(N).q_expansion(n)
    MG=ModularSymbols(N)
    M=MG.cuspidal_subspace()
    genus=M.dimension()/2
    Gens,Mat=hypgens(N)
    Mn = MatrixSpace(QQ,2*genus)
    MatInv=Mat.solve_right(Mn.identity_matrix())
#
# 
#Here we start computing the basis of the de Rham space. First we need a complete q-expansion basis, which we
#find by taking a basis of modular forms and multiplying by a power of the eta product found above. We check
#to determine if this forms a basis for H^1_dR(X) by computing the matrix for the Poincare pairing.
#
#
    print 'Computing basis for H1dR and the matrix for the Poincare pairing...'
    holomorphicforms = [Q(x) for x in CuspForms(N).q_expansion_basis(n)]
    poleord=-Q(u).valuation()
    zeroord=holomorphicforms[len(holomorphicforms)-1].valuation()
    uexp=ceil(QQ(zeroord)/QQ(poleord))
    meromorphicforms = [Q(u)^(uexp)*x for x in holomorphicforms]
    for i in [0..len(M)-1]:
        if M[i].dimension()==2*div:
            if M[i].simple_factors()[0].q_eigenform(2*genus,'a') == f.truncate(2*genus):
                fIndex = i
    H1dR = holomorphicforms+meromorphicforms
    minpolord=max([-x.exponents()[0] for x in meromorphicforms])
    H1dRtrunc=[x.truncate(minpolord+2) for x in H1dR]
    pdmat = pd_matrix(H1dRtrunc)
    while det(pdmat)==0:
        uexp=uexp+1
	meromorphicforms = [Q(u)^(uexp)*x for x in holomorphicforms]
	H1dR=holomorphicforms+meromorphicforms
	minpolord=max([-x.exponents()[0] for x in meromorphicforms])
	H1dRtrunc=[x.truncate(minpolord+2) for x in H1dR]
    	pdmat = pd_matrix(H1dRtrunc)
    homology=[]
    for i in [0..len(M)-1]:
    	homology+=[x.list() for x in M[i].integral_basis()]
    homologymatrix=Matrix(homology)
#
#The echelon basis will be used for computing the alpha corrections later.
#
#
#Next we need bases of the f- and g-isotypic components. This is accomplished using the
#Hecke action on q-expansions and linear algebra.
#
#
    print 'Computing the isotypic components...'
    omega1,eta1=isotypic_bases(H1dR,M[fIndex],pdmat,False)
    omega2,eta2=isotypic_bases(H1dR,M[gIndex],pdmat,False)
    omegaf=[R(x) for x in omega1[0]]
    etaf=[R(x) for x in eta1[0]]
    omegag=[R(x) for x in omega2[0]]
    etag=[R(x) for x in eta2[0]]
    OMEGA=[(x/q).integral() for x in omegag]
    ETA=[(x/q).integral() for x in etag]
    omegaeta=omegag+etag
    OMEGAETA=OMEGA+ETA
    omegaETA=[]
    for i in [0..len(omegag)-1]:
        omegaETA.append([omegag[i]*y for y in OMEGA+ETA])
    for i in [0..len(etag)-1]:
        omegaETA.append([etag[i]*y for y in OMEGA+ETA])
#
#The computation of H^1_dR finishes here. The output is:
#omega,eta are lists which give a basis of H_dR^1(X)
#
#OMEGA,ETA are their primitives, normalized so that have value 0 at hte cusp oo. 
#
#omegaETA are the 1-forms on the upper half plane used for the iterated integral of omega_a·eta_a
#
#Now we compute the Poincare dual gammaf of the extensions of f to level N. This is done by computing the intersection product
#of f against a basis of H_1(X)[f] and the intersection product matrix of this basis, and then using
#linear algebra. Finally, we do a change of basis so that gammaf is in terms of our hypgens basis, with
#small lower-left entries.
#
    print 'Computing gammaf...'    
    fextensions = [R(f).subs({q:q^d}).truncate(n) for d in divisors(N/cond)]
    v1=[]
    fmsbasis=M[fIndex].integral_basis()
    fperiods = []
    for i in [0..len(fmsbasis)-1]:
        fperiods.append([periodS(x,fmsbasis[i],Gens,Mat) for x in omegaf+etaf])
    intersectionmatrix1 = []
    for i in [0..len(fmsbasis)-1]:
    	intersectionmatrix1.append([0 for j in [0..i]])
	for j in [i+1..len(fmsbasis)-1]:
	    intersectionmatrix1[i].append(round(real(sum([fperiods[i][k]*fperiods[j][k+len(omegaf)] - fperiods[i][k+len(omegaf)]*fperiods[j][k] for k in [0..len(omegaf)-1]])/(2*CC(pi)*CC(I)))))
    intersectionmatrix=Matrix(intersectionmatrix1) - Matrix(intersectionmatrix1).transpose()
    for x in fextensions:
    	xperiods = vector([periodS(x,fmsbasis[j],Gens,Mat)/(2*CC(pi)*CC(I)) for j in [0..len(fmsbasis)-1]])
        v1.append(intersectionmatrix^(-1)*xperiods)
    v2 = []
    for x in v1:
    	v3=[]
    	for i in [0..len(M)-1]:
    	    if i==fIndex:
	        v3+=[y for y in x]
	    else:
	        v3+=[0 for k in [0..M[i].dimension()-1]]
	v2.append(v3)
    fduals = [vector(x)*homologymatrix*MatInv for x in v2]
    gammafs=[]
    for x in fduals:
        gammafs.append([[x[i],Gens[i]] for i in [0..len(x)-1]])
#
#Here we finished the conmputation of gamma_f as a vector of elements in Gamma_0(N) with complex coefficients.
#
#
#
#Now we compute the iterated integrals of each simple tensor in H^1_dR(X)[g] \otimes H^1_dR(X)[g]
#with respect to the basis we found earlier. The Chow-Heegner points are linear combinations of
#these iterated integrals depending on how the correspondence in question acts on this basis.
#We take advantadge of the fact that the primitives and products of q-series have been pre-computed.
#
#We may ignore all tensors where both differentials are not holomorphic, since these don't figure into
#the computation of the Chow-Heegner points.
#
    print 'Computing the raw zhang points...'
    RawZh=[]
    for gamma in gammafs:
        RawZh1=[]
	for i in [0..len(omegaETA)-1]:
	    RawZh1.append([])
	    for j in [0..len(omegaETA[i])-1]:
	    	if min(i,j) >= len(omegaETA)/2:
		    RawZh1[i].append(0)
		else:
	            RawZh1[i].append(iteratedintegralv(CC(I)/N,omegaeta[i],OMEGAETA[i],omegaeta[j],OMEGAETA[j],omegaETA[i][j],gamma))
        RawZh.append(RawZh1)

#
#
#Now we need to compute the adjustment factors alpha. Note that alpha is not unique, and we choose alpha
#such that its Poincare pairing with all of the (non-holomorphic) meromorphic forms is 0. Afterwards,
#the adjustments are applied to the previously computed iterated integrals.
#
#
    print 'Computing the alphas...'
    alpha = []
    for i in [0..len(omegaETA)-1]:
	alpha.append([])
	for j in [0..len(omegaETA)-1]:
	    if j <= i or min(i,j)>=len(omegaETA)/2 or max(i,j)<len(omegaETA)/2:
	        alpha[i].append(0)
	    else:
	        pral = omegaETA[i][j].truncate(10)
		alphaholo = [pd(x,pral) for x in holomorphicforms]
		alphamero = [0 for x in meromorphicforms]
		alphaH1dR = alphaholo+alphamero
		H1dRCC = [R(x) for x in H1dR]
		alpha[i].append(R((pdmat^(-1)*vector(alphaH1dR))*vector(H1dRCC)))
    print 'Computing the adjustments...'
    adjustments=[]
    for i in [0..len(gammafs)-1]:
        adjustments1=[]
	for j in [0..len(alpha)-1]:
	    adjustments1.append([])
	    for k in [0..len(alpha[j])-1]:
	        if alpha[j][k]==0:
		    adjustments1[j].append(0)
		else:
		    adjustments1[j].append(periodv(alpha[j][k],gammafs[i]))
        adjustments.append(Matrix(adjustments1) - Matrix(adjustments1).transpose())
    AdjZh=[]
    for i in [0..len(RawZh)-1]:
        AdjZh.append(Matrix(RawZh[i]) - adjustments[i])
#
#Finally, we compute a basis for the g-isotypic component of the Hecke algebra. Then we determine the
#linear combinations of the iterated integrals that compute these points, and finally compute the points
#themselves.
#
#Additionally, as our idempotents are QQ-linear combinations of cycles rather than cycles themselves,
#we must compute the denominators of our points in E(QQ) \otimes QQ, and then multiply through
#to get actual points in E(QQ).
#
    heckegens,heckemat=TtensorQgens(M[gIndex])
    T=HeckeAlgebra(M)
    hecke_denoms=[hecke_denominator(x,M,gIndex,T) for x in heckegens]
    dict={}
    gdim=len(omegag)
    for i in [0..gdim-1]:
        dict[(i,i+gdim)]=1
    symp1=Matrix(2*gdim,dict)
    symp=symp1-symp1.transpose()
    periodmat=poincare_matrix(omegag+etag,M[gIndex].basis(),Mat,Gens)
    hodgeclasses=[symp*periodmat*x.transpose()*periodmat^(-1) for x in heckemat]
    CCpoints=[]
    for k in [0..len(AdjZh)-1]:
        CCpoints.append([])
	for j in [0..len(hodgeclasses)-1]:
	    CCpoints[k].append((heckegens[j],hecke_denoms[j],hecke_denoms[j]*sum([AdjZh[k].rows()[i]*hodgeclasses[j].rows()[i] for i in [0..len(hodgeclasses[j].rows())-1]])))
    P=E.gens()[0]
    z=L.elliptic_logarithm(P)
    fdivs=divisors(N/cond)
    print 'For the generator ',P,' of E, the points are given by:'
    for i in [0..len(CCpoints)-1]:
        print 'For the extension f(q^',fdivs[i],'):'
        for j in [0..len(CCpoints[i])-1]:
    	    lindep = gp.lindep([CCpoints[i][j][2],L.gens()[0],L.gens()[1],z],10)
	    print 'n =',CCpoints[i][j][0],', d_{g,n} =',CCpoints[i][j][1],', d_{g,n}P_{g,f,n} =',-lindep[1]*lindep[4],'*P'
    QQpoints=[]
    for i in [0..len(CCpoints)-1]:
        QQpoints.append([])
	for j in [0..len(CCpoints[i])-1]:
	    QQpoints[i].append((CCpoints[i][j][0],CCpoints[i][j][1],L.elliptic_exponential(CCpoints[i][j][2])))
    return CCpoints, QQpoints