-- Example 8c: a genus 8 reducible surface with four components, each generically nonreduced
S = QQ[x_12,x_13,x_14,x_15,x_16,x_23,x_24,x_25,x_26,x_34,x_35,x_36,x_45,x_46,x_56];
LNonRedSurf8 = {x_12-x_34, x_12-x_13, x_14-x_45, x_16-x_23, x_24-x_35, x_26-x_36};
Gr26 = {x_12*x_34-x_13*x_24+x_14*x_23,
x_12*x_35-x_13*x_25+x_15*x_23,
x_12*x_36-x_13*x_26+x_16*x_23,
x_12*x_45-x_14*x_25+x_15*x_24,
x_12*x_46-x_14*x_26+x_16*x_24,
x_12*x_56-x_15*x_26+x_16*x_25,
x_13*x_45-x_14*x_35+x_15*x_34,
x_13*x_46-x_14*x_36+x_16*x_34,
x_13*x_56-x_15*x_36+x_16*x_35,
x_14*x_56-x_15*x_46+x_16*x_45,
x_23*x_45-x_24*x_35+x_25*x_34,
x_23*x_46-x_24*x_36+x_26*x_34,
x_23*x_56-x_25*x_36+x_26*x_35,
x_24*x_56-x_25*x_46+x_26*x_45,
x_34*x_56-x_35*x_46+x_36*x_45};
I = ideal join(LNonRedSurf8,Gr26);
hilbertPolynomial(I,Projective=>false)
radI = radical(I);
I==radI
pdI = primaryDecomposition(I);
radpdI = apply(pdI, J -> radical J);
apply(#pdI, i -> pdI_i==radpdI_i)
apply(pdI, J -> hilbertPolynomial(J,Projective=>false))
apply(radpdI, J -> hilbertPolynomial(J,Projective=>false))
matrix apply(#pdI, i -> apply(#pdI, j -> dim(radpdI_i+radpdI_j)))