hilb.cc File Reference

#include "kernel/mod2.h"
#include "misc/mylimits.h"
#include "misc/intvec.h"
#include "kernel/combinatorics/hilb.h"
#include "kernel/combinatorics/stairc.h"
#include "kernel/combinatorics/hutil.h"
#include "polys/monomials/ring.h"
#include "polys/monomials/p_polys.h"
#include "polys/simpleideals.h"
#include "kernel/ideals.h"
#include "polys/ext_fields/transext.h"
#include "coeffs/coeffs.h"
#include "kernel/linear_algebra/linearAlgebra.h"
#include "coeffs/numbers.h"
#include <vector>
#include "Singular/ipshell.h"
#include <ctime>
#include <iostream>

Go to the source code of this file.

Macros
#define	omsai   1

Functions
static int	hMinModulweight (intvec *modulweight)
static void	hHilbEst (scfmon stc, int Nstc, varset var, int Nvar)
static int *	hAddHilb (int Nv, int x, int *pol, int *lp)
static void	hLastHilb (scmon pure, int Nv, varset var, int *pol, int lp)
static void	hHilbStep (scmon pure, scfmon stc, int Nstc, varset var, int Nvar, int *pol, int Lpol)
static void	hWDegree (intvec *wdegree)
static void	SortByDeg_p (ideal I, poly p)
	!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
static ideal	SortByDeg (ideal I)
ideal	idQuotMon (ideal Iorig, ideal p)
static void	idAddMon (ideal I, ideal p)
static poly	ChoosePVar (ideal I)
static poly	ChoosePJL (ideal I)
static poly	ChooseP (ideal I)
static poly	SearchP (ideal I)
	searches for a monomial of degree d>=2 and divides it by a variable (result monomial of deg d-1)
static bool	JustVar (ideal I)
static void	eulerchar (ideal I, int variables, mpz_ptr ec)
static poly	SqFree (ideal I)
static bool	IsIn (poly p, ideal I)
static poly	LCMmon (ideal I)
void	rouneslice (ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int *&hilbpower)
void	slicehilb (ideal I)
static intvec *	hSeries (ideal S, intvec *modulweight, int, intvec *wdegree, ideal Q, ring tailRing)
intvec *	hHstdSeries (ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
intvec *	hFirstSeries (ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
intvec *	hSecondSeries (intvec *hseries1)
void	hDegreeSeries (intvec *s1, intvec *s2, int *co, int *mu)
static void	hPrintHilb (intvec *hseries, intvec *modul_weight)
void	hLookSeries (ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
static void	idInsertMonomial (ideal I, poly p)
static int	comapreMonoIdBases (ideal J, ideal Ob)
static int	CountOnIdUptoTruncationIndex (ideal I, int tr)
static int	comapreMonoIdBases_IG_Case (ideal J, int JCount, ideal Ob, int ObCount)
static int	positionInOrbit_IG_Case (ideal I, poly w, std::vector< ideal > idorb, std::vector< poly > polist, int trInd, int trunDegHs)
static int	positionInOrbit_FG_Case (ideal I, poly, std::vector< ideal > idorb, std::vector< poly >, int, int)
static int	positionInOrbitTruncationCase (ideal I, poly w, std::vector< ideal > idorb, std::vector< poly > polist, int, int trunDegHs)
static int	monCompare (const void *m, const void *n)
void	sortMonoIdeal_pCompare (ideal I)
static ideal	minimalMonomialGenSet (ideal I)
static poly	shiftInMon (poly p, int i, int lV, const ring r)
static poly	deleteInMon (poly w, int i, int lV, const ring r)
static void	TwordMap (poly p, poly w, int lV, int d, ideal Jwi, bool &flag)
static ideal	colonIdeal (ideal S, poly w, int lV, ideal Jwi, int trunDegHs)
void	HilbertSeries_OrbitData (ideal S, int lV, bool IG_CASE, bool mgrad, bool odp, int trunDegHs)
ideal	RightColonOperation (ideal S, poly w, int lV)

Variables
STATIC_VAR int **	Qpol
STATIC_VAR int *	Q0
STATIC_VAR int *	Ql
STATIC_VAR int	hLength

Macro Definition Documentation

omsai

#define omsai   1

Definition at line 28 of file hilb.cc.

Function Documentation

ChooseP()

static poly ChooseP ( ideal  I )

static

Definition at line 764 of file hilb.cc.

{
  poly m;
  //  TEST TO SEE WHICH ONE IS BETTER
  //m = ChoosePXL(I);
  //m = ChoosePXF(I);
  //m = ChoosePOL(I);
  //m = ChoosePOF(I);
  //m = ChoosePVL(I);
  //m = ChoosePVF(I);
  m = ChoosePJL(I);
  //m = ChoosePJF(I);
  return(m);
}

ChoosePJL()

static poly ChoosePJL ( ideal  I )

static

Definition at line 705 of file hilb.cc.

{
  int i,j,dummy;
  bool flag = TRUE;
  poly m = p_ISet(1,currRing);
  for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
  {
    flag = TRUE;
    for(j=1;(j<=currRing->N) && (flag);j++)
    {
      dummy = p_GetExp(I->m[i],j,currRing);
      if(dummy >= 2)
      {
        p_SetExp(m,j,dummy-1,currRing);
        p_Setm(m,currRing);
        flag = FALSE;
      }
    }
    if(!p_IsOne(m, currRing))
    {
      return(m);
    }
  }
  p_Delete(&m,currRing);
  m = ChoosePVar(I);
  return(m);
}

ChoosePVar()

static poly ChoosePVar ( ideal  I )

static

Definition at line 479 of file hilb.cc.

{
    bool flag=TRUE;
    int i,j;
    poly res;
    for(i=1;i<=currRing->N;i++)
    {
        flag=TRUE;
        for(j=IDELEMS(I)-1;(j>=0)&&(flag);j--)
        {
            if(p_GetExp(I->m[j], i, currRing)>0)
            {
                flag=FALSE;
            }
        }
 
        if(flag == TRUE)
        {
            res = p_ISet(1, currRing);
            p_SetExp(res, i, 1, currRing);
            p_Setm(res,currRing);
            return(res);
        }
        else
        {
            p_Delete(&res, currRing);
        }
    }
    return(NULL); //i.e. it is the maximal ideal
}

colonIdeal()

static ideal colonIdeal ( ideal  S, poly  w, int  lV, ideal  Jwi, int  trunDegHs  )

static

Definition at line 1939 of file hilb.cc.

{
  /*
   * It computes the right colon ideal of a two-sided ideal S
   * w.r.t. word w and save it in a new object Jwi.
   * It keeps S and w unchanged.
   */
 
  if(idIs0(S))
  {
    return(S);
  }
 
  int i, d;
  d = p_Totaldegree(w, currRing);
  if(trunDegHs !=0 && d >= trunDegHs)
  {
   idInsertMonomial(Jwi, p_One(currRing));
   return(Jwi);
  }
  bool flag = FALSE;
  int SCount = IDELEMS(S);
  for(i = 0; i < SCount; i++)
  {
    TwordMap(S->m[i], w, lV, d, Jwi, flag);
    if(flag)
    {
      break;
    }
  }
 
  Jwi = minimalMonomialGenSet(Jwi);
  return(Jwi);
}

comapreMonoIdBases()

static int comapreMonoIdBases ( ideal  J, ideal  Ob  )

static

Definition at line 1483 of file hilb.cc.

{
  /*
   * Monomials of J and Ob are assumed to
   * be already sorted. J and Ob are
   * represented by the minimal generating set.
   */
  int i, s;
  s = 1;
  int JCount = IDELEMS(J);
  int ObCount = IDELEMS(Ob);
 
  if(idIs0(J))
  {
    return(1);
  }
  if(JCount != ObCount)
  {
    return(0);
  }
 
  for(i = 0; i < JCount; i++)
  {
    if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
    {
      return(0);
    }
  }
  return(s);
}

comapreMonoIdBases_IG_Case()

static int comapreMonoIdBases_IG_Case ( ideal  J, int  JCount, ideal  Ob, int  ObCount  )

static

Definition at line 1541 of file hilb.cc.

{
  /*
   * Monomials of J and Ob are assumed to
   * be already sorted in increasing degrees.
   * J and Ob are represented by the minimal
   * generating set.  It checks if J and Ob have
   * same monomials  up to deg <=tr.
   */
 
  int i, s;
  s = 1;
  //when J is null
  //
  if(JCount != ObCount)
  {
    return(0);
  }
 
  if(JCount == 0)
  {
    return(1);
  }
 
  for(i = 0; i< JCount; i++)
  {
    if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
    {
      return(0);
    }
  }
 
  return(s);
}

CountOnIdUptoTruncationIndex()

static int CountOnIdUptoTruncationIndex ( ideal  I, int  tr  )

static

Definition at line 1514 of file hilb.cc.

{
  /*
   * The ideal I must be sorted in increasing total degree.
   * It counts the number of monomials in I up to
   * degree less than or equal to tr.
   */
 
  //case when I=1;
  if(p_Totaldegree(I->m[0], currRing) == 0)
  {
    return(1);
  }
 
  int count = 0;
  for(int i = 0; i < IDELEMS(I); i++)
  {
    if(p_Totaldegree(I->m[i], currRing) > tr)
    {
      return (count);
    }
    count = count + 1;
  }
 
  return(count);
}

deleteInMon()

static poly deleteInMon ( poly  w, int  i, int  lV, const ring  r  )

static

Definition at line 1844 of file hilb.cc.

{
  /*
   * deletes the variables upto i^th layer of monomial w
   * w remains unchanged
   * creates new poly and returns it for the colon ideal
   */
 
  poly dw = p_One(currRing);
  int *e = (int *)omAlloc((r->N+1)*sizeof(int));
  int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
  p_GetExpV(w, e, r);
  int j, cnt;
  cnt = i*lV;
  /*
  for(j=1;j<=cnt;j++)
  {
    e[j]=0;
  }*/
  for(j = (cnt+1); j < (r->N+1); j++)
  {
    s[j] = e[j];
  }
 
  p_SetExpV(dw, s, currRing);//new exponents
  omFree(e);
  omFree(s);
 
  p_SetComp(dw, p_GetComp(w, currRing), currRing);
  p_Setm(dw, currRing);
 
  return(dw);
}

eulerchar()

static void eulerchar ( ideal  I, int  variables, mpz_ptr  ec  )

static

Definition at line 837 of file hilb.cc.

{
  loop
  {
    mpz_t dummy;
    if(JustVar(I) == TRUE)
    {
      if(IDELEMS(I) == variables)
      {
        mpz_init(dummy);
        if((variables % 2) == 0)
          mpz_set_ui(dummy, 1);
        else
          mpz_set_si(dummy, -1);
        mpz_add(ec, ec, dummy);
        mpz_clear(dummy);
      }
      return;
    }
    ideal p = idInit(1,1);
    p->m[0] = SearchP(I);
    //idPrint(I);
    //idPrint(p);
    //printf("\nNow get in idQuotMon\n");
    ideal Ip = idQuotMon(I,p);
    //idPrint(Ip);
    //Ip = SortByDeg(Ip);
    int i,howmanyvarinp = 0;
    for(i = 1;i<=currRing->N;i++)
    {
      if(p_GetExp(p->m[0],i,currRing)>0)
      {
        howmanyvarinp++;
      }
    }
    eulerchar(Ip, variables-howmanyvarinp, ec);
    id_Delete(&Ip, currRing);
    idAddMon(I,p);
    id_Delete(&p, currRing);
  }
}

hAddHilb()

static int * hAddHilb ( int  Nv, int  x, int *  pol, int *  lp  )

static

Definition at line 104 of file hilb.cc.

{
  int  l = *lp, ln, i;
  int  *pon;
  *lp = ln = l + x;
  pon = Qpol[Nv];
  memcpy(pon, pol, l * sizeof(int));
  if (l > x)
  {/*pon[i] -= pol[i - x];*/
    for (i = x; i < l; i++)
    { int64 t=pon[i];
      int64 t2=pol[i - x];
      t-=t2;
      if ((t>=INT_MIN)&&(t<=INT_MAX)) pon[i]=t;
      else if (!errorreported) WerrorS("int overflow in hilb 1");
    }
    for (i = l; i < ln; i++)
    { /*pon[i] = -pol[i - x];*/
      int64 t= -pol[i - x];
      if ((t>=INT_MIN)&&(t<=INT_MAX)) pon[i]=t;
      else if (!errorreported) WerrorS("int overflow in hilb 2");
    }
  }
  else
  {
    for (i = l; i < x; i++)
      pon[i] = 0;
    for (i = x; i < ln; i++)
      pon[i] = -pol[i - x];
  }
  return pon;
}

hDegreeSeries()

void hDegreeSeries	(	intvec *	s1,
		intvec *	s2,
		int *	co,
		int *	mu
	)

Definition at line 1380 of file hilb.cc.

{
  int i, j, k;
  int m;
  *co = *mu = 0;
  if ((s1 == NULL) || (s2 == NULL))
    return;
  i = s1->length();
  j = s2->length();
  if (j > i)
    return;
  m = 0;
  for(k=j-2; k>=0; k--)
    m += (*s2)[k];
  *mu = m;
  *co = i - j;
}

hFirstSeries()

intvec * hFirstSeries	(	ideal	S,
		intvec *	modulweight,
		ideal	Q,
		intvec *	wdegree,
		ring	tailRing
	)

Definition at line 1335 of file hilb.cc.

{
  id_TestTail(S, currRing, tailRing);
  if (Q!= NULL) id_TestTail(Q, currRing, tailRing);
 
  intvec *hseries1= hSeries(S, modulweight, 1, wdegree, Q, tailRing);
  if (errorreported) { delete hseries1; hseries1=NULL; }
  return hseries1;
}

hHilbEst()

static void hHilbEst ( scfmon  stc, int  Nstc, varset  var, int  Nvar  )

static

Definition at line 63 of file hilb.cc.

{
  int  i, j;
  int  x, y, z = 1;
  int  *p;
  for (i = Nvar; i>0; i--)
  {
    x = 0;
    for (j = 0; j < Nstc; j++)
    {
      y = stc[j][var[i]];
      if (y > x)
        x = y;
    }
    z += x;
    j = i - 1;
    if (z > Ql[j])
    {
      if (z>(MAX_INT_VAL)/2)
      {
       WerrorS("internal arrays too big");
       return;
      }
      p = (int *)omAlloc((unsigned long)z * sizeof(int));
      if (Ql[j]!=0)
      {
        if (j==0)
          memcpy(p, Qpol[j], Ql[j] * sizeof(int));
        omFreeSize((ADDRESS)Qpol[j], Ql[j] * sizeof(int));
      }
      if (j==0)
      {
        for (x = Ql[j]; x < z; x++)
          p[x] = 0;
      }
      Ql[j] = z;
      Qpol[j] = p;
    }
  }
}

hHilbStep()

static void hHilbStep ( scmon  pure, scfmon  stc, int  Nstc, varset  var, int  Nvar, int *  pol, int  Lpol  )

static

Definition at line 177 of file hilb.cc.

{
  int  iv = Nvar -1, ln, a, a0, a1, b, i;
  int  x, x0;
  scmon pn;
  scfmon sn;
  int  *pon;
  if (Nstc==0)
  {
    hLastHilb(pure, iv, var, pol, Lpol);
    return;
  }
  x = a = 0;
  pn = hGetpure(pure);
  sn = hGetmem(Nstc, stc, stcmem[iv]);
  hStepS(sn, Nstc, var, Nvar, &a, &x);
  Q0[iv] = Q0[Nvar];
  ln = Lpol;
  pon = pol;
  if (a == Nstc)
  {
    x = pure[var[Nvar]];
    if (x!=0)
      pon = hAddHilb(iv, x, pon, &ln);
    hHilbStep(pn, sn, a, var, iv, pon, ln);
    return;
  }
  else
  {
    pon = hAddHilb(iv, x, pon, &ln);
    hHilbStep(pn, sn, a, var, iv, pon, ln);
  }
  b = a;
  x0 = 0;
  loop
  {
    Q0[iv] += (x - x0);
    a0 = a;
    x0 = x;
    hStepS(sn, Nstc, var, Nvar, &a, &x);
    hElimS(sn, &b, a0, a, var, iv);
    a1 = a;
    hPure(sn, a0, &a1, var, iv, pn, &i);
    hLex2S(sn, b, a0, a1, var, iv, hwork);
    b += (a1 - a0);
    ln = Lpol;
    if (a < Nstc)
    {
      pon = hAddHilb(iv, x - x0, pol, &ln);
      hHilbStep(pn, sn, b, var, iv, pon, ln);
    }
    else
    {
      x = pure[var[Nvar]];
      if (x!=0)
        pon = hAddHilb(iv, x - x0, pol, &ln);
      else
        pon = pol;
      hHilbStep(pn, sn, b, var, iv, pon, ln);
      return;
    }
  }
}

hHstdSeries()

intvec * hHstdSeries	(	ideal	S,
		intvec *	modulweight,
		intvec *	wdegree,
		ideal	Q,
		ring	tailRing
	)

Definition at line 1328 of file hilb.cc.

{
  id_TestTail(S, currRing, tailRing);
  if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
  return hSeries(S, modulweight, 0, wdegree, Q, tailRing);
}

HilbertSeries_OrbitData()

void HilbertSeries_OrbitData	(	ideal	S,
		int	lV,
		bool	IG_CASE,
		bool	mgrad,
		bool	odp,
		int	trunDegHs
	)

Definition at line 1974 of file hilb.cc.

{
        
        /* new story: 
        no lV is needed, i.e.  it is to be determined
        the rest is extracted from the interface input list in extra.cc and makes the input of this proc
        called from extra.cc
        */
        
  /*
   * This is based on iterative right colon operations on a
   * two-sided monomial ideal of the free associative algebra.
   * The algorithm terminates for those monomial ideals
   * whose monomials define  "regular formal languages",
   * that is, all monomials of the input ideal can be obtained
   * from finite languages by applying finite number of
   * rational operations.
   */
 
  int trInd;
  S = minimalMonomialGenSet(S);
  if( !idIs0(S) && p_Totaldegree(S->m[0], currRing)==0)
  {
    PrintS("Hilbert Series:\n 0\n");
    return;
  }
  int (*POS)(ideal, poly, std::vector<ideal>, std::vector<poly>, int, int);
  if(trunDegHs != 0)
  {
    Print("\nTruncation degree = %d\n",trunDegHs);
    POS=&positionInOrbitTruncationCase;
  }
  else
  {
    if(IG_CASE)
    {
      if(idIs0(S))
      {
        WerrorS("wrong input: it is not an infinitely gen. case");
        return;
      }
      trInd = p_Totaldegree(S->m[IDELEMS(S)-1], currRing);
      POS = &positionInOrbit_IG_Case;
    }
    else
    POS = &positionInOrbit_FG_Case;
  }
  std::vector<ideal > idorb;
  std::vector< poly > polist;
 
  ideal orb_init = idInit(1, 1);
  idorb.push_back(orb_init);
 
  polist.push_back( p_One(currRing));
 
  std::vector< std::vector<int> > posMat;
  std::vector<int> posRow(lV,0);
  std::vector<int> C;
 
  int ds, is, ps;
  unsigned long lpcnt = 0;
 
  poly w, wi;
  ideal Jwi;
 
  while(lpcnt < idorb.size())
  {
    w = NULL;
    w = polist[lpcnt];
    if(lpcnt >= 1 && idIs0(idorb[lpcnt]) == FALSE)
    {
      if(p_Totaldegree(idorb[lpcnt]->m[0], currRing) != 0)
      {
        C.push_back(1);
      }
      else
        C.push_back(0);
    }
    else
    {
      C.push_back(1);
    }
 
    ds = p_Totaldegree(w, currRing);
    lpcnt++;
 
    for(is = 1; is <= lV; is++)
    {
      wi = NULL;
      //make new copy 'wi' of word w=polist[lpcnt]
      //and update it (for the colon operation).
      //if corresponding to wi, right colon operation gives
      //a new (right colon) ideal of S,
      //keep 'wi' in the polist else delete it
 
      wi = pCopy(w);
      p_SetExp(wi, (ds*lV)+is, 1, currRing);
      p_Setm(wi, currRing);
      Jwi = NULL;
      //Jwi stores (right) colon ideal of S w.r.t. word
      //wi if colon operation gives a new ideal place it
      //in the vector of ideals  'idorb'
      //otherwise delete it
 
      Jwi = idInit(1,1);
 
      Jwi = colonIdeal(S, wi, lV, Jwi, trunDegHs);
      ps = (*POS)(Jwi, wi, idorb, polist, trInd, trunDegHs);
 
      if(ps == 0)  // finds a new ideal
      {
        posRow[is-1] = idorb.size();
 
        idorb.push_back(Jwi);
        polist.push_back(wi);
      }
      else // ideal is already there  in the set
      {
        posRow[is-1]=ps-1;
        idDelete(&Jwi);
        pDelete(&wi);
      }
    }
    posMat.push_back(posRow);
    posRow.resize(lV,0);
  }
  int lO = C.size();//size of the orbit
  PrintLn();
  Print("maximal length of words = %ld\n", p_Totaldegree(polist[lO-1], currRing));
  Print("\nlength of the Orbit = %d", lO);
  PrintLn();
 
  if(odp)
  {
    Print("words description of the Orbit: \n");
    for(is = 0; is < lO; is++)
    {
      pWrite0(polist[is]);
      PrintS("    ");
    }
    PrintLn();
    PrintS("\nmaximal degree,  #(sum_j R(w,w_j))");
    PrintLn();
    for(is = 0; is < lO; is++)
    {
      if(idIs0(idorb[is]))
      {
        PrintS("NULL\n");
      }
      else
      {
        Print("%ld,  %d \n",p_Totaldegree(idorb[is]->m[IDELEMS(idorb[is])-1], currRing),IDELEMS(idorb[is]));
      }
    }
  }
 
  for(is = idorb.size()-1; is >= 0; is--)
  {
    idDelete(&idorb[is]);
  }
  for(is = polist.size()-1; is >= 0; is--)
  {
    pDelete(&polist[is]);
  }
 
  idorb.resize(0);
  polist.resize(0);
 
  int adjMatrix[lO][lO];
  memset(adjMatrix, 0, lO*lO*sizeof(int));
  int rowCount, colCount;
  int tm = 0;
  if(!mgrad)
  {
    for(rowCount = 0; rowCount < lO; rowCount++)
    {
      for(colCount = 0; colCount < lV; colCount++)
      {
        tm = posMat[rowCount][colCount];
        adjMatrix[rowCount][tm] = adjMatrix[rowCount][tm] + 1;
      }
    }
  }
 
  ring r = currRing;
  int npar;
  char** tt;
  TransExtInfo p;
  if(!mgrad)
  {
    tt=(char**)omAlloc(sizeof(char*));
    tt[0] = omStrDup("t");
    npar = 1;
  }
  else
  {
    tt=(char**)omalloc(lV*sizeof(char*));
    for(is = 0; is < lV; is++)
    {
      tt[is] = (char*)omAlloc(7*sizeof(char)); //if required enlarge it later
      sprintf (tt[is], "t%d", is+1);
    }
    npar = lV;
  }
 
  p.r = rDefault(0, npar, tt);
  coeffs cf = nInitChar(n_transExt, &p);
  char** xx = (char**)omAlloc(sizeof(char*));
  xx[0] = omStrDup("x");
  ring R = rDefault(cf, 1, xx);
  rChangeCurrRing(R);//rWrite(R);
  /*
   * matrix corresponding to the orbit of the ideal
   */
  matrix mR = mpNew(lO, lO);
  matrix cMat = mpNew(lO,1);
  poly rc;
 
  if(!mgrad)
  {
    for(rowCount = 0; rowCount < lO; rowCount++)
    {
      for(colCount = 0; colCount < lO; colCount++)
      {
        if(adjMatrix[rowCount][colCount] != 0)
        {
          MATELEM(mR, rowCount + 1, colCount + 1) = p_ISet(adjMatrix[rowCount][colCount], R);
          p_SetCoeff(MATELEM(mR, rowCount + 1, colCount + 1), n_Mult(pGetCoeff(mR->m[lO*rowCount+colCount]),n_Param(1, R->cf), R->cf), R);
        }
      }
    }
  }
  else
  {
     for(rowCount = 0; rowCount < lO; rowCount++)
     {
       for(colCount = 0; colCount < lV; colCount++)
       {
          rc=NULL;
          rc=p_One(R);
          p_SetCoeff(rc, n_Mult(pGetCoeff(rc), n_Param(colCount+1, R->cf),R->cf), R);
          MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1)=p_Add_q(rc,MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1), R);
       }
     }
  }
 
  for(rowCount = 0; rowCount < lO; rowCount++)
  {
    if(C[rowCount] != 0)
    {
      MATELEM(cMat, rowCount + 1, 1) = p_ISet(C[rowCount], R);
    }
  }
 
  matrix u;
  unitMatrix(lO, u); //unit matrix
  matrix gMat = mp_Sub(u, mR, R);
 
  char* s;
 
  if(odp)
  {
    PrintS("\nlinear system:\n");
    if(!mgrad)
    {
      for(rowCount = 0; rowCount < lO; rowCount++)
      {
        Print("H(%d) = ", rowCount+1);
        for(colCount = 0; colCount < lV; colCount++)
        {
          StringSetS(""); nWrite(n_Param(1, R->cf));
          s = StringEndS(); PrintS(s);
          Print("*"); omFree(s);
          Print("H(%d) + ", posMat[rowCount][colCount] + 1);
        }
        Print(" %d\n", C[rowCount] );
      }
      PrintS("where H(1) represents the series corresp. to input ideal\n");
      PrintS("and i^th summand in the rhs of an eqn. is according\n");
      PrintS("to the right colon map corresp. to the i^th variable\n");
    }
    else
    {
      for(rowCount = 0; rowCount < lO; rowCount++)
      {
        Print("H(%d) = ", rowCount+1);
        for(colCount = 0; colCount < lV; colCount++)
        {
           StringSetS(""); nWrite(n_Param(colCount+1, R->cf));
           s = StringEndS(); PrintS(s);
           Print("*");omFree(s);
           Print("H(%d) + ", posMat[rowCount][colCount] + 1);
        }
        Print(" %d\n", C[rowCount] );
      }
      PrintS("where H(1) represents the series corresp. to input ideal\n");
    }
  }
  PrintLn();
  posMat.resize(0);
  C.resize(0);
  matrix pMat;
  matrix lMat;
  matrix uMat;
  matrix H_serVec = mpNew(lO, 1);
  matrix Hnot;
 
  //std::clock_t start;
  //start = std::clock();
 
  luDecomp(gMat, pMat, lMat, uMat, R);
  luSolveViaLUDecomp(pMat, lMat, uMat, cMat, H_serVec, Hnot);
 
  //to print system solving time
  //if(odp){
  //std::cout<<"solving time of the system = "<<(std::clock()-start)/(double)(CLOCKS_PER_SEC / 1000)<<" ms"<<std::endl;}
 
  mp_Delete(&mR, R);
  mp_Delete(&u, R);
  mp_Delete(&pMat, R);
  mp_Delete(&lMat, R);
  mp_Delete(&uMat, R);
  mp_Delete(&cMat, R);
  mp_Delete(&gMat, R);
  mp_Delete(&Hnot, R);
  //print the Hilbert series and length of the Orbit
  PrintLn();
  Print("Hilbert series:");
  PrintLn();
  pWrite(H_serVec->m[0]);
  if(!mgrad)
  {
    omFree(tt[0]);
  }
  else
  {
    for(is = lV-1; is >= 0; is--)
 
      omFree( tt[is]);
  }
  omFree(tt);
  omFree(xx[0]);
  omFree(xx);
  rChangeCurrRing(r);
  rKill(R);
}

hLastHilb()

static void hLastHilb ( scmon  pure, int  Nv, varset  var, int *  pol, int  lp  )

static

Definition at line 137 of file hilb.cc.

{
  int  l = lp, x, i, j;
  int  *pl;
  int  *p;
  p = pol;
  for (i = Nv; i>0; i--)
  {
    x = pure[var[i + 1]];
    if (x!=0)
      p = hAddHilb(i, x, p, &l);
  }
  pl = *Qpol;
  j = Q0[Nv + 1];
  for (i = 0; i < l; i++)
  { /* pl[i + j] += p[i];*/
    int64 t=pl[i+j];
    int64 t2=p[i];
    t+=t2;
    if ((t>=INT_MIN)&&(t<=INT_MAX)) pl[i+j]=t;
    else if (!errorreported) WerrorS("int overflow in hilb 3");
  }
  x = pure[var[1]];
  if (x!=0)
  {
    j += x;
    for (i = 0; i < l; i++)
    { /* pl[i + j] -= p[i];*/
      int64 t=pl[i+j];
      int64 t2=p[i];
      t-=t2;
      if ((t>=INT_MIN)&&(t<=INT_MAX)) pl[i+j]=t;
      else if (!errorreported) WerrorS("int overflow in hilb 4");
    }
  }
  j += l;
  if (j > hLength)
    hLength = j;
}

hLookSeries()

void hLookSeries	(	ideal	S,
		intvec *	modulweight,
		ideal	Q,
		intvec *	wdegree,
		ring	tailRing
	)

Definition at line 1424 of file hilb.cc.

{
  id_TestTail(S, currRing, tailRing);
 
  intvec *hseries1 = hFirstSeries(S, modulweight, Q, wdegree, tailRing);
  if (errorreported) return;
 
  hPrintHilb(hseries1,modulweight);
 
  const int l = hseries1->length()-1;
 
  intvec *hseries2 = (l > 1) ? hSecondSeries(hseries1) : hseries1;
 
  int co, mu;
  hDegreeSeries(hseries1, hseries2, &co, &mu);
 
  PrintLn();
  hPrintHilb(hseries2,modulweight);
  if ((l == 1) &&(mu == 0))
    scPrintDegree(rVar(currRing)+1, 0);
  else
    scPrintDegree(co, mu);
  if (l>1)
    delete hseries1;
  delete hseries2;
}

hMinModulweight()

static int hMinModulweight ( intvec *  modulweight )

static

Definition at line 49 of file hilb.cc.

{
  int i,j,k;
 
  if(modulweight==NULL) return 0;
  j=(*modulweight)[0];
  for(i=modulweight->rows()-1;i!=0;i--)
  {
    k=(*modulweight)[i];
    if(k<j) j=k;
  }
  return j;
}

hPrintHilb()

static void hPrintHilb ( intvec *  hseries, intvec *  modul_weight  )

static

Definition at line 1398 of file hilb.cc.

{
  int  i, j, l, k;
  if (hseries == NULL)
    return;
  l = hseries->length()-1;
  k = (*hseries)[l];
  if ((modul_weight!=NULL)&&(modul_weight->compare(0)!=0))
  {
    char *s=modul_weight->ivString(1,0,1);
    Print("module weights:%s\n",s);
    omFree(s);
  }
  for (i = 0; i < l; i++)
  {
    j = (*hseries)[i];
    if (j != 0)
    {
      Print("//  %8d t^%d\n", j, i+k);
    }
  }
}

hSecondSeries()

intvec * hSecondSeries

(

intvec *

hseries1

)

Definition at line 1345 of file hilb.cc.

{
  intvec *work, *hseries2;
  int i, j, k, t, l;
  int s;
  if (hseries1 == NULL)
    return NULL;
  work = new intvec(hseries1);
  k = l = work->length()-1;
  s = 0;
  for (i = k-1; i >= 0; i--)
    s += (*work)[i];
  loop
  {
    if ((s != 0) || (k == 1))
      break;
    s = 0;
    t = (*work)[k-1];
    k--;
    for (i = k-1; i >= 0; i--)
    {
      j = (*work)[i];
      (*work)[i] = -t;
      s += t;
      t += j;
    }
  }
  hseries2 = new intvec(k+1);
  for (i = k-1; i >= 0; i--)
    (*hseries2)[i] = (*work)[i];
  (*hseries2)[k] = (*work)[l];
  delete work;
  return hseries2;
}

hSeries()

static intvec * hSeries ( ideal  S, intvec *  modulweight, int  , intvec *  wdegree, ideal  Q, ring  tailRing  )

static

Definition at line 1172 of file hilb.cc.

{
//  id_TestTail(S, currRing, tailRing);
 
  intvec *work, *hseries1=NULL;
  int  mc;
  int  p0;
  int  i, j, k, l, ii, mw;
  hexist = hInit(S, Q, &hNexist, tailRing);
  if (hNexist==0)
  {
    hseries1=new intvec(2);
    (*hseries1)[0]=1;
    (*hseries1)[1]=0;
    return hseries1;
  }
 
  #if 0
  if (wdegree == NULL)
    hWeight();
  else
    hWDegree(wdegree);
  #else
  if (wdegree != NULL) hWDegree(wdegree);
  #endif
 
  p0 = 1;
  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
  stcmem = hCreate((currRing->N) - 1);
  Qpol = (int **)omAlloc(((currRing->N) + 1) * sizeof(int *));
  Ql = (int *)omAlloc0(((currRing->N) + 1) * sizeof(int));
  Q0 = (int *)omAlloc(((currRing->N) + 1) * sizeof(int));
  *Qpol = NULL;
  hLength = k = j = 0;
  mc = hisModule;
  if (mc!=0)
  {
    mw = hMinModulweight(modulweight);
    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
  }
  else
  {
    mw = 0;
    hstc = hexist;
    hNstc = hNexist;
  }
  loop
  {
    if (mc!=0)
    {
      hComp(hexist, hNexist, mc, hstc, &hNstc);
      if (modulweight != NULL)
        j = (*modulweight)[mc-1]-mw;
    }
    if (hNstc!=0)
    {
      hNvar = (currRing->N);
      for (i = hNvar; i>=0; i--)
        hvar[i] = i;
      //if (notstc) // TODO: no mon divides another
        hStaircase(hstc, &hNstc, hvar, hNvar);
      hSupp(hstc, hNstc, hvar, &hNvar);
      if (hNvar!=0)
      {
        if ((hNvar > 2) && (hNstc > 10))
          hOrdSupp(hstc, hNstc, hvar, hNvar);
        hHilbEst(hstc, hNstc, hvar, hNvar);
        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
        hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
        hLexS(hstc, hNstc, hvar, hNvar);
        Q0[hNvar] = 0;
        hHilbStep(hpure, hstc, hNstc, hvar, hNvar, &p0, 1);
      }
    }
    else
    {
      if(*Qpol!=NULL)
        (**Qpol)++;
      else
      {
        *Qpol = (int *)omAlloc(sizeof(int));
        hLength = *Ql = **Qpol = 1;
      }
    }
    if (*Qpol!=NULL)
    {
      i = hLength;
      while ((i > 0) && ((*Qpol)[i - 1] == 0))
        i--;
      if (i > 0)
      {
        l = i + j;
        if (l > k)
        {
          work = new intvec(l);
          for (ii=0; ii<k; ii++)
            (*work)[ii] = (*hseries1)[ii];
          if (hseries1 != NULL)
            delete hseries1;
          hseries1 = work;
          k = l;
        }
        while (i > 0)
        {
          (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
          (*Qpol)[i - 1] = 0;
          i--;
        }
      }
    }
    mc--;
    if (mc <= 0)
      break;
  }
  if (k==0)
  {
    hseries1=new intvec(2);
    (*hseries1)[0]=0;
    (*hseries1)[1]=0;
  }
  else
  {
    l = k+1;
    while ((*hseries1)[l-2]==0) l--;
    if (l!=k)
    {
      work = new intvec(l);
      for (ii=l-2; ii>=0; ii--)
        (*work)[ii] = (*hseries1)[ii];
      delete hseries1;
      hseries1 = work;
    }
    (*hseries1)[l-1] = mw;
  }
  for (i = 0; i <= (currRing->N); i++)
  {
    if (Ql[i]!=0)
      omFreeSize((ADDRESS)Qpol[i], Ql[i] * sizeof(int));
  }
  omFreeSize((ADDRESS)Q0, ((currRing->N) + 1) * sizeof(int));
  omFreeSize((ADDRESS)Ql, ((currRing->N) + 1) * sizeof(int));
  omFreeSize((ADDRESS)Qpol, ((currRing->N) + 1) * sizeof(int *));
  hKill(stcmem, (currRing->N) - 1);
  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
  hDelete(hexist, hNexist);
  if (hisModule!=0)
    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
  return hseries1;
}

hWDegree()

static void hWDegree ( intvec *  wdegree )

static

Definition at line 245 of file hilb.cc.

{
  int i, k;
  int x;
 
  for (i=(currRing->N); i; i--)
  {
    x = (*wdegree)[i-1];
    if (x != 1)
    {
      for (k=hNexist-1; k>=0; k--)
      {
        hexist[k][i] *= x;
      }
    }
  }
}

idAddMon()

static void idAddMon ( ideal  I, ideal  p  )

static

Definition at line 471 of file hilb.cc.

{
  SortByDeg_p(I,p->m[0]);
  p->m[0]=NULL; // is now in I
  //idSkipZeroes(I);
}

idInsertMonomial()

static void idInsertMonomial ( ideal  I, poly  p  )

static

Definition at line 1456 of file hilb.cc.

{
  /*
   * It adds monomial in I and if required,
   * enlarge the size of poly-set by 16.
   * It does not make copy of  p.
   */
 
  if(I == NULL)
  {
    return;
  }
 
  int j = IDELEMS(I) - 1;
  while ((j >= 0) && (I->m[j] == NULL))
  {
    j--;
  }
  j++;
  if (j == IDELEMS(I))
  {
    pEnlargeSet(&(I->m), IDELEMS(I), 16);
    IDELEMS(I) +=16;
  }
  I->m[j] = p;
}

idQuotMon()

ideal idQuotMon	(	ideal	Iorig,
		ideal	p
	)

Definition at line 409 of file hilb.cc.

{
    if(idIs0(Iorig))
    {
      ideal res = idInit(1,1);
      res->m[0] = poly(0);
      return(res);
    }
    if(idIs0(p))
    {
      ideal res = idInit(1,1);
      res->m[0] = pOne();
      return(res);
    }
    ideal I = id_Head(Iorig,currRing);
    ideal res = idInit(IDELEMS(I),1);
    int i,j;
    int dummy;
    for(i = 0; i<IDELEMS(I); i++)
    {
      res->m[i] = p_Head(I->m[i], currRing);
      for(j = 1; (j<=currRing->N) ; j++)
      {
        dummy = p_GetExp(p->m[0], j, currRing);
        if(dummy > 0)
        {
          if(p_GetExp(I->m[i], j, currRing) < dummy)
          {
            p_SetExp(res->m[i], j, 0, currRing);
          }
          else
          {
            p_SetExp(res->m[i], j, p_GetExp(I->m[i], j, currRing) - dummy, currRing);
          }
        }
      }
      p_Setm(res->m[i], currRing);
      if(p_Totaldegree(res->m[i],currRing) == p_Totaldegree(I->m[i],currRing))
      {
        p_Delete(&res->m[i],currRing);
      }
      else
      {
        p_Delete(&I->m[i],currRing);
      }
    }
    idSkipZeroes(res);
    idSkipZeroes(I);
    if(!idIs0(res))
    {
      for(i = 0; i<=IDELEMS(res)-1; i++)
      {
        SortByDeg_p(I,res->m[i]);
        res->m[i]=NULL; // is now in I
      }
    }
    id_Delete(&res,currRing);
    //idDegSortTest(I);
    return(I);
}

IsIn()

static bool IsIn ( poly  p, ideal  I  )

static

Definition at line 909 of file hilb.cc.

{
  //assumes that I is ordered by degree
  if(idIs0(I))
  {
    if(p==poly(0))
    {
      return(TRUE);
    }
    else
    {
      return(FALSE);
    }
  }
  if(p==poly(0))
  {
    return(FALSE);
  }
  int i,j;
  bool flag;
  for(i = 0;i<IDELEMS(I);i++)
  {
    flag = TRUE;
    for(j = 1;(j<=currRing->N) &&(flag);j++)
    {
      if(p_GetExp(p, j, currRing)<p_GetExp(I->m[i], j, currRing))
      {
        flag = FALSE;
      }
    }
    if(flag)
    {
      return(TRUE);
    }
  }
  return(FALSE);
}

JustVar()

static bool JustVar ( ideal  I )

static

Definition at line 806 of file hilb.cc.

{
    #if 0
    int i,j;
    bool foundone;
    for(i=0;i<=IDELEMS(I)-1;i++)
    {
        foundone = FALSE;
        for(j = 1;j<=currRing->N;j++)
        {
            if(p_GetExp(I->m[i], j, currRing)>0)
            {
                if(foundone == TRUE)
                {
                    return(FALSE);
                }
                foundone = TRUE;
            }
        }
    }
    return(TRUE);
    #else
    if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)>1)
    {
        return(FALSE);
    }
    return(TRUE);
    #endif
}

LCMmon()

static poly LCMmon ( ideal  I )

static

Definition at line 948 of file hilb.cc.

{
  if(idIs0(I))
  {
    return(NULL);
  }
  poly m;
  int dummy,i,j;
  m = p_ISet(1,currRing);
  for(i=1;i<=currRing->N;i++)
  {
    dummy=0;
    for(j=IDELEMS(I)-1;j>=0;j--)
    {
      if(p_GetExp(I->m[j],i,currRing) > dummy)
      {
        dummy = p_GetExp(I->m[j],i,currRing);
      }
    }
    p_SetExp(m,i,dummy,currRing);
  }
  p_Setm(m,currRing);
  return(m);
}

minimalMonomialGenSet()

static ideal minimalMonomialGenSet ( ideal  I )

static

Definition at line 1779 of file hilb.cc.

{
  /*
   * eliminates monomials which
   * can be generated by others in I
   */
  //first sort monomials of the ideal
 
  idSkipZeroes(I);
 
  sortMonoIdeal_pCompare(I);
 
  int i, k;
  int ICount = IDELEMS(I);
 
  for(k = ICount - 1; k >=1; k--)
  {
    for(i = 0; i < k; i++)
    {
 
      if(p_LmDivisibleBy(I->m[i], I->m[k], currRing))
      {
        pDelete(&(I->m[k]));
        break;
      }
    }
  }
 
  idSkipZeroes(I);
  return(I);
}

monCompare()

static int monCompare ( const void *  m, const void *  n  )

static

Definition at line 1759 of file hilb.cc.

{
  /* compares monomials */
 
 return(p_Compare(*(poly*) m, *(poly*)n, currRing));
}

positionInOrbit_FG_Case()

static int positionInOrbit_FG_Case ( ideal  I, poly  , std::vector< ideal >  idorb, std::vector< poly >  , int  , int    )

static

Definition at line 1654 of file hilb.cc.

{
  /*
   * It compares the ideal I with ideals in the set 'idorb'.
   * I and ideals of 'idorb' are sorted.
   *
   * It returns 0 if I is not equal to any ideal of 'idorb'
   * else returns position of the matched ideal.
   */
  int ps = 0;
  int i, s = 0;
  int OrbCount = idorb.size();
 
  if(idIs0(I))
  {
    return(1);
  }
 
  for(i = 1; i < OrbCount; i++)
  {
    s = comapreMonoIdBases(I, idorb[i]);
    if(s)
    {
      ps = i + 1;
      break;
    }
  }
 
  return(ps);
}

positionInOrbit_IG_Case()

static int positionInOrbit_IG_Case ( ideal  I, poly  w, std::vector< ideal >  idorb, std::vector< poly >  polist, int  trInd, int  trunDegHs  )

static

Definition at line 1576 of file hilb.cc.

{
  /*
   * It compares the ideal I with ideals in the set 'idorb'
   * up to total degree =
   * trInd - max(deg of w, deg of word in polist) polynomials.
   *
   * It returns 0 if I is not equal to any ideal in the
   * 'idorb' else returns position of the matched ideal.
   */
 
  int ps = 0;
  int i, s = 0;
  int orbCount = idorb.size();
 
  if(idIs0(I))
  {
    return(1);
  }
 
  int degw = p_Totaldegree(w, currRing);
  int degp;
  int dtr;
  int dtrp;
 
  dtr = trInd - degw;
  int IwCount;
 
   IwCount = CountOnIdUptoTruncationIndex(I, dtr);
 
  if(IwCount == 0)
  {
    return(1);
  }
 
  int ObCount;
 
  bool flag2 = FALSE;
 
  for(i = 1;i < orbCount; i++)
  {
    degp = p_Totaldegree(polist[i], currRing);
    if(degw > degp)
    {
      dtr = trInd - degw;
 
      ObCount = 0;
      ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
      if(ObCount == 0)
      {continue;}
      if(flag2)
      {
        IwCount = 0;
        IwCount = CountOnIdUptoTruncationIndex(I, dtr);
        flag2 = FALSE;
      }
    }
    else
    {
      flag2 = TRUE;
      dtrp = trInd - degp;
      ObCount = 0;
      ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtrp);
      IwCount = 0;
      IwCount = CountOnIdUptoTruncationIndex(I, dtrp);
    }
 
    s = comapreMonoIdBases_IG_Case(I, IwCount, idorb[i], ObCount);
 
    if(s)
    {
      ps = i + 1;
      break;
    }
  }
  return(ps);
}

positionInOrbitTruncationCase()

static int positionInOrbitTruncationCase ( ideal  I, poly  w, std::vector< ideal >  idorb, std::vector< poly >  polist, int  , int  trunDegHs  )

static

Definition at line 1685 of file hilb.cc.

{
  /*
   * It compares the ideal I with ideals in the set 'idorb'.
   * I and ideals in 'idorb' are sorted.
 
   * returns 0 if I is not equal to any ideal of 'idorb'
   * else returns position of the matched ideal.
   */
 
  int ps = 0;
  int i, s = 0;
  int OrbCount = idorb.size();
  int dtr=0; int IwCount, ObCount;
  dtr = trunDegHs - 1 - p_Totaldegree(w, currRing);
 
  if(idIs0(I))
  {
    for(i = 1; i < OrbCount; i++)
    {
      if(p_Totaldegree(w, currRing) == p_Totaldegree(polist[i], currRing))
      {
        if(idIs0(idorb[i]))
          return(i+1);
        ObCount=0;
        ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
        if(ObCount==0)
        {
          ps = i + 1;
          break;
        }
      }
    }
 
    return(ps);
  }
 
  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
 
  if(p_Totaldegree(I->m[0], currRing)==0)
  {
    for(i = 1; i < OrbCount; i++)
    {
       if(idIs0(idorb[i]))
         continue;
       if(p_Totaldegree(idorb[i]->m[0], currRing)==0)
       {
         ps = i + 1;
         break;
       }
     }
     return(ps);
  }
 
  for(i = 1; i < OrbCount; i++)
  {
    if(p_Totaldegree(w, currRing) == p_Totaldegree(polist[i], currRing))
    {
      if(idIs0(idorb[i]))
        continue;
      ObCount=0;
      ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
      s = comapreMonoIdBases_IG_Case(I, IwCount, idorb[i], ObCount);
      if(s)
      {
        ps = i + 1;
        break;
      }
    }
  }
 
  return(ps);
}

RightColonOperation()

ideal RightColonOperation	(	ideal	S,
		poly	w,
		int	lV
	)

Definition at line 2321 of file hilb.cc.

{
  /*
   * This returns right colon ideal of a monomial two-sided ideal of
   * the free associative algebra with respect to a monomial 'w'
   * (S:_R w).
   */
  S = minimalMonomialGenSet(S);
  ideal Iw = idInit(1,1);
  Iw = colonIdeal(S, w, lV, Iw, 0);
  return (Iw);
}

rouneslice()

void rouneslice	(	ideal	I,
		ideal	S,
		poly	q,
		poly	x,
		int &	prune,
		int &	moreprune,
		int &	steps,
		int &	NNN,
		mpz_ptr &	hilbertcoef,
		int *&	hilbpower
	)

Definition at line 974 of file hilb.cc.

{
  loop
  {
    (steps)++;
    int i,j;
    int dummy;
    poly m;
    ideal p;
    //----------- PRUNING OF S ---------------
    //S SHOULD IN THIS POINT BE ORDERED BY DEGREE
    for(i=IDELEMS(S)-1;i>=0;i--)
    {
      if(IsIn(S->m[i],I))
      {
        p_Delete(&S->m[i],currRing);
        prune++;
      }
    }
    idSkipZeroes(S);
    //----------------------------------------
    for(i=IDELEMS(I)-1;i>=0;i--)
    {
      m = p_Head(I->m[i],currRing);
      for(j=1;j<=currRing->N;j++)
      {
        dummy = p_GetExp(m,j,currRing);
        if(dummy > 0)
        {
          p_SetExp(m,j,dummy-1,currRing);
        }
      }
      p_Setm(m, currRing);
      if(IsIn(m,S))
      {
        p_Delete(&I->m[i],currRing);
        //printf("\n Deleted, since pi(m) is in S\n");pWrite(m);
      }
      p_Delete(&m,currRing);
    }
    idSkipZeroes(I);
    //----------- MORE PRUNING OF S ------------
    m = LCMmon(I);
    if(m != NULL)
    {
      for(i=0;i<IDELEMS(S);i++)
      {
        if(!(p_DivisibleBy(S->m[i], m, currRing)))
        {
          S->m[i] = NULL;
          j++;
          moreprune++;
        }
        else
        {
          if(pLmEqual(S->m[i],m))
          {
            S->m[i] = NULL;
            moreprune++;
          }
        }
      }
      idSkipZeroes(S);
    }
    p_Delete(&m,currRing);
    /*printf("\n---------------------------\n");
    printf("\n      I\n");idPrint(I);
    printf("\n      S\n");idPrint(S);
    printf("\n      q\n");pWrite(q);
    getchar();*/
 
    if(idIs0(I))
    {
      id_Delete(&I, currRing);
      id_Delete(&S, currRing);
      break;
    }
    m = LCMmon(I);
    if(!p_DivisibleBy(x,m, currRing))
    {
      //printf("\nx does not divide lcm(I)");
      //printf("\nEmpty set");pWrite(q);
      id_Delete(&I, currRing);
      id_Delete(&S, currRing);
      p_Delete(&m, currRing);
      break;
    }
    p_Delete(&m, currRing);
    m = SqFree(I);
    if(m==NULL)
    {
      //printf("\n      Corner: ");
      //pWrite(q);
      //printf("\n      With the facets of the dual simplex:\n");
      //idPrint(I);
      mpz_t ec;
      mpz_init(ec);
      mpz_ptr ec_ptr = ec;
      eulerchar(I, currRing->N, ec_ptr);
      bool flag = FALSE;
      if(NNN==0)
      {
        hilbertcoef = (mpz_ptr)omAlloc((NNN+1)*sizeof(mpz_t));
        hilbpower = (int*)omAlloc((NNN+1)*sizeof(int));
        mpz_init_set( &hilbertcoef[NNN], ec);
        hilbpower[NNN] = p_Totaldegree(q,currRing);
        NNN++;
      }
      else
      {
        //I look if the power appears already
        for(i = 0;(i<NNN)&&(flag == FALSE)&&(p_Totaldegree(q,currRing)>=hilbpower[i]);i++)
        {
          if((hilbpower[i]) == (p_Totaldegree(q,currRing)))
          {
            flag = TRUE;
            mpz_add(&hilbertcoef[i],&hilbertcoef[i],ec_ptr);
          }
        }
        if(flag == FALSE)
        {
          hilbertcoef = (mpz_ptr)omRealloc(hilbertcoef, (NNN+1)*sizeof(mpz_t));
          hilbpower = (int*)omRealloc(hilbpower, (NNN+1)*sizeof(int));
          mpz_init(&hilbertcoef[NNN]);
          for(j = NNN; j>i; j--)
          {
            mpz_set(&hilbertcoef[j],&hilbertcoef[j-1]);
            hilbpower[j] = hilbpower[j-1];
          }
          mpz_set(  &hilbertcoef[i], ec);
          hilbpower[i] = p_Totaldegree(q,currRing);
          NNN++;
        }
      }
      mpz_clear(ec);
      id_Delete(&I, currRing);
      id_Delete(&S, currRing);
      break;
    }
    else
      p_Delete(&m, currRing);
    m = ChooseP(I);
    p = idInit(1,1);
    p->m[0] = m;
    ideal Ip = idQuotMon(I,p);
    ideal Sp = idQuotMon(S,p);
    poly pq = pp_Mult_mm(q,m,currRing);
    rouneslice(Ip, Sp, pq, x, prune, moreprune, steps, NNN, hilbertcoef,hilbpower);
    idAddMon(S,p);
    p->m[0]=NULL;
    id_Delete(&p, currRing); // p->m[0] was also in S
    p_Delete(&pq,currRing);
  }
}

SearchP()

static poly SearchP ( ideal  I )

static

searches for a monomial of degree d>=2 and divides it by a variable (result monomial of deg d-1)

Definition at line 780 of file hilb.cc.

{
    int i,j,exp;
    poly res;
    if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)<=1)
    {
        res = ChoosePVar(I);
        return(res);
    }
    i = IDELEMS(I)-1;
    res = p_Copy(I->m[i], currRing);
    for(j=1;j<=currRing->N;j++)
    {
        exp = p_GetExp(I->m[i], j, currRing);
        if(exp > 0)
        {
            p_SetExp(res, j, exp - 1, currRing);
            p_Setm(res,currRing);
            break;
        }
    }
    assume( j <= currRing->N );
    return(res);
}

shiftInMon()

static poly shiftInMon ( poly  p, int  i, int  lV, const ring  r  )

static

Definition at line 1811 of file hilb.cc.

{
  /*
   * shifts the varibles of monomial p in the  i^th layer,
   * p remains unchanged,
   * creates new poly and returns it for the colon ideal
   */
  poly smon = p_One(r);
  int j, sh, cnt;
  cnt = r->N;
  sh = i*lV;
  int *e=(int *)omAlloc((r->N+1)*sizeof(int));
  int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
  p_GetExpV(p, e, r);
 
  for(j = 1; j <= cnt; j++)
  {
    if(e[j] == 1)
    {
      s[j+sh] = e[j];
    }
  }
 
  p_SetExpV(smon, s, currRing);
  omFree(e);
  omFree(s);
 
  p_SetComp(smon, p_GetComp(p, currRing), currRing);
  p_Setm(smon, currRing);
 
  return(smon);
}

slicehilb()

void slicehilb

(

ideal

I

)

Definition at line 1130 of file hilb.cc.

{
    //printf("Adi changes are here: \n");
    int i, NNN = 0;
    int steps = 0, prune = 0, moreprune = 0;
    mpz_ptr hilbertcoef;
    int *hilbpower;
    ideal S = idInit(1,1);
    poly q = p_One(currRing);
    ideal X = idInit(1,1);
    X->m[0]=p_One(currRing);
    for(i=1;i<=currRing->N;i++)
    {
      p_SetExp(X->m[0],i,1,currRing);
    }
    p_Setm(X->m[0],currRing);
    I = id_Mult(I,X,currRing);
    ideal Itmp = SortByDeg(I);
    id_Delete(&I,currRing);
    I = Itmp;
    //printf("\n-------------RouneSlice--------------\n");
    rouneslice(I,S,q,X->m[0],prune, moreprune, steps, NNN, hilbertcoef, hilbpower);
    id_Delete(&X,currRing);
    p_Delete(&q,currRing);
    //printf("\nIn total Prune got rid of %i elements\n",prune);
    //printf("\nIn total More Prune got rid of %i elements\n",moreprune);
    //printf("\nSteps of rouneslice: %i\n\n", steps);
    printf("\n//  %8d t^0",1);
    for(i = 0; i<NNN; i++)
    {
      if(mpz_sgn(&hilbertcoef[i])!=0)
      {
        gmp_printf("\n//  %8Zd t^%d",&hilbertcoef[i],hilbpower[i]);
      }
    }
    PrintLn();
    omFreeSize(hilbertcoef, (NNN)*sizeof(mpz_t));
    omFreeSize(hilbpower, (NNN)*sizeof(int));
    //printf("\n-------------------------------------\n");
}

SortByDeg()

static ideal SortByDeg ( ideal  I )

static

Definition at line 388 of file hilb.cc.

{
  if(idIs0(I))
  {
    return id_Copy(I,currRing);
  }
  int i;
  ideal res;
  idSkipZeroes(I);
  res = idInit(1,1);
  for(i = 0; i<=IDELEMS(I)-1;i++)
  {
    SortByDeg_p(res, I->m[i]);
    I->m[i]=NULL; // I->m[i] is now in res
  }
  idSkipZeroes(res);
  //idDegSortTest(res);
  return(res);
}

SortByDeg_p()

static void SortByDeg_p ( ideal  I, poly  p  )

static

!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!

Definition at line 288 of file hilb.cc.

{
  int i,j;
  if(idIs0(I))
  {
    I->m[0] = p;
    return;
  }
  idSkipZeroes(I);
  #if 1
  for(i = 0; (i<IDELEMS(I)) && (p_Totaldegree(I->m[i],currRing)<=p_Totaldegree(p,currRing)); i++)
  {
    if(p_DivisibleBy( I->m[i],p, currRing))
    {
      p_Delete(&p,currRing);
      return;
    }
  }
  for(i = IDELEMS(I)-1; (i>=0) && (p_Totaldegree(I->m[i],currRing)>=p_Totaldegree(p,currRing)); i--)
  {
    if(p_DivisibleBy(p,I->m[i], currRing))
    {
      p_Delete(&I->m[i],currRing);
    }
  }
  if(idIs0(I))
  {
    idSkipZeroes(I);
    I->m[0] = p;
    return;
  }
  #endif
  idSkipZeroes(I);
  //First I take the case when all generators have the same degree
  if(p_Totaldegree(I->m[0],currRing) == p_Totaldegree(I->m[IDELEMS(I)-1],currRing))
  {
    if(p_Totaldegree(p,currRing)<p_Totaldegree(I->m[0],currRing))
    {
      idInsertPoly(I,p);
      idSkipZeroes(I);
      for(i=IDELEMS(I)-1;i>=1; i--)
      {
        I->m[i] = I->m[i-1];
      }
      I->m[0] = p;
      return;
    }
    if(p_Totaldegree(p,currRing)>=p_Totaldegree(I->m[IDELEMS(I)-1],currRing))
    {
      idInsertPoly(I,p);
      idSkipZeroes(I);
      return;
    }
  }
  if(p_Totaldegree(p,currRing)<=p_Totaldegree(I->m[0],currRing))
  {
    idInsertPoly(I,p);
    idSkipZeroes(I);
    for(i=IDELEMS(I)-1;i>=1; i--)
    {
      I->m[i] = I->m[i-1];
    }
    I->m[0] = p;
    return;
  }
  if(p_Totaldegree(p,currRing)>=p_Totaldegree(I->m[IDELEMS(I)-1],currRing))
  {
    idInsertPoly(I,p);
    idSkipZeroes(I);
    return;
  }
  for(i = IDELEMS(I)-2; ;)
  {
    if(p_Totaldegree(p,currRing)==p_Totaldegree(I->m[i],currRing))
    {
      idInsertPoly(I,p);
      idSkipZeroes(I);
      for(j = IDELEMS(I)-1; j>=i+1;j--)
      {
        I->m[j] = I->m[j-1];
      }
      I->m[i] = p;
      return;
    }
    if(p_Totaldegree(p,currRing)>p_Totaldegree(I->m[i],currRing))
    {
      idInsertPoly(I,p);
      idSkipZeroes(I);
      for(j = IDELEMS(I)-1; j>=i+2;j--)
      {
        I->m[j] = I->m[j-1];
      }
      I->m[i+1] = p;
      return;
    }
    i--;
  }
}

sortMonoIdeal_pCompare()

void sortMonoIdeal_pCompare

(

ideal

I

)

Definition at line 1766 of file hilb.cc.

{
  /*
   * sorts monomial ideal in ascending order
   * order must be a total degree
   */
 
  qsort(I->m, IDELEMS(I), sizeof(poly), monCompare);
 
}

SqFree()

static poly SqFree ( ideal  I )

static

Definition at line 880 of file hilb.cc.

{
    int i,j;
    bool flag=TRUE;
    poly notsqrfree = NULL;
    if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)<=1)
    {
        return(notsqrfree);
    }
    for(i=IDELEMS(I)-1;(i>=0)&&(flag);i--)
    {
        for(j=1;(j<=currRing->N)&&(flag);j++)
        {
            if(p_GetExp(I->m[i],j,currRing)>1)
            {
                flag=FALSE;
                notsqrfree = p_ISet(1,currRing);
                p_SetExp(notsqrfree,j,1,currRing);
            }
        }
    }
    if(notsqrfree != NULL)
    {
        p_Setm(notsqrfree,currRing);
    }
    return(notsqrfree);
}

TwordMap()

static void TwordMap ( poly  p, poly  w, int  lV, int  d, ideal  Jwi, bool &  flag  )

static

Definition at line 1878 of file hilb.cc.

{
  /*
   * computes T_w(p) in a new poly object and places it
   * in Jwi which stores elements of colon ideal of I,
   * p and w remain unchanged,
   * the new polys for Jwi are constructed by sub-routines
   * deleteInMon, shiftInMon, p_MDivide,
   * places the result in Jwi and deletes the new polys
   * coming in dw, smon, qmon
   */
  int i;
  poly smon, dw;
  poly qmonp = NULL;
  bool del;
 
  for(i = 0;i <= d - 1; i++)
  {
    dw = deleteInMon(w, i, lV, currRing);
    smon = shiftInMon(p, i, lV, currRing);
    del = TRUE;
 
    if(pLmDivisibleBy(smon, w))
    {
      flag = TRUE;
      del  = FALSE;
 
      pDelete(&dw);
      pDelete(&smon);
 
      //delete all monomials of Jwi
      //and make Jwi =1
 
      for(int j = 0;j < IDELEMS(Jwi); j++)
      {
        pDelete(&Jwi->m[j]);
      }
 
      idInsertMonomial(Jwi, p_One(currRing));
      break;
    }
 
    if(pLmDivisibleBy(dw, smon))
    {
      del = FALSE;
      qmonp = p_MDivide(smon, dw, currRing);
      idInsertMonomial(Jwi, shiftInMon(qmonp, -d, lV, currRing));
      pLmFree(&qmonp);
      pDelete(&dw);
      pDelete(&smon);
    }
    //in case both if are false, delete dw and smon
    if(del)
    {
      pDelete(&dw);
      pDelete(&smon);
    }
  }
 
}

Variable Documentation

hLength

STATIC_VAR int hLength

Definition at line 46 of file hilb.cc.

Q0

STATIC_VAR int* Q0

Definition at line 45 of file hilb.cc.

Ql

STATIC_VAR int * Ql

Definition at line 45 of file hilb.cc.

Qpol

STATIC_VAR int** Qpol

Definition at line 44 of file hilb.cc.