+ M2 --no-readline --print-width 103
Macaulay 2, version 1.2
with packages: Elimination, IntegralClosure, LLLBases, PrimaryDecomposition, ReesAlgebra, SchurRings,
TangentCone

i1 : R=QQ[a..d];

i2 : I = ideal(a*c-b^2,a*d^2-2*b*c*d+c^3);

i3 : loadPackage("StatePolytope");

i4 : toString initialIdeals(I)
     
o4 = {{b^2*d^2, a*d^2, a*c}, {c^4, a*d^2, a*c}, {c^3, a*c, a^2*d^2}, {c^3, b^2*c^2, a*c, a^3*d^2},
      {c^3, b^2*c^2, b^4*c, a*c, a^4*d^2}, {c^3, b^2*c^2, b^4*c, b^6, a*c}, {c^3, b^2}, {b^2, a*d^2},
      {b*c*d, b^2, a*b*d^2, a^2*d^3}, {b*c*d, b*c^4, b^2, a*b*d^2, a^2*d^4}, {c^6, b*c*d, b*c^3, b^2,
      a*b*d^3}, {c^5, b*c*d, b*c^3, b^2}}

--We see that there are 12 initial ideals, and the highest degree of a generator in any initial ideal is 6.

i5 : m6st = statePolytope(6,I)
VERTICES
1 90 61 74 81
1 90 49 98 69
1 86 49 110 61
1 80 55 116 55
1 74 63 118 51
1 70 69 118 49
1 49 111 97 49
1 65 111 49 81
1 58 118 56 74
1 54 120 64 68
1 51 120 73 62
1 49 118 83 56

o5 = {{90, 61, 74, 81}, {90, 49, 98, 69}, {86, 49, 110, 61}, {80, 55, 116, 55}, {74, 63, 118, 51},
     {70, 69, 118, 49}, {49, 111, 97, 49}, {65, 111, 49, 81}, {58, 118, 56, 74}, {54, 120, 64, 68},
     {51, 120, 73, 62}, {49, 118, 83, 56}}

i6 : isStable(6,I)

o6 = false

i7 : isStable(50,I)

o7 = false

--We begin to sense that this ideal may be Hilbert unstable for all finite m >= 6.  Let's look closer at State_6(I) to get a better sense of what is going on.  I used polymake to analyze this polytope.  It has 12 faces.  We expect it to be 3-dimensional and to lie in the plane x + y + z + w = 306, but it is actually 2-dimensional: it also lies in the plane 2*x + y - w  = 160.  In polymake I computed the list of which vertices lie in each facet (vertices in facets, or vinf below):

i8 : vinf = {{9,10},{10,11},{8,9},{7,8},{6,11},{5,6},{4,5},{3,4},{2,3},{1,2},{0,7},{0,1}};

In a two dimensional layout, this is 
             10
          9      11
        8           6 
       7             5
        0           4
          1       3
              2

It is not a regular 12-gon
 
--When m=6, the barycenter is {153/2,153/2,153/2,153/2}.  This is clearly not in the state polytope,  since it lies outside the plane 2*x + y - w = 160. 

Now, I want to compute the proximum from the barycenter to State_6(I).  Start by finding the point in the linear span of State_6(I) that is closest to the barycenter:

Let B = bary - vertex 0 = {-27/2, 31/2, 5/2, -9/2}
Let n be the vector {-3, -1, 1, 3}, which is orthogonal to the plane of the polytope, and to the other normal vector {1, 1, 1, 1}
Then the projection of B onto n is proj_n B = (14/20) {-3, -1, 1, 3}
and its orthogonal complement is orth_n B = {-57/5, 81/5, 9/5, -33/5}
vertex 0 + orth_n B = {393/5, 386/5, 379/5, 372/5}
The distance from bary to this point is 14/sqrt{20} = 3.13....  

Next, I also found the closest point to the barycenter on each face of the polytope.  

--The function distsq computes the square of the distance between two points P and Q
i9 : distsq = (P,Q) -> (sum apply(#P, i-> ((P_i)-(Q_i))^2));

--The function projcoeff computes B dot A over A dot A, a useful coefficient for the calculations below:
i10 : projcoeff = (A,B) -> ( ( sum apply(#A, i-> (A_i)*(B_i))) / ( sum apply(#A, i-> (A_i)*(A_i))) );

--Let st be a list of points.  Let L be a pair of indices.  Let Z be a point.  The function isClptBetween determines whether the closest point to Z on the line through the points st_(L_0) and st_(L_1) is between the two points st_(L_0) and st_(L_1)

i11 : isClptBetween = (st, L , Z) -> (
A:= st_(L_0) - st_(L_1);
B:= Z - st_(L_1);
projcoeff(A,B) <= 1 and projcoeff(A,B) >= 0
);

--Let st be a list of points.  Let L be a pair of indices. The function clpt finds the closest point to Z on the line segment between st_(L_0) and st_(L_1).

i12 : clpt = (st, L , Z) -> (A:={};
B:={};
newP:={};
Ps:=apply(#L, i-> st_(L_i));
if isClptBetween(st,L,Z) == true then 
(A= st_(L_0) - st_(L_1);
B= Z - st_(L_1);
newP = apply(#A, k -> projcoeff(A,B)*(A_k)) + st_(L_1);
Ps = append(Ps, newP)
);
m:=min apply(#Ps, k-> distsq(Z,Ps_k));
clpt:={};
for i from 0 to #Ps-1 do (
if distsq(Z,Ps_i) ==m then clpt = append(clpt, Ps_i));
return flatten clpt
);

--Now we can get the point on each face closest to the barycenter, and the distances:

i13 : bary = {153/2,153/2,153/2,153/2};

i14 : closestPoints = apply(#vinf, k-> clpt(m6st,vinf_k,bary))

o14 = {{54, 120, 64, 68}, {51, 120, 73, 62}, {58, 118, 56, 74}, {65, 111, 49, 81}, {49, 350/3, 257/3,
      164/3}, {185/3, 257/3, 329/3, 49}, {508/7, 456/7, 118, 352/7}, {393/5, 853/15, 1747/15, 811/15},
      {86, 49, 110, 61}, {90, 49, 98, 69}, {83, 75, 67, 81}, {90, 61, 74, 81}}
 
i15 : apply(#closestPoints, k -> distsq(closestPoints_k, bary))

o15 : {2627, 2765, 2491, 2099, 8791/3, 6481/3, 17877/7, 37357/15, 2209, 1457, 155, 449}


--Since 14/sqrt{20} < sqrt{155}, we see that the proximum from State_6(I) to the barycenter is the point computed above, {393/5, 386/5, 379/5, 372/5}.  The barycenter is hovering above the plane of the polytope, over a point in the interior of the polytope.


--Next we analyze Chow stability.  I began by computing the Mumford index for a few 1-ps using Hassett-Hyeon-Lee's MUm function in Macaulay 2.  I very quickly discovered a 1-ps for which the ribbon is Hilbert unstable for any m and Chow semistable.  (It comes from the 1-ps with weights given by the normal vector {-2,-1,0,1} to the state polytope, only I scaled it to make the weights all positive.)

i16 : MUm = (I,w,m) -> (
S = ring I;
N = numgens S;
K = coefficientRing S;
Sw = K[gens S, Weights => w, MonomialOrder => GLex];
W = map(Sw,S,vars Sw);
I = W(I);
P = hilbertPolynomial I;
inI = ideal leadTerm I;
Sbar = Sw/inI;
F = map(Sbar, Sw, vars Sbar);
Bm = flatten entries basis(m,Sw);
Bmbar = basis(m,Sbar);
PSm = #Bm;
monomialWeight = (f) -> (expf = flatten exponents f;
sum(expf, w, times));
e = apply(0..(PSm-1), i -> (if F(Bm_i) ===F(0) then 0 else 1));
TOTALWT = sum for i from 0 to PSm-1 list
product{monomialWeight(Bm_i), e_i};
mu = sum{product{N,-TOTALWT}, product{m,P(m), sum w}}
)

i17 : MUm(I,{1,2,3,4},6)

o17 : -28

i18 : MUm(I,{1,2,3,4},7)

o18 : -34

i19 : MUm(I,{1,2,3,4},8)

o19 : -40

--From this we see that P_{X,rho}(m) = -6*m +8, which means that I is Hilbert unstable for m >= 2 and Chow semistable with respect to this 1-ps rho.

--Finally, I computed the Chow polytope as follows.  My Macaulay 2 package computes the state of a monomial ideal to be the sum of the exponent vectors of the monomials lying INSIDE the ideal (the natural choice for a Grassmannian of subspaces.)  If instead we take the sum of the exponent vectors of the monomials lying OUTSIDE the initial ideal (the natural choice for a Grassmannian of quotients) we will get something closely related.  We can use this second polytope to compute the Chow polytope because for all sufficiently large m the Hilbert scheme embeds in either Grassmannian.  Now, the vertices of this "outer state polytope" are all quadratic functions of m when m >= 6.  So, we can compute the outer state polytopes for m = 6,7,8 and use these to interpolate the vertex functions; then the leading coefficient of these functions give the Chow point.

--Here is an example: Let J be the first initial ideal given in line o4 above, that is, ideal({b^2*d^2, a*d^2, a*c}).  Then 
statePolytope(6,J) = {90,61,74,81}
outerStatePolytope(6,J)  = {binom{m+N-1}{m} * m/N,...,binom{m+N-1}{m} * m/N} - statePolytope(6,J)
So with m=6 and N=4 we get
outerStatePolytope(6,J) = {126,126,126,126} - {90,61,74,81} = {36,65,52,45}

Similarly
outerStatePolytope(7,J) = {49,90,74,60}
outerStatePolytope(8,J) = {64,119,100,77}

We interpolate: outerStatePolytope(m,J) = {m^2, 2*m^2-m-1, 2*m^2-4*m+4, m^2+2*m-3}

and the leading coefficients give the Chow point {1,2,2,1}.

--This procedure is implemented by the functions below.  

i20 : outerStatePolytope = (m,I) -> (
      N:= numgens ring(I);
      totdeg:= ((m+N-1)!*m / (m!*N!));
      st:=statePolytope(m,I);
      apply(#st, j -> apply(N, k -> totdeg - st_j_k))   );

i21 : chowPt = (inid,LL) -> ( 
      M := matrix apply(LL, k-> flatten outerStatePolytope(k,inid));
      vdM:=matrix apply(LL, k-> apply(#LL, l -> k^(#LL-l-1)/1));
      N:=vdM^(-1);
      flatten entries((N*M)^{0})
      );

i22 : maxdeg = (LL) -> ( max apply(#LL, i -> (degree LL_i)_0  )
          );

i23 : chowPolytope = (I) -> (
      inids:=initialIdeals(I);
      r:= max apply(#inids, j -> maxdeg(inids_j));
      d:= dim I;
      LL:=apply(d+1, i -> r+i);
      unique apply(#inids, k -> chowPt(ideal(inids_k),LL))
      );

i24 : chowPolytope(I)

o24 : {{1, 2, 2, 1}, {1, 3, 0, 2}, {3, 0, 0, 3}, {2, 0, 3, 1}}

--Note that this polytope contains the barycenter {3/2,3/2,3/2,3/2} of its ambient plane.

--I also wrote a function which returns the tuple of polynomials describing the outer state of a monomial ideal for all sufficiently large m.  


i25 : outerStatePolynomial = (inid,LL) -> ( R:=ring(inid);
      S:=R[m];
      M := matrix apply(LL, k-> flatten outerStatePolytope(k,inid));
      vdM:=matrix apply(LL, k-> apply(#LL, l -> k^(#LL-l-1)/1));
      N:=vdM^(-1);
      O:=N*M;
      P:=matrix {apply(#LL, l -> m^(#LL-l-1))};
      flatten entries(P*O)     
      )

i26 : inids = initialIdeals(I);

i27 : apply(#inids, i -> outerStatePolynomial(ideal(inids_i),{6,7,8}))

o27 = {{m^2, 2*m^2-m-1, 2*m^2-4*m+4, m^2+2*m-3}, {m^2, 3*m^2-6*m+5, 6*m-8, 2*m^2-3*m+3}, {m^2+m-2,
      3*m^2-6*m+5, 3*m-2, 2*m^2-m-1}, {m^2+3*m-8, 3*m^2-8*m+11, m+4, 2*m^2+m-7}, {m^2+6*m-20,
      3*m^2-12*m+27, 8, 2*m^2+3*m-15}, {3*m^2-12*m+20, 15*m-33, 8, 3*m^2-6*m+5}, {3*m^2-6*m+5, 3*m-3,
      6*m-7, 3*m^2-6*m+5}, {2*m^2-2*m+1, 3*m-3, 3*m^2-6*m+5, m^2+2*m-3}, {2*m^2-4, m+2, 3*m^2-8*m+10,
      m^2+4*m-8}, {2*m^2+2*m-12, 6, 3*m^2-12*m+26, m^2+7*m-20}, {3*m^2-7*m+9, 6, 15*m-37,
      3*m^2-11*m+22}, {3*m^2-6*m+5, m+2, 10*m-17, 3*m^2-8*m+10}}

--From this we can see which vertices of the outer state polytope collapse in the Chow polytope.  Also, we can see that the outer state polytope for any finite m >= 6 lies in a plane perpendicular to {-3,-1,1,3} and {1,1,1,1} (just like the inner state polytope does when m = 6).  Lastly, we can compute the distance from the outer barycenter to the outer state polytope as a function of m.  Is is (3m -16)/(4 sqrt{20}).  Dividing by m^2 and taking the limit as m goes to infinity gives 0.  

--Conclusion: this ribbon is Hilbert unstable for all finite m >= 2, but Chow semistable.