load "newton.magma";

R<x>:=PolynomialRing(Rationals());
R<T>:=PolynomialRing(Rationals());
C:=CremonaDatabase();

function ap(E,p)
	F:=FiniteField(p);
	return(Trace(E,F));
end function;

function MSEven(M,r)
 	return(ModularSymbolEven(M,r)[1]/4);
end function;

function MSOdd(M,r)
 	return(ModularSymbolOdd(M,r)[1]/4);
end function;

function MS_twist(M,r,D)
	answer:=0;
	for a in [0..Abs(D)-1] do
		if KroneckerSymbol(D,a) ne 0 then
			if D ge 0 then
				ms:=MSEven(M[1],r-a/D);
			else
				ms:=MSOdd(M[2],r+a/D);
			end if;
			answer:=answer+ms*KroneckerSymbol(D,a);
		end if;
	end for;

	return(answer);
end function;

Cremona:=CremonaDatabase();

function form_curve_ap(N,a_2,a_3,a_5,a_7)
	R<x>:=PolynomialRing(Rationals());
	M1:=ModularSymbols(N,2,+1);
	M2:=ModularSymbols(N,2,-1);
	I:=[<2,x-a_2>,<3,x-a_3>,<5,x-a_5>,<7,x-a_7>]; 
	V1:=Kernel(I,M1);
	V2:=Kernel(I,M2);
	MSEven(V1,1/2);
	ModularSymbolOdd(V2,1/2);

	return([V1,V2]);
end function;

function form_curve(N,k)
	E:=EllipticCurve(Cremona,N,k,1);
	a_2:=ap(E,2);
	a_3:=ap(E,3);
	a_5:=ap(E,5);
	a_7:=ap(E,7);
	
	return form_curve_ap(N,a_2,a_3,a_5,a_7);
end function;


function alphabet(k)
	case k:
		when 1: return("A");
		when 2: return("B");
		when 3: return("C");
		when 4: return("D");
		when 5: return("E");
		when 6: return("F");
		when 7: return("G");
		when 8: return("H");
		when 9: return("I");
		when 10: return("J");
		when 11: return("K");
		when 12: return("L");
		when 13: return("M");
		when 14: return("N");
		when 15: return("O");
		when 16: return("P");
		when 17: return("Q");
		when 18: return("R");
		when 19: return("S");
		when 20: return("T");
		when 21: return("U");
		when 22: return("V");
		when 23: return("W");
		when 24: return("X");
		when 25: return("Y");
		when 26: return("Z");
	end case;
end function;

function form_filename(N,k,p,sign)
  return(IntegerToString(N) cat alphabet(k) cat "." cat IntegerToString(p) cat "." cat IntegerToString(sign));
end function;

//Here answer = {a} mod p^m
function one_teich(a,p,m)
	local answer,t;

	if p ne 2 then
		answer := a mod p;
		for k in [2..m] do
			t := 0;
			while (answer + t*p^(k-1))^(p-1) mod p^k ne 1 do
				t := t+1;
			end while;
			answer := answer + t*p^(k-1); 
		end for;
	else
		if a mod 4 eq 1 then
			return(1);
		else
			return(-1);
		end if;
	end if;

	return(answer);
end function;
			
//Here answer[a] = {a} (mod p^m) 
function list_of_teich(p,m)
	local answer;

	answer := [];
	if p ne 2 then
		for a in [1..p-1] do
			Append(~answer,one_teich(a,p,m));
		end for;
	else
		for a in [1,3] do
			Append(~answer,one_teich(a,p,m));
		end for;
	end if;		
	return(answer);
end function;

function teich(a,p,TEICH)
	if p ne 2 then
		return(TEICH[a]);
	else
		assert a in [1,3];
		if a eq 1 then
			return(TEICH[1]);
		else
			return(TEICH[2]);
		end if;
	end if;
end function;

function form_alpha(N,p,a_p,D,m)
	assert N mod p^2 ne 0;
	assert a_p mod p ne 0;
	assert D mod p ne 0;

	if N mod p eq 0 then
		return(a_p*KroneckerSymbol(D,p));
	else
		a_p := a_p * KroneckerSymbol(D,p);
		temp := a_p mod p;
		for k in [2..m] do
			t := (-(temp^2-a_p*temp+p) div p^(k-1) mod p) * InverseMod(2*temp - a_p,p);
			temp := temp + t*p^(k-1);
		end for;
		return(temp);
	end if;
end function;

function Pn(V,p,n,D)
//	print "Computing at level ",n;
	R<T>:=PolynomialRing(Rationals());
	TEICH := list_of_teich(p,n+1);
	answer:=0;
	oneplusTpower:=1;
	gammapower:=1;
	if p ne 2 then
		gamma:=1+p;
		if p ne 3 then
			for j in [0..p^n-1] do  // if j mod 1000 eq 0 then print(j); end if;
				for a in [1..((p-1) div 2)] do 
					ms := MS_twist(V,gammapower * teich(a,p,TEICH)/p^(n+1),D); 
					answer := (answer + ms * oneplusTpower) ;
				end for;
				gammapower := gammapower * gamma;
				oneplusTpower := (oneplusTpower * (1+T));
			end for;
		else 
			for j in [0..p^n-1] do
				// if j mod 1000 eq 0 then print j; end if;
				ms := MS_twist(V,gammapower/p^(n+1),D);
				answer := answer + ms * oneplusTpower;
				gammapower := gammapower * gamma;
				oneplusTpower := oneplusTpower * (1+T);
			end for;
		end if;
	else
		gamma:=5;
		for j in [0..p^n-1] do
			ms := 2*MS_twist(V,gammapower/p^(n+2),D);
			answer := answer + ms * oneplusTpower;
			gammapower := gammapower * gamma;
			oneplusTpower := oneplusTpower * (1+T);
		end for;
	end if;

	return(answer);
end function;

function cyclotomic(p,n)
	if n gt 0 then
		answer:=0;
		for a in [0..p-1] do
			answer:=answer+(1+T)^(a*p^(n-1));
		end for;
		return(answer);
	else
		return(T);
	end if;
end function;

function Lfunc_ord(V,p,a_p,n,D)
        N:=Level(V[1]);
        assert N mod p^2 ne 0;
	assert a_p mod p ne 0;
        alpha:=form_alpha(N,p,a_p,D,11);
        P:=Pn(V,p,n,D);
        Q:=Pn(V,p,n-1,D);
        L:=1/alpha^n*P-1/alpha^(n+1)*Q*cyclotomic(p,n);

	return L;
end function;

function valuation_of_coefficients(f,p)
	answer := [];
	f:=f+0*T;
	if f eq 0 then
		return([10000]);
 	else
		for j in [1..Degree(f)+1] do
			if Coefficient(f,j-1) eq 0 then
				answer := answer cat [10000];
			else
				answer := answer cat [Valuation(Coefficient(f,j-1),p)];
			end if;
		end for;
	end if;

	return(answer);
end function;
		
function iwasawa_invariants_poly(f,p)
	m,l:=Min(valuation_of_coefficients(f,p));
	return m,l-1;
end function;

//num only counts in ss case where it
//denotes + or -
function max_level(a_p,p,num)
	if a_p mod p ne 0 then
		if p eq 2 then
			return(6);  //(ok to 31)
		else
			if p eq 3 then
				return(4); //(ok to 26)
			else 
				if p eq 5 then
					return(3); //(ok to 24)
				else
					return(2); //(ok to p-1)
				end if;
			end if;
		end if;
	end if;
	if num eq 1 then
		if p eq 2 then
			return(5); //(ok to 85)
		else
			return(3); //(ok to 20 resp. 104)
		end if;
	else
		if p eq 2 then
			return(4); //(ok to 42)
		else
			if p eq 3 then	
				return(3); //(ok to 60)
			else 
				return(2); //(p=5 ok to 20)
			end if;
		end if;
	end if;
end function;

function split(N,p,a_p,D)
	return((N mod p eq 0) and (KroneckerSymbol(D,p)*a_p eq 1));
end function;

function twisted_sign(e,N,D)
	return(e*KroneckerSymbol(D,N)*KroneckerSymbol(D,-1));
end function;

function iwasawa_invariants_ord(V,p,a_p,D,e)
	local cycl,N,P,Q,r,L,m,l;
	cycl:=[];
	N:=Level(V[1]);
	assert N mod p^2 ne 0;
	alpha:=form_alpha(N,p,a_p,D,11); 
	P:=Pn(V,p,0,D); 
	Q:=Pn(V,p,1,D); 
	lmw:=0;
	r:=0;
	L:=Q-1/alpha*P*cyclotomic(p,1); 
	if (twisted_sign(e,N,D) eq -1) or (not split(N,p,a_p,D) and (L mod T eq 0)) then
		cycl:=cycl cat [0];
		lmw:=1;
		r:=1; 
	end if;
	m,l:=iwasawa_invariants_poly(L,p); 
	if (m eq 0) then
		if(#cycl gt 0) and (cycl[1] eq 0) and (twisted_sign(e,N,D) eq 1) then
			lmw := lmw+1;
			r := 2;
		end if;
		if split(N,p,a_p,D) then
			return(<D,r,m,l,lmw,"split">);
		else
			return(<D,r,m,l,lmw>);
		end if;
	else
		n:=2;  
		while (n le max_level(a_p,p,1)) do //print(n);
			P:=Q;
			Q:=Pn(V,p,n,D);
			if (N mod p eq 0) then
				L:=Q;
			else;
				L:=Q-1/alpha*P*cyclotomic(p,n);  
			end if; 
			if L mod cyclotomic(p,n-1) eq 0 then
				cycl:=cycl cat [n-1];
				lmw := lmw + p^(n-1) - p^(n-2);
			end if;
			m,l:=iwasawa_invariants_poly(L,p); //print([m,l]);
			if (m eq 0) then
				if (#cycl gt 0) and (cycl[1] eq 0) and (twisted_sign(e,N,D) eq 1) then
					lmw := lmw+1;
					r := 2;
				end if;
				if not split(N,p,a_p,D) then
					return(<D,r,m,l,lmw>);
				else
					return(<D,r,m,l,lmw,"split">);
				end if;
			else
				n:=n+1;
			end if;
		end while;
	end if;
	if not split(N,p,a_p,D) then
		return(<D,r,m,l,lmw,"non-zero mu">);
	else
		return(<D,r,m,l,lmw,"split","non-zero mu">);
	end if;
end function;

function q(n,p)
	answer:=0;
	while n gt 1 do
		answer:=answer + p^(n-1)-p^(n-2);
		n:=n-2;
	end while;
	return(answer);
end function;

function iwasawa_invariants_ss(V,p,a_p,D,e)
	cycl:=[];
	err:=false;
	N:=Level(V[1]);
	L:=Pn(V,p,0,D);
	r:=0;
	pluslmw:=0;
	minuslmw:=0;
	minusm:=0;
	minusl:=0;
	plusm:=0;
	plusl:=0;
	if L mod cyclotomic(p,0) eq 0 then
		cycl:=cycl cat [0];
		r:=1;
		pluslmw:=pluslmw+1;
		minuslmw:=minuslmw+1;
	end if;
	m,l:=iwasawa_invariants_poly(L,p); //print(0); print([m,l]);
	if (m eq 0) then  
		plusm:=0;
		plusl:=0;
		pluslmw:=0;
		minusm:=0;
		minusl:=0;
		minuslmw:=0;
	else
		n:=1;
		while (n le max_level(a_p,p,-1)) and (m ne 0) do 
			L:=Pn(V,p,2*n,D);
			if L mod cyclotomic(p,2*n) eq 0 then
				cycl:=cycl cat [2*n];
				minuslmw := minuslmw + p^(2*n) - p^(2*n-1);
			end if;			
			m,l:=iwasawa_invariants_poly(L,p); //print(n); print([m,l]);
			if (m eq 0) then
				minusm:=m;
				minusl:=l-q(2*n,p);
			end if;
			n:=n+1;
		end while;
		if (m ne 0) then
			err:=true;
			minusm:=m;
			minusl:=l-q(2*(n-1),p);
		end if;
	end if;

	L:=Pn(V,p,1,D);
	if L mod cyclotomic(p,1) eq 0 then
		cycl:=cycl cat [1];			
		pluslmw:=pluslmw+p-1;
	end if;			
	
	m,l:=iwasawa_invariants_poly(L,p); //print(1);print([m,l]);
	if (m eq 0) then  
		plusm:=m;
		plusl:=l;
	else
		n:=2;
		while (n le max_level(a_p,p,1)) and (m ne 0) do
			L:=Pn(V,p,2*n-1,D);
			if L mod cyclotomic(p,2*n-1) eq 0 then
				cycl:=cycl cat [2*n-1];
				pluslmw:=pluslmw+p^(2*n-1)-p^(2*n-2);
			end if;			
			m,l:=iwasawa_invariants_poly(L,p); //print(n);print([m,l]);
			if (m eq 0) then
				plusm:=m;
				plusl:=l-q(2*n-1,p);
			end if;
			n:=n+1;
		end while;
		if (m ne 0) then
			err:=true;
			plusm:=m;
			plusl:=l-q(2*n-3,p);
		end if;
	end if;

	if (r gt 0) and (twisted_sign(e,N,D) eq 1) then
		r:=2;
		pluslmw:=pluslmw+1;
		minuslmw:=minuslmw+1;
	end if;
	if p ne 2 then
		if not err then
			return(<D,r,plusm,plusl,pluslmw,minusm,minusl,minuslmw>);
		else
			return(<D,r,plusm,plusl,pluslmw,minusm,minusl,minuslmw,"error">);
		end if;
	else
		if not err then
			return(<D,r,minusm,minusl,minuslmw,plusm,plusl,pluslmw>);
			else
			return(<D,r,minusm,minusl,minuslmw,plusm,plusl,pluslmw,"error">);
		end if;
	end if;
end function;
	

//This function computes the Iwasawa invariants of the curve V twisted
//by the quadratic character of conductor D at the prime p
//
//Input:	V   - space of plus/minus modular symbols (produced by form_curve)
//		p   - prime
//		a_p - Fourier coefficient
//		D   - discriminant of qudratic twist
//		e   - sign of functional equation (e = 1 means plus; e = -1 means minus)
//
//Output:	if p is ordinary this function returns
//			<D,rank,mu,lambda,lambda_MW>
//		if p is ss it returns
//			<D,rank,mu^+,lambda^+,lambda_MW^+,mu^-,lambda^-,lambda_MW>
function iwasawa_invariants(V,p,a_p,D,e)
//	print [D,p];
	assert D mod p ne 0;
	assert Level(V[1]) mod p^2 ne 0;
	if e eq 0 then 
		if MSEven(V[1],0) eq 0 then
			e:=-1;
		else
			e:=1;
		end if;
	end if;
	if a_p mod p eq 0 then
		return(iwasawa_invariants_ss(V,p,a_p,D,e));
	else
		return(iwasawa_invariants_ord(V,p,a_p,D,e));
	end if;
end function;

function form_list_of_good_twists(N,p,bottom,top)
	answer:=[];
	for d in [bottom..top] do
		if (Gcd(d,p) eq 1) and (Gcd(d,N) eq 1) and (KroneckerSymbol(d,-1) eq 1) then
			if (d eq Squarefree(d)) and (d mod 4 eq 1) then
				answer := answer cat [d];
			end if;
			if (d mod 4 eq 0) and (d div 4 eq SquareFree(d div 4)) and ((d div 4) mod 4 ne 1) then
				answer := answer cat [d];
			end if;
		end if;
	end for;
	return(answer);
end function;

function form_list_of_good_twists_bound(N,p,bottom,top,start,sign)
	answer:=[];
	if sign eq -1 then
		return([]);
	else 
		for d in [start..top] do
			if (Gcd(d,p) eq 1) and (Gcd(d,N) eq 1) and (KroneckerSymbol(d,-1) eq 1) then
				if (d eq Squarefree(d)) and (d mod 4 eq 1) then
					answer := answer cat [d];
				end if;
				if (d mod 4 eq 0) and (d div 4 eq SquareFree(d div 4)) and ((d div 4) mod 4 ne 1) then
					answer := answer cat [d];
				end if;
			end if;
		end for;
		return(answer);
	end if;
end function;

function form_list_of_good_twists_neg(N,p,bottom,top)
	answer:=[];
	for e in [bottom..top] do
		d := -e;
		if (Gcd(d,p) eq 1) and (Gcd(d,N) eq 1) then
			if (d eq Squarefree(d)) and (d mod 4 eq 1) then
				answer := answer cat [d];
			end if;
			if (d mod 4 eq 0) and (d div 4 eq SquareFree(d div 4)) and ((d div 4) mod 4 ne 1) then
				answer := answer cat [d];
			end if;
		end if;
	end for;
	return(answer);
end function;

function form_list_of_good_twists_neg_bound(N,p,bottom,top,start,sign)
	answer:=[];
	if sign eq 1 then
		return(form_list_of_good_twists_neg(N,p,bottom,top));
	else
		for e in [start..top] do
			d := -e;
			if (Gcd(d,p) eq 1) and (Gcd(d,N) eq 1) then
				if (d eq Squarefree(d)) and (d mod 4 eq 1) then
					answer := answer cat [d];
				end if;
				if (d mod 4 eq 0) and (d div 4 eq SquareFree(d div 4)) and ((d div 4) mod 4 ne 1) then
					answer := answer cat [d];
				end if;
			end if;
		end for;
		return(answer);
	end if;
end function;

function display(data)
	if #data lt 8 then
		info:=[];
		for k in [1..5] do
			info := info cat [data[k]];
		end for;
		if #data eq 5 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ";
		end if;
		if #data eq 6 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],"|";
		end if;
		if #data eq 7 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],"|",data[7],"|";
		end if;
	else
		info:=[];
		for k in [1..8] do
			info := info cat [data[k]];
		end for;
		if #data eq 8 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],data[7],data[8]," |";
		end if;
		if #data eq 9 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],data[7],data[8]," |",data[9],"|";
		end if;
	end if;	
	return(1);
end function;

function display_single_curve(data,curve)
	if #data lt 8 then
		info:=[];
		for k in [1..5] do
			info := info cat [data[k]];
		end for;
		if #data eq 5 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ";
		end if;
		if #data eq 6 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],"|";
		end if;
		if #data eq 7 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],"|",data[7],"|";
		end if;
	else
		info:=[];
		for k in [1..8] do
			info := info cat [data[k]];
		end for;
		if #data eq 8 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],data[7],data[8]," |";
		end if;
		if #data eq 9 then
			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],data[7],data[8]," |",data[9],"|";
		end if;
	end if;	
	return(1);
end function;

function display_title(N,p,a_p,curve,sign)
	print curve;
	print "p =",p;
	if (N mod p eq 0) then
		print "Multiplicative";
	else
		if (a_p mod p ne 0) then
			print "Good ordinary ( a_p =",a_p,")";
		else
			print "Good supersingular ( a_p =",a_p,")";
		end if;
	end if;
	if sign eq 1 then
		print "Positive twists";
	else
		print "Negative twists";
	end if;
	print "";
	return(1);
end function;

function compute_invariants_of_twists(V,p,a_p,minD,maxD,sign,filename)
	assert Level(V[1]) mod p^2 ne 0;
	if MSEven(V[1],0) eq 0 then
		e := -1;
	else
		e := 1;
	end if;
	if sign gt 0 then
		twists:=form_list_of_good_twists(Level(V[1]),p,minD,maxD);
	else
		twists:=form_list_of_good_twists_neg(Level(V[1]),p,minD,maxD);
	end if;
	k:=1;
	for k in [1..#twists] do
		D:=twists[k];
		data:=iwasawa_invariants(V,p,a_p,D,e);
		SetOutputFile(filename);
		a:=display(data);
		UnsetOutputFile();
	end for;
	return(1);
end function;

function compute_invariants(V,p,a_p)
	assert Level(V[1]) mod p^2 ne 0;
	data:=iwasawa_invariants(V,p,a_p,1,0);
	a:=display_single_curve(data);
	return(1);
end function;
			
function run_thru_primes(V,Q,minp,maxp)
	if MSEven(V[1],0) eq 0 then
		e:=-1;
	else
		e:=1;
	end if;
	N:=Level(V[1]);
	for p in [minp..maxp] do
		if IsPrime(p) then
			if N mod p^2 ne 0 then
				a_p :=Numerator(Coefficient(Q,p));
				data:=iwasawa_invariants(V,p,a_p,1,e);
				if a_p mod p ne 0 then
					if split(N,p,a_p,1) then
						print p," | ",data[3],data[4],data[5],"*|";
					else
						print p," | ",data[3],data[4],data[5]," |";
					end if;
				else
					print  p," | ",data[3],data[4],data[5]," | ",data[6],data[7],data[8];
				end if;
			else
				print p,"-";
			end if;
		end if;
	end for;
	return(1);
end function;
	
//This function produces the plus and minus spaces of modular symbols 
//for the k-th curve of conductor N in Cremona's tables
//
//Input:	N - conductor
//		k - isogeny class
function form_curve(N,k)
	curve:=IntegerToString(N) cat alphabet(k); 
//	print "Forming curve ",curve;
	E:=EllipticCurve(C,curve);
	M1:=ModularSymbols(N,2,1);
	p:=2;
	while Dimension(M1) gt 1 do
		M1:=Kernel([<p,x-ap(E,p)>],M1);
		p:=NextPrime(p);
	end while;
	M2:=ModularSymbols(N,2,-1);
	p:=2;
	while Dimension(M2) gt 1 do
		M2:=Kernel([<p,x-ap(E,p)>],M2);
		p:=NextPrime(p);
	end while;
	V:=[M1,M2];
	return(V);
end function;

function run_thru_curves(minp,maxp,minN,maxN,minD,maxD)
	for N in [minN..maxN] do
		for k in [1..NumberOfIsogenyClasses(C,N)] do			
			curve:=IntegerToString(N) cat alphabet(k); 
			print "Forming curve ",curve;
			M:=ModularSymbols(curve);
			Q:=qEigenform(M,50);
			M1:=ModularSymbols(N,2,1);
			p:=2;
			while Dimension(M1) gt 1 do
				M1:=Kernel([<p,x-Coefficient(Q,p)>],M1);
				p:=NextPrime(p);
			end while;
			M2:=ModularSymbols(N,2,-1);
			p:=2;
			while Dimension(M2) gt 1 do
				M2:=Kernel([<p,x-Coefficient(Q,p)>],M2);
				p:=NextPrime(p);
			end while;
			V:=[M1,M2];
			for p in [minp..maxp] do
				if IsPrime(p) then
					if N mod p^2 ne 0 then
						a_p:=Numerator(Coefficient(Q,p));

						filename:=form_filename(N,k,p,1);
						curve:=IntegerToString(N) cat alphabet(k);
						SetOutputFile(filename);
						a:=display_title(N,p,a_p,curve,1);
						UnsetOutputFile();
						a:=compute_invariants_of_twists(V,p,a_p,minD,maxD,1,filename);

						filename:=form_filename(N,k,p,-1);
						SetOutputFile(filename);
						a:=display_title(N,p,a_p,curve,-1);
						UnsetOutputFile();
						a:=compute_invariants_of_twists(V,p,a_p,minD,maxD,-1,filename);
					end if;
				end if;
			end for;
			DisownChildren(M);
			DisownChildren(M1);
			DisownChildren(M2);
		end for;
	end for;
	return(1);
end function;

//startD is the first discriminant to start with on the first run
//and sign is the sign of the discriminant
function run_thru_twists(minp,maxp,minN,maxN,minD,maxD,startD,startsign)
	local p,N,k,M,Q,M1,V,L,e,sign,list1,list2,D,first;
	first:=1;
	for N in [minN..maxN] do
		for k in [1..NumberOfIsogenyClasses(C,N)] do			
			V:=form_curve(N,k);
			curve:=IntegerToString(N) cat alphabet(k); 
			E:=EllipticCurve(C,curve);
			L:=LRatio(V[1],1);
			if L eq 0 then
				e := -1;
			else
				e := 1;
			end if;
			done := false;
			for sign in [1,-1] do
			if first eq 1 then
				if sign gt 0 then
					twists:=form_list_of_good_twists_bound(Level(V[1]),-1,minD,maxD,startD,startsign);
				else
					twists:=form_list_of_good_twists_neg_bound(Level(V[1]),-1,minD,maxD,startD,startsign);
					first:=0;
				end if;
			else
				if sign gt 0 then
					twists:=form_list_of_good_twists(Level(V[1]),-1,minD,maxD);
				else
					twists:=form_list_of_good_twists_neg(Level(V[1]),-1,minD,maxD);
				end if;
			end if;
			for k in [1..#twists] do
				D:=twists[k];
					
				list1:=[];
				list2:=[];
				L:=MS_twist(V,0,D);
			for p in [minp..maxp] do
				if IsPrime(p) then
					a_p:=ap(E,p);
					if (N mod p^2 ne 0) and (D mod p ne 0) then
						if (N mod p eq 0) or (a_p*KroneckerSymbol(D,p) mod p eq 1) or 
						   (Valuation(L,p) ne 0) or (p eq 2) then
							done := true;
							data:=iwasawa_invariants(V,p,a_p,D,e); print(data);
							if ((#data eq 6) or (#data eq 7) or (#data eq 9)) and 
								(data[#data] eq "error") then
								print "error!";
							end if;
							if #data lt 8 then
								list1 := list1 cat [[data[3],data[4],data[5]]];
								if (#data gt 5) and (data[6] eq "split") then
									list1 := list1 cat [[-1,-1,-1]];
								end if;
								list2 := list2 cat [[-2,-2,-2]];
							else
								list1 := list1 cat [[data[3],data[4],data[5]]];
								list2 := list2 cat [[data[6],data[7],data[8]]];
							end if;
						else
							list1 := list1 cat [[0,0,0]];
							if a_p mod p eq 0 then
								list2 := list2 cat [[0,0,0]];
							else
								list2 := list2 cat [[-2,-2,-2]];
							end if;
						end if;
					else
						list1 := list1 cat [[-2,-2,-2]];
						list2 := list2 cat [[-2,-2,-2]];
						data := [0,0,0,0,0];
					end if;			
				end if;
			end for;
			if not done then
				data := [0,0];
			end if;
			line:=IntegerToString(D);
			if Abs(D) lt 10 then
				line := line cat "  ";
			else
				if Abs(D) lt 100 then
					line := line cat " ";
				end if;
			end if;
			line := line cat " | " cat IntegerToString(data[2]) cat " | ";
			for k in [1..#list1] do
				if list1[k][1] ge 0 then
					if ((k lt #list1) and (list1[k+1][1] ne -1)) or (k eq #list1) then
						if list1[k][2] lt 10 then
							line := line cat " ";
						end if;
						line := line cat IntegerToString(list1[k][1]) cat " "
								cat IntegerToString(list1[k][2]) cat " "
								  cat IntegerToString(list1[k][3]);
						if list1[k][3] lt 10 then
							line := line cat " ";
						end if;
						line := line cat " | ";
					else
						if list1[k][2] lt 10 then
							line := line cat " ";
						end if;
						line := line cat IntegerToString(list1[k][1]) cat " "
								cat IntegerToString(list1[k][2]) cat " "
								  cat IntegerToString(list1[k][3]);
						if list1[k][3] lt 10 then
							line := line cat " ";
						end if;
						line := line cat "*| ";
					end if;
				else
					if list1[k][1] eq -2 then
						line := line cat "        | ";
					end if;
				end if;
			end for;
			SetOutputFile(curve cat "." cat IntegerToString(sign));
			print line;
			line := "";
			if D gt 0 then
				line := "   ";
			else 
				line := "    ";
			end if;
			line := line cat " |   | ";
			for k in [1..#list2] do
				if list2[k][1] ne -2 then
					if list1[k][2] lt 10 then
						line := line cat " ";
					end if;
					line := line cat IntegerToString(list2[k][1]) cat " "
							cat IntegerToString(list2[k][2]) cat " "
							  cat IntegerToString(list2[k][3]);
					if list1[k][3] lt 10 then
						line := line cat " ";
					end if;
					line := line cat " | ";
				else
					line := line cat "        | ";
				end if;
			end for;
			print line;				
			UnsetOutputFile();
			end for;
			end for;
			DisownChildren(V[1]);
			DisownChildren(V[2]);
		end for;
	end for;
	return(1);
end function;


function run_thru_cremona(minp,maxp,minN,maxN)
	SetOutputFile("curves1-5000");
	for N in [minN..maxN] do 
		for k in [1..NumberOfIsogenyClasses(C,N)] do			
			curve:=IntegerToString(N) cat alphabet(k); 
			V:=form_curve(N,k);
			E:=EllipticCurve(C,curve);
			list1:=[];
			list2:=[];
			L:=LRatio(V[1],1);
			if L eq 0 then
				e := -1;
			else
				e := 1;
			end if;
			done := false;
			for p in [minp..maxp] do 
				if IsPrime(p) then //print(p); print(ap(E,p));
					a_p:=ap(E,p);
					if N mod p^2 ne 0 then
						if (N mod p eq 0) or (a_p mod p eq 1) or (Valuation(L,p) ne 0) or (p eq 2) then
							done := true;
							data:=iwasawa_invariants(V,p,a_p,1,e); //print(data);
							if ((#data eq 6) or (#data eq 7) or (#data eq 9)) and 
								(data[#data] eq "error") then
								print "error!";
							end if;
							if #data lt 8 then
								list1 := list1 cat [[data[3],data[4],data[5]]];
								if (#data gt 5) and (data[6] eq "split") then
									list1 := list1 cat [[-1,-1,-1]];
								end if;
								list2 := list2 cat [[-2,-2,-2]];
							else
								list1 := list1 cat [[data[3],data[4],data[5]]];
								list2 := list2 cat [[data[6],data[7],data[8]]];
							end if;
						else
							list1 := list1 cat [[0,0,0]];
							if a_p mod p eq 0 then
								list2 := list2 cat [[0,0,0]];
							else
								list2 := list2 cat [[-2,-2,-2]];
							end if;
						end if;
					else
						list1 := list1 cat [[-2,-2,-2]];
						list2 := list2 cat [[-2,-2,-2]];
					end if;			
				end if;
			end for;
			if not done then
				data := [0,0];
			end if;
			line:=curve;
			for j in [1..5-#line] do
				line:= line cat " ";
			end for;
			line:=line cat " | " cat IntegerToString(data[2]) cat " | ";
			for k in [1..#list1] do
				if list1[k][1] ge 0 then
					if ((k lt #list1) and (list1[k+1][1] ne -1)) or (k eq #list1) then
						if list1[k][2] lt 10 then
							line := line cat " ";
						end if;
						line := line cat IntegerToString(list1[k][1]) cat " "
								cat IntegerToString(list1[k][2]) cat " "
								  cat IntegerToString(list1[k][3]);
						if list1[k][3] lt 10 then
							line := line cat " ";
						end if;
						line := line cat " | ";
					else
						if list1[k][2] lt 10 then
							line := line cat " ";
						end if;
						line := line cat IntegerToString(list1[k][1]) cat " "
								cat IntegerToString(list1[k][2]) cat " "
								  cat IntegerToString(list1[k][3]);
						if list1[k][3] lt 10 then
							line := line cat " ";
						end if;
						line := line cat "*| ";
					end if;
				else
					if list1[k][1] eq -2 then
						line := line cat "        | ";
					end if;
				end if;
			end for;
			print line;
			line := "     ";
			line := line cat " |   | ";
			for k in [1..#list2] do
				if list2[k][1] ne -2 then
					if list2[k][2] lt 10 then
						line := line cat " ";
					end if;
					line := line cat IntegerToString(list2[k][1]) cat " "
							cat IntegerToString(list2[k][2]) cat " "
							  cat IntegerToString(list2[k][3]);
					if list2[k][3] lt 10 then
						line := line cat " ";
					end if;
					line := line cat " | ";
				else
					line := line cat "        | ";
				end if;
			end for;
			print line;				
			UnsetOutputFile();
			SetOutputFile("curves1-5000");
			DisownChildren(V[1]);
		end for;
	end for;
	UnsetOutputFile();
	return(1);
end function;



function ranktwo(V,e,bottom,top,bottom_neg,top_neg)
	answer:=[];
	list:=form_list_of_good_twists(Level(V[1]),1,bottom,top);
	for k in [1..#list] do print(list[k]);
		if (twisted_sign(e,Level(V[1]),list[k]) eq 1) then 
			if (MS_twist(V,0,list[k]) eq 0) then
				answer := answer cat [list[k]];
			end if;
		end if;
	end for;
	answer_neg:=[];
	list:=form_list_of_good_twists_neg(Level(V[1]),1,bottom_neg,top_neg);
	for k in [1..#list] do print(list[k]);
		if (twisted_sign(e,Level(V[1]),list[k]) eq 1) then 
			if (MS_twist(V,0,list[k]) eq 0) then
				answer_neg := answer_neg cat [list[k]];
			end if;
		end if;
	end for;
	return(<answer,answer_neg>);
end function;

function frac_mod(r,p,n)

        a:=Numerator(r) mod p^n;
        b:=InverseMod(Denominator(r),p^n);
        return(a*b mod p^n);

end function;

//This function computes the value of f/g evaluated at zero.
//
//Input: 	V - space of modular symbols
//		p - prime
//		n - accuracy
//		r - order of vanishing of f and g
//		D - discriminant of twist
function ratio_constant_terms_ss(V,p,n,r,D)

        P:=Pn(V,p,2*n+1,D);
	for i in [1..2*n+1] do
		if i mod 2 eq 0 then
			P := P div cyclotomic(p,i);
		end if;
	end for;
        Q:=Pn(V,p,2*n,D);
	for i in [1..2*n] do
		if i mod 2 eq 1 then
			Q := Q div cyclotomic(p,i);
		end if;
	end for;

        a:=Coefficient(P,r);
	if n ne 0 then
	        b:=Coefficient(Q,r);
	else
		b:=Q;
	end if;

	print "Plus constant term:",a;
	print "Minus constant term:",b;
	print "Valuation of plus term:",Valuation(a,p);
	print "Valuation of minus term:",Valuation(b,p);
	done:=false;
	if (Valuation(a,p) eq Valuation(b,p)) and (n gt Valuation(a,p)) then
		print "The ratio a/b is:",frac_mod(a/b,p,n-Valuation(a,p)),"(mod ",p,"^",n-Valuation(a,p),")";
	else
		print "Accuracy too low.";
	end if;

	print "";

	if ((Valuation(a,p) gt Valuation(b,p)) and (Valuation(b,p) lt n)) or
	   ((Valuation(b,p) gt Valuation(a,p)) and (Valuation(a,p) lt n)) or
	   ((Valuation(a,p) eq Valuation(b,p)) and (Valuation(a,p) lt n) and (frac_mod(a/b,p,n-Valuation(a,p)) ne ((p-1) div 2))) then
		return [*true,a,b,Valuation(a,p),Valuation(b,p),n*];
	else
		return [*false,a,b,Valuation(a,p),Valuation(b,p),n*];
	end if;

end function;

function smooth_valuations(C)
	for j in [2..#C] do
		if C[j] gt C[j-1] then
			C[j]:=C[j-1];
		end if;
	end for;

	return C;
end function;	

function accurate_coefficients(f,p,n)
	w:=T;
	for i in [1..n] do
		if (n-i) mod 2 eq 0 then 
			w := w * cyclotomic(p,i);
		end if;
	end for;
	D:=smooth_valuations(valuation_of_coefficients(w,p));
	E:=valuation_of_coefficients(f,p);
	C:=[(j lt #D) and (E[j] lt D[j]) : j in [1..#E]];

	if E[1] gt 1000 then
		C[1]:=true;
	end if;

	if (#E gt 1) and (E[2] gt 1000) then
		C[2]:=true;
	end if;

	return C;
			
end function;

function accuracy_of_coefficient(f,p,n,r)
	w:=1;
	for i in [1..n] do
		if (n-i) mod 2 eq 0 then 
			w := w * cyclotomic(p,i);
		end if;
	end for;
	D:=smooth_valuations(valuation_of_coefficients(w,p));

	return D[r+1];
end function;

function newton_polygon(C)
        local k,here,next,min_slope,answer,okay,prev_slope;

        answer:= [];
        okay:=true;

        here:=1;
        next:=1;

        while here lt #C do

                min_slope:=100000;

                for k in [here+1..#C] do
                        if (C[k]-C[here])*(1.0)/(k-here) le 1.0*min_slope then
                                min_slope := (C[k]-C[here])/(k-here);
                                next := k;
                        end if;
                end for;

                answer := answer cat [[next-here,-min_slope]];
                here := next;

        end while;

	C:=answer;
	return [C[j] : j in [1..#C] | C[j][2] gt 0];

end function;


function newton(f,p)
	C:=newton_polygon(valuation_of_coefficients(f,p));
	return [C[j] : j in [1..#C] | C[j][2] gt 0];
end function;

function break_points(C)
	count:=1;
	D:=[count];
	for j in [1..#C] do
		count:=count + Integers()!C[j][1];
		D:=D cat [count];
	end for;

	return D;
end function;

function is_newton_accurate(f,p,n);
	if Min(valuation_of_coefficients(f,p)) eq 0 then
		w:=T;
		for i in [1..n] do
			if (n-i) mod 2 eq 0 then 
				w := w * cyclotomic(p,i);
			end if;
		end for;
		D:=smooth_valuations(valuation_of_coefficients(w,p));
		C:=valuation_of_coefficients(f,p);
		for j in [1..#C] do
			if (C[j] gt D[j]) and (C[j] lt 100) then
				C[j]:=D[j];
			end if;
		end for;

		B:=break_points(newton_polygon(C));
		acc:=accurate_coefficients(f,p,n);
		correct:=true;
		for j in [1..#B] do
			correct := correct and acc[B[j]];
			if not correct then
				print "         not enough accuracy at coefficient",B[j];
			end if;
		end for;

		if correct then
			print "         Done!";
		end if;
	else
		print "         Mu is positive.";
		correct:=false;
	end if;

	return correct;
end function;

function signed_newton(V,p,D,sgn,e)
	e:=twisted_sign(e,Level(V[1]),D);

	done:=false;

	if sgn eq 1 then 
		del:=1;
	else
		del:=0;
	end if;	

	j:=del;
	while not done do
		print "Trying j =",j;
		P:=Pn(V,p,j,D);
		if P ne 0 then
			for i in [1..j] do
				if (i mod 2) eq ((del+1) mod 2) then
					P := P div cyclotomic(p,i);
				end if;
			end for;
			if (Coefficient(P,0) eq 0) and (e eq 1) then
				P:=P-Coefficient(P,1)*T;
			end if;
			print newton(P,p);
			done:=is_newton_accurate(P,p,j);
		else 
			print "        zero";
		end if;
		j:=j+2;
	end while;
	

	return newton(P,p);
end function;

function pm_newton(V,p,D,e)

	print "Plus newton..."; 
	print "";
	plus:=signed_newton(V,p,D,1,e);
	print plus;
	print "";
	print "";

	print "Minus newton..."; print "";
	minus:=signed_newton(V,p,D,-1,e);
	print minus;

	print "";
	print "";
	print "The plus newton polygon is:";
	print plus;
	print "";
	print "The minus newton polygon is:";
	print minus;

	return 1;
end function;

function newton_pm(V,p,D,e)
	e:=twisted_sign(e,Level(V[1]),D);

	done_plus_np:=false;
	done_minus_np:=false;

	j:=0;
	while (not done_plus_np) or (not done_minus_np) do 
		print "";
		print "";
		print "Trying j =",j,",",j+1;
		print "   For j =",j;
		if (not done_minus_np) then	
			print "      Minus newton polygon equals:";
			P:=Pn(V,p,j,D);
			if P ne 0 then
				for i in [1..j] do
					if (i mod 2) eq 1 then
						P := P div cyclotomic(p,i);
					end if;
				end for;
				if (j ne 0) and (Coefficient(P,0) eq 0) and (e eq 1) then
					P:=P-Coefficient(P,1)*T;
				end if;
				minus_np:=newton(P,p);
				print minus_np;
				done_minus_np:=is_newton_accurate(P,p,j);
			else 
				print "        zero";
			end if;
		else
			print "      Minus newton polygon already computed.";
		end if;


		print "   For j =",j+1;
		j:=j+1;
		if (not done_plus_np) then	
			print "      Plus newton polygon equals:";
			P:=Pn(V,p,j,D);
			if P ne 0 then
				for i in [1..j] do
					if (i mod 2) eq 0 then
						P := P div cyclotomic(p,i);
					end if;
				end for;
				if (Coefficient(P,0) eq 0) and (e eq 1) then
					P:=P-Coefficient(P,1)*T;
				end if;
				plus_np:=newton(P,p);
				print plus_np;
				done_plus_np:=is_newton_accurate(P,p,j);
			else 
				print "        zero";
			end if;
		else
			print "      Plus newton polygon already computed.";
		end if;
						
		j:=j+1;
	end while;

	print "";
	print "";
	print "";
	
	print "Plus newton:";
	print plus_np;

	print "Minus newton:";
	print minus_np;

	return "Done!";
end function;

function do_it_all(V,p,list,e)

	for k in [1..#list] do
		D:=list[k][1];
		r:=list[k][2];
		print "WORKING ON D=",D;
		done:=false;
		n:=1;
		info:=ratio_constant_terms_ss(V,p,n,r,D);
		while (not info[1]) and (n lt 4) do
			n:=n+1;
			info:=ratio_constant_terms_ss(V,p,n,r,D);
		end while;
		print info;
	end for;

	return "Done!";

end function;