Changeset 2fce0e in git


Ignore:
Timestamp:
Sep 9, 2011, 4:02:38 PM (12 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
e48172b11af08005e2d025ff5a132809d88b6ea5
Parents:
9d5ba2ee7540ee4eb709cfef410629f44c00e75b
git-author:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2011-09-09 16:02:38+02:00
git-committer:
Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 16:12:39+01:00
Message:
FIX: exposed mp_Wedge
FIX: moving duplicates from linalg_from_matpol.cc to matpol.cc
FIX: commented out unused(?!) stuff from linalg_from_matpol.cc in matpol.h
TODO: match&clean up!
Location:
libpolys/polys
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • libpolys/polys/linalg_from_matpol.cc

    r9d5ba2 r2fce0e  
    167167  }
    168168}
    169 
    170 
    171 /*
    172 /// entries of a are minors and go to result (only if not in R)
    173 void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
    174                      ideal R, const ring R)
    175 {
    176   poly *q1;
    177   int e=IDELEMS(result);
    178   int i,j;
    179 
    180   if (R != NULL)
    181   {
    182     for (i=r-1;i>=0;i--)
    183     {
    184       q1 = &(a->m)[i*a->ncols];
    185       for (j=c-1;j>=0;j--)
    186       {
    187         if (q1[j]!=NULL) q1[j] = kNF(R,currQuotient,q1[j]);
    188       }
    189     }
    190   }
    191   for (i=r-1;i>=0;i--)
    192   {
    193     q1 = &(a->m)[i*a->ncols];
    194     for (j=c-1;j>=0;j--)
    195     {
    196       if (q1[j]!=NULL)
    197       {
    198         if (elems>=e)
    199         {
    200           if(e<SIZE_OF_SYSTEM_PAGE)
    201           {
    202             pEnlargeSet(&(result->m),e,e);
    203             e += e;
    204           }
    205           else
    206           {
    207             pEnlargeSet(&(result->m),e,SIZE_OF_SYSTEM_PAGE);
    208             e += SIZE_OF_SYSTEM_PAGE;
    209           }
    210           IDELEMS(result) =e;
    211         }
    212         result->m[elems] = q1[j];
    213         q1[j] = NULL;
    214         elems++;
    215       }
    216     }
    217   }
    218 }
    219 
    220 /// produces recursively the ideal of all arxar-minors of a
    221 void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
    222               poly barDiv, ideal R, const ring R)
    223 {
    224   int k;
    225   int kr=lr-1,kc=lc-1;
    226   matrix nextLevel=mpNew(kr,kc);
    227 
    228   loop
    229   {
    230 // --- look for an optimal row and bring it to last position ------------
    231     if(mpPrepareRow(a,lr,lc)==0) break;
    232 // --- now take all pivots from the last row ------------
    233     k = lc;
    234     loop
    235     {
    236       if(mpPreparePiv(a,lr,k)==0) break;
    237       mpElimBar(a,nextLevel,barDiv,lr,k);
    238       k--;
    239       if (ar>1)
    240       {
    241         mpRecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R);
    242         mpPartClean(nextLevel,kr,k);
    243       }
    244       else mpMinorToResult(result,elems,nextLevel,kr,k,R);
    245       if (ar>k-1) break;
    246     }
    247     if (ar>=kr) break;
    248 // --- now we have to take out the last row...------------
    249     lr = kr;
    250     kr--;
    251   }
    252   mpFinalClean(nextLevel);
    253 }
    254 */
    255 
    256 
    257169
    258170/*
    … …  
    322234/*2
    323235*returns the determinant of the matrix m;
    324 *uses Bareiss algorithm
    325 */
    326 poly mp_DetBareiss (matrix a, const ring R)
    327 {
    328   int s;
    329   poly div, res;
    330   if (MATROWS(a) != MATCOLS(a))
    331   {
    332     Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
    333     return NULL;
    334   }
    335   matrix c = mp_Copy(a, R);
    336   mp_permmatrix *Bareiss = new mp_permmatrix(c, R);
    337   row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
    338 
    339   /* Bareiss */
    340   div = NULL;
    341   while(Bareiss->mpPivotBareiss(&w))
    342   {
    343     Bareiss->mpElimBareiss(div);
    344     div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
    345   }
    346   Bareiss->mpRowReorder();
    347   Bareiss->mpColReorder();
    348   Bareiss->mpSaveArray();
    349   s = Bareiss->mpGetSign();
    350   delete Bareiss;
    351 
    352   /* result */
    353   res = MATELEM(c,1,1);
    354   MATELEM(c,1,1) = NULL;
    355   id_Delete((ideal *)&c, R);
    356   if (s < 0)
    357     res = p_Neg(res, R);
    358   return res;
    359 }
    360 
    361 
    362 
    363 /*2
    364 *returns the determinant of the matrix m;
    365236*uses Newtons formulea for symmetric functions
    366237*/
    … …  
    447318}
    448319
    449 /*2
    450 * compute all ar-minors of the matrix a
    451 */
    452 matrix mp_Wedge(matrix a, int ar, const ring R)
    453 {
    454   int     i,j,k,l;
    455   int *rowchoise,*colchoise;
    456   BOOLEAN rowch,colch;
    457   matrix result;
    458   matrix tmp;
    459   poly p;
    460 
    461   i = binom(a->nrows,ar);
    462   j = binom(a->ncols,ar);
    463 
    464   rowchoise=(int *)omAlloc(ar*sizeof(int));
    465   colchoise=(int *)omAlloc(ar*sizeof(int));
    466   result =mpNew(i,j);
    467   tmp=mpNew(ar,ar);
    468   l = 1; /* k,l:the index in result*/
    469   idInitChoise(ar,1,a->nrows,&rowch,rowchoise);
    470   while (!rowch)
    471   {
    472     k=1;
    473     idInitChoise(ar,1,a->ncols,&colch,colchoise);
    474     while (!colch)
    475     {
    476       for (i=1; i<=ar; i++)
    477       {
    478         for (j=1; j<=ar; j++)
    479         {
    480           MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
    481         }
    482       }
    483       p = mp_DetBareiss(tmp, R);
    484       if ((k+l) & 1) p=p_Neg(p, R);
    485       MATELEM(result,l,k) = p;
    486       k++;
    487       idGetNextChoise(ar,a->ncols,&colch,colchoise);
    488     }
    489     idGetNextChoise(ar,a->nrows,&rowch,rowchoise);
    490     l++;
    491   }
    492 
    493   /*delete the matrix tmp*/
    494   for (i=1; i<=ar; i++)
    495   {
    496     for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
    497   }
    498   id_Delete((ideal *) &tmp, R);
    499   return (result);
    500 }
    501320
    502321///*2

  • libpolys/polys/matpol.cc

    r9d5ba2 r2fce0e  
    14541454}
    14551455
    1456 /*2
    1457 * entries of a are minors and go to result (only if not in R)
    1458 */
     1456/*2*/
     1457/// entries of a are minors and go to result (only if not in R)
    14591458void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
    14601459                     ideal R, const ring ri)
    … …  
    14951494  }
    14961495}
     1496/*
     1497// from  linalg_from_matpol.cc: TODO: compare with above & remove...
     1498void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
     1499                     ideal R, const ring R)
     1500{
     1501  poly *q1;
     1502  int e=IDELEMS(result);
     1503  int i,j;
     1504
     1505  if (R != NULL)
     1506  {
     1507    for (i=r-1;i>=0;i--)
     1508    {
     1509      q1 = &(a->m)[i*a->ncols];
     1510      for (j=c-1;j>=0;j--)
     1511      {
     1512        if (q1[j]!=NULL) q1[j] = kNF(R,currQuotient,q1[j]);
     1513      }
     1514    }
     1515  }
     1516  for (i=r-1;i>=0;i--)
     1517  {
     1518    q1 = &(a->m)[i*a->ncols];
     1519    for (j=c-1;j>=0;j--)
     1520    {
     1521      if (q1[j]!=NULL)
     1522      {
     1523        if (elems>=e)
     1524        {
     1525          if(e<SIZE_OF_SYSTEM_PAGE)
     1526          {
     1527            pEnlargeSet(&(result->m),e,e);
     1528            e += e;
     1529          }
     1530          else
     1531          {
     1532            pEnlargeSet(&(result->m),e,SIZE_OF_SYSTEM_PAGE);
     1533            e += SIZE_OF_SYSTEM_PAGE;
     1534          }
     1535          IDELEMS(result) =e;
     1536        }
     1537        result->m[elems] = q1[j];
     1538        q1[j] = NULL;
     1539        elems++;
     1540      }
     1541    }
     1542  }
     1543}
     1544*/
    14971545
    14981546static void mpFinalClean(matrix a)
    … …  
    15021550}
    15031551
    1504 /*2
    1505 * produces recursively the ideal of all arxar-minors of a
    1506 */
     1552/*2*/
     1553/// produces recursively the ideal of all arxar-minors of a
    15071554void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
    15081555              poly barDiv, ideal R, const ring r)
    … …  
    15381585  mpFinalClean(nextLevel);
    15391586}
    1540 
    1541 /*2
    1542 *returns the determinant of the matrix m;
    1543 *uses Bareiss algorithm
    1544 */
     1587/*
     1588// from  linalg_from_matpol.cc: TODO: compare with above & remove...
     1589void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
     1590              poly barDiv, ideal R, const ring R)
     1591{
     1592  int k;
     1593  int kr=lr-1,kc=lc-1;
     1594  matrix nextLevel=mpNew(kr,kc);
     1595
     1596  loop
     1597  {
     1598// --- look for an optimal row and bring it to last position ------------
     1599    if(mpPrepareRow(a,lr,lc)==0) break;
     1600// --- now take all pivots from the last row ------------
     1601    k = lc;
     1602    loop
     1603    {
     1604      if(mpPreparePiv(a,lr,k)==0) break;
     1605      mpElimBar(a,nextLevel,barDiv,lr,k);
     1606      k--;
     1607      if (ar>1)
     1608      {
     1609        mpRecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R);
     1610        mpPartClean(nextLevel,kr,k);
     1611      }
     1612      else mpMinorToResult(result,elems,nextLevel,kr,k,R);
     1613      if (ar>k-1) break;
     1614    }
     1615    if (ar>=kr) break;
     1616// --- now we have to take out the last row...------------
     1617    lr = kr;
     1618    kr--;
     1619  }
     1620  mpFinalClean(nextLevel);
     1621}
     1622*/
     1623
     1624/*2*/
     1625/// returns the determinant of the matrix m;
     1626/// uses Bareiss algorithm
    15451627poly mp_DetBareiss (matrix a, const ring r)
    15461628{
    … …  
    15771659  return res;
    15781660}
     1661/*
     1662// from  linalg_from_matpol.cc: TODO: compare with above & remove...
     1663poly mp_DetBareiss (matrix a, const ring R)
     1664{
     1665  int s;
     1666  poly div, res;
     1667  if (MATROWS(a) != MATCOLS(a))
     1668  {
     1669    Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
     1670    return NULL;
     1671  }
     1672  matrix c = mp_Copy(a, R);
     1673  mp_permmatrix *Bareiss = new mp_permmatrix(c, R);
     1674  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
     1675
     1676  // Bareiss
     1677  div = NULL;
     1678  while(Bareiss->mpPivotBareiss(&w))
     1679  {
     1680    Bareiss->mpElimBareiss(div);
     1681    div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
     1682  }
     1683  Bareiss->mpRowReorder();
     1684  Bareiss->mpColReorder();
     1685  Bareiss->mpSaveArray();
     1686  s = Bareiss->mpGetSign();
     1687  delete Bareiss;
     1688
     1689  // result
     1690  res = MATELEM(c,1,1);
     1691  MATELEM(c,1,1) = NULL;
     1692  id_Delete((ideal *)&c, R);
     1693  if (s < 0)
     1694    res = p_Neg(res, R);
     1695  return res;
     1696}
     1697*/
     1698
     1699/*2
     1700* compute all ar-minors of the matrix a
     1701*/
     1702matrix mp_Wedge(matrix a, int ar, const ring R)
     1703{
     1704  int     i,j,k,l;
     1705  int *rowchoise,*colchoise;
     1706  BOOLEAN rowch,colch;
     1707  matrix result;
     1708  matrix tmp;
     1709  poly p;
     1710
     1711  i = binom(a->nrows,ar);
     1712  j = binom(a->ncols,ar);
     1713
     1714  rowchoise=(int *)omAlloc(ar*sizeof(int));
     1715  colchoise=(int *)omAlloc(ar*sizeof(int));
     1716  result = mpNew(i,j);
     1717  tmp    = mpNew(ar,ar);
     1718  l = 1; /* k,l:the index in result*/
     1719  idInitChoise(ar,1,a->nrows,&rowch,rowchoise);
     1720  while (!rowch)
     1721  {
     1722    k=1;
     1723    idInitChoise(ar,1,a->ncols,&colch,colchoise);
     1724    while (!colch)
     1725    {
     1726      for (i=1; i<=ar; i++)
     1727      {
     1728        for (j=1; j<=ar; j++)
     1729        {
     1730          MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
     1731        }
     1732      }
     1733      p = mp_DetBareiss(tmp, R);
     1734      if ((k+l) & 1) p=p_Neg(p, R);
     1735      MATELEM(result,l,k) = p;
     1736      k++;
     1737      idGetNextChoise(ar,a->ncols,&colch,colchoise);
     1738    }
     1739    idGetNextChoise(ar,a->nrows,&rowch,rowchoise);
     1740    l++;
     1741  }
     1742
     1743  /*delete the matrix tmp*/
     1744  for (i=1; i<=ar; i++)
     1745  {
     1746    for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
     1747  }
     1748  id_Delete((ideal *) &tmp, R);
     1749  return (result);
     1750}

  • libpolys/polys/matpol.h

    r9d5ba2 r2fce0e  
    5656poly TraceOfProd ( matrix a, matrix b, int n, const ring r);
    5757
    58 poly mp_Det (matrix m, const ring r);
     58// poly mp_Det (matrix m, const ring r);
    5959matrix mp_Wedge(matrix a, int ar, const ring r);
    6060
    6161// BOOLEAN mpJacobi(leftv res,leftv a);
    6262// BOOLEAN mpKoszul(leftv res,leftv b/*in*/, leftv c/*ip*/, leftv id=NULL);
     63
    6364poly mp_DetBareiss (matrix a, const ring r);
    6465
Note: See TracChangeset for help on using the changeset viewer.