source: git/libpolys/polys/monomials/p_polys.h @ 6caad65

/****************************************
*  Computer Algebra System SINGULAR     *
****************************************/
/***************************************************************
 *  File:    p_polys.h
 *  Purpose: declaration of poly stuf which are independent of
 *           currRing
 *  Author:  obachman (Olaf Bachmann)
 *  Created: 9/00
 *******************************************************************/
/***************************************************************
 *  Purpose: implementation of poly procs which iter over ExpVector
 *  Author:  obachman (Olaf Bachmann)
 *  Created: 8/00
 *******************************************************************/
#ifndef P_POLYS_H
#define P_POLYS_H

#include <omalloc/omalloc.h>

#include <misc/mylimits.h>
#include <misc/intvec.h>
#include <coeffs/coeffs.h>

#include <polys/monomials/monomials.h>
#include <polys/monomials/ring.h>

#include <polys/templates/p_MemAdd.h>
#include <polys/templates/p_MemCmp.h>
#include <polys/templates/p_Procs.h>

#include <polys/sbuckets.h>

#ifdef HAVE_PLURAL
#include <polys/nc/nc.h>
#endif

/***************************************************************
 *
 * Divisiblity tests, args must be != NULL, except for
 * pDivisbleBy
 *
 ***************************************************************/
unsigned long p_GetShortExpVector(poly a, const ring r);

#ifdef HAVE_RINGS
/*! divisibility check over ground ring (which may contain zero divisors);
   TRUE iff LT(f) divides LT(g), i.e., LT(f)*c*m = LT(g), for some
   coefficient c and some monomial m;
   does not take components into account
 */
BOOLEAN p_DivisibleByRingCase(poly f, poly g, const ring r);
#endif

/***************************************************************
 *
 * Misc things on polys
 *
 ***************************************************************/

poly p_One(const ring r);

int p_MinDeg(poly p,intvec *w, const ring R);

long p_DegW(poly p, const short *w, const ring R);

/// return TRUE if all monoms have the same component
BOOLEAN   p_OneComp(poly p, const ring r);

/// return i, if head depends only on var(i)
int       p_IsPurePower(const poly p, const ring r);

/// return i, if poly depends only on var(i)
int       p_IsUnivariate(poly p, const ring r);

/// set entry e[i] to 1 if var(i) occurs in p, ignore var(j) if e[j]>0
/// return #(e[i]>0)
int      p_GetVariables(poly p, int * e, const ring r);

/// returns the poly representing the integer i
poly      p_ISet(long i, const ring r);

/// returns the poly representing the number n, destroys n
poly      p_NSet(number n, const ring r);

void  p_Vec2Polys(poly v, poly**p, int *len, const ring r);

/***************************************************************
 *
 * Copying/Deletion of polys: args may be NULL
 *
 ***************************************************************/

// simply deletes monomials, does not free coeffs
void p_ShallowDelete(poly *p, const ring r);



/***************************************************************
 *
 * Copying/Deleteion of polys: args may be NULL
 *  - p/q as arg mean a poly
 *  - m a monomial
 *  - n a number
 *  - pp (resp. qq, mm, nn) means arg is constant
 *  - p (resp, q, m, n)     means arg is destroyed
 *
 ***************************************************************/

poly      p_Sub(poly a, poly b, const ring r);

poly      p_Power(poly p, int i, const ring r);


/***************************************************************
 *
 * PDEBUG stuff
 *
 ***************************************************************/
#ifdef PDEBUG
// Returns TRUE if m is monom of p, FALSE otherwise
BOOLEAN pIsMonomOf(poly p, poly m);
// Returns TRUE if p and q have common monoms
BOOLEAN pHaveCommonMonoms(poly p, poly q);

// p_Check* routines return TRUE if everything is ok,
// else, they report error message and return false

// check if Lm(p) is from ring r
BOOLEAN p_LmCheckIsFromRing(poly p, ring r);
// check if Lm(p) != NULL, r != NULL and initialized && Lm(p) is from r
BOOLEAN p_LmCheckPolyRing(poly p, ring r);
// check if all monoms of p are from ring r
BOOLEAN p_CheckIsFromRing(poly p, ring r);
// check r != NULL and initialized && all monoms of p are from r
BOOLEAN p_CheckPolyRing(poly p, ring r);
// check if r != NULL and initialized
BOOLEAN p_CheckRing(ring r);
// only do check if cond


#define pIfThen(cond, check) do {if (cond) {check;}} while (0)

BOOLEAN _p_Test(poly p, ring r, int level);
BOOLEAN _p_LmTest(poly p, ring r, int level);
BOOLEAN _pp_Test(poly p, ring lmRing, ring tailRing, int level);

#define p_Test(p,r)     _p_Test(p, r, PDEBUG)
#define p_LmTest(p,r)   _p_LmTest(p, r, PDEBUG)
#define pp_Test(p, lmRing, tailRing)    _pp_Test(p, lmRing, tailRing, PDEBUG)

#else // ! PDEBUG

#define pIsMonomOf(p, q)        (TRUE)
#define pHaveCommonMonoms(p, q) (TRUE)
#define p_LmCheckIsFromRing(p,r)  ((void)0)
#define p_LmCheckPolyRing(p,r)    ((void)0)
#define p_CheckIsFromRing(p,r)  ((void)0)
#define p_CheckPolyRing(p,r)    ((void)0)
#define p_CheckRing(r)          ((void)0)
#define P_CheckIf(cond, check)  ((void)0)

#define p_Test(p,r)     ((void) 1)
#define p_LmTest(p,r)   ((void) 1)
#define pp_Test(p, lmRing, tailRing) ((void) 1)

#endif

/***************************************************************
 *
 * Misc stuff
 *
 ***************************************************************/
/*2
* returns the length of a polynomial (numbers of monomials)
*/
static inline int pLength(poly a)
{
  int l = 0;
  while (a!=NULL)
  {
    pIter(a);
    l++;
  }
  return l;
}

void      p_Norm(poly p1, const ring r);
void      p_Normalize(poly p,const ring r);

void      p_Content(poly p, const ring r);
#if 1
// currently only used by Singular/janet
void      p_SimpleContent(poly p, int s, const ring r);
#endif

poly      p_Cleardenom(poly p, const ring r);
void      p_Cleardenom_n(poly p, const ring r,number &c);
number    p_GetAllDenom(poly ph, const ring r);

int       p_Size( poly p, const ring r );

// homogenizes p by multiplying certain powers of the varnum-th variable
poly      p_Homogen (poly p, int varnum, const ring r);

BOOLEAN   p_IsHomogeneous (poly p, const ring r);

static inline void p_Setm(poly p, const ring r);
p_SetmProc p_GetSetmProc(ring r);

poly      p_Subst(poly p, int n, poly e, const ring r);

// TODO:
#define p_SetmComp  p_Setm

// component
static inline  unsigned long p_SetComp(poly p, unsigned long c, ring r)
{
  p_LmCheckPolyRing2(p, r);
  pAssume2(rRing_has_Comp(r));
  __p_GetComp(p,r) = c;
  return c;
}
// sets component of poly a to i
static inline   void p_SetCompP(poly p, int i, ring r)
{
  if (p != NULL)
  {
#ifdef PDEBUG
    p_Test(p, r);
#endif
    if (rOrd_SetCompRequiresSetm(r))
    {
      do
      {
        p_SetComp(p, i, r);
        p_SetmComp(p, r);
        pIter(p);
      }
      while (p != NULL);
    }
    else
    {
      do
      {
        p_SetComp(p, i, r);
        pIter(p);
      }
      while(p != NULL);
    }
  }
}

static inline   void p_SetCompP(poly p, int i, ring lmRing, ring tailRing)
{
  if (p != NULL)
  {
    p_SetComp(p, i, lmRing);
    p_SetmComp(p, lmRing);
    p_SetCompP(pNext(p), i, tailRing);
  }
}

// returns maximal column number in the modul element a (or 0)
static inline long p_MaxComp(poly p, ring lmRing, ring tailRing)
{
  long result,i;

  if(p==NULL) return 0;
  result = p_GetComp(p, lmRing);
  if (result != 0)
  {
    loop
    {
      pIter(p);
      if(p==NULL) break;
      i = p_GetComp(p, tailRing);
      if (i>result) result = i;
    }
  }
  return result;
}

static inline long p_MaxComp(poly p,ring lmRing) {return p_MaxComp(p,lmRing,lmRing);}

static inline   long p_MinComp(poly p, ring lmRing, ring tailRing)
{
  long result,i;

  if(p==NULL) return 0;
  result = p_GetComp(p,lmRing);
  if (result != 0)
  {
    loop
    {
      pIter(p);
      if(p==NULL) break;
      i = p_GetComp(p,tailRing);
      if (i<result) result = i;
    }
  }
  return result;
}

static inline long p_MinComp(poly p,ring lmRing) {return p_MinComp(p,lmRing,lmRing);}


static inline poly pReverse(poly p)
{
  if (p == NULL || pNext(p) == NULL) return p;

  poly q = pNext(p), // == pNext(p)
    qn;
  pNext(p) = NULL;
  do
  {
    qn = pNext(q);
    pNext(q) = p;
    p = q;
    q = qn;
  }
  while (qn != NULL);
  return p;
}
void      pEnlargeSet(poly**p, int length, int increment);


/***************************************************************
 *
 * I/O
 *
 ***************************************************************/
/// print p according to ShortOut in lmRing & tailRing
char*     p_String0(poly p, ring lmRing, ring tailRing);
char*     p_String(poly p, ring lmRing, ring tailRing);
void      p_Write(poly p, ring lmRing, ring tailRing);
void      p_Write0(poly p, ring lmRing, ring tailRing);
void      p_wrp(poly p, ring lmRing, ring tailRing);

/// print p in a short way, if possible
char* p_String0Short(const poly p, ring lmRing, ring tailRing);

/// print p in a long way
char* p_String0Long(const poly p, ring lmRing, ring tailRing);


/***************************************************************
 *
 * Degree stuff -- see p_polys.cc for explainations
 *
 ***************************************************************/

static inline long  p_FDeg(const poly p, const ring r)  { return r->pFDeg(p,r); }
static inline long  p_LDeg(const poly p, int *l, const ring r)  { return r->pLDeg(p,l,r); }

long p_WFirstTotalDegree(poly p, ring r);
long p_WTotaldegree(poly p, const ring r);
long p_WDegree(poly p,const ring r);
long pLDeg0(poly p,int *l, ring r);
long pLDeg0c(poly p,int *l, ring r);
long pLDegb(poly p,int *l, ring r);
long pLDeg1(poly p,int *l, ring r);
long pLDeg1c(poly p,int *l, ring r);
long pLDeg1_Deg(poly p,int *l, ring r);
long pLDeg1c_Deg(poly p,int *l, ring r);
long pLDeg1_Totaldegree(poly p,int *l, ring r);
long pLDeg1c_Totaldegree(poly p,int *l, ring r);
long pLDeg1_WFirstTotalDegree(poly p,int *l, ring r);
long pLDeg1c_WFirstTotalDegree(poly p,int *l, ring r);

BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r);

/// same as the usual p_EqualPolys for polys belonging to *equal* rings
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r1, const ring r2);

long p_Deg(poly a, const ring r);


/***************************************************************
 *
 * Primitives for accessing and setting fields of a poly
 *
 ***************************************************************/

static inline number p_SetCoeff(poly p, number n, ring r)
{
  p_LmCheckPolyRing2(p, r);
  n_Delete(&(p->coef), r->cf);
  (p)->coef=n;
  return n;
}

// order
static inline long p_GetOrder(poly p, ring r)
{
  p_LmCheckPolyRing2(p, r);
  if (r->typ==NULL) return ((p)->exp[r->pOrdIndex]);
  int i=0;
  loop
  {
    switch(r->typ[i].ord_typ)
    {
      case ro_am:
      case ro_wp_neg:
        return (((long)((p)->exp[r->pOrdIndex]))-POLY_NEGWEIGHT_OFFSET);
      case ro_syzcomp:
      case ro_syz:
      case ro_cp:
        i++;
        break;
      //case ro_dp:
      //case ro_wp:
      default:
        return ((p)->exp[r->pOrdIndex]);
    }
  }
}

// Setm
static inline void p_Setm(poly p, const ring r)
{
  p_CheckRing2(r);
  r->p_Setm(p, r);
}


static inline unsigned long p_AddComp(poly p, unsigned long v, ring r)
{
  p_LmCheckPolyRing2(p, r);
  pAssume2(rRing_has_Comp(r));
  return __p_GetComp(p,r) += v;
}
static inline unsigned long p_SubComp(poly p, unsigned long v, ring r)
{
  p_LmCheckPolyRing2(p, r);
  pAssume2(rRing_has_Comp(r));
  _pPolyAssume2(__p_GetComp(p,r) >= v,p,r);
  return __p_GetComp(p,r) -= v;
}

#ifndef HAVE_EXPSIZES

/// get a single variable exponent
/// @Note:
/// the integer VarOffset encodes:
/// 1. the position of a variable in the exponent vector p->exp (lower 24 bits)
/// 2. number of bits to shift to the right in the upper 8 bits (which takes at most 6 bits for 64 bit)
/// Thus VarOffset always has 2 zero higher bits!
static inline long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
{
  pAssume2((VarOffset >> (24 + 6)) == 0);
#if 0
  int pos=(VarOffset & 0xffffff);
  int bitpos=(VarOffset >> 24);
  unsigned long exp=(p->exp[pos] >> bitmask) & iBitmask;
  return exp;
#else
  return (long)
         ((p->exp[(VarOffset & 0xffffff)] >> (VarOffset >> 24))
          & iBitmask);
#endif
}


/// set a single variable exponent
/// @Note:
/// VarOffset encodes the position in p->exp @see p_GetExp
static inline unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
{
  pAssume2(e>=0);
  pAssume2(e<=iBitmask);
  pAssume2((VarOffset >> (24 + 6)) == 0);

  // shift e to the left:
  register int shift = VarOffset >> 24;
  unsigned long ee = e << shift /*(VarOffset >> 24)*/;
  // find the bits in the exponent vector
  register int offset = (VarOffset & 0xffffff);
  // clear the bits in the exponent vector:
  p->exp[offset]  &= ~( iBitmask << shift );
  // insert e with |
  p->exp[ offset ] |= ee;
  return e;
}


#else // #ifdef HAVE_EXPSIZES // EXPERIMENTAL!!!

static inline unsigned long BitMask(unsigned long bitmask, int twobits)
{
  // bitmask = 00000111111111111
  // 0 must give bitmask!
  // 1, 2, 3 - anything like 00011..11
  pAssume2((twobits >> 2) == 0);
  static const unsigned long _bitmasks[4] = {-1, 0x7fff, 0x7f, 0x3};
  return bitmask & _bitmasks[twobits];
}


/// @Note: we may add some more info (6 ) into VarOffset and thus encode
static inline long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
{
  int pos  =(VarOffset & 0xffffff);
  int hbyte= (VarOffset >> 24); // the highest byte
  int bitpos = hbyte & 0x3f; // last 6 bits
  long bitmask = BitMask(iBitmask, hbyte >> 6);

  long exp=(p->exp[pos] >> bitpos) & bitmask;
  return exp;

}

static inline long p_SetExp(poly p, const long e, const unsigned long iBitmask, const int VarOffset)
{
  pAssume2(e>=0);
  pAssume2(e <= BitMask(iBitmask, VarOffset >> 30));

  // shift e to the left:
  register int hbyte = VarOffset >> 24;
  int bitmask = BitMask(iBitmask, hbyte >> 6);
  register int shift = hbyte & 0x3f;
  long ee = e << shift;
  // find the bits in the exponent vector
  register int offset = (VarOffset & 0xffffff);
  // clear the bits in the exponent vector:
  p->exp[offset]  &= ~( bitmask << shift );
  // insert e with |
  p->exp[ offset ] |= ee;
  return e;
}

#endif // #ifndef HAVE_EXPSIZES


static inline long p_GetExp(const poly p, const ring r, const int VarOffset)
{
  p_LmCheckPolyRing2(p, r);
  pAssume2(VarOffset != -1);
  return p_GetExp(p, r->bitmask, VarOffset);
}

static inline long p_SetExp(poly p, const long e, const ring r, const int VarOffset)
{
  p_LmCheckPolyRing2(p, r);
  pAssume2(VarOffset != -1);
  return p_SetExp(p, e, r->bitmask, VarOffset);
}



/// get v^th exponent for a monomial
static inline long p_GetExp(const poly p, const int v, const ring r)
{
  p_LmCheckPolyRing2(p, r);
  pAssume2(v>0 && v <= r->N);
  pAssume2(r->VarOffset[v] != -1);
  return p_GetExp(p, r->bitmask, r->VarOffset[v]);
}


/// set v^th exponent for a monomial
static inline long p_SetExp(poly p, const int v, const long e, const ring r)
{
  p_LmCheckPolyRing2(p, r);
  pAssume2(v>0 && v <= r->N);
  pAssume2(r->VarOffset[v] != -1);
  return p_SetExp(p, e, r->bitmask, r->VarOffset[v]);
}

// the following should be implemented more efficiently
static inline  long p_IncrExp(poly p, int v, ring r)
{
  p_LmCheckPolyRing2(p, r);
  int e = p_GetExp(p,v,r);
  e++;
  return p_SetExp(p,v,e,r);
}
static inline  long p_DecrExp(poly p, int v, ring r)
{
  p_LmCheckPolyRing2(p, r);
  int e = p_GetExp(p,v,r);
  pAssume2(e > 0);
  e--;
  return p_SetExp(p,v,e,r);
}
static inline  long p_AddExp(poly p, int v, long ee, ring r)
{
  p_LmCheckPolyRing2(p, r);
  int e = p_GetExp(p,v,r);
  e += ee;
  return p_SetExp(p,v,e,r);
}
static inline  long p_SubExp(poly p, int v, long ee, ring r)
{
  p_LmCheckPolyRing2(p, r);
  long e = p_GetExp(p,v,r);
  pAssume2(e >= ee);
  e -= ee;
  return p_SetExp(p,v,e,r);
}
static inline  long p_MultExp(poly p, int v, long ee, ring r)
{
  p_LmCheckPolyRing2(p, r);
  long e = p_GetExp(p,v,r);
  e *= ee;
  return p_SetExp(p,v,e,r);
}

static inline long p_GetExpSum(poly p1, poly p2, int i, ring r)
{
  p_LmCheckPolyRing2(p1, r);
  p_LmCheckPolyRing2(p2, r);
  return p_GetExp(p1,i,r) + p_GetExp(p2,i,r);
}
static inline long p_GetExpDiff(poly p1, poly p2, int i, ring r)
{
  return p_GetExp(p1,i,r) - p_GetExp(p2,i,r);
}

static inline int p_Comp_k_n(poly a, poly b, int k, ring r)
{
  if ((a==NULL) || (b==NULL) ) return FALSE;
  p_LmCheckPolyRing2(a, r);
  p_LmCheckPolyRing2(b, r);
  pAssume2(k > 0 && k <= r->N);
  int i=k;
  for(;i<=r->N;i++)
  {
    if (p_GetExp(a,i,r) != p_GetExp(b,i,r)) return FALSE;
    //    if (a->exp[(r->VarOffset[i] & 0xffffff)] != b->exp[(r->VarOffset[i] & 0xffffff)]) return FALSE;
  }
  return TRUE;
}


/***************************************************************
 *
 * Allocation/Initalization/Deletion
 *
 ***************************************************************/
#if PDEBUG > 2
static inline poly p_New(const ring r, omBin bin)
#else
static inline poly p_New(const ring r, omBin bin)
#endif
{
  p_CheckRing2(r);
  pAssume2(bin != NULL && omSizeWOfBin(r->PolyBin) == omSizeWOfBin(bin));
  poly p;
  omTypeAllocBin(poly, p, bin);
  p_SetRingOfLm(p, r);
  return p;
}

static inline poly p_New(ring r)
{
  return p_New(r, r->PolyBin);
}

#if PDEBUG > 2
static inline void p_LmFree(poly p, ring r)
#else
static inline void p_LmFree(poly p, ring)
#endif
{
  p_LmCheckPolyRing2(p, r);
  omFreeBinAddr(p);
}
#if PDEBUG > 2
static inline void p_LmFree(poly *p, ring r)
#else
static inline void p_LmFree(poly *p, ring)
#endif
{
  p_LmCheckPolyRing2(*p, r);
  poly h = *p;
  *p = pNext(h);
  omFreeBinAddr(h);
}
#if PDEBUG > 2
static inline poly p_LmFreeAndNext(poly p, ring r)
#else
static inline poly p_LmFreeAndNext(poly p, ring)
#endif
{
  p_LmCheckPolyRing2(p, r);
  poly pnext = pNext(p);
  omFreeBinAddr(p);
  return pnext;
}
static inline void p_LmDelete(poly p, const ring r)
{
  p_LmCheckPolyRing2(p, r);
  n_Delete(&pGetCoeff(p), r->cf);
  omFreeBinAddr(p);
}
static inline void p_LmDelete(poly *p, const ring r)
{
  p_LmCheckPolyRing2(*p, r);
  poly h = *p;
  *p = pNext(h);
  n_Delete(&pGetCoeff(h), r->cf);
  omFreeBinAddr(h);
}
static inline poly p_LmDeleteAndNext(poly p, const ring r)
{
  p_LmCheckPolyRing2(p, r);
  poly pnext = pNext(p);
  n_Delete(&pGetCoeff(p), r->cf);
  omFreeBinAddr(p);
  return pnext;
}

/***************************************************************
 *
 * Misc routines
 *
 ***************************************************************/

/// return the maximal exponent of p in form of the maximal long var
unsigned long p_GetMaxExpL(poly p, const ring r, unsigned long l_max = 0);

/// return monomial r such that GetExp(r,i) is maximum of all
/// monomials in p; coeff == 0, next == NULL, ord is not set
poly p_GetMaxExpP(poly p, ring r);

static inline unsigned long p_GetMaxExp(const unsigned long l, const ring r)
{
  unsigned long bitmask = r->bitmask;
  unsigned long max = (l & bitmask);
  unsigned long j = r->ExpPerLong - 1;

  if (j > 0)
  {
    unsigned long i = r->BitsPerExp;
    long e;
    loop
    {
      e = ((l >> i) & bitmask);
      if ((unsigned long) e > max)
        max = e;
      j--;
      if (j==0) break;
      i += r->BitsPerExp;
    }
  }
  return max;
}

static inline unsigned long p_GetMaxExp(const poly p, const ring r)
{
  return p_GetMaxExp(p_GetMaxExpL(p, r), r);
}

static inline unsigned long
p_GetTotalDegree(const unsigned long l, const ring r, const int number_of_exps)
{
  const unsigned long bitmask = r->bitmask;
  unsigned long sum = (l & bitmask);
  unsigned long j = number_of_exps - 1;

  if (j > 0)
  {
    unsigned long i = r->BitsPerExp;
    loop
    {
      sum += ((l >> i) & bitmask);
      j--;
      if (j==0) break;
      i += r->BitsPerExp;
    }
  }
  return sum;
}

static inline unsigned long
p_GetTotalDegree(const unsigned long l, const ring r)
{
  return p_GetTotalDegree(l, r, r->ExpPerLong);
}

/***************************************************************
 *
 * Dispatcher to r->p_Procs, they do the tests/checks
 *
 ***************************************************************/
// returns a copy of p
static inline poly p_Copy(poly p, const ring r)
{
#ifdef PDEBUG
  poly pp= r->p_Procs->p_Copy(p, r);
  p_Test(pp,r);
  return pp;
#else
  return r->p_Procs->p_Copy(p, r);
#endif
}

static inline poly p_Head(poly p, const ring r)
{
  if (p == NULL) return NULL;
  p_LmCheckPolyRing1(p, r);
  poly np;
  omTypeAllocBin(poly, np, r->PolyBin);
  p_SetRingOfLm(np, r);
  memcpy(np->exp, p->exp, r->ExpL_Size*sizeof(long));
  pNext(np) = NULL;
  pSetCoeff0(np, n_Copy(pGetCoeff(p), r->cf));
  return np;
}

// returns a copy of p with Lm(p) from lmRing and Tail(p) from tailRing
static inline poly p_Copy(poly p, const ring lmRing, const ring tailRing)
{
#ifndef PDEBUG
  if (tailRing == lmRing)
    return tailRing->p_Procs->p_Copy(p, tailRing);
#endif
  if (p != NULL)
  {
    poly pres = p_Head(p, lmRing);
    pNext(pres) = tailRing->p_Procs->p_Copy(pNext(p), tailRing);
    return pres;
  }
  else
    return NULL;
}

// deletes *p, and sets *p to NULL
static inline void p_Delete(poly *p, const ring r)
{
  r->p_Procs->p_Delete(p, r);
}

static inline void p_Delete(poly *p,  const ring lmRing, const ring tailRing)
{
#ifndef PDEBUG
  if (tailRing == lmRing)
  {
    tailRing->p_Procs->p_Delete(p, tailRing);
    return;
  }
#endif
  if (*p != NULL)
  {
    if (pNext(*p) != NULL)
      tailRing->p_Procs->p_Delete(&pNext(*p), tailRing);
    p_LmDelete(p, lmRing);
  }
}

// copys monomials of p, allocates new monomials from bin,
// deletes monomoals of p
static inline poly p_ShallowCopyDelete(poly p, const ring r, omBin bin)
{
  p_LmCheckPolyRing2(p, r);
  pAssume2(omSizeWOfBin(r->PolyBin) == omSizeWOfBin(bin));
  return r->p_Procs->p_ShallowCopyDelete(p, r, bin);
}

// returns p+q, destroys p and q
static inline poly p_Add_q(poly p, poly q, const ring r)
{
  assume( (p != q) || (p == NULL && q == NULL) );
  int shorter;
  return r->p_Procs->p_Add_q(p, q, shorter, r);
}

/// like p_Add_q, except that if lp == pLength(lp) lq == pLength(lq) then lp == pLength(p+q)
static inline poly p_Add_q(poly p, poly q, int &lp, int lq, const ring r)
{
  assume( (p != q) || (p == NULL && q == NULL) );
  int shorter;
  poly res = r->p_Procs->p_Add_q(p, q, shorter, r);
  lp = (lp + lq) - shorter;
  return res;
}

// returns p*n, destroys p
static inline poly p_Mult_nn(poly p, number n, const ring r)
{
  if (n_IsOne(n, r->cf))
    return p;
  else
    return r->p_Procs->p_Mult_nn(p, n, r);
}

static inline poly p_Mult_nn(poly p, number n, const ring lmRing,
                        const ring tailRing)
{
#ifndef PDEBUG
  if (lmRing == tailRing)
  {
    return p_Mult_nn(p, n, tailRing);
  }
#endif
  poly pnext = pNext(p);
  pNext(p) = NULL;
  p = lmRing->p_Procs->p_Mult_nn(p, n, lmRing);
  pNext(p) = tailRing->p_Procs->p_Mult_nn(pnext, n, tailRing);
  return p;
}

// returns p*n, does not destroy p
static inline poly pp_Mult_nn(poly p, number n, const ring r)
{
  if (n_IsOne(n, r->cf))
    return p_Copy(p, r);
  else
    return r->p_Procs->pp_Mult_nn(p, n, r);
}

// test if the monomial is a constant as a vector component
// i.e., test if all exponents are zero
static inline BOOLEAN p_LmIsConstantComp(const poly p, const ring r)
{
  //p_LmCheckPolyRing(p, r);
  int i = r->VarL_Size - 1;

  do
  {
    if (p->exp[r->VarL_Offset[i]] != 0)
      return FALSE;
    i--;
  }
  while (i >= 0);
  return TRUE;
}

// test if monomial is a constant, i.e. if all exponents and the component
// is zero
static inline BOOLEAN p_LmIsConstant(const poly p, const ring r)
{
  if (p_LmIsConstantComp(p, r))
    return (p_GetComp(p, r) == 0);
  return FALSE;
}

// returns Copy(p)*m, does neither destroy p nor m
static inline poly pp_Mult_mm(poly p, poly m, const ring r)
{
  if (p_LmIsConstant(m, r))
    return pp_Mult_nn(p, pGetCoeff(m), r);
  else
  {
    poly last;
    return r->p_Procs->pp_Mult_mm(p, m, r, last);
  }
}

// returns p*m, destroys p, const: m
static inline poly p_Mult_mm(poly p, poly m, const ring r)
{
  if (p_LmIsConstant(m, r))
    return p_Mult_nn(p, pGetCoeff(m), r);
  else
    return r->p_Procs->p_Mult_mm(p, m, r);
}

// return p - m*Copy(q), destroys p; const: p,m
static inline poly p_Minus_mm_Mult_qq(poly p, poly m, poly q, const ring r)
{
#ifdef HAVE_PLURAL
  if (rIsPluralRing(r))
  {
    int lp, lq = 0;
    poly spNoether = NULL;
    return nc_p_Minus_mm_Mult_qq(p, m, q, lp, lq, spNoether, r);
  }
#endif

  int shorter;
  poly last;

  return r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, NULL, r, last); // !!!
}

// like p_Minus_mm_Mult_qq, except that if lp == pLength(lp) lq == pLength(lq)
// then lp == pLength(p -m*q)
static inline poly p_Minus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq,
                                 poly spNoether, const ring r)
{
#ifdef HAVE_PLURAL
  if (rIsPluralRing(r))
     return nc_p_Minus_mm_Mult_qq(p, m, q, lp, lq, spNoether, r);
#endif

  int shorter;
  poly last,res;
  res = r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, spNoether, r, last);
  lp = (lp + lq) - shorter;
  return res;
}

// returns p*Coeff(m) for such monomials pm of p, for which m is divisble by pm
static inline poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r)
{
  int shorter;
  return r->p_Procs->pp_Mult_Coeff_mm_DivSelect(p, m, shorter, r);
}

// returns p*Coeff(m) for such monomials pm of p, for which m is divisble by pm
// if lp is length of p on input then lp is length of returned poly on output
static inline poly pp_Mult_Coeff_mm_DivSelect(poly p, int &lp, const poly m, const ring r)
{
  int shorter;
  poly pp = r->p_Procs->pp_Mult_Coeff_mm_DivSelect(p, m, shorter, r);
  lp -= shorter;
  return pp;
}

// returns -p, destroys p
static inline poly p_Neg(poly p, const ring r)
{
  return r->p_Procs->p_Neg(p, r);
}

extern poly  _p_Mult_q(poly p, poly q, const int copy, const ring r);
// returns p*q, destroys p and q
static inline poly p_Mult_q(poly p, poly q, const ring r)
{
  assume( (p != q) || (p == NULL && q == NULL) );

  if (p == NULL)
  {
    r->p_Procs->p_Delete(&q, r);
    return NULL;
  }
  if (q == NULL)
  {
    r->p_Procs->p_Delete(&p, r);
    return NULL;
  }

  if (pNext(p) == NULL)
  {
#ifdef HAVE_PLURAL
    if (rIsPluralRing(r))
      q = nc_mm_Mult_p(p, q, r);
    else
#endif /* HAVE_PLURAL */
      q = r->p_Procs->p_Mult_mm(q, p, r);

    r->p_Procs->p_Delete(&p, r);
    return q;
  }

  if (pNext(q) == NULL)
  {
  // NEEDED
#ifdef HAVE_PLURAL
/*    if (rIsPluralRing(r))
      p = gnc_p_Mult_mm(p, q, r); // ???
    else*/
#endif /* HAVE_PLURAL */
      p = r->p_Procs->p_Mult_mm(p, q, r);

    r->p_Procs->p_Delete(&q, r);
    return p;
  }
#ifdef HAVE_PLURAL
  if (rIsPluralRing(r))
    return _nc_p_Mult_q(p, q, r);
  else
#endif
  return _p_Mult_q(p, q, 0, r);
}

// returns p*q, does neither destroy p nor q
static inline poly pp_Mult_qq(poly p, poly q, const ring r)
{
  poly last;
  if (p == NULL || q == NULL) return NULL;

  if (pNext(p) == NULL)
  {
#ifdef HAVE_PLURAL
    if (rIsPluralRing(r))
      return nc_mm_Mult_pp(p, q, r);
#endif
    return r->p_Procs->pp_Mult_mm(q, p, r, last);
  }

  if (pNext(q) == NULL)
  {
    return r->p_Procs->pp_Mult_mm(p, q, r, last);
  }

  poly qq = q;
  if (p == q)
    qq = p_Copy(q, r);

  poly res;
#ifdef HAVE_PLURAL
  if (rIsPluralRing(r))
    res = _nc_pp_Mult_qq(p, qq, r);
  else
#endif
    res = _p_Mult_q(p, qq, 1, r);

  if (qq != q)
    p_Delete(&qq, r);
  return res;
}

// returns p + m*q destroys p, const: q, m
static inline poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq,
                                const ring r)
{
#ifdef HAVE_PLURAL
  if (rIsPluralRing(r))
    return nc_p_Plus_mm_Mult_qq(p, m, q, lp, lq, r);
#endif

// this should be implemented more efficiently
  poly res, last;
  int shorter;
  number n_old = pGetCoeff(m);
  number n_neg = n_Copy(n_old, r->cf);
  n_neg = n_Neg(n_neg, r->cf);
  pSetCoeff0(m, n_neg);
  res = r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, NULL, r, last);
  lp = (lp + lq) - shorter;
  pSetCoeff0(m, n_old);
  n_Delete(&n_neg, r->cf);
  return res;
}

static inline poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, const ring r)
{
  int lp = 0, lq = 0;
  return p_Plus_mm_Mult_qq(p, m, q, lp, lq, r);
}

// returns merged p and q, assumes p and q have no monomials which are equal
static inline poly p_Merge_q(poly p, poly q, const ring r)
{
  assume( (p != q) || (p == NULL && q == NULL) );
  return r->p_Procs->p_Merge_q(p, q, r);
}

// like p_SortMerge, except that p may have equal monimals
static inline poly p_SortAdd(poly p, const ring r, BOOLEAN revert= FALSE)
{
  if (revert) p = pReverse(p);
  return sBucketSortAdd(p, r);
}

// sorts p using bucket sort: returns sorted poly
// assumes that monomials of p are all different
// reverses it first, if revert == TRUE, use this if input p is "almost" sorted
// correctly
static inline poly p_SortMerge(poly p, const ring r, BOOLEAN revert= FALSE)
{
  if (revert) p = pReverse(p);
  return sBucketSortMerge(p, r);
}

/***************************************************************
 *
 * I/O
 *
 ***************************************************************/
static inline char*     p_String(poly p, ring p_ring)
{
  return p_String(p, p_ring, p_ring);
}
static inline char*     p_String0(poly p, ring p_ring)
{
  return p_String0(p, p_ring, p_ring);
}
static inline void      p_Write(poly p, ring p_ring)
{
  p_Write(p, p_ring, p_ring);
}
static inline void      p_Write0(poly p, ring p_ring)
{
  p_Write0(p, p_ring, p_ring);
}
static inline void      p_wrp(poly p, ring p_ring)
{
  p_wrp(p, p_ring, p_ring);
}


#if PDEBUG > 0

#define _p_LmCmpAction(p, q, r, actionE, actionG, actionS)  \
do                                                          \
{                                                           \
  int _cmp = p_LmCmp(p,q,r);                                \
  if (_cmp == 0) actionE;                                   \
  if (_cmp == 1) actionG;                                   \
  actionS;                                                  \
}                                                           \
while(0)

#else

#define _p_LmCmpAction(p, q, r, actionE, actionG, actionS)                      \
 p_MemCmp_LengthGeneral_OrdGeneral(p->exp, q->exp, r->CmpL_Size, r->ordsgn,    \
                                   actionE, actionG, actionS)

#endif

#define pDivAssume(x)   ((void)0)



/***************************************************************
 *
 * Allocation/Initalization/Deletion
 *
 ***************************************************************/
// adjustments for negative weights
static inline void p_MemAdd_NegWeightAdjust(poly p, const ring r)
{
  if (r->NegWeightL_Offset != NULL)
  {
    for (int i=r->NegWeightL_Size-1; i>=0; i--)
    {
      p->exp[r->NegWeightL_Offset[i]] -= POLY_NEGWEIGHT_OFFSET;
    }
  }
}
static inline void p_MemSub_NegWeightAdjust(poly p, const ring r)
{
  if (r->NegWeightL_Offset != NULL)
  {
    for (int i=r->NegWeightL_Size-1; i>=0; i--)
    {
      p->exp[r->NegWeightL_Offset[i]] += POLY_NEGWEIGHT_OFFSET;
    }
  }
}
// ExpVextor(d_p) = ExpVector(s_p)
static inline void p_ExpVectorCopy(poly d_p, poly s_p, const ring r)
{
  p_LmCheckPolyRing1(d_p, r);
  p_LmCheckPolyRing1(s_p, r);
  memcpy(d_p->exp, s_p->exp, r->ExpL_Size*sizeof(long));
}

static inline poly p_Init(const ring r, omBin bin)
{
  p_CheckRing1(r);
  pAssume1(bin != NULL && omSizeWOfBin(r->PolyBin) == omSizeWOfBin(bin));
  poly p;
  omTypeAlloc0Bin(poly, p, bin);
  p_MemAdd_NegWeightAdjust(p, r);
  p_SetRingOfLm(p, r);
  return p;
}
static inline poly p_Init(const ring r)
{
  return p_Init(r, r->PolyBin);
}

static inline poly p_LmInit(poly p, const ring r)
{
  p_LmCheckPolyRing1(p, r);
  poly np;
  omTypeAllocBin(poly, np, r->PolyBin);
  p_SetRingOfLm(np, r);
  memcpy(np->exp, p->exp, r->ExpL_Size*sizeof(long));
  pNext(np) = NULL;
  pSetCoeff0(np, NULL);
  return np;
}
static inline poly p_LmInit(poly s_p, const ring s_r, const ring d_r, omBin d_bin)
{
  p_LmCheckPolyRing1(s_p, s_r);
  p_CheckRing(d_r);
  pAssume1(d_r->N <= s_r->N);
  poly d_p = p_Init(d_r, d_bin);
  for (int i=d_r->N; i>0; i--)
  {
    p_SetExp(d_p, i, p_GetExp(s_p, i,s_r), d_r);
  }
  if (rRing_has_Comp(d_r))
  {
    p_SetComp(d_p, p_GetComp(s_p,s_r), d_r);
  }
  p_Setm(d_p, d_r);
  return d_p;
}
static inline poly p_LmInit(poly s_p, const ring s_r, const ring d_r)
{
  pAssume1(d_r != NULL);
  return p_LmInit(s_p, s_r, d_r, d_r->PolyBin);
}

// set all exponents l..k to 0, assume exp. k+1..n and 1..l-1 are in
// different blocks
// set coeff to 1
static inline poly p_GetExp_k_n(poly p, int l, int k, const ring r)
{
  if (p == NULL) return NULL;
  p_LmCheckPolyRing1(p, r);
  poly np;
  omTypeAllocBin(poly, np, r->PolyBin);
  p_SetRingOfLm(np, r);
  memcpy(np->exp, p->exp, r->ExpL_Size*sizeof(long));
  pNext(np) = NULL;
  pSetCoeff0(np, n_Init(1, r->cf));
  int i;
  for(i=l;i<=k;i++)
  {
    //np->exp[(r->VarOffset[i] & 0xffffff)] =0;
    p_SetExp(np,i,0,r);
  }
  p_Setm(np,r);
  return np;
}

// simialar to p_ShallowCopyDelete but does it only for leading monomial
static inline poly p_LmShallowCopyDelete(poly p, const ring r)
{
  p_LmCheckPolyRing1(p, r);
  pAssume1(omSizeWOfBin(bin) == omSizeWOfBin(r->PolyBin));
  poly new_p = p_New(r);
  memcpy(new_p->exp, p->exp, r->ExpL_Size*sizeof(long));
  pSetCoeff0(new_p, pGetCoeff(p));
  pNext(new_p) = pNext(p);
  omFreeBinAddr(p);
  return new_p;
}

/***************************************************************
 *
 * Operation on ExpVectors
 *
 ***************************************************************/
// ExpVector(p1) += ExpVector(p2)
static inline void p_ExpVectorAdd(poly p1, poly p2, const ring r)
{
  p_LmCheckPolyRing1(p1, r);
  p_LmCheckPolyRing1(p2, r);
#if PDEBUG >= 1
  for (int i=1; i<=r->N; i++)
    pAssume1((unsigned long) (p_GetExp(p1, i, r) + p_GetExp(p2, i, r)) <= r->bitmask);
  pAssume1(p_GetComp(p1, r) == 0 || p_GetComp(p2, r) == 0);
#endif

  p_MemAdd_LengthGeneral(p1->exp, p2->exp, r->ExpL_Size);
  p_MemAdd_NegWeightAdjust(p1, r);
}
// ExpVector(pr) = ExpVector(p1) + ExpVector(p2)
static inline void p_ExpVectorSum(poly pr, poly p1, poly p2, const ring r)
{
  p_LmCheckPolyRing1(p1, r);
  p_LmCheckPolyRing1(p2, r);
  p_LmCheckPolyRing1(pr, r);
#if PDEBUG >= 1
  for (int i=1; i<=r->N; i++)
    pAssume1((unsigned long) (p_GetExp(p1, i, r) + p_GetExp(p2, i, r)) <= r->bitmask);
  pAssume1(p_GetComp(p1, r) == 0 || p_GetComp(p2, r) == 0);
#endif

  p_MemSum_LengthGeneral(pr->exp, p1->exp, p2->exp, r->ExpL_Size);
  p_MemAdd_NegWeightAdjust(pr, r);
}
// ExpVector(p1) -= ExpVector(p2)
static inline void p_ExpVectorSub(poly p1, poly p2, const ring r)
{
  p_LmCheckPolyRing1(p1, r);
  p_LmCheckPolyRing1(p2, r);
#if PDEBUG >= 1
  for (int i=1; i<=r->N; i++)
    pAssume1(p_GetExp(p1, i, r) >= p_GetExp(p2, i, r));
  pAssume1(p_GetComp(p1, r) == 0 || p_GetComp(p2, r) == 0 ||
          p_GetComp(p1, r) == p_GetComp(p2, r));
#endif

  p_MemSub_LengthGeneral(p1->exp, p2->exp, r->ExpL_Size);
  p_MemSub_NegWeightAdjust(p1, r);

}
// ExpVector(p1) += ExpVector(p2) - ExpVector(p3)
//static inline void p_ExpVectorAddSub(poly p1, poly p2, poly p3, const ring r)
//{
//  p_LmCheckPolyRing1(p1, r);
//  p_LmCheckPolyRing1(p2, r);
//  p_LmCheckPolyRing1(p3, r);
//#if PDEBUG >= 1
//  for (int i=1; i<=r->N; i++)
//    pAssume1(p_GetExp(p1, i, r) + p_GetExp(p2, i, r) >= p_GetExp(p3, i, r));
//  pAssume1(p_GetComp(p1, r) == 0 ||
//           (p_GetComp(p2, r) - p_GetComp(p3, r) == 0) ||
//           (p_GetComp(p1, r) == p_GetComp(p2, r) - p_GetComp(p3, r)));
//#endif
//
//  p_MemAddSub_LengthGeneral(p1->exp, p2->exp, p3->exp, r->ExpL_Size);
//  // no need to adjust in case of NegWeights
//}

// ExpVector(pr) = ExpVector(p1) - ExpVector(p2)
static inline void p_ExpVectorDiff(poly pr, poly p1, poly p2, const ring r)
{
  p_LmCheckPolyRing1(p1, r);
  p_LmCheckPolyRing1(p2, r);
  p_LmCheckPolyRing1(pr, r);
#if PDEBUG >= 2
  for (int i=1; i<=r->N; i++)
    pAssume1(p_GetExp(p1, i, r) >= p_GetExp(p2, i, r));
  pAssume1(!rRing_has_Comp(r) || p_GetComp(p1, r) == p_GetComp(p2, r));
#endif

  p_MemDiff_LengthGeneral(pr->exp, p1->exp, p2->exp, r->ExpL_Size);
  p_MemSub_NegWeightAdjust(pr, r);
}

static inline BOOLEAN p_ExpVectorEqual(poly p1, poly p2, const ring r)
{
  p_LmCheckPolyRing1(p1, r);
  p_LmCheckPolyRing1(p2, r);

  int i = r->ExpL_Size;
  unsigned long *ep = p1->exp;
  unsigned long *eq = p2->exp;

  do
  {
    i--;
    if (ep[i] != eq[i]) return FALSE;
  }
  while (i);
  return TRUE;
}

static inline long p_Totaldegree(poly p, const ring r)
{
  p_LmCheckPolyRing1(p, r);
  unsigned long s = p_GetTotalDegree(p->exp[r->VarL_Offset[0]],
                                     r,
                                     r->MinExpPerLong);
  for (int i=r->VarL_Size-1; i>0; i--)
  {
    s += p_GetTotalDegree(p->exp[r->VarL_Offset[i]], r);
  }
  return (long)s;
}

static inline void p_GetExpV(poly p, int *ev, const ring r)
{
  p_LmCheckPolyRing1(p, r);
  for (int j = r->N; j; j--)
      ev[j] = p_GetExp(p, j, r);

  ev[0] = p_GetComp(p, r);
}
static inline void p_SetExpV(poly p, int *ev, const ring r)
{
  p_LmCheckPolyRing1(p, r);
  for (int j = r->N; j; j--)
      p_SetExp(p, j, ev[j], r);

  p_SetComp(p, ev[0],r);
  p_Setm(p, r);
}

/***************************************************************
 *
 * Comparison w.r.t. monomial ordering
 *
 ***************************************************************/

static inline int p_LmCmp(poly p, poly q, const ring r)
{
  p_LmCheckPolyRing1(p, r);
  p_LmCheckPolyRing1(q, r);

  const unsigned long* _s1 = ((unsigned long*) p->exp);
  const unsigned long* _s2 = ((unsigned long*) q->exp);
  register unsigned long _v1;
  register unsigned long _v2;
  const unsigned long _l = r->CmpL_Size;

  register unsigned long _i=0;

  LengthGeneral_OrdGeneral_LoopTop:
  _v1 = _s1[_i];
  _v2 = _s2[_i];
  if (_v1 == _v2)
  {
    _i++;
    if (_i == _l) return 0;
    goto LengthGeneral_OrdGeneral_LoopTop;
  }
  const long* _ordsgn = (long*) r->ordsgn;
  if (_v1 > _v2)
  {
    if (_ordsgn[_i] == 1) return 1;
    return -1;
  }
  if (_ordsgn[_i] == 1) return -1;
  return 1;

}

/// returns TRUE if p1 is a skalar multiple of p2
/// assume p1 != NULL and p2 != NULL
BOOLEAN p_ComparePolys(poly p1,poly p2, const ring r);


/***************************************************************
 *
 * Comparisons: they are all done without regarding coeffs
 *
 ***************************************************************/
#define p_LmCmpAction(p, q, r, actionE, actionG, actionS) \
  _p_LmCmpAction(p, q, r, actionE, actionG, actionS)

// returns 1 if ExpVector(p)==ExpVector(q): does not compare numbers !!
#define p_LmEqual(p1, p2, r) p_ExpVectorEqual(p1, p2, r)

// pCmp: args may be NULL
// returns: (p2==NULL ? 1 : (p1 == NULL ? -1 : p_LmCmp(p1, p2)))
static inline int p_Cmp(poly p1, poly p2, ring r)
{
  if (p2==NULL)
    return 1;
  if (p1==NULL)
    return -1;
  return p_LmCmp(p1,p2,r);
}


/***************************************************************
 *
 * divisibility
 *
 ***************************************************************/
/// return: FALSE, if there exists i, such that a->exp[i] > b->exp[i]
///         TRUE, otherwise
/// (1) Consider long vars, instead of single exponents
/// (2) Clearly, if la > lb, then FALSE
/// (3) Suppose la <= lb, and consider first bits of single exponents in l:
///     if TRUE, then value of these bits is la ^ lb
///     if FALSE, then la-lb causes an "overflow" into one of those bits, i.e.,
///               la ^ lb != la - lb
static inline BOOLEAN _p_LmDivisibleByNoComp(poly a, poly b, const ring r)
{
  int i=r->VarL_Size - 1;
  unsigned long divmask = r->divmask;
  unsigned long la, lb;

  if (r->VarL_LowIndex >= 0)
  {
    i += r->VarL_LowIndex;
    do
    {
      la = a->exp[i];
      lb = b->exp[i];
      if ((la > lb) ||
          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
      {
        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
        return FALSE;
      }
      i--;
    }
    while (i>=r->VarL_LowIndex);
  }
  else
  {
    do
    {
      la = a->exp[r->VarL_Offset[i]];
      lb = b->exp[r->VarL_Offset[i]];
      if ((la > lb) ||
          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
      {
        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
        return FALSE;
      }
      i--;
    }
    while (i>=0);
  }
#ifdef HAVE_RINGS
  pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == nDivBy(p_GetCoeff(b, r), p_GetCoeff(a, r)));
  return (!rField_is_Ring(r)) || n_DivBy(p_GetCoeff(b, r), p_GetCoeff(a, r), r->cf);
#else
  pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == TRUE);
  return TRUE;
#endif
}

static inline BOOLEAN _p_LmDivisibleByNoComp(poly a, const ring r_a, poly b, const ring r_b)
{
  int i=r_a->N;
  pAssume1(r_a->N == r_b->N);

  do
  {
    if (p_GetExp(a,i,r_a) > p_GetExp(b,i,r_b))
      return FALSE;
    i--;
  }
  while (i);
#ifdef HAVE_RINGS
  return n_DivBy(p_GetCoeff(b, r_b), p_GetCoeff(a, r_a), r_a->cf);
#else
  return TRUE;
#endif
}

#ifdef HAVE_RATGRING
static inline BOOLEAN _p_LmDivisibleByNoCompPart(poly a, const ring r_a, poly b, const ring r_b,const int start, const int end)
{
  int i=end;
  pAssume1(r_a->N == r_b->N);

  do
  {
    if (p_GetExp(a,i,r_a) > p_GetExp(b,i,r_b))
      return FALSE;
    i--;
  }
  while (i>=start);
#ifdef HAVE_RINGS
  return nDivBy(p_GetCoeff(b, r), p_GetCoeff(a, r));
#else
  return TRUE;
#endif
}
static inline BOOLEAN _p_LmDivisibleByPart(poly a, const ring r_a, poly b, const ring r_b,const int start, const int end)
{
  if (p_GetComp(a, r_a) == 0 || p_GetComp(a,r_a) == p_GetComp(b,r_b))
    return _p_LmDivisibleByNoCompPart(a, r_a, b, r_b,start,end);
  return FALSE;
}
static inline BOOLEAN p_LmDivisibleByPart(poly a, poly b, const ring r,const int start, const int end)
{
  p_LmCheckPolyRing1(b, r);
  pIfThen1(a != NULL, p_LmCheckPolyRing1(b, r));
  if (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r))
    return _p_LmDivisibleByNoCompPart(a, r, b, r,start, end);
  return FALSE;
}
#endif
static inline BOOLEAN _p_LmDivisibleBy(poly a, poly b, const ring r)
{
  if (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r))
    return _p_LmDivisibleByNoComp(a, b, r);
  return FALSE;
}
static inline BOOLEAN _p_LmDivisibleBy(poly a, const ring r_a, poly b, const ring r_b)
{
  if (p_GetComp(a, r_a) == 0 || p_GetComp(a,r_a) == p_GetComp(b,r_b))
    return _p_LmDivisibleByNoComp(a, r_a, b, r_b);
  return FALSE;
}
static inline BOOLEAN p_LmDivisibleByNoComp(poly a, poly b, const ring r)
{
  p_LmCheckPolyRing1(a, r);
  p_LmCheckPolyRing1(b, r);
  return _p_LmDivisibleByNoComp(a, b, r);
}

static inline BOOLEAN p_LmDivisibleByNoComp(poly a, const ring ra, poly b, const ring rb)
{
  p_LmCheckPolyRing1(a, ra);
  p_LmCheckPolyRing1(b, rb);
  return _p_LmDivisibleByNoComp(a, ra, b, rb);
}

static inline BOOLEAN p_LmDivisibleBy(poly a, poly b, const ring r)
{
  p_LmCheckPolyRing1(b, r);
  pIfThen1(a != NULL, p_LmCheckPolyRing1(b, r));
  if (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r))
    return _p_LmDivisibleByNoComp(a, b, r);
  return FALSE;
}

static inline BOOLEAN p_DivisibleBy(poly a, poly b, const ring r)
{
  pIfThen1(b!=NULL, p_LmCheckPolyRing1(b, r));
  pIfThen1(a!=NULL, p_LmCheckPolyRing1(a, r));

  if (a != NULL && (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r)))
      return _p_LmDivisibleByNoComp(a,b,r);
  return FALSE;
}
static inline BOOLEAN p_DivisibleBy(poly a, const ring r_a, poly b, const ring r_b)
{
  pIfThen1(b!=NULL, p_LmCheckPolyRing1(b, r_b));
  pIfThen1(a!=NULL, p_LmCheckPolyRing1(a, r_a));
  if (a != NULL) {
      return _p_LmDivisibleBy(a, r_a, b, r_b);
  }
  return FALSE;
}
static inline BOOLEAN p_LmDivisibleBy(poly a, const ring r_a, poly b, const ring r_b)
{
  p_LmCheckPolyRing(a, r_a);
  p_LmCheckPolyRing(b, r_b);
  return _p_LmDivisibleBy(a, r_a, b, r_b);
}

static inline BOOLEAN p_LmShortDivisibleBy(poly a, unsigned long sev_a,
                                    poly b, unsigned long not_sev_b, const ring r)
{
  p_LmCheckPolyRing1(a, r);
  p_LmCheckPolyRing1(b, r);
#ifndef PDIV_DEBUG
  _pPolyAssume2(p_GetShortExpVector(a, r) == sev_a, a, r);
  _pPolyAssume2(p_GetShortExpVector(b, r) == ~ not_sev_b, b, r);

  if (sev_a & not_sev_b)
  {
    pAssume1(p_LmDivisibleByNoComp(a, b, r) == FALSE);
    return FALSE;
  }
  return p_LmDivisibleBy(a, b, r);
#else
  return pDebugLmShortDivisibleBy(a, sev_a, r, b, not_sev_b, r);
#endif
}

static inline BOOLEAN p_LmShortDivisibleByNoComp(poly a, unsigned long sev_a,
                                           poly b, unsigned long not_sev_b, const ring r)
{
  p_LmCheckPolyRing1(a, r);
  p_LmCheckPolyRing1(b, r);
#ifndef PDIV_DEBUG
  _pPolyAssume2(p_GetShortExpVector(a, r) == sev_a, a, r);
  _pPolyAssume2(p_GetShortExpVector(b, r) == ~ not_sev_b, b, r);

  if (sev_a & not_sev_b)
  {
    pAssume1(p_LmDivisibleByNoComp(a, b, r) == FALSE);
    return FALSE;
  }
  return p_LmDivisibleByNoComp(a, b, r);
#else
  return pDebugLmShortDivisibleByNoComp(a, sev_a, r, b, not_sev_b, r);
#endif
}

static inline BOOLEAN p_LmShortDivisibleBy(poly a, unsigned long sev_a, const ring r_a,
                                      poly b, unsigned long not_sev_b, const ring r_b)
{
  p_LmCheckPolyRing1(a, r_a);
  p_LmCheckPolyRing1(b, r_b);
#ifndef PDIV_DEBUG
  _pPolyAssume2(p_GetShortExpVector(a, r_a) == sev_a, a, r_a);
  _pPolyAssume2(p_GetShortExpVector(b, r_b) == ~ not_sev_b, b, r_b);

  if (sev_a & not_sev_b)
  {
    pAssume1(_p_LmDivisibleByNoComp(a, r_a, b, r_b) == FALSE);
    return FALSE;
  }
  return _p_LmDivisibleBy(a, r_a, b, r_b);
#else
  return pDebugLmShortDivisibleBy(a, sev_a, r_a, b, not_sev_b, r_b);
#endif
}

/***************************************************************
 *
 * Misc things on Lm
 *
 ***************************************************************/


// like the respective p_LmIs* routines, except that p might be empty
static inline BOOLEAN p_IsConstantComp(const poly p, const ring r)
{
  if (p == NULL) return TRUE;
  return (pNext(p)==NULL) && p_LmIsConstantComp(p, r);
}

static inline BOOLEAN p_IsConstant(const poly p, const ring r)
{
  if (p == NULL) return TRUE;
  return (pNext(p)==NULL) && p_LmIsConstant(p, r);
}

static inline BOOLEAN p_IsConstantPoly(const poly p, const ring r)
{
  poly pp=p;
  while(pp!=NULL)
  {
    if (! p_LmIsConstantComp(pp, r))
      return FALSE;
    pIter(pp);
  }
  return TRUE;
}

static inline BOOLEAN p_IsUnit(const poly p, const ring r)
{
  if (p == NULL) return FALSE;
#ifdef HAVE_RINGS
  if (rField_is_Ring(r))
    return (p_LmIsConstant(p, r) && n_IsUnit(pGetCoeff(p),r->cf));
#endif
  return p_LmIsConstant(p, r);
}

static inline BOOLEAN p_LmExpVectorAddIsOk(const poly p1, const poly p2,
                                      const ring r)
{
  p_LmCheckPolyRing(p1, r);
  p_LmCheckPolyRing(p2, r);
  unsigned long l1, l2, divmask = r->divmask;
  int i;

  for (i=0; i<r->VarL_Size; i++)
  {
    l1 = p1->exp[r->VarL_Offset[i]];
    l2 = p2->exp[r->VarL_Offset[i]];
    // do the divisiblity trick
    if ( (l1 > ULONG_MAX - l2) ||
         (((l1 & divmask) ^ (l2 & divmask)) != ((l1 + l2) & divmask)))
      return FALSE;
  }
  return TRUE;
}
void      p_Split(poly p, poly * r);   /*p => IN(p), r => REST(p) */
BOOLEAN p_HasNotCF(poly p1, poly p2, const ring r);
poly      p_mInit(const char *s, BOOLEAN &ok, const ring r); /* monom s -> poly, interpreter */
const char *    p_Read(const char *s, poly &p,const ring r); /* monom -> poly */
poly      p_Divide(poly a, poly b, const ring r);
poly      p_DivideM(poly a, poly b, const ring r);
poly      p_Div_nn(poly p, const number n, const ring r);
void      p_Lcm(poly a, poly b, poly m, const ring r);
poly      p_Diff(poly a, int k, const ring r);
poly      p_DiffOp(poly a, poly b,BOOLEAN multiply, const ring r);
int       p_Weight(int c, const ring r);

/* assumes that p and divisor are univariate polynomials in r,
   mentioning the same variable;
   assumes divisor != NULL;
   p may be NULL;
   assumes a global monomial ordering in r;
   performs polynomial division of p by divisor:
     - afterwards p contains the remainder of the division, i.e.,
       p_before = result * divisor + p_afterwards;
     - if needResult == TRUE, then the method computes and returns 'result',
       otherwise NULL is returned (This parametrization can be used when
       one is only interested in the remainder of the division. In this
       case, the method will be slightly faster.)
   leaves divisor unmodified */
poly      p_PolyDiv(poly &p, poly divisor, BOOLEAN needResult, ring r);

/* returns NULL if p == NULL, otherwise makes p monic by dividing
   by its leading coefficient (only done if this is not already 1);
   this assumes that we are over a ground field so that division
   is well-defined;
   modifies p */
void      p_Monic(poly p, const ring r);

/* assumes that p and q are univariate polynomials in r,
   mentioning the same variable;
   assumes a global monomial ordering in r;
   assumes that not both p and q are NULL;
   returns the gcd of p and q;
   leaves p and q unmodified */
poly      p_Gcd(poly p, poly q, ring r);

/* assumes that p and q are univariate polynomials in r,
   mentioning the same variable;
   assumes a global monomial ordering in r;
   assumes that not both p and q are NULL;
   returns the gcd of p and q;
   moreover, afterwards pFactor and qFactor contain appropriate
   factors such that gcd(p, q) = p * pFactor + q * qFactor;
   leaves p and q unmodified */
poly      p_ExtGcd(poly p, poly &pFactor, poly q, poly &qFactor, ring r);

/* syszygy stuff */
BOOLEAN   p_VectorHasUnitB(poly p, int * k, const ring r);
void      p_VectorHasUnit(poly p, int * k, int * len, const ring r);
poly      p_TakeOutComp1(poly * p, int k, const ring r);
// Splits *p into two polys: *q which consists of all monoms with
// component == comp and *p of all other monoms *lq == pLength(*q)
// On return all components pf *q == 0
void p_TakeOutComp(poly *p, long comp, poly *q, int *lq, const ring r);

// This is something weird -- Don't use it, unless you know what you are doing
poly      p_TakeOutComp(poly * p, int k, const ring r);

void      p_DeleteComp(poly * p,int k, const ring r);

/*-------------ring management:----------------------*/

// resets the pFDeg and pLDeg: if pLDeg is not given, it is
// set to currRing->pLDegOrig, i.e. to the respective LDegProc which
// only uses pFDeg (and not pDeg, or pTotalDegree, etc).
// If you use this, make sure your procs does not make any assumptions
// on ordering and/or OrdIndex -- otherwise they might return wrong results
// on strat->tailRing
void pSetDegProcs(ring r, pFDegProc new_FDeg, pLDegProc new_lDeg = NULL);
// restores pFDeg and pLDeg:
void pRestoreDegProcs(ring r, pFDegProc old_FDeg, pLDegProc old_lDeg);

/*-------------pComp for syzygies:-------------------*/
void p_SetModDeg(intvec *w, ring r);

/*------------ Jet ----------------------------------*/
poly pp_Jet(poly p, int m, const ring R);
poly p_Jet(poly p, int m,const ring R);
poly pp_JetW(poly p, int m, short *w, const ring R);
poly p_JetW(poly p, int m, short *w, const ring R);

poly n_PermNumber(const number z, const int *par_perm, const int OldPar, const ring src, const ring dst);

poly p_PermPoly (poly p, const int * perm,const ring OldRing, const ring dst,
                     nMapFunc nMap, const int *par_perm=NULL, int OldPar=0);

/*----------------------------------------------------*/
poly p_Series(int n,poly p,poly u, intvec *w, const ring R);
poly p_Invers(int n,poly u,intvec *w, const ring R);



/*----------------------------------------------------*/
int   p_Var(poly mi, const ring r);
/// the minimal index of used variables - 1
int   p_LowVar (poly p, const ring r);

/*----------------------------------------------------*/
// returns the length of a polynomial (numbers of monomials) and the last mon.
// respect syzComp
poly p_Last(poly a, int &l, const ring r);

/// shifts components of the vector p by i
void p_Shift (poly * p,int i, const ring r);
#endif // P_POLYS_H