-- Elizabeth Gasparim and Irena Swanson
-- (modified by Thomas Koeppe)
-- April 22, 2007. University of Edinburgh.
--
-- Program to compute the size of the space of holomorphic
-- sections of rank-2 vector bundles E with splitting type
-- (j, -j) and extension class p, i.e. with transition
--                         / z^j   p  \
-- functions given by M = |            | .
--                         \  0  z^-j /
--
-- On O(-1) --> CP^1, fix charts U = {(z,u)}
--                           and V = {(1/z, z u)} .
-- Then p is a polynomial in z, z^-1 and u, and sections of E
-- are functions (a,b) : C -> C^2 holomorphic in z and u
-- such that M(a,b) is holomorphic in 1/z, z u .
--
--
-- The main functions are:
--	iWidth(p, j)
--	iHeight(p, j)
-- where p is a polynomial in x and y with rational/integer coefficients
-- and j is a nonnegative integer specifying the above splitting type.
--
-- The two functions compute the instanton widths and heights
-- of the vector bundle E on the total space of O(-1) whose
-- transition matrix in canonical coordinates is of the form
--               /  z^j     p'  \
--              |                | ,
--               \   0     z^-j /
-- where p' is the sum of those terms u^i z^l in p(u,zu) for which
-- 1 <= i <= 2j-2 and 0 <= l <= j-1.


-------------------------
-- Program begins here --
-------------------------


--
-- Set up the basic environment
--

basering = QQ
S        = basering[x, y]
Rbase    = basering[z, u, MonomialOrder => GroupLex=> 1, Degrees => { {1,0}, {0,1} } ]

--
-- Here are some test polynomials
--
uuu = u^4-u^6*z^5
vvv = u^4-u^5*z^3-u^7*z^5-u^8*z^8
p1  = z^-2 * u^3 + 4 * z^-1 * u^4 + z * u^2 + 4 * z^2 * u
p2  = z^-2 * u
p3  = z^3  * u

--
-- For the "canonical form" of the polynomial p, certain coefficients
-- suffice. The function "polynomialParameters" returns (m, N),
-- where m is the smallest degree of u in the terms in p
-- and N is the maximal degree of u in the terms in p, plus 2j-2.
--
polynomialParameters = (p, j) ->
(
  m := 0;
  N := j;
  if p != 0 then
  (
    pp := terms(p);
    m   = min apply(pp, i -> (exponents i)_0_1);
    N   = max apply(pp, i -> (exponents i)_0_1) + 2*j-2;
  );

  return {m, N};
)


--
-- For a given extension p, "doRing" computes a list
-- of monomials in z and u that generate sections of E.
--
-- Monomials of excessive degree in u (>2j-2) are marked to be thrown out right away
-- by being stored in "irrelvars", but to get a basis of the space of sections we also
-- need to find relations between the generators, which is done by the next function, "getRelations".
--
-- Input: the extension class p and the splitting type j.
--
-- Return value: { list of coefficients, polynomial in u and z, a generic section (a, b), ideal of discardible coefficients }
--
doRing = (p, j) ->
(
  -- We only need to consider monomials up to certain degrees:
  m := (polynomialParameters(p, j))_0;
  N := (polynomialParameters(p, j))_1;

  -- A section is of the form
  --          __      / u^r z^s \        /    0    \
  -- (a, b) = \  a   |           | + b  |           | .
  --          /_  rs  \    0    /     rs \ u^r z^s /
  --
  aIndex1 := flatten apply(m, i -> apply(i - j + 1, l -> (i, l)));
  aIndex2 := toList flatten apply(m..N, i -> apply(i + j + 1, l -> (i, l)));
  aIndex  := aIndex1 | aIndex2;
  aVars   := apply(aIndex, i -> a_i);

  bIndex  := flatten apply(N+1, i -> toList apply(0..i+j, l -> (i, l)));
  bVars   := apply(bIndex, i -> b_i);

  -- Construct the actual section from the generating elements.
  R := basering[z, u, reverse(bVars | aVars), MonomialOrder => GroupLex=> 1];
  aPoly := sum apply(aIndex, i -> a_i * u^(i_0) *z^(i_1) );
  bPoly := sum apply(bIndex, i -> b_i * u^(i_0) * z^(i_1) );

  -- We extract the variables bVars|aVars from the new ring R again,
  -- but this time they are elements of the ring R, not just symbols.
  allVars := apply(numgens R - 2, i -> R_(i+2));

  -- All coefficients with degrees of u greater than 2j-2 are zero,
  -- so we remember those in order to discard them later.
  -- (They are remembered as ideals in R, called "irrelevantVars".)
  aIndexSel := select(aIndex, i -> (i_0 > 2*j-2));
  bIndexSel := select(bIndex, i -> (i_0 > 2*j-2));
  irrelevantVars := apply(aIndexSel, i -> a_i) | apply(bIndexSel, i -> b_i);

  F := map(R, Rbase, {z, u} );


  -- The return value is a list of
  -- { all coefficients,
  -- { the generic section (a, b),
  -- { the ideal of discardible coefficients }.
  return {allVars, F(p), aPoly, bPoly, ideal (irrelevantVars | {0_R})};
)


--
-- The function "getRelations" computes the relations between the monomial
-- generators of sections. The input is the first coordinate of the section
-- on the V-chart, i.e. M(a,b)[1] = fTv; the output is a list of expressions
-- in the coefficients expressing their relations, so they have to be set to
-- zero. (E.g. "a_10 + b_02" means "a_12 + b_02 = 0", i.e. b_02 = -a_12.)
--
-- Return value: { ideal of non-free coefficients, ideal of relations between those coefficients }
--
getRelations = (N, fTv) ->
(
  R := ring fTv;
  nonfree := ideal(0_R);
  rels := ideal(0_R);

  for k when k <= N do
  (

    -- truncate fTv to u-degrees <= k:
    tempoly := sum select(terms fTv, i -> (apply(exponents i, j -> j_{1}))_0_0 <= k);

    while tempoly != 0 do
    (

      -- The aim is for "partp" to hold the coefficient of the leading term in u and z of "tempoly" in terms of a..., b... :

      --
      -- First get the maximum total degree:
      maxSum  := max apply(terms tempoly, i -> sum((flatten exponents i)_{0,1}));
      --
      -- Then select the terms of that total degree:
      maxTerm := select(terms tempoly, i -> sum((flatten exponents i)_{0,1}) == maxSum);
      --
      -- And finally pick only those of maximal degree in u:
      maxU    := max apply(maxTerm, j -> (flatten exponents j)_1);
      maxTerm  = sum select(maxTerm, i -> (flatten exponents i)_1 == maxU);

      -- Now extract the coefficient in a..., b... from the leading term:
      partp   := maxTerm // (u^maxU * z^(maxSum-maxU));

      tempoly = tempoly - maxTerm;           -- Remove the leading term from tempoly ...
      fTv     = fTv     - maxTerm;           -- ... and also from the input polynomial.

      -- If deg(z) > deg(u), store the coefficient of z^deg(z) u^deg(u) as a relation:
      if maxSum > (2 * maxU) then
      (
        depTerm := leadTerm(partp);          -- "leadTerm" here in terms of a_rs, b_rs
        nonfree  = nonfree + ideal(depTerm);
        partp    = partp // leadCoefficient(depTerm);  -- and this is the actual relation!
        rels     = rels + ideal(partp);
      );

    );

  );

  return {nonfree, rels};
)


--
-- "RtoS" converts polynomials in u and z into polynomials in x=u and y=zu.
--
RtoS = (p) ->
(
  g := 0_S;

  for i in (terms p) do
  (
    degz := (exponents i)_0_0;
    degu := (exponents i)_0_1;
    if (degz <= degu) then g = g + y^degz * x^(degu - degz) * leadCoefficient(i);
  );

  return g;
)


--
-- "setVectors" creates the S-module from the generators of sections.
--
setVectors = (aPoly, bPoly, changeables, allVars) ->
(
  Smodule := image matrix{{0_S}, {0_S}};
  tozero  := apply(allVars, i -> i => 0);

  for k when k < #changeables do
  (
    taPoly := substitute(aPoly, changeables_k => 1);
    tbPoly := substitute(bPoly, changeables_k => 1);
    taPoly  = substitute(taPoly, tozero);
    tbPoly  = substitute(tbPoly, tozero);

    vecab  := matrix{{RtoS(taPoly)}, {RtoS(tbPoly)}};
    Smodule = Smodule + image vecab;
  );

  return (presentation image mingens Smodule);
)


--
-- Given a sequence 0 -> A -> A** -> Q -> 0,
-- Q is a skyscraper sheaf, and "qLength" computes
-- the dimension of its stalk:
--
qLength = (A) ->
(
  B := kernel transpose A;
  Q := (kernel transpose presentation B) / (image transpose gens B);
  return (degree Q);
)


--
-- Given a bundle E as above with extension class p and splitting type (j, -j),
-- there is an S-module A of holomorphic sections of E, and the qLength of that
-- is calculated by the following function "iWidth":
--
-- The optional flag "Trunc => false" can be specified to prevent the polynomial being truncated to its equivalent form.
--
iWidth = {Trunc => true} >> opts -> (p, j) ->
(
  -- The bundle with extension class p is holomorphically equivalent
  -- to the bundle with extension class ptrimmed, which is obtained
  -- from p by throwing away all terms with
  -- 1 <= i <= 2j-2  and i-j+1 <= l <= j-1 ,
  -- where the general term is z^l u^i.
  ptrimmed := sum select(terms p,
                         i -> 1 <= (exponents i)_0_1                         and
                              (exponents i)_0_1 <= 2 * j - 2                 and
                              (exponents i)_0_1 - j + 1 <= (exponents i)_0_0 and
                              (exponents i)_0_0 <= j - 1
                        );
  if opts.Trunc then ( << "Working with the equivalent function p' = " << ptrimmed; print ""; p = ptrimmed; )
  else << "Should really work with the equivalent function p' = " << ptrimmed << ", but will not at your request."; print "";


  -- Get a generating set of sections and a list 
  -- or irrelevant coefficients
  N       := (polynomialParameters(p, j))_1;
  r       := doRing(p, j);
  allVars := r_0;
  pp      := r_1;
  aPoly   := r_2;
  bPoly   := r_3;
  irrel   := r_4;

  -- "fTv" is the first coordinate of the section on V
  fTv     := z^j * aPoly + pp * bPoly;

  -- Get the non-trivial relations between basic sections
  relsRes := getRelations(N, fTv);
  nonfree := relsRes_0 + irrel;
  rels    := relsRes_1;

  -- "templist" contains all the expressions that are zero
  relsList := flatten entries mingens rels;

  -- From that we make a list of substitutions...
  dosubst := apply(relsList, i -> { leadTerm i => (leadTerm i) - i } );

  print "Considering the following coefficients:";
  print allVars; print "";
  print "Found the following relations:";
  print relsList; print "";
  print "Will enact the following substitutions:";
  print dosubst; print "";
  
  -- ... and apply it to the polynomials:
  for i when i < #dosubst do
  (
    aPoly = substitute(aPoly, dosubst_i);
    bPoly = substitute(bPoly, dosubst_i);
  );

  -- Now form a list of generators from "allVars" by removing
  -- all zero coefficients:
  gb nonfree;
  changeables := apply(allVars, i -> i % nonfree);
  changeables  = (entries mingens ideal changeables)_0;

  -- "A" is a presentation of the module M according to Lemma 2.2
  A = setVectors(u^j * aPoly, u^j * bPoly, changeables, allVars);
  print "The module is:";
  print A;
  return(qLength A);
)


iHeight = (p, j) ->
(
  if (p == 0_Rbase) then return j*(j-1)/2;

  m := (polynomialParameters(p, j))_0;
  M := j * (j - 1) / 2;

  if (j > m + 1) then M = M - (j-m) * (j-m-1) / 2;

  return M;
)