One Hat Cyber Team

Edit File: totaro.m2

--  From: Burt Totaro <B.J.Totaro@dpmms.cam.ac.uk>
--  Date: Wed, 15 Jan 2003 19:02:53 +0000

--  Here is a more detailed description of Grobner bases for commutative
--  Z/p^r-algebras R, in particular how to find a Grobner basis.
--  This note ends with some programs I wrote, for reducing modulo
--  such a Grobner basis (with various restrictions).

--  I hope that you and Mike Stillman may want to implement
--  some of this. I wrote about Grobner bases before when R is assumed to be flat
--  over Z/p^r, but here I won't assume that.

--  Briefly, everything works as it does for Grobner bases over fields,
--  except that one treats "p" as an extra variable. For the algorithms
--  I describe, the basic example of a monomial ordering is reverse lexicographic
--  based on the order of "variables" x_1,...,x_n,p: that is, we always
--  try to reduce any expression to one divisible by a high power of p,
--  if possible.

--  Here are some more detailed definitions. Fix a monomial ordering
--  on x_1,...,x_n. Given a nonzero homogeneous polynomial f over Z/p^r
--  (or over Z localized at p), define the "leading monomial" of f
--  to be the element p^i.m, with m a monomial in x_1,...,x_n,
--  such that i is the least number with f=0 (mod p^i) and
--  m is the leading monomial of the Z/p-polynomial (f/p^i) (mod p).
--  For any homogeneous ideal I in Z/p^r[x_1,...,x_n], define the "initial ideal"
--  in(I) to be the ideal generated by the leading monomials of all
--  nonzero elements of I. (It's convenient to think of in(I) as an ideal
--  in Z[x_1,...x_n] which contains the element p^r.)
--  The initial ideal in(I) has a unique minimal
--  set of generators of the form p^i.m with m a monomial in x_1,...,x_n;
--  call these elements p^i.m the "basic non-reduced monomials" of I.
--  As usual, by the Hilbert basis theorem, there are only finitely many
--  basic non-reduced monomials, p^(i_j).m_j for 1<=j<=N.

--  We say that a Z/p^r-polynomial in x_1,...,x_n is "reduced" with respect to I
--  if the coefficient c of each monomial m in x_1,...,x_n satisfies 0 <= c < p^i,
--  where i is the least number such that p^i.m is a multiple of some
--  basic non-reduced monomial. Clearly i<=r for all monomials m.
--  It is standard to check that every element of R = Z/p^r[x_1,...,x_n]/I
--  is represented by a unique
--  reduced polynomial in x_1,...,x_n.

--  One could define a "Grobner basis" of I to be any set of elements of
--  I < Z/p^r[x_1,...,x_n] whose leading monomials generate
--  in(I) < Z/p^r[x_1,...,x_n]. Or one could make a stronger definition,
--  which makes the Grobner basis of I uniquely defined (given the
--  ordering of monomials in x_1,...,x_n): the Grobner basis of I
--  is the set of Z/p^r-polynomials g_j, 1<=j<=N (one for each basic
--  non-reduced monomial p^(i_j).m_j), such that g_j is in I,
--  g_j has leading monomial p^(i_j).m_j, and the polynomial p^(i_j)m_j - g_j
--  is reduced. Using either definition, once we have found a Grobner basis
--  for I, there is a simple algorithm to compute the reduced polynomial
--  representing any given element f(x_1,...,x_n) of R: one first makes
--  f reduced mod p, then mod p^2, and so on up to p^r. (It may save
--  time in this process to compute as much as possible over Z/p and
--  then "lift" the results; this is what I do in my program "red" to reduce
--  modulo a given Grobner basis, as discussed below.)

--  It remains to say how to compute a Grobner basis for the ideal
--  I in Z/p^r[x_1,...,x_n] generated by a given set of polynomials
--  f_1,...,f_s. Briefly, Buchberger's algorithm works, using the above
--  notion of "leading monomial". That is, one looks at the leading monomials
--  of f_1,...,f_s, for example  p.xy, p^4.yz, and so on, and then one
--  has to consider "overlaps" such as p^4.xyz which can be reduced in
--  two different ways, leading to a new element of I with a possibly
--  interesting leading monomial. As usual, this process stops
--  after finitely many steps.

--  To save time, both in finding a Grobner basis and in reducing a polynomial
--  modulo a known Grobner basis, each step of the calculation should be done
--  modulo Z/p^i for i as small as possible. For example, once we find
--  the Grobner basis over Z/p^r, it is probably a good idea to reduce the Grobner
--  basis modulo Z/p^{r-1},..., and Z/p, and to store all this information,
--  for use in the reduction algorithm.

--  ---------------

--  For an application I had in mind, I wrote programs to reduce a given
--  polynomial modulo a Grobner basis for a commutative flat Z/2^B-algebra R.
--  This is less general than what is suggested above in several ways:
--  I take p=2, I don't implement Buchberger's algorithm for finding
--  a Grobner basis, and I assume that R is flat over Z/2^B, which implies that the
--  basic non-reduced monomials are simply monomials in x_1,...,x_n
--  (rather than a power of p times such a monomial).

--  Still, I send these programs in case you want to take a look,
--  at the end of this note.
--  These programs include a particular Grobner basis over Z/2^B, where I
--  take B=7, for the ring H^*(SO(2l+1)/U(l),Z/2^B)
--  with l=5. The generators of this ring are called e_1,...,e_5,
--  in degrees 1,...,5, and e_1...e_5 is a basis for the ring in the top
--  dimension, which is 15.
--  So for example, one can compute the degree modulo 2^7 of the "spinor variety"
--  SO(11)/U(5) in its basic projective embedding (corresponding
--  to the ample divisor class e_1), by typing

--  load "filename"   (the file below)
--  red(e_1^15,B)

--  The answer is: 30e_1...e_5, which is correct. Indeed, integrally
--  one has e_1^15 = 286 e_1...e_5: that is, this spinor variety has degree 286;
--  and 286 (mod 2^7) is indeed 30.

--  Here are some aspects of these programs that might be worth noting
--  if you write similar things. First, it would be handy to have
--  Z/p^B available as a basic ring, with multiplication and division by
--  units in this ring. I just implemented these things ad hoc, by thinking
--  of polynomials over Z/2^B as integral polynomials with coefficients
--  in 0,...,2^B-1. Other useful operations one would want in dealing with
--  Z/p^B or over Z are "ord" (the p-adic valuation) and "ordpol" (minimum "ord"
--  of the coefficients of a polynomial); again, I just implemented these
--  by hand. Finally, it might be nice to have special routines for arithmetic
--  and polynomial rings over Z/p^B when p=2; computers should be able to work
--  much faster in that case.

--  It would be tempting to think of Z/p[x_1,...,x_n] as a quotient ring
--  of Z[x_1,...,x_n], but I'm not sure if that's allowed by Macaulay 2.
--  So in these programs, I used different variable names in the two rings,
--  Z[e_1,...e_5] and Z/2[E_1,...,E_5]. This is only a minor complication,
--  but it led to something interesting: I defined "redmap" (below) as the obvious
--  ring map from the Z-algebra to the Z/2-algebra, and "liftmap" as
--  the map back. Of course, there is actually no ring map back, but
--  fortunately Macaulay 2 did not give me an error message: "liftmap" gives
--  the unique lift of any Z/2-polynomial to the corresponding Z-polynomial
--  with coefficients 0 or 1. This is a very convenient operation, so
--  you might want to describe it "officially" in the documentation.
--  Likewise, if you implement polynomial rings over Z/p^B, there
--  should be a built-in way to reduce them to Z/p-polynomials, or to lift them
--  to Z-polynomials with coefficients in 0,...,p^B-1.

--  The main program below is red(f,B), which reduces a polynomial f
--  modulo the given Grobner basis (which is stored in several global variables,
--  grob0,grob1,Gmat,gmat) and modulo 2^B. After several improvements,
--  described in the comments below, this program is fairly fast. It
--  relies heavily on Macaulay 2's ability to reduce modulo a Grobner basis
--  over Z/2, and does as little arithmetic over Z/2^B as possible.

--  Finally, in my applications, I often found myself multiplying
--  two elements of R which are divisible by high powers of 2 and then
--  reducing the result. Clearly, one can save a lot of work by
--  being aware of these powers of 2 when doing the multiplication.
--  As a result, I introduced an alternate notation for elements of R,
--  namely as a pair z = (j,f), representing (2^j)f, where the polynomial f
--  only needs to be stored modulo 2^(B-j). The last two programs below
--  are red2(z,B), which reduces an expression z of this type modulo
--  the given Grobner basis, and mult(x,y,B), which multiplies two
--  expressions x and y of this type and reduces the result.
--  These programs are slightly more specialized, so I'm not
--  sure whether programs like these should be included
--  in a general-purpose set of Grobner basis routines over Z/p^B.
--  Maybe so, since they should speed things up for a user who knows
--  that he will be have to multiply elements of R which are divisible
--  by high powers of p.

--  Burt
--  To: Dan Grayson <dan@math.uiuc.edu>

-- These programs work modulo 2^B, where B is a global variables.
B=7;

-- We define a sequence  two, with  two#i = 2^i  for i=0,...,B.
two=(); i=0; while (i<=B) do (two=append(two,2^i); i=i+1);

-- recip(x,C) computes 1/x mod 2^C, where C is at most the global
-- variable B. If x is even,
-- we stop with an error.
recip = (x,C) -> (if (x%2==0) then (
    error ("Tried to compute 1/",x," mod 2^",C))
  else (C=two#C;
    y:=(1-x)%C; tot:=1; pow:=y;
-- Then 1/x = 1/(1-y) = 1+y+y^2+....
    while (pow != 0) do (tot=tot+pow; pow=(pow*y)%C);
    tot%C));

-- ord(i) computes the 2-adic order of an integer i. The program
-- gives ord(0) = infty.
ord = i -> (tot:=0; if i==0 then (tot=infty) else
      (while (i%2==0) do (tot=tot+1; i=i//2));
      tot)

-- ordpol(f) computes the 2-adic order of a polynomial over the integers.
ordpol = f -> (tot:=0; coeff:=coefficients f; sizef:= size f; i:=0;
  if (f==0) then (tot=infty) else (
    tot=ord(lift(coeff#1_(0,0),ZZ)); i=1;
    while (i<sizef) do (
      tot=min(tot,ord(lift(coeff#1_(0,i),ZZ))); i=i+1));
  tot);

R0=ZZ[e_1..e_5,Degrees=>{1..5}];

-- Here grob1 is a Grobner basis for H^*(SO(11)/U(5),Z/2^B), and grob0
-- is the corresponding list of basic non-reduced monomials.
-- This is based on the reverse lexicographic monomial ordering, as usual.
-- It happens that the usual relations (see for example Mimura and Toda's book
-- The Topology of Lie Groups) form a Grobner basis over Z/2^C for any C:
--      H^*(SO(2l+1)/U(l),Z) = 
--Z[e_1,...,e_l]/(e_i^2-2e_(i-1)e_(i+1)+2e_(i-2)e_(i+2)-...+(-1)^i e_(2i) = 0),
-- where e_i=0 for i>l.
--                         
grob0=(e_1^2,e_2^2,e_3^2,e_4^2,e_5^2);
grob1=(e_1^2-e_2, e_2^2-2*e_1*e_3+e_4, e_3^2-2*e_2*e_4+2*e_1*e_5,
e_4^2-2*e_3*e_5, e_5^2);

R1=ZZ/2[E_1..E_5,Degrees=>{1..5}];
redmap=map(R1,R0,{E_1,E_2,E_3,E_4,E_5});
liftmap=map(R0,R1,{e_1,e_2,e_3,e_4,e_5});
-- In defining "liftmap", I am using "map" in an undocumented way,
-- but it seems to work.
-- There is no ring homomorphism from the Z/2-algebra R1 to the Z-algebra R0,
-- but "liftmap" seems to provide the natural lift from Z/2-polynomials to
-- Z-polynomials with coefficients 0 or 1.

-- We define a global variable Gmat. It is a sequence such that
-- Gmat#i is a row matrix containing our Grobner basis grob1 over 2^B
-- reduced modulo 2^i, for 1<=i<=B. Here Gmat#0 should not be used,
-- so we make it an empty sequence () rather than a matrix. Having this
-- sequence Gmat saves time in our reduction procedures.
Gmat=();
  i=B; prev=grob1; len=#grob1; nxt=(); j=0;
  while i>=1 do (
    nxt=();
    j=0; while j<len do (nxt=append(nxt,(prev#j)%(two#i)); j=j+1);
    prev=nxt; Gmat=prepend(matrix {toList prev},Gmat); i=i-1);
  Gmat=prepend((),Gmat);

gmat=redmap(Gmat#1);
-- It may not be necessary to ask Macaulay 2 to compute a Grobner basis
-- for the image of the row matrix "gmat" of Z/2-polynomials, since
-- Macaulay 2 will do so automatically, whenever it first needs to. Still,
-- it seems natural to get it done now. 
-- It should be fast, since the matrix "gmat" really forms
-- a Grobner basis already, if all is well.
gb gmat;

-- red(f,C) uses a mod-2^C Grobner basis, where C is at most the
-- global variable B, to reduce f to a reduced polynomial
-- over Z/2^r. Here f must be an element of a polynomial ring over Z,
-- the same as the ring containing grob0 and grob1; in this exposition,
-- I call it R0, but any name would do.
--
-- There is a lot of choice about what order to do things in; the wrong
-- choice may involve a lot of extra work.
-- For red, I try to save time by using Macaulay 2 to do Grobner
-- basis calculations over Z/2. We begin by reducing f mod 2
-- using the Grobner basis (g_j) = (G_j mod 2). That is, we find
-- Z/2-polynomials e_j such that
--       f mod 2 = r + sum_j g_j e_j
-- where r is a reduced Z/2-polynomial. Let E_j be any lifts
-- of the polynomials e_j to Z/2^C. (It might seem preferable to make
-- a more intelligent choice of lift, but I don't see a simple way
-- to do that. The program actually takes the lifts E_j to be the unique
-- lifts to Z/2^r-polynomials whose coefficients are 0 or 1. This turns
-- out to work very well, in that the following multiplication becomes very
-- fast.) Then we replace f by
--           f - sum_j G_j E_j,
-- which is reduced mod 2.
-- 
-- To repeat the process, we need to subtract off some reduced
-- Z/2^C-polynomial to make f zero mod 2. Let us just do this in
-- a crude way: let R_0 be the lift of r to Z/2^C which has coefficients
-- 0 or 1. (We could instead subtract all monomials in f that occur in r
-- with their precise coefficient in f, but this crude way seems
-- only marginally slower.) We remember R_0,
-- and replace f by (f-R)/2 and C by C-1. When we're done, we have to
-- add up R_0 + 2R_1+...+2^(C-1)R_(C-1).
--
-- The Grobner basis must be in global variables grob0 and Gmat,
-- the latter being computed earlier from grob1.
-- Here grob0 is the sequence of basic non-reduced monomials, as elements
-- of the integral ring R0, in any order.
-- And grob1 is the corresponding sequence of Grobner relations.
-- They are polynomials over Z/2^B, represented as integer polynomials, in R0,
-- with coefficients in 0,...,2^B-1, such that the coefficient
-- of the corresponding monomial in grob0 is 1. Finally, Gmat#i,
-- for 1<=i<=B, is the sequence grob1 reduced modulo 2^i and viewed
-- as a row matrix, and gmat is the row matrix Gmat#1 viewed as a
-- row matrix over the Z/2-polynomial ring R1.
-- We also use the ring maps redmap and liftmap, defined above.
--
-- The program "red", to reduce modulo a given Grobner basis over Z/2^i,
-- is significantly faster than several earlier versions. In one test,
-- I tried a big calculation using an earlier version of "red",
-- which does this reduction "by hand", without using Macaulay 2's
-- Grobner basis routines over Z/2, but otherwise using some natural tricks
-- to save calculation, such as storing the given Grobner basis modulo
-- 2^B, ..., and 2^1, where I had B=7. 
-- This took 196 seconds. Then I did the same
-- calculation using another version of "red", which used Macaulay 2's
-- Grobner basis routines over Z/2, but without bothering to store
-- the given Grobner basis over different powers of 2 (in effect,
-- it used just Gmat#B instead of Gmat#(B-j), in "red", below).
-- This took 80 seconds. Finally, I did the same calculation using
-- the version of red below, which uses both tricks: it uses Macaulay 2's
-- Grobner basis routines over Z/2 and also stores the given Grobner basis over
-- 2^B,...,and 2^1, where I had B=7. This took 51 seconds,
-- a convincing success.
--
red= (f,C) -> (j:=0; h:=f; tot:=0_(ring f); hred:=0;
  Emat:=0; R:=0;
-- We will be looking at a polynomial h (which changes) modulo 2^(C-j).
-- At the end, the answer will be tot.
  while (j<C) do (
    h=h%(two#(C-j));
    hred=redmap(h);
    R=liftmap(hred%gmat);
    Emat=liftmap(hred//gmat);
    h=h-(((Gmat#(C-j))*Emat)_(0,0))-R;
-- At this point, h is zero modulo 2.
    h=h//2;
    tot=tot+(two#j)*R;
    j=j+1);
-- Here I will use that liftmap yields integer polynomials with coefficients
-- 0 or 1. It follows that tot = R_0 + 2R_1 + ... +2^(C-1)R_(C-1),
-- in the above notation, has coefficients in the range 0,...,2^C-1,
-- and so I don't need to reduce it mod 2^C.
  tot);

-- red2(z,C) uses a Grobner basis modulo 2^C (maybe at first defined
-- modulo a higher power of 2) to reduce z to a reduced polynomial over Z/2^C.
-- Here z is a sequence (i,f) where 0<=i<=C-1 and f is an element
-- of a polynomial ring over Z, the same as the ring containing grob0
-- and grob1; in exposition, I call it R0, but any name would do.
-- Here z means (2^i)f mod 2^C, so in particular f is only meaningful
-- mod 2^(C-i). Also, we allow the possibility z=(infty, 0_R0),
-- representing 0. The output of red2(z,C) is the same kind of pair.
--
red2 = (z,C) -> (tot:=(); tot0:=0; tot1:=0;
  if (#z != 2) then (error("First input to red2 should be a pair")) else (
    if (z#0===infty) then (tot=(infty,0_(ring z#1))) else (
      if (z#0>=C) then (tot=(infty,0_(ring z#1))) else (
	tot1=red(z#1,C-z#0);
	if (tot1 == 0) then (tot=(infty,0_(ring z#1))) else (
	  tot0=ordpol(tot1);
	  tot=(z#0+tot0,tot1//(two#tot0)))));
    tot));

-- mult(x,y,C) computes x*y (mod 2^C) in reduced form. Here x and y
-- are pairs of the form (i,f), meaning (2^i)f (mod 2^C), and
-- the output of mult is in the same form. Zero is denoted by
-- (infty, 0_R0).
--
mult = (x,y,C) -> (tot:=(); tot0:=0; tot1:=0; ex:=0;
  if ((x#0===infty) or (y#0===infty)) then (tot=(infty, 0_(ring x#1))) else (
    if (x#0+y#0 >= C) then (tot=(infty, 0_(ring x#1))) else (
      ex=C-x#0-y#0;
      tot1=((x#1)%(two#ex)) * ((y#1)%(two#ex));
      tot=red2((x#0+y#0,tot1),C)));
  tot);

-- a quick test!
assert( red(e_1^15,B) == 30*e_1*e_2*e_3*e_4*e_5 )