attachment:pLseries.m of days22/wuthrich

   1 load "newton.magma";
   2 
   3 R<x>:=PolynomialRing(Rationals());
   4 R<T>:=PolynomialRing(Rationals());
   5 C:=CremonaDatabase();
   6 
   7 function ap(E,p)
   8 	F:=FiniteField(p);
   9 	return(Trace(E,F));
  10 end function;
  11 
  12 function MSEven(M,r)
  13  	return(ModularSymbolEven(M,r)[1]/4);
  14 end function;
  15 
  16 function MSOdd(M,r)
  17  	return(ModularSymbolOdd(M,r)[1]/4);
  18 end function;
  19 
  20 function MS_twist(M,r,D)
  21 	answer:=0;
  22 	for a in [0..Abs(D)-1] do
  23 		if KroneckerSymbol(D,a) ne 0 then
  24 			if D ge 0 then
  25 				ms:=MSEven(M[1],r-a/D);
  26 			else
  27 				ms:=MSOdd(M[2],r+a/D);
  28 			end if;
  29 			answer:=answer+ms*KroneckerSymbol(D,a);
  30 		end if;
  31 	end for;
  32 
  33 	return(answer);
  34 end function;
  35 
  36 Cremona:=CremonaDatabase();
  37 
  38 function form_curve_ap(N,a_2,a_3,a_5,a_7)
  39 	R<x>:=PolynomialRing(Rationals());
  40 	M1:=ModularSymbols(N,2,+1);
  41 	M2:=ModularSymbols(N,2,-1);
  42 	I:=[<2,x-a_2>,<3,x-a_3>,<5,x-a_5>,<7,x-a_7>]; 
  43 	V1:=Kernel(I,M1);
  44 	V2:=Kernel(I,M2);
  45 	MSEven(V1,1/2);
  46 	ModularSymbolOdd(V2,1/2);
  47 
  48 	return([V1,V2]);
  49 end function;
  50 
  51 function form_curve(N,k)
  52 	E:=EllipticCurve(Cremona,N,k,1);
  53 	a_2:=ap(E,2);
  54 	a_3:=ap(E,3);
  55 	a_5:=ap(E,5);
  56 	a_7:=ap(E,7);
  57 	
  58 	return form_curve_ap(N,a_2,a_3,a_5,a_7);
  59 end function;
  60 
  61 
  62 function alphabet(k)
  63 	case k:
  64 		when 1: return("A");
  65 		when 2: return("B");
  66 		when 3: return("C");
  67 		when 4: return("D");
  68 		when 5: return("E");
  69 		when 6: return("F");
  70 		when 7: return("G");
  71 		when 8: return("H");
  72 		when 9: return("I");
  73 		when 10: return("J");
  74 		when 11: return("K");
  75 		when 12: return("L");
  76 		when 13: return("M");
  77 		when 14: return("N");
  78 		when 15: return("O");
  79 		when 16: return("P");
  80 		when 17: return("Q");
  81 		when 18: return("R");
  82 		when 19: return("S");
  83 		when 20: return("T");
  84 		when 21: return("U");
  85 		when 22: return("V");
  86 		when 23: return("W");
  87 		when 24: return("X");
  88 		when 25: return("Y");
  89 		when 26: return("Z");
  90 	end case;
  91 end function;
  92 
  93 function form_filename(N,k,p,sign)
  94   return(IntegerToString(N) cat alphabet(k) cat "." cat IntegerToString(p) cat "." cat IntegerToString(sign));
  95 end function;
  96 
  97 //Here answer = {a} mod p^m
  98 function one_teich(a,p,m)
  99 	local answer,t;
 100 
 101 	if p ne 2 then
 102 		answer := a mod p;
 103 		for k in [2..m] do
 104 			t := 0;
 105 			while (answer + t*p^(k-1))^(p-1) mod p^k ne 1 do
 106 				t := t+1;
 107 			end while;
 108 			answer := answer + t*p^(k-1); 
 109 		end for;
 110 	else
 111 		if a mod 4 eq 1 then
 112 			return(1);
 113 		else
 114 			return(-1);
 115 		end if;
 116 	end if;
 117 
 118 	return(answer);
 119 end function;
 120 			
 121 //Here answer[a] = {a} (mod p^m) 
 122 function list_of_teich(p,m)
 123 	local answer;
 124 
 125 	answer := [];
 126 	if p ne 2 then
 127 		for a in [1..p-1] do
 128 			Append(~answer,one_teich(a,p,m));
 129 		end for;
 130 	else
 131 		for a in [1,3] do
 132 			Append(~answer,one_teich(a,p,m));
 133 		end for;
 134 	end if;		
 135 	return(answer);
 136 end function;
 137 
 138 function teich(a,p,TEICH)
 139 	if p ne 2 then
 140 		return(TEICH[a]);
 141 	else
 142 		assert a in [1,3];
 143 		if a eq 1 then
 144 			return(TEICH[1]);
 145 		else
 146 			return(TEICH[2]);
 147 		end if;
 148 	end if;
 149 end function;
 150 
 151 function form_alpha(N,p,a_p,D,m)
 152 	assert N mod p^2 ne 0;
 153 	assert a_p mod p ne 0;
 154 	assert D mod p ne 0;
 155 
 156 	if N mod p eq 0 then
 157 		return(a_p*KroneckerSymbol(D,p));
 158 	else
 159 		a_p := a_p * KroneckerSymbol(D,p);
 160 		temp := a_p mod p;
 161 		for k in [2..m] do
 162 			t := (-(temp^2-a_p*temp+p) div p^(k-1) mod p) * InverseMod(2*temp - a_p,p);
 163 			temp := temp + t*p^(k-1);
 164 		end for;
 165 		return(temp);
 166 	end if;
 167 end function;
 168 
 169 function Pn(V,p,n,D)
 170 //	print "Computing at level ",n;
 171 	R<T>:=PolynomialRing(Rationals());
 172 	TEICH := list_of_teich(p,n+1);
 173 	answer:=0;
 174 	oneplusTpower:=1;
 175 	gammapower:=1;
 176 	if p ne 2 then
 177 		gamma:=1+p;
 178 		if p ne 3 then
 179 			for j in [0..p^n-1] do  // if j mod 1000 eq 0 then print(j); end if;
 180 				for a in [1..((p-1) div 2)] do 
 181 					ms := MS_twist(V,gammapower * teich(a,p,TEICH)/p^(n+1),D); 
 182 					answer := (answer + ms * oneplusTpower) ;
 183 				end for;
 184 				gammapower := gammapower * gamma;
 185 				oneplusTpower := (oneplusTpower * (1+T));
 186 			end for;
 187 		else 
 188 			for j in [0..p^n-1] do
 189 				// if j mod 1000 eq 0 then print j; end if;
 190 				ms := MS_twist(V,gammapower/p^(n+1),D);
 191 				answer := answer + ms * oneplusTpower;
 192 				gammapower := gammapower * gamma;
 193 				oneplusTpower := oneplusTpower * (1+T);
 194 			end for;
 195 		end if;
 196 	else
 197 		gamma:=5;
 198 		for j in [0..p^n-1] do
 199 			ms := 2*MS_twist(V,gammapower/p^(n+2),D);
 200 			answer := answer + ms * oneplusTpower;
 201 			gammapower := gammapower * gamma;
 202 			oneplusTpower := oneplusTpower * (1+T);
 203 		end for;
 204 	end if;
 205 
 206 	return(answer);
 207 end function;
 208 
 209 function cyclotomic(p,n)
 210 	if n gt 0 then
 211 		answer:=0;
 212 		for a in [0..p-1] do
 213 			answer:=answer+(1+T)^(a*p^(n-1));
 214 		end for;
 215 		return(answer);
 216 	else
 217 		return(T);
 218 	end if;
 219 end function;
 220 
 221 function Lfunc_ord(V,p,a_p,n,D)
 222         N:=Level(V[1]);
 223         assert N mod p^2 ne 0;
 224 	assert a_p mod p ne 0;
 225         alpha:=form_alpha(N,p,a_p,D,11);
 226         P:=Pn(V,p,n,D);
 227         Q:=Pn(V,p,n-1,D);
 228         L:=1/alpha^n*P-1/alpha^(n+1)*Q*cyclotomic(p,n);
 229 
 230 	return L;
 231 end function;
 232 
 233 function valuation_of_coefficients(f,p)
 234 	answer := [];
 235 	f:=f+0*T;
 236 	if f eq 0 then
 237 		return([10000]);
 238  	else
 239 		for j in [1..Degree(f)+1] do
 240 			if Coefficient(f,j-1) eq 0 then
 241 				answer := answer cat [10000];
 242 			else
 243 				answer := answer cat [Valuation(Coefficient(f,j-1),p)];
 244 			end if;
 245 		end for;
 246 	end if;
 247 
 248 	return(answer);
 249 end function;
 250 		
 251 function iwasawa_invariants_poly(f,p)
 252 	m,l:=Min(valuation_of_coefficients(f,p));
 253 	return m,l-1;
 254 end function;
 255 
 256 //num only counts in ss case where it
 257 //denotes + or -
 258 function max_level(a_p,p,num)
 259 	if a_p mod p ne 0 then
 260 		if p eq 2 then
 261 			return(6);  //(ok to 31)
 262 		else
 263 			if p eq 3 then
 264 				return(4); //(ok to 26)
 265 			else 
 266 				if p eq 5 then
 267 					return(3); //(ok to 24)
 268 				else
 269 					return(2); //(ok to p-1)
 270 				end if;
 271 			end if;
 272 		end if;
 273 	end if;
 274 	if num eq 1 then
 275 		if p eq 2 then
 276 			return(5); //(ok to 85)
 277 		else
 278 			return(3); //(ok to 20 resp. 104)
 279 		end if;
 280 	else
 281 		if p eq 2 then
 282 			return(4); //(ok to 42)
 283 		else
 284 			if p eq 3 then	
 285 				return(3); //(ok to 60)
 286 			else 
 287 				return(2); //(p=5 ok to 20)
 288 			end if;
 289 		end if;
 290 	end if;
 291 end function;
 292 
 293 function split(N,p,a_p,D)
 294 	return((N mod p eq 0) and (KroneckerSymbol(D,p)*a_p eq 1));
 295 end function;
 296 
 297 function twisted_sign(e,N,D)
 298 	return(e*KroneckerSymbol(D,N)*KroneckerSymbol(D,-1));
 299 end function;
 300 
 301 function iwasawa_invariants_ord(V,p,a_p,D,e)
 302 	local cycl,N,P,Q,r,L,m,l;
 303 	cycl:=[];
 304 	N:=Level(V[1]);
 305 	assert N mod p^2 ne 0;
 306 	alpha:=form_alpha(N,p,a_p,D,11); 
 307 	P:=Pn(V,p,0,D); 
 308 	Q:=Pn(V,p,1,D); 
 309 	lmw:=0;
 310 	r:=0;
 311 	L:=Q-1/alpha*P*cyclotomic(p,1); 
 312 	if (twisted_sign(e,N,D) eq -1) or (not split(N,p,a_p,D) and (L mod T eq 0)) then
 313 		cycl:=cycl cat [0];
 314 		lmw:=1;
 315 		r:=1; 
 316 	end if;
 317 	m,l:=iwasawa_invariants_poly(L,p); 
 318 	if (m eq 0) then
 319 		if(#cycl gt 0) and (cycl[1] eq 0) and (twisted_sign(e,N,D) eq 1) then
 320 			lmw := lmw+1;
 321 			r := 2;
 322 		end if;
 323 		if split(N,p,a_p,D) then
 324 			return(<D,r,m,l,lmw,"split">);
 325 		else
 326 			return(<D,r,m,l,lmw>);
 327 		end if;
 328 	else
 329 		n:=2;  
 330 		while (n le max_level(a_p,p,1)) do //print(n);
 331 			P:=Q;
 332 			Q:=Pn(V,p,n,D);
 333 			if (N mod p eq 0) then
 334 				L:=Q;
 335 			else;
 336 				L:=Q-1/alpha*P*cyclotomic(p,n);  
 337 			end if; 
 338 			if L mod cyclotomic(p,n-1) eq 0 then
 339 				cycl:=cycl cat [n-1];
 340 				lmw := lmw + p^(n-1) - p^(n-2);
 341 			end if;
 342 			m,l:=iwasawa_invariants_poly(L,p); //print([m,l]);
 343 			if (m eq 0) then
 344 				if (#cycl gt 0) and (cycl[1] eq 0) and (twisted_sign(e,N,D) eq 1) then
 345 					lmw := lmw+1;
 346 					r := 2;
 347 				end if;
 348 				if not split(N,p,a_p,D) then
 349 					return(<D,r,m,l,lmw>);
 350 				else
 351 					return(<D,r,m,l,lmw,"split">);
 352 				end if;
 353 			else
 354 				n:=n+1;
 355 			end if;
 356 		end while;
 357 	end if;
 358 	if not split(N,p,a_p,D) then
 359 		return(<D,r,m,l,lmw,"non-zero mu">);
 360 	else
 361 		return(<D,r,m,l,lmw,"split","non-zero mu">);
 362 	end if;
 363 end function;
 364 
 365 function q(n,p)
 366 	answer:=0;
 367 	while n gt 1 do
 368 		answer:=answer + p^(n-1)-p^(n-2);
 369 		n:=n-2;
 370 	end while;
 371 	return(answer);
 372 end function;
 373 
 374 function iwasawa_invariants_ss(V,p,a_p,D,e)
 375 	cycl:=[];
 376 	err:=false;
 377 	N:=Level(V[1]);
 378 	L:=Pn(V,p,0,D);
 379 	r:=0;
 380 	pluslmw:=0;
 381 	minuslmw:=0;
 382 	minusm:=0;
 383 	minusl:=0;
 384 	plusm:=0;
 385 	plusl:=0;
 386 	if L mod cyclotomic(p,0) eq 0 then
 387 		cycl:=cycl cat [0];
 388 		r:=1;
 389 		pluslmw:=pluslmw+1;
 390 		minuslmw:=minuslmw+1;
 391 	end if;
 392 	m,l:=iwasawa_invariants_poly(L,p); //print(0); print([m,l]);
 393 	if (m eq 0) then  
 394 		plusm:=0;
 395 		plusl:=0;
 396 		pluslmw:=0;
 397 		minusm:=0;
 398 		minusl:=0;
 399 		minuslmw:=0;
 400 	else
 401 		n:=1;
 402 		while (n le max_level(a_p,p,-1)) and (m ne 0) do 
 403 			L:=Pn(V,p,2*n,D);
 404 			if L mod cyclotomic(p,2*n) eq 0 then
 405 				cycl:=cycl cat [2*n];
 406 				minuslmw := minuslmw + p^(2*n) - p^(2*n-1);
 407 			end if;			
 408 			m,l:=iwasawa_invariants_poly(L,p); //print(n); print([m,l]);
 409 			if (m eq 0) then
 410 				minusm:=m;
 411 				minusl:=l-q(2*n,p);
 412 			end if;
 413 			n:=n+1;
 414 		end while;
 415 		if (m ne 0) then
 416 			err:=true;
 417 			minusm:=m;
 418 			minusl:=l-q(2*(n-1),p);
 419 		end if;
 420 	end if;
 421 
 422 	L:=Pn(V,p,1,D);
 423 	if L mod cyclotomic(p,1) eq 0 then
 424 		cycl:=cycl cat [1];			
 425 		pluslmw:=pluslmw+p-1;
 426 	end if;			
 427 	
 428 	m,l:=iwasawa_invariants_poly(L,p); //print(1);print([m,l]);
 429 	if (m eq 0) then  
 430 		plusm:=m;
 431 		plusl:=l;
 432 	else
 433 		n:=2;
 434 		while (n le max_level(a_p,p,1)) and (m ne 0) do
 435 			L:=Pn(V,p,2*n-1,D);
 436 			if L mod cyclotomic(p,2*n-1) eq 0 then
 437 				cycl:=cycl cat [2*n-1];
 438 				pluslmw:=pluslmw+p^(2*n-1)-p^(2*n-2);
 439 			end if;			
 440 			m,l:=iwasawa_invariants_poly(L,p); //print(n);print([m,l]);
 441 			if (m eq 0) then
 442 				plusm:=m;
 443 				plusl:=l-q(2*n-1,p);
 444 			end if;
 445 			n:=n+1;
 446 		end while;
 447 		if (m ne 0) then
 448 			err:=true;
 449 			plusm:=m;
 450 			plusl:=l-q(2*n-3,p);
 451 		end if;
 452 	end if;
 453 
 454 	if (r gt 0) and (twisted_sign(e,N,D) eq 1) then
 455 		r:=2;
 456 		pluslmw:=pluslmw+1;
 457 		minuslmw:=minuslmw+1;
 458 	end if;
 459 	if p ne 2 then
 460 		if not err then
 461 			return(<D,r,plusm,plusl,pluslmw,minusm,minusl,minuslmw>);
 462 		else
 463 			return(<D,r,plusm,plusl,pluslmw,minusm,minusl,minuslmw,"error">);
 464 		end if;
 465 	else
 466 		if not err then
 467 			return(<D,r,minusm,minusl,minuslmw,plusm,plusl,pluslmw>);
 468 			else
 469 			return(<D,r,minusm,minusl,minuslmw,plusm,plusl,pluslmw,"error">);
 470 		end if;
 471 	end if;
 472 end function;
 473 	
 474 
 475 //This function computes the Iwasawa invariants of the curve V twisted
 476 //by the quadratic character of conductor D at the prime p
 477 //
 478 //Input:	V   - space of plus/minus modular symbols (produced by form_curve)
 479 //		p   - prime
 480 //		a_p - Fourier coefficient
 481 //		D   - discriminant of qudratic twist
 482 //		e   - sign of functional equation (e = 1 means plus; e = -1 means minus)
 483 //
 484 //Output:	if p is ordinary this function returns
 485 //			<D,rank,mu,lambda,lambda_MW>
 486 //		if p is ss it returns
 487 //			<D,rank,mu^+,lambda^+,lambda_MW^+,mu^-,lambda^-,lambda_MW>
 488 function iwasawa_invariants(V,p,a_p,D,e)
 489 //	print [D,p];
 490 	assert D mod p ne 0;
 491 	assert Level(V[1]) mod p^2 ne 0;
 492 	if e eq 0 then 
 493 		if MSEven(V[1],0) eq 0 then
 494 			e:=-1;
 495 		else
 496 			e:=1;
 497 		end if;
 498 	end if;
 499 	if a_p mod p eq 0 then
 500 		return(iwasawa_invariants_ss(V,p,a_p,D,e));
 501 	else
 502 		return(iwasawa_invariants_ord(V,p,a_p,D,e));
 503 	end if;
 504 end function;
 505 
 506 function form_list_of_good_twists(N,p,bottom,top)
 507 	answer:=[];
 508 	for d in [bottom..top] do
 509 		if (Gcd(d,p) eq 1) and (Gcd(d,N) eq 1) and (KroneckerSymbol(d,-1) eq 1) then
 510 			if (d eq Squarefree(d)) and (d mod 4 eq 1) then
 511 				answer := answer cat [d];
 512 			end if;
 513 			if (d mod 4 eq 0) and (d div 4 eq SquareFree(d div 4)) and ((d div 4) mod 4 ne 1) then
 514 				answer := answer cat [d];
 515 			end if;
 516 		end if;
 517 	end for;
 518 	return(answer);
 519 end function;
 520 
 521 function form_list_of_good_twists_bound(N,p,bottom,top,start,sign)
 522 	answer:=[];
 523 	if sign eq -1 then
 524 		return([]);
 525 	else 
 526 		for d in [start..top] do
 527 			if (Gcd(d,p) eq 1) and (Gcd(d,N) eq 1) and (KroneckerSymbol(d,-1) eq 1) then
 528 				if (d eq Squarefree(d)) and (d mod 4 eq 1) then
 529 					answer := answer cat [d];
 530 				end if;
 531 				if (d mod 4 eq 0) and (d div 4 eq SquareFree(d div 4)) and ((d div 4) mod 4 ne 1) then
 532 					answer := answer cat [d];
 533 				end if;
 534 			end if;
 535 		end for;
 536 		return(answer);
 537 	end if;
 538 end function;
 539 
 540 function form_list_of_good_twists_neg(N,p,bottom,top)
 541 	answer:=[];
 542 	for e in [bottom..top] do
 543 		d := -e;
 544 		if (Gcd(d,p) eq 1) and (Gcd(d,N) eq 1) then
 545 			if (d eq Squarefree(d)) and (d mod 4 eq 1) then
 546 				answer := answer cat [d];
 547 			end if;
 548 			if (d mod 4 eq 0) and (d div 4 eq SquareFree(d div 4)) and ((d div 4) mod 4 ne 1) then
 549 				answer := answer cat [d];
 550 			end if;
 551 		end if;
 552 	end for;
 553 	return(answer);
 554 end function;
 555 
 556 function form_list_of_good_twists_neg_bound(N,p,bottom,top,start,sign)
 557 	answer:=[];
 558 	if sign eq 1 then
 559 		return(form_list_of_good_twists_neg(N,p,bottom,top));
 560 	else
 561 		for e in [start..top] do
 562 			d := -e;
 563 			if (Gcd(d,p) eq 1) and (Gcd(d,N) eq 1) then
 564 				if (d eq Squarefree(d)) and (d mod 4 eq 1) then
 565 					answer := answer cat [d];
 566 				end if;
 567 				if (d mod 4 eq 0) and (d div 4 eq SquareFree(d div 4)) and ((d div 4) mod 4 ne 1) then
 568 					answer := answer cat [d];
 569 				end if;
 570 			end if;
 571 		end for;
 572 		return(answer);
 573 	end if;
 574 end function;
 575 
 576 function display(data)
 577 	if #data lt 8 then
 578 		info:=[];
 579 		for k in [1..5] do
 580 			info := info cat [data[k]];
 581 		end for;
 582 		if #data eq 5 then
 583 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ";
 584 		end if;
 585 		if #data eq 6 then
 586 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],"|";
 587 		end if;
 588 		if #data eq 7 then
 589 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],"|",data[7],"|";
 590 		end if;
 591 	else
 592 		info:=[];
 593 		for k in [1..8] do
 594 			info := info cat [data[k]];
 595 		end for;
 596 		if #data eq 8 then
 597 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],data[7],data[8]," |";
 598 		end if;
 599 		if #data eq 9 then
 600 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],data[7],data[8]," |",data[9],"|";
 601 		end if;
 602 	end if;	
 603 	return(1);
 604 end function;
 605 
 606 function display_single_curve(data,curve)
 607 	if #data lt 8 then
 608 		info:=[];
 609 		for k in [1..5] do
 610 			info := info cat [data[k]];
 611 		end for;
 612 		if #data eq 5 then
 613 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ";
 614 		end if;
 615 		if #data eq 6 then
 616 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],"|";
 617 		end if;
 618 		if #data eq 7 then
 619 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],"|",data[7],"|";
 620 		end if;
 621 	else
 622 		info:=[];
 623 		for k in [1..8] do
 624 			info := info cat [data[k]];
 625 		end for;
 626 		if #data eq 8 then
 627 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],data[7],data[8]," |";
 628 		end if;
 629 		if #data eq 9 then
 630 			print data[1],"|",data[2],"| ",data[3],data[4],data[5]," | ",data[6],data[7],data[8]," |",data[9],"|";
 631 		end if;
 632 	end if;	
 633 	return(1);
 634 end function;
 635 
 636 function display_title(N,p,a_p,curve,sign)
 637 	print curve;
 638 	print "p =",p;
 639 	if (N mod p eq 0) then
 640 		print "Multiplicative";
 641 	else
 642 		if (a_p mod p ne 0) then
 643 			print "Good ordinary ( a_p =",a_p,")";
 644 		else
 645 			print "Good supersingular ( a_p =",a_p,")";
 646 		end if;
 647 	end if;
 648 	if sign eq 1 then
 649 		print "Positive twists";
 650 	else
 651 		print "Negative twists";
 652 	end if;
 653 	print "";
 654 	return(1);
 655 end function;
 656 
 657 function compute_invariants_of_twists(V,p,a_p,minD,maxD,sign,filename)
 658 	assert Level(V[1]) mod p^2 ne 0;
 659 	if MSEven(V[1],0) eq 0 then
 660 		e := -1;
 661 	else
 662 		e := 1;
 663 	end if;
 664 	if sign gt 0 then
 665 		twists:=form_list_of_good_twists(Level(V[1]),p,minD,maxD);
 666 	else
 667 		twists:=form_list_of_good_twists_neg(Level(V[1]),p,minD,maxD);
 668 	end if;
 669 	k:=1;
 670 	for k in [1..#twists] do
 671 		D:=twists[k];
 672 		data:=iwasawa_invariants(V,p,a_p,D,e);
 673 		SetOutputFile(filename);
 674 		a:=display(data);
 675 		UnsetOutputFile();
 676 	end for;
 677 	return(1);
 678 end function;
 679 
 680 function compute_invariants(V,p,a_p)
 681 	assert Level(V[1]) mod p^2 ne 0;
 682 	data:=iwasawa_invariants(V,p,a_p,1,0);
 683 	a:=display_single_curve(data);
 684 	return(1);
 685 end function;
 686 			
 687 function run_thru_primes(V,Q,minp,maxp)
 688 	if MSEven(V[1],0) eq 0 then
 689 		e:=-1;
 690 	else
 691 		e:=1;
 692 	end if;
 693 	N:=Level(V[1]);
 694 	for p in [minp..maxp] do
 695 		if IsPrime(p) then
 696 			if N mod p^2 ne 0 then
 697 				a_p :=Numerator(Coefficient(Q,p));
 698 				data:=iwasawa_invariants(V,p,a_p,1,e);
 699 				if a_p mod p ne 0 then
 700 					if split(N,p,a_p,1) then
 701 						print p," | ",data[3],data[4],data[5],"*|";
 702 					else
 703 						print p," | ",data[3],data[4],data[5]," |";
 704 					end if;
 705 				else
 706 					print  p," | ",data[3],data[4],data[5]," | ",data[6],data[7],data[8];
 707 				end if;
 708 			else
 709 				print p,"-";
 710 			end if;
 711 		end if;
 712 	end for;
 713 	return(1);
 714 end function;
 715 	
 716 //This function produces the plus and minus spaces of modular symbols 
 717 //for the k-th curve of conductor N in Cremona's tables
 718 //
 719 //Input:	N - conductor
 720 //		k - isogeny class
 721 function form_curve(N,k)
 722 	curve:=IntegerToString(N) cat alphabet(k); 
 723 //	print "Forming curve ",curve;
 724 	E:=EllipticCurve(C,curve);
 725 	M1:=ModularSymbols(N,2,1);
 726 	p:=2;
 727 	while Dimension(M1) gt 1 do
 728 		M1:=Kernel([<p,x-ap(E,p)>],M1);
 729 		p:=NextPrime(p);
 730 	end while;
 731 	M2:=ModularSymbols(N,2,-1);
 732 	p:=2;
 733 	while Dimension(M2) gt 1 do
 734 		M2:=Kernel([<p,x-ap(E,p)>],M2);
 735 		p:=NextPrime(p);
 736 	end while;
 737 	V:=[M1,M2];
 738 	return(V);
 739 end function;
 740 
 741 function run_thru_curves(minp,maxp,minN,maxN,minD,maxD)
 742 	for N in [minN..maxN] do
 743 		for k in [1..NumberOfIsogenyClasses(C,N)] do			
 744 			curve:=IntegerToString(N) cat alphabet(k); 
 745 			print "Forming curve ",curve;
 746 			M:=ModularSymbols(curve);
 747 			Q:=qEigenform(M,50);
 748 			M1:=ModularSymbols(N,2,1);
 749 			p:=2;
 750 			while Dimension(M1) gt 1 do
 751 				M1:=Kernel([<p,x-Coefficient(Q,p)>],M1);
 752 				p:=NextPrime(p);
 753 			end while;
 754 			M2:=ModularSymbols(N,2,-1);
 755 			p:=2;
 756 			while Dimension(M2) gt 1 do
 757 				M2:=Kernel([<p,x-Coefficient(Q,p)>],M2);
 758 				p:=NextPrime(p);
 759 			end while;
 760 			V:=[M1,M2];
 761 			for p in [minp..maxp] do
 762 				if IsPrime(p) then
 763 					if N mod p^2 ne 0 then
 764 						a_p:=Numerator(Coefficient(Q,p));
 765 
 766 						filename:=form_filename(N,k,p,1);
 767 						curve:=IntegerToString(N) cat alphabet(k);
 768 						SetOutputFile(filename);
 769 						a:=display_title(N,p,a_p,curve,1);
 770 						UnsetOutputFile();
 771 						a:=compute_invariants_of_twists(V,p,a_p,minD,maxD,1,filename);
 772 
 773 						filename:=form_filename(N,k,p,-1);
 774 						SetOutputFile(filename);
 775 						a:=display_title(N,p,a_p,curve,-1);
 776 						UnsetOutputFile();
 777 						a:=compute_invariants_of_twists(V,p,a_p,minD,maxD,-1,filename);
 778 					end if;
 779 				end if;
 780 			end for;
 781 			DisownChildren(M);
 782 			DisownChildren(M1);
 783 			DisownChildren(M2);
 784 		end for;
 785 	end for;
 786 	return(1);
 787 end function;
 788 
 789 //startD is the first discriminant to start with on the first run
 790 //and sign is the sign of the discriminant
 791 function run_thru_twists(minp,maxp,minN,maxN,minD,maxD,startD,startsign)
 792 	local p,N,k,M,Q,M1,V,L,e,sign,list1,list2,D,first;
 793 	first:=1;
 794 	for N in [minN..maxN] do
 795 		for k in [1..NumberOfIsogenyClasses(C,N)] do			
 796 			V:=form_curve(N,k);
 797 			curve:=IntegerToString(N) cat alphabet(k); 
 798 			E:=EllipticCurve(C,curve);
 799 			L:=LRatio(V[1],1);
 800 			if L eq 0 then
 801 				e := -1;
 802 			else
 803 				e := 1;
 804 			end if;
 805 			done := false;
 806 			for sign in [1,-1] do
 807 			if first eq 1 then
 808 				if sign gt 0 then
 809 					twists:=form_list_of_good_twists_bound(Level(V[1]),-1,minD,maxD,startD,startsign);
 810 				else
 811 					twists:=form_list_of_good_twists_neg_bound(Level(V[1]),-1,minD,maxD,startD,startsign);
 812 					first:=0;
 813 				end if;
 814 			else
 815 				if sign gt 0 then
 816 					twists:=form_list_of_good_twists(Level(V[1]),-1,minD,maxD);
 817 				else
 818 					twists:=form_list_of_good_twists_neg(Level(V[1]),-1,minD,maxD);
 819 				end if;
 820 			end if;
 821 			for k in [1..#twists] do
 822 				D:=twists[k];
 823 					
 824 				list1:=[];
 825 				list2:=[];
 826 				L:=MS_twist(V,0,D);
 827 			for p in [minp..maxp] do
 828 				if IsPrime(p) then
 829 					a_p:=ap(E,p);
 830 					if (N mod p^2 ne 0) and (D mod p ne 0) then
 831 						if (N mod p eq 0) or (a_p*KroneckerSymbol(D,p) mod p eq 1) or 
 832 						   (Valuation(L,p) ne 0) or (p eq 2) then
 833 							done := true;
 834 							data:=iwasawa_invariants(V,p,a_p,D,e); print(data);
 835 							if ((#data eq 6) or (#data eq 7) or (#data eq 9)) and 
 836 								(data[#data] eq "error") then
 837 								print "error!";
 838 							end if;
 839 							if #data lt 8 then
 840 								list1 := list1 cat [[data[3],data[4],data[5]]];
 841 								if (#data gt 5) and (data[6] eq "split") then
 842 									list1 := list1 cat [[-1,-1,-1]];
 843 								end if;
 844 								list2 := list2 cat [[-2,-2,-2]];
 845 							else
 846 								list1 := list1 cat [[data[3],data[4],data[5]]];
 847 								list2 := list2 cat [[data[6],data[7],data[8]]];
 848 							end if;
 849 						else
 850 							list1 := list1 cat [[0,0,0]];
 851 							if a_p mod p eq 0 then
 852 								list2 := list2 cat [[0,0,0]];
 853 							else
 854 								list2 := list2 cat [[-2,-2,-2]];
 855 							end if;
 856 						end if;
 857 					else
 858 						list1 := list1 cat [[-2,-2,-2]];
 859 						list2 := list2 cat [[-2,-2,-2]];
 860 						data := [0,0,0,0,0];
 861 					end if;			
 862 				end if;
 863 			end for;
 864 			if not done then
 865 				data := [0,0];
 866 			end if;
 867 			line:=IntegerToString(D);
 868 			if Abs(D) lt 10 then
 869 				line := line cat "  ";
 870 			else
 871 				if Abs(D) lt 100 then
 872 					line := line cat " ";
 873 				end if;
 874 			end if;
 875 			line := line cat " | " cat IntegerToString(data[2]) cat " | ";
 876 			for k in [1..#list1] do
 877 				if list1[k][1] ge 0 then
 878 					if ((k lt #list1) and (list1[k+1][1] ne -1)) or (k eq #list1) then
 879 						if list1[k][2] lt 10 then
 880 							line := line cat " ";
 881 						end if;
 882 						line := line cat IntegerToString(list1[k][1]) cat " "
 883 								cat IntegerToString(list1[k][2]) cat " "
 884 								  cat IntegerToString(list1[k][3]);
 885 						if list1[k][3] lt 10 then
 886 							line := line cat " ";
 887 						end if;
 888 						line := line cat " | ";
 889 					else
 890 						if list1[k][2] lt 10 then
 891 							line := line cat " ";
 892 						end if;
 893 						line := line cat IntegerToString(list1[k][1]) cat " "
 894 								cat IntegerToString(list1[k][2]) cat " "
 895 								  cat IntegerToString(list1[k][3]);
 896 						if list1[k][3] lt 10 then
 897 							line := line cat " ";
 898 						end if;
 899 						line := line cat "*| ";
 900 					end if;
 901 				else
 902 					if list1[k][1] eq -2 then
 903 						line := line cat "        | ";
 904 					end if;
 905 				end if;
 906 			end for;
 907 			SetOutputFile(curve cat "." cat IntegerToString(sign));
 908 			print line;
 909 			line := "";
 910 			if D gt 0 then
 911 				line := "   ";
 912 			else 
 913 				line := "    ";
 914 			end if;
 915 			line := line cat " |   | ";
 916 			for k in [1..#list2] do
 917 				if list2[k][1] ne -2 then
 918 					if list1[k][2] lt 10 then
 919 						line := line cat " ";
 920 					end if;
 921 					line := line cat IntegerToString(list2[k][1]) cat " "
 922 							cat IntegerToString(list2[k][2]) cat " "
 923 							  cat IntegerToString(list2[k][3]);
 924 					if list1[k][3] lt 10 then
 925 						line := line cat " ";
 926 					end if;
 927 					line := line cat " | ";
 928 				else
 929 					line := line cat "        | ";
 930 				end if;
 931 			end for;
 932 			print line;				
 933 			UnsetOutputFile();
 934 			end for;
 935 			end for;
 936 			DisownChildren(V[1]);
 937 			DisownChildren(V[2]);
 938 		end for;
 939 	end for;
 940 	return(1);
 941 end function;
 942 
 943 
 944 function run_thru_cremona(minp,maxp,minN,maxN)
 945 	SetOutputFile("curves1-5000");
 946 	for N in [minN..maxN] do 
 947 		for k in [1..NumberOfIsogenyClasses(C,N)] do			
 948 			curve:=IntegerToString(N) cat alphabet(k); 
 949 			V:=form_curve(N,k);
 950 			E:=EllipticCurve(C,curve);
 951 			list1:=[];
 952 			list2:=[];
 953 			L:=LRatio(V[1],1);
 954 			if L eq 0 then
 955 				e := -1;
 956 			else
 957 				e := 1;
 958 			end if;
 959 			done := false;
 960 			for p in [minp..maxp] do 
 961 				if IsPrime(p) then //print(p); print(ap(E,p));
 962 					a_p:=ap(E,p);
 963 					if N mod p^2 ne 0 then
 964 						if (N mod p eq 0) or (a_p mod p eq 1) or (Valuation(L,p) ne 0) or (p eq 2) then
 965 							done := true;
 966 							data:=iwasawa_invariants(V,p,a_p,1,e); //print(data);
 967 							if ((#data eq 6) or (#data eq 7) or (#data eq 9)) and 
 968 								(data[#data] eq "error") then
 969 								print "error!";
 970 							end if;
 971 							if #data lt 8 then
 972 								list1 := list1 cat [[data[3],data[4],data[5]]];
 973 								if (#data gt 5) and (data[6] eq "split") then
 974 									list1 := list1 cat [[-1,-1,-1]];
 975 								end if;
 976 								list2 := list2 cat [[-2,-2,-2]];
 977 							else
 978 								list1 := list1 cat [[data[3],data[4],data[5]]];
 979 								list2 := list2 cat [[data[6],data[7],data[8]]];
 980 							end if;
 981 						else
 982 							list1 := list1 cat [[0,0,0]];
 983 							if a_p mod p eq 0 then
 984 								list2 := list2 cat [[0,0,0]];
 985 							else
 986 								list2 := list2 cat [[-2,-2,-2]];
 987 							end if;
 988 						end if;
 989 					else
 990 						list1 := list1 cat [[-2,-2,-2]];
 991 						list2 := list2 cat [[-2,-2,-2]];
 992 					end if;			
 993 				end if;
 994 			end for;
 995 			if not done then
 996 				data := [0,0];
 997 			end if;
 998 			line:=curve;
 999 			for j in [1..5-#line] do
1000 				line:= line cat " ";
1001 			end for;
1002 			line:=line cat " | " cat IntegerToString(data[2]) cat " | ";
1003 			for k in [1..#list1] do
1004 				if list1[k][1] ge 0 then
1005 					if ((k lt #list1) and (list1[k+1][1] ne -1)) or (k eq #list1) then
1006 						if list1[k][2] lt 10 then
1007 							line := line cat " ";
1008 						end if;
1009 						line := line cat IntegerToString(list1[k][1]) cat " "
1010 								cat IntegerToString(list1[k][2]) cat " "
1011 								  cat IntegerToString(list1[k][3]);
1012 						if list1[k][3] lt 10 then
1013 							line := line cat " ";
1014 						end if;
1015 						line := line cat " | ";
1016 					else
1017 						if list1[k][2] lt 10 then
1018 							line := line cat " ";
1019 						end if;
1020 						line := line cat IntegerToString(list1[k][1]) cat " "
1021 								cat IntegerToString(list1[k][2]) cat " "
1022 								  cat IntegerToString(list1[k][3]);
1023 						if list1[k][3] lt 10 then
1024 							line := line cat " ";
1025 						end if;
1026 						line := line cat "*| ";
1027 					end if;
1028 				else
1029 					if list1[k][1] eq -2 then
1030 						line := line cat "        | ";
1031 					end if;
1032 				end if;
1033 			end for;
1034 			print line;
1035 			line := "     ";
1036 			line := line cat " |   | ";
1037 			for k in [1..#list2] do
1038 				if list2[k][1] ne -2 then
1039 					if list2[k][2] lt 10 then
1040 						line := line cat " ";
1041 					end if;
1042 					line := line cat IntegerToString(list2[k][1]) cat " "
1043 							cat IntegerToString(list2[k][2]) cat " "
1044 							  cat IntegerToString(list2[k][3]);
1045 					if list2[k][3] lt 10 then
1046 						line := line cat " ";
1047 					end if;
1048 					line := line cat " | ";
1049 				else
1050 					line := line cat "        | ";
1051 				end if;
1052 			end for;
1053 			print line;				
1054 			UnsetOutputFile();
1055 			SetOutputFile("curves1-5000");
1056 			DisownChildren(V[1]);
1057 		end for;
1058 	end for;
1059 	UnsetOutputFile();
1060 	return(1);
1061 end function;
1062 
1063 
1064 
1065 function ranktwo(V,e,bottom,top,bottom_neg,top_neg)
1066 	answer:=[];
1067 	list:=form_list_of_good_twists(Level(V[1]),1,bottom,top);
1068 	for k in [1..#list] do print(list[k]);
1069 		if (twisted_sign(e,Level(V[1]),list[k]) eq 1) then 
1070 			if (MS_twist(V,0,list[k]) eq 0) then
1071 				answer := answer cat [list[k]];
1072 			end if;
1073 		end if;
1074 	end for;
1075 	answer_neg:=[];
1076 	list:=form_list_of_good_twists_neg(Level(V[1]),1,bottom_neg,top_neg);
1077 	for k in [1..#list] do print(list[k]);
1078 		if (twisted_sign(e,Level(V[1]),list[k]) eq 1) then 
1079 			if (MS_twist(V,0,list[k]) eq 0) then
1080 				answer_neg := answer_neg cat [list[k]];
1081 			end if;
1082 		end if;
1083 	end for;
1084 	return(<answer,answer_neg>);
1085 end function;
1086 
1087 function frac_mod(r,p,n)
1088 
1089         a:=Numerator(r) mod p^n;
1090         b:=InverseMod(Denominator(r),p^n);
1091         return(a*b mod p^n);
1092 
1093 end function;
1094 
1095 //This function computes the value of f/g evaluated at zero.
1096 //
1097 //Input: 	V - space of modular symbols
1098 //		p - prime
1099 //		n - accuracy
1100 //		r - order of vanishing of f and g
1101 //		D - discriminant of twist
1102 function ratio_constant_terms_ss(V,p,n,r,D)
1103 
1104         P:=Pn(V,p,2*n+1,D);
1105 	for i in [1..2*n+1] do
1106 		if i mod 2 eq 0 then
1107 			P := P div cyclotomic(p,i);
1108 		end if;
1109 	end for;
1110         Q:=Pn(V,p,2*n,D);
1111 	for i in [1..2*n] do
1112 		if i mod 2 eq 1 then
1113 			Q := Q div cyclotomic(p,i);
1114 		end if;
1115 	end for;
1116 
1117         a:=Coefficient(P,r);
1118 	if n ne 0 then
1119 	        b:=Coefficient(Q,r);
1120 	else
1121 		b:=Q;
1122 	end if;
1123 
1124 	print "Plus constant term:",a;
1125 	print "Minus constant term:",b;
1126 	print "Valuation of plus term:",Valuation(a,p);
1127 	print "Valuation of minus term:",Valuation(b,p);
1128 	done:=false;
1129 	if (Valuation(a,p) eq Valuation(b,p)) and (n gt Valuation(a,p)) then
1130 		print "The ratio a/b is:",frac_mod(a/b,p,n-Valuation(a,p)),"(mod ",p,"^",n-Valuation(a,p),")";
1131 	else
1132 		print "Accuracy too low.";
1133 	end if;
1134 
1135 	print "";
1136 
1137 	if ((Valuation(a,p) gt Valuation(b,p)) and (Valuation(b,p) lt n)) or
1138 	   ((Valuation(b,p) gt Valuation(a,p)) and (Valuation(a,p) lt n)) or
1139 	   ((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
1140 		return [*true,a,b,Valuation(a,p),Valuation(b,p),n*];
1141 	else
1142 		return [*false,a,b,Valuation(a,p),Valuation(b,p),n*];
1143 	end if;
1144 
1145 end function;
1146 
1147 function smooth_valuations(C)
1148 	for j in [2..#C] do
1149 		if C[j] gt C[j-1] then
1150 			C[j]:=C[j-1];
1151 		end if;
1152 	end for;
1153 
1154 	return C;
1155 end function;	
1156 
1157 function accurate_coefficients(f,p,n)
1158 	w:=T;
1159 	for i in [1..n] do
1160 		if (n-i) mod 2 eq 0 then 
1161 			w := w * cyclotomic(p,i);
1162 		end if;
1163 	end for;
1164 	D:=smooth_valuations(valuation_of_coefficients(w,p));
1165 	E:=valuation_of_coefficients(f,p);
1166 	C:=[(j lt #D) and (E[j] lt D[j]) : j in [1..#E]];
1167 
1168 	if E[1] gt 1000 then
1169 		C[1]:=true;
1170 	end if;
1171 
1172 	if (#E gt 1) and (E[2] gt 1000) then
1173 		C[2]:=true;
1174 	end if;
1175 
1176 	return C;
1177 			
1178 end function;
1179 
1180 function accuracy_of_coefficient(f,p,n,r)
1181 	w:=1;
1182 	for i in [1..n] do
1183 		if (n-i) mod 2 eq 0 then 
1184 			w := w * cyclotomic(p,i);
1185 		end if;
1186 	end for;
1187 	D:=smooth_valuations(valuation_of_coefficients(w,p));
1188 
1189 	return D[r+1];
1190 end function;
1191 
1192 function newton_polygon(C)
1193         local k,here,next,min_slope,answer,okay,prev_slope;
1194 
1195         answer:= [];
1196         okay:=true;
1197 
1198         here:=1;
1199         next:=1;
1200 
1201         while here lt #C do
1202 
1203                 min_slope:=100000;
1204 
1205                 for k in [here+1..#C] do
1206                         if (C[k]-C[here])*(1.0)/(k-here) le 1.0*min_slope then
1207                                 min_slope := (C[k]-C[here])/(k-here);
1208                                 next := k;
1209                         end if;
1210                 end for;
1211 
1212                 answer := answer cat [[next-here,-min_slope]];
1213                 here := next;
1214 
1215         end while;
1216 
1217 	C:=answer;
1218 	return [C[j] : j in [1..#C] | C[j][2] gt 0];
1219 
1220 end function;
1221 
1222 
1223 function newton(f,p)
1224 	C:=newton_polygon(valuation_of_coefficients(f,p));
1225 	return [C[j] : j in [1..#C] | C[j][2] gt 0];
1226 end function;
1227 
1228 function break_points(C)
1229 	count:=1;
1230 	D:=[count];
1231 	for j in [1..#C] do
1232 		count:=count + Integers()!C[j][1];
1233 		D:=D cat [count];
1234 	end for;
1235 
1236 	return D;
1237 end function;
1238 
1239 function is_newton_accurate(f,p,n);
1240 	if Min(valuation_of_coefficients(f,p)) eq 0 then
1241 		w:=T;
1242 		for i in [1..n] do
1243 			if (n-i) mod 2 eq 0 then 
1244 				w := w * cyclotomic(p,i);
1245 			end if;
1246 		end for;
1247 		D:=smooth_valuations(valuation_of_coefficients(w,p));
1248 		C:=valuation_of_coefficients(f,p);
1249 		for j in [1..#C] do
1250 			if (C[j] gt D[j]) and (C[j] lt 100) then
1251 				C[j]:=D[j];
1252 			end if;
1253 		end for;
1254 
1255 		B:=break_points(newton_polygon(C));
1256 		acc:=accurate_coefficients(f,p,n);
1257 		correct:=true;
1258 		for j in [1..#B] do
1259 			correct := correct and acc[B[j]];
1260 			if not correct then
1261 				print "         not enough accuracy at coefficient",B[j];
1262 			end if;
1263 		end for;
1264 
1265 		if correct then
1266 			print "         Done!";
1267 		end if;
1268 	else
1269 		print "         Mu is positive.";
1270 		correct:=false;
1271 	end if;
1272 
1273 	return correct;
1274 end function;
1275 
1276 function signed_newton(V,p,D,sgn,e)
1277 	e:=twisted_sign(e,Level(V[1]),D);
1278 
1279 	done:=false;
1280 
1281 	if sgn eq 1 then 
1282 		del:=1;
1283 	else
1284 		del:=0;
1285 	end if;	
1286 
1287 	j:=del;
1288 	while not done do
1289 		print "Trying j =",j;
1290 		P:=Pn(V,p,j,D);
1291 		if P ne 0 then
1292 			for i in [1..j] do
1293 				if (i mod 2) eq ((del+1) mod 2) then
1294 					P := P div cyclotomic(p,i);
1295 				end if;
1296 			end for;
1297 			if (Coefficient(P,0) eq 0) and (e eq 1) then
1298 				P:=P-Coefficient(P,1)*T;
1299 			end if;
1300 			print newton(P,p);
1301 			done:=is_newton_accurate(P,p,j);
1302 		else 
1303 			print "        zero";
1304 		end if;
1305 		j:=j+2;
1306 	end while;
1307 	
1308 
1309 	return newton(P,p);
1310 end function;
1311 
1312 function pm_newton(V,p,D,e)
1313 
1314 	print "Plus newton..."; 
1315 	print "";
1316 	plus:=signed_newton(V,p,D,1,e);
1317 	print plus;
1318 	print "";
1319 	print "";
1320 
1321 	print "Minus newton..."; print "";
1322 	minus:=signed_newton(V,p,D,-1,e);
1323 	print minus;
1324 
1325 	print "";
1326 	print "";
1327 	print "The plus newton polygon is:";
1328 	print plus;
1329 	print "";
1330 	print "The minus newton polygon is:";
1331 	print minus;
1332 
1333 	return 1;
1334 end function;
1335 
1336 function newton_pm(V,p,D,e)
1337 	e:=twisted_sign(e,Level(V[1]),D);
1338 
1339 	done_plus_np:=false;
1340 	done_minus_np:=false;
1341 
1342 	j:=0;
1343 	while (not done_plus_np) or (not done_minus_np) do 
1344 		print "";
1345 		print "";
1346 		print "Trying j =",j,",",j+1;
1347 		print "   For j =",j;
1348 		if (not done_minus_np) then	
1349 			print "      Minus newton polygon equals:";
1350 			P:=Pn(V,p,j,D);
1351 			if P ne 0 then
1352 				for i in [1..j] do
1353 					if (i mod 2) eq 1 then
1354 						P := P div cyclotomic(p,i);
1355 					end if;
1356 				end for;
1357 				if (j ne 0) and (Coefficient(P,0) eq 0) and (e eq 1) then
1358 					P:=P-Coefficient(P,1)*T;
1359 				end if;
1360 				minus_np:=newton(P,p);
1361 				print minus_np;
1362 				done_minus_np:=is_newton_accurate(P,p,j);
1363 			else 
1364 				print "        zero";
1365 			end if;
1366 		else
1367 			print "      Minus newton polygon already computed.";
1368 		end if;
1369 
1370 
1371 		print "   For j =",j+1;
1372 		j:=j+1;
1373 		if (not done_plus_np) then	
1374 			print "      Plus newton polygon equals:";
1375 			P:=Pn(V,p,j,D);
1376 			if P ne 0 then
1377 				for i in [1..j] do
1378 					if (i mod 2) eq 0 then
1379 						P := P div cyclotomic(p,i);
1380 					end if;
1381 				end for;
1382 				if (Coefficient(P,0) eq 0) and (e eq 1) then
1383 					P:=P-Coefficient(P,1)*T;
1384 				end if;
1385 				plus_np:=newton(P,p);
1386 				print plus_np;
1387 				done_plus_np:=is_newton_accurate(P,p,j);
1388 			else 
1389 				print "        zero";
1390 			end if;
1391 		else
1392 			print "      Plus newton polygon already computed.";
1393 		end if;
1394 						
1395 		j:=j+1;
1396 	end while;
1397 
1398 	print "";
1399 	print "";
1400 	print "";
1401 	
1402 	print "Plus newton:";
1403 	print plus_np;
1404 
1405 	print "Minus newton:";
1406 	print minus_np;
1407 
1408 	return "Done!";
1409 end function;
1410 
1411 function do_it_all(V,p,list,e)
1412 
1413 	for k in [1..#list] do
1414 		D:=list[k][1];
1415 		r:=list[k][2];
1416 		print "WORKING ON D=",D;
1417 		done:=false;
1418 		n:=1;
1419 		info:=ratio_constant_terms_ss(V,p,n,r,D);
1420 		while (not info[1]) and (n lt 4) do
1421 			n:=n+1;
1422 			info:=ratio_constant_terms_ss(V,p,n,r,D);
1423 		end while;
1424 		print info;
1425 	end for;
1426 
1427 	return "Done!";
1428 
1429 end function;