primes := [ p : p in [1..10000] | IsPrime(p) ];
print "primes = list of primes less than 10000";

/*** The product of all Tamagawa Numbers ***/

prodtamagawa := function(E)
			SOL := 1;
			tama := TamagawaNumbers(E);
			for i:=1 to #tama do
				SOL := SOL*tama[i];
			end for;
			return SOL;
end function;

/*** AnalyticSha returns the conjectural value of Sha, using the BS-D Conjecture ***/

AnalyticSha := function(E : points := [[0, 1 ,0]], preci := 10, quiet := 1)
			M := MinimalModel(E);
			points := [E!P : P in points];
			O:=E![0,1,0];
			if points eq [O] then 
					reg := Regulator(M); 
				else 
					reg := Regulator(points);
			end if;
			disc := Discriminant(M);
			if disc lt 0 then 
					c:=1; 
				else 
					c:=2;
			end if;
			om := RealPeriod(M)*c;
			tama := prodtamagawa(M);
			anr, lval :=AnalyticRank(M : Precision:= preci);
			tors := (#TorsionSubgroup(M))^2; 
			sha := lval*tors/(om*reg*tama);
			if quiet eq 0 then
					print E;
					print "(# Torsion Subgroup of E)^2 = ", tors;
					print "Regulator = ", reg;
					print "Omega = ", om;
					print "Tamagawa numbers = ", tama;
					print "Analytic rank = ", anr;
					print "Residue of L at s=1 = ", lval;
					print "Conjectural value of Sha = ", sha;
			end if;
			return sha, Round(sha);
end function;