This notebook computes a family of Dolgachev surfaces with  a double and a triple fiber.
The main idea is to compute the ring \oplus H^0(aF+bS) where F is the canonical class and S
is one of the 6-sections.

(* we start by computing the generating series *)
GenSeries=Factor[1/(1-t^2)/(1-t^3)+Sum[Sum[(((2 a -1 )b -3 b^2)/2+1)t^a s^b,{a,3b,Infinity}] ,{b,1,Infinity}]+ Sum[ (((6b -1 )b -3 b^2)/2+1)  t^(3b) s^b s/(1-s), {b,1,Infinity}]+ Sum[ (((6b+2 -1 )b -3 b^2)/2+1)  t^(3b+1)s^b s/(1-s), {b,1,Infinity}]+ Sum[ (((6b+4 -1 )b -3 b^2)/2+1)  t^(3b+2)s^b s/(1-s), {b,1,Infinity}]
+ (1+t^2)s/(1-s)
]

We are going to have variables u0,u1,v1,v2,v3,v4,v5,v6 of weights 
(2,0),(3,0),(0,1),(3,1),(4,1),(4,1),(5,1),(5,1) respectively.

(* we see that there are three equations of degrees (2,8),(2,9) and (2,10) *)
Series[SeriesCoefficient[1/((1-t^2)(1-t^3)(1-s)(1-s t^3) (1- s t^4)^2 (1 - s t^5)^2) - GenSeries,{s,0,2}],{t,0,21}]

We will set up the quadrics.

Coeffs=Table[ToExpression[StringJoin["c",ToString[k]]],{k,1,500}];
uvs={u0,u1,v1,v2,v3,v4,v5,v6};
degs={{2,0},{3,0},{0,1},{3,1},{4,1},{4,1},{5,1},{5,1}};
as={a1,a2,a3,a4,a5,a6,a7,a8};

dega=as.degs;
count=0;
deg28=Solve[Join[{dega[[1]]==8,dega[[2]]==2},{a1≥0,a2≥0,a3≥0,a4≥0,a5≥0,a6≥0,a7≥0,a8≥0}],Integers];
deg29=Solve[Join[{dega[[1]]==9,dega[[2]]==2},{a1≥0,a2≥0,a3≥0,a4≥0,a5≥0,a6≥0,a7≥0,a8≥0}],Integers];
deg210=Solve[Join[{dega[[1]]==10,dega[[2]]==2},{a1≥0,a2≥0,a3≥0,a4≥0,a5≥0,a6≥0,a7≥0,a8≥0}],Integers];

Quadrics81=0;Quadrics82=0;Quadrics83=0;
For[i=1,i≤Length[deg28],i++,count++;Quadrics81+=Coeffs[[count]] Product[uvs[[k]]^(as[[k]]/.deg28[[i]]),{k,1,8}]];
For[i=1,i≤Length[deg28],i++,count++;Quadrics82+=Coeffs[[count]] Product[uvs[[k]]^(as[[k]]/.deg28[[i]]),{k,1,8}]];
For[i=1,i≤Length[deg28],i++,count++;Quadrics83+=Coeffs[[count]] Product[uvs[[k]]^(as[[k]]/.deg28[[i]]),{k,1,8}]];

Quadrics91=0;Quadrics92=0;Quadrics93=0;
For[i=1,i≤Length[deg29],i++,count++;Quadrics91+=Coeffs[[count]] Product[uvs[[k]]^(as[[k]]/.deg29[[i]]),{k,1,8}]];
For[i=1,i≤Length[deg29],i++,count++;Quadrics92+=Coeffs[[count]] Product[uvs[[k]]^(as[[k]]/.deg29[[i]]),{k,1,8}]];
For[i=1,i≤Length[deg29],i++,count++;Quadrics93+=Coeffs[[count]] Product[uvs[[k]]^(as[[k]]/.deg29[[i]]),{k,1,8}]];

Quadrics101=0;Quadrics102=0;Quadrics103=0;
For[i=1,i≤Length[deg210],i++,count++;Quadrics101+=Coeffs[[count]] Product[uvs[[k]]^(as[[k]]/.deg210[[i]]),{k,1,8}]];
For[i=1,i≤Length[deg210],i++,count++;Quadrics102+=Coeffs[[count]] Product[uvs[[k]]^(as[[k]]/.deg210[[i]]),{k,1,8}]];
For[i=1,i≤Length[deg210],i++,count++;Quadrics103+=Coeffs[[count]] Product[uvs[[k]]^(as[[k]]/.deg210[[i]]),{k,1,8}]];

count

Now we simplify the quadrics using various symmetries of the problem.

(* fix coefficients of Quadrics81,82,83 at v3^2, v3 v4, v4^2 *)
sls=Solve[{Coefficient[Quadrics81,v3^2]==1,Coefficient[Quadrics81,v3 v4]==0, Coefficient[Quadrics81,v4^2]==0,
  Coefficient[Quadrics82,v3^2]==0,Coefficient[Quadrics82,v3 v4]==1, Coefficient[Quadrics82,v4^2]==0,
  Coefficient[Quadrics83,v3^2]==0,Coefficient[Quadrics83,v3 v4]==0, Coefficient[Quadrics83,v4^2]==1}][[1]];
Quadrics81=Quadrics81/.sls;
Quadrics82=Quadrics82/.sls;
Quadrics83=Quadrics83/.sls;

(* Set some other coefficients using v3 -> b1 v3 + b2 v4 +  b3 v1 u0^2, etc., *)
Eqs={};
VarChange={v3->b1 v3 + b2 v4 + b3 v1 u0^2, v4->b4 v3 + b5 v4 + b6 v1 u0^2};
{Quadrics81n,Quadrics82n,Quadrics83n}={Quadrics81,Quadrics82,Quadrics83}/.VarChange;
lincomb= {a1,a2,a3}.{Quadrics81n,Quadrics82n,Quadrics83n};
sls=Solve[{Coefficient[lincomb,v3^2]==1,Coefficient[lincomb,v3 v4]==0, Coefficient[lincomb,v4^2]==0},{a1,a2,a3}][[1]];
AppendTo[Eqs,Coefficient[lincomb/.sls,v3 v1 u0^2]==1];
AppendTo[Eqs,Coefficient[lincomb/.sls,v4 v1 u0^2]==0];
lincomb= {a1,a2,a3}.{Quadrics81n,Quadrics82n,Quadrics83n};
sls=Solve[{Coefficient[lincomb,v3^2]==0,Coefficient[lincomb,v3 v4]==1, Coefficient[lincomb,v4^2]==0},{a1,a2,a3}][[1]];
AppendTo[Eqs,Coefficient[lincomb/.sls,v3 v1 u0^2]==0];
AppendTo[Eqs,Coefficient[lincomb/.sls,v4 v1 u0^2]==0];
lincomb= {a1,a2,a3}.{Quadrics81n,Quadrics82n,Quadrics83n};
sls=Solve[{Coefficient[lincomb,v3^2]==0,Coefficient[lincomb,v3 v4]==0, Coefficient[lincomb,v4^2]==1},{a1,a2,a3}][[1]];
AppendTo[Eqs,Coefficient[lincomb/.sls,v3 v1 u0^2]==0];
AppendTo[Eqs,Coefficient[lincomb/.sls,v4 v1 u0^2]==1];
(*
sls=Solve[Eqs,{b1,b2,b3,b4,b5,b6}];
Table[Length[sls[[l]]],{l,1,Length[sls]}]
*)
(* this works but takes some time, which means we will not do it every time *)
sls=Solve[{Coefficient[Quadrics81,v3 v1 u0^2]==1,Coefficient[Quadrics81,v4 v1 u0^2]==0,Coefficient[Quadrics82,v3 v1 u0^2]==0,Coefficient[Quadrics82,v4 v1 u0^2]==0,Coefficient[Quadrics83,v3 v1 u0^2]==0,Coefficient[Quadrics83,v4 v1 u0^2]==1}][[1]];
Quadrics81=Quadrics81/.sls;
Quadrics82=Quadrics82/.sls;
Quadrics83=Quadrics83/.sls;

(* we will have the (v5,v6) a dual basis to (v3,v4) in the pairing in the Gorenstein ring. Thus the v3 v5, v3 v6, v4 v5, v4 v6 coefficients of the (2,9) quadrics can be set *)
sls=Solve[{Coefficient[Quadrics91,v3 v5]==0,Coefficient[Quadrics91,v3 v6]==1, Coefficient[Quadrics91,v4 v5]==0,Coefficient[Quadrics91,v4 v6]==0,
Coefficient[Quadrics92,v3 v5]==1,Coefficient[Quadrics92,v3 v6]==0, Coefficient[Quadrics92,v4 v5]==0,Coefficient[Quadrics92,v4 v6]==-1,
Coefficient[Quadrics93,v3 v5]==0,Coefficient[Quadrics93,v3 v6]==0, Coefficient[Quadrics93,v4 v5]==1,Coefficient[Quadrics93,v4 v6]==0 
}][[1]];
Quadrics91=Quadrics91/.sls;
Quadrics92=Quadrics92/.sls;
Quadrics93=Quadrics93/.sls;

(* we also have symmetries v5 -> v5 + b1 u0 v2 + b2 u0 u1 v1 and similar for v6 which allows us to get rid of u0 v3 v2, u0 u1 v1 v2 in Quad91 and similar in Quad93 *)
VarChange={v5->v5 + b1 u0 v2 + b2 u0 u1 v1 , v6->v6 + b3 u0 v2 + b4 u0 u1 v1};
Quadrics91n=Quadrics91/.VarChange;
Quadrics93n=Quadrics93/.VarChange;
sls=Solve[{Coefficient[Quadrics91n,v3 u0 v2]==0,Coefficient[Quadrics91n,v3 u0 u1 v1]==0, Coefficient[Quadrics93n,v4 u0 v2]==0,Coefficient[Quadrics93n,v4 u0 u1 v1]==0},{b1,b2,b3,b4}];
(* it works *)
sls=Solve[{Coefficient[Quadrics91,v3 u0 v2]==0,Coefficient[Quadrics91,v3 u0 u1 v1]==0, Coefficient[Quadrics93,v4 u0 v2]==0,Coefficient[Quadrics93,v4 u0 u1 v1]==0}][[1]];
Quadrics91=Quadrics91/.sls;
Quadrics93=Quadrics93/.sls;

(* now for the last (2,10) quadrics we can make sure that we don't have any combinations of u0 Quad8i, 
and the standard coefficients at v5^2, v5 v6, v6^2*)

sls=Solve[{Coefficient[Quadrics101,u0 v3^2]==0,Coefficient[Quadrics101,u0 v3 v4]==0, Coefficient[Quadrics101,u0 v4^2]==0,
Coefficient[Quadrics102,u0 v3^2]==0,Coefficient[Quadrics102,u0 v3 v4]==0, Coefficient[Quadrics102,u0 v4^2]==0,
Coefficient[Quadrics103,u0 v3^2]==0,Coefficient[Quadrics103,u0 v3 v4]==0, Coefficient[Quadrics103,u0 v4^2]==0
}][[1]];
Quadrics101=Quadrics101/.sls;
Quadrics102=Quadrics102/.sls;
Quadrics103=Quadrics103/.sls;
sls=Solve[{Coefficient[Quadrics101,v5^2]==1,Coefficient[Quadrics101,v5 v6]==0, Coefficient[Quadrics101,v6^2]==0,
Coefficient[Quadrics102,v5^2]==0,Coefficient[Quadrics102,v5 v6]==1, Coefficient[Quadrics102,v6^2]==0,
Coefficient[Quadrics103,v5^2]==0,Coefficient[Quadrics103,v5 v6]==0, Coefficient[Quadrics103,v6^2]==1
}][[1]];
Quadrics101=Quadrics101/.sls;
Quadrics102=Quadrics102/.sls;
Quadrics103=Quadrics103/.sls;

(* Finally, there are symmetries v2-> a1 v2 + a1 v1 u1, v1 -> a3 v1. We will use them to set the coeefficients of  *)
VarChange={v2->a1 v2+ a2 v1 u1, v1-> a3 v1};
sls=Solve[{Coefficient[Quadrics91/.VarChange, u1 v2^2]==1,Coefficient[Quadrics91/.VarChange, u1^2 v1 v2]==1,Coefficient[Quadrics91/.VarChange, u1^3 v1^2]==0},{a1,a2,a3}];
(* it works *)
sls=Solve[{Coefficient[Quadrics91, u1 v2^2]==1,Coefficient[Quadrics91, u1^2 v1 v2]==1,Coefficient[Quadrics91, u1^3 v1^2]==0}][[1]];
Quadrics91=Quadrics91/.sls;

We will  use the quadrics to construct the multiplication table of the free module with respect to u0, u1, v1, v2. The basis will be denoted by 1, v3, v4, v5, v6, w.

MultTable=Table[0,{i,1,6},{j,1,6}];
MultTable[[1]]={1,v3, v4, v5, v6, w};
For[k=1,k≤6,k++,MultTable[[k]][[1]]=MultTable[[1]][[k]]];
(* we will use v3 v5 = w *)
MultTable[[2]][[4]]=w;MultTable[[4]][[2]]=w;
MultTable//MatrixForm

(* we build products of v3, v4, v5, v6 using the quadric relations *)
MultTable[[2]][[2]] = v3^2 - Quadrics81;
MultTable[[2]][[3]]= v3 v4 - Quadrics82;MultTable[[3]][[2]]=MultTable[[2]][[3]];
MultTable[[3]][[3]] = v4^2 - Quadrics83;
MultTable[[2]][[5]] = v3 v6  - Quadrics91;MultTable[[5]][[2]] =MultTable[[2]][[5]];
MultTable[[3]][[4]] = v4 v5 - Quadrics93;MultTable[[4]][[3]] =MultTable[[3]][[4]];
MultTable[[3]][[5]] =v4 v6 + Quadrics92 + w - v3 v5;MultTable[[5]][[3]]=MultTable[[3]][[5]];
MultTable[[4]][[4]] = v5^2 - Quadrics101;
MultTable[[4]][[5]] = v5 v6- Quadrics102; MultTable[[5]][[4]] =MultTable[[4]][[5]] ;
MultTable[[5]][[5]] = v6^2 - Quadrics103;

(* nNw we construct multiplication with w. The idea is that v3 w = v3 v3 v5 = MT[[2]][[2]] v5, v4 w = v3 v4 v5, v5 w = v5 v5 v3, v6 w = v3 v6 v5 *)
v3w=MultTable[[2]][[2]] v5;
MultTable[[2]][[6]] =Expand[ v3w+
Coefficient[v3w, v5 v3] (MultTable[[2]][[4]]-v5 v3) + Coefficient[v3w, v5 v4 ]  (MultTable[[3]][[4]]-v5 v4)   + Coefficient[v3w, v5^2 ] (MultTable[[4]][[4]] -v5^2)+
Coefficient[v3w, v5 v6] (MultTable[[4]][[5]]-v5 v6) ];MultTable[[6]][[2]] =MultTable[[2]][[6]] ;
v4w=MultTable[[2]][[3]] v5;
MultTable[[3]][[6]] =Expand[ v4w+
Coefficient[v4w, v5 v3] (MultTable[[2]][[4]]-v5 v3) + Coefficient[v4w, v5 v4 ]  (MultTable[[3]][[4]]-v5 v4)   + Coefficient[v4w, v5^2 ] (MultTable[[4]][[4]] -v5^2)+
Coefficient[v4w, v5 v6] (MultTable[[4]][[5]]-v5 v6) ];MultTable[[6]][[3]] =MultTable[[3]][[6]];
v6w = MultTable[[2]][[5]] v5;
MultTable[[5]][[6]] =Expand[ v6w+
Coefficient[v6w, v5 v3] (MultTable[[2]][[4]]-v5 v3) + Coefficient[v6w, v5 v4 ]  (MultTable[[3]][[4]]-v5 v4)   + Coefficient[v6w, v5^2 ] (MultTable[[4]][[4]] -v5^2)+
Coefficient[v6w, v5 v6] (MultTable[[4]][[5]]-v5 v6) ];MultTable[[6]][[5]] =MultTable[[5]][[6]];
v5w = MultTable[[4]][[4]]v3;
MultTable[[4]][[6]] =Expand[v5w +
Coefficient[v5w, v3^2](MultTable[[2]][[2]]-v3^2) + Coefficient[v5w, v3 v4](MultTable[[2]][[3]]-v3 v4) + Coefficient[v5w, v3 v5](MultTable[[2]][[4]]-v3 v5) +Coefficient[v5w, v3 v6](MultTable[[2]][[5]]-v3 v6)  
];MultTable[[6]][[4]] =MultTable[[4]][[6]];


(* now we need to compute w w = v3 v3 v5 v5 = (v3 w) v5 = MT[[2]][[6] v5*)
ww = MultTable[[2]][[6]] v5;
MultTable[[6]][[6]] = Expand[ ww+Coefficient[ww, v3 v3] (MultTable[[2]][[2]]-v3 v3) + Coefficient[ww, v3 v4] (MultTable[[2]][[3]]-v3 v4)
+ Coefficient[ww, v3 v5] (MultTable[[2]][[4]]-v3 v5) + Coefficient[ww, v3 v6] (MultTable[[2]][[5]]-v3 v6)
+ Coefficient[ww, v4 v4] (MultTable[[3]][[3]]-v4 v4)+ Coefficient[ww, v4 v5] (MultTable[[3]][[4]]-v4 v5)+ Coefficient[ww, v4 v6] (MultTable[[3]][[5]]-v4 v6)
+ Coefficient[ww, v5 v5] (MultTable[[4]][[4]]-v5 v5)+ Coefficient[ww, v5 v6] (MultTable[[4]][[5]]-v5 v6) + Coefficient[ww, v6 v6] (MultTable[[5]][[5]]-v6 v6)
+Coefficient[ww,v5 w](MultTable[[4]][[6]]-v5 w)
];

Table[ByteCount[MultTable[[k]][[l]]],{k,1,6},{l,1,6}]//MatrixForm

Let us now compute the associativity relations.

(* some basic sanity check: we want to make sure that the terms of the multiplication table are linear *)
For[i=1,i≤6,i++,For[j=1,j≤6,j++,
entry = MultTable[[i]][[j]];
If[0≠ tocheck - Coefficient[entry,v3] v3 -  Coefficient[entry,v4] v4 -  Coefficient[entry,v5] v5 -  Coefficient[entry,v6] v6 - Coefficient[entry,w] w - entry/.{v3->0,v4->0,v5->0,v6->0,w->0},Print["Waah ",{i,j}]]
]];

(* now convert the table into a more useful format *)
MultTableVector=Table[0,{i,1,6},{j,1,6}];
For[i=1,i≤6,i++,For[j=1,j≤6,j++,
entry=MultTable[[i]][[j]];
MultTableVector[[i]][[j]]={entry/.{v3->0,v4->0,v5->0,v6->0,w->0},Coefficient[entry,v3] ,  Coefficient[entry,v4],Coefficient[entry,v5] , Coefficient[entry,v6] , Coefficient[entry,w] }; 
]]

(* set up associativity relations v_i (v_j v_k ) = v_j (v_i v_k); we may assume that i<j and k > 1*)
Eqs={};
For[k=2,k≤6,k++,
For[i=1,i≤6,i++,For[j=1,j<i,j++,
tbz =Expand[Sum[(MultTableVector[[j]][[k]] )[[m]] MultTableVector[[m]][[i]],{m,1,6}]-Sum[(MultTableVector[[i]][[k]] )[[m]] MultTableVector[[m]][[j]],{m,1,6}]];
(* we want polynomials in u0, u1, v1, v2 to be zero *)
For[n=1,n≤6,n++,
Eqs=Join[Eqs,Flatten[CoefficientList[tbz[[n]],{u0,u1,v1,v2}]]];
];
]];
];
Eqs=Complement[DeleteDuplicates[Eqs],{0}];

We will now work to solve these equations. We will keep track of the quadrics as we do so. This calculation has been attempted many times and there are some idiosyncratic choices made in the process that allowed it to work. There is no claim that we get the only nine-dimensional component of the solution locus, but it may well be the case. 

The calculation runs under a minute. In the process it keeps printing out the number of equations left, their byte count and the number of variables left.

Quadrics={Quadrics81,Quadrics82,Quadrics83,Quadrics91,Quadrics92,Quadrics93,Quadrics101,Quadrics102,Quadrics103};
{Length[Eqs],ByteCount[Eqs],Length[Variables[Eqs]]}

GoodEqsList={};
For[i=1,i≤Length[Eqs],i++,
If[ByteCount[Eqs[[i]]]≤160,AppendTo[GoodEqsList,i];]
];
GoodEqsList;
sols=Solve[Table[0== Eqs[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}],{c102,c115,c116,c117,c118,c119,c124,c129,c13,c131,c136,c139,c143,c18,c20,c26,c31,c33,c39,c49,c51,c52,c62,c64,c66,c67,c81,c82,c88,c90,c95,c96,c97,c98,c99}][[1]]
Quadrics1=Quadrics/.sols;
Eqs1=Complement[DeleteDuplicates[Eqs/.sols],{0}];
{Length[Eqs1],ByteCount[Eqs1],Length[Variables[Eqs1]]}

GoodEqsList={};
For[i=1,i≤Length[Eqs1],i++,
If[ByteCount[Eqs1[[i]]]≤160,AppendTo[GoodEqsList,i];]
];
GoodEqsList;
sols1=Solve[Table[0== Eqs1[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}],{c104,c109,c111,c122,c144,c17,c6}][[1]]
Quadrics2=Quadrics1/.sols1;
Eqs2=Complement[DeleteDuplicates[Eqs1/.sols1],{0}];
{Length[Eqs2],ByteCount[Eqs2],Length[Variables[Eqs2]]}

ts=Sort[Table[{ByteCount[Eqs2[[k]]],Eqs2[[k]],k},{k,1,Length[Eqs2]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={4,15,442,78,678}
sols2=Solve[Join[{c110≠0},Table[0== Eqs2[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]]][[1]]
Quadrics3=Quadrics2/.sols2;
Eqs3=Complement[DeleteDuplicates[Eqs2/.sols2],{0}];
{Length[Eqs3],ByteCount[Eqs3],Length[Variables[Eqs3]]}

ts=Sort[Table[{Length[Eqs3[[k]]],Eqs3[[k]],k},{k,1,Length[Eqs3]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={62,38,3,8};
sols3=Solve[Join[{c17≠0,c19≠0,c123≠0},Table[0== Eqs3[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]]][[1]]
Quadrics4=Quadrics3/.sols3;
Eqs4=Complement[DeleteDuplicates[Factor[Eqs3/.sols3]],{0}];
{Length[Eqs4],ByteCount[Eqs4],Length[Variables[Eqs4]]}

ts=Sort[Table[{ByteCount[Eqs4[[k]]],Eqs4[[k]],k},{k,1,Length[Eqs4]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={40,113,193,54,87,341,58};
sols4=Solve[Join[{c17≠0,c19≠0,c123≠0},Table[0== Eqs4[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]]][[1]]
Quadrics5=Quadrics4/.sols4;
Eqs5=Complement[DeleteDuplicates[Factor[Eqs4/.sols4]],{0}];
{Length[Eqs5],ByteCount[Eqs5],Length[Variables[Eqs5]]}

ts=Sort[Table[{ByteCount[Eqs5[[k]]],Eqs5[[k]],k},{k,1,Length[Eqs5]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={2,46,503,158};
sols5=Solve[Join[{c17≠0,c19≠0,c123≠0,c110≠0,c108≠0},Table[0== Eqs5[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]],{c108,c128,c130,c5}][[1]]
Quadrics6=Quadrics5/.sols5;
Eqs6=Complement[DeleteDuplicates[Factor[Eqs5/.sols5]],{0}];
{Length[Eqs6],ByteCount[Eqs6],Length[Variables[Eqs6]]}

ts=Sort[Table[{ByteCount[Eqs6[[k]]],Eqs6[[k]],k},{k,1,Length[Eqs6]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={28,436,117};
sols6=Solve[Join[{c17≠0,c19≠0,c123≠0,c110≠0,c108≠0},Table[0== Eqs6[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]],{c19,c78,c91}][[1]]
Quadrics7=Quadrics6/.sols6;
Eqs7=Complement[DeleteDuplicates[Factor[Eqs6/.sols6]],{0}];
{Length[Eqs7],ByteCount[Eqs7],Length[Variables[Eqs7]]}

ts=Sort[Table[{ByteCount[Eqs7[[k]]],Eqs7[[k]],k},{k,1,Length[Eqs7]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={304,190,312};
sols7=Solve[Join[{c17≠0,c19≠0,c123≠0,c110≠0,c108≠0},Table[0== Eqs7[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]],{c89,c68,c75}][[1]]
Quadrics8=Quadrics7/.sols7;
Eqs8=Complement[DeleteDuplicates[Factor[Eqs7/.sols7]],{0}];
{Length[Eqs8],ByteCount[Eqs8],Length[Variables[Eqs8]]}

ts=Sort[Table[{ByteCount[Eqs8[[k]]],Eqs8[[k]],k},{k,1,Length[Eqs8]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={246,118};
sols8=Solve[Join[{c17≠0,c19≠0,c123≠0,c110≠0,c108≠0},Table[0== Eqs8[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]],{c69,c65}][[1]]
Quadrics9=Quadrics8/.sols8;
Eqs9=Complement[DeleteDuplicates[Factor[Eqs8/.sols8]],{0}];
{Length[Eqs9],ByteCount[Eqs9],Length[Variables[Eqs9]]}

ts=Sort[Table[{ByteCount[Eqs9[[k]]],Eqs9[[k]],k},{k,1,Length[Eqs9]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={4,21,1};
sols9=Solve[Join[{c17≠0,c19≠0,c123≠0,c110≠0,c108≠0},Table[0== Eqs9[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]],{c101,c36,c121}][[1]]
Quadrics10=Quadrics9/.sols9;
Eqs10=Complement[DeleteDuplicates[Factor[Eqs9/.sols9]],{0}];
{Length[Eqs10],ByteCount[Eqs10],Length[Variables[Eqs10]]}

ts=Sort[Table[{ByteCount[Eqs10[[k]]],Eqs10[[k]],k},{k,1,Length[Eqs10]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={5,66};
sols10=Solve[Join[{c17≠0,c19≠0,c123≠0,c110≠0,c108≠0},Table[0== Eqs10[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]],{c100,c60}][[1]]
Quadrics11=Quadrics10/.sols10;
Eqs11=Complement[DeleteDuplicates[Factor[Eqs10/.sols10]],{0}];
{Length[Eqs11],ByteCount[Eqs11],Length[Variables[Eqs11]]}

ts=Sort[Table[{ByteCount[Eqs11[[k]]],Eqs11[[k]],k},{k,1,Length[Eqs11]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={60,3};
sols11=Solve[Join[{c17≠0,c19≠0,c123≠0,c110≠0,c108≠0},Table[0== Eqs11[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]],{c63,c22}][[1]](*it's an Occam's razor choice *)
Quadrics12=Quadrics11/.sols11;
Eqs12=Complement[DeleteDuplicates[Factor[Eqs11/.sols11]],{0}];
{Length[Eqs12],ByteCount[Eqs12],Length[Variables[Eqs12]]}

ts=Sort[Table[{ByteCount[Eqs12[[k]]],Eqs12[[k]],k},{k,1,Length[Eqs12]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={123,195};
sols12=Solve[Join[{c17≠0,c19≠0,c123≠0,c110≠0,c108≠0},Table[0== Eqs12[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]],{c8,c9}][[1]]
Quadrics13=Quadrics12/.sols12;
Eqs13=Complement[DeleteDuplicates[Factor[Eqs12/.sols12]],{0}];
{Length[Eqs13],ByteCount[Eqs13],Length[Variables[Eqs13]]}

ts=Sort[Table[{ByteCount[Eqs13[[k]]],k},{k,1,Length[Eqs13]}]];
Table[ts[[k]],{k,1,50}];
GoodEqsList={4};
sols13=Solve[Join[{c17≠0,c19≠0,c123≠0,c110≠0,c108≠0},Table[0== Eqs13[[GoodEqsList[[k]]]],{k,1,Length[GoodEqsList]}]]][[1]]
Quadrics14=Quadrics13/.sols13;
Eqs14=Complement[DeleteDuplicates[Factor[Eqs13/.sols13]],{0}];
{Length[Eqs14],ByteCount[Eqs14],Length[Variables[Eqs14]]}

sols14=Solve[Table[0==Eqs14[[k]],{k,1,Length[Eqs14]}],{c120,c140,c141}][[1]]
Complement[DeleteDuplicates[Factor[Eqs14/.sols14]],{0}] (* we are done *)
Quadrics15=Simplify[Quadrics14/.sols14];

We can verify that a generic member of this nine-parameter family is a (2,3) Dolgachev surface. The next step is to impose a condition that the fiber at u0=u1=1 contains a line. This leads to a seven-parameter family.

(* we define a line in P^5 *)
vs={v1,v2,v3,v4,v5,v6};
LineInP5={v1->1,v2->t,v3->b1+b5 t,v4->b2+b6 t,v5->b3+b7 t,v6->b4+b8 t};
Eqs=Flatten[Table[0==SeriesCoefficient[Quadrics15[[k]]/.Join[LineInP5,{u0->1,u1->1}],{t,0,l}],{k,1,9},{l,0,2}]];
(* The choice of variables to solve for is a consequence of some trial and error approach, similar to solving for the nine-parameter family. *)
containsline=Solve[Eqs,{b3,b7,b8,c30,c32,c103,c123,c135,c137,c142}][[1]];
QuadricsLine=Quadrics15/.containsline;
Complement[Variables[QuadricsLine],uvs]

We will now consider a condition of having two lines in this fiber. This is tricky and there are likely several cases. We have found at least two different 5-parameter families. In one family the degrees of the components of the u0=u1=1 fiber are 1,1,1,3.

Line2InP5={v1->1,v2->t,v3->r1+r5 t,v4->r2+r6 t,v5->r3+r7 t,v6->r4+r8 t};
containsline2=containsline/.{b1->r1,b2->r2,b3->r3,b4->r4,b5->r5,b6->r6,b7->r7,b8->r8};
Eqs=Factor[Flatten[Table[SeriesCoefficient[Quadrics15[[k]]/.Join[Line2InP5,{u0->1,u1->1}],{t,0,l}],{k,1,9},{l,0,2}]]];
Table[containsline[[k]][[1]],{k,1,Length[containsline]}]
For[i=4,i≤10,i++,
AppendTo[Eqs,Factor[containsline[[i]][[2]]-containsline2[[i]][[2]]]];
];

Sort[Table[{ByteCount[Eqs[[k]]],k},{k,1,Length[Eqs]}]]
sls1=Solve[Eqs[[29]]==0,b4][[1]]
Timing[Eqs1=Complement[Factor[Eqs/.sls1],{0}];]
Length[Eqs1]

Sort[Table[{ByteCount[Eqs1[[k]]],k},{k,1,33}]]
sls2=Solve[Eqs1[[6]]==0,r4][[1]]
Timing[Eqs2=Complement[Factor[Eqs1/.sls2],{0}];]
Length[Eqs2]

Sort[Table[{ByteCount[Eqs2[[k]]],k},{k,1,32}]]
sls3=Solve[{Eqs2[[16]]==0,Eqs2[[6]]==0,Eqs2[[29]]==0},{r3,r7,r8}][[1]]
Timing[Eqs3=Complement[Factor[Eqs2/.sls3],{0}];]
Length[Eqs3]

Sort[Table[{ByteCount[Eqs3[[k]]],k},{k,1,29}]]
sls4=Solve[Eqs3[[1]]==0,r2][[1]]
Timing[Eqs4=Complement[Factor[Eqs3/.sls4],{0}];]
Length[Eqs4]

Sort[Table[{ByteCount[Eqs4[[k]]],k},{k,1,27}]]
sls5=Solve[{Eqs4[[22]]==0,Eqs4[[26]]==0},{c103,c135}][[1]]
Timing[Eqs5=Complement[Factor[Eqs4/.sls5],{0}];]
Length[Eqs5]

Sort[Table[{ByteCount[Eqs5[[k]]],k},{k,1,20}]]
sls6=Factor[Solve[Eqs5[[8]]==0,c137][[1]]]
Timing[Eqs6=Complement[Factor[Eqs5/.sls6],{0}];]
Length[Eqs6]

Sort[Table[{ByteCount[Eqs6[[k]]],k},{k,1,15}]]
sls7=Solve[Eqs6[[1]]==0,c30][[1]]
Timing[Eqs7=Complement[Factor[Eqs6/.sls7],{0}];]
Length[Eqs7]

Sort[Table[{ByteCount[Eqs7[[k]]],k},{k,1,14}]]
sls8=Solve[Eqs7[[2]]==0][[1]]
Timing[Eqs8=Complement[Factor[Eqs7/.sls8],{0}];]
Length[Eqs8]

Sort[Table[{ByteCount[Eqs8[[k]]],k},{k,1,6}]]
sls9=Factor[Solve[{Eqs8[[1]]==0,Eqs8[[2]]==0},{c123,c142}][[1]]]
Eqs9={};
For[k=1,k≤6,k++,
Print[k];
AppendTo[Eqs9,Factor[Eqs8[[k]]/.sls9]];
];
Eqs9=Complement[Eqs9,{0}];
Length[Eqs9]

Sort[Table[{ByteCount[Eqs9[[k]]],k},{k,1,4}]]
vars=Variables[Numerator[Eqs9[[1]]]]
Exponent[Numerator[Eqs9[[1]]][[2]],vars]
Exponent[Numerator[Eqs9[[1]]][[3]],vars]
sls10=Solve[Numerator[Eqs9[[1]]][[2]]==0,b6][[1]]
Eqs10={};
For[k=1,k≤4,k++,
Print[k];
AppendTo[Eqs10,Factor[Eqs9[[k]]/.sls10]];
];
Eqs10=Complement[Eqs10,{0}];
Length[Eqs10]

vars
Exponent[Numerator[Eqs10[[1]]],vars]
Exponent[Numerator[Eqs10[[2]]][[3]],vars]
Exponent[Numerator[Eqs10[[3]]],vars]
sls11=Solve[Numerator[Eqs10[[1]]]==0,c110][[1]];
Timing[Eqs11=Complement[Factor[Eqs10/.sls11],{0}];]

Variables[Eqs11]

sls12=Solve[Numerator[Eqs11[[1]]][[3]]==0,c32][[1]];
Timing[Eqs12=Factor[Eqs11[[2]]/.sls12];]

ByteCount[Numerator[Eqs12]]
Variables[Numerator[Eqs12]]
Exponent[Numerator[Eqs12], Variables[Numerator[Eqs12]]]
sls13=Solve[Eqs12==0,c4][[1]];

We will now try to write down the quadrics in these five variables. It is somewhat cumbersome.

(* this was developed based on some empirical computations of the connecting line. *)
btos={b1->s1-s4,r1->s1,b5->((-s1 s2+s3+s1 s4+s2 s4-(s4^2))/s4),r5->((-s1 s2+s3+s1 s4)/s4)}
Eqs12s=Numerator[Factor[Eqs12/.btos]];(* takes about a minute *)

Length[Eqs12s]
Eqs12s[[1]]
Timing[sls13s=Join[btos,Factor[Solve[0== Eqs12s[[2]],c4][[1]]]];]
ByteCount[sls13s]

Variables[Together[Quadrics15/.containsline2]]

transf={sls1,sls2,sls3,sls4,sls5,sls6,sls7,sls8,sls9,sls10,sls11,sls12,sls13s};
new={c110,c4,r1,r2,r4,r5,r6};
For[k=1,k≤13,k++,
new=Together[new/.transf[[k]]];
Print[{k, ByteCount[new]}];
];

Table[ByteCount[new[[k]]],{k,1,Length[new]}]

Quadrics3111={};
vs={v1,v2,v3,v4,v5,v6};
For[i=1,i≤6,i++, (* only the first six are easy to compute, we ran into trouble with the seventh *)
Print[i];
q=0;
For[j=1,j≤6,j++,For[l=1,l≤j,l++,
Print[{j,l}];
q+=Together[(Coefficient[Quadrics15[[i]],vs[[j]]vs[[l]]]/.containsline2)/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]] vs[[j]]vs[[l]];
]];
Print[ByteCount[q]];
AppendTo[Quadrics3111,q];
];

(* we will treat last three quadrics separately. In fact the 7-th quadric gave the most trouble, we will do it last. *)
Length[Quadrics3111]
mons210=Table[Product[uvs[[i]]^as[[i]]/.deg210[[k]],{i,1,8}],{k,1,Length[deg210]}]
Length[mons210]
For[i=8,i≤9,i++,
Print[i];
q=0;
For[j=1,j≤ Length[mons210],j++,
Print["j=",j];
q+=Together[Coefficient[Quadrics15[[i]],mons210[[j]]]/.containsline2/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]] mons210[[j]];
];
Print[ByteCount[q]];
AppendTo[Quadrics3111,q];
];

(* now we will work on the seventh quadric; start with the first 14 monomials *)
q=0;
For[j=1,j≤14,j++,
Print["j=",j];
q+=Together[Coefficient[Quadrics15[[7]],mons210[[j]]]/.containsline2/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]] mons210[[j]];
];

(* this somehow works for the 15th monomial *)
q+=Together[Factor[Factor[Coefficient[Quadrics15[[7]],mons210[[15]]]]/.containsline2]/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]]mons210[[15]];

q16=Together[Coefficient[Quadrics15[[7]],mons210[[16]]]/.containsline2/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]]mons210[[16]];

q17=Together[Coefficient[Quadrics15[[7]],mons210[[17]]]/.containsline2/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]]mons210[[17]];

q18=Together[Coefficient[Quadrics15[[7]],mons210[[18]]]/.containsline2/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]]mons210[[18]];

q19=Together[Coefficient[Quadrics15[[7]],mons210[[19]]]/.containsline2/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]]mons210[[19]];

q20=Together[Coefficient[Quadrics15[[7]],mons210[[20]]]/.containsline2/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]]mons210[[20]];

Quadric7=q+q16+q17+q18+q19+q20;
ByteCount[Quadric7]

Quadrics3111=AppendTo[Quadrics3111,Quadric7];

Table[ByteCount[Quadrics3111[[k]]],{k,1,9}]

Our next goal is to find the equations of the lines and the four intersection points of the components of the special fiber.

Line25params=Join[{u0->1,u1->1},Factor[Line2InP5/.containsline2/.Table[{c110,c4,r1,r2,r4,r5,r6}[[k]]->new[[k]],{k,1,7}]]];

transf={containsline,sls1,sls2,sls3,sls4,sls5,sls6,sls7,sls8,sls9,sls10,sls11,sls12,sls13s};
bs={b1,b2,b3,b4,b5,b6,b7,b8};
For[k=1,k≤14,k++,
bs=Together[bs/.transf[[k]]];
Print[{k, ByteCount[bs]}];
];

Variables[bs]
Line15params=Join[{u0->1,u1->1},{v1->1,v2->t,v3->bs[[1]]+bs[[5]] t,v4->bs[[2]]+bs[[6]] t,v5->bs[[3]]+bs[[7]] t,v6->bs[[4]]+bs[[8]] t}];

Sanity check: lines are in the quadrics.

Rand5=Table[{s1,s2,s3,s4,b2}[[k]]->RandomInteger[30]+1,{k,1,5}]
Timing[Quadrics3111rand=Together[Quadrics3111/.Rand5];]
Timing[Line15paramsrand=Together[Line15params/.Rand5];]
Timing[Line25paramsrand=Together[Line25params/.Rand5];]
Together[Quadrics3111rand/.Line15paramsrand]
Together[Quadrics3111rand/.Line25paramsrand]

We will now work to compute the four intersection points of the components of the special fiber. We do so by looking at vanishing of some minors of the Jacobian matrix of the fiber.

Timing[JacMat1=Together[Table[D[{Quadrics3111[[1]],Quadrics3111[[6]],Quadrics3111[[7]],Quadrics3111[[4]]}/.{u0->1,u1->1,v1->1},vs[[k]]],{k,2,5}]/.Line15params];]
Timing[DetJacMat1=Together[Det[JacMat1]];]

Timing[JacMat2=Together[Table[D[{Quadrics3111[[1]],Quadrics3111[[3]],Quadrics3111[[7]],Quadrics3111[[4]]}/.{u0->1,u1->1,v1->1},vs[[k]]],{k,2,5}]/.Line15params];]
Timing[DetJacMat2=Together[Det[JacMat2]];]

Timing[tbz1=Together[Numerator[DetJacMat1]/(t+1)];]
ByteCount[tbz1]
Timing[tbz2=Together[Numerator[DetJacMat2]/(t+1)];]
ByteCount[tbz2]

We can guess that one of the solutions has numerator and denominator which is of degree two in b2. So we could figure them out by looking at special values of b2 (and s4, which it should be of degree 4 in). The difficulty lies in normalizing the fractions.

AllNums={};
AllDens={};
Allbs={};

solst1={};
values={-3,-2,-1,1,2,4}; (* didn't work for other small values *)

For[k=1,k≤6,k++,
trythis={b2->2,s4->values[[k]]};
Print[k];
AppendTo[solst1,Together[t/.Solve[{0==tbz1/.trythis,0==tbz2/.trythis,t≠-1},t]]];
];

Table[Length[solst1[[k]]],{k,1,6}]
Table[ByteCount[solst1[[k]][[l]]],{k,1,6},{l,1,Length[solst1[[k]]]}]
nums=Table[Numerator[solst1[[k]][[1]]],{k,1,6}];
Coefficient[nums,s1^7] (* it looks like this is s4 s1^7, well normalized. *)

numsolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[nums[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]
dens=Table[Denominator[solst1[[k]][[1]]],{k,1,6}];
densolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[dens[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]

AppendTo[AllNums,numsolved];
AppendTo[AllDens,densolved];
AppendTo[Allbs,2];

solst1={};
values={-3,-2,-4,1,2,4}; (* didn't work for other small values *)

For[k=1,k≤6,k++,
trythis={b2->-2,s4->values[[k]]};
Print[k];
AppendTo[solst1,Together[t/.Solve[{0==tbz1/.trythis,0==tbz2/.trythis,t≠-1},t]]];
];

Table[Length[solst1[[k]]],{k,1,6}]
Table[ByteCount[solst1[[k]][[l]]],{k,1,6},{l,1,Length[solst1[[k]]]}]
nums=Table[Numerator[solst1[[k]][[1]]],{k,1,6}];
Coefficient[nums,s1^7] (* it looks like this is s4 s1^7, well normalized. *)

numsolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[nums[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]
dens=Table[Denominator[solst1[[k]][[1]]],{k,1,6}];
densolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[dens[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]

AppendTo[AllNums,numsolved]
AppendTo[AllDens,densolved];
AppendTo[Allbs,-2];solst1={};

solst1={};
values={-3,-4,-1,1,2,4}; (* didn't work for other small values *)

For[k=1,k≤6,k++,
trythis={b2->-3,s4->values[[k]]};
Print[k];
AppendTo[solst1,Together[t/.Solve[{0==tbz1/.trythis,0==tbz2/.trythis,t≠-1},t]]];
];

Table[Length[solst1[[k]]],{k,1,6}]
Table[ByteCount[solst1[[k]][[l]]],{k,1,6},{l,1,Length[solst1[[k]]]}]
nums=Table[Numerator[solst1[[k]][[1]]],{k,1,6}];
Coefficient[nums,s1^7] (* it looks like this is s4 s1^7, well normalized. *)

numsolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[nums[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]
dens=Table[Denominator[solst1[[k]][[1]]],{k,1,6}];
densolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[dens[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]

AppendTo[AllNums,numsolved];
AppendTo[AllDens,densolved];
AppendTo[Allbs,-3];solst1={};

(* now that we have it for three values of b we can figure it out *)
valb={2,-2,-3};
numb={ AllNums[[1]],AllNums[[2]], AllNums[[3]]}; (* normalization was required, and it was a guess that it was already good *)
Coefficient[numb,s2^4 s1^4]
Coefficient[numb,s1^7]
Coefficient[numb,s1^4s2^4]

numbsolved=Factor[(a0 + a1 b2 + a2 b2^2)/.Solve[Table[numb[[k]]== (a0 + a1 b2 + a2 b2^2)/.b2->valb[[k]],{k,1,3}],{a0,a1,a2}][[1]]]
denb={ AllDens[[1]],AllDens[[2]],AllDens[[3]]};
denbsolved=Factor[(a0 + a1 b2 + a2 b2^2)/.Solve[Table[denb[[k]]== (a0 + a1 b2 + a2 b2^2)/.b2->valb[[k]],{k,1,3}],{a0,a1,a2}][[1]]]

(* sanity check: it works! *)
Solve[0== tbz2/.{s1->13,s2->30,s3->7,s4->7,b2->3},t]
(numbsolved/denbsolved)/.{s1->13,s2->30,s3->7,s4->7,b2->3}

(* Now the moment of truth: trying to get rid of one factor *)
Timing[tbz1root2=Together[tbz1/(numbsolved - t denbsolved)];]
Timing[tbz2root2=Together[tbz2/(numbsolved - t denbsolved)];]
ByteCount[tbz1root2]
ByteCount[tbz2root2]

This sometimes gets stuck.

Timing[nicetry=Factor[tbz2root2];]

Length[nicetry]

Rand5={s1->13,s2->30,s3->5,s4->7,b2->5};

Solve[0== nicetry[[1]]/.Rand5,t]
Solve[0== nicetry[[2]]/.Rand5,t]
Solve[0== nicetry[[3]]/.Rand5,t]

Quadrand=(Quadrics3111/.Rand5)/.{u0->1,u1->1};
Solve[0== PolynomialGCD[tbz1root2/.Rand5,tbz2root2/.Rand5]]

FirstRoot=Factor[numbsolved/denbsolved]
SecondRoot=Factor[t/.Solve[0==nicetry[[3]],t][[1]]]

So now we have two roots. We then compute the corresponding points.

pt1=Factor[Line15params/.t->FirstRoot];
pt2=Factor[Line15params/.t->SecondRoot];
Table[ByteCount[pt1[[k]]],{k,1,Length[pt1]}]
Table[ByteCount[pt2[[k]]],{k,1,Length[pt1]}]

(* sanity check, should be singular pts in the fiber, so rank 3 *)
Rand5={s1->13,s2->30,s3->5,s4->7,b2->5};
Quadrand=(Quadrics3111/.Rand5)/.{u0->1,u1->1};
MatrixRank[Table[D[Quadrand[[k]],vs[[l]]]/.pt1/.Rand5,{k,1,9},{l,1,6}]]
MatrixRank[Table[D[Quadrand[[k]],vs[[l]]]/.pt2/.Rand5,{k,1,9},{l,1,6}]]

We will now repeat this process for the points on the second line. Hopefully it goes reasonably smoothly.

Timing[JacMat1=Together[Table[D[{Quadrics3111[[1]],Quadrics3111[[6]],Quadrics3111[[7]],Quadrics3111[[4]]}/.{u0->1,u1->1,v1->1},vs[[k]]],{k,2,5}]/.Line25params];]
Timing[DetJacMat1=Together[Det[JacMat1]];]

Timing[JacMat2=Together[Table[D[{Quadrics3111[[1]],Quadrics3111[[3]],Quadrics3111[[7]],Quadrics3111[[4]]}/.{u0->1,u1->1,v1->1},vs[[k]]],{k,2,5}]/.Line25params];]
Timing[DetJacMat2=Together[Det[JacMat2]];]

Timing[tbz1=Together[Numerator[DetJacMat1]/(t+1)];]
ByteCount[tbz1]
Timing[tbz2=Together[Numerator[DetJacMat2]/(t+1)];]
ByteCount[tbz2]

AllNums={};
AllDens={};
Allbs={};

solst1={};
values={-3,-2,-1,1,5,4}; (* didn't work for other small values *)

For[k=1,k≤6,k++,
trythis={b2->1,s4->values[[k]]};
Print[k];
AppendTo[solst1,Together[t/.Solve[{0==tbz1/.trythis,0==tbz2/.trythis,t≠-1},t]]];
];

Table[Length[solst1[[k]]],{k,1,6}]
Table[ByteCount[solst1[[k]][[l]]],{k,1,6},{l,1,Length[solst1[[k]]]}]
nums=Table[Numerator[solst1[[k]][[1]]],{k,1,6}];
Coefficient[nums,s1^7] (* it looks like this is s4 s1^7, well normalized. *)

numsolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[nums[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]
dens=Table[Denominator[solst1[[k]][[1]]],{k,1,6}];
densolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[dens[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]

AppendTo[AllNums,numsolved];
AppendTo[AllDens,densolved];
AppendTo[Allbs,1];

solst1={};
values={-4,-2,-1,1,5,4}; (* didn't work for other small values *)

For[k=1,k≤6,k++,
trythis={b2->-1,s4->values[[k]]};
Print[k];
AppendTo[solst1,Together[t/.Solve[{0==tbz1/.trythis,0==tbz2/.trythis,t≠-1},t]]];
];

Table[Length[solst1[[k]]],{k,1,6}]
Table[ByteCount[solst1[[k]][[l]]],{k,1,6},{l,1,Length[solst1[[k]]]}]
nums=Table[Numerator[solst1[[k]][[1]]],{k,1,6}];
Coefficient[nums,s1^7] (* it looks like this is s4 s1^7, well normalized. *)

numsolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[nums[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]
dens=Table[Denominator[solst1[[k]][[1]]],{k,1,6}];
densolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[dens[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]

AppendTo[AllNums,numsolved];
AppendTo[AllDens,densolved];
AppendTo[Allbs,-1];

solst1={};
values={-4,-2,-1,1,5,3}; (* didn't work for other small values *)

For[k=1,k≤6,k++,
trythis={b2->3,s4->values[[k]]};
Print[k];
AppendTo[solst1,Together[t/.Solve[{0==tbz1/.trythis,0==tbz2/.trythis,t≠-1},t]]];
];

Table[Length[solst1[[k]]],{k,1,6}]
Table[ByteCount[solst1[[k]][[l]]],{k,1,6},{l,1,Length[solst1[[k]]]}]
nums=Table[Numerator[solst1[[k]][[1]]],{k,1,6}];
Coefficient[nums,s1^7] (* it looks like this is s4 s1^7, well normalized. *)

numsolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[nums[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]
dens=Table[Denominator[solst1[[k]][[1]]],{k,1,6}];
densolved=Factor[(a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.Solve[Table[dens[[k]]== (a0 + a1 s4 + a2 s4^2 + a3 s4^3 + a4 s4^4 )/.s4->values[[k]],{k,1,6}],{a0,a1,a2,a3,a4}][[1]]]

AppendTo[AllNums,numsolved];
AppendTo[AllDens,densolved];
AppendTo[Allbs,3];

(* now that we have it for three values of b we can figure it out *)
valb={1,-1,3};
numb={ AllNums[[1]],AllNums[[2]], AllNums[[3]]}; (* normalization was required, and it was a guess that it was already good *)
Coefficient[numb,s2^4 s1^4]
Coefficient[numb,s1^7]
Coefficient[numb,s1^4s2^4]

numbsolved=Factor[(a0 + a1 b2 + a2 b2^2)/.Solve[Table[numb[[k]]== (a0 + a1 b2 + a2 b2^2)/.b2->valb[[k]],{k,1,3}],{a0,a1,a2}][[1]]]
denb={ AllDens[[1]],AllDens[[2]],AllDens[[3]]};
denbsolved=Factor[(a0 + a1 b2 + a2 b2^2)/.Solve[Table[denb[[k]]== (a0 + a1 b2 + a2 b2^2)/.b2->valb[[k]],{k,1,3}],{a0,a1,a2}][[1]]]

(* sanity check: it works! *)
Rand5={s1->13,s2->30,s3->7,s4->7,b2->5};
Solve[0== tbz2/.Rand5,t]
(numbsolved/denbsolved)/.Rand5

(* Now the moment of truth: trying to get rid of one factor *)
Timing[tbz1root2=Together[tbz1/(numbsolved - t denbsolved)];]
Timing[tbz2root2=Together[tbz2/(numbsolved - t denbsolved)];]
ByteCount[tbz1root2]
ByteCount[tbz2root2]

Timing[nicetry=Factor[tbz2root2];]

Length[nicetry]

Rand5={s1->13,s2->30,s3->5,s4->7,b2->5};

Solve[0== nicetry[[1]]/.Rand5,t]
Solve[0== nicetry[[2]]/.Rand5,t]
Solve[0== nicetry[[3]]/.Rand5,t]
Solve[0== nicetry[[4]]/.Rand5,t]

Quadrand=(Quadrics3111/.Rand5)/.{u0->1,u1->1};
Solve[0== PolynomialGCD[tbz1root2/.Rand5,tbz2root2/.Rand5]]

FirstRoot=Factor[numbsolved/denbsolved]
SecondRoot=Factor[t/.Solve[0==nicetry[[4]],t][[1]]]
pt3=Factor[Line25params/.t->FirstRoot];
pt4=Factor[Line25params/.t->SecondRoot];
Table[ByteCount[pt3[[k]]],{k,1,Length[pt1]}]
Table[ByteCount[pt4[[k]]],{k,1,Length[pt1]}]

We will now check the points.

(* We want to make sure that the third point is on one of the lines in the fiber; and that we have the singular points of the fiber right *)
Rand5={s1->13,s2->30,s3->5,s4->7,b2->5};
Quadrand=(Quadrics3111/.Rand5)/.{u0->1,u1->1};
MatrixRank[Table[D[Quadrand[[k]],vs[[l]]]/.pt1/.Rand5,{k,1,9},{l,1,6}]]
MatrixRank[Table[D[Quadrand[[k]],vs[[l]]]/.pt2/.Rand5,{k,1,9},{l,1,6}]]
MatrixRank[Table[D[Quadrand[[k]],vs[[l]]]/.pt3/.Rand5,{k,1,9},{l,1,6}]]
MatrixRank[Table[D[Quadrand[[k]],vs[[l]]]/.pt4/.Rand5,{k,1,9},{l,1,6}]]


(Quadrics3111/.Rand5)/.Table[uvs[[k]]-> (1/2)((uvs[[k]]/.pt1/.Rand5)+(uvs[[k]]/.pt3/.Rand5)),{k,1,8}]
(* the first and the third points are connected by a line *)

We want A_2 singularities at pt2 and pt4 and nodes at pt1, pt2, pt4.

(* clear denominators of pts *)
pt1vs=vs/.pt1;
For[i=1,i≤6,i++,
den=Denominator[pt1vs[[i]]];
pt1vs = Factor[pt1vs den];
];
pt1poly=Simplify[Join[{u0->1,u1->1},Table[vs[[k]]-> pt1vs[[k]],{k,1,6}]]];
Print["here"];

pt2vs=vs/.pt2;
For[i=1,i≤6,i++,
den=Denominator[pt2vs[[i]]];
pt2vs = Factor[pt2vs den];
];
pt2poly=Simplify[Join[{u0->1,u1->1},Table[vs[[k]]-> pt2vs[[k]],{k,1,6}]]];
Print["here"];

pt4vs=vs/.pt4;
For[i=1,i≤6,i++,
den=Denominator[pt4vs[[i]]];
pt4vs = Factor[pt4vs den];
];
pt4poly=Simplify[Join[{u0->1,u1->1},Table[vs[[k]]-> pt4vs[[k]],{k,1,6}]]];
Print["here"];

ByteCount[pt1poly]
ByteCount[pt2poly]
ByteCount[pt4poly]

This writes out preliminary computations, which can then be used to write code in C.

fname ="Desktop/MumfordFPP2021/KeumFPP/Quadricscontainsline2";
OpenWrite[fname];
Write[fname, Factor[Quadrics15/.containsline2]];
Close[fname];

fname ="Desktop/MumfordFPP2021/KeumFPP/transf";
OpenWrite[fname];
Write[fname,transf];
Close[fname];

fname ="Desktop/MumfordFPP2021/KeumFPP/Quadrics3111";
OpenWrite[fname];
Write[fname, Quadrics3111];
Close[fname];

fname ="Desktop/MumfordFPP2021/KeumFPP/ThreePointsPoly";
OpenWrite[fname];
Write[fname, {pt1poly,pt2poly,pt4poly}];
Close[fname];

fname ="Desktop/MumfordFPP2021/KeumFPP/LastPoint";
OpenWrite[fname];
Write[fname,pt3];
Close[fname];