ppinitialReduction.cc File Reference

#include "polys/monomials/p_polys.h"
#include "Singular/ipid.h"
#include "singularWishlist.h"
#include "ppinitialReduction.h"

Go to the source code of this file.

Functions
bool	isOrderingLocalInT (const ring r)
void	divideByCommonGcd (poly &g, const ring r)
void	pReduce (poly &g, const number p, const ring r)
bool	p_xLeadmonomDivisibleBy (const poly g, const poly f, const ring r)
void	pReduceInhomogeneous (poly &g, const number p, const ring r)
void	ptNormalize (poly *gStar, const number p, const ring r)
void	ptNormalize (ideal I, const number p, const ring r)
BOOLEAN	ptNormalize (leftv res, leftv args)
void	pReduce (ideal &I, const number p, const ring r)
bool	ppreduceInitially (poly *hStar, const poly g, const ring r)
	reduces h initially with respect to g, returns false if h was initially reduced in the first place, returns true if reductions have taken place. More...
bool	ppreduceInitially (ideal I, const number p, const ring r)
int	ppreduceInitially (ideal I, const number p, const poly g, const ring r)
static poly	ppNext (poly p, int l)
static void	sortMarks (const ideal H, const ring r, std::vector< mark > &T)
static poly	getTerm (const ideal H, const mark ab)
static void	adjustMarks (std::vector< mark > &T, const int newEntry)
static void	cleanupMarks (const ideal H, std::vector< mark > &T)
bool	ppreduceInitially (ideal &H, const number p, const ideal G, const ring r)
bool	ppreduceInitially (ideal I, const ring r, const number p)
	reduces I initially with respect to itself. More...

Function Documentation

adjustMarks()

static void adjustMarks ( std::vector< mark > &  T, const int  newEntry  )

static

Definition at line 478 of file ppinitialReduction.cc.

{
  for (unsigned i=0; i<T.size(); i++)
  {
    if (T[i].first>=newEntry)
      T[i].first = T[i].first+1;
  }
  return;
}

cleanupMarks()

static void cleanupMarks ( const ideal  H, std::vector< mark > &  T  )

static

Definition at line 489 of file ppinitialReduction.cc.

{
  for (unsigned i=0; i<T.size();)
  {
    if (getTerm(H,T[i])==NULL)
      T.erase(T.begin()+i);
    else
      i++;
  }
  return;
}

divideByCommonGcd()

void divideByCommonGcd	(	poly &	g,
		const ring	r
	)

Definition at line 21 of file ppinitialReduction.cc.

{
  number commonGcd = n_Copy(p_GetCoeff(g,r),r->cf);
  for (poly gCache=pNext(g); gCache; pIter(gCache))
  {
    number commonGcdCache = n_Gcd(commonGcd,p_GetCoeff(gCache,r),r->cf);
    n_Delete(&commonGcd,r->cf);
    commonGcd = commonGcdCache;
    if (n_IsOne(commonGcd,r->cf))
    {
      n_Delete(&commonGcd,r->cf);
      return;
    }
  }
  for (poly gCache=g; gCache; pIter(gCache))
  {
    number oldCoeff = p_GetCoeff(gCache,r);
    number newCoeff = n_Div(oldCoeff,commonGcd,r->cf);
    p_SetCoeff(gCache,newCoeff,r);
  }
  p_Test(g,r);
  n_Delete(&commonGcd,r->cf);
  return;
}

getTerm()

static poly getTerm ( const ideal  H, const mark  ab  )

static

Definition at line 470 of file ppinitialReduction.cc.

{
  int a = ab.first;
  int b = ab.second;
  return ppNext(H->m[a],b);
}

isOrderingLocalInT()

bool isOrderingLocalInT

(

const ring

r

)

Definition at line 8 of file ppinitialReduction.cc.

{
  poly one = p_One(r);
  poly t = p_One(r);
  p_SetExp(t,1,1,r);
  p_Setm(t,r);
  int s = p_LmCmp(one,t,r);
  p_Delete(&one,r);
  p_Delete(&t,r);
  return (s==1);
}

p_xLeadmonomDivisibleBy()

bool p_xLeadmonomDivisibleBy	(	const poly	g,
		const poly	f,
		const ring	r
	)

Definition at line 118 of file ppinitialReduction.cc.

{
  poly gx = p_Head(g,r);
  poly fx = p_Head(f,r);
  p_SetExp(gx,1,0,r);
  p_SetExp(fx,1,0,r);
  p_Setm(gx,r);
  p_Setm(fx,r);
  bool b = p_LeadmonomDivisibleBy(gx,fx,r);
  p_Delete(&gx,r);
  p_Delete(&fx,r);
  return b;
}

ppNext()

static poly ppNext ( poly  p, int  l  )

static

Definition at line 430 of file ppinitialReduction.cc.

{
  poly q = p;
  for (int i=0; i<l; i++)
  {
    if (q==NULL)
      break;
    pIter(q);
  }
  return q;
}

ppreduceInitially() [1/5]

bool ppreduceInitially	(	ideal &	H,
		const number	p,
		const ideal	G,
		const ring	r
	)

Definition at line 508 of file ppinitialReduction.cc.

{
  /***
   * Step 1: reduce H initially with respect to itself and with respect to p-t
   **/
  if (ppreduceInitially(H,p,r)) return true;
 
  /***
   * Step 2: initialize an ideal I in which the reductions will take place-
   *   along the reduction it will be enlarged with elements that will be discarded at the end
   *   initialize a working list T which keeps track.
   *   the working list T is a vector of pairs of integer.
   *   if T contains a pair (i,j) then that means that in the i-th element of H
   *   term j and subsequent terms need to be checked for reduction.
   *   T is sorted by the ordering on the temrs the pairs correspond to.
   **/
  int m=IDELEMS(H);
  ideal I = idInit(m);
  std::vector<mark> T;
  for (int i=0; i<m; i++)
  {
    if(H->m[i]!=NULL)
    {
      I->m[i]=H->m[i];
      if (pNext(I->m[i])!=NULL)
        T.push_back(std::pair<int,int>(i,1));
    }
  }
 
  /***
   * Step 3: as long as the working list is not empty, successively reduce terms in it
   *   by adding suitable elements to I and reducing it initially with respect to itself
   **/
  int k=IDELEMS(G);
  while (T.size()>0)
  {
    sortMarks(I,r,T);
    int i=0;
    for (; i<k; i++)
    {
      if(G->m[i]!=NULL)
      {
        if (p_LeadmonomDivisibleBy(G->m[i],getTerm(I,T[0]),r)) break;
      }
    }
    if (i<k)
    {
      poly g = p_One(r); poly h0 = getTerm(I,T[0]);
      assume(h0!=NULL);
      for (int j=2; j<=r->N; j++)
        p_SetExp(g,j,p_GetExp(h0,j,r)-p_GetExp(G->m[i],j,r),r);
      p_Setm(g,r);
      g = p_Mult_q(g,p_Copy(G->m[i],r),r);
      int newEntry = ppreduceInitially(I,p,g,r);
      adjustMarks(T,newEntry);
    }
    else
      T[0].second = T[0].second+1;
    cleanupMarks(I,T);
  }
 
  /***
   * Step 4: cleanup, delete all polynomials in I which have been added in Step 3
   **/
  k=IDELEMS(I);
  for (int i=0; i<k; i++)
  {
    if(I->m[i]!=NULL)
    {
      for (int j=0; j<m; j++)
      {
        if ((H->m[j]!=NULL)
        && (p_LeadmonomDivisibleBy(H->m[j],I->m[i],r)))
        {
          I->m[i]=NULL;
          break;
        }
      }
    }
  }
  id_Delete(&I,r);
  return false;
}

ppreduceInitially() [2/5]

int ppreduceInitially	(	ideal	I,
		const number	p,
		const poly	g,
		const ring	r
	)

Definition at line 372 of file ppinitialReduction.cc.

{
  id_Test(I,r);
  p_Test(g,r);
  idInsertPoly(I,g);
  idSkipZeroes(I);
  int n=IDELEMS(I);
  int j;
  for (j=n-1; j>0; j--)
  {
    if (p_LmCmp(I->m[j], I->m[j-1],r)>0) /*p_LmCmp(p,q) requires: p,q!=NULL */
    {
      poly cache = I->m[j];
      I->m[j] = I->m[j-1];
      I->m[j-1] = cache;
    }
    else
      break;
  }
 
  /***
   * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
   * removing terms with the same monomials in x as lt(g_j) out of g_k for j<k
   **/
  for (int i=0; i<j; i++)
    if (ppreduceInitially(&I->m[j], I->m[i], r))
      pReduce(I->m[j],p,r);
  for (int k=j+1; k<n; k++)
    if (ppreduceInitially(&I->m[k], I->m[j], r))
    {
      pReduce(I->m[k],p,r);
      for (int l=j+1; l<k; l++)
        if (ppreduceInitially(&I->m[k], I->m[l], r))
          pReduce(I->m[k],p,r);
    }
 
  /***
   * the second pass. removing terms divisible by lt(g_j) and lt(g_k) out of g_i for i<j<k
   * removing terms divisible by lt(g_k) out of g_j for j<k
   **/
  for (int i=0; i<j; i++)
    for (int k=j; k<n; k++)
      if (ppreduceInitially(&I->m[i], I->m[k], r))
        pReduce(I->m[i],p,r);
  for (int k=j; k<n-1; k++)
    for (int l=k+1; l<n; l++)
      if (ppreduceInitially(&I->m[k], I->m[l], r))
        pReduce(I->m[k],p,r);
 
  /***
   * removes the elements of I which have been reduced to 0 in the previous two passes
   **/
  idSkipZeroes(I);
  id_Test(I,r);
  return j;
}

ppreduceInitially() [3/5]

bool ppreduceInitially	(	ideal	I,
		const number	p,
		const ring	r
	)

Definition at line 321 of file ppinitialReduction.cc.

{
  idSkipZeroes(I);
  int m=IDELEMS(I),n=m; poly cache;
  do
  {
    int j=0;
    for (int i=1; i<n; i++)
    {
      if (p_LmCmp(I->m[i-1],I->m[i],r)<0) /*p_LmCmp(p,q): requires: p,q!=NULL*/
      {
        cache=I->m[i-1];
        I->m[i-1]=I->m[i];
        I->m[i]=cache;
        j = i;
      }
    }
    n=j;
  } while(n);
  for (int i=0; i<m; i++)
    pReduce(I->m[i],p,r);
 
  /***
   * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
   **/
  for (int i=0; i<m-1; i++)
    for (int j=i+1; j<m; j++)
      if (ppreduceInitially(&I->m[j], I->m[i], r))
        pReduce(I->m[j],p,r);
 
  /***
   * the second pass. removing terms divisible by lt(g_j) out of g_i for i<j
   **/
  for (int i=0; i<m-1; i++)
    for (int j=i+1; j<m; j++)
      if (ppreduceInitially(&I->m[i], I->m[j],r))
        pReduce(I->m[i],p,r);
 
  /***
   * removes the elements of I which have been reduced to 0 in the previous two passes
   **/
  idSkipZeroes(I);
  return false;
}

ppreduceInitially() [4/5]

bool ppreduceInitially	(	ideal	I,
		const ring	r,
		const number	p
	)

reduces I initially with respect to itself.

assumes that the generators of I are homogeneous in x and that p-t is in I.

sorts Hi according to degree in t in descending order (lowest first, highest last)

Definition at line 597 of file ppinitialReduction.cc.

{
  assume(!n_IsUnit(p,r->cf));
 
  /***
   * Step 1: split up I into components of same degree in x
   *  the lowest component should only contain p-t
   **/
  std::map<long,ideal> H; int n = IDELEMS(I);
  for (int i=0; i<n; i++)
  {
    if(I->m[i]!=NULL)
    {
      I->m[i] = p_Cleardenom(I->m[i],r);
      long d = 0;
      for (int j=2; j<=r->N; j++)
        d += p_GetExp(I->m[i],j,r);
      std::map<long,ideal>::iterator it = H.find(d);
      if (it != H.end())
        idInsertPoly(it->second,I->m[i]);
      else
      {
        std::pair<long,ideal> Hd(d,idInit(1));
        Hd.second->m[0] = I->m[i];
        H.insert(Hd);
      }
    }
  }
 
  std::map<long,ideal>::iterator it=H.begin();
  ideal Hi = it->second;
  idShallowDelete(&Hi);
  it++;
  Hi = it->second;
 
  /***
   * Step 2: reduce each component initially with respect to itself
   *  and all lower components
   **/
  if (ppreduceInitially(Hi,p,r)) return true;
  idSkipZeroes(Hi);
  id_Test(Hi,r);
  id_Test(I,r);
 
  ideal G = idInit(n); int m=0;
  ideal GG = (ideal) omAllocBin(sip_sideal_bin);
  GG->nrows = 1; GG->rank = 1; GG->m=NULL;
 
  for (it++; it!=H.end(); it++)
  {
    int l=IDELEMS(Hi); int k=l; poly cache;
    /**
     * sorts Hi according to degree in t in descending order
     * (lowest first, highest last)
     */
    do
    {
      int j=0;
      for (int i=1; i<k; i++)
      {
        if (p_GetExp(Hi->m[i-1],1,r)<p_GetExp(Hi->m[i],1,r))
        {
          cache=Hi->m[i-1];
          Hi->m[i-1]=Hi->m[i];
          Hi->m[i]=cache;
          j = i;
        }
      }
      k=j;
    } while(k);
    int kG=n-m, kH=0;
    for (int i=n-m-l; i<n; i++)
    {
      if (kG==n)
      {
        memcpy(&(G->m[i]),&(Hi->m[kH]),(n-i)*sizeof(poly));
        break;
      }
      if (kH==l)
        break;
      if (p_GetExp(G->m[kG],1,r)>p_GetExp(Hi->m[kH],1,r))
        G->m[i] = G->m[kG++];
      else
        G->m[i] = Hi->m[kH++];
    }
    m += l; IDELEMS(GG) = m; GG->m = &G->m[n-m];
    id_Test(it->second,r);
    id_Test(GG,r);
    if (ppreduceInitially(it->second,p,GG,r)) return true;
    id_Test(it->second,r);
    id_Test(GG,r);
    idShallowDelete(&Hi); Hi = it->second;
  }
  idShallowDelete(&Hi);
 
  ptNormalize(I,p,r);
  omFreeBin((ADDRESS)GG, sip_sideal_bin);
  idShallowDelete(&G);
  return false;
}

ppreduceInitially() [5/5]

bool ppreduceInitially	(	poly *	hStar,
		const poly	g,
		const ring	r
	)

reduces h initially with respect to g, returns false if h was initially reduced in the first place, returns true if reductions have taken place.

assumes that h and g are in pReduced form and homogeneous in x of the same degree

Definition at line 282 of file ppinitialReduction.cc.

{
  poly h = *hStar;
  if (h==NULL || g==NULL)
    return false;
  p_Test(h,r);
  p_Test(g,r);
  poly hCache;
  for (hCache=h; hCache; pIter(hCache))
    if (p_LeadmonomDivisibleBy(g,hCache,r)) break;
  if (hCache)
  {
    number gAlpha = p_GetCoeff(g,r);
    poly hAlphaT = p_Init(r);
    p_SetCoeff(hAlphaT,n_Copy(p_GetCoeff(hCache,r),r->cf),r);
    p_SetExp(hAlphaT,1,p_GetExp(hCache,1,r)-p_GetExp(g,1,r),r);
    for (int i=2; i<=r->N; i++)
      p_SetExp(hAlphaT,i,0,r);
    p_Setm(hAlphaT,r); p_Test(hAlphaT,r);
    poly q1 = p_Mult_nn(h,gAlpha,r); p_Test(q1,r);
    poly q2 = p_Mult_q(p_Copy(g,r),hAlphaT,r); p_Test(q2,r);
    q2 = p_Neg(q2,r); p_Test(q2,r);
    h = p_Add_q(q1,q2,r);
    p_Test(h,r);
    p_Test(g,r);
    *hStar = h;
    return true;
  }
  p_Test(h,r);
  p_Test(g,r);
  return false;
}

pReduce() [1/2]

void pReduce	(	ideal &	I,
		const number	p,
		const ring	r
	)

Definition at line 261 of file ppinitialReduction.cc.

{
  int k = IDELEMS(I);
  for (int i=0; i<k; i++)
  {
    if (I->m[i]!=NULL)
    {
      number c = p_GetCoeff(I->m[i],r);
      if (!n_DivBy(p,c,r->cf))
        pReduce(I->m[i],p,r);
    }
  }
}

pReduce() [2/2]

void pReduce	(	poly &	g,
		const number	p,
		const ring	r
	)

Definition at line 54 of file ppinitialReduction.cc.

{
  if (g==NULL)
    return;
  p_Test(g,r);
 
  poly toBeChecked = pNext(g);
  pNext(g) = NULL; poly gEnd = g;
  poly gCache;
 
  number coeff, pPower; int power; poly subst;
  while(toBeChecked)
  {
    for (gCache = g; gCache; pIter(gCache))
      if (p_LeadmonomDivisibleBy(gCache,toBeChecked,r)) break;
    if (gCache)
    {
      n_Power(p,p_GetExp(toBeChecked,1,r)-p_GetExp(gCache,1,r),&pPower,r->cf);
      coeff = n_Mult(p_GetCoeff(toBeChecked,r),pPower,r->cf);
      p_SetCoeff(gCache,n_Add(p_GetCoeff(gCache,r),coeff,r->cf),r);
      n_Delete(&pPower,r->cf); n_Delete(&coeff,r->cf);
      toBeChecked=p_LmDeleteAndNext(toBeChecked,r);
    }
    else
    {
      if (n_DivBy(p_GetCoeff(toBeChecked,r),p,r->cf))
      {
        power=1;
        coeff=n_Div(p_GetCoeff(toBeChecked,r),p,r->cf);
        while (n_DivBy(coeff,p,r->cf))
        {
          power++;
          number coeff0 = n_Div(coeff,p,r->cf);
          n_Delete(&coeff,r->cf);
          coeff = coeff0;
          coeff0 = NULL;
          if (power<1)
          {
            WerrorS("pReduce: overflow in exponent");
            throw 0;
          }
        }
        subst=p_LmInit(toBeChecked,r);
        p_AddExp(subst,1,power,r);
        p_SetCoeff(subst,coeff,r);
        p_Setm(subst,r); p_Test(subst,r);
        toBeChecked=p_LmDeleteAndNext(toBeChecked,r);
        toBeChecked=p_Add_q(toBeChecked,subst,r);
        p_Test(toBeChecked,r);
      }
      else
      {
        pNext(gEnd)=toBeChecked;
        pIter(gEnd); pIter(toBeChecked);
        pNext(gEnd)=NULL;
        p_Test(g,r);
      }
    }
  }
  p_Test(g,r);
  divideByCommonGcd(g,r);
  return;
}

pReduceInhomogeneous()

void pReduceInhomogeneous	(	poly &	g,
		const number	p,
		const ring	r
	)

Definition at line 132 of file ppinitialReduction.cc.

{
  if (g==NULL)
    return;
  p_Test(g,r);
 
  poly toBeChecked = pNext(g);
  pNext(g) = NULL; poly gEnd = g;
  poly gCache;
 
  number coeff, pPower; int power; poly subst;
  while(toBeChecked)
  {
    for (gCache = g; gCache; pIter(gCache))
      if (p_xLeadmonomDivisibleBy(gCache,toBeChecked,r)) break;
    if (gCache)
    {
      n_Power(p,p_GetExp(toBeChecked,1,r)-p_GetExp(gCache,1,r),&pPower,r->cf);
      coeff = n_Mult(p_GetCoeff(toBeChecked,r),pPower,r->cf);
      p_SetCoeff(gCache,n_Add(p_GetCoeff(gCache,r),coeff,r->cf),r);
      n_Delete(&pPower,r->cf); n_Delete(&coeff,r->cf);
      toBeChecked=p_LmDeleteAndNext(toBeChecked,r);
    }
    else
    {
      if (n_DivBy(p_GetCoeff(toBeChecked,r),p,r->cf))
      {
        power=1;
        coeff=n_Div(p_GetCoeff(toBeChecked,r),p,r->cf);
        while (n_DivBy(coeff,p,r->cf))
        {
          power++;
          number coeff0 = n_Div(coeff,p,r->cf);
          n_Delete(&coeff,r->cf);
          coeff = coeff0;
          coeff0 = NULL;
          if (power<1)
          {
            WerrorS("pReduce: overflow in exponent");
            throw 0;
          }
        }
        subst=p_LmInit(toBeChecked,r);
        p_AddExp(subst,1,power,r);
        p_SetCoeff(subst,coeff,r);
        p_Setm(subst,r); p_Test(subst,r);
        toBeChecked=p_LmDeleteAndNext(toBeChecked,r);
        toBeChecked=p_Add_q(toBeChecked,subst,r);
        p_Test(toBeChecked,r);
      }
      else
      {
        pNext(gEnd)=toBeChecked;
        pIter(gEnd); pIter(toBeChecked);
        pNext(gEnd)=NULL;
        p_Test(g,r);
      }
    }
  }
  p_Test(g,r);
  divideByCommonGcd(g,r);
  return;
}

ptNormalize() [1/3]

void ptNormalize	(	ideal	I,
		const number	p,
		const ring	r
	)

Definition at line 230 of file ppinitialReduction.cc.

{
  for (int i=0; i<IDELEMS(I); i++)
    ptNormalize(&(I->m[i]),p,r);
  return;
}

ptNormalize() [2/3]

BOOLEAN ptNormalize	(	leftv	res,
		leftv	args
	)

Definition at line 238 of file ppinitialReduction.cc.

{
  leftv u = args;
  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
  {
    leftv v = u->next;
    if ((v!=NULL) && (v->Typ()==NUMBER_CMD))
    {
      omUpdateInfo();
      Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
      ideal I = (ideal) u->CopyD();
      number p = (number) v->CopyD();
      ptNormalize(I,p,currRing);
      n_Delete(&p,currRing->cf);
      res->rtyp = IDEAL_CMD;
      res->data = (char*) I;
      return FALSE;
    }
  }
  return TRUE;
}

ptNormalize() [3/3]

void ptNormalize	(	poly *	gStar,
		const number	p,
		const ring	r
	)

Definition at line 196 of file ppinitialReduction.cc.

{
  poly g = *gStar;
  if (g==NULL || n_DivBy(p_GetCoeff(g,r),p,r->cf))
    return;
  p_Test(g,r);
 
  // create p-t
  poly pt = p_Init(r);
  p_SetCoeff(pt,n_Copy(p,r->cf),r);
 
  pNext(pt) = p_Init(r);
  p_SetExp(pNext(pt),1,1,r);
  p_Setm(pNext(pt),r);
  p_SetCoeff(pNext(pt),n_Init(-1,r->cf),r);
 
  // make g monic with the help of p-t
  number a,b;
  number gcd = n_ExtGcd(p_GetCoeff(g,r),p,&a,&b,r->cf);
  assume(n_IsUnit(gcd,r->cf));
  // now a*leadcoef(g)+b*p = gcd with gcd being a unit
  // so a*g+b*(p-t)*leadmonom(g) should have a unit as leading coefficient
  poly m = p_Head(g,r);
  p_SetCoeff(m,n_Init(1,r->cf),r);
  g = p_Add_q(p_Mult_nn(g,a,r),p_Mult_nn(p_Mult_mm(pt,m,r),b,r),r);
  n_Delete(&a,r->cf);
  n_Delete(&b,r->cf);
  n_Delete(&gcd,r->cf);
  p_Delete(&m,r);
 
  p_Test(g,r);
  return;
}

sortMarks()

static void sortMarks ( const ideal  H, const ring  r, std::vector< mark > &  T  )

static

Definition at line 443 of file ppinitialReduction.cc.

{
  std::pair<int,int> pointerToTerm;
  int k=T.size();
  do
  {
    int j=0;
    for (int i=1; i<k-1; i++)
    {
      int generatorA = T[i-1].first;
      int termA = T[i-1].second;
      int generatorB = T[i].first;
      int termB = T[i].second;
      if (p_LmCmp(ppNext(H->m[generatorA],termA),ppNext(H->m[generatorB],termB),r)<0)
      {
        mark cache=T[i-1];
        T[i-1]=T[i];
        T[i]=cache;
        j = i;
      }
    }
    k=j;
  } while(k);
  return;
}