-- -- Free Resolutions of Ideals of Points -- -- DE+SP 96--97, last updated 02/00 -- Routines: -- points -- randomPoints -- canonicalPoints -- expectedBettiMRC -- minmaxBetti -- selfAssociated -- macaulayRep -- macaulayRepToInteger -- numLexSeg -- lexSeg -- galeTransform points = method() points Matrix := (M) -> ( V := vars ring M; if not isPolynomialRing(ring M) then error "expected a matrix over a polynomial ring"; if rank source V =!= rank target M then error "given matrix has not expected number of rows"; Points := ideal(V*(gens ker transpose M_{0})); scan(rank(source M)-1, i -> Points = intersect(Points,ideal(V*(gens ker transpose M_{i+1})))); Points ) document {points, TT "points(M)", " -- returns the ideal of the set of points in Proj(R) whose coordinates are the columns of M, a matrix over the polynomial ring R."} randomPoints = method() randomPoints (Ring, ZZ) := (R,n) -> ( if not isPolynomialRing(R) then error "expected a polynomial ring"; if (n<=0) then error "expected a positive integer"; d := numgens R; points(random(R^d,R^n)) ) randomPoints (ZZ, ZZ) := (r,n) -> ( if (n < 1) or (r < 1) then error "expected a positive integers"; R := ZZ/31991[x_0..x_r]; randomPoints(R,n) ) randomPoints (Ring,ZZ,ZZ) := (R,n,whatever) -> ( if not isPolynomialRing(R) then error "expected a polynomial ring"; if (n<=0) then error "expected a positive integer"; d := numgens R; M := id_(R^d) | transpose matrix(R,{toList(d:1)}) | random(R^d,R^(n-d-1)); points(M) ) document{randomPoints, TT "randomPoints(R,n)"," -- given a ring R, with r+1 variables, and a number n the script returns the homogeneous ideal of n random points.", BR,NOINDENT, TT "randomPoints(r,n)","-- given positive integers r and n, returns the homogeneous ideal of n random points over a polynomial ring in r+1 variables. Characteristic of the base field is assumed to be 31991", BR,NOINDENT, TT "randomPoints(r,n, whatever)", "--given a ring R, with r+1 variables, and a number n >= r+1, the script returns an r+1 x n matrix whose columns represent n random points including the standard simplex." } EXAMPLE /// R=ZZ/101[x_0..x_6] I=randomPoints(R,11) hilbertPolynomial(R/I) betti res I -- total: 1 17 46 46 30 18 4 -The first counterexample to the MRC -- 0: 1 . . . . . . -- 1: . 17 46 45 5 . . -- 2: . . . 1 25 18 4 /// EXAMPLE /// R=ZZ/101[x_0..x_8] I=randomPoints(R,13,1) hilbertPolynomial(R/I) betti res I -- total: 1 32 136 266 280 150 51 24 4 -Another counterexample to the MRC -- 0: 1 . . . . . . . . -- 1: . 32 136 266 280 140 2 . . -- 2: . . . . . 10 49 24 4 /// canonicalPoints = method() canonicalPoints Matrix := (M) -> ( R := ring M; if not isPolynomialRing(R) then error "expected a polynomial ring"; r := (rank target M)-1; if (numgens R =!= r+1) then error "given matrix has wrong size"; V := vars R; dualmat := gens ker M; s := (rank source dualmat)-1; n := rank source M; mult := matrix(R, table(n,n^2,(i,j)->(if ((j//n==j%n) and (j%n==i)) then 1 else 0))); prod := mult*((transpose M)**dualmat); V = id_(R^{1})**V; (V**id_(R^(s+1)))*(gens ker prod) ) document {canonicalPoints, TT "canonicalPoints(M)"," -- given an (r+1) x n matrix M over a ring with r+1 variables, interpreted as a set of n points in P^r, the script produces the linear part of the presentation matrix of w_{>=-1}, where w is the canonical module of the affine cone over the points. It is necessary for this to assume that no subset of n+1 of the points is linearly dependent. The presentation is actually a presentation of w if the points do not lie on a rational normal curve (so there are no quadratic relations on w_{>=-1}) and impose independent conditions on quadrics (so the homogeneous coordinate ring is 3-regular, and w is generated in degree -1)."} EXAMPLE /// R = ZZ/101[x_0..x_6] M = random(R^7,R^(11)) omega = canonicalPoints(M) betti res coker omega -- total: 4 18 30 46 46 17 1 The first counterexample to MRC -- -1: 4 18 25 1 . . . -- 0: . . 5 45 46 17 . -- 1: . . . . . . 1 /// expectedBettiMRC = (r,n) -> ( -- define the integer d = e-1 e := 1; while binomial(r+e,e)<= n do e=e+1; d := e-1; a := n-binomial(r+d,d); toprow := apply(1..r,i-> max((binomial(d+i-1,i-1))*(binomial(r+d,d+i))-a*(binomial(r,i-1)),0)); bottomrow := apply(1..r, i-> max(a*(binomial(r,i))-(binomial(d+i,i))*(binomial(r+d,d+i+1)),0)); R := ZZ/2[t]; C := new ChainComplex; C.ring = R; C#0 = R^{0}; apply(1..r, i->C#i = R^({toprow#(i-1):-d-i,bottomrow#(i-1):-d-i-1})); betti C ) document {expectedBettiMRC, TT "expectedBettiMRC(r,n)"," -- computes the expected Betti numbers of n points in P^r (cf. Minimal Resolution Conjecture)", PARA, "These are displayed as graded Betti numbers (Net) for an appropriate ", TO "ChainComplex", " over the ring ZZ/2[t]", PARA, "Here is a sample display:", EXAMPLE { "R=ZZ/101[x_0..x_3]", "I=randomPoints(R,10)", "betti res I", "expectedBettiMRC(3,10)", }, PARA, "The first counterexample to the MRC:", EXAMPLE /// R=ZZ/101[x_0..x_6] I=randomPoints(R,11) hilbertPolynomial(R/I) betti res I -- total: 1 17 46 46 30 18 4 -- 0: 1 . . . . . . -- 1: . 17 46 45 5 . . -- 2: . . . 1 25 18 4 expectedBettiMRC(6,11) --gives -- total: 1 17 46 46 30 18 4 -- 0: 1 . . . . . . -- 1: . 17 46 45 4 . . -- 2: . . . . 25 18 4 /// } selfAssociated = method() selfAssociated (Ring, Matrix) := (Q,V) -> ( if not (isAffineRing(Q) and (dim Q==0)) then error "expected a finite dimensional algebra"; B := (basis Q)**Q; n := rank source B; k := coefficientRing Q; S := k[A_0..A_(n-1)]; G := S^n; -- Make the multiplication table, over Q use ring B; m1 := 1*(B**B); -- and the row of elements of V*V vv1 := 1*(V**V); -- (the operation ** does not automatically reduce its entries; -- mult by 1 corrects this.) -- -- The columns of m2 are the coefficients of the -- elements of the row matrix m1 in terms of the basis B m2 := transpose substitute(contract(transpose m1,B),0); m2 = lift(m2,k); -- m3 is the same, expressed as a row matrix of linear forms m3 := (vars S)*m2; -- and finally the table m := adjoint(m3,G,G); -- Now do the same for VV vv2 := transpose substitute(contract(transpose vv1,B),0); vv2 = lift(vv2,k); vv := (vars S)*vv2; mbar := m**(S/ideal(vv)); (kernel mbar)==0 ) selfAssociated Ring := (Q) -> ( if not (isAffineRing(Q) and (dim Q==1) and isHomogeneous(Q)) then error "expected the homogeneous coordinate ring of a finite scheme"; u := ideal(1-(vars Q) * random(source vars Q, Q^{-1})); Qbar := Q/u; selfAssociated(Qbar,(vars Q)**Qbar) ) document {selfAssociated, TT "selfAssociated(Q,V)", "-- Given a finite scheme in a projective space represented by a FINITE DIMENSIONAL algebra Q, and a row matrix V giving generators for the subspace of linear forms, the function returns TRUE if the finite scheme is self-associated.", BR,NOINDENT, "The idea is to form the multiplication table of Q as a matrix of linear forms over a new polynomial ring S = Sym Q whose variables correspond to a basis of Q. We then tensor the matrix with S/(VV), where VV is the space of linear forms in S coming from the subspace V*V of Q. This reduction represents the generic pairing on Q that annihilates VV. Self-association means that it becomes an isomorphism on generic specialization to a symmetric matrix over the ground field, that is, if its rank is equal to the dimension of Q.", BR,NOINDENT, TT "selfAssoc(R)", "-- Given the homogeneous coordinate ring R of a finite scheme in projective space, this function returns TRUE if the scheme is self-associated.", BR,NOINDENT, "We first reduce mod a generic linear form, then call the first version of", TO "selfAssociated", "to deal with the inhomogeneous case. See the first aprt of the documentation for more information.", BR,NOINDENT, "-- CAVEATS: The linear form must not vanish on the scheme; if the field is too small this could be a problem. Could probably do it a little faster by replacing one of the variables by 1 in the matrix representing the linear forms. Also, if one of the variables does not vanish at any point of x, using it as the linear form would probably make the computation faster.", PARA, EXAMPLE, /// -- the set consisting of 4 points on each of -- two skew lines in P3 is self-associated R=ZZ/101[a,b,c,d] I= saturate (intersect(ideal(a,b), ideal(c,d))+ideal((a+c)^4)) hilbertPolynomial(R/I) x=random(R^1, R^{-1}) Q=R/(I+ideal(x-1)) V=matrix{{1,b,c,d}} selfAssociated (Q,V) ///, EXAMPLE, /// -- While 5 points on one and three points on the other --is not. R=ZZ/101[a,b,c,d] use R I=intersect( ideal(a,b,c), ideal(a,b,d), ideal(a,b,c+d), ideal(a,b,c-d), ideal(a,b,c+7*d), ideal(c,d,a), ideal(c,d,b), ideal(c,d,a+b) ) hilbertPolynomial(R/I) x=random(R^1, R^{-1}) Q=R/(I+ideal(x-1)) V=matrix{{1,b,c,d}} selfAssociated (Q,V) ///, EXAMPLE, /// -- the set consisting of 4 points on each of -- two skew lines in P3 is self-associated S=ZZ/101[a,b,c,d] I= saturate (intersect(ideal(a,b), ideal(c,d))+ ideal((a+c)^4)) selfAssociated (S/I) ///, EXAMPLE, /// --- The base locus of a net of quadrics in P^3 is --- self associated R=ZZ/101[x_0..x_3] I=ideal(random(R^1, R^{3:-2})) selfAssociated (R/I) /// } macaulayRep = method() macaulayRep (ZZ,ZZ) :=(a,d) -> ( if (a < 0 or d < 0) then error "expected positive integers"; if a == 0 then {} else if d == 1 then {a} else ( b := d; while (binomial(b,d) <= a) do b = b+1; prepend(b-1, macaulayRep(a-binomial(b-1,d),d-1))) ) document{macaulayRep, TT "macaulayRep(a,d)"," -- returns the list {m_0,..,m_k} of positive integers such that a = binomial(m_0,d)+..+binomial(m_k,d-k)", PARA, "Here is an example:", EXAMPLE { "macaulayRep(6,2)", "macaulayRep(111,3)", } } macaulayRepToInteger =method() macaulayRepToInteger (List, ZZ) := (L,d) -> ( if d<0 then error "expected positive integers"; sum(#L, i -> binomial(L_i,d-i)) ) document{macaulayRepToInteger, TT "macaulayRepToInteger({a_0,..,a_s},d)","-- returns the integer binomial(a_0,d)+...+binomial(a_s, d-s). Thus it is inverse to macaulayRep(a,d)."} TEST /// assert(macaulayRepToInteger({19, 18, 14, 11, 10, 6, 5}, 8) == 111111) /// numLexSeg = method() numLexSeg (ZZ,ZZ,ZZ) := (a,d,e) -> ( if (e i+e-d),e) ) document{numLexSeg, TT "numLexSeg(a,d,e)", " Returns the minimal codimension in degree e of an ideal of codimension a in degree d. It is attained for the lexicographic segment ideal. Here we must have e>=d>=1."} makeHilbertFcn = L -> ( b := 1; apply( #L , i -> if i==0 or i == 1 then b = L_i else b = min(L_i,numLexSeg(b,i-1,i)) ) ) document{makeHilbertFcn, TT "makeHilbertFcn {1,a1,..,as}"," -- Returns the first s terms of the termwise largest Hilbert function of a cyclic module S/I that is termwise <= {1,a1,..,as}. The a_i's should all be non-negative integers"} lexSeg = method() lexSeg (Ring, List) := (R, L) -> ( if not isPolynomialRing(R) then error "expected a polynomial ring"; C := makeHilbertFcn L; I := matrix(R,{{}}); mm := vars R; sizemm := numgens R; scan(# C , i -> ( if i == 1 then ( I = submatrix(mm,{0},{0..sizemm-C_1-1})); if i > 1 then ( d := numLexSeg(C_(i-1),i-1,i); if C_i < d then ( mm = symmetricPower(i,vars R); sizemm = rank source mm; I = I | submatrix(mm,{0},{sizemm-d..sizemm-C_i-1}); ) ) ) ); ideal I) document{lexSeg, TT "lexSeg(R, {1,a1,..,a_s})"," -- Given a list {1,a1,..,as} of nonnegative integers, and a ring R with at least a1 variables, the script returns the smallest lexicographic segment ideal I over R such that the Hilbert function of R/I begins with 1,a1,b2,..,bs and bi <= ai for each i. (Thus if the given sequence is a possible Hilbert function, bi = ai for all i, and the ideal I is generated in degrees <=s."} TEST /// assert(macaulayRep(27,5) == {7,5,3}) assert(macaulayRepToInteger(macaulayRep(250,10),10)==250) assert(makeHilbertFcn{3,4,3,2,3,2,1} == {3, 4, 3, 2, 2, 2, 1}) R = ZZ/31991[x_1..x_3]; assert(lexSeg(R,{1,3,2}) == ideal(x_1^2, x_1*x_2, x_1*x_3, x_2^2)) /// minmaxBetti = method() minmaxBetti (ZZ, ZZ) := (r,n) -> ( if (r<1 or n<1) then error "expected positive integers"; print expectedBettiMRC(r,n); e := 0; while (binomial(r+e,e) <= n) do e = e+1; -- e is the minimal degree of a generator of the ideal H := apply(e, i -> binomial(r+i-1,i)); H = append(H,n-binomial(r+e-1,e-1)); if last H =!= 0 then H = append(H,0); R := ZZ/2[x_0 .. x_(r-1)]; betti res lexSeg(R,H) ) minmaxBetti (Ring, ZZ) := (R,n) -> ( if not isPolynomialRing(R) then error "expected a polynomial ring"; r := (numgens R) - 1; minmaxBetti(r,n) ) document{minmaxBetti, TT "minmaxBetti(r,n)"," -- given integers r,n, produces the Betti diagrams of the expected (MRC) resolution of n points in P^r, and the Betti diagram of the lexicographic segment ideal with the same hilbert function as the general hyperplane section of the ideal of the general points (in one fewer variables); this is the maximum possible resolution for a purely one-dimensional ideal with the same Hilbert function as the points.", PARA, EXAMPLE /// -- The first two are the smallest counterexamples to -- the Minimal Resolution Conjecture: -- minmaxBetti(6,11) -- total: 1 17 46 45 29 18 4 -- 0: 1 . . . . . . -- 1: . 17 46 45 4 . . -- 2: . . . . 25 18 4 -- -- total: 1 22 75 115 94 40 7 -- 0: 1 . . . . . . -- 1: . 17 51 69 50 19 3 -- 2: . 5 24 46 44 21 4 -- minmaxBetti(7,12) minmaxBetti(9,20) /// } galeTransform = method() galeTransform (Matrix,Ring) := (M, S) -> ( V := vars ring M; if not isPolynomialRing(ring M) then error "expected a matrix over a polynomial ring"; if coefficientRing(ring M)=!=coefficientRing(S) then error "expected rings with the same coefficient rings"; if rank source V =!= rank target M then error "given matrix has not expected number of rows"; P:=transpose gens ker M; points(substitute(P,S))) galeTransform (Matrix) := (M) -> ( V := vars ring M; if not isPolynomialRing(ring M) then error "expected a matrix over a polynomial ring"; if rank source V =!= rank target M then error "given matrix has not expected number of rows"; transpose gens ker M) galeTransform (Ideal, Ring) := (I,S) -> ( V := vars ring I; r := numgens ring I; if not isPolynomialRing(ring I) then error "expected a matrix over a polynomial ring"; lCan := transpose((res I).dd_(r-1)); Big := R**S; M := substitute(vars S, Big)*substitute(lCan, Big); Madj := diff(transpose substitute(V,Big),M); ideal(transpose (res coker(substitute(Madj,S))).dd_((numgens S)-1))) document {galeTransform, TT "galeTransform(M,S)", " -- given the matrix of coordinates M for a set of points the script returns the homogeneous ideal of the Gale transform in Proj(S). This is well defined only up to PGL ...", BR,NOINDENT, TT "galeTransform(M)", "-- given the matrix of coordinates M for a set of points the script returns the kernel of the given matrix", BR,NOINDENT, TT "galeTransform(I,S)", "-- given the homogeneous ideal I of a set of points the script returns the homogeneous ideal of the Gale transform in Proj(S). Again this is well defined only up to PGL ...", PARA, "The Gale transform of 6 points in P^2:", EXAMPLE /// R=ZZ/101[x_0..x_2] M=transpose matrix(R,{{1,0,0},{0,1,0},{0,0,1},{1,1,1},{1,2,4},{1,7,11}}) S=ZZ/101[y_0..y_2] I=galeTransform(M,S) betti res I ///, PARA, "The Gale transform of 7 points in P^2:", EXAMPLE /// R=ZZ/101[x_0..x_2] M=transpose matrix(R,{{1,0,0},{0,1,0},{0,0,1},{1,1,1},{1,2,4},{1,7,11},{1,-1,19}}) J=points(M) betti res J S=ZZ/101[y_0..y_2] I=galeTransform(J,S) betti res I /// }