powers := [[0, 0, 0, 3],
		   [0, 0, 1, 2],
		   [0, 0, 2, 1],
		   [0, 0, 3, 0],
	       [0, 1, 0, 2],
	       [0, 1, 1, 1],
		   [0, 1, 2, 0],
	       [0, 2, 0, 1],
	       [0, 2, 1, 0],
	       [0, 3, 0, 0],
	       [1, 0, 0, 2],
		   [1, 0, 1, 1],
		   [1, 0, 2, 0],
		   [1, 1, 0, 1],
		   [1, 1, 1, 0],
		   [1, 2, 0, 0],
		   [2, 0, 0, 1],
		   [2, 0, 1, 0],
		   [2, 1, 0, 0],
		   [3, 0, 0, 0]];




function randomSmoothCubicSurface(N,p)
	
	_<x,y,z,w> := PolynomialRing(Rationals(),4);
	P3<x,y,z,w> := ProjectiveSpace(Rationals(),3);
	
	listCoef := [];
	for i in [0..N] do
		Append(~listCoef,p^i);
	end for;	
	
	bool := true; 
	while bool do
	
		randSurface := 0*x;		
		for part in powers do
			coef := Random(listCoef);
			monomial := x^part[1] * y^part[2] * z^part[3] * w^part[4];
			randSurface := randSurface + coef * monomial;
		end for;

		S := Surface(P3,randSurface);	
		bool := IsSingular(S);

	end while;

	return randSurface;

end function;



function howManyLines(surface,p)

	f := surface;
		
	QQabcd<a,b,c,d> := PolynomialRing(Rationals(),4);
	_<s,t> := PolynomialRing(QQabcd,2);

	params := [[a*s+b*t,c*s+d*t,s,t],
			   [a*s+b*t,s,c*t,t],
			   [a*s+b*t,s,t,0],
			   [s,a*t,b*t,t],
			   [s,a*t,t,0],
			   [s,t,0,0]];

	pt0 := [1,0];
	pt1 := [1,1];
	ptmin1 := [1,-1]; 
	ptinf := [0,1];

	polys := [];

	for par in params do
		
		g0 := Evaluate(f,[Evaluate(par[i],pt0) : i in [1..4]]);
		g1 := Evaluate(f,[Evaluate(par[i],pt1) : i in [1..4]]);
		gmin1 := Evaluate(f,[Evaluate(par[i],ptmin1) : i in [1..4]]);
		ginf := Evaluate(f,[Evaluate(par[i],ptinf) : i in [1..4]]);
	  
		I := Ideal([g0,g1,gmin1,ginf]);
		GB := GroebnerBasis(I);
		g := GB[#GB];
		Append(~polys,g);

	end for;

	K := pAdicField(p,750);

	nsol := 0;

	for g in polys do

		newcoeff := [];
		for i in [0..27] do
			Append(~newcoeff,K!Rationals()!Coefficient(g,4,i));
		end for;

		_<x> := PolynomialRing(K);
		gK := 0*x;
		for i in [1..28] do
			gK := gK + newcoeff[i]*x^(i-1);
		end for;

		try nsol := nsol + #Roots(gK);
		catch e return -1;
		end try;

	end for;

	return nsol;

end function;



function lineCounts(N,p,M)

	counts := [0 : j in [-1..27]]; 
	
	for i in [1..M] do
		print i;
		
		f := randomSmoothCubicSurface(N,p);
		count := howManyLines(f,p);
		counts[count+2] := counts[count+2] + 1; 
	end for;

	return counts;

end function;

///////////////////////////////////////////
////////////////// TESTS //////////////////
///////////////////////////////////////////

p := 7; // prime
M := 1000; // number of trials
N := 7; // precision

counts :=  lineCounts(N,p,M);

RR := RealField(5);
 
print("--------------------------------------------------------------------------------------------");
print "Measure: Tropical  |  Prime: ",p," |  Number of trials: ",M,"| Range of coeffs: [0,p^",N+1,")";
print("--------------------------------------------------------------------------------------------");

for j in [-1,0,1,2,3,5,7,9,15,27] do
	print j," lines: ", RR!(counts[j+2]/M);
end for;

 

