/* This script verifies the primary decomposition in section 5 of the paper * "The Hilbert scheme of a diagonal in a product of projective spaces" by * Dustin Cartwright and Bernd Sturmfels. It also tests all symmetric classes of * monomial ideals for containment in the 3 symmetric classes of components. */ int ch = 0; option(redSB); // This is essential for the function tohilb /*** Utility functions */ proc linear(p, v) { /* Returns the portion of the polynomial which is linear in the variable, i.e. dp/dv evaluated at v = 0 */ matrix m = coeffs(p, v); if (nrows(m) < 2) { return (0); } else { return (m[2,1]); } } proc c(p, v1, v2) { return (-linear(linear(p, v1), v2)/leadcoef(p)); } proc c3(p, v1, v2, v3) { return (-linear(linear(linear(p, v1), v2), v3)/leadcoef(p)); } /* Returns the polynomial in i whose leading coefficient is l, or 0 otherwise. */ proc p(i, l) { int ndx; for (ndx = 1; ndx <= size(i); ndx++) { if (leadmonom(i[ndx]) == l) { return(i[ndx]); } } printf("Couldn't find leading coefficient %s in %s", l, i); return(0); } /*** Smaller generating set for Hilbert scheme */ ring hsmin = ch, (//a12, a13, a21, a23, a31, a32, //c12, c13, c23, ap12, ap13, ap21, ap23, ap31, ap32, bt12, bt23, bt31, cp12, cp13, cp21, cp23, cp31, cp32, f1, f2, f3), wp(/*2,2,2,2,2,2, 2,2,2,*/ 1,1,1,1,1,1, 2,2,2, 1,1,1,1,1,1, 1,1,1); /* Constructs the map corresponding to the map to the Hilbert shceme. */ proc tohilb(ideal i) { map g = hsmin, //c(p(i, x1*x2), x1, z2), c(p(i, x1*x3), x1, z3), /* a */ //c(p(i, x1*x2), x2, z1), c(p(i, x2*x3), x2, z3), //c(p(i, x1*x3), x3, z1), c(p(i, x2*x3), x3, z2), //c(p(i, x1*x2), y1, y2), c(p(i, x1*x3), y1, y3), /* c */ //c(p(i, x2*x3), y2, y3), c(p(i, x1*y2), x1, z2), c(p(i, x1*y3), x1, z3), /* a prime */ c(p(i, x2*y1), x2, z1), c(p(i, x2*y3), x2, z3), c(p(i, x3*y1), x3, z1), c(p(i, x3*y2), x3, z2), // c(p(i, x1*y2), z1, x2), c(p(i, x1*y3), z1, x3), /* a tilde */ // c(p(i, x2*y1), z2, x1), c(p(i, x2*y3), z2, x3), // c(p(i, x3*y1), z3, x1), c(p(i, x3*y2), z3, x2), c(p(i, x1*y2), z1, y2), c(p(i, x2*y3), z2, y3), /* b tilde */ c(p(i, x3*y1), z3, y1), c(p(i, x1*y2), y1, y2), c(p(i, x1*y3), y1, y3), /* c prime */ c(p(i, x2*y1), y2, y1), c(p(i, x2*y3), y2, y3), c(p(i, x3*y1), y3, y1), c(p(i, x3*y2), y3, y2), c3(p(i, y1*y2*y3), z1, y2, y3), /* f */ c3(p(i, y1*y2*y3), y1, z2, y3), c3(p(i, y1*y2*y3), y1, y2, z3); return (g); } /*** parameter space of extra components */ ring s = ch, (x1, x2, x3, y1, y2, y3, z1, z2, z3, a12, a13, a21, a23, a31, a32, // parameters for first extra comp a, b2, b3, b1, // parameter for second extra comp d1, e1, f1, d2, e2, f2, d3, e3, f3), (dp(9), dp(15)); /* Change of coordinates */ map f = s, x1 + d1*y1 + e1*z1, x2 + d2*y2 + e2*z2, x3 + d3*y3 + e3*z3, y1 + f1*z1, y2 + f2*z2, y3 + f3*z3, z1, z2, z3, a12, a13, a21, a23, a31, a32, a, b2, b3, b1, d1, e1, f1, d2, e2, f2, d3, e3, f3; ideal j = x1*x2, x2*x3, x1*x3, x1*y2 - a12*x1*z2 - a21*y1*y2 - (a31*a32-a12*a31-a21*a32)*y1*z2, x1*y3 - a13*x1*z3 - a31*y1*y3 - (a21*a23-a13*a21-a31*a23)*y1*z3, x2*y1, x2*y3, x3*y1, x3*y2, y1*y2*y3 - a23*y1*y2*z3 - a32*y1*z2*y3 - (a12*a13-a23*a12-a32*a13)*y1*z2*z3, d1; // This portion of the change of variables is unnecessary ideal i = std(f(j)); /* Second set of extra components */ ideal jj = x1*x2, x1*x3, x2*x3, x1*y2, x1*y2, x1*y3, x2*y1, x3*y1, x3*y2, x2*y3 - a*x2*z3, y1*y2*y3 - b3*y1*y2*z3 - b2*y1*z2*y3 - b1*z1*y2*(y3-a*z3); ideal ii = std(f(jj)); ideal zero = 0; /* Permutation of the 3 sets of variables. */ map cycle = s, x2, x3, x1, y2, y3, y1, z2, z3, z1, a12, a13, a21, a23, a31, a32, a, b2, b3, b1, d1, e1, f1, d2, e2, f2, d3, e3, f3; map g1 = tohilb(i); map g2 = tohilb(cycle(i)); map g3 = tohilb(preimage(s, cycle, i)); map gg1 = tohilb(ii); map gg2 = tohilb(cycle(ii)); map gg3 = tohilb(preimage(s, cycle, ii)); /* Implicitize the parametrizations of the extra components. */ setring hsmin; ideal extra1 = std(preimage(s, g1, zero)); ideal extra2 = std(preimage(s, g2, zero)); ideal extra3 = std(preimage(s, g3, zero)); ideal extra11 = std(preimage(s, gg1, zero)); ideal extra22 = std(preimage(s, gg2, zero)); ideal extra33 = std(preimage(s, gg3, zero)); /*** Parameter space of main component */ ring mainparam = ch, (x1, x2, x3, y1, y2, y3, z1, z2, z3, m11, m12, m13, m21, m22, m23, m31, m32, m33, n11, n12, n13, n21, n22, n23, n31, n32, n33, in31, im31, in11, im11, im1, in1, i1, icx2y3, icy2x3, icy1y2y3), (dp(9), dp(28)); matrix cx2y3[3][3] = m11, m21, m31, n11, n21, n31, n12, n22, n32; matrix cy2x3[3][3] = m11, m21, m31, m12, m22, m32, n11, n21, n31; ideal i = n31*in31 - 1, /* y1x3 */ m31*im31 - 1, /* y1x2 */ n11*in11 - 1, /* x1x3 */ m11*im11 - 1, /* x1x2 */ (m22*m31-m21*m32)*im1 -1, /* x1y2 */ (n22*n31-n21*n32)*in1 - 1, /* x1y3 */ (n11*m21 - n21*m11)*i1 - 1, /* x2x3 */ det(cx2y3)*icx2y3 - 1, /* x2y3 */ det(cy2x3)*icy2x3 - 1, /* y2x3 */ ((m21*m32-m31*m22)*(n12*n31-n32*n11) - /* y1y2y3 */ (m12*m31-m32*m11)*(n21*n32-n31*n22)) * icy1y2y3 - 1; i = std(i); /* This is the transpose of the usual matrix. */ matrix a[3][3] = x1, y1, z1, m11*x2 + m12*y2 + m13*z2, m21*x2 + m22*y2 + m23*z2, m31*x2 + m32*y2 + m33*z2, n11*x3 + n12*y3 + n13*z3, n21*x3 + n22*y3 + n23*z3, n31*x3 + n32*y3 + n33*z3; ideal family = std(i + minor(a, 2)); map g = tohilb(family); /* Implicitize main component */ setring hsmin; print("Beginning implicitization of main component"); ideal main = preimage(mainparam, g, i); print("End implicitization"); /*** H_{3,3} Hilbert scheme */ /* Note that at* are zero in the Hilbert scheme and can be moved to the * parameters which are eliminated when passing to hsmin. */ ring hs = ch, (x1, x2, x3, y1, y2, y3, z1, z2, z3, b12, b13, b21, b23, b31, b32, d12, d13, d23, at12, at13, at21, at23, at31, at32, bp12, bp13, bp21, bp23, bp31, bp32, bt32, bt13, bt21, dp12, dp13, dp21, dp23, dp31, dp32, e1, e2, e3, g1, g2, g3, h, a12, a13, a21, a23, a31, a32, c12, c13, c23, ap12, ap13, ap21, ap23, ap31, ap32, bt12, bt23, bt31, cp12, cp13, cp21, cp23, cp31, cp32, f1, f2, f3), (dp(9), dp(37), wp(2,2,2,2,2,2, 2,2,2, 1,1,1,1,1,1, 2,2,2, 1,1,1,1,1,1, 1,1,1)); ideal family = x1*x2 - a12*x1*z2 - a21*z1*x2 - b12*y1*z2 - b21*z1*y2 - c12*y1*y2 - d12*z1*z2, x1*x3 - a13*x1*z3 - a31*z1*x3 - b13*y1*z3 - b31*z1*y3 - c13*y1*y3 - d13*z1*z3, x2*x3 - a23*x2*z3 - a32*z2*x3 - b23*y2*z3 - b32*z2*y3 - c23*y2*y3 - d23*z2*z3, x1*y2 - ap12*x1*z2 - at12*z1*x2 - bp12*y1*z2 - bt12*z1*y2 - cp12*y1*y2 - dp12*z1*z2, x1*y3 - ap13*x1*z3 - at13*z1*x3 - bp13*y1*z3 - bt13*z1*y3 - cp13*y1*y3 - dp13*z1*z3, x2*y1 - ap21*x2*z1 - at21*z2*x1 - bp21*y2*z1 - bt21*z2*y1 - cp21*y2*y1 - dp21*z2*z1, x2*y3 - ap23*x2*z3 - at23*z2*x3 - bp23*y2*z3 - bt23*z2*y3 - cp23*y2*y3 - dp23*z2*z3, x3*y1 - ap31*x3*z1 - at31*z3*x1 - bp31*y3*z1 - bt31*z3*y1 - cp31*y3*y1 - dp31*z3*z1, x3*y2 - ap32*x3*z2 - at32*z3*x2 - bp32*y3*z2 - bt32*z3*y2 - cp32*y3*y2 - dp32*z3*z2, y1*y2*y3 - e1*x1*z2*z3 - e2*z1*x2*z3 - e3*z1*z2*x3 - f1*z1*y2*y3 - f2*y1*z2*y3 - f3*y1*y2*z3 - g1*y1*z2*z3 - g2*z1*y2*z3 - g3*z1*z2*y3 - h*z1*z2*z3; /* Generates relations among the parameters by reading off coefficients of a * normal form. This generates a warning since i is not a standard basis. */ proc rel(expr, i) { matrix m = coef(reduce(expr, i), x1*x2*x3*y1*y2*y3*z1*z2*z3); return(ideal(m[2, 1..ncols(m)])); } print("Ignore the following warnings about not being a standard basis"); /* The ideal h contains the relations of the Hilbert scheme H_{3,3}. The calls * to rel correspond to a generating set of the syzygies of the Borel-fixed * ideal Z. */ ideal h = rel(x3*family[1] - x2*family[2], family) + /* b/t xx and xx */ rel(x3*family[1] - x1*family[2], family) + rel(y3*family[1] - x2*family[5], family) + /* b/t xx and xy */ rel(y2*family[2] - x3*family[4], family) + rel(y1*family[3] - x3*family[6], family) + rel(y2*family[1] - x2*family[4], family) + rel(y3*family[2] - x3*family[5], family) + rel(y1*family[2] - x1*family[8], family) + rel(y2*family[3] - x2*family[9], family) + rel(y3*family[3] - x3*family[7], family) + rel(y1*family[1] - x1*family[6], family) + rel(x1*family[7] - x2*family[5], family) + /* b/t xy and xy */ rel(x3*family[4] - x1*family[9], family) + rel(x2*family[8] - x3*family[6], family) + rel(y2*family[5] - y3*family[4], family) + rel(y1*family[9] - y2*family[8], family) + rel(y3*family[6] - y1*family[7], family) + rel(y2*y3*family[6] - x2*family[10], family) + /* b/t xy & yyy */ rel(y2*y3*family[8] - x3*family[10], family) + rel(y1*y3*family[9] - x3*family[10], family) + rel(y1*y3*family[4] - x1*family[10], family) + rel(y1*y2*family[5] - x1*family[10], family) + rel(y1*y2*family[7] - x2*family[10], family); print("End ignore"); h = std(h); map f = hsmin, //a12, a13, a21, a23, a31, a32, //c12, c13, c23, ap12, ap13, ap21, ap23, ap31, ap32, bt12, bt23, bt31, cp12, cp13, cp21, cp23, cp31, cp32, f1, f2, f3; setring hsmin; ideal h = preimage(hs, f, h); h = std(h); if (0) { /* Verify that g parametrizes a set inside the Hilbert scheme */ poly hh; int ndx; int szh = size(h); for (ndx = 1; ndx <= szh; ndx++) { setring hsmin; hh = h[ndx]; setring mainparam; if (reduce(g(hh), i) == 0) { printf("Checked relation %s/%s", ndx, szh); } else { printf("FAILED: %s", ndx); } } } /* Now check that the Hilbert scheme is the intersection of the 7 ideals * constructed above */ setring hsmin; ideal inter = std(intersect(main, extra1, extra2, extra3, extra11, extra22, extra33)); if (std(reduce(h, inter)) == 0) { if (std(reduce(inter, h)) == 0) { print("Hilbert scheme is the union of the 7 components"); } else { print("Hilbert scheme is NOT the union of the 7 components"); } } else { print("Components are not even contained in Hilbert scheme"); } /** Now determine which monomial ideals are in which component */ ring mon = ch, (x1, x2, x3, y1, y2, y3, z1, z2, z3), dp; list l = intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 1 ideal(x2, y2, x3, z3), ideal(x1, x2, y3, z3), ideal(z3, x1, x2, y2), ideal(x2, y3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, z1, x3, y3), // 2 ideal(x2, y2, x3, z3), ideal(x1, x2, x3, y3), ideal(x3, x1, x2, y2), ideal(x2, y3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 3 ideal(x2, y2, x3, z3), ideal(x1, y2, x3, y3), ideal(x3, x1, x2, y2), ideal(x2, y3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 4 ideal(x2, y2, x3, y3), ideal(x1, y2, x3, y3), ideal(x3, y1, x2, y2), ideal(x2, y3, x1, y1)), intersect(ideal(x1, y1, x2, z2), ideal(x1, z1, x3, y3), // 5 ideal(x2, y2, x3, z3), ideal(x1, x2, x3, z3), ideal(x3, x1, x2, z2), ideal(x2, x3, x1, z1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, z1, x3, y3), // 6 ideal(x2, y2, x3, z3), ideal(x1, y2, x3, y3), ideal(x3, x1, x2, y2), ideal(x2, x3, x1, z1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 7 ideal(x2, y2, x3, z3), ideal(x1, x2, x3, z3), ideal(x3, x1, x2, z2), ideal(x2, x3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 8 ideal(x2, y2, x3, z3), ideal(x1, x2, x3, z3), ideal(z3, x1, x2, y2), ideal(x2, x3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 9 ideal(x2, y2, x3, z3), ideal(x1, x2, x3, z3), ideal(x3, y1, x2, y2), ideal(x2, x3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 10 ideal(x2, y2, x3, y3), ideal(x1, x2, x3, y3), ideal(x3, y1, x2, y2), ideal(x2, y3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, z1, x3, y3), // 11 ideal(x2, y2, x3, z3), ideal(x1, x2, x3, z3), ideal(x3, x1, x2, y2), ideal(x2, x3, x1, z1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, z1, x3, y3), // 12 ideal(x2, y2, x3, z3), ideal(x1, x2, x3, y3), ideal(x3, x1, x2, y2), ideal(x2, x3, x1, z1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 13 ideal(x2, y2, x3, y3), ideal(x1, x2, x3, z3), ideal(x3, x1, x2, y2), ideal(x2, x3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 14 ideal(x2, y2, x3, y3), ideal(x1, x2, x3, y3), ideal(x3, x1, x2, y2), ideal(x2, y3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 15 ideal(x2, y2, x3, z3), ideal(x1, x2, x3, z3), ideal(x3, x1, x2, y2), ideal(x2, x3, x1, y1)), intersect(ideal(x1, y1, x2, y2), ideal(x1, y1, x3, y3), // 16 ideal(x2, y2, x3, y3), ideal(x1, x2, x3, y3), ideal(x3, x1, x2, y2), ideal(x2, x3, x1, y1)); map cov = mon, x1, x2, x3, y1+x1, y2+x2, y3+x3, z1+2*y1+x1, z2+2*y2+x2, z3+2*y3+x3; map f; ideal i; ideal zero = 0; setring hsmin; ideal monpt; for (int ndx = 1; ndx <= 16; ndx++) { setring mon; i = l[ndx]; f = tohilb(std(cov(i))); setring hsmin; monpt = std(preimage(mon, f, zero)); printf("Monomial ideal %s:", ndx); if (std(reduce(main, monpt)) == 0) { printf(" is in the main component"); } if (std(reduce(extra1, monpt)) == 0 || std(reduce(extra2, monpt)) == 0 || std(reduce(extra3, monpt)) == 0) { printf(" is in a 14-dim. extra component"); } if (std(reduce(extra11, monpt)) == 0 || std(reduce(extra22, monpt)) == 0 || std(reduce(extra33, monpt)) == 0) { printf(" is in a 13-dim. extra component"); } }