///////////////////////////////////////////////////////////////////////////////
version="version schreyer.lib 4.1.2.0 Feb_2019 "; // $Id$
category="General purpose";
info="
LIBRARY: schreyer.lib helpers for derham.lib
AUTHOR:  Oleksandr Motsak <U@D>, where U={motsak}, D={mathematik.uni-kl.de}
KEYWORDS: Schreyer ordering; Schreyer resolution; syzygy
OVERVIEW:
The library contains several procedures for computing a/part of Schreyer
resolution (cf. [SFO]), and some helpers for derham.lib (which requires
resolutions over the homogenized Weyl algebra) for that purpose.
The input for any resolution computation is a set of vectors M in form of a
module over some basering R. The helpers works both in the commutative and
non-commutative setting (cf. [MO]), that is the ring R may be non-commutative,
in which case the ring ordering over it must be global. They produce/work with
partial Schreyer resolutions of (R^rank(M))/M in form of a specially constructed
ring (endowed with a special ring ordering that will be extended in the course
of a resolution computation) containing the following objects:
@* RES: the list of modules contains the images of maps (also called syzygy
modules) substituting the computed beginning of a Schreyer resolution, that is,
each syzygy module is given by a Groebner basis with respect to the
corresponding Schreyer ordering. RES starts with a zero map given by rank(M)
zero generators indicating that the image of the first differential map is
zero. The second map RES[2] is given by M, which indicates that the resolution
of (R^rank(M))/M is being computed.
@* MRES: the module is a direct sum of modules from RES and thus comprises all
computed differentials. Syzygies are shifted so that gen(i) is mapped to MRES[i]
under the differential map.
@* Here, we call a free resolution a Schreyer resolution if each syzygy
module is given by a Groebner basis with respect to the corresponding Schreyer
ordering. A Schreyer resolution can be much bigger than a minimal resolution of
the same module, but may be easier to construct. The Schreyer ordering
succesively extends the starting module ordering on M (defined in Singular by
the basering R) and is extended to higher syzygies using the following
definition:
@* a < b if and only if (d(a)<d(b)) OR ( (d(a)=d(b) AND (comp(a)<comp(b)) ),
@* where d(a) is the image of an under the differential (given by MRES), and
comp(a) is the module component, for any module terms a and b from the same
higher syzygy module.

NOTE: Since most comutations require the module syzextra.so, please be make sure
 to build it into Singular on Windows.

REFERENCES:
@*
[BMSS] Burcin, E., Motsak, O., Schreyer, F.-O., Steenpass, A.:
Refined algorithms to compute syzygies, 2015 (to appear).
@*
[SFO]  Schreyer, F.O.: Die Berechnung von Syzygien mit dem verallgemeinerten
Weierstrassschen Divisionssatz, Master's thesis, Univ. Hamburg, 1980.
@*
[MO]   Motsak, O.: Non-commutative Computer Algebra with applications:
Graded commutative algebra and related structures in Singular with applications,
Ph.D. thesis, TU Kaiserslautern, 2010.

PROCEDURES:
  Sres(M,l)    helper for computing Schreyer resolution
  Ssyz(M)      helper for computing Schreyer resolution of module M of length 1
  Scontinue(l) helper for extending currently active resolution

SEE ALSO: syz, sres, lres, res, fres
";

static proc prepareSyz( module I, list # )
{
  int i;
  int k = 0;
  int r = nrows(I);
  int c = ncols(I);


  if( size(#) > 0 )
  {
    if( typeof(#[1]) == "int" || typeof(#[1]) == "bigint" )
    {
      k = #[1];
    }
  }

  if( k < r )
  {
    "// *** Wrong k: ", k, " < nrows: ", r, " => setting k = r = ", r;
    k = r;
  }

//   "k: ", k;  "c: ", c;   "I: ", I;

  for( i = c; i > 0; i-- )
  {
    I[i] = I[i] + gen(k + i);
  }

  return(I);
}

static proc separateSyzGB( module J, int c )
{
  module II, G; vector v; int i;

  J = simplify(J, 2);

  for( i = ncols(J); i > 0; i-- )
  {
    v = J[i];
    if(   Syzextra::leadcomp(v) > c )
    {
      II[i] = v;
    } else
    {
      G[i] = v; // leave only gen(i): i <= c
    }
  }

  II = simplify(II, 2);
  G = simplify(G, 2);

  return (list(G, II));
}

static proc Sinit(module M)
{
  def @save = basering;

  int @DEBUG = 0; // !system("with", "ndebug");

  int @RANK = nrows(M); int @SIZE = ncols(M);

  int @IS_A_SB = attrib(M, "isSB"); // ??? only if all weights were zero?!

  if( !@IS_A_SB )
  {
    M = std(M); // this should be faster than computing std in S (later on)
  }

  def S =   Syzextra::MakeInducedSchreyerOrdering(1); // 1 puts history terms to the back
  // TODO: NOTE: +1 causes trouble to Singular interpreter!!!???
  setring S; // a new ring with a Schreyer ordering

  // Setup the leading syzygy^{-1} module to zero:
  module Z = 0; Z[@RANK] = 0; attrib(Z, "isHomog", intvec(0));

  module MRES = Z;

  list RES; RES[1] = Z;

  module F = freemodule(@RANK);
  intvec @V = deg(F[1..@RANK]);

  module M = imap(@save, M);

  attrib(M, "isHomog", @V);
  attrib(M, "isSB", 1);

  RES[size(RES)+1] = M; // list of all syzygy modules
  MRES = MRES, M;

  attrib(MRES, "isHomog", @V);

  attrib(S, "InducionLeads", lead(M));
  attrib(S, "InducionStart", @RANK);

  export RES;
  export MRES;
  return (S);
}

static proc Sstep()
{
  int @DEBUG = 0; // !system("with", "ndebug");

  // syzygy step:

/*
  // is initial weights are all zeroes!
  def L =  lead(M);
  intvec @V = deg(M[1..ncols(M)]);  @W;  @V;  @W = @V;  attrib(L, "isHomog", @W);
    Syzextra::SetInducedReferrence(L, @RANK, 0);
*/

//  def L =  lead(MRES);
//  @W = @W, @V;
//  attrib(L, "isHomog", @W);


  // General setting:
//    Syzextra::SetInducedReferrence(MRES, 0, 0); // limit: 0!
  int @l = size(RES);

  module M = RES[@l];

  module L = attrib(basering, "InducionLeads");
  int limit = attrib(basering, "InducionStart");

//  L;  limit;

  int @RANK = ncols(MRES) - ncols(M); // nrows(M); // what if M is zero?!

/*
  if( @RANK !=  nrows(M) )
  {
    type(MRES);
    @RANK;
    type(M);
    pause();
  }
*/

  intvec @W = attrib(M, "isHomog");
  intvec @V = deg(M[1..ncols(M)]);
  @V = @W, @V;

  Syzextra::SetInducedReferrence(L, limit, 0);

  def K = prepareSyz(M, @RANK);
//  K;

//  pause();

  K =   Syzextra::idPrepare(K, @RANK); // std(K); // ?
  K = simplify(K, 2);

//  K;

  module N = separateSyzGB(K, @RANK)[2]; // 1^st syz. module: vectors which start in lower part (comp >= @RANK)

//  basering; print(@V); type(N);
//  attrib(N, "isHomog", @V);  // TODO: fix "wrong weights"!!!? deg is wrong :(((
  N = std(N);
  attrib(N, "isHomog", @V);

  RES[@l + 1] = N; // list of all syzygy modules

  MRES = MRES, N;
  attrib(MRES, "isHomog", @V);


  L = L, lead(N);
  attrib(basering, "InducionLeads", L);

}

proc Scontinue(int l)
"USAGE:  Scontinue(int len)
RETURN:  nothing, instead it changes the currently active resolution
PURPOSE: extends the currently active resolution by at most len syzygies
ASSUME:  must be used within a ring returned by Sres or Ssyz
EXAMPLE: example Scontinue; shows an example
"
{
  def data =   Syzextra::GetInducedData();

  if( (!defined(RES)) || (!defined(MRES)) || (typeof(data) != "list") || (size(data) != 2) )
  {
    ERROR("Sorry, but basering does not seem to be returned by Sres or Ssyz");
  }
  for (;  (l != 0) && (size(RES[size(RES)]) > 0); l-- )
  {
    Sstep();
  }
}
example
{ "EXAMPLE:"; echo = 2;
  ring r;
  module M = maxideal(1); M;
  def S = Ssyz(M); setring S; S;
  "Only the first syzygy: ";
  RES; MRES;
  "More syzygies: ";
  Scontinue(10);
  RES; MRES;
}

proc Sres(module M, int l)
"USAGE:  Sres(module M, int len)
RETURN:  ring, containing a Schreyer resolution
PURPOSE: computes a Schreyer resolution of M of length at most len (see the library overview)
NOTE:    If given len is zero then nvars(basering) + 1 is used instead.
SEE ALSO: Ssyz
EXAMPLE: example Sres; shows an example
"
{
  def S = Sinit(M); setring S;

  if (l == 0)
  {
    l = nvars(basering) + 1; // not really an estimate...?!
  }

  Sstep(); l = l - 1;

  Scontinue(l);

  return (S);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r;
  module M = maxideal(1); M;
  def S = Sres(M, 0); setring S; S;
  RES;
  MRES;
  kill S;
  setring r; kill M;

  def A = nc_algebra(-1,0); setring A;
  ideal Q = var(1)^2, var(2)^2, var(3)^2;
  qring SCA = twostd(Q);
  basering;

  module M = maxideal(1);
  def S = Sres(M, 2); setring S; S;
  RES;
  MRES;
}

// ================================================================== //

static proc mod_init()
{
  load("syzextra.so");
}