Changeset bf183f in git


Ignore:
Timestamp:
Nov 10, 2010, 12:09:08 PM (14 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'c1c65551e4b3a0de0b13ddf30446f6e9569681b5')
Children:
8d1d30c2bf4983fc7e94265f6fac94f09b3385fb
Parents:
041dfe18affe0cfc2ab3f798958a8e596f6bad6f
git-author:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2010-11-10 12:09:08+01:00
git-committer:
Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 11:55:37+01:00
Message:
degree stuff and p_Power, p_Sub
Location:
polys
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • polys/monomials/p_polys.cc

    r041dfe1 rbf183f  
    413413
    414414/* comptible with ordering */
    415 long pDeg(poly a, const ring r)
     415long p_Deg(poly a, const ring r)
    416416{
    417417  p_LmCheckPolyRing(a, r);
     
    516516}
    517517
    518 int pWeight(int i, const ring r)
     518int p_Weight(int i, const ring r)
    519519{
    520520  if ((r->firstwv==NULL) || (i>r->firstBlockEnds))
     
    525525}
    526526
    527 long pWDegree(poly p, const ring r)
     527long p_WDegree(poly p, const ring r)
    528528{
    529529  if (r->firstwv==NULL) return p_Totaldegree(p, r);
     
    14581458  return result;
    14591459}
     1460/*2
     1461* subtract p2 from p1, p1 and p2 are destroyed
     1462* do not put attention on speed: the procedure is only used in the interpreter
     1463*/
     1464poly p_Sub(poly p1, poly p2, const ring r)
     1465{
     1466  return p_Add_q(p1, p_Neg(p2,r),r);
     1467}
     1468
     1469/*3
     1470* compute for a monomial m
     1471* the power m^exp, exp > 1
     1472* destroys p
     1473*/
     1474static poly p_MonPower(poly p, int exp, const ring r)
     1475{
     1476  int i;
     1477
     1478  if(!n_IsOne(pGetCoeff(p),r))
     1479  {
     1480    number x, y;
     1481    y = pGetCoeff(p);
     1482    n_Power(y,exp,&x,r);
     1483    n_Delete(&y,r);
     1484    pSetCoeff0(p,x);
     1485  }
     1486  for (i=rVar(r); i!=0; i--)
     1487  {
     1488    p_MultExp(p,i, exp,r);
     1489  }
     1490  p_Setm(p,r);
     1491  return p;
     1492}
     1493
     1494/*3
     1495* compute for monomials p*q
     1496* destroys p, keeps q
     1497*/
     1498static void p_MonMult(poly p, poly q, const ring r)
     1499{
     1500  number x, y;
     1501  int i;
     1502
     1503  y = pGetCoeff(p);
     1504  x = n_Mult(y,pGetCoeff(q),r);
     1505  n_Delete(&y,r);
     1506  pSetCoeff0(p,x);
     1507  //for (i=pVariables; i!=0; i--)
     1508  //{
     1509  //  pAddExp(p,i, pGetExp(q,i));
     1510  //}
     1511  //p->Order += q->Order;
     1512  p_ExpVectorAdd(p,q,r);
     1513}
     1514
     1515/*3
     1516* compute for monomials p*q
     1517* keeps p, q
     1518*/
     1519static poly p_MonMultC(poly p, poly q, const ring rr)
     1520{
     1521  number x;
     1522  int i;
     1523  poly r = p_Init(rr);
     1524
     1525  x = n_Mult(pGetCoeff(p),pGetCoeff(q),rr);
     1526  pSetCoeff0(r,x);
     1527  p_ExpVectorSum(r,p, q, rr);
     1528  return r;
     1529}
     1530
     1531/*
     1532*  compute for a poly p = head+tail, tail is monomial
     1533*          (head + tail)^exp, exp > 1
     1534*          with binomial coef.
     1535*/
     1536static poly p_TwoMonPower(poly p, int exp, const ring r)
     1537{
     1538  int eh, e;
     1539  long al;
     1540  poly *a;
     1541  poly tail, b, res, h;
     1542  number x;
     1543  number *bin = pnBin(exp);
     1544
     1545  tail = pNext(p);
     1546  if (bin == NULL)
     1547  {
     1548    p_MonPower(p,exp,r);
     1549    p_MonPower(tail,exp,r);
     1550#ifdef PDEBUG
     1551    p_Test(p,r);
     1552#endif
     1553    return p;
     1554  }
     1555  eh = exp >> 1;
     1556  al = (exp + 1) * sizeof(poly);
     1557  a = (poly *)omAlloc(al);
     1558  a[1] = p;
     1559  for (e=1; e<exp; e++)
     1560  {
     1561    a[e+1] = p_MonMultC(a[e],p,r);
     1562  }
     1563  res = a[exp];
     1564  b = p_Head(tail,r);
     1565  for (e=exp-1; e>eh; e--)
     1566  {
     1567    h = a[e];
     1568    x = n_Mult(bin[exp-e],pGetCoeff(h),r);
     1569    p_SetCoeff(h,x,r);
     1570    p_MonMult(h,b,r);
     1571    res = pNext(res) = h;
     1572    p_MonMult(b,tail,r);
     1573  }
     1574  for (e=eh; e!=0; e--)
     1575  {
     1576    h = a[e];
     1577    x = n_Mult(bin[e],pGetCoeff(h),r);
     1578    p_SetCoeff(h,x,r);
     1579    p_MonMult(h,b,r);
     1580    res = pNext(res) = h;
     1581    p_MonMult(b,tail,r);
     1582  }
     1583  p_LmDelete(&tail,r);
     1584  pNext(res) = b;
     1585  pNext(b) = NULL;
     1586  res = a[exp];
     1587  omFreeSize((ADDRESS)a, al);
     1588  pnFreeBin(bin, exp);
     1589//  tail=res;
     1590// while((tail!=NULL)&&(pNext(tail)!=NULL))
     1591// {
     1592//   if(nIsZero(pGetCoeff(pNext(tail))))
     1593//   {
     1594//     pLmDelete(&pNext(tail));
     1595//   }
     1596//   else
     1597//     pIter(tail);
     1598// }
     1599#ifdef PDEBUG
     1600  p_Test(res,r);
     1601#endif
     1602  return res;
     1603}
     1604
     1605static poly p_Pow(poly p, int i, const ring r)
     1606{
     1607  poly rc = p_Copy(p,r);
     1608  i -= 2;
     1609  do
     1610  {
     1611    rc = p_Mult_q(rc,p_Copy(p,r),r);
     1612    p_Normalize(rc,r);
     1613    i--;
     1614  }
     1615  while (i != 0);
     1616  return p_Mult_q(rc,p,r);
     1617}
     1618
     1619/*2
     1620* returns the i-th power of p
     1621* p will be destroyed
     1622*/
     1623poly p_Power(poly p, int i, const ring r)
     1624{
     1625  poly rc=NULL;
     1626
     1627  if (i==0)
     1628  {
     1629    p_Delete(&p,r);
     1630    return p_One(r);
     1631  }
     1632
     1633  if(p!=NULL)
     1634  {
     1635    if ( (i > 0) && ((unsigned long ) i > (r->bitmask)))
     1636    {
     1637      Werror("exponent %d is too large, max. is %ld",i,r->bitmask);
     1638      return NULL;
     1639    }
     1640    switch (i)
     1641    {
     1642// cannot happen, see above
     1643//      case 0:
     1644//      {
     1645//        rc=pOne();
     1646//        pDelete(&p);
     1647//        break;
     1648//      }
     1649      case 1:
     1650        rc=p;
     1651        break;
     1652      case 2:
     1653        rc=p_Mult_q(p_Copy(p,r),p,r);
     1654        break;
     1655      default:
     1656        if (i < 0)
     1657        {
     1658          p_Delete(&p,r);
     1659          return NULL;
     1660        }
     1661        else
     1662        {
     1663#ifdef HAVE_PLURAL
     1664          if (rIsPluralRing(r)) /* in the NC case nothing helps :-( */
     1665          {
     1666            int j=i;
     1667            rc = p_Copy(p,r);
     1668            while (j>1)
     1669            {
     1670              rc = p_Mult_q(p_Copy(p,r),rc,r);
     1671              j--;
     1672            }
     1673            p_Delete(&p,r);
     1674            return rc;
     1675          }
     1676#endif
     1677          rc = pNext(p);
     1678          if (rc == NULL)
     1679            return p_MonPower(p,i,r);
     1680          /* else: binom ?*/
     1681          int char_p=rChar(r);
     1682          if ((pNext(rc) != NULL)
     1683#ifdef HAVE_RINGS
     1684             || rField_is_Ring(r)
     1685#endif
     1686             )
     1687            return p_Pow(p,i,r);
     1688          if ((char_p==0) || (i<=char_p))
     1689            return p_TwoMonPower(p,i,r);
     1690          poly p_p=p_TwoMonPower(p_Copy(p,r),char_p,r);
     1691          return p_Mult_q(p_Power(p,i-char_p,r),p_p,r);
     1692        }
     1693      /*end default:*/
     1694    }
     1695  }
     1696  return rc;
     1697}
    14601698/***************************************************************
    14611699 *
  • polys/monomials/p_polys.h

    r041dfe1 rbf183f  
    249249static inline poly p_Mult_mm(poly p, poly m, const ring r);
    250250
    251 // returns p+q, destroys p and q
     251/// returns p+q, destroys p and q
    252252static inline poly p_Add_q(poly p, poly q, const ring r);
    253 // like p_Add_q, except that if lp == pLength(lp) lq == pLength(lq)
    254 // then lp == pLength(p+q)
     253/// like p_Add_q, except that if lp == pLength(lp) lq == pLength(lq) then lp == pLength(p+q)
    255254static inline poly p_Add_q(poly p, poly q, int &lp, int lq, const ring r);
     255
     256poly      p_Sub(poly a, poly b, const ring r);
    256257
    257258// return p - m*q, destroys p; const: q,m
     
    290291static inline poly p_SortAdd(poly p, const ring r, BOOLEAN revert = FALSE);
    291292
     293poly      p_Power(poly p, int i, const ring r);
    292294/***************************************************************
    293295 *
     
    345347extern pLDegProc pLDeg;
    346348extern pFDegProc pFDeg;
    347 int  pWeight(int i, ring r);
    348 long pDeg(poly p, ring r);
    349349long p_WFirstTotalDegree(poly p, ring r);
    350350long p_WTotaldegree(poly p, const ring r);
    351 long pWDegree(poly p, ring r);
     351long p_WDegree(poly p,const ring r);
    352352long pLDeg0(poly p,int *l, ring r);
    353353long pLDeg0c(poly p,int *l, ring r);
     
    17601760poly      p_Diff(poly a, int k, const ring r);
    17611761poly      p_DiffOp(poly a, poly b,BOOLEAN multiply, const ring r);
     1762int       p_Weight(int c, const ring r);
     1763static inline long pTotaldegree(poly p) { return p_Totaldegree(p,currRing); }
     1764
    17621765#endif // P_POLYS_H
    17631766
  • polys/polys.h

    r041dfe1 rbf183f  
    258258extern pLDegProc pLDeg;
    259259extern pFDegProc pFDeg;
    260 int  pWeight(int c, const ring r = currRing);
    261 long pDeg(poly p,const ring r = currRing);
     260#define pWeight(c) p_Weight(c,currRing)
     261#define pDeg(p)    p_Deg(p,currRing)
    262262static inline long pTotaldegree(poly p) { return p_Totaldegree(p,currRing); }
    263 long pWTotaldegree(poly p,const ring r = currRing);
    264 long pWDegree(poly p,const ring r = currRing);
    265 long pLDeg0(poly p,int *l,const ring r = currRing);
    266 long pLDeg0c(poly p,int *l,const ring r = currRing);
    267 long pLDegb(poly p,int *l,const ring r = currRing);
    268 long pLDeg1(poly p,int *l,const ring r = currRing);
    269 long pLDeg1c(poly p,int *l,const ring r = currRing);
    270 long pLDeg1_Deg(poly p,int *l,const ring r = currRing);
    271 long pLDeg1c_Deg(poly p,int *l,const ring r = currRing);
    272 long pLDeg1_Totaldegree(poly p,int *l,const ring r = currRing);
    273 long pLDeg1c_Totaldegree(poly p,int *l,const ring r = currRing);
    274 long pLDeg1_WFirstTotalDegree(poly p,int *l,const ring r=currRing);
    275 long pLDeg1c_WFirstTotalDegree(poly p,int *l,const ring r=currRing);
     263#define pWTotaldegree(p) p_WTotaldegree(p,currRing)
     264#define pWDegree(poly p) p_WDegree(p,currRing)
    276265
    277266/*-------------pComp for syzygies:-------------------*/
     
    283272
    284273/*-------------operations on polynomials:------------*/
    285 poly      pSub(poly a, poly b);
    286 poly      p_Power(poly p, int i, const ring r);
     274#define   pSub(a,b) p_Sub(a,b,currRing)
     275
    287276#define pmInit(a,b) p_mInit(a,b,currRing)
    288277
  • polys/polys1.cc

    r041dfe1 rbf183f  
    7171
    7272
    73 /*2
    74 * subtract p2 from p1, p1 and p2 are destroyed
    75 * do not put attention on speed: the procedure is only used in the interpreter
    76 */
    77 poly pSub(poly p1, poly p2)
    78 {
    79   return pAdd(p1, pNeg(p2));
    80 }
    81 
    8273/*3
    8374*  create binomial coef.
     
    124115  }
    125116  omFreeSize((ADDRESS)bin, h*sizeof(number));
    126 }
    127 
    128 /*3
    129 * compute for a monomial m
    130 * the power m^exp, exp > 1
    131 * destroys p
    132 */
    133 static poly p_MonPower(poly p, int exp, const ring r)
    134 {
    135   int i;
    136 
    137   if(!n_IsOne(pGetCoeff(p),r))
    138   {
    139     number x, y;
    140     y = pGetCoeff(p);
    141     n_Power(y,exp,&x,r);
    142     n_Delete(&y,r);
    143     pSetCoeff0(p,x);
    144   }
    145   for (i=rVar(r); i!=0; i--)
    146   {
    147     p_MultExp(p,i, exp,r);
    148   }
    149   p_Setm(p,r);
    150   return p;
    151 }
    152 
    153 /*3
    154 * compute for monomials p*q
    155 * destroys p, keeps q
    156 */
    157 static void p_MonMult(poly p, poly q, const ring r)
    158 {
    159   number x, y;
    160 
    161   y = pGetCoeff(p);
    162   x = n_Mult(y,pGetCoeff(q),r);
    163   n_Delete(&y,r);
    164   pSetCoeff0(p,x);
    165   //for (int i=pVariables; i!=0; i--)
    166   //{
    167   //  pAddExp(p,i, pGetExp(q,i));
    168   //}
    169   //p->Order += q->Order;
    170   p_ExpVectorAdd(p,q,r);
    171 }
    172 
    173 /*3
    174 * compute for monomials p*q
    175 * keeps p, q
    176 */
    177 static poly p_MonMultC(poly p, poly q, const ring rr)
    178 {
    179   number x;
    180   poly r = p_Init(rr);
    181 
    182   x = n_Mult(pGetCoeff(p),pGetCoeff(q),rr);
    183   pSetCoeff0(r,x);
    184   p_ExpVectorSum(r,p, q, rr);
    185   return r;
    186 }
    187 
    188 /*
    189 *  compute for a poly p = head+tail, tail is monomial
    190 *          (head + tail)^exp, exp > 1
    191 *          with binomial coef.
    192 */
    193 static poly p_TwoMonPower(poly p, int exp, const ring r)
    194 {
    195   int eh, e;
    196   long al;
    197   poly *a;
    198   poly tail, b, res, h;
    199   number x;
    200   number *bin = pnBin(exp);
    201 
    202   tail = pNext(p);
    203   if (bin == NULL)
    204   {
    205     p_MonPower(p,exp,r);
    206     p_MonPower(tail,exp,r);
    207 #ifdef PDEBUG
    208     p_Test(p,r);
    209 #endif
    210     return p;
    211   }
    212   eh = exp >> 1;
    213   al = (exp + 1) * sizeof(poly);
    214   a = (poly *)omAlloc(al);
    215   a[1] = p;
    216   for (e=1; e<exp; e++)
    217   {
    218     a[e+1] = p_MonMultC(a[e],p,r);
    219   }
    220   res = a[exp];
    221   b = p_Head(tail,r);
    222   for (e=exp-1; e>eh; e--)
    223   {
    224     h = a[e];
    225     x = n_Mult(bin[exp-e],pGetCoeff(h),r);
    226     p_SetCoeff(h,x,r);
    227     p_MonMult(h,b,r);
    228     res = pNext(res) = h;
    229     p_MonMult(b,tail,r);
    230   }
    231   for (e=eh; e!=0; e--)
    232   {
    233     h = a[e];
    234     x = n_Mult(bin[e],pGetCoeff(h),r);
    235     p_SetCoeff(h,x,r);
    236     p_MonMult(h,b,r);
    237     res = pNext(res) = h;
    238     p_MonMult(b,tail,r);
    239   }
    240   p_LmDelete(&tail,r);
    241   pNext(res) = b;
    242   pNext(b) = NULL;
    243   res = a[exp];
    244   omFreeSize((ADDRESS)a, al);
    245   pnFreeBin(bin, exp);
    246 //  tail=res;
    247 // while((tail!=NULL)&&(pNext(tail)!=NULL))
    248 // {
    249 //   if(nIsZero(pGetCoeff(pNext(tail))))
    250 //   {
    251 //     pLmDelete(&pNext(tail));
    252 //   }
    253 //   else
    254 //     pIter(tail);
    255 // }
    256 #ifdef PDEBUG
    257   p_Test(res,r);
    258 #endif
    259   return res;
    260 }
    261 
    262 static poly p_Pow(poly p, int i, const ring r)
    263 {
    264   poly rc = p_Copy(p,r);
    265   i -= 2;
    266   do
    267   {
    268     rc = p_Mult_q(rc,p_Copy(p,r),r);
    269     p_Normalize(rc,r);
    270     i--;
    271   }
    272   while (i != 0);
    273   return p_Mult_q(rc,p,r);
    274 }
    275 
    276 /*2
    277 * returns the i-th power of p
    278 * p will be destroyed
    279 */
    280 poly p_Power(poly p, int i, const ring r)
    281 {
    282   poly rc=NULL;
    283 
    284   if (i==0)
    285   {
    286     p_Delete(&p,r);
    287     return p_One(r);
    288   }
    289 
    290   if(p!=NULL)
    291   {
    292     if ( (i > 0) && ((unsigned long ) i > (r->bitmask)))
    293     {
    294       Werror("exponent %d is too large, max. is %ld",i,r->bitmask);
    295       return NULL;
    296     }
    297     switch (i)
    298     {
    299 // cannot happen, see above
    300 //      case 0:
    301 //      {
    302 //        rc=pOne();
    303 //        pDelete(&p);
    304 //        break;
    305 //      }
    306       case 1:
    307         rc=p;
    308         break;
    309       case 2:
    310         rc=p_Mult_q(p_Copy(p,r),p,r);
    311         break;
    312       default:
    313         if (i < 0)
    314         {
    315           p_Delete(&p,r);
    316           return NULL;
    317         }
    318         else
    319         {
    320 #ifdef HAVE_PLURAL
    321           if (rIsPluralRing(r)) /* in the NC case nothing helps :-( */
    322           {
    323             int j=i;
    324             rc = p_Copy(p,r);
    325             while (j>1)
    326             {
    327               rc = p_Mult_q(p_Copy(p,r),rc,r);
    328               j--;
    329             }
    330             p_Delete(&p,r);
    331             return rc;
    332           }
    333 #endif
    334           rc = pNext(p);
    335           if (rc == NULL)
    336             return p_MonPower(p,i,r);
    337           /* else: binom ?*/
    338           int char_p=rChar(r);
    339           if ((pNext(rc) != NULL)
    340 #ifdef HAVE_RINGS
    341              || rField_is_Ring(r)
    342 #endif
    343              )
    344             return p_Pow(p,i,r);
    345           if ((char_p==0) || (i<=char_p))
    346             return p_TwoMonPower(p,i,r);
    347           poly p_p=p_TwoMonPower(p_Copy(p,r),char_p,r);
    348           return p_Mult_q(p_Power(p,i-char_p,r),p_p,r);
    349         }
    350       /*end default:*/
    351     }
    352   }
    353   return rc;
    354117}
    355118
Note: See TracChangeset for help on using the changeset viewer.