Changeset e44b149 in git


Ignore:
Timestamp:
Jul 3, 2013, 12:04:32 PM (11 years ago)
Author:
Yue Ren <ren@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
f80a5305e1e0960855d4bd8702e48d1c2bbf8710
Parents:
a667138296109b29e11fe1654025d4581abfc50a
git-author:
Yue Ren <ren@mathematik.uni-kl.de>2013-07-03 12:04:32+02:00
git-committer:
Yue Ren <ren@mathematik.uni-kl.de>2013-07-18 14:56:50+02:00
Message:
fix: missed assume for allowing empty bigintmats
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/gitfan.lib

    ra667138 re44b149  
    2727
    2828LIB "parallel.lib"; // for parallelWaitAll
     29LIB "general.lib";
    2930
    3031////////////////////////////////////////////////////
    31 proc mod_init()
    32 {
    33   LIB"gfanlib.so";
    34 }
    35 
    36 static proc int2face(int n, int r)
     32
     33proc int2face(int n, int r)
    3734{
    3835  int k = r-1;
    … …  
    4340    while(2^k > n0){
    4441      k--;
    45       //v[size(v)+1] = 0;
    4642    }
    4743
    … …  
    5652/////////////////////////////////
    5753
    58 proc isAface(ideal a, intvec gam0)
     54proc isAface(ideal a, intvec gam0, int n)
    5955"USAGE:  isAface(a,gam0); a: ideal, gam0:intvec
    6056PURPOSE: Checks whether the face of the positive orthant,
    … …  
    6662"
    6763{
    68   poly pz;
    69 
    7064  // special case: gam0 is the zero-cone:
    7165  if (size(gam0) == 1 and gam0[1] == 0){
    72     ideal G;
     66    poly pz; ideal G;
    7367
    7468    // is an a-face if and only if RL0(0,...,0) = const
    … …  
    9387  }
    9488
    95 
    96   // ring is too big: Switch to KK[T_i; e_i\in gam0] instead:
     89  // ring is too big: Switch to KK[T_i | e_i\in gam0] instead:
     90  intvec w=ringlist(basering)[3][1][2];
     91  intvec w0;
    9792  string initNewRing = "ring Rgam0 = 0,(";
    9893
    9994  for (int i=1; i<size(gam0); i++){
    10095    initNewRing = initNewRing + string(var(gam0[i])) + ",";
    101   }
    102 
     96    w0[i]=w[gam0[i]];
     97  }
     98
     99  w0 = w0,w[gam0[i]],1;
     100  initNewRing = initNewRing + string(var(gam0[size(gam0)])) + ",U),Wp("+string(w0)+");";
     101  def R = basering;
     102  execute(initNewRing);
     103
     104  ideal agam0 = imap(R,a);
     105
     106  for (i = 1; i<=size(agam0); i=i+1)
     107  {
     108    if (size(agam0[i]) == 1)
     109    { return(0,0); }
     110  }
     111
     112  poly p = var(1); // first entry of g; p = prod T_i with i element of g
     113  for ( i = 2; i <= nvars(basering); i++ ) {
     114      p = p * var(i);
     115    }
     116  // p is now the product over all T_i, with e_i in gam0
     117
     118  agam0 = agam0, p - 1; // rad-membership
     119  agam0 = timeStd(agam0,5);
     120  // "timestd finished after: "+string(timer-t);
     121  // int timeout = 0;
     122  if (attrib(agam0,"isSB") < 1)
     123  {
     124    kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R;
     125    return(0,1);
     126    // // "timestd failed in "+string(gam0)+", falling back to saturation!";
     127    // setring R; kill Rgam0;
     128
     129    // initNewRing = "ring Rgam0 = 0,(";
     130
     131    // w0 = 0;
     132    // for (i=1; i<size(gam0); i++){
     133    //   initNewRing = initNewRing + string(var(gam0[i])) + ",";
     134    //   w0[i]=w[gam0[i]];
     135    // }
     136
     137    // w0 = w0,w[gam0[i]];
     138    // initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),Wp("+string(w0)+");";
     139    // execute(initNewRing);
     140
     141    // ideal G = imap(R,a);
     142
     143    // timeout = 1;
     144    // int t=rtimer;
     145    // for(int k=nvars(basering); k>0; k--) { G=sat(G,var(k))[1]; }
     146    // t = (rtimer - t) div 1000;
     147    // string(n)+": saturation successful after "+string(t)+" with: "+string(G<>1);
     148  }
     149
     150  // does G contain 1?, i.e. is G = 1?
     151  if(agam0 <> 1) {
     152    kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R;
     153    return(1,0); // true
     154  }
     155
     156  kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R;
     157  return(0,0); // false
     158}
     159example
     160{
     161  echo = 2;
     162
     163  ring R = 0,(T(1..4)),dp;
     164  ideal I = T(1)*T(2)-T(4);
     165
     166  intvec w = 1,4;
     167  intvec v = 1,2,4;
     168
     169  isAface(I,w); // should be 0
     170  "-----------";
     171  isAface(I,v); // should be 1
     172}
     173
     174
     175proc isAfaceNonZero(ideal a, intvec gam0)
     176{
     177  string initNewRing = "ring Rgam0 = 0,(";
     178  for (int i=1; i<size(gam0); i++)
     179    { initNewRing = initNewRing + string(var(gam0[i])) + ","; }
    103180  initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),dp;";
    104181  def R = basering;
    … …  
    107184  ideal agam0 = imap(R,a);
    108185
    109   poly p = var(1); // first entry of g; p = prod T_i with i element of g
    110   for ( i = 2; i <= nvars(basering); i++ ) {
    111       p = p * var(i);
    112     }
    113   // p is now the product over all T_i, with e_i in gam0
    114 
    115   agam0 = agam0, p - 1; // rad-membership
     186  for ( i = 1; i<=size(agam0); i=i+1)
     187    { if (size(agam0[i]) == 1) { return(0); } }
     188
     189  poly p = var(1);
     190  for ( i = 2; i <= nvars(basering); i++ )
     191    { p = p * var(i); }
     192
     193  agam0 = agam0, p - 1;
    116194  ideal G = std(agam0);
    117195
    118   // does G contain 1?, i.e. is G = 1?
    119   if(G <> 1) {
    120     return(1); // true
    121   }
    122 
    123   return(0); // false
    124 }
    125 example
    126 {
    127   echo = 2;
    128 
    129   ring R = 0,(T(1..4)),dp;
    130   ideal I = T(1)*T(2)-T(4);
    131 
    132   intvec w = 1,4;
    133   intvec v = 1,2,4;
    134 
    135   isAface(I,w); // should be 0
    136   "-----------";
    137   isAface(I,v); // should be 1
    138 }
    139 
     196  if(G <> 1)
     197    { return(1); }
     198
     199  return(0);
     200}
    140201////////////////////////////////////////////////////
    141202
    … …  
    145206  list AF;
    146207
    147   for(i = start; i <= end; i++){
     208  for(i = start; i <= end; i=i+1){
    148209    if(i < 2^r){
     210      string(i)+"    "+string(size(AF));
    149211      gam0 = int2face(i,r);
    150212
    … …  
    165227
    166228proc afaces(ideal a, list #)
    167 "USAGE:  afaces(a, b, c); a: ideal, d: int, c: int
     229"USAGE:  afaces(a, b, c); a: ideal, b: int, c: int
    168230PURPOSE: Returns a list of all a-faces (represented by intvecs).
    169231         Moreover, it is possible to specify a dimensional bound b,
    … …  
    210272
    211273  // do remaining ones:
    212   for(i = ncores * step +1; i < 2^r; i++){
     274  for(i = ncores * step +1; i < 2^r; i=i+1){
    213275    "another one needed";
    214276    gam0 = int2face(i,r);
    … …  
    223285  }
    224286
     287  int l;
    225288  // read out afaces of out into AF:
    226   for(i = 1; i <= size(out); i++){
    227     AF = AF + out[i];
     289  for(l = 1; l <= size(out); l++){
     290    AF = AF + out[l];
    228291  }
    229292

  • dyn_modules/callgfanlib/Makefile.am

    ra667138 re44b149  
    77  module_LTLIBRARIES=gfanlib.la
    88
    9 # if WANT_DEBUG
    10 #   module_LTLIBRARIES+=gfanlib_g.la
    11 # endif
    12  
     9if WANT_DEBUG
     10  module_LTLIBRARIES+=gfanlib_g.la
     11endif
     12
    1313endif
    1414endif
    1515
    16 # CXXTEMPLFLAGS = ${PIPE}
     16CXXTEMPLFLAGS   = ${PIPE} -lcddgmp
    1717## -fno-implicit-templates
    1818## --no-exceptions
    1919
    2020# DEBUGCXXFLAGS = -g -Wextra -Wall -pedantic -Wno-long-long ${CXXTEMPLFLAGS}
    21 # 
     21#
    2222# NODEBUGCFLAGS   = -O3 -w -fomit-frame-pointer ${PIPE}
    2323# NODEBUGCXXFLAGS = -O3 -w -fomit-frame-pointer ${CXXTEMPLFLAGS}
    … …  
    4545P_PROCS_CPPFLAGS_COMMON = -DDYNAMIC_VERSION
    4646
    47 gfanlib_la_CPPFLAGS   = ${MYINCLUDES} ${P_PROCS_CPPFLAGS_COMMON} 
     47gfanlib_la_CPPFLAGS   = ${MYINCLUDES} ${P_PROCS_CPPFLAGS_COMMON}
    4848# ${NODEBUGDEFS}
    4949# gfanlib_g_la_CPPFLAGS = ${MYINCLUDES} ${P_PROCS_CPPFLAGS_COMMON}

  • dyn_modules/callgfanlib/bbcone.cc

    ra667138 re44b149  
    111111  bigintmat* bim = zMatrixToBigintmat(zm);
    112112  char* s = bim->StringAsPrinted();
     113  if (s==NULL)
     114    s = (char*) omAlloc0(sizeof(char));
    113115  delete bim;
    114116  return s;
    … …  
    127129    s<<"INEQUALITIES"<<std::endl;
    128130  s<<toString(i)<<std::endl;
     131  // std::cout << toString(i) << std::endl;
    129132  if (c->areImpliedEquationsKnown())
    130133    s<<"LINEAR_SPAN"<<std::endl;
    … …  
    132135    s<<"EQUATIONS"<<std::endl;
    133136  s<<toString(e)<<std::endl;
     137  // std::cout << toString(e) << std::endl;
    134138  return s.str();
    135139}
    … …  
    15591563  return TRUE;
    15601564}
     1565
     1566gfan::ZVector intStar2ZVector(const int d, const int* i)
     1567{
     1568  gfan::ZVector zv(d);
     1569  for(int j=0; j<d; j++)
     1570    zv[j]=i[j+1];
     1571  return zv;
     1572}
     1573
     1574BOOLEAN maximalGroebnerCone(leftv res, leftv args)
     1575{
     1576  leftv u = args;
     1577  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
     1578  {
     1579    leftv v = u->next;
     1580    if (v == NULL)
     1581    {
     1582      int n = currRing->N;
     1583      ideal I = (ideal) u->Data();
     1584      poly g = NULL;
     1585      int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
     1586      int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
     1587      gfan::ZVector leadexpw = gfan::ZVector(n);
     1588      gfan::ZVector tailexpw = gfan::ZVector(n);
     1589      gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
     1590      for (int i=0; i<IDELEMS(I); i++)
     1591      {
     1592        g = (poly) I->m[i]; pGetExpV(g,leadexpv);
     1593        leadexpw = intStar2ZVector(n, leadexpv);
     1594        pIter(g);
     1595        while (g != NULL)
     1596        {
     1597          pGetExpV(g,tailexpv);
     1598          tailexpw = intStar2ZVector(n, tailexpv);
     1599          inequalities.appendRow(leadexpw-tailexpw);
     1600          pIter(g);
     1601        }
     1602      }
     1603      gfan::ZCone* gCone = new gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
     1604      omFreeSize(leadexpv,(n+1)*sizeof(int));
     1605      omFreeSize(tailexpv,(n+1)*sizeof(int));
     1606
     1607      res->rtyp = coneID;
     1608      res->data = (void*) gCone;
     1609      return FALSE;
     1610    }
     1611  }
     1612  WerrorS("maximalGroebnerCone: unexpected parameters");
     1613  return TRUE;
     1614}
     1615
     1616
     1617poly initial(poly p)
     1618{
     1619  poly g = p;
     1620  poly h = p_Head(g,currRing);
     1621  poly f = h;
     1622  long d = p_Deg(g,currRing);
     1623  pIter(g);
     1624  while ((g != NULL) && (p_Deg(g,currRing) == d))
     1625  {
     1626    pNext(h) = p_Head(g,currRing);
     1627    pIter(h);
     1628    pIter(g);
     1629  }
     1630  return(f);
     1631}
     1632
     1633
     1634BOOLEAN initial(leftv res, leftv args)
     1635{
     1636  leftv u = args;
     1637  if ((u != NULL) && (u->Typ() == POLY_CMD))
     1638  {
     1639    leftv v = u->next;
     1640    if (v == NULL)
     1641    {
     1642      poly p = (poly) u->Data();
     1643      res->rtyp = POLY_CMD;
     1644      res->data = (void*) initial(p);
     1645      return FALSE;
     1646    }
     1647  }
     1648  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
     1649  {
     1650    leftv v = u->next;
     1651    if (v == NULL)
     1652    {
     1653      ideal I = (ideal) u->Data();
     1654      ideal inI = idInit(IDELEMS(I));
     1655      poly g; poly h; long d;
     1656      for (int i=0; i<IDELEMS(I); i++)
     1657      {
     1658        g = (poly) I->m[i];
     1659        inI->m[i]=initial(g);
     1660      }
     1661      res->rtyp = IDEAL_CMD;
     1662      res->data = (void*) inI;
     1663      return FALSE;
     1664    }
     1665  }
     1666  WerrorS("initial: unexpected parameters");
     1667  return TRUE;
     1668}
     1669
     1670
     1671BOOLEAN homogeneitySpace(leftv res, leftv args)
     1672{
     1673  leftv u = args;
     1674  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
     1675  {
     1676    leftv v = u->next;
     1677    if (v == NULL)
     1678    {
     1679      int n = currRing->N;
     1680      ideal I = (ideal) u->Data();
     1681      poly g;
     1682      int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
     1683      int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
     1684      gfan::ZVector leadexpw = gfan::ZVector(n);
     1685      gfan::ZVector tailexpw = gfan::ZVector(n);
     1686      gfan::ZMatrix equations = gfan::ZMatrix(0,n);
     1687      for (int i=0; i<IDELEMS(I); i++)
     1688      {
     1689        g = (poly) I->m[i]; pGetExpV(g,leadexpv);
     1690        leadexpw = intStar2ZVector(n, leadexpv);
     1691        pIter(g);
     1692        while (g != NULL)
     1693        {
     1694          pGetExpV(g,tailexpv);
     1695          tailexpw = intStar2ZVector(n, tailexpv);
     1696          equations.appendRow(leadexpw-tailexpw);
     1697          pIter(g);
     1698        }
     1699      }
     1700      gfan::ZCone* gCone = new gfan::ZCone(gfan::ZMatrix(0, equations.getWidth()),equations);
     1701      omFreeSize(leadexpv,(n+1)*sizeof(int));
     1702      omFreeSize(tailexpv,(n+1)*sizeof(int));
     1703
     1704      res->rtyp = coneID;
     1705      res->data = (void*) gCone;
     1706      return FALSE;
     1707    }
     1708  }
     1709  WerrorS("homogeneitySpace: unexpected parameters");
     1710  return TRUE;
     1711}
     1712
     1713
     1714BOOLEAN groebnerCone(leftv res, leftv args)
     1715{
     1716  leftv u = args;
     1717  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
     1718  {
     1719    leftv v = u->next;
     1720    if (v == NULL)
     1721    {
     1722      int n = currRing->N;
     1723      ideal I = (ideal) u->Data();
     1724      poly g = NULL;
     1725      int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
     1726      int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
     1727      gfan::ZVector leadexpw = gfan::ZVector(n);
     1728      gfan::ZVector tailexpw = gfan::ZVector(n);
     1729      gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
     1730      gfan::ZMatrix equations = gfan::ZMatrix(0,n);
     1731      long d;
     1732      for (int i=0; i<IDELEMS(I); i++)
     1733      {
     1734        g = (poly) I->m[i]; pGetExpV(g,leadexpv);
     1735        leadexpw = intStar2ZVector(n, leadexpv);
     1736        pIter(g);
     1737        d = p_Deg(g,currRing);
     1738        while ((g != NULL) && (p_Deg(g,currRing) == d))
     1739        {
     1740          pGetExpV(g,tailexpv);
     1741          tailexpw = intStar2ZVector(n, tailexpv);
     1742          equations.appendRow(leadexpw-tailexpw);
     1743          pIter(g);
     1744        }
     1745
     1746        if (g != NULL)
     1747        {
     1748          while (g != NULL)
     1749          {
     1750            pGetExpV(g,tailexpv);
     1751            tailexpw = intStar2ZVector(n, tailexpv);
     1752            inequalities.appendRow(leadexpw-tailexpw);
     1753            pIter(g);
     1754          }
     1755        }
     1756      }
     1757      gfan::ZCone* gCone = new gfan::ZCone(inequalities,equations);
     1758      omFreeSize(leadexpv,(n+1)*sizeof(int));
     1759      omFreeSize(tailexpv,(n+1)*sizeof(int));
     1760
     1761      res->rtyp = coneID;
     1762      res->data = (void*) gCone;
     1763      return FALSE;
     1764    }
     1765  }
     1766  WerrorS("groebnerCone: unexpected parameters");
     1767  return TRUE;
     1768}
     1769
     1770/***
     1771 * Given a cone and a point in its boundary,
     1772 * returns the inner normal vector of a facet
     1773 * containing the point.
     1774 * In case no facet contains the point,
     1775 * then 0 is returned.
     1776 **/
     1777gfan::ZVector* facetContaining(gfan::ZCone* zc, gfan::ZVector* zv)
     1778{
     1779  gfan::ZMatrix facets = zc->getFacets();
     1780  for (int i=0; i<facets.getHeight(); i++)
     1781  {
     1782    gfan::ZVector facet = facets[i];
     1783    if (dot(facet,*zv) == (long) 0)
     1784      return new gfan::ZVector(facet);
     1785  }
     1786  return new gfan::ZVector(zc->ambientDimension());
     1787}
     1788
     1789
     1790BOOLEAN facetContaining(leftv res, leftv args)
     1791{
     1792  leftv u = args;
     1793  if ((u != NULL) && (u->Typ() == coneID))
     1794  {
     1795    leftv v = u->next;
     1796    if ((v != NULL) && ((v->Typ() == BIGINTMAT_CMD) || (v->Typ() == INTVEC_CMD)))
     1797    {
     1798      gfan::ZCone* zc = (gfan::ZCone*) u->Data();
     1799
     1800      bigintmat* point1;
     1801      if (v->Typ() == INTVEC_CMD)
     1802      {
     1803        intvec* point0 = (intvec*) v->Data();
     1804        point1 = iv2bim(point0,coeffs_BIGINT)->transpose();
     1805      }
     1806      else
     1807        point1 = (bigintmat*) v->Data();
     1808
     1809      gfan::ZVector* point = bigintmatToZVector(*point1);
     1810      gfan::ZVector* facet = facetContaining(zc, point);
     1811
     1812      res->rtyp = BIGINTMAT_CMD;
     1813      res->data = (void*) zVectorToBigintmat(*facet);
     1814
     1815      delete facet;
     1816      delete point;
     1817      if (v->Typ() == INTVEC_CMD)
     1818        delete point1;
     1819      return FALSE;
     1820    }
     1821  }
     1822  WerrorS("facetContaining: unexpected parameters");
     1823  return TRUE;
     1824}
     1825
    15611826
    15621827void bbcone_setup(SModulFunctions* p)
    … …  
    16121877  p->iiAddCproc("","uniquePoint",FALSE,uniquePoint);
    16131878  p->iiAddCproc("","listContainsCone",FALSE,containsCone);
     1879  p->iiAddCproc("","maximalGroebnerCone",FALSE,maximalGroebnerCone);
     1880  p->iiAddCproc("","groebnerCone",FALSE,groebnerCone);
     1881  p->iiAddCproc("","initial",FALSE,initial);
     1882  p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
     1883  p->iiAddCproc("","facetContaining",FALSE,facetContaining);
    16141884  coneID=setBlackboxStuff(b,"cone");
    16151885}

  • dyn_modules/callgfanlib/bbcone.h

    ra667138 re44b149  
    2424gfan::ZVector* bigintmatToZVector(const bigintmat &bim);
    2525
    26 gfan::ZVector* intStar2ZVector(const int d, const int* i);
     26gfan::ZVector intStar2ZVector(const int d, const int* i);
    2727char* toString(gfan::ZMatrix const &m);
    2828std::string toString(const gfan::ZCone* const c);

  • dyn_modules/callgfanlib/bbfan.cc

    ra667138 re44b149  
    279279            if (d-ld>=0)
    280280            {
    281               int n = zf->numberOfConesOfDimension(d,oo,mm);
     281              int n = zf->numberOfConesOfDimension(d-ld,oo,mm);
    282282              res->rtyp = INT_CMD;
    283283              res->data = (void*) (long) n;
    284284              return FALSE;
    285285            }
    286             else
    287             {
    288               res->rtyp = INT_CMD;
    289               res->data = (void*) (long) 0;
    290               return FALSE;
    291             }
     286            res->rtyp = INT_CMD;
     287            res->data = (void*) (long) 0;
     288            return FALSE;
    292289          }
    293290        }
    … …  
    309306
    310307      for (int i=0; i<=d; i++)
    311         { n = n + zf->numberOfConesOfDimension(i,0,0); }
     308        n = n + zf->numberOfConesOfDimension(i,0,0);
    312309
    313310      res->rtyp = INT_CMD;
    … …  
    396393
    397394      leftv w=v->next;
    398       int n = 1;
     395      int n;
    399396      if ((w != NULL) && (w->Typ() == INT_CMD))
    400         int n = (int)(long) w;
     397        n = (int)(long) w;
    401398
    402399      if (n != 0)
    … …  
    416413    }
    417414  }
    418   else
    419   {
    420     WerrorS("insertCone: unexpected parameters");
    421     return TRUE;
    422   }
     415  WerrorS("insertCone: unexpected parameters");
     416  return TRUE;
    423417}
    424418
    … …  
    601595    }
    602596  }
    603   else
    604   {
    605     WerrorS("getCone: unexpected parameters");
    606     return TRUE;
    607   }
     597  WerrorS("getCone: unexpected parameters");
     598  return TRUE;
    608599}
    609600
    … …  
    642633          {
    643634            lists L = (lists)omAllocBin(slists_bin);
    644             int n = zf->numberOfConesOfDimension(d,oo,mm);
     635            int n = zf->numberOfConesOfDimension(d-ld,oo,mm);
    645636            L->Init(n);
    646637            for (int i=0; i<n; i++)
    … …  
    672663    }
    673664  }
    674   else
    675   {
    676     WerrorS("getCones: unexpected parameters");
    677     return TRUE;
    678   }
     665  WerrorS("getCones: unexpected parameters");
     666  return TRUE;
    679667}
    680668
    … …  
    732720gfan::ZMatrix rays(gfan::ZFan* zf)
    733721{
    734   int linDim = zf->getLinealityDimension();
    735722  gfan::ZMatrix rays(0,zf->getAmbientDimension());
    736   for (int i=0; i<zf->numberOfConesOfDimension(linDim+1,0,0); i++)
    737   {
    738     gfan::ZCone zc = zf->getCone(linDim+1, i, 0, 0);
     723  for (int i=0; i<zf->numberOfConesOfDimension(1,0,0); i++)
     724  {
     725    gfan::ZCone zc = zf->getCone(1, i, 0, 0);
    739726    rays.append(zc.extremeRays());
    740727  }
    … …  
    924911  }
    925912  WerrorS("listOfFacets: unexpected parameters");
     913  return TRUE;
     914}
     915
     916BOOLEAN tropicalVariety(leftv res, leftv args)
     917{
     918  leftv u=args;
     919  if ((u != NULL) && (u->Typ() == POLY_CMD))
     920  {
     921    int n = rVar(currRing);
     922    gfan::ZFan* zf = new gfan::ZFan(n);
     923    int* expv1 = (int*)omAlloc((n+1)*sizeof(int));
     924    int* expv2 = (int*)omAlloc((n+1)*sizeof(int));
     925    int* expvr = (int*)omAlloc((n+1)*sizeof(int));
     926    gfan::ZVector expw1 = gfan::ZVector(n);
     927    gfan::ZVector expw2 = gfan::ZVector(n);
     928    gfan::ZVector expwr = gfan::ZVector(n);
     929    gfan::ZMatrix eq, ineq;
     930    for (poly s1=(poly)u->Data(); s1!=NULL; pIter(s1))
     931    {
     932      pGetExpV(s1,expv1);
     933      expw1 = intStar2ZVector(n,expv1);
     934      for (poly s2=pNext(s1); s2!=NULL; pIter(s2))
     935      {
     936        pGetExpV(s2,expv2);
     937        expw2 = intStar2ZVector(n,expv2);
     938        eq = gfan::ZMatrix(0,n);
     939        eq.appendRow(expw1-expw2);
     940        ineq = gfan::ZMatrix(0,n);
     941        for (poly r=(poly)u->Data(); r!=NULL; pIter(r))
     942        {
     943          pGetExpV(r,expvr);
     944          expwr = intStar2ZVector(n,expvr);
     945          if ((r!=s1) && (r!=s2))
     946          {
     947            ineq.appendRow(expw1-expwr);
     948          }
     949        }
     950        gfan::ZCone zc = gfan::ZCone(ineq,eq);
     951        zf->insert(zc);
     952      }
     953    }
     954    omFreeSize(expv1,(n+1)*sizeof(int));
     955    omFreeSize(expv2,(n+1)*sizeof(int));
     956    omFreeSize(expvr,(n+1)*sizeof(int));
     957    res->rtyp = fanID;
     958    res->data = (void*) zf;
     959    return FALSE;
     960  }
     961  WerrorS("tropicalVariety: unexpected parameters");
     962  return TRUE;
     963}
     964
     965gfan::ZFan* commonRefinement(gfan::ZFan* zf, gfan::ZFan* zg)
     966{
     967  assume(zf->getAmbientDimension() == zg->getAmbientDimension());
     968
     969  // gather all maximal cones of f and g
     970  std::list<gfan::ZCone> maximalConesOfF;
     971  for (int d=0; d<=zf->getAmbientDimension(); d++)
     972  {
     973    for (int i=0; i<zf->numberOfConesOfDimension(d,0,1); i++)
     974    {
     975      maximalConesOfF.push_back(zf->getCone(d,i,0,1));
     976    }
     977  }
     978
     979  std::list<gfan::ZCone> maximalConesOfG;
     980  for (int d=0; d<=zg->getAmbientDimension(); d++)
     981  {
     982    for (int i=0; i<zg->numberOfConesOfDimension(d,0,1); i++)
     983    {
     984      maximalConesOfG.push_back(zg->getCone(d,i,0,1));
     985    }
     986  }
     987
     988  // construct a new fan out of their intersections
     989  gfan::ZFan* zr = new gfan::ZFan(zf->getAmbientDimension());
     990  for (std::list<gfan::ZCone>::iterator itf=maximalConesOfF.begin();
     991       itf != maximalConesOfF.end(); itf++)
     992  {
     993    for (std::list<gfan::ZCone>::iterator itg=maximalConesOfG.begin();
     994         itg != maximalConesOfG.end(); itg++)
     995    {
     996      zr->insert(intersection(*itf,*itg));
     997    }
     998  }
     999
     1000  return zr;
     1001}
     1002
     1003BOOLEAN commonRefinement(leftv res, leftv args)
     1004{
     1005  leftv u=args;
     1006  if ((u != NULL) && (u->Typ() == fanID))
     1007  {
     1008    leftv v=u->next;
     1009    if ((v != NULL) && (v->Typ() == fanID))
     1010    {
     1011      gfan::ZFan* zf = (gfan::ZFan*) u->Data();
     1012      gfan::ZFan* zg = (gfan::ZFan*) v->Data();
     1013      gfan::ZFan* zr = commonRefinement(zf,zg);
     1014      res->rtyp = fanID;
     1015      res->data = (void*) zr;
     1016      return FALSE;
     1017    }
     1018  }
     1019  WerrorS("commonRefinement: unexpected parameters");
    9261020  return TRUE;
    9271021}
    … …  
    10121106  p->iiAddCproc("","fVector",FALSE,fVector);
    10131107  p->iiAddCproc("","containsInCollection",FALSE,containsInCollection);
     1108  p->iiAddCproc("","tropicalVariety",FALSE,tropicalVariety);
     1109  p->iiAddCproc("","commonRefinement",FALSE,commonRefinement);
    10141110  // iiAddCproc("","grFan",FALSE,grFan);
    10151111  fanID=setBlackboxStuff(b,"fan");

  • dyn_modules/callgfanlib/gitfan.cc

    ra667138 re44b149  
    164164  {
    165165    leftv v=u->next;
    166     if ((u != NULL) && (v->Typ() == BIGINTMAT_CMD))
     166    if ((v != NULL) && (v->Typ() == BIGINTMAT_CMD))
    167167    {
    168168      lists cones = (lists) u->Data();
    … …  
    257257}
    258258
     259
     260static int binomial(int n, int k)
     261{
     262  if (n<k)
     263    return(0);
     264  gfan::Integer num = 1;
     265  gfan::Integer den = 1;
     266  for (int i=1; i<=k; i++)
     267    den = den*i;
     268  for (int j=n-k+1; j<=n; j++)
     269    num = num*j;
     270  gfan::Integer bin = num/den;
     271  return(bin.toInt());
     272}
     273
     274
     275intvec* intToAface(unsigned int v0, int n, int k)
     276{
     277  intvec* v = new intvec(k);
     278  int j = 0;
     279  for (int i=0; i<n; i++)
     280  {
     281    if (v0 & (1<<i))
     282      (*v)[j++] = i+1;
     283  }
     284  return v;
     285}
     286
     287
     288BOOLEAN listOfAfacesToCheck(leftv res, leftv args)
     289{
     290  leftv u = args;
     291  if ((u != NULL) && (u->Typ() == INT_CMD))
     292  {
     293    leftv v = u->next;
     294    if ((v != NULL) && (v->Typ() == INT_CMD))
     295    {
     296      int n = (int)(long) u->Data();
     297      int k = (int)(long) v->Data();
     298      unsigned int v = 0;
     299      for (int i=0; i<k; i++)
     300        v |= 1<<i;  // sets the first k bits of v as 1
     301
     302      lists L = (lists)omAllocBin(slists_bin);
     303      int count = (int) binomial(n,k); L->Init(count);
     304      unsigned int t;
     305      while (!(v & (1<<n)))
     306      {
     307        L->m[--count].rtyp = INTVEC_CMD;
     308        L->m[count].data = (void*) intToAface(v,n,k);
     309
     310        // t gets v's least significant 0 bits set to 1
     311        t = v | (v - 1);
     312        // Next set to 1 the most significant bit to change,
     313        // set to 0 the least significant ones, and add the necessary 1 bits.
     314        v = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1));
     315      }
     316      res->rtyp = LIST_CMD;
     317      res->data = (void*) L;
     318      return FALSE;
     319    }
     320  }
     321  WerrorS("listOfAfacesToCheck: unexpected parameter");
     322  return TRUE;
     323}
     324
     325
     326BOOLEAN nextAfaceToCheck(leftv res, leftv args)
     327{
     328  leftv u = args;
     329  if ((u != NULL) && (u->Typ() == INTVEC_CMD))
     330  {
     331    leftv v = u->next;
     332    if ((v != NULL) && (v->Typ() == INT_CMD))
     333    {
     334      leftv w = v->next;
     335      if ((w != NULL) && (w->Typ() == INT_CMD))
     336      {
     337        intvec* aface = (intvec*) u->Data();
     338        int ambientDimension = (int)(long) v->Data();
     339        int dimension = (int)(long) w->Data();
     340
     341        unsigned int af = 0;
     342        for (int i=0; i<aface->length(); i++)
     343          af |= 1<<((*aface)[i]-1);
     344
     345        unsigned int t = af | (af - 1);
     346        af = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(af) + 1));
     347
     348        if (af & (1<<ambientDimension))
     349        {
     350          res->rtyp = INTVEC_CMD;
     351          res->data = (void*) new intvec(1);
     352          return FALSE;
     353        }
     354
     355        res->rtyp = INTVEC_CMD;
     356        res->data = (void*) intToAface(af,ambientDimension,dimension);
     357        return FALSE;
     358      }
     359    }
     360  }
     361  WerrorS("nextAfaceToCheck: unexpected parameter");
     362  return TRUE;
     363}
     364
     365
    259366void gitfan_setup(SModulFunctions* p)
    260367{
    261368  p->iiAddCproc("","refineCones",FALSE,refineCones);
    262 }
     369  p->iiAddCproc("","listOfAfacesToCheck",FALSE,listOfAfacesToCheck);
     370  p->iiAddCproc("","nextAfaceToCheck",FALSE,nextAfaceToCheck);
     371}

  • libpolys/coeffs/bigintmat.h

    ra667138 re44b149  
    3030    bigintmat(int r, int c, const coeffs n): m_coeffs(n), v(NULL), row(r), col(c)
    3131    {
    32       assume (rows() > 0);
    33       assume (cols() > 0);
     32      assume (rows() >= 0);
     33      assume (cols() >= 0);
    3434
    3535      const int l = r*c;
    … …  
    5353      if (l > 0)
    5454      {
    55         assume (rows() > 0);
    56         assume (cols() > 0);
     55        assume (rows() >= 0);
     56        assume (cols() >= 0);
    5757
    5858        assume (m->v != NULL);
Note: See TracChangeset for help on using the changeset viewer.