\\ new eulerchar calculation for the fine selmer group

if(fineisinitialized != 1,\
 print("First initalisation of fine.gp. I will kill the variables t,w");\
 fineisinitialized =1;\
 read("~/global/pari/lseriesderi.gp");\
 kill(t); kill(w); kill(z); w;t;);

\\bigO is a global variable for the precision

bigO=20;

debug_verbose=2;
output_verbose=1;

   \\drops the first coordinates of the vector v that equal x, or v itself
   drop(v,x) =
    {local(ind);
   		ind=1;
		while(ind<=length(v) & v[ind] != x,ind++);
		if(ind<=length(v),return(vecextract(v,concat("^",ind))),return(v))} 
	
	redissmooth(el,P)= (P==[0])||(el.a1*P[2]!=3*P[1]^2+2*el.a2*P[1]+el.a4)||(2*P[2]+el.a1*P[1]+el.a3!=0)


\\ ANALYTIC FUNCTIONS =================================================

\\ elliptic logarithm
\\ el is an elliptic curve and bigO the series precision
\\ the output is a series in t with coefficients in the field of
\\ definition of el

elllog(el) =
 { local(f,s,xs);
 		if(debug_verbose>2,print("start elllog"));
	f = Mod(t^3,t^bigO) + Mod(el.a1 * t  + el.a2 * t^2,t^bigO) * w + Mod(el.a3 + el.a4 * t,t^bigO) * w^2 + Mod(el.a6,t^bigO)* w^3;
	s = f;
	while(poldegree(s,w)>0,
		s=truncate(subst(s,w,f)+O(w)^bigO));
	s = (lift(s)+O(t)^bigO);
	xs = t / s;
	return(intformal(1/(-2/s+el.a1*xs+el.a3)*deriv(xs,t),t)); }

\\the sigma function of Bernardi
\\ it is a series in t of prec 2* \ps- 

ells0(el) =
 { local(sz);
 		if(debug_verbose>2,print("start ells0"));
    sz = z*exp(-intformal(intformal(ellwp(el,z)-1/z^2))); 
	return(subst(sz,z,elllog(el)));}
 
\\ HEIGHTS =================================================================

\\the p-adic height of a point P in the formal group with
\\good reduction at all primes

\\using Bernardi's sigma function s0

h_sigma(el,P,p,s0) = log(ee(P)/subst(s0,t,tt(P)+O(p)^bigO));

hacc_sigma(el,P,Q,p,s0)=   h_sigma(el,elladd(el,P,Q),p,s0) - h_sigma(el,P,p,s0) - h_sigma(el,Q,p,s0);

\\the height matrix associated to a basis ba

hmat_sigma(el,ba,p,s0) =
 {local(rank,hm,ind,jnd);
 	if(debug_verbose>1,print("start hmat_sigma"));
  rank=length(ba);
  hm=matrix(rank,rank);
  for(ind=1,rank,for(jnd=1,rank,
   hm[ind,jnd]=hacc_sigma(el,ba[ind],ba[jnd],p,s0)));
  return(hm)}

\\------------------------------------------------------------------------------
\\this is the approximation given by sigma =t + O(t)^***
  
  prec_naive(a1,a2,a3,p,vt)=if(p==3,\
  		if(a1==0,if(a2==0,if(a3==0,5*vt-2,if(vt==1,3,4*vt)),3*vt-1),2*vt),\
		if(a1==0,if(a2==0,if(a3==0,5*vt,4*vt),3*vt),2*vt))

  
h_naive(el,P,p) = log(ee(P)/(tt(P)+O(p)^prec_naive(el[1],el[2],el[3],p,valuation(tt(P),p))))

hacc_naive(el,P,Q,p)=  h_naive(el,elladd(el,P,Q),p) - h_naive(el,P,p) - h_naive(el,Q,p);

\\the height matrix associated to a basis ba
hmat_naive(el,ba,p) =
 {local(rank,hm,ind,jnd);
 	if(debug_verbose>1,print("start hmat_naive"));
  rank=length(ba);
  hm=matrix(rank,rank);
  for(ind=1,rank,for(jnd=1,rank,
   hm[ind,jnd]=hacc_naive(el,ba[ind],ba[jnd],p)));
  return(hm)}

\\Tate's trick to take the limit of naive heights
\\turns out to be much slower.

h_limit(el,P,p,n)=1/p^(2*n)*h_naive(el,ellpow(el,P,p^n),p) 



\\ THE KERNEL OF LOG ===============================================

\\given en elliptic curve el and a set ba=[Q_1,...,Q_r] 
\\in B = E^(Q_p) n E(Q)

\\this gives the vector of p-adic elliptic logarithms of Q_i 
\\only over Q

logsbasis(el,ba,p)=
 {local(ind2,zs,rank);
 		if(debug_verbose>1,print("start logsbasis"));
 	rank=length(ba);
	zs=truncate(elllog(el));
	return(vector(rank,ind2,subst(zs,t,tt(ba[ind2])+O(p)^bigO)));}
	
\\ this is a list of vectors f_1 = [x_1,  .... x_r] put in a matrix
\\ such that x1 * Q_1 + .... + x_r *Q_r is in the kernel of log.
\\ these vectors form a basis of it.
\\ in perticular if ba is a Z-basis of B, the is a Z_p-basis of M
\\ and if ba is a Z-basis of B_S, then it is of M_S

finebasis(el,ba,p) =
 { local(rank,ind,logs,mini,holder,resvec);
 		if(debug_verbose>1,print("start finebasis"));
	rank=length(ba);
	if(rank<1,return([]),
		logs=logsbasis(el,ba,p);
	\\pivot
		mini=1000;holder=0;
		for(ind=1,rank,
			if(valuation(logs[ind],p)<mini,
				holder=ind;
				mini=valuation(logs[ind],p)));
		resvec=matrix(rank,rank-1);
		for(ind=1,rank,
			if(ind<holder,
				resvec[ind,ind]=1;
				resvec[holder,ind]=-logs[ind]/logs[holder],
				if(ind>holder,
					resvec[ind,ind-1]=1;
					resvec[holder,ind-1]=-logs[ind]/logs[holder])));
		return(resvec))}


\\ BASIS OF THE FORMAL B = E(Q)  n E^ ==================================
										
\\find a basis of B. 

\\    !!! stupid programming !!!

\\given a basis ba of the Mordell-Weil, hopefully of small height, this produces
\\ [M,nba]
\\ where nba is a vector containing a basis of B
\\ and M is the transformation matrix.


epbasis(el,ba,p)=
 {local(cp,nba,nmat);
 	if(length(ba)==0,return([1,[]]),
	cp=tama(el,p);
	if(cp!=1,
		nba=epbasis_bad(el,ba,p);nmat=nba[1];nba=nba[2],
		nba=ba;nmat=1);
	nba=epbasis_good(el,nba,p);
	return([nmat*nba[1],nba[2]]) ) }

epindex(el,ba,p)=abs(matdet(epbasis(el,ba,p)[1]))


\\first we suppose that none of the vectors in ba map to the singular point.

\\this is very slow ! 

epbasis_good(el,ba,p)=
 {local(np,rank,bared,ind,dvec,d,P,count,amat,aredmat,jnd,newba,v,image,boo,bro);
	np=ellnp(el,p); \\ order of E~ns
	rank=length(ba);
	bared=vector(rank);
	for(ind=1,rank,
		if(ee(ba[ind])%p,
			bared[ind]=Mod(ba[ind],p),
			bared[ind]=[0]));
	 \\ find the order of the points in ba, when reduced
	dvec=vector(rank,ind,1);
	for(ind=1,rank,
		if(ee(ba[ind])%p, \\if it is in the formal group, d_i=1, otherwise
			fordiv(np,d,
				if(ellpow(el,bared[ind],d)==[0],dvec[ind]=d;break)))); 
		if(debug_verbose>1,print(" found orders :"dvec));
	 \\ find all combinations a_i P_i +... in the formal group
	 \\ count the points found in the reduction
	 image=listcreate(np-1); bro=0;
	 count=rank; amat=matrix(rank,20+rank);
	 for(ind=1,rank,
	 	amat[ind,ind]=dvec[ind]); 
	 forvec(avec=vector(rank,ind,[0,dvec[ind]-1]),
	 	P=ellpow(el,bared[1],avec[1]);
		for(jnd=2,rank,
			P=elladd(el,P,ellpow(el,bared[jnd],avec[jnd])));
		if(P==[0],
			count++;
			amat[,count]=avec~;
			  \\reduce the matrix every 20th time
			if(count==rank+1,    \\20,
				aredmat=amat*qflll(amat,1);
				amat=matrix(rank,20+rank);
				for(jnd=1,length(aredmat),amat[,jnd]=aredmat[,jnd]);
				count=length(aredmat);
				if(debug_verbose>2,
					print("reduced to "amat),
					if(debug_verbose>1,print1("r")));
				\\check if we found as many points in the image as the index found so far
				if(length(image)+1==abs(matdet(aredmat)),
					bro=1;
					if(debug_verbose > 1,print(" found enough points :"length(image)+1));
					break);),
			\\the reduction is not 0, append it to the list image, if new
			boo=0;
			for(knd=1,length(image),
				if(P==image[knd],boo=1;break) );
			if(boo==0,
				listput(image,P));
			));
	\\ reduce a last time
	aredmat=amat*qflll(amat,1);
	   if(debug_verbose>1,print1(" reduced basis of the lattice :"aredmat);print(" of det "matdet(aredmat)" compared with N_p = "np));
	newba=vector(rank);
	for(ind=1,rank,
		v=aredmat[,ind]~;
	 	P=ellpow(el,ba[1],v[1]);
		for(jnd=2,rank,
			P=elladd(el,P,ellpow(el,ba[jnd],v[jnd])););
		newba[ind]=P;);
	return([aredmat,newba]);}


\\ this gives a basis of W n E^0(Q_p). if W is given by a basis ba
\\ and p is any prime.

	
epbasis_bad(el,ba,p)=
 {local(cp,rank,count,amat,ind,avec,jnd,P,aredmat,newba,v);
	cp=tama(el,p); \\ order of Phi_p
	rank=length(ba);
	 \\ find all cobinations a_i P_i +... that lie in E^0(Q_p)
	 count=rank; amat=matrix(rank,20+rank);
	 for(ind=1,rank,
	 	amat[ind,ind]=cp); 
	 forvec(avec=vector(rank,ind,[0,cp-1]),
	 	P=ellpow(el,ba[1],avec[1]);
		for(jnd=2,rank,
			P=elladd(el,P,ellpow(el,ba[jnd],avec[jnd])));
		if(issmooth(el,P,p)==1,
			count++;
			amat[,count]=avec~;
			  \\reduce the matrix every 20th time
			if(count==20,
				aredmat=amat*qflll(amat,1);
				amat=matrix(rank,20+rank);
				for(jnd=1,length(aredmat),amat[,jnd]=aredmat[,jnd]);
				count=length(aredmat);
				if(debug_verbose>2,
					print("reduced to "amat),
					if(debug_verbose>1,print1("r")));)));
	\\ reduce a last time
	  if(debug_verbose>2,print(" unreduced generators : "amat));
	aredmat=amat*qflll(amat,1);
	   if(debug_verbose>1,print1(" reduced basis of the lattice :"aredmat);print(" of det "matdet(aredmat)" compared with c_p = "cp));
	newba=vector(rank);
	for(ind=1,rank,
		v=aredmat[,ind]~;
	 	P=ellpow(el,ba[1],v[1]);
		for(jnd=2,rank,
			P=elladd(el,P,ellpow(el,ba[jnd],v[jnd])););
		newba[ind]=P;);
	 if(debug_verbose> 2,print("new basis of good reduction : "newba));
	return([aredmat,newba]);}
	
	
\\ INDEX OF M_S IN M =================================================
\\ this should calculate the index of the finest mordell-weil group
\\ inside the fine mordell-weil group

\\ it is the order of the image of M in the product of E(Q_v)*
\\ for bad v different from p.

findex(el,ba,p)=
 {local(epba,bp,c,ind,logs);
     if(debug_verbose>1,print("start findex "));
 	bp=drop(badprimes(el),p);
	c=p^valuation(ellglobalred(el)[3]/tama(el,p),p);
	c=c*p^(sum(ind=1,length(bp),valuation(ellnp(el,bp[ind]),p)));
	 if(debug_verbose>1,print(" upper bound :"c));
	 \\c is an upper bound for the index
	if(c==1,return(1),
		epba=epbasis(el,ba,p)[2];
		logs=finebasis(el,epba,p);
		return(findex_1(el,epba,logs,p)))}

findex_2(el,epba,p)=
 {local(bp,c,ind,logs);
     if(debug_verbose>1,print("start findex _2 "));
 	bp=drop(badprimes(el),p);
	c=p^valuation(ellglobalred(el)[3]/tama(el,p),p);
	c=c*p^(sum(ind=1,length(bp),valuation(ellnp(el,bp[ind]),p)));
	 if(debug_verbose>1,print(" upper bound :"c));
	 \\c is an upper bound for the index
	if(c==1,return(1),
		logs=finebasis(el,epba,p);
		return(findex_1(el,epba,logs,p)))}

findex_1(el,epba,logs,p)=
 {local(rank,bp,bpp,np,c,fba,nba,nmba,nnba,knd,count,amat,avec,P,ind,jnd,boo,aredmat);
 	 if(debug_verbose>1,print("start findex_1"));
	rank=length(ba);
	bp=drop(badprimes(el),p);
	c=p^valuation(ellglobalred(el)[3]/tama(el,p),p);
	c=c*p^(sum(ind=1,length(bp),valuation(ellnp(el,bp[ind]),p)));
	 \\c is an upper bound for the index
	 if(debug_verbose>1,print(" upper bound :"c));
	if(c==1,return(1),  \\should happen almost always
	if(c>p,
		if(p>20,print("findex not yet implemented for c="c" and p="p);return(0),
		 \\ this is the slow programming
		 \\ calculate a Z/cZ-basis of M represented by points on E(Q)
	     if(debug_verbose>1,print("enter complicated case")); 
		fba=lift(Mod(logs,c));
		fba=fba*qflll(fba,1);
		 if(debug_verbose>2,print("lifted f-basis is: "fba));
		nba=vector(rank-1);
		for(ind=1,rank-1,
			P=ellpow(el,epba[1],fba[1,ind]);
			for(jnd=2,rank,
				P=elladd(el,P,ellpow(el,epba[jnd],fba[jnd,ind])));
			nba[ind]=P;);;
		\\ now nba represents a Z/cZ-basis od M/cM
		find=1;
		for(knd=1,length(bp),
			if(tama(el,bp[knd])!=1,
				nnba=epbasis_bad(el,nba,bp[knd]);
				nba=nnba[2];
				find=find*p^valuation(abs(matdet(nnba[1])),p);));
		\\ now all elements in nba are smooth !
		\\ so we can work modulo bad primes whose reduction has p-primary parts
		bpp=[];	
		for(knd=1,length(bp),
			if(Mod(ellnp(el,bp[knd]),p)==0,bpp=concat(bpp,bp[knd])) );
		 if(debug_verbose>1,print(" index so far ::"find", only primes to consider now :"bpp));
		if(length(bpp)==0,
			return(find), \\done
			 \\	
			nmba=vector(rank-1,ind,vector(length(bpp)));
			for(ind=1,rank-1,
				for(knd=1,length(bpp),
					if(valuation(nba[ind][1],bpp[knd])<0,
						nmba[ind][knd]=[0],
						nmba[ind][knd]=Mod(nba[ind],bpp[knd]) )));
			  if(debug_verbose>1,print(" basis ",nmba));
			  \\ multiply with non-p-primary parts
			c=1;
			for(knd=1,length(bpp),
				for(ind=1,rank-1,
				np=ellnp(el,bpp[knd]);
				c=c*p^valuation(np,p);
				nmba[ind][knd]=ellpow(el,nmba[ind][knd],np/p^valuation(np,p)) ));
			  if(debug_verbose>1,print(" basis ",nmba));
			 \\ construct all combination
		 	count=rank-1; amat=matrix(rank-1,20+rank);
			for(ind=1,rank-1,amat[ind,ind]=c); 
			forvec(avec=vector(rank-1,ind,[0,c-1]),
				boo=1;knd=0;
				while(knd<length(bpp) && boo==1,
		 			knd++;
					P=ellpow(el,nmba[1][knd],avec[1]); 
					for(jnd=2,rank-1,
						P=elladd(el,P,ellpow(el,nmba[jnd][knd],avec[jnd])) );
					if(P!=[0],boo=0) );
				if(boo==1,
					count++;
					amat[,count]=avec~;
					  \\reduce the matrix every 20th time
					if(count==20,
						aredmat=amat*qflll(amat,1);
						amat=matrix(rank-1,20+rank);
						for(jnd=1,length(aredmat),amat[,jnd]=aredmat[,jnd]);
						count=length(aredmat);
						if(debug_verbose>2,
							print("reduced to "aredmat),
							if(debug_verbose>1,print1("r")) );)));
			\\ reduce a last time
			aredmat=amat*qflll(amat,1);
			 if(debug_verbose>1,print1(" reduced basis of the lattice :"aredmat);print(" of det "matdet(aredmat)));
		    return(find*abs(matdet(aredmat))) )),
	\\case when c=p
		if(debug_verbose>1,print("enter easy case"));
		 \\ calculate a Z/cZ-basis of M represented by points on E(Q)
		fba=lift(Mod(logs,c));
		fba=fba*qflll(fba,1);
		 if(debug_verbose>2,print("lifted f-basis is: "fba));
		vu=0;
		for(knd=1,length(bp),if(Mod(tama(el,bp[knd]),p)==0,vu=bp[knd]));
		if(vu!=0,
		  \\ in this case we know that Phi(Q_vu)(p)=Z/pZ + factors prime to p
		  \\ all we need to check is if the basis fba maps
		  \\ surjectively on the Z/pZ or if all have order prime to p
		  \\ in Phi(Q_v)(p)
		  	if(debug_verbose>1,print(" "p"-primary Phi at "vu));
			boo=1;
			for(ind=1,rank-1,
				P=ellpow(el,epba[1],fba[1,ind]);
				for(jnd=2,rank,
					P=elladd(el,P,ellpow(el,epba[jnd],fba[jnd,ind])));
				if(issmooth(el,P,vu)==0,
					P=ellpow(el,P,tama(el,vu)/p);
					if(issmooth(el,P,vu)==0,boo=p;break;)));
			return(boo),
			\\
			for(knd=1,length(bp),if(Mod(ellnp(el,bp[knd]),p)==0,vu=bp[knd]));
			 \\ in this case, we know that E~_ns(F_vu) is cyclic
		  	 \\ of order dividing p exactly once. But Phi is prime to p.
			boo=1;if(debug_verbose>1,print(" "p"-primary reduction at "vu));
		 	epbared=vector(rank);
			for(ind=1,rank,
				if(issmooth(el,epba[ind],vu)==0,
						\\sortir de la merde.
						\\ ATTENTION. Mod(ellpow(el,P,n),p) =!= ellpow(el,Mod(P,p),n)
						\\ quand P n'est pas lisse.
					P=ellpow(el,epba[ind],tama(el,vu));
					if(valuation(P[1],vu)<0,
						epbared[ind]=[0],
						epbared[ind]=Mod(P,vu) ),
					if(valuation(epba[ind][1],vu)<0,
						epbared[ind]=[0],
						epbared[ind]=Mod(epba[ind],vu) )));
			for(ind=1,rank-1,
				P=ellpow(el,epbared[1],fba[1,ind]);
				for(jnd=2,rank,
					P=elladd(el,P,ellpow(el,epbared[jnd],fba[jnd,ind])));
				P=ellpow(el,P,ellnp(el,vu)/p);
				if(P!=[0],boo=p;break));
			return(boo))))}




\\ REGULATORS ======================================================================

\\ here we give the valuation of the p-adic regulator of the
\\ fine subgroup in this
\\ case using sigma-functions.

\\determines the valuation of the regulator of the p-adic height pairing on the
\\fine Mordell-Weil.
\\ el is an elliptic curve in minimal form,
\\ p any ODD prime and 
\\ ba a basis of E(Q) (or at least of index prime to p)

regfine(el,ba,p)=
 {local(s0,epba);
	s0=truncate(ells0(el));
	nba=epbasis(el,ba,p)[2];
	return(regfine_sigma(el,nba,p,s0)) }

regfine_sigma(el,nba,p,s0) =
 {local(rank,c,bp,mba,index,ind,hm,fba);
 	rank=length(ba);
	 \\ we calculate a basis of B= E n E^
	\\nba=epbasis(el,ba,p)[2];
	mba=vector(rank);
	 \\ shift bad points
	bp=drop(badprimes(el),p);
	if(length(bp)==0,c=1,
		c=tama(el,bp[1]);
		for(ind=2,length(bp),
			c=lcm(c,tama(el,bp[ind])) ;););
		\\c=lcm(tamvec); \\is smaller than ellglobalred[3].
	if(debug_verbose>1,print("factor to multiply into E0 :"c));
	index=vector(rank,ind,1);
	for(ind=1,rank,
		if(isallsmooth(el,nba[ind])==0,
			mba[ind]=ellpow(el,nba[ind],c);
			index[ind]=1/c,
			mba[ind]=nba[ind]));
	 \\ calcualte the height matrix with respect to the sigma fn of bernardi
	index=matdiagonal(index);
	hm=hmat_sigma(el,mba,p,s0); 
	hm=index*hm*index;
	 \\ calculate now the kernel of log
	fba=finebasis(el,nba,p);
	return(valuation(matdet(fba~ * hm * fba ),p))}

\\==== COKERNEL OF LOCALISATION IN E(Q_p)* ================

\\ given a basis ba of the free part of E(Q), this calculate the valuation #D(p) 
\\ of the torsion part of the
\\ index of the subgroup genereated by E(Q)* in E(Q_p)* 

\\ Here p>2 is important.

 Dindex(el,ba,p) =
  { local(epba);
  	if(length(ba)>0,epba=epbasis(el,ba,p)[2],epba=[]);
	return(Dindex_ba(el,ba,epba,p))}


 Dindex_ba(el,ba,epba,p) =
 { local(rank,torspts,qts,ind,lpts,c,dind,boo,vvec);
	rank=length(ba);
	if(rank==0,
		return(valuation(localtors(el,p)/elltors(el)[1],p)),
		 \\rank >0
		torspts=elltors(el)[3];
 		qts=vector(rank+length(torspts));
		for(ind=1,rank,qts[ind]=ba[ind]);
		for(ind=1,length(torspts),qts[rank+ind]=torspts[ind]);
		qts=qts*(1+O(p)^30);
		lpts=length(qts);
			\\ surjective on Phi* ?
		c=tama(el,p)/p^valuation(tama(el,p),p);
		for(ind=1,lpts, qts[ind]=ellpow(el,qts[ind],c));
		if(Mod(tama(el,p),p) == 0,
			if(valuation(tama(el,p),p)>1,
				print(" *** Tamagawa :"tama(el,p)" too high valuation at "p", not yet implemented ***");
				break(); );
			dind=1;
			for(ind=1,lpts,
				if(localissmooth(el,qts[ind],p) ==0,
					dind=0;
					qts[ind]=ellpow(el,qts[ind],p) )););
			  \\qts are now in E^0(Q_p), surjectivity on N_p
		if(Mod(ellnp(el,p),p)==0,
			boo=1;
			c=ellnp(el,p)/p^valuation(ellnp(el,p),p);  \\only for p=3 and N_3=6, then c=2
			for(ind=1,lpts, qts[ind]=ellpow(el,qts[ind],c));
			for(ind=1,lpts,
				if(qts[ind]!=[0] && valuation(qts[ind][1],p)>-1 , \\not in the formal group
					boo=0; ););
			dind=dind+boo; );
			   \\qts are now in the formal group, index here is the min of valutations of the non-torsion
			   if(debug_verbose>2,print(" into the formal group : "qts*(1+O(p)^2)));
		\\epba=epbasis(el,vecextract(qts,concat("1..",rank)),p)[2];
		vvec=vector(rank);
		for(ind=1,rank,vvec[ind]=-valuation(epba[ind][1],p)/2);
		dind=dind+vecmin(vvec)-1;
		return(dind); ) }
		
				
\\ LOCAL TORSION ====================================================
\\
\\ check if E(Q_p) has p-torsion.
\\
\\ 0 -> E(Q_p)(p) -> E~(F_p)(p) -d-> Q_p/Z_p -> E^0(Q_p) ox Q_p/Z_p -> 0
\\
\\ 0 -> E^0(Q_p)(p) -> E(Q_p)(p) -> Phi(Q_p)(p) -d'-> E^0(Q_p) ox Q_p/Z_p ->  E(Q_p) ox Q_p/Z_p -> 0

\\ el must be an elliptic curve over Q in minimal form
\\ p is any ODD prime
\\ the result is the order of E(Q_p)(p)

localtors(el,p) =
 {local(cp,ci,sp,boo,x0,y0,ellf,Q);
	cp=p^valuation(tama(el,p),p);
	if(cp>1,
		\\if(cp>p,print(" **** high valuation of c_p ::: possibly wrong **** "));
		  \\ find singular point
		  if(debug_verbose>2,print("search for a singular point"));
		 sp=lift(singpoint(el,p));
		  \\ lift to a point in Q_p with singular reduction
		 Q=findlift(el,sp,p); if(debug_verbose>2,print("...found "Q));
		 \\ check if Q is a generator of the cyclic group Phi(Q_p)(p),
		 \\ otherwise search for another point
		 if(cp>p,
		 	while(localissmooth(el,ellpow(el,Q,tama(el,p)/p),p),
				sp=sp+[p,0];
				Q=findlift(el,sp,p);
				if(debug_verbose>2,print("changed to "Q)); ));
		 Q=ellpow(el,Q,tama(el,p));
		 Q=ellpow(el,Q,ellnp(el,p));
		 if(-valuation(Q[1],p)/2 > valuation(cp,p) || Q==[0],
		 	\\d' is trivial
			ci=cp,
			ci=p^(-valuation(Q[1],p)/2-1) ),
		ci=1;);
     if(Mod(ellnp(el,p),p)!=0, return(ci),
	 	 	\\ Choose any point in E^0 not in E^.
			if(debug_verbose>2,print("find a non-trivial element in E^0 mod E^"));
			boo=0;x0=O(p)^5;
			while(boo==0,
				ellf=Y^2+(el.a1*x0+el.a3)*Y-(x0^3+el.a2*x0^2+el.a4*x0+el.a6);
				y0=polrootspadic(ellf,p,5);
				if(length(y0)==0,x0=x0+1,
					y0=y0[1];
					if(localissmooth(el,[x0,y0],p)==1, \\is [x0,y0] in E^0
						boo=1,
						x0=x0+1 )));
			 \\ [x0,y0] is the lift of a non-trivial point
			 \\ in E~ = Z/pZ. The connecting hom d is obtained
			 \\ by multiplying with p
			 \\ d is trivial if Q is in E^(p^2) and injective if not
			 if(debug_verbose>2,print("...found ",[x0+O(p)^3,y0+O(p)^3]));
			Q=ellpow(el,[x0,y0],ellnp(el,p));
			if(debug_verbose>2,print("multiplies to ",Q*(1+O(p)^2)));
			if(Q==[0]||valuation(Q[1],p)<-2,
				\\ d is trivial
				if(debug_verbose>2,print(" hence delta is trivial "));
				return(ci*p),
				\\ d is injective
				if(debug_verbose>2,print(" hence delta is non-trivial "));
				return(ci) )); }


\\checks is a local point or a global point is smooth at p
localissmooth(el, P, p) =
 { local(xp, yp);
	 if(P==[0],return(1),
   		if(valuation(P[1],p)<0,return(1),
      		yp=Mod(P[2],p);
	   		xp=Mod(P[1],p);
	   		if((el.a1*yp==3*xp^2+2*el.a2*xp+el.a4)&&(2*yp+el.a1*xp+el.a3==0),
	   			return(0),
				return(1)))) }

\\ find the singular point in the reduction.. yes this is a stupid code, try to make it better.
singpoint(el,p)=
 { local(xx,xsing,ysing,boo);
 	if(Mod(e.disc,p)!=0 || p==2, return([]),
	boo=0;xx=0;
 	while(boo==0&& xx<p,
		if(Mod(xx,p)^3+el.a2*xx^2+el.a4*xx+el.a6+(el.a1*xx+el.a3)^2/4 == 0,
		xsing=Mod(xx,p);
		ysing=-(el.a3+el.a1*xsing)/2;
		if(el.a1*ysing == 3*xsing^2+2*el.a2*xsing+el.a4,boo=1); );
		xx++; );
	return([xsing,ysing]) );}

\\ellordinate for p-adics
elly(el,x,p,prec)= polrootspadic(Y^2+(el.a1*x+el.a3)*Y-(x^3+el.a2*x^2+el.a4*x+el.a6),p,prec)

\\ this should find à la Hensel a point close to sp. used only in localtors !
findlift(el,sp,p) =
{local(x0,y0,boof,preci,x1);
	if(debug_verbose >2 ,print(" enter findlift "));
	if(ellisoncurve(e,Mod(sp,p))==0,
		if(debug_verbose >2 ,print(" point "sp" is not on the cruve mod "p));
		return([]),
		if(Mod(2*sp[2]+(el.a1*sp[1]+el.a3),p)!=0,
			if(debug_verbose >2 ,print(" del_Y !=0 :: lift with x="sp[1]));
			y0=elly(el,sp[1],p,10);	return([sp[1],y0[1]]),
			if(Mod(-el.a1*sp[2]+3*sp[1]^2+2*el.a2*sp[1]+el.a4,p)!=0,
				if(debug_verbose >2 ,print(" del_X !=0 :: lift with y="sp[2]));
				x0=polrootspadic(sp[2]^2+(el.a1*X+el.a3)*sp[2]-(X^3+el.a2*X^2+el.a4*X+el.a6),p,20);
				return([x0[1],sp[2]]),
				\\ otherwise the reduction of sp is singular
				if(elllocalred(e,p)[4]==1,
					if(debug_verbose >2 ,print(" singular point and c_p=1 :: no lift can exist. "));
					return([]), \\can't have any solutions
					boof=0; preci=0;
					if(debug_verbose >2 ,print(" singular point and c_p="elllocalred(e,p)[4]));
					while(boof==0,
						preci++;
						if(debug_verbose>2 && preci>1, print(" ... increase preci to "preci));
						if(preci>10,error("problem in findlift"));
						x0=sp[1]-p;
						while(boof==0 && x0<p^preci,
							x0=x0+p;
							y0=elly(e,x0,p,20);
							if(debug_verbose >4 ,print1(" "x0" : "length(y0)));
							if(length(y0)>0,boof=1);
						);
					);
					return([x0,y0[1]]);
				);			
			);
		);
	);			
} 

				
\\ EULER-CHARACTERISTIC ====================================================
\\

\\ Euler characteristic bound


eulerfine(el,ba,p)=
 {local(bp,epba,s0,reg,rank,tamag,dind,fshabound,upperbound);
 	if(p==2,
		print("you will be punished in hell for doing Iwasawa-theory for p=2");
		return([]),
 		if(valuation(el.j,p)>=0, 
			\\	potential good reduction
		 	rank=length(ba);
			bp=badprimes(el);
			epba=epbasis(el,ba,p)[2];
			if(rank>1,	
				s0=truncate(ells0(el));
				reg = regfine_sigma(el,epba,p,s0)-(rank-1),
				reg = 0 );
			tamag=valuation(ellglobalred(el)[3]/tama(el,p),p);
			dind=Dindex_ba(el,ba,epba,p);
			
			fshabound=valuation(round(ellSha(el,ba)),p);
	
			upperbound= reg +dind +fshabound+tamag;

			return([upperbound,reg,dind,fshabound]),
			  \\
			  \\ not potential good reductionrank=length(ba);
			rank=length(ba);
			bp=badprimes(el);
			epba=epbasis(el,ba,p)[2];
			if(rank>1,	
				s0=truncate(ells0(el));
				reg = regfine_sigma(el,epba,p,s0)-(rank-1),
				reg = 0 );
			tamag=valuation(ellglobalred(el)[3]/tama(el,p),p);
			dind=Dindex_ba(el,ba,epba,p);
			
			
			fshabound=valuation(round(ellSha(el,ba)),p);
	
			upperbound= reg +dind +fshabound+tamag + valuation(elltors(el)[1],p);

			return([upperbound,reg,dind,fshabound]) )) }


\\=========== COMPLETE INFORMATION ================

fine(el,ba,p) = 
 {local(rank,bp,epba,s0,reg,tamag,dind,fshabound,lt,gt,upperbound,find,lred,lredstr);
 	rank=length(ba);
	bp=badprimes(el);
	if(rank>0,
		epba=epbasis(el,ba,p)[2];
		if(rank>1,	
			s0=truncate(ells0(el));
			reg = regfine_sigma(el,epba,p,s0)-(rank-1),
			reg = 0 ),
		epba=[];
		reg=0; );
	tamag=valuation(ellglobalred(el)[3]/tama(el,p),p);
	dind=Dindex_ba(el,ba,epba,p);
	fshabound=valuation(round(ellSha(el,ba)),p);
	lt=valuation(localtors(el,p),p);
	gt=valuation(elltors(el)[1],p);
	upperbound= reg +dind +fshabound+tamag + if(valuation(el.j,p)>=0,0,gt);
	if(rank>1,find=valuation(findex_2(el,epba,p),p),find=0); 
	lred=elllocalred(el,p); 
	  if(lred[2]==1, lredstr=concat(" good",if(ellnp(el,p)%p == 0,", anomalous "," ")));
	  if(lred[2]==1 && isordinary(el,p) ==0, lredstr=concat(lredstr,", supersingular "));
	  if(lred[2]==2, lredstr="II");
	  if(lred[2]==3, lredstr="III");
	  if(lred[2]==4, lredstr="IV");
	  if(lred[2]==-1, lredstr="I_0 *");
	  if(lred[2]==-2, lredstr="II*");
	  if(lred[2]==-3, lredstr="III*");
	  if(lred[2]==-4, lredstr="IV*");
	  if(lred[2]>4, lredstr=concat("I_",lred[2]-4));
	  if(lred[2]<-4, lredstr=concat(concat("I_",- lred[2] -4)," *"));
	  if(lred[2]!=1,lredstr=concat(lredstr,concat(" c = ",lred[4]))); 

	 \\output
	if(output_verbose >0,
		print(" curve of conductor     "conductor(el));
		print(" reduction at p         "lredstr);
		print(" global torsion         "p"^"gt);
		print(" local torsion          "p"^"lt);
		print(" tamagawa outside p     "p"^"tamag);
		print(" regulator              "p"^"reg);
		print(" analytic Sha(p)        "p"^"fshabound);
		print(" # D(p)                 "p"^"dind);
		print(" Index M : M_S          "p"^"find);
		print("");
		print(" bound on chi           "p"^"upperbound);
		); 
	return([upperbound,reg,dind,fshabound,find,lt]);}