source: git/Singular/p_polys.cc @ fdca1c0

/****************************************
*  Computer Algebra System SINGULAR     *
****************************************/
/***************************************************************
 *  File:    p_polys.cc
 *  Purpose: implementation of currRing independent poly procedures
 *  Author:  obachman (Olaf Bachmann)
 *  Created: 8/00
 *  Version: $Id: p_polys.cc,v 1.8 2000-11-09 16:32:53 obachman Exp $
 *******************************************************************/

#include "mod2.h"
#include "structs.h"
#include "tok.h"
#include "p_polys.h"
#include "ring.h"

/***************************************************************
 *
 * Completing what needs to be set for the monomial
 *
 ***************************************************************/
// this is special for the syz stuff
static int* _Components = NULL;
static long* _ShiftedComponents = NULL;
static int _ExternalComponents = 0;

void p_Setm_General(poly p, ring r)
{
  p_LmCheckPolyRing(p, r);
  int pos=0;
  if (r->typ!=NULL)
  {
    loop
    {
      long ord=0;
      sro_ord* o=&(r->typ[pos]);
      switch(o->ord_typ)
      {
        case ro_dp:
        {
          int a,e;
          a=o->data.dp.start;
          e=o->data.dp.end;
          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r);
          p->exp[o->data.dp.place]=ord;
          break;
        }
        case ro_wp_neg:
          ord=POLY_NEGWEIGHT_OFFSET;
          // no break;
        case ro_wp:
        {
          int a,e;
          a=o->data.wp.start;
          e=o->data.wp.end;
          int *w=o->data.wp.weights;
          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r)*w[i-a];
          p->exp[o->data.wp.place]=ord;
          break;
        }
        case ro_cp:
        {
          int a,e;
          a=o->data.cp.start;
          e=o->data.cp.end;
          int pl=o->data.cp.place;
          for(int i=a;i<=e;i++) { p->exp[pl]=p_GetExp(p,i,r); pl++; }
          break;
        }
        case ro_syzcomp:
        {
          int c=p_GetComp(p,r);
          long sc = c;
          int* Components = (_ExternalComponents ? _Components : 
                             o->data.syzcomp.Components);
          long* ShiftedComponents = (_ExternalComponents ? _ShiftedComponents: 
                                     o->data.syzcomp.ShiftedComponents);
          if (ShiftedComponents != NULL)
          {
            assume(Components != NULL);
            assume(c == 0 || Components[c] != 0);
            sc = ShiftedComponents[Components[c]];
            assume(c == 0 || sc != 0);
          }
          p->exp[o->data.syzcomp.place]=sc;
          break;
        }
        case ro_syz:
        {
          int c=p_GetComp(p, r);
          if (c > o->data.syz.limit)
            p->exp[o->data.syz.place] = o->data.syz.curr_index;
          else if (c > 0)
            p->exp[o->data.syz.place]= o->data.syz.syz_index[c];
          else
          {
            assume(c == 0);
            p->exp[o->data.syz.place]= 0;
          }
          break;
        }
        default:
          dReportError("wrong ord in rSetm:%d\n",o->ord_typ);
          return;
      }
      pos++;
      if (pos == r->OrdSize) return;
    }
  }
}

void p_Setm_Syz(poly p, ring r, int* Components, long* ShiftedComponents)
{
  _Components = Components;
  _ShiftedComponents = ShiftedComponents;
  _ExternalComponents = 1;
  p_Setm_General(p, r);
  _ExternalComponents = 0;
}

// dummy for lp, ls, etc
void p_Setm_Dummy(poly p, ring r)
{
  p_LmCheckPolyRing(p, r);
}

// for dp, Dp, ds, etc
void p_Setm_TotalDegree(poly p, ring r)
{
  p_LmCheckPolyRing(p, r);
  p->exp[r->pOrdIndex] = p_ExpVectorQuerSum(p, r);
}

// for wp, Wp, ws, etc
void p_Setm_WFirstTotalDegree(poly p, ring r)
{
  p_LmCheckPolyRing(p, r);
  p->exp[r->pOrdIndex] = pWFirstTotalDegree(p, r);
}

p_SetmProc p_GetSetmProc(ring r)
{
  // covers lp, rp, ls, 
  if (r->typ == NULL) return p_Setm_Dummy;

  if (r->OrdSize == 1)
  {
    if (r->typ[0].ord_typ == ro_dp && 
        r->typ[0].data.dp.start == 1 &&
        r->typ[0].data.dp.end == r->N &&
        r->typ[0].data.dp.place == r->pOrdIndex)
      return p_Setm_TotalDegree;
    if (r->typ[0].ord_typ == ro_wp && 
        r->typ[0].data.wp.start == 1 &&
        r->typ[0].data.wp.end == r->N &&
        r->typ[0].data.wp.place == r->pOrdIndex &&
        r->typ[0].data.wp.weights == r->firstwv)
      return p_Setm_WFirstTotalDegree;
  }
  return p_Setm_General;
}


/* -------------------------------------------------------------------*/
/* several possibilities for pFDeg: the degree of the head term       */
/*2
* compute the degree of the leading monomial of p
* the ordering is compatible with degree, use a->order
*/
long pDeg(poly a, ring r)
{
  p_LmCheckPolyRing(a, r);
  assume(p_GetOrder(a, r) == pWTotaldegree(a, r));
  return p_GetOrder(a, r);
}

/*2
* compute the degree of the leading monomial of p
* with respect to weigths 1
* (all are 1 so save multiplications or they are of different signs)
* the ordering is not compatible with degree so do not use p->Order
*/
long pTotaldegree(poly p, ring r)
{
  p_LmCheckPolyRing(p, r);
  return (long) p_ExpVectorQuerSum(p, r);
}


// pWTotalDegree for weighted orderings which cover all variables
// whose first block covers all variables
long pWFirstTotalDegree(poly p, ring r)
{
  int i;
  long sum = 0;
  
  for (i=1; i<= r->firstBlockEnds; i++)
  {
    sum += p_GetExp(p, i, r)*r->firstwv[i-1];
  }
  return sum;
}


/*2
* compute the degree of the leading monomial of p
* with respect to weigths from the ordering
* the ordering is not compatible with degree so do not use p->Order
*/
long pWTotaldegree(poly p, ring r)
{
  p_LmCheckPolyRing(p, r);
  int i, k;
  long j =0;

  // iterate through each block:
  for (i=0;r->order[i]!=0;i++)
  {
    switch(r->order[i])
    {
      case ringorder_wp:
      case ringorder_ws:
      case ringorder_Wp:
      case ringorder_Ws:
        for (k=r->block0[i];k<=r->block1[i];k++)
        { // in jedem block:
          j+= p_GetExp(p,k,r)*r->wvhdl[i][k - r->block0[i]];
        }
        break;
      case ringorder_M:
      case ringorder_lp:
      case ringorder_dp:
      case ringorder_ds:
      case ringorder_Dp:
      case ringorder_Ds:
        for (k=r->block0[i];k<=r->block1[i];k++)
        {
          j+= p_GetExp(p,k,r);
        }
        break;
      case ringorder_c:
      case ringorder_C:
      case ringorder_S:
      case ringorder_s:
        break;
      case ringorder_a:
        for (k=r->block0[i];k<=r->block1[i];k++)
        { // only one line
          j+= p_GetExp(p,k, r)*r->wvhdl[i][ k- r->block0[i]];
        }
        return j;
    }
  }
  return  j;
}

int pWeight(int i, ring r)
{
  if ((r->firstwv==NULL) || (i>r->firstBlockEnds))
  {
    return 1;
  }
  return r->firstwv[i-1];
}

long pWDegree(poly p, ring r)
{
  if (r->firstwv==NULL) return pTotaldegree(p, r);
  p_LmCheckPolyRing(p, r);
  int i, k;
  long j =0;

  for(i=1;i<=r->firstBlockEnds;i++)
    j+=p_GetExp(p, i, r)*r->firstwv[i-1];

  for (;i<=r->N;i++)
    j+=p_GetExp(p,i, r)*pWeight(i, r);

  return j;
}


/* ---------------------------------------------------------------------*/
/* several possibilities for pLDeg: the maximal degree of a monomial in p*/
/*  compute in l also the pLength of p                                   */

/*2
* compute the length of a polynomial (in l)
* and the degree of the monomial with maximal degree: the last one
*/
long pLDeg0(poly p,int *l, ring r)
{
  p_CheckPolyRing(p, r);
  Exponent_t k= p_GetComp(p, r);
  int ll=1;

  while ((pNext(p)!=NULL) && (p_GetComp(pNext(p), r)==k))
  {
    pIter(p);
    ll++;
  }
  *l=ll;
  return p_GetOrder(p, r);
}

/*2
* compute the length of a polynomial (in l)
* and the degree of the monomial with maximal degree: the last one
* but search in all components before syzcomp
*/
#if 0
long pLDeg0c(poly p,int *l, ring r)
{
  p_CheckPolyRing(p, r);
  long o=pFDeg(p, r);
  int ll=1;

  if (! rIsSyzIndexRing(r))
  {
    while ((p=pNext(p))!=NULL)
    {
      o=pFDeg(p, r);
      ll++;
    }
  }
  else
  {
    int curr_limit = rGetCurrSyzLimit(r);
    while ((p=pNext(p))!=NULL)
    {
      if (p_GetComp(p, r)<=curr_limit/*syzComp*/)
      {
        o=pFDeg(p, r);
        ll++;
      }
      else break;
    }
  }
  *l=ll;
  return o;
}
#else
long pLDeg0c(poly p,int *l, ring r)
{
  p_CheckPolyRing(p, r);
  long o;
  int ll=1;

  if (! rIsSyzIndexRing(r))
  {
    while (pNext(p) != NULL) 
    {
      pIter(p);
      ll++;
    }
    o = pFDeg(p, r);
  }
  else
  {
    int curr_limit = rGetCurrSyzLimit(r);
    poly pp = p;
    while ((p=pNext(p))!=NULL)
    {
      if (p_GetComp(p, r)<=curr_limit/*syzComp*/)
        ll++;
      else break;
      pp = p;
    }
    o = pFDeg(pp, r);
  }
  *l=ll;
  return o;
}
#endif

/*2
* compute the length of a polynomial (in l)
* and the degree of the monomial with maximal degree: the first one
* this works for the polynomial case with degree orderings
* (both c,dp and dp,c)
*/
long pLDegb(poly p,int *l, ring r)
{
  p_CheckPolyRing(p, r);
  Exponent_t k= p_GetComp(p, r);
  long o = pFDeg(p, r);
  int ll=1;

  while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
  {
    ll++;
  }
  *l=ll;
  return o;
}

/*2
* compute the length of a polynomial (in l)
* and the degree of the monomial with maximal degree:
* this is NOT the last one, we have to look for it
*/
long pLDeg1(poly p,int *l, ring r)
{
  p_CheckPolyRing(p, r);
  Exponent_t k= p_GetComp(p, r);
  int ll=1;
  long  t,max;

  max=pFDeg(p);
  while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
  {
    t=pFDeg(p, r);
    if (t>max) max=t;
    ll++;
  }
  *l=ll;
  return max;
}

/*2
* compute the length of a polynomial (in l)
* and the degree of the monomial with maximal degree:
* this is NOT the last one, we have to look for it
* in all components
*/
long pLDeg1c(poly p,int *l, ring r)
{
  p_CheckPolyRing(p, r);
  int ll=1;
  long  t,max;

  max=pFDeg(p, r);
  while ((p=pNext(p))!=NULL)
  {
    if (! rIsSyzIndexRing(r) ||
        (p_GetComp(p, r)<=rGetCurrSyzLimit(r)))
    {
       if ((t=pFDeg(p, r))>max) max=t;
       ll++;
    }
    else break;
  }
  *l=ll;
  return max;
}

/***************************************************************
 *
 * Maximal Exponent business
 *
 ***************************************************************/

static inline unsigned long 
p_GetMaxExpL2(unsigned long l1, unsigned long l2, ring r, 
              unsigned long number_of_exp)
{
  const unsigned long bitmask = r->bitmask;
  unsigned long ml1 = l1 & bitmask;
  unsigned long ml2 = l2 & bitmask;
  unsigned long max = (ml1 > ml2 ? ml1 : ml2);
  unsigned long j = number_of_exp - 1;

  if (j > 0)
  {
    unsigned long mask = bitmask << r->BitsPerExp;
    while (1)
    {
      ml1 = l1 & mask;
      ml2 = l2 & mask;
      max |= ((ml1 > ml2 ? ml1 : ml2) & mask);
      j--;
      if (j == 0) break;
      mask = mask << r->BitsPerExp;
    }
  }
  return max;
}

static inline unsigned long
p_GetMaxExpL2(unsigned long l1, unsigned long l2, ring r)
{
  return p_GetMaxExpL2(l1, l2, r, r->ExpPerLong);
}

poly p_GetMaxExpP(poly p, ring r)
{
  p_CheckPolyRing(p, r);
  if (p == NULL) return p_Init(r);
  poly max = p_LmInit(p, r);
  pIter(p);
  if (p == NULL) return max;
  int i, offset;
  unsigned long l_p, l_max;
  unsigned long divmask = r->divmask;
  
  do
  {
    offset = r->VarL_Offset[0];
    l_p = p->exp[offset];
    l_max = max->exp[offset];
    // do the divisibility trick to find out whether l has an exponent
    if (l_p > l_max || 
        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
      max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);

    for (i=1; i<r->VarL_Size; i++)
    {
      offset = r->VarL_Offset[i];
      l_p = p->exp[offset];
      l_max = max->exp[offset];
      // do the divisibility trick to find out whether l has an exponent
      if (l_p > l_max || 
          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
        max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
    }
    pIter(p);
  }
  while (p != NULL);
  return max;
}

unsigned long p_GetMaxExpL(poly p, ring r, unsigned long l_max)
{
  unsigned long l_p, divmask = r->divmask;
  int i;
  
  while (p != NULL)
  {
    l_p = p->exp[r->VarL_Offset[0]];
    if (l_p > l_max ||
        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
      l_max = p_GetMaxExpL2(l_max, l_p, r);
    for (i=1; i<r->VarL_Size; i++)
    {
      l_p = p->exp[r->VarL_Offset[i]];
      // do the divisibility trick to find out whether l has an exponent
      if (l_p > l_max || 
          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
        l_max = p_GetMaxExpL2(l_max, l_p, r);
    }
    pIter(p);
  }
  return l_max;
}