Changeset ff2859d in git


Ignore:
Timestamp:
Oct 26, 2018, 8:21:54 PM (6 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b4f17ed1d25f93d46dbe29e4b499baecc2fd51bb')
Children:
7e94f3d724c603bff225ef75ac25ce4722b500ca
Parents:
5c73a12e2a76e9452bc04071e2ff07ca8dc0912539ccee4a37dbeac2e9321f1ab6ad47ceb1df2855
Message:
Merge branch 'spielwiese' of github.com:levandov/Sources into spielwiese
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/ncHilb.lib

    r39ccee rff2859d  
    11//////////////////////////////////////////////////////////////////////////////
    22version="version nc_hilb.lib 4.1.1.0 Dec_2017 "; // $Id$
    3 category="Noncommutative algebra";
     3category="Noncommutative";
    44info="
    5 LIBRARY:  fpahilb.lib: A library for computing graded and multi-graded Hilbert
    6                       series of general non-commutative algebras (available via Letterplace).
    7                                           It also computes the truncated Hilbert series of the algebras whose Hilbert
    8                       series are possibly not rational or unknown.
     5LIBRARY:  fpahilb.lib: Computation of graded and multi-graded Hilbert series of non-commutative algebras (Letterplace).
    96
    107AUTHOR:   Sharwan K. Tiwari   shrawant@gmail.com
    … …  
    1815          La Scala R., Tiwari Sharwan K.: Multigraded Hilbert Series of noncommutative modules, https://arxiv.org/abs/1705.01083.
    1916
     17KEYWORDS: finitely presented algebra; infinitely presented algebra; graded Hilbert series; multi-graded Hilbert series;
     18                 
    2019PROCEDURES:
    2120          fpahilb(L,d,#); Hilbert series of a non-commutative algebra
    … …  
    3635@* N>2: computation of truncated Hilbert series up to degree N-1, and
    3736@* nonempty string: the details about the orbit and system of equations will be printed.
    38 
    3937Let the orbit be @code{O_I = {T_{w_1}(I),...,T_{w_r}(I)} ($w_i\in W$)}, where we assume that if @code{T_{w_i}(I)=T_{w_i'}(I)}
    4038for some @code{w'_i\in W}, then @code{deg(w_i)\leq deg(w'_i)} follows.
    … …  
    4341where @code{{b_j}} is a basis of I.
    4442Moreover, it also prints the linear system (for the information about adjacency matrix).
    45 
    4643EXAMPLE: example fpahilb; shows an example "
    4744{
    48         if (attrib(@R, "isLetterplaceRing")==0)
     45        if (attrib(basering, "isLetterplaceRing")==0)
    4946        {
    5047      ERROR("Basering should be Letterplace ring");
    … …  
    5855    string odp="";
    5956    int i;
    60     for(i=sz ;i >= 1; i--)
     57    for(i=sz; i >= 1; i--)
    6158    {
    6259      if(typeof(#[i])=="string")
    … …  
    7269    }
    7370    i=1;
    74 //only one optional parameter (for printing the details) is allowed as a string.
    75     while(typeof(#[i])=="int" && i<=sz)
    76     {
    77       if(#[i] == 1 && ig==0)
     71// VL: changing the old "only one optional parameter (for printing the details) is allowed as a string."
     72    while( (typeof(#[i])=="int") && (i<=sz) )
     73    {
     74      if( (#[i] == 1) && (ig==0) )
    7875      {
    7976        ig = 1;
    … …  
    8178      else
    8279      {
    83         if(#[i] == 2 && mgrad==0)
     80        if ( (#[i] == 2) && (mgrad==0) )
    8481        {
    8582          mgrad = 2;
    … …  
    8784        else
    8885        {
    89           if(#[i] > 2 && tdeg==0)
     86          if ( (#[i] > 2) && (tdeg==0) )
    9087          {
    9188            tdeg = #[i];
    … …  
    103100      ERROR("error:only int 1,2, >2, and a string are allowed as optional parameters");
    104101    }
    105     /* new: truncation should be < than degbound */
    106         if (tdeg> lV)
     102    /* new: truncation should be < than degbound/2 */
     103        if (tdeg > lV div 2 + 1)
    107104        {
    108        ERROR("unappropriate degree bound");
     105       ERROR("degree bound is too high");
    109106        }
    110107    ideal J_lm = I;
    … …  
    130127  def R = makeLetterplaceRing(10); setring R;
    131128  ideal I = Y*Z, Y*Z*X, Y*Z*Z*X, Y*Z*Z*Z*X, Y*Z*Z*Z*Z*X, Y*Z*Z*Z*Z*Z*X, Y*Z*Z*Z*Z*Z*Z*X, Y*Z*Z*Z*Z*Z*Z*Z*X;
    132   nchilb(I,10);
     129  fpahilb(I,5);
    133130  kill R;
    134131       
    … …  
    139136  // inspecting J we see that this is a homogeneous Groebner basis
    140137  // which is potentially infinite, i.e. J is not finitely generated
    141   nchilb(I,6,1); //third argument '1' is for non-finitely generated case
    142 
    143   ring r=0,(a,b),dp;
    144   ideal I = a^3, a*b^2, a^2*b;
    145   nchilb(l3,5,2); //third argument '2' is to compute multi-graded Hilbert series
     138  fpahilb(J,5,1,2,"p"); // '1' i for non-finitely generated case, string to print details
     139  //'5' here is to compute the truncated HS up to degree 5.
     140  //'2' is to compute multi-graded Hilbert series
     141  fpahilb(J,3,1,"p");
     142}
     143
     144// overhauled by VL
     145proc rcolon(ideal I, poly W)
     146"USAGE:  rcolon(I,w); ideal I, poly w
     147RETURNS: ideal
     148ASSUME: - basering is a Letterplace ring
     149- W is a monomial
     150- I is a monomial ideal
     151PURPOSE: compute a right colon ideal of I by a monomial w
     152NOTE: Output is the set of generators, which should be added to I
     153EXAMPLE: example rcolon; shows an example"
     154{
     155        int lV =  attrib(R,"isLetterplaceRing"); //nvars(save);
     156        if (lV == 0)
     157        {
     158      ERROR("Basering must be a Letterplace ring");
     159        }
     160    poly wc = leadmonom(W);
     161        ideal II = lead(I);
     162    ideal J = system("rcolon", II, wc, lV);
     163    // VL: printlevel here? before only printing was there, no output
     164        return(J);
     165 }
     166example
     167{
     168"EXAMPLE:"; echo = 2;
     169  ring r=0,(X,Y,Z),dp;
     170  def R = makeLetterplaceRing(10); setring R;
     171  ideal I = Y*Z, Y*Z*X, Y*Z*Z*X, Y*Z*Z*Z*X, Y*Z*Z*Z*Z*X, Y*Z*Z*Z*Z*Z*X, Y*Z*Z*Z*Z*Z*Z*X, Y*Z*Z*Z*Z*Z*Z*Z*X;
     172  poly w = Y;
     173  ideal J = rcolon(I,w);
     174  J; // new generators, which need to be added to I
     175}
     176
     177
     178/*
     179proc nchilb(list L_wp, int d, list #)
     180"USAGE: fpahilb(I, d[, L]), ideal I, int d, optional list L
     181PURPOSE: compute Hilbert series of a non-commutative algebra
     182ASSUME:
     183NOTE: d is an integer for the degree bound (maximal total degree of
     184         polynomials of the generating set of the input ideal),
     185#[]=1, computation for non-finitely generated regular ideals,
     186#[]=2, computation of multi-graded Hilbert series,
     187#[]=tdeg, for obtaining the truncated Hilbert series up to the total degree tdeg-1 (tdeg should be > 2), and
     188#[]=string(p), to print the details about the orbit and system of equations.
     189Let the orbit is O_I = {T_{w_1}(I),...,T_{w_r}(I)} ($w_i\in W$), where we assume that if T_{w_i}(I)=T_{w_i'}(I)$
     190for some $w'_i\in W$, then $deg(w_i)\leq deg(w'_i)$.
     191Then, it prints words description of orbit: w_1,...,w_r.
     192It also prints the maximal degree and the cardinality of \sum_j R(w_i, b_j) corresponding to each w_i,
     193where {b_j} is a basis of I.
     194Moreover, it also prints the linear system (for the information about adjacency matrix) and its solving time.
     195
     196NOTE  : A Groebner basis of two-sided ideal of the input should be given in a
     197        special form. This form is a list of modules, where each generator
     198        of every module represents a monomial times a coefficient in the free
     199        associative algebra. The first entry, in each generator, represents a
     200        coefficient and every next entry is a variable.
     201
     202        Ex: module p1=[1,y,z],[-1,z,y], represents the poly y*z-z*y;
     203            module p2=[1,x,z,x],[-1,z,x,z], represents the poly x*z*x-z*x*z
     204        for more details about the input, see examples.
     205EXAMPLE: example nchilb; shows an example "
     206{
     207    if (d<1)
     208    {
     209      ERROR("bad degree bound");
     210    }
     211
     212    def save = basering;
     213    int sz=size(#);
     214    int lV=nvars(save);
     215    int ig=0;
     216    int mgrad=0;
     217    int tdeg=0;
     218    string odp="";
     219    int i;
     220    for(i=sz ;i >= 1; i--)
     221    {
     222      if(typeof(#[i])=="string")
     223      {
     224        if(#[i]!="")
     225        {
     226          odp = "p";
     227        }
     228        # = delete(#,i);
     229        sz = sz-1;
     230        break;
     231      }
     232    }
     233    i=1;
     234//only one optional parameter (for printing the details) is allowed as a string.
     235    while(typeof(#[i])=="int" && i<=sz)
     236    {
     237      if(#[i] == 1 && ig==0)
     238      {
     239        ig = 1;
     240      }
     241      else
     242      {
     243        if(#[i] == 2 && mgrad==0)
     244        {
     245          mgrad = 2;
     246        }
     247        else
     248        {
     249          if(#[i] > 2 && tdeg==0)
     250          {
     251            tdeg = #[i];
     252          }
     253          else
     254          {
     255            ERROR("error: only int 1,2 and >2 are allowed as
     256            optional parameters");
     257          }
     258        }
     259      }
     260      i = i + 1;
     261    }
     262    if( i <= sz)
     263    {
     264      ERROR("error:only int 1,2, >2, and a string are allowed as
     265       optional parameters");
     266    }
     267    if(tdeg==0)
     268    {def R = makeLetterplaceRing(2*d);}
     269    else
     270    {def R = makeLetterplaceRing(2*(tdeg-1));}
     271    setring R;
     272    ideal I;
     273    poly p;
     274    poly q=0;
     275    // convert list L_wp of free-poly to letterPlace-poly format
     276    setring save;
     277    module M;
     278    int j,k,sw,sm,slm;
     279    vector w;
     280    poly pc=0;
     281    intvec v;
     282    slm = size(L_wp);              // number of polys in the given ideal
     283    for (i=1; i<=slm; i++)
     284    {
     285        M  = L_wp[i];
     286        sm = ncols(M);            // number of words in the free-poly M
     287        for (j=1; j<=sm; j++)
     288        {
     289            w  = M[j];
     290            sw = size(w);
     291            for (k=2; k<=sw; k++)
     292            {
     293              v[k-1]=rvar(w[k]);
     294            }
     295            pc=w[1];
     296            setring R;
     297            p=imap(save,pc);
     298            for (k=2; k<=sw; k++)
     299            {
     300              p=p*var(v[k-1]+(k-2)*lV);
     301            }
     302            q=q+p;
     303            setring save;
     304        }
     305        setring R;
     306        I = I,q; //lp-polynomial added to I
     307        q=0;   //ready for the next polynomial
     308        setring save;
     309    }
     310    setring R;
     311    I=simplify(I,2);
     312    ideal J_lm;
     313    for(i=1;i<=size(I);i++)
     314    {
     315        J_lm[i]=leadmonom(I[i]);
     316    }
     317    //compute the Hilbert series
     318    if(odp == "")
     319    {system("nc_hilb", J_lm, lV, ig, mgrad,tdeg);}
     320    else
     321    {system("nc_hilb", J_lm, lV, ig, mgrad,tdeg, odp);}
     322}
     323
     324example
     325{
     326"EXAMPLE:"; echo = 2;
     327
     328    ring r=0,(X,Y,Z),dp;
     329    module p1 =[1,Y,Z];             //represents the poly Y*Z
     330    module p2 =[1,Y,Z,X];          //represents the poly Y*Z*X
     331    module p3 =[1,Y,Z,Z,X,Z];
     332    module p4 =[1,Y,Z,Z,Z,X,Z];
     333    module p5 =[1,Y,Z,Z,Z,Z,X,Z];
     334    module p6 =[1,Y,Z,Z,Z,Z,Z,X,Z];
     335    module p7 =[1,Y,Z,Z,Z,Z,Z,Z,X,Z];
     336    module p8 =[1,Y,Z,Z,Z,Z,Z,Z,Z,X,Z];
     337    list l1=list(p1,p2,p3,p4,p5,p6,p7,p8);
     338    nchilb(l1,10);
     339
     340    ring r=0,(x,y,z),dp;
     341
     342    module p1=[1,y,z],[-1,z,y];               //y*z-z*y
     343    module p2=[1,x,z,x],[-1,z,x,z];           // x*z*x-z*x*z
     344    module p3=[1,x,z,z,x,z],[-1,z,x,z,z,x];   // x*z^2*x*z-z*x*z^2*x
     345    module p4=[1,x,z,z,z,x,z];[-1,z,x,z,z,x,x]; // x*z^3*x*z-z*x*z^2*x^2
     346    list l2=list(p1,p2,p3,p4);
     347
     348    nchilb(l2,6,1); //third argument '1' is for non-finitely generated case
     349
     350    ring r=0,(a,b),dp;
     351    module p1=[1,a,a,a];
     352    module p2=[1,a,b,b];
     353    module p3=[1,a,a,b];
     354
     355    list l3=list(p1,p2,p3);
     356    nchilb(l3,5,2);//third argument '2' is to compute multi-graded HS
    146357
    147358    ring r=0,(x,y,z),dp;
    … …  
    187398    //'11' is to compute the truncated HS up to degree 10.
    188399}
    189 
    190 // overhauled by VL
    191 proc rcolon(ideal I, poly W)
    192 "USAGE:  rcolon(I,w,d); ideal I, poly w
    193 RETURNS: ideal
    194 ASSUME: - basering is a Letterplace ring
    195 - W is a monomial
    196 - I is a monomial ideal
    197 PURPOSE: compute a right colon ideal of I by a monomial w
    198 NOTE: Output is the set of generators, which should be added to I
    199 EXAMPLE: example rcolon; shows an example"
    200 {
    201     if (d<1)
    202     {
    203        ERROR("bad degree bound");
    204     }
    205         int lV =  attrib(R,"isLetterplaceRing"); //nvars(save);
    206         if (lV == 0)
    207         {
    208       ERROR("Basering must be a Letterplace ring");
    209         }
    210     poly wc = leadmonom(W);
    211         ideal II = lead(I);
    212     ideal J = system("rcolon", II, wc, lV);
    213     // VL: printlevel here? before only printing was there, no output
    214         return(J);
    215  }
    216 example
    217 {
    218 "EXAMPLE:"; echo = 2;
    219   ring r=0,(X,Y,Z),dp;
    220   poly w = Y;
    221   ideal I = Y*Z, Y*Z*X, Y*Z*Z*X, Y*Z*Z*Z*X, Y*Z*Z*Z*Z*X, Y*Z*Z*Z*Z*Z*X, Y*Z*Z*Z*Z*Z*Z*X, Y*Z*Z*Z*Z*Z*Z*Z*X;
    222   ideal J = rcolon(I,w,10);
    223   J; // new generators, which need to be added to I
    224 }
    225 
    226 
    227 /*
    228 proc nchilb(list L_wp, int d, list #)
    229 "USAGE: fpahilb(I, d[, L]), ideal I, int d, optional list L
    230 PURPOSE: compute Hilbert series of a non-commutative algebra
    231 ASSUME:
    232 NOTE: d is an integer for the degree bound (maximal total degree of
    233          polynomials of the generating set of the input ideal),
    234 #[]=1, computation for non-finitely generated regular ideals,
    235 #[]=2, computation of multi-graded Hilbert series,
    236 #[]=tdeg, for obtaining the truncated Hilbert series up to the total degree tdeg-1 (tdeg should be > 2), and
    237 #[]=string(p), to print the details about the orbit and system of equations.
    238 Let the orbit is O_I = {T_{w_1}(I),...,T_{w_r}(I)} ($w_i\in W$), where we assume that if T_{w_i}(I)=T_{w_i'}(I)$
    239 for some $w'_i\in W$, then $deg(w_i)\leq deg(w'_i)$.
    240 Then, it prints words description of orbit: w_1,...,w_r.
    241 It also prints the maximal degree and the cardinality of \sum_j R(w_i, b_j) corresponding to each w_i,
    242 where {b_j} is a basis of I.
    243 Moreover, it also prints the linear system (for the information about adjacency matrix) and its solving time.
    244 
    245 NOTE  : A Groebner basis of two-sided ideal of the input should be given in a
    246         special form. This form is a list of modules, where each generator
    247         of every module represents a monomial times a coefficient in the free
    248         associative algebra. The first entry, in each generator, represents a
    249         coefficient and every next entry is a variable.
    250 
    251         Ex: module p1=[1,y,z],[-1,z,y], represents the poly y*z-z*y;
    252             module p2=[1,x,z,x],[-1,z,x,z], represents the poly x*z*x-z*x*z
    253         for more details about the input, see examples.
    254 EXAMPLE: example nchilb; shows an example "
    255 {
    256     if (d<1)
    257     {
    258       ERROR("bad degree bound");
    259     }
    260 
    261     def save = basering;
    262     int sz=size(#);
    263     int lV=nvars(save);
    264     int ig=0;
    265     int mgrad=0;
    266     int tdeg=0;
    267     string odp="";
    268     int i;
    269     for(i=sz ;i >= 1; i--)
    270     {
    271       if(typeof(#[i])=="string")
    272       {
    273         if(#[i]!="")
    274         {
    275           odp = "p";
    276         }
    277         # = delete(#,i);
    278         sz = sz-1;
    279         break;
    280       }
    281     }
    282     i=1;
    283 //only one optional parameter (for printing the details) is allowed as a string.
    284     while(typeof(#[i])=="int" && i<=sz)
    285     {
    286       if(#[i] == 1 && ig==0)
    287       {
    288         ig = 1;
    289       }
    290       else
    291       {
    292         if(#[i] == 2 && mgrad==0)
    293         {
    294           mgrad = 2;
    295         }
    296         else
    297         {
    298           if(#[i] > 2 && tdeg==0)
    299           {
    300             tdeg = #[i];
    301           }
    302           else
    303           {
    304             ERROR("error: only int 1,2 and >2 are allowed as
    305             optional parameters");
    306           }
    307         }
    308       }
    309       i = i + 1;
    310     }
    311     if( i <= sz)
    312     {
    313       ERROR("error:only int 1,2, >2, and a string are allowed as
    314        optional parameters");
    315     }
    316     if(tdeg==0)
    317     {def R = makeLetterplaceRing(2*d);}
    318     else
    319     {def R = makeLetterplaceRing(2*(tdeg-1));}
    320     setring R;
    321     ideal I;
    322     poly p;
    323     poly q=0;
    324     // convert list L_wp of free-poly to letterPlace-poly format
    325     setring save;
    326     module M;
    327     int j,k,sw,sm,slm;
    328     vector w;
    329     poly pc=0;
    330     intvec v;
    331     slm = size(L_wp);              // number of polys in the given ideal
    332     for (i=1; i<=slm; i++)
    333     {
    334         M  = L_wp[i];
    335         sm = ncols(M);            // number of words in the free-poly M
    336         for (j=1; j<=sm; j++)
    337         {
    338             w  = M[j];
    339             sw = size(w);
    340             for (k=2; k<=sw; k++)
    341             {
    342               v[k-1]=rvar(w[k]);
    343             }
    344             pc=w[1];
    345             setring R;
    346             p=imap(save,pc);
    347             for (k=2; k<=sw; k++)
    348             {
    349               p=p*var(v[k-1]+(k-2)*lV);
    350             }
    351             q=q+p;
    352             setring save;
    353         }
    354         setring R;
    355         I = I,q; //lp-polynomial added to I
    356         q=0;   //ready for the next polynomial
    357         setring save;
    358     }
    359     setring R;
    360     I=simplify(I,2);
    361     ideal J_lm;
    362     for(i=1;i<=size(I);i++)
    363     {
    364         J_lm[i]=leadmonom(I[i]);
    365     }
    366     //compute the Hilbert series
    367     if(odp == "")
    368     {system("nc_hilb", J_lm, lV, ig, mgrad,tdeg);}
    369     else
    370     {system("nc_hilb", J_lm, lV, ig, mgrad,tdeg, odp);}
    371 }
    372 
    373 example
    374 {
    375 "EXAMPLE:"; echo = 2;
    376 
    377     ring r=0,(X,Y,Z),dp;
    378     module p1 =[1,Y,Z];             //represents the poly Y*Z
    379     module p2 =[1,Y,Z,X];          //represents the poly Y*Z*X
    380     module p3 =[1,Y,Z,Z,X,Z];
    381     module p4 =[1,Y,Z,Z,Z,X,Z];
    382     module p5 =[1,Y,Z,Z,Z,Z,X,Z];
    383     module p6 =[1,Y,Z,Z,Z,Z,Z,X,Z];
    384     module p7 =[1,Y,Z,Z,Z,Z,Z,Z,X,Z];
    385     module p8 =[1,Y,Z,Z,Z,Z,Z,Z,Z,X,Z];
    386     list l1=list(p1,p2,p3,p4,p5,p6,p7,p8);
    387     nchilb(l1,10);
    388 
    389     ring r=0,(x,y,z),dp;
    390 
    391     module p1=[1,y,z],[-1,z,y];               //y*z-z*y
    392     module p2=[1,x,z,x],[-1,z,x,z];           // x*z*x-z*x*z
    393     module p3=[1,x,z,z,x,z],[-1,z,x,z,z,x];   // x*z^2*x*z-z*x*z^2*x
    394     module p4=[1,x,z,z,z,x,z];[-1,z,x,z,z,x,x]; // x*z^3*x*z-z*x*z^2*x^2
    395     list l2=list(p1,p2,p3,p4);
    396 
    397     nchilb(l2,6,1); //third argument '1' is for non-finitely generated case
    398 
    399     ring r=0,(a,b),dp;
    400     module p1=[1,a,a,a];
    401     module p2=[1,a,b,b];
    402     module p3=[1,a,a,b];
    403 
    404     list l3=list(p1,p2,p3);
    405     nchilb(l3,5,2);//third argument '2' is to compute multi-graded HS
    406 
    407     ring r=0,(x,y,z),dp;
    408     module p1=[1,x,z,y,z,x,z];
    409     module p2=[1,x,z,x];
    410     module p3=[1,x,z,y,z,z,x,z];
    411     module p4=[1,y,z];
    412     module p5=[1,x,z,z,x,z];
    413 
    414     list l4=list(p1,p2,p3,p4,p5);
    415     nchilb(l4,7,"p"); //third argument "p" is to print the details
    416                      // of the orbit and system
    417 
    418     ring r=0,(x,y,z),dp;
    419 
    420     module p1=[1,y,z,z];
    421     module p2=[1,y,y,z];
    422     module p3=[1,x,z,z];
    423     module p4=[1,x,z,y];
    424     module p5=[1,x,y,z];
    425     module p6=[1,x,y,y];
    426     module p7=[1,x,x,z];
    427     module p8=[1,x,x,y];
    428     module p9=[1,y,z,y,z];
    429     module p10=[1,y,z,x,z];
    430     module p11=[1,y,z,x,y];
    431     module p12=[1,x,z,x,z];
    432     module p13=[1,x,z,x,y];
    433     module p14=[1,x,y,x,z];
    434     module p15=[1,x,y,x,y];
    435     module p16=[1,y,z,y,x,z];
    436     module p17=[1,y,z,y,x,y];
    437     module p18=[1,y,z,y,y,x,z];
    438     module p19=[1,y,z,y,y,x,y];
    439     module p20=[1,y,z,y,y,y,x,z];
    440     module p21=[1,y,z,y,y,y,x,y];
    441 
    442     list l5=list(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,
    443     p14,p15,p16,p17,p18,p19,p20,p21);
    444     nchilb(l5,7,1,2,"p");
    445 
    446     nchilb(l5,7,1,2,11,"p");
    447     //'11' is to compute the truncated HS up to degree 10.
    448 }
    449400*/
    450401

  • Singular/dyn_modules/cohomo/cohomo.cc

    r5c73a12 rff2859d  
    88#include "kernel/mod2.h"
    99
     10#include "omalloc/omalloc.h"
    1011#include "misc/mylimits.h"
    1112#include "libpolys/misc/intvec.h"
    … …  
    2526#include "kernel/linear_algebra/linearAlgebra.h"//for printNumber
    2627#include "kernel/GBEngine/kstd1.h"//for gb
    27 
    28 
    29 // #include <kernel/structs.h>
    30 // #include <kernel/polys.h>
    31 //ADICHANGES:
    3228#include <kernel/ideals.h>
    33 // #include <kernel/GBEngine/kstd1.h>
    34 // #include<gmp.h>
    35 // #include<vector>
    36 
    37 //# define omsai 1
    38 //#if omsai
    39 //#if HAVE_STANLEYREISNER
    4029#if 1
    4130
    … …  
    4736#include <Singular/ipshell.h>
    4837#include <Singular/libsingular.h>
    49 
    50 
    51 
    52 
    53 
    54 
    55 
    56 /******************************************************************************/
    57 
    58 
    59 
    60 
    61 
    62 /***************************print***********************************************/
    63 //print vector
     38#include <time.h>
     39
     40
     41
     42
     43
     44
     45
     46/***************************print(only for debugging)***********************************************/
     47//print vector of integers.
    6448void listprint(std::vector<int> vec)
    6549{
    … …  
    7862
    7963
    80 //print vector of vectors
     64//print vector of vectors of integers.
    8165void listsprint(std::vector<std::vector<int> > posMat)
    8266{
    … …  
    9781
    9882
    99 
     83//print ideal.
    10084void id_print(ideal h)
    10185{
    … …  
    10690    pWrite(h->m[i]);
    10791    PrintLn();
    108     //PrintS(",");
    109   }
    110 }
    111 
    112 
    113 
    114 
    115 /************************only for T^2**********************************/
     92  }
     93}
     94
     95
     96
     97
     98//only for T^2,
     99//print vector of polynomials.
    116100void lpprint( std::vector<poly> pv)
    117101{
    … …  
    130114
    131115
    132 
     116//print vector of vectors of polynomials.
    133117void lpsprint(std::vector<std::vector<poly> > pvs)
    134118{
    … …  
    152136
    153137
    154 /*****************************About simplicial complex and stanly-reisner ring**************************/
    155 
    156 
    157 
    158 
    159 
    160 
    161 
    162 
    163 
     138
     139
     140
     141
     142/*************operations for vectors (regard vectors as sets)*********/
     143
     144//returns true if integer n is in vector vec,
     145//otherwise, returns false
     146bool IsinL(int a, std::vector<int> vec)
     147{
     148  int i;
     149  for(i=0;i<vec.size();i++)
     150  {
     151    if(a==vec[i])
     152    {
     153      return true;
     154    }
     155  }
     156  return false;
     157}
     158
     159
     160
     161
     162
     163//returns intersection of vectors p and q,
     164//returns empty if they are disjoint
     165std::vector<int> vecIntersection(std::vector<int> p, std::vector<int> q)
     166{
     167  int i;
     168  std::vector<int> inte;
     169  for(i=0;i<p.size();i++)
     170  {
     171    if(IsinL(p[i],q))
     172      inte.push_back(p[i]);
     173  }
     174  return inte;
     175}
     176
     177
     178
     179
     180
     181
     182
     183
     184
     185//returns true if vec1 is equal to vec2 (strictly equal, including the order)
     186//is not used
     187bool vEv(std::vector<int> vec1,std::vector<int> vec2)
     188{
     189  int i,j, lg1=vec1.size(),lg2=vec2.size();
     190  if(lg1!=lg2)
     191  {
     192    return false;
     193  }
     194  else
     195  {
     196    for(j=0;j<vec1.size();j++)
     197    {
     198      if(vec1[j]!=vec2[j])
     199        return false;
     200    }
     201  }
     202  return true;
     203}
     204
     205
     206
     207
     208//returns true if vec1 is contained in vec2
     209bool vsubset(std::vector<int> vec1, std::vector<int> vec2)
     210{
     211  int i;
     212  if(vec1.size()>vec2.size())
     213    return false;
     214  for(i=0;i<vec1.size();i++)
     215  {
     216    if(!IsinL(vec1[i],vec2))
     217      return false;
     218  }
     219  return true;
     220}
     221
     222//not strictly equal(order doesn't matter)
     223bool vEvl(std::vector<int> vec1, std::vector<int> vec2)
     224{
     225  if(vec1.size()==0 && vec2.size()==0)
     226    return true;
     227  if(vsubset(vec1,vec2)&&vsubset(vec2,vec1))
     228    return true;
     229  return false;
     230}
     231
     232
     233//the length of vec must be same to it of the elements of vecs
     234//returns true if vec is as same as some element of vecs ((not strictly same))
     235//returns false if vec is not in vecs
     236bool vInvsl(std::vector<int> vec, std::vector<std::vector<int> > vecs)
     237{
     238  int i;
     239  for(i=0;i<vecs.size();i++)
     240  {
     241    if(vEvl(vec,vecs[i]))
     242    {
     243      return true;
     244    }
     245  }
     246  return false;
     247}
     248
     249
     250//the length of vec must be same to it of the elements of vecs (strictly same)
     251//returns the position of vec in vecs,
     252//returns -1 if vec is not in vecs
     253//actrually is not used.
     254int vInvs(std::vector<int> vec, std::vector<std::vector<int> > vecs)
     255{
     256  int i;
     257  for(i=0;i<vecs.size();i++)
     258  {
     259    if(vEv(vec,vecs[i]))
     260    {
     261      return i+1;
     262    }
     263  }
     264  return -1;
     265}
     266
     267
     268
     269//returns the union of two vectors(as the union of sets)
     270std::vector<int> vecUnion(std::vector<int> vec1, std::vector<int> vec2)
     271{
     272  std::vector<int> vec=vec1;
     273  int i;
     274  for(i=0;i<vec2.size();i++)
     275  {
     276    if(!IsinL(vec2[i],vec))
     277      vec.push_back(vec2[i]);
     278  }
     279  return vec;
     280}
     281
     282
     283
     284std::vector<int> vecMinus(std::vector<int> vec1,std::vector<int> vec2)
     285{
     286  std::vector<int> vec;
     287  for(int i=0;i<vec1.size();i++)
     288  {
     289    if(!IsinL(vec1[i],vec2))
     290    {
     291      vec.push_back(vec1[i]);
     292    }
     293  }
     294  return vec;
     295}
     296
     297
     298
     299
     300
     301
     302std::vector<std::vector<int> > vsMinusv(std::vector<std::vector<int> > vecs, std::vector<int> vec)
     303{
     304  int i;
     305  std::vector<std::vector<int> > rem;
     306  for(i=0;i<vecs.size();i++)
     307  {
     308    if(!vEvl(vecs[i],vec))
     309    {
     310      rem.push_back(vecs[i]);
     311    }
     312  }
     313  return (rem);
     314}
     315
     316
     317std::vector<std::vector<int> > vsUnion(std::vector<std::vector<int> > vs1, std::vector<std::vector<int> > vs2)
     318{
     319  int i;
     320  std::vector<std::vector<int> > vs=vs1;
     321  for(i=0;i<vs2.size();i++)
     322  {
     323    if(!vInvsl(vs2[i],vs))
     324    {
     325      vs.push_back(vs2[i]);
     326    }
     327  }
     328  return vs;
     329}
     330
     331
     332
     333
     334
     335
     336std::vector<std::vector<int> > vsIntersection(std::vector<std::vector<int> > vs1, std::vector<std::vector<int> > vs2)
     337{
     338  int i;
     339  std::vector<std::vector<int> > vs;
     340  for(i=0;i<vs2.size();i++)
     341  {
     342    if(vInvsl(vs2[i],vs1))
     343    {
     344      vs.push_back(vs2[i]);
     345    }
     346  }
     347  return vs;
     348}
     349
     350
     351
     352
     353
     354/*************************************for transition between ideal and vectors******************************************/
     355
     356//P should be monomial,
     357// vector version of poly support(poly p)
     358std::vector<int> support1(poly p)
     359{
     360  int j;
     361  std::vector<int> supset;
     362  if(p==0) return supset;
     363  for(j=1;j<=rVar(currRing);j++)
     364  {
     365    if(pGetExp(p,j)>0)
     366    {
     367      supset.push_back(j);
     368    }
     369  }
     370  return (supset);
     371}
     372
     373
     374
     375
     376
     377
     378//simplicial complex(the faces set is ideal h)
     379std::vector<std::vector<int> >  supports(ideal h)
     380{
     381  std::vector<std::vector<int> > vecs;
     382  std::vector<int> vec;
     383  if(!idIs0(h))
     384  {
     385    for(int s=0;s<IDELEMS(h);s++)
     386    {
     387      vec=support1(h->m[s]);
     388      vecs.push_back(vec);
     389    }
     390  }
     391  return vecs;
     392}
     393
     394
     395
     396
     397// only for eqsolve1
     398// p could be any polynomial
     399std::vector<int> support2(poly p)
     400{
     401  int j;
     402  poly q;
     403  std::vector<int> supset;
     404  for(j=1;j<=rVar(currRing);j++)
     405  {
     406    q=pCopy(p);
     407    while (q!=NULL)
     408    {
     409      if(p_GetExp(q,j,currRing)!=0)
     410      {
     411        supset.push_back(j);
     412        break;
     413      }
     414      q=pNext(q);
     415    }
     416  }
     417  return (supset);
     418}
     419
     420
     421
     422//the supports of ideal
     423std::vector<std::vector<int> >  supports2(ideal h)
     424{
     425  std::vector<std::vector<int> > vecs;
     426  std::vector<int> vec;
     427  if(!idIs0(h))
     428  {
     429    for(int s=0;s<IDELEMS(h);s++)
     430    {
     431      vec=support2(h->m[s]);
     432      vecs.push_back(vec);
     433    }
     434  }
     435  return vecs;
     436}
     437//convert the vector(vbase[i] are the coefficients of x_{i+1}) to a polynomial w.r.t. current ring
     438//vector vbase has length of currRing->N.
     439poly pMake(std::vector<int> vbase)
     440{
     441  int n=vbase.size(); poly p,q=0;
     442  for(int i=0;i<n;i++)
     443  {
     444    if(vbase[i]!=0)
     445    {
     446      p = pOne();pSetExp(p, i+1, 1);pSetm(p);pSetCoeff(p, nInit(vbase[i]));
     447      q = pAdd(q, p);
     448    }
     449
     450  }
     451  return q;
     452}
     453
     454
     455
     456
     457//convert the vectors to a ideal(for T^1)
     458ideal idMake(std::vector<std::vector<int> > vecs)
     459{
     460  int lv=vecs.size(), i, j;
     461  poly p;
     462  ideal id_re=idInit(1,1);
     463  for(i=0;i<lv;i++)
     464  {
     465    p=pMake(vecs[i]);
     466    idInsertPoly(id_re, p);
     467  }
     468  idSkipZeroes(id_re);
     469  return id_re;
     470}
     471
     472
     473
     474/*****************************quotient ring of two ideals*********************/
     475
     476//the quotient ring of h1 respect to h2
     477ideal idmodulo(ideal h1,ideal h2)
     478{
     479  int i;
     480  ideal gb=kStd(h2,NULL,testHomog,NULL,NULL,0,0,NULL);
     481  idSkipZeroes(gb);
     482  ideal idq=kNF(gb,NULL,h1);
     483  idSkipZeroes(idq);
     484  return idq;
     485}
     486
     487//returns the coeff of the monomial of polynomial p which involves the mth varialbe
     488//assume the polynomial p has form of y1+y2+...
     489int pcoef(poly p, int m)
     490{
     491  int i,j,co; poly q=pCopy(p);
     492  for(i=1;i<=currRing->N;i++)
     493  {
     494    if(p_GetExp(q,m,currRing)!=0)
     495    {
     496      co=n_Int(pGetCoeff(q),currRing->cf);
     497      return co;
     498    }
     499    else
     500      q=pNext(q);
     501  }
     502  if(q!=NULL)
     503    co=0;
     504  return co;
     505}
     506
     507//returns true if p involves the mth variable of the current ring
     508bool vInp(int m,poly p)
     509{
     510  int i;
     511  poly q=pCopy(p);
     512  while (q!=NULL)
     513  {
     514    if(p_GetExp(q,m,currRing)!=0)
     515    {
     516      return true;
     517    }
     518    q=pNext(q);
     519  }
     520  return false;
     521}
     522
     523
     524
     525//returns the vector w.r.t. polynomial p
     526std::vector<int> vMake(poly p)
     527{
     528  int i; poly q=pCopy(p);
     529  std::vector<int> vbase;
     530  for(i=1;i<=currRing->N;i++)
     531  {
     532    if(vInp(i,p))
     533    {
     534      vbase.push_back(pcoef(p,i));
     535    }
     536    else
     537    {
     538      vbase.push_back(0);
     539    }
     540  }
     541  return (vbase);
     542}
     543
     544
     545//returns the vectors w.r.t. ideal h
     546std::vector<std::vector<int> > vsMake(ideal h)
     547{
     548  std::vector<int> vec;
     549  std::vector<std::vector<int> > vecs;
     550  int i;
     551  for(i=0;i<IDELEMS(h);i++)
     552  {
     553    vec=vMake(h->m[i]);
     554    vecs.push_back(vec);
     555  }
     556  return vecs;
     557}
     558
     559
     560//the quotient ring of two ideals which are represented by vectors,
     561//the result is also represented by vector.
     562std::vector<std::vector<int> > vecqring(std::vector<std::vector<int> > vec1, std::vector<std::vector<int> > vec2)
     563{
     564  int i,j;
     565  ideal h1=idMake(vec1), h2=idMake(vec2);
     566  ideal h=idmodulo(h1,h2);
     567  std::vector<std::vector<int> > vecs= vsMake(h);
     568  return vecs;
     569}
     570
     571
     572
     573/****************************************************************/
     574//construct a monomial based on the support of it
     575//returns a squarefree monomial
     576poly pMaken(std::vector<int> vbase)
     577{
     578  int n=vbase.size();
     579  poly p,q=pOne();
     580  for(int i=0;i<n;i++)
     581  {
     582    p = pOne();pSetExp(p, vbase[i], 1);pSetm(p);pSetCoeff(p, nInit(1));
     583    //pWrite(p);
     584    q=pp_Mult_mm(q,p,currRing);
     585  }
     586  return q;
     587}
     588
     589// returns a ideal according to a set of supports
     590ideal idMaken(std::vector<std::vector<int> > vecs)
     591{
     592  ideal id_re=idInit(1,1);
     593  poly p;
     594  int i,lv=vecs.size();
     595  for(i=0;i<lv;i++)
     596  {
     597    p=pMaken(vecs[i]);
     598    idInsertPoly(id_re, p);
     599  }
     600  idSkipZeroes(id_re);
     601  //id_print(id_re);
     602  return id_re;
     603}
     604
     605
     606
     607/********************************new version for stanley reisner ideal ***********************************************/
     608
     609
     610std::vector<std::vector<int> > b_subsets(std::vector<int> vec)
     611{
     612  int i,j;
     613  std::vector<int> bv;
     614  std::vector<std::vector<int> > vecs;
     615  for(i=0;i<vec.size();i++)
     616  {
     617    bv.push_back(vec[i]);
     618    vecs.push_back(bv);
     619    bv.clear();
     620  }
     621  //listsprint(vecs);
     622  for(i=0;i<vecs.size();i++)
     623  {
     624    for(j=i+1;j<vecs.size();j++)
     625    {
     626      bv=vecUnion(vecs[i], vecs[j]);
     627      if(!vInvsl(bv,vecs))
     628        vecs.push_back(bv);
     629    }
     630  }
     631  //listsprint(vecs);
     632  return(vecs);
     633}
     634
     635
     636//the number of the variables
     637int idvert(ideal h)
     638{
     639  int i, j, vert=0;
     640  if(idIs0(h))
     641    return vert;
     642  for(i=currRing->N;i>0;i--)
     643  {
     644    for(j=0;j<IDELEMS(h);j++)
     645    {
     646      if(pGetExp(h->m[j],i)>0)
     647      {
     648        vert=i;
     649        return vert;
     650      }
     651    }
     652  }
     653  return vert;
     654}
     655
     656
     657
     658
     659int pvert(poly p)
     660{
     661  int i, j, vert=0;
     662  for(i=currRing->N;i>0;i--)
     663  {
     664      if(pGetExp(p,i)>0)
     665      {
     666        vert=i;
     667        return vert;
     668      }
     669  }
     670  return vert;
     671}
     672
     673
     674/*
     675//full complex
     676std::vector<std::vector<int> > fullcomplex(ideal h)
     677{
     678  int vert=vertnum(h), i, j;
     679  //Print("there are %d vertices\n", vert);
     680  std::vector<std::vector<int> > fmons;
     681  std::vector<int> pre;
     682  for(i=1;i<=vert;i++)
     683  {
     684    pre.push_back(i);
     685  }
     686  fmons=b_subsets(pre);
     687  return fmons;
     688
     689}*/
     690
     691
     692/*
     693//all the squarefree monomials whose degree is less or equal to n
     694std::vector<std::vector<int> > sfrmons(ideal h, int n)
     695{
     696  int vert=vertnum(h), i, j, time=0;
     697  std::vector<std::vector<int> > fmons, pres, pres0, pres1;
     698  std::vector<int> pre;
     699  for(i=1;i<=vert;i++)
     700  {
     701    pre.push_back(i);
     702    pres0.push_back(pre);
     703  }
     704  pres=pres0;
     705  for(i=0;i<size(pres),time<=n;i++)
     706  {
     707    time++;
     708    pre=pres[i];
     709    for(j=0;j<size(pres0);j++)
     710    {
     711      pre=vecUnion(pre, pres0[j]);
     712      if(pre.)
     713    }
     714  }
     715  return fmons;
     716
     717}
     718*/
     719
     720/*
     721ideal id_complement(ideal h)
     722{
     723  int i,j;
     724  std::vector<std::vector<int> > full=fullcomplex(h), hvs=supports(h), res;
     725  for(i=0;i<full.size();i++)
     726  {
     727    if(!vInvsl(full[i], hvs))
     728    {
     729      res.push_back(full[i]);
     730    }
     731  }
     732  return idMaken(res);
     733}*/
     734
     735
     736/*****************About simplicial complex and stanley-reisner ideal and ring **************************/
    164737
    165738//h1 minus h2
    … …  
    173746    for(j=0;j<IDELEMS(h2);j++)
    174747    {
    175       if(p_EqualPolys(h1->m[i],h2->m[j], currRing))
     748      if(p_EqualPolys(pCopy(h1->m[i]),pCopy(h2->m[j]), currRing))
    176749      {
    177750        eq=1;
    … …  
    181754    if(eq==0)
    182755    {
    183       idInsertPoly(h, h1->m[i]);
     756      idInsertPoly(h, pCopy(h1->m[i]));
    184757    }
    185758  }
    … …  
    206779}
    207780
    208 /*find all squarefree monomials of degree deg in ideal h*/
     781
     782
     783//returns the set of all squarefree monomials of degree deg in ideal h
    209784ideal sfreemon(ideal h,int deg)
    210785{
    … …  
    232807
    233808
    234 
    235 
    236809//full simplex represented by ideal.
    237810//(all the squarefree monomials over the polynomial ring)
    238 ideal id_sfmon()
     811ideal id_sfmon(ideal h)
    239812{
    240813  ideal asfmons,sfmons,mons,p;
    241   int j;
     814  int j, vert=idvert(h);
    242815  mons=id_MaxIdeal(1, currRing);
    243816  asfmons=sfreemon(mons,1);
    244   for(j=2;j<=rVar(currRing);j++)
     817  for(j=2;j<=vert;j++)
    245818  {
    246819    mons=id_MaxIdeal(j, currRing);
    … …  
    256829
    257830
    258 
    259831//if the input ideal is simplicial complex, returns the stanley-reisner ideal,
    260832//if the input ideal is stanley-reisner ideal, returns the monomial ideal according to simplicial complex.
    261833//(nonfaces and faces).
    262 //returns the complement of the ideal h
     834//returns the complement of the ideal h (consisting of only squarefree polynomials)
    263835ideal id_complement(ideal h)
    264836{
    265   int j;
    266   ideal i1=id_sfmon();
    267   ideal i2=idMinus(i1,h);
     837  int j, vert=idvert(h);
     838  ideal i1=id_sfmon(h);
     839  ideal i3=idInit(1,1);
     840  poly p;
     841  for(j=0;j<IDELEMS(i1);j++)
     842  {
     843    p=pCopy(i1->m[j]);
     844    if(pvert(p)<=vert)
     845    {
     846      idInsertPoly(i3, p);
     847    }
     848  }
     849  ideal i2=idMinus(i3,h);
     850  idSkipZeroes(i2);
    268851  return (i2);
    269852}
    … …  
    272855
    273856
    274 
    275 
    276 
    277 
    278 
    279 
    280 
    281 
    282 //returns the monomials in the quotient ring R/(h1+h2) which have degree deg,
    283 //and R is currRing
    284 ideal qringadd(ideal h1, ideal h2,int deg)
     857//Returns true if p is one of the generators of ideal X
     858//returns false otherwise
     859bool IsInX(poly p,ideal X)
     860{
     861  int i,j;
     862  for(i=0;i<IDELEMS(X);i++)
     863  {
     864    if(pEqualPolys(p,X->m[i]))
     865    {
     866      //PrintS("yes\n");
     867      return(true);
     868    }
     869  }
     870  //PrintS("no\n");
     871  return(false);
     872}
     873
     874
     875
     876
     877
     878
     879//returns the monomials in the quotient ring R/(h1+h2) which have degree deg.
     880ideal qringadd(ideal h1, ideal h2, int deg)
    285881{
    286882  ideal h,qrh;
    … …  
    301897  for(i=1;i<IDELEMS(h);i++)
    302898  {
    303     if(pTotaldegree(h->m[i])>max)
    304     max=pTotaldegree(h->m[i]);
     899    if(pTotaldegree(h->m[i]) > max)
     900      max=pTotaldegree(h->m[i]);
    305901  }
    306902  return (max);
    307903}
    308904
    309 //input ideal h is the ideal associated to simplicial complex,
     905
     906
     907
     908
     909
     910
     911//input ideal h (a squarefree monomial ideal) is the ideal associated to simplicial complex,
    310912//and returns the Stanley-Reisner ideal(minimal generators)
    311913ideal idsrRing(ideal h)
    312914{
    313915  int max,i,j,n;
    314   ideal hc=idCopy(h);
    315   //Print("This is the complement generators\n");
    316   //id_print(hc);
    317   ideal pp,qq,rsr,ppp;
     916  ideal pp,qq,rsr,ppp,hc=idCopy(h);
    318917  for(i=1;i<=rVar(currRing);i++)
    319918  {
    … …  
    324923      pp=sfreemon(pp,i);
    325924      rsr=pp;
    326   //Print("This is the first quotient generators %d:\n",i);
    327   //id_print(rsr);
     925      //Print("This is the first quotient generators %d:\n",i);
     926      //id_print(rsr);
    328927      break;
    329928    }
    … …  
    332931  {
    333932    qq=sfreemon(hc,n);
    334     if(!idIs0(qq))
    335     {
    336       pp=qringadd(qq,rsr,n);
    337       ppp=sfreemon(pp,n);
    338   //Print("This is the quotient generators %d:\n",n);
    339   //id_print(ppp);
    340       rsr=idAdd(rsr,ppp);
    341   //Print("This is the current minimal set\n");
    342   //id_print(rsr);
    343     }
     933    pp=qringadd(qq,rsr,n);
     934    ppp=sfreemon(pp,n);
     935    rsr=idAdd(rsr,ppp);
    344936  }
    345937  idSkipZeroes(rsr);
    346   //PrintS("This is the minimal generators:\n");
    347   //id_print(rsr);
    348938  return rsr;
    349939}
    … …  
    351941
    352942
    353 //returns  the set of subset of p
     943//returns the set of all the polynomials could divide p
    354944ideal SimFacset(poly p)
    355945{
    356   int i,j;
    357   ideal h1,mons;
    358   int max=pTotaldegree(p);
    359   ideal id_re=idInit(1,1);
     946  int i,j,max=pTotaldegree(p);
     947  ideal h1,mons,id_re=idInit(1,1);
    360948  for(i=1;i<max;i++)
    361949  {
    362950    mons=id_MaxIdeal(i, currRing);
    363951    h1=sfreemon(mons,i);
     952
    364953    for(j=0;j<IDELEMS(h1);j++)
    365954    {
    366955      if(p_DivisibleBy(h1->m[j],p,currRing))
     956      {
    367957        idInsertPoly(id_re, h1->m[j]);
    368     }
     958      }
     959    }
     960
    369961  }
    370962  idSkipZeroes(id_re);
    371   //PrintS("This is the facset\n");
    372   //id_print(id_re);
    373963  return id_re;
    374964}
    375965
    376 //complicated version(print the polynomial set which can make the input to be a simplex)
    377 //input h is at least part of faces
    378 bool IsSimplex(ideal h)
    379 {
    380   int i,j,ifbreak=0;
    381   ideal id_re;
    382   ideal id_so=idCopy(h);
    383   int max=id_maxdeg(h);
     966
     967
     968ideal idadda(ideal h1, ideal h2)
     969{
     970  ideal h=idInit(1,1);
     971  for(int i=0;i<IDELEMS(h1);i++)
     972  {
     973    if(!IsInX(h1->m[i],h))
     974    {
     975      idInsertPoly(h, h1->m[i]);
     976    }
     977  }
     978  for(int i=0;i<IDELEMS(h2);i++)
     979  {
     980    if(!IsInX(h2->m[i],h))
     981    {
     982      idInsertPoly(h, h2->m[i]);
     983    }
     984  }
     985  idSkipZeroes(h);
     986  return h;
     987}
     988
     989
     990//complicated version
     991//(returns false if it is not a simplicial complex and print the simplex)
     992//input h is need to be at least part of faces
     993ideal IsSimplex(ideal h)
     994{
     995  int i,j,ifbreak=0,max=id_maxdeg(h);
     996  poly e=pOne();
     997  ideal id_re, id_so=idCopy(h);
    384998  for(i=0;i<IDELEMS(h);i++)
    385999  {
    3861000    id_re=SimFacset(h->m[i]);
    387     //id_print(id_re);
    3881001    if(!idIs0(id_re))
    3891002    {
    390       id_so=idAdd(id_so, id_re);
    391     }
    392   }
     1003      id_so=idadda(id_so, id_re);//idAdd(id_so,id_re);
     1004    }
     1005  }
     1006  idInsertPoly(id_so,e);
    3931007  idSkipZeroes(id_so);
    394   if(!idIs0(idMinus(id_so,id_re)))
    395   {
    396     PrintS("It is not simplex.\n");
    397     PrintS("This is the simplicial complex:\n");
    398     id_print(id_so);
    399     return false;
    400   }
    401   PrintS("It is simplex.\n");
    402   return true;
    403 }
    404 
    405 
    406 //input is the subset of the Stainly-Reisner ideal
     1008  return (idMinus(id_so,h));
     1009}
     1010
     1011
     1012//input is the subset of the Stainley-Reisner ideal
    4071013//returns the faces
     1014//is not used
    4081015ideal complementsimplex(ideal h)
    4091016{
    410   int i,j;poly p;
    411   ideal h1=idInit(1,1);
    412   ideal pp;
    413   ideal h3=idInit(1,1);
     1017  int i,j;poly p,e=pOne();
     1018  ideal h1=idInit(1,1), pp, h3;
    4141019  for(i=1;i<=rVar(currRing);i++)
    4151020  {
    4161021    p = pOne(); pSetExp(p, i, 2); pSetm(p); pSetCoeff(p, nInit(1));
    417     //pWrite(p);
    4181022    idInsertPoly(h1, p);
    4191023  }
    … …  
    4271031    h3=idAdd(h3,pp);
    4281032  }
    429   PrintS("This is the simplicial complex:\n");
    430   id_print(h3);
     1033  idInsertPoly(h3, e);
     1034  idSkipZeroes(h3);
    4311035  return (h3);
    4321036}
    … …  
    4341038
    4351039
    436 
    437 /*************operations for vectors(sets)*********/
    438 
    439 //returns true if integer n is in vector vec,
    440 //otherwise returns false
    441 bool IsinL(int a,std::vector<int> badset)
    442 {
    443   int i;
    444   for(i=0;i<badset.size();i++)
    445   {
    446     if(a==badset[i])
    447     {
    448       return true;
    449     }
    450   }
    451   return false;
    452 }
    453 
    454 
    455 
    456 
    457 
    458 //intersection of vectors p and q, returns empty if they are disjoint
    459 std::vector<int> vecIntersection(std::vector<int> p,std::vector<int> q)
    460 {
    461   int i;
    462   std::vector<int> inte;
    463   for(i=0;i<p.size();i++)
    464   {
    465     if(IsinL(p[i],q))
    466       inte.push_back(p[i]);
    467   }
    468   //listprint(inte);
    469   return inte;
    470 }
    471 
    472 
    473 
    474 
    475 
    476 
    477 
    478 
    479 
    480 //returns true if vec1 is equal to vec2 (including the order)
    481 bool vEv(std::vector<int> vec1,std::vector<int> vec2)
    482 {
    483   int i,j;
    484   int lg1=vec1.size(),lg2=vec2.size();
    485   if(lg1!=lg2)
    486   {
    487     return false;
    488   }
    489   else
    490   {
    491     for(j=0;j<vec1.size();j++)
    492     {
    493       if(vec1[j]!=vec2[j])
    494         return false;
    495     }
    496   }
    497   return true;
    498 }
    499 
    500 
    501 
    502 
    503 //returns true if vec1 is contained in vec2
    504 bool vsubset(std::vector<int> vec1,std::vector<int> vec2)
    505 {
    506   int i;
    507   if(vec1.size()>vec2.size())
    508     return false;
    509   for(i=0;i<vec1.size();i++)
    510   {
    511     if(!IsinL(vec1[i],vec2))
    512       return false;
    513   }
    514   return true;
    515 }
    516 
    517 //not strictly equal(order doesn't matter)
    518 bool vEvl(std::vector<int> vec1,std::vector<int> vec2)
    519 {
    520   if(vec1.size()==0 && vec2.size()==0)
    521     return true;
    522   if(vsubset(vec1,vec2)&&vsubset(vec2,vec1))
    523     return true;
    524   return false;
    525 }
    526 
    527 
    528 //the length of vec must be same to it of the elements of vecs (not strictly same)
    529 //returns false if vec is not in vecs
    530 bool vInvsl(std::vector<int> vec, std::vector<std::vector<int> > vecs)
    531 {
    532   int i;
    533   for(i=0;i<vecs.size();i++)
    534   {
    535     if(vEvl(vec,vecs[i]))
    536     {
    537       return true;
    538     }
    539   }
    540   return false;
    541 }
    542 
    543 
    544 //the length of vec must be same to it of the elements of vecs (strictly same)
    545 //returns the position of vec in vecs,
    546 //returns -1 if vec is not in vecs
    547 int vInvs(std::vector<int> vec, std::vector<std::vector<int> > vecs)
    548 {
    549   int i;
    550   for(i=0;i<vecs.size();i++)
    551   {
    552     if(vEv(vec,vecs[i]))
    553     {
    554       //Print("This is the %dth variable\n",i+1);
    555       return i+1;
    556     }
    557   }
    558   //Print("This is not the new variable\n");
    559   //listprint(vec);
    560   return -1;
    561 }
    562 
    563 
    564 
    565 //returns the union of two vectors(like the union of sets)
    566 std::vector<int> vecUnion(std::vector<int> vec1, std::vector<int> vec2)
    567 {
    568   std::vector<int> vec=vec1;
    569   int i;
    570   for(i=0;i<vec2.size();i++)
    571   {
    572     if(!IsinL(vec2[i],vec))
    573       vec.push_back(vec2[i]);
    574   }
    575   return vec;
    576 }
    577 
    578 
    579 
    580 std::vector<int> vecMinus(std::vector<int> vec1,std::vector<int> vec2)
    581 {
    582   std::vector<int> vec;
    583   for(int i=0;i<vec1.size();i++)
    584   {
    585     if(!IsinL(vec1[i],vec2))
    586     {
    587       vec.push_back(vec1[i]);
    588     }
    589   }
    590   return vec;
    591 }
    592 
     1040int dim_sim(ideal h)
     1041{
     1042  int dim=pTotaldegree(h->m[0]), i;
     1043  for(i=1; i<IDELEMS(h);i++)
     1044  {
     1045    if(dim<pTotaldegree(h->m[i]))
     1046    {
     1047      dim=pTotaldegree(h->m[i]);
     1048    }
     1049  }
     1050  return dim;
     1051}
     1052
     1053
     1054int num4dim(ideal h, int n)
     1055{
     1056  int num=0;
     1057  for(int i=0; i<IDELEMS(h); i++)
     1058  {
     1059    if(pTotaldegree(h->m[i])==n)
     1060    {
     1061      num++;
     1062    }
     1063  }
     1064  return num;
     1065}
    5931066
    5941067
    … …  
    6001073
    6011074
    602 
    603 
    604 
    605 
    606 //P should be monomial(not used) vector version of poly support(poly p)
    607 std::vector<int> support1(poly p)
    608 {
    609   int j;
    610   std::vector<int> supset;
    611   for(j=1;j<=rVar(currRing);j++)
    612   {
    613     if(pGetExp(p,j)>0)
    614     {
    615       supset.push_back(j);
    616     }
    617   }
    618   return (supset);
    619 }
    620 
    621 //simplicial complex(the faces set is ideal h)
    622 std::vector<std::vector<int> >  supports(ideal h)
    623 {
    624   std::vector<std::vector<int> > vecs;
    625   std::vector<int> vec;
    626   if(!idIs0(h))
    627   {
    628     for(int s=0;s<IDELEMS(h);s++)
    629     {
    630       //pWrite(h->m[s]);
    631       vec=support1(h->m[s]);
    632       vecs.push_back(vec);
    633     }
    634   }
    635   return vecs;
    636 }
    637 
    638 
    639 
    640 
    641 
    6421075//h is ideal( monomial ideal) associated to simplicial complex
    643 //returns the all the monomials x^b (x^b must be able to divide at least one monomial in Stanley-Reisner ring)
     1076//returns the all the monomials x^b (x^b must be able to divide
     1077//at least one monomial in Stanley-Reisner ring)
     1078//not so efficient
    6441079ideal findb(ideal h)
    6451080{
    646   ideal ib=id_sfmon();
    647   ideal nonf=id_complement(h);
    648   ideal bset=idInit(1,1);
     1081  ideal ib=id_sfmon(h), nonf=id_complement(h), bset=idInit(1,1);
    6491082  poly e=pOne();
    6501083  int i,j;
    … …  
    6751108ideal finda(ideal h,poly S,int ddeg)
    6761109{
    677   ideal aset=idInit(1,1);
    6781110  poly e=pOne();
    679   ideal h2=id_complement(h);
    680   poly p;
     1111  ideal h2=id_complement(h), aset=idInit(1,1);
    6811112  int i,j,deg1=pTotaldegree(S);
    6821113  int tdeg=deg1+ddeg;
    … …  
    6881119    for(i=0;i<IDELEMS(ia);i++)
    6891120    {
    690     //p=support(ia->m[i]);
    6911121      v=support1(ia->m[i]);
    6921122      in=vecIntersection(v,bv);
    … …  
    6961126      }
    6971127    }
    698     //idInsertPoly(aset,e);
    6991128    idSkipZeroes(aset);
    7001129  }
    … …  
    7121141//returns true if support(p) union support(a) minus support(b) is face,
    7131142//otherwise returns false
    714 //not be used yet, (the vector version of mabcondition)
     1143//(the vector version of mabcondition)
    7151144bool mabconditionv(std::vector<std::vector<int> > hvs,std::vector<int> pv,std::vector<int> av,std::vector<int> bv)
    7161145{
    7171146  std::vector<int> uv=vecUnion(pv,av);
    7181147  uv=vecMinus(uv,bv);
    719 //Print("this is the support of p union a minus b\n");
    720 //listprint(uv);
    721   //Print("this is the support of ideal\n");
    722   //listsprint(hvs);
    7231148  if(vInvsl(uv,hvs))
    7241149  {
    … …  
    7291154
    7301155
     1156// returns the set of nonfaces p where mabconditionv(h, p, a, b) is true
    7311157std::vector<std::vector<int> > Mabv(ideal h,poly a,poly b)
    7321158{
    733   std::vector<int> pv;
    734   std::vector<std::vector<int> > vecs;
    735   //Print("this is the support of p\n");
    736   //listprint(pv);
    737   std::vector<int> av=support1(a);
    738 //Print("this is the support of a\n");
    739 //listprint(av);
    740   std::vector<int> bv=support1(b);
    741 //Print("this is the support of b\n");
    742 //listprint(bv);
     1159  std::vector<int> av=support1(a), bv=support1(b), pv, vec;
    7431160  ideal h2=id_complement(h);
    744   std::vector<std::vector<int> > hvs=supports(h);
    745   std::vector<std::vector<int> > h2v=supports(h2);
    746   std::vector<int> vec;
     1161  std::vector<std::vector<int> > hvs=supports(h), h2v=supports(h2), vecs;
    7471162  for(int i=0;i<h2v.size();i++)
    7481163  {
    … …  
    7671182
    7681183/***************************************************************************/
    769 //For solving the equations which has form of x_i-x_j.(equations got by T_1)
     1184//For solving the equations which has form of x_i-x_j.(equations got from T_1)
    7701185/***************************************************************************/
    7711186
    … …  
    7901205}
    7911206
    792 
     1207/*
    7931208//get triangular form(eqs.size()>0)
    7941209std::vector<std::vector<int> > soleli1( std::vector<std::vector<int> > eqs)
    7951210{
    796   int i,j;std::vector<std::vector<int> > re;
    797   std::vector<std::vector<int> >  pre=eqs,ppre;
     1211  int i,j;
     1212  std::vector<int> yaya;
     1213  std::vector<std::vector<int> >  pre=eqs, ppre, re;
    7981214  if(eqs.size()>0)
    7991215  {
    … …  
    8011217    pre.erase(pre.begin());
    8021218  }
    803   //listsprint(pre);
    804   std::vector<int> yaya;
    805   for(i=0;(i<re.size()) && (pre.size()>0);i++)
     1219  for(i=0;i<re.size(),pre.size()>0;i++)
    8061220  {
    8071221    yaya=eli1(re[i],pre[0]);
    808     //listprint(yaya);
    8091222    re.push_back(yaya);
    8101223    for(j=1;j<pre.size();j++)
    8111224    {
    812       //listprint(pre[j]);
    8131225      ppre.push_back(eli1(re[i],pre[j]));
    8141226    }
    … …  
    8171229  }
    8181230  return re;
    819 }
    820 
    821 
    822 // input is a set of equations who is of triangular form(every equations has a form of x_i-x_j) n is the number of variables
     1231}*/
     1232//make sure the first element is smaller that the second one
     1233std::vector<int> keeporder(  std::vector<int> vec)
     1234{
     1235  std::vector<int> yaya;
     1236  int n;
     1237  if(vec[0]>vec[1])
     1238  {
     1239    n=vec[0];
     1240    vec[0]=vec[1];
     1241    vec[1]=n;
     1242  }
     1243  return vec;
     1244}
     1245
     1246
     1247std::vector<std::vector<int> > soleli1( std::vector<std::vector<int> > eqs)
     1248{
     1249  int i,j;
     1250  std::vector<int> yaya;
     1251  std::vector<std::vector<int> >  pre=eqs, ppre, re;
     1252  if(eqs.size()>0)
     1253  {
     1254    re.push_back(eqs[0]);
     1255    pre.erase(pre.begin());
     1256  }
     1257  while(pre.size()>0)
     1258  {
     1259    yaya=keeporder(eli1(re[0],pre[0]));
     1260    for(i=1;i<re.size();i++)
     1261    {
     1262      if(!vInvsl(yaya, re))
     1263      {
     1264        yaya=eli1(re[i],yaya);
     1265        yaya=keeporder(yaya);
     1266      }
     1267    }
     1268    if(!vInvsl(yaya, re))
     1269    {
     1270      re.push_back(yaya);
     1271    }
     1272    pre.erase(pre.begin());
     1273  }
     1274  return re;
     1275}
     1276
     1277
     1278
     1279// input is a set of equations who is of triangular form(every equations has a form of x_i-x_j)
     1280// n is the number of variables
    8231281//get the free variables and the dimension
    8241282std::vector<int> freevars(int n,  std::vector<int> bset, std::vector<std::vector<int> > gset)
    8251283{
    826   int ql=gset.size();
    827   int bl=bset.size();
    828   std::vector<int> mvar;
    829   std::vector<int> fvar;
     1284  int ql=gset.size(), bl=bset.size(), i;
     1285  std::vector<int> mvar, fvar;
     1286  for(i=0;i<bl;i++)
     1287  {
     1288    mvar.push_back(bset[i]);
     1289  }
     1290  for(i=0;i<ql;i++)
     1291  {
     1292    mvar.push_back(gset[i][0]);
     1293  }
     1294  for(i=1;i<=n;i++)
     1295  {
     1296    if(!IsinL(i,mvar))
     1297    {
     1298      fvar.push_back(i);
     1299    }
     1300  }
     1301    return fvar;
     1302}
     1303
     1304
     1305//return the set of free variables except the vnum one
     1306std::vector<int> fvarsvalue(int vnum, std::vector<int> fvars)
     1307{
    8301308  int i;
    831     for(i=0;i<bl;i++)
    832     {
    833       mvar.push_back(bset[i]);
    834     }
    835     for(i=0;i<ql;i++)
    836     {
    837       mvar.push_back(gset[i][0]);
    838     }
    839     for(i=1;i<=n;i++)
    840     {
    841       if(!IsinL(i,mvar))
    842       {
    843         fvar.push_back(i);
    844       }
    845     }
    846     return fvar;
    847 }
    848 
    849 
    850 //return the set of free variables except vnum
    851 std::vector<int> fvarsvalue(int vnum, std::vector<int> fvars)
    852 {
    853   int i;std::vector<int> fset=fvars;
     1309  std::vector<int> fset=fvars;
    8541310  for(i=0;i<fset.size();i++)
    8551311  {
    … …  
    8651321
    8661322
    867 
     1323//returns the simplified bset and gset
     1324//enlarge bset, simplify gset
    8681325std::vector<std::vector<int> > vAbsorb( std::vector<int> bset,std::vector<std::vector<int> > gset)
    8691326{
    870   int i,j,m;
    8711327  std::vector<int> badset=bset;
    872   int bl=badset.size();
    873   int gl=gset.size();
     1328  int i,j,m, bl=bset.size(), gl=gset.size();
    8741329  for(i=0;i<bl;i++)
    8751330  {
    … …  
    9151370
    9161371
    917 //the new variables are started from 1
     1372//the labels of new variables are started with 1
     1373//returns a vector of solution space according to index
    9181374std::vector<int> vecbase1(int num, std::vector<int> oset)
    9191375{
    … …  
    9661422std::vector<int> ofindbases1(int num, int vnum, std::vector<int> bset,std::vector<std::vector<int> > gset)
    9671423{
    968   int i,j,m;std::vector<std::vector<int> > goodset;
    969   std::vector<int> fvars=freevars(num,   bset,  gset);
     1424  int i,j,m;
     1425  std::vector<std::vector<int> > goodset;
     1426  std::vector<int> fvars=freevars(num,   bset,  gset), oset, base;
    9701427  std::vector<int> zset=fvarsvalue(vnum, fvars);
    9711428  zset=vecUnion(zset,bset);
    972   std::vector<int> oset;
    9731429  oset.push_back(vnum);
    9741430  goodset=vAbsorb(oset, gset);
    9751431  oset=goodset[goodset.size()-1];
    9761432  goodset.erase(goodset.end());
    977   //PrintS("This is the 1 set:\n");
    978   //listprint(oset);
    979   //goodset=vAbsorb(zset, goodset);
    980   //zset=goodset[goodset.size()-1];
    981   //goodset.erase(goodset.end());
    982   //PrintS("This is the 0 set:\n");
    983   //listprint(zset);
    984   //Print("thet size of goodset is %ld\n", goodset.size());
    985   std::vector<int> base= vecbase1(num,  oset);
     1433  base= vecbase1(num,  oset);
    9861434  return base;
    9871435}
    … …  
    9991447{
    10001448  int i,j,m;
    1001   std::vector<int> base1;
    10021449  std::vector<std::vector<int> > bases;
    1003   std::vector<int> fvars=freevars(num,   bset,  gset);
     1450  std::vector<int> fvars=freevars(num,   bset,  gset), base1;
    10041451  if (fvars.size()==0)
    10051452  {
    … …  
    10311478//num is the number of varialbes also the length of the set which we need to consider
    10321479//output is trigular form of gset and badset where x_i=0
    1033 std::vector<std::vector<int> > eli2(int num,std::vector<int> bset,std::vector<std::vector<int> > gset)
     1480std::vector<std::vector<int> > eli2(int num, std::vector<int> bset,std::vector<std::vector<int> > gset)
    10341481{
    10351482  int i,j;
    10361483  std::vector<int> badset;
    1037   std::vector<std::vector<int> > goodset;
    1038   std::vector<std::vector<int> > solve;
     1484  std::vector<std::vector<int> > goodset, solve;
     1485//PrintS("This is the input bset\n");listprint(bset);
     1486//PrintS("This is the input gset\n");listsprint(gset);
    10391487  if(gset.size()!=0)//gset is not empty
    10401488  {
    1041 //find all the variables which are zeroes
     1489   //find all the variables which are zeroes
    10421490
    10431491    if(bset.size()!=0)//bset is not empty
    … …  
    10521500      goodset=gset;//badset is empty
    10531501    }//goodset is already the set which doesn't contain zero variables
    1054 
     1502//PrintS("This is the badset after absorb \n");listprint(badset);
     1503//PrintS("This is the goodset after absorb \n");listsprint(goodset);
    10551504    goodset=soleli1(goodset);//get the triangular form of goodset
    1056     //PrintS("the triangular form of the good set is:\n");
     1505//PrintS("This is the goodset after triangulization \n");listsprint(goodset);
    10571506    solve=ofindbases(num,badset,goodset);
    1058    // goodset.push_back(badset);
    1059    //listsprint(goodset);
    10601507  }
    10611508  else
    … …  
    10631510    solve=ofindbases(num,bset,gset);
    10641511  }
     1512//PrintS("This is the solution\n");listsprint(solve);
    10651513  return solve;
    10661514}
    … …  
    10701518
    10711519
    1072 //convert the vector to a polynomial w.r.t. current ring
    1073 //vector vbase has length of currRing->N.
    1074 poly pMake(std::vector<int> vbase)
    1075 {
    1076   int n=vbase.size();poly p,q=0;
    1077   for(int i=0;i<n;i++)
    1078   {
    1079     if(vbase[i]!=0)
    1080     {
    1081       p = pOne();pSetExp(p, i+1, 1);pSetm(p);pSetCoeff(p, nInit(vbase[i]));
    1082       q = pAdd(q, p);
    1083     }
    1084 
    1085   }
    1086   return q;
    1087 }
    1088 
    1089 
    1090 
    1091 
    1092 //convert the vectors to a ideal(for T^1)
    1093 ideal idMake(std::vector<std::vector<int> > vecs)
    1094 {
    1095   int lv=vecs.size();poly p;
    1096   int i,j;
    1097   ideal id_re=idInit(1,1);
    1098   for(i=0;i<lv;i++)
    1099   {
    1100     p=pMake(vecs[i]);
    1101     idInsertPoly(id_re, p);
    1102   }
    1103   idSkipZeroes(id_re);
    1104   return id_re;
    1105 }
    1106 
    1107 
    1108 
    1109 /*****************************quotient ring of two ideals*********************/
    1110 
    1111 //the quotient ring of h1 respect to h2
    1112 ideal idmodulo(ideal h1,ideal h2)
    1113 {
    1114   int i;
    1115   ideal gb=kStd(h2,NULL,testHomog,NULL,NULL,0,0,NULL);
    1116   ideal idq=kNF(gb,NULL,h1);
    1117   idSkipZeroes(idq);
    1118   return idq;
    1119 }
    1120 
    1121 //returns the coeff of the monomial of polynomial p which involves the mth varialbe
    1122 //assume the polynomial p has form of y1+y2+...
    1123 int pcoef(poly p, int m)
    1124 {
    1125   int i,j,co;poly q=pCopy(p);
    1126   for(i=1;i<=currRing->N;i++)
    1127   {
    1128     if(p_GetExp(q,m,currRing)!=0)
    1129     {
    1130       co=n_Int(pGetCoeff(q),currRing->cf);
    1131       return co;
    1132     }
    1133     else
    1134       q=pNext(q);
    1135   }
    1136   if(q!=NULL)
    1137     co=0;
    1138   return co;
    1139 }
    1140 
    1141 //returns true if p involves the mth variable of the current ring
    1142 bool vInp(int m,poly p)
    1143 {
    1144   int i;
    1145   poly q=pCopy(p);
    1146   while (q!=NULL)
    1147   {
    1148     if(p_GetExp(q,m,currRing)!=0)
    1149     {
    1150       return true;
    1151     }
    1152     q=pNext(q);
    1153   }
    1154   return false;
    1155 }
    1156 
    1157 
    1158 
    1159 //returns the vector w.r.t. polynomial p
    1160 std::vector<int> vMake(poly p)
    1161 {
    1162   int i;poly q=pCopy(p);
    1163   std::vector<int> vbase;
    1164   for(i=1;i<=currRing->N;i++)
    1165   {
    1166     if(vInp(i,p))
    1167     {
    1168       vbase.push_back(pcoef(p,i));
    1169     }
    1170     else
    1171     {
    1172       vbase.push_back(0);
    1173     }
    1174   }
    1175   return (vbase);
    1176 }
    1177 
    1178 
    1179 //returns the vectors w.r.t. ideal h
    1180 std::vector<std::vector<int> > vsMake(ideal h)
    1181 {
    1182   std::vector<int> vec;
    1183   std::vector<std::vector<int> > vecs;
    1184   int i;
    1185   for(i=0;i<IDELEMS(h);i++)
    1186   {
    1187     vec=vMake(h->m[i]);
    1188     vecs.push_back(vec);
    1189   }
    1190   return vecs;
    1191 }
    1192 
    1193 
    1194 //the quotient ring of two ideals which are represented by vectors,
    1195 //the result is also represented by vector.
    1196 std::vector<std::vector<int> > vecqring(std::vector<std::vector<int> > vec1, std::vector<std::vector<int> > vec2)
    1197 {
    1198   int i,j;
    1199   ideal h1=idMake(vec1);
    1200   //id_print(h1);
    1201   ideal h2=idMake(vec2);
    1202   //id_print(h2);
    1203   ideal h=idmodulo(h1,h2);
    1204   std::vector<std::vector<int> > vecs= vsMake(h);
    1205   return vecs;
    1206 }
    1207 
    1208 /***********************************************************/
    1209 
    1210 
    1211 
    1212 
    1213 //returns the link of face a in simplicial complex X
     1520
     1521
     1522
     1523
     1524
     1525/************************links***********************************/
     1526
     1527
     1528//returns the links of face a in simplicial complex X
    12141529std::vector<std::vector<int> > links(poly a, ideal h)
    12151530{
    … …  
    12201535  {
    12211536    U=vecUnion(av,X[i]);
    1222     //PrintS("**********************\n");
    1223     //listprint(X[i]);
    1224     //PrintS("The union of them is:\n");
    1225     //listprint(U);
    12261537    In=vecIntersection(av,X[i]);
    1227     //PrintS("The intersection of them is:\n");
    1228     //listprint(In);
    12291538    if( In.size()==0 && vInvsl(U,X))
    12301539    {
    … …  
    12421551
    12431552
    1244 
    1245 /*vector version
    1246 //returns the link of face a in simplicial complex X
    1247 std::vector<std::vector<int> > links(std::vector<int> a, std::vector<std::vector<int> > X)
    1248 {
    1249   int i;
    1250   std::vector<std::vector<int> > lk;
    1251   std::vector<int> U;
    1252   std::vector<int> In;
    1253   for(i=0;i<X.size();i++)
    1254   {
    1255     U=vecUnion(a,X[i]);
    1256     //PrintS("**********************\n");
    1257     //listprint(X[i]);
    1258     //PrintS("The union of them is:\n");
    1259     //listprint(U);
    1260     In=vecIntersection(a,X[i]);
    1261     //PrintS("The intersection of them is:\n");
    1262     //listprint(In);
    1263     if( In.size()==0 && vInvsl(U,X))
    1264     {
    1265       //PrintS("The union of them is FACE and intersection is EMPTY!\n");
    1266       lk.push_back(X[i]);
     1553int redefinedeg(poly p, int  num)
     1554{
     1555  int deg=0, deg0;
     1556  for(int i=1;i<=currRing->N;i++)
     1557  {
     1558    deg0=pGetExp(p, i);
     1559    if(i>num)
     1560    {
     1561      deg= deg+2*deg0;
    12671562    }
    12681563    else
    12691564    {
    1270       ;
    1271     }
    1272   }
    1273   return lk;
    1274 }
    1275 */
     1565      deg=deg+deg0;
     1566    }
     1567  }
     1568  //Print("the new degree is: %d\n", deg);
     1569  return (deg);
     1570}
     1571
     1572
     1573// the degree of variables should be same
     1574ideal p_a(ideal h)
     1575{
     1576  poly e=pOne(), p;
     1577  int i,j,deg=0,deg0;
     1578  ideal aset=idCopy(h),ia,h1=idsrRing(h);
     1579//PrintS("idsrRing is:\n");id_print(h1);
     1580  std::vector<int> as;
     1581  std::vector<std::vector<int> > hvs=supports(h);
     1582  for(i=0;i<IDELEMS(h1);i++)
     1583  {
     1584    deg0=pTotaldegree(h1->m[i]);
     1585    if(deg < deg0)
     1586      deg=deg0;
     1587  }
     1588  for(i=2;i<=deg;i++)
     1589  {
     1590    ia=id_MaxIdeal(i, currRing);
     1591    for(j=0;j<IDELEMS(ia);j++)
     1592    {
     1593      p=pCopy(ia->m[j]);
     1594      if(!IsInX(p,h))
     1595      {
     1596        as=support1(p);
     1597        if(vInvsl(as,hvs))
     1598        {
     1599          idInsertPoly(aset, p);
     1600        }
     1601      }
     1602    }
     1603  }
     1604  idSkipZeroes(aset);
     1605  return(aset);
     1606}
     1607
     1608
     1609/*only for the exampels whose variables has degree more than 1*/
     1610/*ideal p_a(ideal h)
     1611{
     1612  poly e=pOne(), p;
     1613  int i,j,deg=0,deg0, ord=4;
     1614  ideal aset=idCopy(h),ia,h1=idsrRing(h);
     1615//PrintS("idsrRing is:\n");id_print(h1);
     1616  std::vector<int> as;
     1617  std::vector<std::vector<int> > hvs=supports(h);
     1618  for(i=0;i<IDELEMS(h1);i++)
     1619  {
     1620    deg0=redefinedeg(h1->m[i],ord);
     1621    if(deg < deg0)
     1622      deg=deg0;
     1623  }
     1624  for(i=2;i<=deg;i++)
     1625  {
     1626    ia=id_MaxIdeal(i, currRing);
     1627    for(j=0;j<IDELEMS(ia);j++)
     1628    {
     1629      p=pCopy(ia->m[j]);
     1630      if(!IsInX(p,h))
     1631      {
     1632        as=support1(p);
     1633        if(vInvsl(as,hvs))
     1634        {
     1635          idInsertPoly(aset, p);
     1636        }
     1637      }
     1638    }
     1639  }
     1640  idSkipZeroes(aset);
     1641  return(aset);
     1642}*/
     1643
     1644
     1645
     1646
     1647std::vector<std::vector<int> > id_subsets(std::vector<std::vector<int> > vecs)
     1648{
     1649  int i,j;
     1650  std::vector<std::vector<int> > vvs, res;
     1651  for(i=0;i<vecs.size();i++)
     1652  {
     1653    vvs=b_subsets(vecs[i]);
     1654    //listsprint(vvs);
     1655    for(j=0;j<vvs.size();j++)
     1656    {
     1657      if(!vInvsl(vvs[j],res))
     1658        res.push_back(vvs[j]);
     1659    }
     1660  }
     1661  //listsprint(res);
     1662  return (res);
     1663}
     1664
     1665
     1666
     1667
     1668std::vector<int> vertset(std::vector<std::vector<int> > vecs)
     1669{
     1670  int i,j;
     1671  std::vector<int> vert;
     1672  std::vector<std::vector<int> > vvs;
     1673  for(i=1;i<=currRing->N;i++)
     1674  {
     1675    for(j=0;j<vecs.size();j++)
     1676    {
     1677      if(IsinL(i, vecs[j]))
     1678      {
     1679        if(!IsinL(i , vert))
     1680        {
     1681          vert.push_back(i);
     1682        }
     1683        break;
     1684      }
     1685    }
     1686  }
     1687  return (vert);
     1688}
     1689
     1690//smarter way
     1691ideal p_b(ideal h, poly a)
     1692{
     1693  std::vector<std::vector<int> > pbv,lk=links(a,h), res;
     1694  std::vector<int> vert=vertset(lk), bv;
     1695  res=b_subsets(vert);
     1696  int i, j, nu=res.size(), adg=pTotaldegree(a);
     1697  poly e=pOne();
     1698  ideal idd=idInit(1,1);
     1699  for(i=0;i<res.size();i++)
     1700  {
     1701    if(res[i].size()==adg)
     1702      pbv.push_back(res[i]);
     1703  }
     1704  if(pEqualPolys(a,e))
     1705  {
     1706    idInsertPoly(idd, e);
     1707    idSkipZeroes(idd);
     1708    return (idd);
     1709  }
     1710  idd=idMaken(pbv);
     1711  return(idd);
     1712}
     1713
     1714/*//dump way to get pb
     1715// the degree of variables should be same
     1716ideal p_b(ideal h, poly a)
     1717{
     1718  std::vector<std::vector<int> > pbv,lk=links(a,h),res;
     1719// PrintS("Its links are :\n");id_print(idMaken(lk));
     1720  res=id_subsets(lk);
     1721  //PrintS("res is :\n");listsprint(res);
     1722  std::vector<int> bv;
     1723  ideal bset=findb(h);
     1724  int i,j,nu=res.size(),adg=pTotaldegree(a);
     1725  poly e=pOne();ideal idd=idInit(1,1);
     1726  for(i=0;i<res.size();i++)
     1727  {
     1728    if(res[i].size()==adg)
     1729      pbv.push_back(res[i]);
     1730  }
     1731  if(pEqualPolys(a,e)){idInsertPoly(idd, e); idSkipZeroes(idd); return (idd);}
     1732  for(i=0;i<nu;i++)
     1733  {
     1734    for(j=i+1;j<nu;j++)
     1735    {
     1736      if(res[i].size()!=0 && res[j].size()!=0)
     1737      {
     1738        bv = vecUnion(res[i], res[j]);
     1739        if(IsInX(pMaken(bv),bset)  && bv.size()==adg && !vInvsl(bv,pbv))
     1740          {pbv.push_back(bv);}
     1741      }
     1742    }
     1743  }
     1744  idd=idMaken(pbv);
     1745  //id_print(idd);
     1746  return(idd);
     1747}*/
     1748
     1749// also only for the examples whose variables have degree more than 1(ndegreeb and p_b)
     1750/*int ndegreeb(std::vector<int> vec, int num)
     1751{
     1752  int deg, deg0=0;
     1753  for(int i=0;i<vec.size();i++)
     1754  {
     1755    if(vec[i]>num)
     1756    {
     1757      deg0++;
     1758    }
     1759  }
     1760  deg=vec.size()+deg0;
     1761  return(deg);
     1762}
     1763
     1764ideal p_b(ideal h, poly a)
     1765{
     1766  std::vector<std::vector<int> > pbv,lk=links(a,h),res;
     1767// PrintS("Its links are :\n");id_print(idMaken(lk));
     1768  res=id_subsets(lk);
     1769  //PrintS("res is :\n");listsprint(res);
     1770  std::vector<int> bv;
     1771  ideal bset=findb(h);
     1772  int i,j,nu=res.size(),ord=4,adg=redefinedeg(a, ord);
     1773  poly e=pOne();ideal idd=idInit(1,1);
     1774  for(i=0;i<res.size();i++)
     1775  {
     1776    if(ndegreeb(res[i],ord)==adg)
     1777      pbv.push_back(res[i]);
     1778  }
     1779  if(pEqualPolys(a,e)){idInsertPoly(idd, e); idSkipZeroes(idd); return (idd);}
     1780  for(i=0;i<nu;i++)
     1781  {
     1782    for(j=i+1;j<nu;j++)
     1783    {
     1784      if(res[i].size()!=0 && res[j].size()!=0)
     1785      {
     1786        bv = vecUnion(res[i], res[j]);
     1787  //PrintS("bv is :\n");listprint(bv);
     1788 //Print("bv's degree is : %d\n", ndegreeb(bv,ord));
     1789        if(IsInX(pMaken(bv),bset)  && ndegreeb(bv,ord)==adg && !vInvsl(bv,pbv))
     1790        {
     1791          pbv.push_back(bv);
     1792        }
     1793      }
     1794    }
     1795  }
     1796  idd=idMaken(pbv);
     1797  //id_print(idd);
     1798  return(idd);
     1799}*/
     1800
    12761801
    12771802
    … …  
    12811806ideal psubset(poly p)
    12821807{
    1283   int i,j;
    1284   ideal h1,mons;
    1285   int max=pTotaldegree(p);
    1286   ideal id_re=idInit(1,1);
     1808  int i,j,max=pTotaldegree(p);
     1809  ideal h1,mons, id_re=idInit(1,1);
    12871810  for(i=1;i<max;i++)
    12881811  {
    … …  
    13371860poly pMake3(std::vector<int> vbase)
    13381861{
    1339   int n=vbase.size();poly p,q=0;
    1340   int co=1;
    1341   //equmab(n);
     1862  int n=vbase.size(),co=1;
     1863  poly p,q=0;
    13421864  for(int i=0;i<3;i++)
    13431865  {
    … …  
    13461868      if(i==1) co=-1;
    13471869      p = pOne();pSetExp(p, vbase[i], 1);pSetm(p);pSetCoeff(p, nInit(co));
    1348       //Print("attention! ");pWrite(p);
    13491870    }
    13501871    else p=0;
    … …  
    13521873    co=1;
    13531874  }
    1354   //PrintS("the polynomial according to the metrix M is:\n");
    1355   //listprint(vbase);
    1356   //pWrite(q);
    13571875  return q;
    13581876}
    … …  
    13661884  for(i=0;i<lv;i++)
    13671885  {
    1368     //Print("The vector is:  ");
    1369     //listprint(vecs[i]);
    13701886    p=pMake3(vecs[i]);
    1371     //Print("The polynomial is:  ");
    1372     //pWrite(p);
    13731887    idInsertPoly(id_re, p);
    13741888  }
    1375   //PrintS("This is the metrix M:\n");
    1376   //listsprint(vecs);
    1377   //PrintS("the ideal according to metrix M is:\n");
    13781889  idSkipZeroes(id_re);
    13791890  return id_re;
    … …  
    13861897{
    13871898  int i,j;
    1388   //Print("There are %d variables in total.\n",num);
     1899  //Print("There are %d new variables for equations solving.\n",num);
    13891900  ring r=currRing;
    13901901  char** tt;
    … …  
    13961907    sprintf (tt[i], "t(%d)", i+1);
    13971908    tt[i]=omStrDup(tt[i]);
    1398     //Print("%s\n",tt[i]);
    13991909  }
    14001910  ring R=rDefault(cf,num,tt,ringorder_lp);
    … …  
    14091919std::vector<int> subspace1(std::vector<std::vector<int> > mv, std::vector<int> bv)
    14101920{
    1411   int i;
    1412   int num=mv.size();
     1921  int i, num=mv.size();
    14131922  std::vector<int> base;
    1414   std::vector<int> pv;
    14151923  for(i=0;i<num;i++)
    14161924  {
    … …  
    14281936
    14291937
    1430 /****************************************************************/
    1431 //construct a monomial based on the support of it
    1432 //returns a squarefree monomial
    1433 poly pMaken(std::vector<int> vbase)
    1434 {
    1435   int n=vbase.size();
    1436   poly p,q=pOne();
    1437   //equmab(n);
    1438   for(int i=0;i<n;i++)
    1439   {
    1440       p = pOne();pSetExp(p, vbase[i], 1);pSetm(p);pSetCoeff(p, nInit(1));
    1441       //Print("attention! ");pWrite(p);
    1442       q=pp_Mult_mm(q,p,currRing);
    1443   }
    1444   //PrintS("the polynomial according to the metrix M is:\n");
    1445   //listprint(vbase);
    1446   //pWrite(q);
    1447   return q;
    1448 }
    1449 
    1450 // returns a ideal according to a set of supports
    1451 ideal idMaken(std::vector<std::vector<int> > vecs)
    1452 {
    1453   ideal id_re=idInit(1,1);
    1454   poly p;
    1455   int i,lv=vecs.size();
    1456   for(i=0;i<lv;i++)
    1457   {
    1458     //Print("The vector is:  ");
    1459     //listprint(vecs[i]);
    1460 
    1461     p=pMaken(vecs[i]);
    1462     //Print("The polynomial is:  ");
    1463     //pWrite(p);
    1464     idInsertPoly(id_re, p);
    1465   }
    1466   //PrintS("This is the metrix M:\n");
    1467   //listsprint(vecs);
    1468   //PrintS("the ideal according to metrix M is:\n");
    1469   idSkipZeroes(id_re);
    1470   //id_print(id_re);
    1471   return id_re;
    1472 }
    1473 
    1474 
    14751938
    14761939
    … …  
    14931956}
    14941957
     1958
     1959
    14951960// returns a ideal according to a set of supports
    14961961 std::vector<std::vector<poly> > idMakei(std::vector<std::vector<int> > mv,std::vector<std::vector<int> > vecs)
    … …  
    15201985
    15211986//return the graded pieces of cohomology T^1 according to a,b
     1987//original method (only for debugging)
    15221988void gradedpiece1(ideal h,poly a,poly b)
    15231989{
    1524   int i,j;
    1525   std::vector<std::vector<int> > hvs=supports(h);
    1526   std::vector<int> av=support1(a);
    1527   std::vector<int> bv=support1(b);
     1990  int i,j,m;
    15281991  ideal sub=psubset(b);
    1529   std::vector<std::vector<int> > sbv=supports(sub);
    1530   //ideal M=Mab(h,a,b);
    1531   std::vector<std::vector<int> > mv=Mabv(h,a,b);
    1532   PrintS("The homophisim is map onto the set:\n");
    1533   id_print(idMaken(mv));
    1534   int m=mv.size();
    1535   std::vector<std::vector<int> > good;
    1536   std::vector<int> bad,vv;
     1992  std::vector<int> av=support1(a), bv=support1(b), bad, vv;
     1993  std::vector<std::vector<int> > hvs=supports(h), sbv=supports(sub), mv=Mabv(h,a,b),good;
     1994  m=mv.size();
    15371995  ring r=currRing;
    15381996  if( m > 0 )
    … …  
    15642022    if(bv.size()!=1)
    15652023    {
    1566       PrintS("This is the solution of coefficients:\n");
     2024      //PrintS("This is the solution of coefficients:\n");
    15672025      listsprint(solve);
    15682026    }
    … …  
    15762034      equmab(solve[0].size());
    15772035      std::vector<std::vector<int> > solves=vecqring(solve,suu);
    1578       PrintS("This is the solution of coefficients:\n");
     2036      //PrintS("This is the solution of coefficients:\n");
    15792037      listsprint(solves);
    15802038      rChangeCurrRing(r);
    … …  
    16972155std::vector<int> numfree(ideal h)
    16982156{
    1699   int i,j,num=0;poly p;
     2157  int i,j,num=0;
    17002158  std::vector<int> fvar;
    17012159  for(j=1;j<=currRing->N;j++)
    … …  
    17512209    ideal h1=getpresolve(h2);
    17522210    poly q,e=pOne();
    1753     int lg=IDELEMS(h1);
     2211    int lg=IDELEMS(h1),n,i,j,t;
    17542212    std::vector<int> fvar=numfree(h1);
    1755     int n=fvar.size();
    1756     //Print("*************&&&&&&&&&&&*******************There are %d free variables in total\n",n);
    1757     int i,j,t;
    1758     //Print("lg is %d\n",lg);
     2213    n=fvar.size();
    17592214    if(n==0)
    17602215    {
    … …  
    17682223        for(i=0;i<lg;i++)
    17692224        {
    1770           //Print("The polynomial is the %dth one:\n",i+1);
    17712225          q=pCopy(h1->m[i]);
    17722226          //pWrite(q);
    … …  
    17922246            else
    17932247            {
    1794               //Print("aiyamaya is %d \n",n_Int(pGetCoeff(q),currRing->cf));
    17952248              vec.push_back(n_Int(pGetCoeff(q),currRing->cf));
    17962249            }
    … …  
    18452298std::vector<int> subspacet1(int num, std::vector<std::vector<int> > ntvs)
    18462299{
    1847   int i,j,t;
    1848   int n=ntvs.size();
     2300  int i, j, t, n=ntvs.size();
    18492301  std::vector<int> subase;
    18502302  for(t=0;t<n;t++)
    … …  
    18772329{
    18782330  int i,j;
    1879   std::vector<int> alset=findalpha(mv,bv);
    1880   std::vector<int> subase;
     2331  std::vector<int> alset=findalpha(mv,bv), subase;
    18812332  std::vector<std::vector<int> > subases;
    18822333  for(i=0;i<alset.size();i++)
    … …  
    19192370
    19202371//fix the problem of the number of the new variables
     2372//original method for T^2(only for debugging)
    19212373void gradedpiece2(ideal h,poly a,poly b)
    19222374{
    1923   int t0,t1,t2,i,j,t;
     2375  int t0,t1,t2,i,j,t,m;
     2376  ideal sub=psubset(b);
    19242377  ring r=rCopy(currRing);
    1925   std::vector<std::vector<int> > hvs=supports(h);
    1926   std::vector<int> av=support1(a);
    1927   std::vector<int> bv=support1(b);
    1928   ideal sub=psubset(b);
    1929   std::vector<std::vector<int> > mv=Mabv(h,a,b);
    1930   std::vector<std::vector<int> > mts=mabtv(hvs,mv,av,bv);
     2378  std::vector<std::vector<int> > hvs=supports(h), mv=Mabv(h,a,b), mts, vecs,vars;
     2379  std::vector<int> av=support1(a), bv=support1(b), vec,var;
     2380  mts=mabtv(hvs,mv,av,bv);
    19312381  PrintS("The homomorphism should map onto:\n");
    19322382  lpsprint(idMakei(mv,mts));
    1933   int m=mv.size();
    1934   //ideal M=Mab(h,a,b);
    1935   std::vector<std::vector<int> > vecs,vars;
    1936   std::vector<int> vec,var;
     2383  m=mv.size();
    19372384  if(m > 0)
    19382385  {
    … …  
    20312478
    20322479
    2033 
    2034 
    2035 
    2036 
    2037 
    2038 
    2039 
    2040 
    2041 
    2042 
    2043 
    2044 
    2045 
    2046 
    20472480/**********************************************************************/
    20482481//For the method of N_{a-b}
    … …  
    20542487bool nabconditionv(std::vector<std::vector<int> > hvs, std::vector<int> pv,  std::vector<int> av,  std::vector<int> bv)
    20552488{
    2056   std::vector<int> vec1=vecIntersection(pv,bv);
    2057   std::vector<int> vec2=vecUnion(pv,bv);
     2489  std::vector<int> vec1=vecIntersection(pv,bv), vec2=vecUnion(pv,bv);
    20582490  int s1=vec1.size();
    20592491  if(!vInvsl(vec2,hvs) && s1==0 && vsubset(av,pv))
    … …  
    20792511  {
    20802512    if(nabconditionv(hvs,hvs[i],av,bv))
     2513    {
     2514      //PrintS("satisfy:\n");
    20812515      vecs.push_back(hvs[i]);
     2516    }
    20822517  }
    20832518  return vecs;
    … …  
    20972532  if(vInvsl(v1,hvs))
    20982533  {
    2099     //PrintS("They are in Nab2\n");
    21002534    return (true);
    21012535  }
    2102   //PrintS("They are not in Nab2\n");
    21032536  return (false);
    21042537}
    … …  
    21462579    if(!vInvsl(sv,hvs))
    21472580    {
    2148       //PrintS("TRUE! It is in tilde Nab!\n");
    21492581      return true;
    21502582    }
    21512583  }
    2152   //PrintS("FALSE! It is not in tilde Nab!\n");
    21532584  return false;
    21542585}
    … …  
    21622593std::vector<int>  tnab(std::vector<std::vector<int> > hvs,std::vector<std::vector<int> > nvs,std::vector<std::vector<int> > bvs)
    21632594{
    2164   std::vector<int> pv;
    2165   std::vector<int> vec;
     2595  std::vector<int> pv, vec;
    21662596  for(int j=0;j<nvs.size();j++)
    21672597  {
    … …  
    21962626std::vector<std::vector<int> > value1(std::vector<std::vector<int> > mvs, std::vector<std::vector<int> > nvs, std::vector<std::vector<int> > vecs,std::vector<int> av, std::vector<int> bv)
    21972627{
    2198   std::vector<int> pv;
    2199   std::vector<int> base;
     2628  int j;
     2629  std::vector<int> pv, base;
    22002630  std::vector<std::vector<int> > bases;
    22012631  for(int t=0;t<vecs.size();t++)
    … …  
    22042634    {
    22052635      pv=phimage(mvs[i],av,bv);
    2206       for(int j=0;j<nvs.size();j++)
     2636      for( j=0;j<nvs.size();j++)
    22072637      {
    22082638        if(vEvl(pv,nvs[j]))
     2639        {
    22092640          base.push_back(vecs[t][j]);
     2641          break;
     2642        }
    22102643      }
     2644      if(j==nvs.size())
     2645      {
     2646        base.push_back(0);
     2647      }
    22112648    }
    22122649    if(base.size()!=mvs.size())
    22132650    {
    2214       WerrorS("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1");
     2651      //WerrorS("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1");
     2652      WerrorS("Errors in Equations solving (Values Finding)!");
    22152653      usleep(1000000);
    22162654      assert(false);
    2217       WerrorS("Errors in Nab set!");
     2655
    22182656    }
    22192657
    … …  
    22382676 //vtm(solve);
    22392677  intvec *m;
    2240   int i,j;
    2241   int a=vecs.size();
     2678  int i,j, a=vecs.size();
    22422679  if(a==0)
    22432680  {
    … …  
    22642701
    22652702
    2266 /*void ppppp(int l)
    2267 {
    2268   int i;
    2269 std::vector<std::vector<int> > vecs;
    2270 std::vector<int> vec;
    2271   for(i=0;i<l;i++)
    2272   {
    2273     for(int j=i*l;j<i*l+l;j++)
    2274     {
    2275       vec.push_back(j);
    2276     }
    2277     vecs.push_back(vec);
    2278     vec.clear();
    2279   }
    2280 PrintS("May I have your attention please!\n");
    2281 listsprint(vecs);
    2282     intvec *m=T1mat(vecs);
    2283 PrintS("May I have your attention again!\n");
    2284   m->show(0,0);
    2285 }*/
    2286 
    2287 
    2288 
    2289 
    2290 
     2703
     2704
     2705
     2706//returns the set of position number of minimal gens in M
    22912707std::vector<int> gensindex(ideal M, ideal ids)
    22922708{
    … …  
    23112727{
    23122728  int i;
    2313   ideal hi=idInit(1,1);
    23142729  std::vector<std::vector<int> > mv=Mabv(h,a,b);
    2315   ideal M=idMaken(mv);
    2316 //PrintS("mab\n");
    2317   //id_print(M);
    2318   //PrintS("idsrRing\n");
    2319   //id_print(idsrRing(h));
     2730  ideal M=idMaken(mv), hi=idInit(1,1);
    23202731  std::vector<int> index = gensindex(M, idsrRing(h));
    2321 //PrintS("index\n");
    2322   //listprint(index);
    23232732  for(i=0;i<index.size();i++)
    23242733  {
    … …  
    23262735  }
    23272736  idSkipZeroes(hi);
    2328 //PrintS("over\n");
    2329  // id_print(hi);
    23302737  return (hi);
    23312738}
    … …  
    23362743{
    23372744  int i,j;
    2338   std::vector<int> vec;
     2745  std::vector<int> vec,solm;
    23392746  std::vector<std::vector<int> > solsm;
    2340   std::vector<int> solm;
    23412747  for(i=0;i<solve.size();i++)
    23422748  {
    … …  
    23552761
    23562762//T_1 graded piece(N method)
     2763//frame of the most efficient version
     2764//regardless of links
    23572765
    23582766intvec * gradedpiece1n(ideal h,poly a,poly b)
    23592767{
    2360   int i,j,co;
    2361   std::vector<std::vector<int> > hvs=supports(h);
    2362   std::vector<int> av=support1(a);
    2363   //listprint(av);
    2364   std::vector<int> bv=support1(b);
    2365   //listprint(bv);
    2366   ideal sub=psubset(b);
    2367 //id_print(sub);
    2368   std::vector<std::vector<int> > sbv=supports(sub);
    2369 //listsprint(sbv);
    2370   std::vector<std::vector<int> > nv=Nabv(hvs,av,bv);
    2371   //PrintS("The N set is:\n");
    2372   //listsprint(nv);
    2373   std::vector<std::vector<int> > mv=Mabv(h,a,b);
    2374   //listsprint(mv);
    2375   ideal M=idMaken(mv);
    2376   std::vector<int> index = gensindex(M, idsrRing(h));
    2377   //ideal gens=mingens(M,index);
    2378   int n=nv.size();
    2379   //PrintS("The homophisim is map onto the set:\n");
    2380   //id_print(M);
    2381   std::vector<std::vector<int> > good,solve;
    2382   std::vector<int> bad;
     2768  int i,j,co,n;
     2769  std::vector<std::vector<int> > hvs=supports(h),mv=Mabv(h,a,b),sbv,nv,good,solve;
     2770  std::vector<int> av=support1(a), bv=support1(b), bad, tnv, index;
     2771  ideal sub=psubset(b),M;
     2772  sbv=supports(sub);
     2773  nv=Nabv(hvs,av,bv);
     2774  M=idMaken(mv);
     2775  index = gensindex(M, idsrRing(h));
     2776  n=nv.size();
    23832777  ring r=currRing;
    2384   std::vector<int> tnv;
    23852778  if(n > 0)
    23862779  {
    … …  
    24392832
    24402833
    2441 /*
    2442 void vtm(std::vector<std::vector<int> > vecs)
    2443 {
    2444   int i,j;
    2445 
    2446   intvec *m;
    2447 
    2448   int r=vecs.size();
    2449   Print("r is %d\n",r);
    2450  // c=(vecs[0]).size();
    2451 
    2452   PrintS("no error yet:\n");
    2453   m=new intvec(r);
    2454 
    2455 for(i=0;i<r;i++)
    2456   {
    2457 
    2458     for(j=0;j<r;j++)
    2459     {
    2460    (*m)[i]=*(vecs[i]);
    2461        Print("%d",IMATELEM(*m,i,j));
    2462     }
    2463   }
    2464  // return matt;
    2465 }
    2466 */
    2467 
    2468 
    2469 
    2470 
    2471 
     2834
     2835
     2836
     2837//for debugging
    24722838void T1(ideal h)
    24732839{
    … …  
    24792845  for(int i=0;i<IDELEMS(bi);i++)
    24802846  {
    2481     PrintS("This is aset according to:");
     2847    //PrintS("This is aset according to:");
    24822848    b=pCopy(bi->m[i]);
    24832849    pWrite(b);
    … …  
    25472913{
    25482914  int i,j;
    2549   std::vector<int> alset=findalphan(N,tN);
    2550   std::vector<int> subase;
     2915  std::vector<int> alset=findalphan(N,tN), subase;
    25512916  std::vector<std::vector<int> > subases;
    25522917  for(i=0;i<alset.size();i++)
    … …  
    25682933std::vector<std::vector<int> > value2(std::vector<std::vector<int> > mvs, std::vector<std::vector<int> > nvs, std::vector<std::vector<int> > mts, std::vector<std::vector<int> > nts, std::vector<std::vector<int> > vecs,std::vector<int> av,   std::vector<int> bv)
    25692934{
    2570   std::vector<int> pv,qv;
    2571   std::vector<int> base;
    2572   int row,col;
     2935  int row,col,j;
     2936  std::vector<int> pv,qv, base;
    25732937  std::vector<std::vector<int> > bases;
    25742938  //PrintS("This is the nabt:\n");
    … …  
    25782942  //listsprint(mts);
    25792943  //PrintS("mabt ends:\n");
    2580 
    25812944  for(int t=0;t<vecs.size();t++)
    25822945  {
    … …  
    25892952      if(vEvl(pv,qv))
    25902953        base.push_back(0);
    2591       //PrintS("This is image of p and q:\n");
    2592       //listprint(pv);  PrintS("*********************\n");listprint(qv);
    2593       //PrintS("nabt ends:\n");
    25942954      else
    25952955      {
    2596         for(int j=0;j<nts.size();j++)
     2956        for(j=0;j<nts.size();j++)
    25972957        {
    25982958          row=nts[j][0];
    25992959          col=nts[j][1];
    2600           //PrintS("This is nvs:\n");
    2601           //listprint(nvs[row]); PrintS("*********************\n");listprint(nvs[col]);
    2602           //PrintS("nabt ends:\n");
    26032960          if(vEvl(pv,nvs[row])&&vEvl(qv,nvs[col]))
    26042961          {
    26052962            base.push_back(vecs[t][j]);break;
    2606             //PrintS("This is nvs,they are the same:\n");
    2607             //listprint(nvs[row]); listprint(nvs[col]);
    2608             //PrintS("nabt ends:\n");
    26092963          }
    2610           else
     2964          else if(vEvl(pv,nvs[col])&&vEvl(qv,nvs[row]))
    26112965          {
    26122966            base.push_back(-vecs[t][j]);break;
    2613             //PrintS("This is nvs,they are the same:\n");
    2614             //listprint(nvs[row]); listprint(nvs[col]);
    2615             //PrintS("nabt ends:\n");
    26162967          }
    26172968        }
     2969        if(j==nts.size()) {base.push_back(0);}
    26182970      }
    26192971    }
    26202972    if(base.size()!=mts.size())
    26212973    {
    2622        WerrorS("Errors in Nab set!");
    2623         //WerrorS("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1");
     2974      WerrorS("Errors in Values Finding(value2)!");
     2975       //WerrorS("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1");
    26242976      usleep(1000000);
    26252977      assert(false);
    … …  
    26372989{
    26382990  int i,j;
    2639   std::vector<std::vector<int> > hvs=supports(h);
     2991  std::vector<std::vector<int> > hvs=supports(h),mv,mts;
    26402992  std::vector<int> av=support1(a), bv=support1(b);
    2641   std::vector<std::vector<int> > mv=Mabv(h,a,b), mts=mabtv(hvs,mv,av,bv);
     2993  mv=Mabv(h,a,b);
     2994  mts=mabtv(hvs,mv,av,bv);
    26422995  std::vector<std::vector<poly> > pvs=idMakei(mv,mts);
    26432996  ideal gens=idInit(1,1);
    … …  
    26483001  }
    26493002  idSkipZeroes(gens);
    2650 //PrintS("This is the mix set of mab2!\n");
    2651 //id_print(gens);
    26523003  return (gens);
    26533004}
    … …  
    26623013intvec * gradedpiece2n(ideal h,poly a,poly b)
    26633014{
    2664   int i,j,t;
    2665   std::vector<std::vector<int> > hvs=supports(h);
    2666   std::vector<int> av=support1(a);
    2667   std::vector<int> bv=support1(b);
     3015  int i,j,t,n;
     3016  std::vector<std::vector<int> > hvs=supports(h),nv,mv,mts,sbv,vecs,vars,ntvs,solve;
     3017  std::vector<int> av=support1(a), bv=support1(b),tnv,vec,var;
    26683018  ideal sub=psubset(b);
    2669   std::vector<std::vector<int> > sbv=supports(sub);
    2670   std::vector<std::vector<int> > nv=Nabv(hvs,av,bv);
    2671   int n=nv.size();
    2672   std::vector<int> tnv=tnab(hvs,nv,sbv);
     3019  sbv=supports(sub);
     3020  nv=Nabv(hvs,av,bv);
     3021  n=nv.size();
     3022  tnv=tnab(hvs,nv,sbv);
    26733023  ring r=currRing;
    2674   std::vector<std::vector<int> > mv=Mabv(h,a,b);
    2675   std::vector<std::vector<int> > mts=mabtv(hvs,mv,av,bv);
     3024  mv=Mabv(h,a,b);
     3025  mts=mabtv(hvs,mv,av,bv);
    26763026  //PrintS("The relations are:\n");
    26773027  //listsprint(mts);
    26783028  //PrintS("The homomorphism should map onto:\n");
    26793029  //lpsprint(idMakei(mv,mts));
    2680   std::vector<std::vector<int> > vecs,vars,ntvs;
    2681   std::vector<int> vec,var;
    2682   std::vector<std::vector<int> > solve;
    26833030  if(n>0)
    26843031  {
    … …  
    26933040      if(tNab(hvs,nv[i],sbv)&&tNab(hvs,nv[j],sbv))//condition 1
    26943041      {
    2695         //pWrite(pMaken(nv[i]));pWrite(pMaken(nv[j]));
    2696         //PrintS("They are both in tilde N.\n");
    2697         //PrintS("tilde N is:\n");  listprint(tnv);
    26983042        vec=makeequation(t0+1,0,0);
    26993043        vecs.push_back(vec);
    … …  
    27503094
    27513095
    2752 
     3096//for debugging
    27533097void T2(ideal h)
    27543098{
    … …  
    27943138
    27953139
    2796 
     3140/*****************************for links*******************************************/
     3141//the image phi(pv)=pv minus av minus bv
     3142std::vector<int> phimagel(std::vector<int> fv,  std::vector<int> av, std::vector<int> bv)
     3143{
     3144  std::vector<int> nv;
     3145  nv=vecMinus(fv,bv);
     3146  nv=vecMinus(nv,av);
     3147  return nv;
     3148}
     3149
     3150
     3151
     3152//mvs and nvs are the supports of ideal Mab and Nab
     3153//vecs is the solution of nab
     3154std::vector<std::vector<int> > value1l(std::vector<std::vector<int> > mvs, std::vector<std::vector<int> > lks, std::vector<std::vector<int> > vecs,std::vector<int> av, std::vector<int> bv)
     3155{
     3156  int j;
     3157  std::vector<int> pv;
     3158  std::vector<int> base;
     3159  std::vector<std::vector<int> > bases;
     3160  for(int t=0;t<vecs.size();t++)
     3161  {
     3162    for(int i=0;i<mvs.size();i++)
     3163    {
     3164      pv=phimagel(mvs[i], av, bv);
     3165      for(j=0;j<lks.size();j++)
     3166      {
     3167        if(vEvl(pv,lks[j]))
     3168        {
     3169          base.push_back(vecs[t][j]);break;
     3170        }
     3171      }
     3172      //if(j==lks.size()) {base.push_back(0);}
     3173    }
     3174    if(base.size()!=mvs.size())
     3175    {
     3176      WerrorS("Errors in Values Finding(value1l)!");
     3177      usleep(1000000);
     3178      assert(false);
     3179
     3180    }
     3181
     3182    bases.push_back(base);
     3183    base.clear();
     3184  }
     3185  return bases;
     3186}
     3187
     3188/***************************************************/
     3189clock_t t_begin, t_mark, t_start, t_construct=0, t_solve=0, t_value=0, t_total=0;
     3190/**************************************************/
     3191
     3192
     3193static void TimeShow(clock_t t_construct, clock_t t_solve, clock_t t_value ,clock_t t_total)
     3194{
     3195  Print("The time of value matching for first order deformation:   %.2f sec ;\n", ((double) t_value)/CLOCKS_PER_SEC);
     3196  Print("The total time of fpiece:  %.2f sec ;\n", ((double) t_total)/CLOCKS_PER_SEC);
     3197  Print("The time of equations construction for fpiece:   %.2f sec ;\n", ((double) t_construct)/CLOCKS_PER_SEC);
     3198  Print("The total time of equations solving for fpiece:  %.2f sec ;\n", ((double) t_solve)/CLOCKS_PER_SEC);
     3199  PrintS("__________________________________________________________\n");
     3200}
     3201
     3202
     3203
     3204std::vector<std::vector<int> > gpl(ideal h,poly a,poly b)
     3205{
     3206  int i,j,co;
     3207  std::vector<std::vector<int> > hvs=supports(h),sbv,nv,mv,good,solve;
     3208  std::vector<int> av=support1(a), bv=support1(b),index,bad,tnv;
     3209  ideal sub=psubset(b);
     3210  sbv=supports(sub);
     3211  nv=Nabv(hvs,av,bv);
     3212  mv=Mabv(h,a,b);
     3213  ideal M=idMaken(mv);
     3214  index = gensindex(M, idsrRing(h));
     3215  int n=nv.size();
     3216  ring r=currRing;
     3217  t_begin=clock();
     3218  if(n > 0)
     3219  {
     3220    tnv=tnab(hvs,nv,sbv);
     3221    for(i=0;i<tnv.size();i++)
     3222    {
     3223      co=tnv[i];
     3224      bad.push_back(co+1);
     3225    }
     3226    for(i=0;i<n;i++)
     3227    {
     3228      for(j=i+1;j<n;j++)
     3229      {
     3230        if(nabtconditionv(hvs,nv[i],nv[j],av,bv))
     3231        {
     3232          good=listsinsertlist(good,i+1,j+1);
     3233        }
     3234        else
     3235        {
     3236          ;
     3237        }
     3238      }
     3239    }
     3240    t_construct=t_construct+clock()-t_begin;
     3241    t_begin=clock();
     3242    solve=eli2(n,bad,good);
     3243    t_solve=t_solve+clock()-t_begin;
     3244    if(bv.size()!=1)
     3245    {;
     3246    }
     3247    else
     3248    {
     3249      std::vector<int> su=make1(n);
     3250      std::vector<std::vector<int> > suu;
     3251      suu.push_back(su);
     3252      equmab(n);
     3253      solve=vecqring(solve,suu);
     3254      rChangeCurrRing(r);
     3255    }
     3256  }
     3257  else
     3258  {
     3259    solve.clear();
     3260  }
     3261  //listsprint(solve);
     3262  //sl->show(0,0);
     3263  return solve;
     3264}
     3265
     3266
     3267//T^1
     3268//only need to consider the links of a, and reduce a to empty set
     3269intvec * gradedpiece1nl(ideal h,poly a,poly b, int set)
     3270{
     3271  t_start=clock();
     3272  int i,j,co;
     3273  poly e=pOne();
     3274  std::vector<int> av=support1(a),bv=support1(b),index, em;
     3275  std::vector<std::vector<int> > solve, hvs=supports(h), lks=links(a,h),  mv=Mabv(h,a,b), nvl;
     3276  ideal id_links=idMaken(lks);
     3277  ideal M=idMaken(mv);
     3278  index = gensindex(M, idsrRing(h));
     3279  solve=gpl(id_links,e,b);
     3280  t_mark=clock();
     3281  nvl=Nabv(lks,em,bv);
     3282  solve=value1l(mv, nvl , solve, av, bv);
     3283  if(set==1)
     3284  {
     3285    solve=minisolve(solve,index);
     3286  }
     3287  intvec *sl=Tmat(solve);
     3288  t_value=t_value+clock()-t_mark;
     3289  t_total=t_total+clock()-t_start;
     3290  return sl;
     3291}
     3292
     3293
     3294
     3295
     3296//for finding values of T^2
     3297std::vector<std::vector<int> > value2l(std::vector<std::vector<int> > mvs, std::vector<std::vector<int> > lks, std::vector<std::vector<int> > mts, std::vector<std::vector<int> > lkts, std::vector<std::vector<int> > vecs,std::vector<int> av,   std::vector<int> bv)
     3298{
     3299  std::vector<int> pv,qv,base;
     3300  int row,col,j;
     3301  std::vector<std::vector<int> > bases;
     3302  if(vecs.size()==0)
     3303  {
     3304
     3305  }
     3306  for(int t=0;t<vecs.size();t++)
     3307  {
     3308    for(int i=0;i<mts.size();i++)
     3309    {
     3310      row=mts[i][0];
     3311      col=mts[i][1];
     3312      pv=phimagel(mvs[row],av,bv);
     3313      qv=phimagel(mvs[col],av,bv);
     3314      if(vEvl(pv,qv))
     3315        base.push_back(0);
     3316      else
     3317      {
     3318        for(j=0;j<lkts.size();j++)
     3319        {
     3320          row=lkts[j][0];
     3321          col=lkts[j][1];
     3322          if(vEvl(pv,lks[row])&&vEvl(qv,lks[col]))
     3323          {
     3324            base.push_back(vecs[t][j]);break;
     3325          }
     3326          else if(vEvl(qv,lks[row])&&vEvl(pv,lks[col]))
     3327          {
     3328            base.push_back(-vecs[t][j]);break;
     3329          }
     3330        }
     3331        //if(j==lkts.size())
     3332        //{
     3333          //base.push_back(0);
     3334        //}
     3335      }
     3336    }
     3337    if(base.size()!=mts.size())
     3338    {
     3339       WerrorS("Errors in Values Finding!");
     3340        //WerrorS("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1");
     3341      usleep(1000000);
     3342      assert(false);
     3343    }
     3344    bases.push_back(base);
     3345    base.clear();
     3346  }
     3347  return bases;
     3348}
     3349
     3350
     3351std::vector<std::vector<int> > gpl2(ideal h,poly a,poly b)
     3352{
     3353  int i,j,t,n;
     3354  std::vector<std::vector<int> > hvs=supports(h),sbv,nv,mv,mts,vecs,vars,ntvs,solve;
     3355  std::vector<int> av=support1(a), bv=support1(b),vec,var,tnv;
     3356  ideal sub=psubset(b);
     3357  sbv=supports(sub);
     3358  nv=Nabv(hvs,av,bv);
     3359  n=nv.size();
     3360  tnv=tnab(hvs,nv,sbv);
     3361  ring r=currRing;
     3362  mv=Mabv(h,a,b);
     3363  mts=mabtv(hvs,mv,av,bv);
     3364  if(n>0)
     3365  {
     3366    ntvs=nabtv( hvs, nv, av, bv);
     3367    int l=ntvs.size();
     3368    if(l>0)
     3369    {
     3370      for(int t0=0;t0<l;t0++)
     3371      {
     3372        i=ntvs[t0][0];
     3373        j=ntvs[t0][1];
     3374        if(tNab(hvs,nv[i],sbv)&&tNab(hvs,nv[j],sbv))//condition 1
     3375        {
     3376          vec=makeequation(t0+1,0,0);
     3377          vecs.push_back(vec);
     3378          vec.clear();
     3379        }
     3380        for(int t1=t0+1;t1<ntvs.size();t1++)
     3381        {
     3382          for(int t2=t1+1;t2<ntvs.size();t2++)
     3383          {
     3384            if(ntvs[t0][0]==ntvs[t1][0]&&ntvs[t1][1]==ntvs[t2][1]&&ntvs[t0][1]==ntvs[t2][0])
     3385            {
     3386              i=ntvs[t0][0];
     3387              j=ntvs[t0][1];
     3388              t=ntvs[t1][1];
     3389              if(condition2for2nv(hvs,nv[i],nv[j],nv[t]))
     3390              {
     3391                vec=makeequation(t0+1,t1+1,t2+1);
     3392                vecs.push_back(vec);
     3393                vec.clear();
     3394              }
     3395            }
     3396          }
     3397        }
     3398      }
     3399      if(n==1) {l=1;}
     3400      equmab(l);
     3401      ideal id_re=idMake3(vecs);
     3402      std::vector<std::vector<int> > re=getvector(id_re,l);
     3403      rChangeCurrRing(r);
     3404      std::vector<std::vector<int> > sub=subspacetn(nv, tnv,ntvs);
     3405      equmab(l);
     3406      solve=vecqring(re, sub);
     3407      rChangeCurrRing(r);
     3408    }
     3409    else
     3410    {
     3411      solve.clear();
     3412    }
     3413  }
     3414  else
     3415    solve.clear();
     3416  return solve;
     3417}
     3418
     3419
     3420
     3421
     3422
     3423
     3424intvec * gradedpiece2nl(ideal h,poly a,poly b)
     3425{
     3426  int i,j,t;
     3427  poly e=pOne();
     3428  std::vector<int> av=support1(a), bv=support1(b), em;
     3429  std::vector<std::vector<int> > hvs=supports(h), mv=Mabv(h,a,b),mts,solve,lks,nvl,ntsl;
     3430  mts=mabtv(hvs,mv,av,bv);
     3431  lks=links(a,h);
     3432  ideal id_links=idMaken(lks);
     3433//PrintS("This is the links of a:\n"); id_print(id_links);
     3434  nvl=Nabv(lks,em,bv);
     3435//PrintS("This is the N set:\n"); id_print(idMaken(nvl));
     3436  ntsl=nabtv(lks,nvl,em,bv);
     3437//PrintS("This is N^2:\n"); listsprint(ntsl);
     3438  solve=gpl2(id_links,e,b);
     3439//PrintS("This is pre solution of N:\n"); listsprint(solve);
     3440  if(solve.size() > 0)
     3441  {
     3442    solve=value2l(mv, nvl, mts, ntsl, solve, av, bv);
     3443  }
     3444//PrintS("This is solution of N:\n"); listsprint(solve);
     3445  intvec *sl=Tmat(solve);
     3446  return sl;
     3447}
     3448
     3449
     3450
     3451//for debugging
     3452/*
    27973453void Tlink(ideal h,poly a,poly b,int n)
    27983454{
    … …  
    28203476    gradedpiece2n(li,p,b);
    28213477}
    2822 
    2823 
    2824 
    2825 
    2826 
    2827 
    2828 
    2829 
    2830 
    2831 
    2832 
    2833 
    2834 
    2835 
     3478*/
     3479
     3480
     3481
     3482/******************************for triangulation***********************************/
     3483
     3484
     3485
     3486//returns all the faces which are triangles
     3487ideal trisets(ideal h)
     3488{
     3489  int i;
     3490  ideal ids=idInit(1,1);
     3491  std::vector<int> pv;
     3492  for(i=0;i<IDELEMS(h);i++)
     3493  {
     3494    pv= support1(h->m[i]);
     3495    if(pv.size()==3)
     3496      idInsertPoly(ids, pCopy(h->m[i]));
     3497  }
     3498  idSkipZeroes(ids);
     3499  return ids;
     3500}
     3501
     3502
     3503
     3504
     3505// case 1 new faces
     3506std::vector<std::vector<int> > triface(poly p, int vert)
     3507{
     3508  int i;
     3509  std::vector<int> vec, fv=support1(p);
     3510  std::vector<std::vector<int> > fvs0, fvs;
     3511  vec.push_back(vert);
     3512  fvs.push_back(vec);
     3513  fvs0=b_subsets(fv);
     3514  fvs0=vsMinusv(fvs0,fv);
     3515  for(i=0;i<fvs0.size();i++)
     3516  {
     3517    vec=fvs0[i];
     3518    vec.push_back(vert);
     3519    fvs.push_back(vec);
     3520  }
     3521  return (fvs);
     3522}
     3523
     3524
     3525
     3526
     3527
     3528
     3529
     3530// the size of p's support must be 3
     3531//returns the new complex which is a triangulation based on the face p
     3532ideal triangulations1(ideal h, poly p, int vert)
     3533{
     3534  std::vector<int> vec, pv=support1(p);
     3535  std::vector<std::vector<int> > vecs=supports(h),vs,vs0;
     3536  vs0=triface(p,vert);
     3537  vecs=vsMinusv(vecs, pv);
     3538  vecs=vsUnion(vecs,vs0);
     3539  //PrintS("This is the new simplicial complex according to the face \n"); pWrite(p);
     3540  //PrintS("is:\n");
     3541  //listsprint(vecs);
     3542
     3543  ideal re=idMaken(vecs);
     3544
     3545  return re;
     3546}
     3547
     3548
     3549
     3550
     3551/*
     3552ideal triangulations1(ideal h)
     3553{
     3554  int i,vert=currRing->N+1;
     3555  std::vector<int> vec;
     3556  std::vector<std::vector<int> > vecs=supports(h),vs,vs0;
     3557  for (i=0;i<vecs.size();i++)
     3558  {
     3559    if((vecs[i]).size()==3)
     3560    {
     3561      vs0=triface(vecs[i],vert);
     3562      vs=vsMinusv(vecs,vecs[i]);
     3563      vs=vsUnion(vs,vs0);
     3564      PrintS("This is the new simplicial complex according to the face \n");listprint(vecs[i]);
     3565      PrintS("is:\n");
     3566      listsprint(vs);
     3567    }
     3568    //else if((vecs[i]).size()==4)
     3569      //tetraface(vecs[i]);
     3570  }
     3571  //ideal hh=idMaken(vs);
     3572  return h;
     3573}*/
     3574
     3575
     3576std::vector<int> commonedge(poly p, poly q)
     3577{
     3578  int i,j;
     3579  std::vector<int> ev, fv1= support1(p), fv2= support2(q);
     3580  for(i=0;i<fv1.size();i++)
     3581  {
     3582    if(IsinL(fv1[i], fv2))
     3583      ev.push_back(fv1[i]);
     3584  }
     3585  return ev;
     3586}
     3587
     3588
     3589intvec *edgemat(poly p, poly q)
     3590{
     3591  intvec *m;
     3592  int i,j;
     3593  std::vector<int> dg=commonedge(p, q);
     3594  int lg=dg.size();
     3595  m=new intvec(lg);
     3596  if(lg!=0)
     3597  {
     3598    m=new intvec(lg);
     3599    for(i=0;i<lg;i++)
     3600    {
     3601        (*m)[i]=dg[i];
     3602    }
     3603  }
     3604  return (m);
     3605}
     3606
     3607// case 2 the new face
     3608std::vector<std::vector<int> > tetraface(poly p, poly q, int vert)
     3609{
     3610  int i;
     3611  std::vector<int> ev=commonedge(p, q), vec, fv1=support1(p), fv2=support1(q);
     3612  std::vector<std::vector<int> > fvs1, fvs2, fvs;
     3613  vec.push_back(vert);
     3614  fvs.push_back(vec);
     3615  fvs1=b_subsets(fv1);
     3616  fvs2=b_subsets(fv2);
     3617  fvs1=vsMinusv(fvs1, fv1);
     3618  fvs2=vsMinusv(fvs2, fv2);
     3619  fvs2=vsUnion(fvs1, fvs2);
     3620  fvs2=vsMinusv(fvs2, ev);
     3621  for(i=0;i<fvs2.size();i++)
     3622  {
     3623    vec=fvs2[i];
     3624    vec.push_back(vert);
     3625    fvs.push_back(vec);
     3626  }
     3627  return (fvs);
     3628}
     3629
     3630
     3631//if p and q have  a common edge
     3632ideal triangulations2(ideal h, poly p, poly q, int vert)
     3633{
     3634  int i,j;
     3635  std::vector<int> ev, fv1=support1(p), fv2=support1(q);
     3636  std::vector<std::vector<int> > vecs=supports(h), vs1;
     3637  ev=commonedge(p, q);
     3638  vecs=vsMinusv(vecs, ev);
     3639  vecs=vsMinusv(vecs,fv1);
     3640  vecs=vsMinusv(vecs,fv2);
     3641  vs1=tetraface(p, q, vert);
     3642  vecs=vsUnion(vecs,vs1);
     3643  ideal hh=idMaken(vecs);
     3644  return hh;
     3645}
     3646
     3647
     3648
     3649
     3650// case 2 the new face
     3651std::vector<std::vector<int> > penface(poly p, poly q, poly g, int vert)
     3652{
     3653  int i, en=0;
     3654  std::vector<int> ev1=commonedge(p, q), ev2=commonedge(p, g), ev3=commonedge(q, g), ind, vec, fv1=support1(p), fv2=support1(q), fv3=support1(g);
     3655  std::vector<std::vector<int> > fvs1, fvs2, fvs3, fvs, evec;
     3656  evec.push_back(ev1);
     3657  evec.push_back(ev2);
     3658  evec.push_back(ev3);
     3659  for(i=0;i<evec.size();i++)
     3660  {
     3661    if(evec[i].size()==2)
     3662    {
     3663      en++;
     3664    }
     3665  }
     3666  if(en==2)
     3667  {
     3668    vec.push_back(vert);
     3669    fvs.push_back(vec);
     3670    fvs1=b_subsets(fv1);
     3671    fvs2=b_subsets(fv2);
     3672    fvs3=b_subsets(fv3);
     3673    fvs1=vsMinusv(fvs1, fv1);
     3674    fvs2=vsMinusv(fvs2, fv2);
     3675    fvs3=vsMinusv(fvs3, fv3);
     3676    fvs3=vsUnion(fvs3, fvs2);
     3677    fvs3=vsUnion(fvs3, fvs1);
     3678    for(i=0;i<evec.size();i++)
     3679    {
     3680      if(evec[i].size()==2)
     3681      {
     3682        fvs3=vsMinusv(fvs3, evec[i]);
     3683      }
     3684    }
     3685    for(i=0;i<fvs3.size();i++)
     3686    {
     3687      vec=fvs3[i];
     3688      vec.push_back(vert);
     3689      fvs.push_back(vec);
     3690    }
     3691  }
     3692  return (fvs);
     3693}
     3694
     3695
     3696
     3697ideal triangulations3(ideal h, poly p, poly q, poly g, int vert)
     3698{
     3699  int i,j;
     3700  std::vector<int> ev1=commonedge(p, q), ev2=commonedge(p, g), ev3=commonedge(q, g), fv1=support1(p), fv2=support1(q), fv3=support1(g);
     3701  std::vector<std::vector<int> > vecs=supports(h), vs1, evec;
     3702  evec.push_back(ev1);
     3703  evec.push_back(ev2);
     3704  evec.push_back(ev3);
     3705  for(i=0;i<evec.size();i++)
     3706  {
     3707    if(evec[i].size()==2)
     3708    {
     3709      vecs=vsMinusv(vecs, evec[i]);
     3710    }
     3711  }
     3712  vecs=vsMinusv(vecs,fv1);
     3713  vecs=vsMinusv(vecs,fv2);
     3714  vecs=vsMinusv(vecs,fv3);
     3715  vs1=penface(p, q, g, vert);
     3716  vecs=vsUnion(vecs,vs1);
     3717  ideal hh=idMaken(vecs);
     3718  return hh;
     3719}
     3720
     3721
     3722//returns p's valency in h
     3723//p must be a vertex
     3724int valency(ideal h, poly p)
     3725{
     3726  int i, val=0;
     3727  std::vector<int> ev=support1(pCopy(p));
     3728  int ver=ev[0];
     3729//PrintS("the vertex is :\n"); listprint(p);
     3730  std::vector<std::vector<int> > vecs=supports(idCopy(h));
     3731  for(i=0;i<vecs.size();i++)
     3732  {
     3733    if(vecs[i].size()==2 && IsinL(ver, vecs[i]))
     3734      val++;
     3735  }
     3736  return (val);
     3737}
     3738
     3739/*ideal triangulations2(ideal h)
     3740{
     3741  int i,j,vert=currRing->N+1;
     3742  std::vector<int> ev;
     3743  std::vector<std::vector<int> > vecs=supports(h),vs,vs0,vs1;
     3744  vs0=tetrasets(h);
     3745  for (i=0;i<vs0.size();i++)
     3746  {
     3747    for(j=i;j<vs0.size();j++)
     3748    {
     3749      ev=commonedge(vs0[i],vs0[j]);
     3750      if(ev.size()==2)
     3751      {
     3752        vecs=vsMinusv(vecs, ev);
     3753        vs=vsMinusv(vecs,vs0[i]);
     3754        vs=vsMinusv(vecs,vs0[j]);
     3755        vs1=tetraface(vs0[i],vs0[j],vert);
     3756        vs=vsUnion(vs,vs1);
     3757        PrintS("This is the new simplicial complex according to the face 1 \n");listprint(vecs[i]);
     3758PrintS("face 2: \n");
     3759        PrintS("is:\n");
     3760        listsprint(vs);
     3761      }
     3762    }
     3763
     3764    //else if((vecs[i]).size()==4)
     3765      //tetraface(vecs[i]);
     3766  }
     3767  //ideal hh=idMaken(vs);
     3768  return h;
     3769}*/
     3770
     3771
     3772
     3773/*********************************For computation of X_n***********************************/
     3774std::vector<std::vector<int> > vsMinusvs(std::vector<std::vector<int> > vs1, std::vector<std::vector<int> > vs2)
     3775{
     3776  int i;
     3777  std::vector<std::vector<int> > vs=vs1;
     3778  for(i=0;i<vs2.size();i++)
     3779  {
     3780    vs=vsMinusv(vs, vs2[i]);
     3781  }
     3782  return vs;
     3783}
     3784
     3785
     3786std::vector<std::vector<int> > vs_subsets(std::vector<std::vector<int> > vs)
     3787{
     3788  std::vector<std::vector<int> >  sset, bv;
     3789  for(int i=0;i<vs.size();i++)
     3790  {
     3791    bv=b_subsets(vs[i]);
     3792    sset=vsUnion(sset, bv);
     3793  }
     3794  return sset;
     3795}
     3796
     3797
     3798
     3799std::vector<std::vector<int> > p_constant(ideal Xo,  ideal Sigma)
     3800{
     3801  std::vector<std::vector<int> > xs=supports(idCopy(Xo)), ss=supports(idCopy(Sigma)), fvs1;
     3802  fvs1=vs_subsets(ss);
     3803  fvs1=vsMinusvs(xs, fvs1);
     3804  return fvs1;
     3805}
     3806
     3807
     3808std::vector<std::vector<int> > p_change(ideal Sigma)
     3809{
     3810  std::vector<std::vector<int> > ss=supports(idCopy(Sigma)), fvs;
     3811  fvs=vs_subsets(ss);
     3812  return (fvs);
     3813}
     3814
     3815
     3816
     3817std::vector<std::vector<int> > p_new(ideal Xo, ideal Sigma)
     3818{
     3819  int vert=0;
     3820  std::vector<std::vector<int> > ss=supports(idCopy(Sigma)), fvs;
     3821  for(int i=1;i<=currRing->N;i++)
     3822  {
     3823    for(int j=0;j<IDELEMS(Xo);j++)
     3824    {
     3825      if(pGetExp(Xo->m[j],i)>0)
     3826      {
     3827        vert=i+1;
     3828        break;
     3829      }
     3830    }
     3831  }
     3832  int typ=ss.size();
     3833  if(typ==1)
     3834  {
     3835    fvs=triface(Sigma->m[0], vert);
     3836  }
     3837  else if(typ==2)
     3838  {
     3839     fvs=tetraface(Sigma->m[0], Sigma->m[1], vert);
     3840  }
     3841  else
     3842  {
     3843     fvs=penface(Sigma->m[0], Sigma->m[1], Sigma->m[2], vert);
     3844  }
     3845  return (fvs);
     3846}
     3847
     3848
     3849
     3850
     3851ideal c_New(ideal Io, ideal sig)
     3852{
     3853  poly p, q, g;
     3854  std::vector<std::vector<int> > vs1=p_constant(Io, sig), vs2=p_change(sig), vs3=p_new(Io, sig), vsig=supports(sig), vs;
     3855  std::vector<int> ev;
     3856  int ednum=vsig.size();
     3857  if(ednum==2)
     3858  {
     3859    vsig.push_back(commonedge(sig->m[0], sig->m[1]));
     3860  }
     3861  else if(ednum==3)
     3862  {
     3863    for(int i=0;i<IDELEMS(sig);i++)
     3864    {
     3865      for(int j=i+1;j<IDELEMS(sig);j++)
     3866      {
     3867        ev=commonedge(sig->m[i], sig->m[j]);
     3868        if(ev.size()==2)
     3869        {
     3870          vsig.push_back(ev);
     3871        }
     3872      }
     3873    }
     3874  }
     3875//PrintS("the first part is:\n");id_print(idMaken(vs1));
     3876//PrintS("the second part is:\n");id_print(idMaken(vsig));
     3877//PrintS("the third part is:\n");id_print(idMaken(vs3));
     3878  vs2=vsMinusvs(vs2, vsig);
     3879//PrintS("the constant part2 is:\n");id_print(idMaken(vs2));
     3880  vs=vsUnion(vs2, vs1);
     3881//PrintS("the constant part is:\n");id_print(idMaken(vs));
     3882  vs=vsUnion(vs, vs3);
     3883//PrintS("the whole part is:\n");id_print(idMaken(vs));
     3884  return(idMaken(vs));
     3885}
     3886
     3887
     3888
     3889
     3890std::vector<std::vector<int> > phi1(poly a,  ideal Sigma)
     3891{
     3892  std::vector<std::vector<int> > ss=supports(idCopy(Sigma)), fvs;
     3893  std::vector<int> av=support1(a), intvec, vv;
     3894  for(int i=0;i<ss.size();i++)
     3895  {
     3896    intvec=vecIntersection(ss[i], av);
     3897    if(intvec.size()==av.size())
     3898    {
     3899      vv=vecMinus(ss[i], av);
     3900      fvs.push_back(vv);
     3901    }
     3902  }
     3903  return fvs;
     3904}
     3905
     3906
     3907
     3908std::vector<std::vector<int> > phi2(poly a, ideal Xo, ideal Sigma, int vert)
     3909{
     3910
     3911  std::vector<std::vector<int> > ss=p_new(Sigma, Xo), fvs;
     3912  std::vector<int> av=support1(a), intvec, vv;
     3913  for(int i=0;i<ss.size();i++)
     3914  {
     3915    intvec=vecIntersection(ss[i], av);
     3916    if(intvec.size()==av.size())
     3917    {
     3918      vv=vecMinus(ss[i], av);
     3919      fvs.push_back(vv);
     3920    }
     3921  }
     3922  return fvs;
     3923}
     3924
     3925
     3926std::vector<std::vector<int> > links_new(poly a, ideal Xo, ideal Sigma, int vert, int ord)
     3927{
     3928  std::vector<int> av=support1(a);
     3929  std::vector<std::vector<int> > lko, lkn, lk1, lk2;
     3930  lko=links(a, Xo);
     3931  if(ord==1)
     3932    return lko;
     3933  if(ord==2)
     3934  {
     3935
     3936    lk1=phi1(a, Sigma);
     3937    lk2=phi2(a, Xo, Sigma, vert);
     3938    lkn=vsMinusvs(lko, lk1);
     3939    lkn=vsUnion(lkn, lk2);
     3940    return lkn;
     3941  }
     3942  if(ord==3)
     3943  {
     3944    lkn=phi2(a, Xo, Sigma, vert);
     3945    return lkn;
     3946  }
     3947  WerrorS("Cannot find the links smartly!");
     3948}
     3949
     3950
     3951
     3952
     3953//returns 1 if there is a real divisor of b not in Xs
     3954int existIn(poly b, ideal Xs)
     3955{
     3956  std::vector<int> bv=support1(pCopy(b));
     3957  std::vector<std::vector<int> > xvs=supports(idCopy(Xs)), bs=b_subsets(bv);
     3958  bs=vsMinusv(bs, bv);
     3959  for(int i=0;i<bs.size();i++)
     3960  {
     3961    if(!vInvsl(bs[i], xvs))
     3962    {
     3963      return 1;
     3964    }
     3965  }
     3966  return 0;
     3967}
     3968
     3969
     3970int isoNum(poly p, ideal I, poly a, poly b)
     3971{
     3972  int i;
     3973  std::vector<std::vector<int> > vs=supports(idCopy(I));
     3974  std::vector<int> v1=support1(a), v2=support1(b), v=support1(p);
     3975  std::vector<int>  vp, iv=phimagel(v, v1, v2);
     3976  for(i=0;i<IDELEMS(I);i++)
     3977  {
     3978    vp=support1(pCopy(I->m[i]));
     3979    if(vEvl(iv, phimagel(vp, v1, v2)))
     3980    {
     3981      return (i+1);
     3982    }
     3983  }
     3984  return (0);
     3985}
     3986
     3987
     3988
     3989
     3990int ifIso(poly p, poly q, poly f, poly g, poly a, poly b)
     3991{
     3992  int i;
     3993  std::vector<int> va=support1(a), vb=support1(b), vp=support1(p),  vq=support1(q), vf=support1(f), vg=support1(g);
     3994  std::vector<int>   v1=phimagel(vp, va, vb), v2=phimagel(vq, va, vb), v3=phimagel(vf, va, vb), v4=phimagel(vg, va, vb);
     3995  if((vEvl(v1, v3)&& vEvl(v2,v4))||(vEvl(v1, v4)&& vEvl(v2,v3)) )
     3996  {
     3997    return (1);
     3998  }
     3999  return (0);
     4000}
     4001
     4002
     4003
     4004
     4005ideal idMinusp(ideal I, poly p)
     4006{
     4007  ideal h=idInit(1,1);
     4008  int i,j,eq=0;
     4009  for(i=0;i<IDELEMS(I);i++)
     4010  {
     4011    if(!p_EqualPolys(I->m[i], p, currRing))
     4012    {
     4013      idInsertPoly(h, pCopy(I->m[i]));
     4014    }
     4015  }
     4016  idSkipZeroes(h);
     4017  return h;
     4018}
    28364019
    28374020
    … …  
    28674050  }
    28684051  std::vector<int> vec=v_minus(av,bv);
    2869 //PrintS("the degree is:\n");
    2870 //listprint(vec);
     4052  //PrintS("The degree is:\n");
     4053  //listprint(vec);
    28714054  return vec;
    28724055}
    … …  
    28754058
    28764059
     4060
     4061
     4062/********************************for stellar subdivision******************************/
     4063
     4064
     4065std::vector<std::vector<int> > star(poly a, ideal h)
     4066{
     4067  int i;
     4068  std::vector<std::vector<int> > st,X=supports(h);
     4069  std::vector<int> U,av=support1(a);
     4070  for(i=0;i<X.size();i++)
     4071  {
     4072    U=vecUnion(av,X[i]);
     4073    if(vInvsl(U,X))
     4074    {
     4075      st.push_back(X[i]);
     4076    }
     4077  }
     4078  return st;
     4079}
     4080
     4081
     4082std::vector<std::vector<int> > boundary(poly a)
     4083{
     4084  std::vector<int> av=support1(a), vec;
     4085  std::vector<std::vector<int> > vecs;
     4086  vecs=b_subsets(av);
     4087  vecs.push_back(vec);
     4088  vecs=vsMinusv(vecs, av);
     4089  return vecs;
     4090}
     4091
     4092
     4093
     4094
     4095
     4096
     4097std::vector<std::vector<int> > stellarsub(poly a, ideal h)
     4098{
     4099  std::vector<std::vector<int> > vecs_minus, vecs_plus, lk=links(a,h), hvs=supports(h), sub, bys=boundary(a);
     4100  std::vector<int> av=support1(a), vec, vec_n;
     4101  int i,j,vert=0;
     4102  for(i=1;i<=currRing->N;i++)
     4103  {
     4104    for(j=0;j<IDELEMS(h);j++)
     4105    {
     4106      if(pGetExp(h->m[j],i)>0)
     4107      {
     4108        vert=i+1;
     4109        break;
     4110      }
     4111    }
     4112  }
     4113  vec_n.push_back(vert);
     4114  for(i=0;i<lk.size();i++)
     4115  {
     4116    vec=vecUnion(av, lk[i]);
     4117    vecs_minus.push_back(vec);
     4118    for(j=0;j<bys.size();j++)
     4119    {
     4120      vec=vecUnion(lk[i], vec_n);
     4121      vec=vecUnion(vec, bys[j]);
     4122      vecs_plus.push_back(vec);
     4123    }
     4124  }
     4125  sub=vsMinusvs(hvs, vecs_minus);
     4126  sub=vsUnion(sub, vecs_plus);
     4127  return(sub);
     4128}
     4129
     4130
     4131std::vector<std::vector<int> > bsubsets_1(poly b)
     4132{
     4133  std::vector<int>  bvs=support1(b), vs;
     4134  std::vector<std::vector<int> > bset;
     4135  for(int i=0;i<bvs.size();i++)
     4136  {
     4137    for(int j=0;j<bvs.size(), j!=i; j++)
     4138    {
     4139      vs.push_back(bvs[j]);
     4140    }
     4141    bset.push_back(vs);
     4142    vs.resize(0);
     4143  }
     4144  return bset;
     4145}
     4146
     4147
     4148
     4149/***************************for time testing******************************/
     4150ideal T_1h(ideal h)
     4151{
     4152  int i, j;
     4153  //std::vector < intvec > T1;
     4154  ideal ai=p_a(h), bi;
     4155  //intvec *L;
     4156  for(i=0;i<IDELEMS(ai);i++)
     4157  {
     4158    bi=p_b(h,ai->m[i]);
     4159    if(!idIs0(bi))
     4160    {
     4161      for(j=0;j<IDELEMS(bi);j++)
     4162      {
     4163        //PrintS("This is for:\n");pWrite(ai->m[i]); pWrite(bi->m[j]);
     4164        gradedpiece1nl(h,ai->m[i],bi->m[j], 0);
     4165        //PrintS("Succeed!\n");
     4166        //T1.push_back(L);
     4167      }
     4168    }
     4169  }
     4170  TimeShow(t_construct, t_solve, t_value, t_total);
     4171  return h;
     4172
     4173}
    28774174/**************************************interface T1****************************************/
    28784175/*
    … …  
    29054202
    29064203
     4204
     4205BOOLEAN SRideal(leftv res, leftv args)
     4206{
     4207  leftv h=args;
     4208   if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4209   {
     4210     ideal hh=(ideal)h->Data();
     4211     res->rtyp =IDEAL_CMD;
     4212     res->data =idsrRing(hh);
     4213   }
     4214  return false;
     4215}
     4216
     4217
     4218
     4219
     4220
     4221
    29074222BOOLEAN idcomplement(leftv res, leftv args)
    29084223{
    … …  
    29204235
    29214236
     4237
     4238
     4239BOOLEAN t1h(leftv res, leftv args)
     4240{
     4241  leftv h=args;
     4242   if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4243   {
     4244     ideal hh=(ideal)h->Data();
     4245     res->rtyp =IDEAL_CMD;
     4246     res->data =T_1h(hh);
     4247   }
     4248  return false;
     4249}
     4250
     4251
    29224252BOOLEAN idsr(leftv res, leftv args)
    29234253{
    … …  
    29264256  {
    29274257     ideal h1= (ideal)h->Data();
    2928 //T1(h1);
    29294258     h   = h->next;
    29304259     if((h != NULL)&&(h->Typ() == POLY_CMD))
    … …  
    29484277  int i,j;
    29494278  std::vector<int> dg=gdegree(a,b);
    2950 //PrintS("This is the degree of a and b\n");
    2951 //listprint(dg);
    29524279  int lg=dg.size();
    29534280  m=new intvec(lg);
    … …  
    29584285    {
    29594286        (*m)[i]=dg[i];
    2960         //Print("This is the %dth degree of a and b: %d, and %d is copied\n",i,dg[i],(*m)[i]);
    2961     }
    2962   }
    2963   /*for(j=0;j<lg;j++)
    2964   {
    2965     Print("[%d]: %d\n",j+1,(*m)[j]);
    2966   }*/
    2967   //(m)->show(1,1);
    2968 return (m);
     4287    }
     4288  }
     4289  return (m);
    29694290}
    29704291
    … …  
    29904311
    29914312
     4313BOOLEAN comedg(leftv res, leftv args)
     4314{
     4315  leftv h=args;
     4316  if((h != NULL)&&(h->Typ() == POLY_CMD))
     4317  {
     4318     poly p= (poly)h->Data();
     4319     h   = h->next;
     4320     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4321     {
     4322       poly q= (poly)h->Data();
     4323       res->rtyp =INTVEC_CMD;
     4324       res->data =edgemat(p,q);
     4325     }
     4326  }
     4327  return false;
     4328}
     4329
     4330
     4331
     4332
    29924333BOOLEAN fb(leftv res, leftv args)
    29934334{
    … …  
    29984339     res->rtyp =IDEAL_CMD;
    29994340     res->data =findb(h1);
     4341  }
     4342  return false;
     4343}
     4344
     4345
     4346BOOLEAN pa(leftv res, leftv args)
     4347{
     4348  leftv h=args;
     4349  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4350  {
     4351     ideal h1= (ideal)h->Data();
     4352     res->rtyp =IDEAL_CMD;
     4353     res->data =p_a(h1);
     4354  }
     4355  return false;
     4356}
     4357
     4358
     4359
     4360BOOLEAN makeSimplex(leftv res, leftv args)
     4361{
     4362  leftv h=args;
     4363  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4364  {
     4365     ideal h1= (ideal)h->Data();
     4366     res->rtyp =IDEAL_CMD;
     4367     res->data =complementsimplex(h1);
     4368  }
     4369  return false;
     4370}
     4371
     4372
     4373BOOLEAN pb(leftv res, leftv args)
     4374{
     4375  leftv h=args;
     4376  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4377  {
     4378     ideal h1= (ideal)h->Data();
     4379     h   = h->next;
     4380     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4381     {
     4382       poly p= (poly)h->Data();
     4383       res->rtyp =IDEAL_CMD;
     4384       res->data =p_b(h1,p);
     4385     }
    30004386  }
    30014387  return false;
    … …  
    30504436
    30514437
     4438BOOLEAN fgpl(leftv res, leftv args)
     4439{
     4440  leftv h=args;
     4441  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4442  {
     4443     ideal h1= (ideal)h->Data();
     4444     h   = h->next;
     4445     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4446     {
     4447       poly p= (poly)h->Data();
     4448       h   = h->next;
     4449       if((h != NULL)&&(h->Typ() == POLY_CMD))
     4450       {
     4451         poly q= (poly)h->Data();
     4452         h   = h->next;
     4453         if((h != NULL)&&(h->Typ() == INT_CMD))
     4454         {
     4455           int d= (int)(long)h->Data();
     4456           res->rtyp =INTVEC_CMD;
     4457           res->data =gradedpiece1nl(h1,p,q,d);
     4458         }
     4459       }
     4460     }
     4461  }
     4462  return false;
     4463}
     4464
     4465
    30524466
    30534467BOOLEAN genstt(leftv res, leftv args)
    … …  
    30974511
    30984512
     4513BOOLEAN sgpl(leftv res, leftv args)
     4514{
     4515  leftv h=args;
     4516  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4517  {
     4518     ideal h1= (ideal)h->Data();
     4519     h   = h->next;
     4520     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4521     {
     4522       poly p= (poly)h->Data();
     4523       h   = h->next;
     4524       if((h != NULL)&&(h->Typ() == POLY_CMD))
     4525       {
     4526         poly q= (poly)h->Data();
     4527         res->rtyp =INTVEC_CMD;
     4528         res->data =gradedpiece2nl(h1,p,q);
     4529       }
     4530     }
     4531  }
     4532  return false;
     4533}
     4534
     4535
    30994536BOOLEAN Links(leftv res, leftv args)
    31004537{
    … …  
    31154552}
    31164553
     4554BOOLEAN isSim(leftv res, leftv args)
     4555{
     4556  leftv h=args;
     4557  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4558  {
     4559     ideal h1= (ideal)h->Data();
     4560     res->rtyp =IDEAL_CMD;
     4561     res->data =IsSimplex(h1);
     4562  }
     4563  return false;
     4564}
     4565
     4566
     4567BOOLEAN nfaces1(leftv res, leftv args)
     4568{
     4569  leftv h=args;
     4570  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4571  {
     4572     ideal h1= (ideal)h->Data();
     4573     h   = h->next;
     4574     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4575     {
     4576       poly p= (poly)h->Data();
     4577       h   = h->next;
     4578       if((h != NULL)&&(h->Typ() == INT_CMD))
     4579       {
     4580         int d= (int)(long)h->Data();
     4581         res->rtyp =IDEAL_CMD;
     4582         res->data =triangulations1(h1, p, d);
     4583       }
     4584     }
     4585  }
     4586  return false;
     4587}
     4588
     4589
     4590BOOLEAN nfaces2(leftv res, leftv args)
     4591{
     4592  leftv h=args;
     4593  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4594  {
     4595     ideal h1= (ideal)h->Data();
     4596     h   = h->next;
     4597     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4598     {
     4599       poly p= (poly)h->Data();
     4600       h   = h->next;
     4601       if((h != NULL)&&(h->Typ() == POLY_CMD))
     4602       {
     4603         poly q= (poly)h->Data();
     4604         h   = h->next;
     4605         if((h != NULL)&&(h->Typ() == INT_CMD))
     4606         {
     4607           int d= (int)(long)h->Data();
     4608           res->rtyp =IDEAL_CMD;
     4609           res->data =triangulations2(h1,p,q,d);
     4610         }
     4611       }
     4612     }
     4613  }
     4614  return false;
     4615}
     4616
     4617
     4618BOOLEAN nfaces3(leftv res, leftv args)
     4619{
     4620  leftv h=args;
     4621  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4622  {
     4623     ideal h1= (ideal)h->Data();
     4624     h   = h->next;
     4625     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4626     {
     4627       poly p= (poly)h->Data();
     4628       h   = h->next;
     4629       if((h != NULL)&&(h->Typ() == POLY_CMD))
     4630       {
     4631         poly q= (poly)h->Data();
     4632         h   = h->next;
     4633         if((h != NULL)&&(h->Typ() == POLY_CMD))
     4634         {
     4635           poly g= (poly)h->Data();
     4636           h   = h->next;
     4637           if((h != NULL)&&(h->Typ() == INT_CMD))
     4638           {
     4639             int d= (int)(long)h->Data();
     4640             res->rtyp =IDEAL_CMD;
     4641             res->data =triangulations3(h1,p,q,g,d);
     4642           }
     4643         }
     4644       }
     4645     }
     4646  }
     4647  return false;
     4648}
     4649
     4650
     4651
     4652
     4653
     4654BOOLEAN eqsolve1(leftv res, leftv args)
     4655{
     4656  leftv h=args;int i;
     4657  std::vector<int> bset,bs;
     4658  std::vector<std::vector<int> > gset;
     4659  if((h != NULL)&&(h->Typ() == INT_CMD))
     4660  {
     4661     int n= (int)(long)h->Data();
     4662     h   = h->next;
     4663     if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4664     {
     4665       ideal bi= (ideal)h->Data();
     4666       h   = h->next;
     4667       if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4668       {
     4669         ideal gi= (ideal)h->Data();
     4670         for(i=0;i<IDELEMS(bi);i++)
     4671         {
     4672           bs=support1(bi->m[i]);
     4673           if(bs.size()==1)
     4674             bset.push_back(bs[0]);
     4675           else if(bs.size()==0)
     4676             ;
     4677           else
     4678           {
     4679             WerrorS("Errors in T^1 Equations Solving!");
     4680             usleep(1000000);
     4681             assert(false);
     4682           }
     4683
     4684         }
     4685         gset=supports2(gi);
     4686         res->rtyp =INTVEC_CMD;
     4687         std::vector<std::vector<int> > vecs=eli2(n,bset,gset);
     4688         res->data =Tmat(vecs);
     4689       }
     4690     }
     4691  }
     4692  return false;
     4693}
     4694
     4695
     4696BOOLEAN tsets(leftv res, leftv args)
     4697{
     4698  leftv h=args;
     4699  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4700  {
     4701     ideal h1= (ideal)h->Data();
     4702     res->rtyp =IDEAL_CMD;
     4703     res->data =trisets(h1);
     4704  }
     4705  return false;
     4706}
     4707
     4708
     4709
     4710
     4711
     4712BOOLEAN Valency(leftv res, leftv args)
     4713{
     4714  leftv h=args;
     4715  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4716  {
     4717     ideal h1= (ideal)h->Data();
     4718     h   = h->next;
     4719     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4720     {
     4721       poly p= (poly)h->Data();
     4722       res->rtyp =INT_CMD;
     4723       res->data =(void *)(long)valency(h1,p);
     4724     }
     4725  }
     4726  return false;
     4727}
     4728
     4729
     4730
     4731
     4732BOOLEAN nabvl(leftv res, leftv args)
     4733{
     4734  leftv h=args;
     4735  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4736  {
     4737     ideal h1= (ideal)h->Data();
     4738     h   = h->next;
     4739     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4740     {
     4741       poly p= (poly)h->Data();
     4742       h   = h->next;
     4743       if((h != NULL)&&(h->Typ() == POLY_CMD))
     4744       {
     4745         poly q= (poly)h->Data();
     4746         res->rtyp =IDEAL_CMD;
     4747         std::vector<std::vector<int> > vecs=supports(h1);
     4748         std::vector<int> pv=support1(p), qv=support1(q);
     4749         res->data =idMaken(Nabv(vecs,pv,qv));
     4750       }
     4751     }
     4752  }
     4753  return false;
     4754}
     4755
     4756
     4757
     4758BOOLEAN tnabvl(leftv res, leftv args)
     4759{
     4760  leftv h=args;
     4761  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4762  {
     4763     ideal h1= (ideal)h->Data();
     4764     h   = h->next;
     4765     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4766     {
     4767       poly p= (poly)h->Data();
     4768       h   = h->next;
     4769       if((h != NULL)&&(h->Typ() == POLY_CMD))
     4770       {
     4771         poly q= (poly)h->Data();
     4772         res->rtyp =IDEAL_CMD;
     4773         std::vector<std::vector<int> > vecs=supports(h1), sbv,tnbr;
     4774         std::vector<int> pv=support1(p), qv=support1(q);
     4775         std::vector<std::vector<int> > nvs=Nabv(vecs, pv, qv);
     4776         ideal sub=psubset(q);
     4777         sbv=supports(sub);
     4778         std::vector<int> tnv =tnab(vecs,nvs,sbv);
     4779         for(int i=0;i<tnv.size();i++)
     4780         {
     4781           tnbr.push_back(nvs[tnv[i]]);
     4782         }
     4783         res->data =idMaken(tnbr);
     4784       }
     4785     }
     4786  }
     4787  return false;
     4788}
     4789
     4790
     4791BOOLEAN vsIntersec(leftv res, leftv args)
     4792{
     4793  leftv h=args;
     4794  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4795  {
     4796     ideal h1= (ideal)h->Data();
     4797     h   = h->next;
     4798     if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4799     {
     4800       ideal h2= (ideal)h->Data();
     4801       res->rtyp =INT_CMD;
     4802       std::vector<std::vector<int> > vs1=supports(h1), vs2=supports(h2);
     4803       res->data =(void *)(long)(vsIntersection(vs1, vs2).size());
     4804     }
     4805  }
     4806  return false;
     4807}
     4808
     4809
     4810BOOLEAN mabvl(leftv res, leftv args)
     4811{
     4812  leftv h=args;
     4813  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4814  {
     4815     ideal h1= (ideal)h->Data();
     4816     h   = h->next;
     4817     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4818     {
     4819       poly p= (poly)h->Data();
     4820       h   = h->next;
     4821       if((h != NULL)&&(h->Typ() == POLY_CMD))
     4822       {
     4823         poly q= (poly)h->Data();
     4824         res->rtyp =IDEAL_CMD;
     4825         res->data =idMaken(Mabv(h1,p,q));
     4826       }
     4827     }
     4828  }
     4829  return false;
     4830}
     4831
     4832
     4833
     4834BOOLEAN nabtvl(leftv res, leftv args)
     4835{
     4836  leftv h=args;
     4837  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4838  {
     4839     ideal h1= (ideal)h->Data();
     4840     h   = h->next;
     4841     if((h != NULL)&&(h->Typ() == POLY_CMD))
     4842     {
     4843       poly p= (poly)h->Data();
     4844       h   = h->next;
     4845       if((h != NULL)&&(h->Typ() == POLY_CMD))
     4846       {
     4847         poly q= (poly)h->Data();
     4848         std::vector<std::vector<int> > hvs=supports(h1), nv, ntvs;
     4849         std::vector<int> av=support1(p), bv=support1(q);
     4850         nv=Nabv(hvs,av,bv);
     4851         ntvs=nabtv( hvs, nv, av, bv);
     4852         std::vector<std::vector<poly> > pvs=idMakei(nv,ntvs);
     4853         ideal gens=idInit(1,1);
     4854         for(int i=0;i<pvs.size();i++)
     4855         {
     4856           idInsertPoly(gens,pvs[i][0]);
     4857           idInsertPoly(gens,pvs[i][1]);
     4858         }
     4859         idSkipZeroes(gens);
     4860         res->rtyp =IDEAL_CMD;
     4861         res->data =gens;
     4862       }
     4863     }
     4864  }
     4865  return false;
     4866}
     4867
     4868
     4869BOOLEAN linkn(leftv res, leftv args)
     4870{
     4871  leftv h=args;
     4872  if((h != NULL)&&(h->Typ() == POLY_CMD))
     4873  {
     4874    poly a= (poly)h->Data();
     4875    h   = h->next;
     4876    if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4877    {
     4878      ideal Xo= (ideal)h->Data();
     4879      h   = h->next;
     4880      if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4881      {
     4882        ideal Sigma= (ideal)h->Data();
     4883        h   = h->next;
     4884        if((h != NULL)&&(h->Typ() == INT_CMD))
     4885        {
     4886          int vert= (int)(long)h->Data();
     4887          h   = h->next;
     4888          if((h != NULL)&&(h->Typ() == INT_CMD))
     4889          {
     4890            int ord= (int)(long)h->Data();
     4891            res->rtyp =IDEAL_CMD;
     4892            res->data =idMaken(links_new(a,  Xo,  Sigma,  vert,  ord));
     4893          }
     4894        }
     4895      }
     4896    }
     4897  }
     4898  return false;
     4899}
     4900
     4901
     4902
     4903BOOLEAN existsub(leftv res, leftv args)
     4904{
     4905  leftv h=args;
     4906  if((h != NULL)&&(h->Typ() == POLY_CMD))
     4907  {
     4908     poly p= (poly)h->Data();
     4909     h   = h->next;
     4910     if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4911     {
     4912       ideal h1= (ideal)h->Data();
     4913       res->rtyp =INT_CMD;
     4914       res->data =(void *)(long)existIn(p, h1);
     4915     }
     4916   }
     4917  return false;
     4918}
     4919
     4920
     4921BOOLEAN pConstant(leftv res, leftv args)
     4922{
     4923  leftv h=args;
     4924  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4925  {
     4926     ideal h1= (ideal)h->Data();
     4927     h   = h->next;
     4928     if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4929     {
     4930       ideal h2= (ideal)h->Data();
     4931       res->rtyp =IDEAL_CMD;
     4932       res->data =idMaken(p_constant(h1,h2));
     4933     }
     4934  }
     4935  return false;
     4936}
     4937
     4938BOOLEAN pChange(leftv res, leftv args)
     4939{
     4940  leftv h=args;
     4941  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4942  {
     4943     ideal h1= (ideal)h->Data();
     4944     res->rtyp =IDEAL_CMD;
     4945     res->data =idMaken(p_change(h1));
     4946  }
     4947  return false;
     4948}
     4949
     4950
     4951
     4952BOOLEAN p_New(leftv res, leftv args)
     4953{
     4954  leftv h=args;
     4955  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4956  {
     4957     ideal h1= (ideal)h->Data();
     4958     h   = h->next;
     4959     if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     4960     {
     4961       ideal h2= (ideal)h->Data();
     4962       res->rtyp =IDEAL_CMD;
     4963       res->data =idMaken(p_new(h1,h2));
     4964     }
     4965  }
     4966  return false;
     4967}
     4968
     4969
     4970
     4971
     4972BOOLEAN support(leftv res, leftv args)
     4973{
     4974  leftv h=args;
     4975  if((h != NULL)&&(h->Typ() == POLY_CMD))
     4976  {
     4977     poly p= (poly)h->Data();
     4978     res->rtyp =INT_CMD;
     4979     res->data =(void *)(long)(support1(p).size());
     4980  }
     4981  return false;
     4982}
     4983
     4984
     4985
     4986
     4987
     4988
     4989BOOLEAN bprime(leftv res, leftv args)
     4990{
     4991  leftv h=args;
     4992  if((h != NULL)&&(h->Typ() == POLY_CMD))
     4993  {
     4994     poly p= (poly)h->Data();
     4995     res->rtyp =IDEAL_CMD;
     4996     res->data =idMaken(bsubsets_1(p));
     4997  }
     4998  return false;
     4999}
     5000
     5001
     5002
     5003BOOLEAN psMinusp(leftv res, leftv args)
     5004{
     5005  leftv h=args;
     5006  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5007  {
     5008     ideal h1= (ideal)h->Data();
     5009     h   = h->next;
     5010     if((h != NULL)&&(h->Typ() == POLY_CMD))
     5011     {
     5012       poly p= (poly)h->Data();
     5013       res->rtyp =IDEAL_CMD;
     5014       res->data =idMinusp(h1, p);
     5015     }
     5016  }
     5017  return false;
     5018}
     5019
     5020
     5021
     5022BOOLEAN stellarremain(leftv res, leftv args)
     5023{
     5024  leftv h=args;
     5025  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5026  {
     5027     ideal h1= (ideal)h->Data();
     5028     h   = h->next;
     5029     if((h != NULL)&&(h->Typ() == POLY_CMD))
     5030     {
     5031       poly p= (poly)h->Data();
     5032       std::vector<std::vector<int> > st=star(p, h1);
     5033       std::vector<std::vector<int> > hvs=supports(h1);
     5034       std::vector<std::vector<int> > re= vsMinusvs(hvs, st);
     5035       res->rtyp =IDEAL_CMD;
     5036       res->data =idMaken(re);
     5037     }
     5038  }
     5039  return false;
     5040}
     5041
     5042
     5043BOOLEAN cNew(leftv res, leftv args)
     5044{
     5045  leftv h=args;
     5046  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5047  {
     5048     ideal h1= (ideal)h->Data();
     5049     h   = h->next;
     5050     if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5051     {
     5052       ideal h2= (ideal)h->Data();
     5053       res->rtyp =IDEAL_CMD;
     5054       res->data =c_New(h1, h2);
     5055     }
     5056  }
     5057  return false;
     5058}
     5059
     5060
     5061
     5062
     5063BOOLEAN stars(leftv res, leftv args)
     5064{
     5065  leftv h=args;
     5066  if((h != NULL)&&(h->Typ() == POLY_CMD))
     5067  {
     5068     poly p= (poly)h->Data();
     5069     h   = h->next;
     5070     if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5071     {
     5072       ideal h1= (ideal)h->Data();
     5073       res->rtyp =IDEAL_CMD;
     5074       res->data =idMaken(star(p, h1));
     5075     }
     5076   }
     5077  return false;
     5078}
     5079
     5080
     5081
     5082
     5083BOOLEAN stellarsubdivision(leftv res, leftv args)
     5084{
     5085  leftv h=args;
     5086  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5087  {
     5088     ideal h2= (ideal)h->Data();
     5089     h   = h->next;
     5090     if((h != NULL)&&(h->Typ() == POLY_CMD))
     5091     {
     5092       poly p= (poly)h->Data();
     5093       res->rtyp =IDEAL_CMD;
     5094       res->data =idMaken(stellarsub(p, h2));
     5095     }
     5096  }
     5097  return false;
     5098}
     5099
     5100
     5101
     5102BOOLEAN idModulo(leftv res, leftv args)
     5103{
     5104  leftv h=args;
     5105  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5106  {
     5107     ideal h1= (ideal)h->Data();
     5108     h   = h->next;
     5109     if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5110     {
     5111       ideal h2= (ideal)h->Data();
     5112       res->rtyp =IDEAL_CMD;
     5113       res->data =idmodulo(h1, h2);
     5114     }
     5115  }
     5116  return false;
     5117}
     5118
     5119
     5120BOOLEAN idminus(leftv res, leftv args)
     5121{
     5122  leftv h=args;
     5123  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5124  {
     5125     ideal h1= (ideal)h->Data();
     5126     h   = h->next;
     5127     if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5128     {
     5129       ideal h2= (ideal)h->Data();
     5130       res->rtyp =IDEAL_CMD;
     5131       res->data =idMinus(h1, h2);
     5132     }
     5133  }
     5134  return false;
     5135}
     5136
     5137
     5138
     5139BOOLEAN isoNumber(leftv res, leftv args)
     5140{
     5141  leftv h=args;
     5142  if((h != NULL)&&(h->Typ() == POLY_CMD))
     5143  {
     5144    poly p= (poly)h->Data();
     5145    h   = h->next;
     5146    if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5147    {
     5148      ideal h1= (ideal)h->Data();
     5149      h   = h->next;
     5150      if((h != NULL)&&(h->Typ() == POLY_CMD))
     5151      {
     5152        poly a= (poly)h->Data();
     5153        h   = h->next;
     5154        if((h != NULL)&&(h->Typ() == POLY_CMD))
     5155        {
     5156          poly b= (poly)h->Data();
     5157          res->rtyp =INT_CMD;
     5158          res->data =(void *)(long)isoNum(p, h1, a, b);
     5159        }
     5160      }
     5161    }
     5162  }
     5163  return false;
     5164}
     5165
     5166
     5167
     5168BOOLEAN ifIsomorphism(leftv res, leftv args)
     5169{
     5170  leftv h=args;
     5171  if((h != NULL)&&(h->Typ() == POLY_CMD))
     5172  {
     5173     poly p= (poly)h->Data();
     5174     h   = h->next;
     5175     if((h != NULL)&&(h->Typ() == POLY_CMD))
     5176     {
     5177       poly q= (poly)h->Data();
     5178       h   = h->next;
     5179       if((h != NULL)&&(h->Typ() == POLY_CMD))
     5180       {
     5181         poly f= (poly)h->Data();
     5182         h   = h->next;
     5183         if((h != NULL)&&(h->Typ() == POLY_CMD))
     5184         {
     5185           poly g= (poly)h->Data();
     5186           h   = h->next;
     5187           if((h != NULL)&&(h->Typ() == POLY_CMD))
     5188           {
     5189            poly a= (poly)h->Data();
     5190             h   = h->next;
     5191             if((h != NULL)&&(h->Typ() == POLY_CMD))
     5192             {
     5193               poly b= (poly)h->Data();
     5194               res->rtyp =INT_CMD;
     5195               res->data =(void *)(long)ifIso(p,q,f,g, a, b);
     5196             }
     5197           }
     5198         }
     5199       }
     5200     }
     5201  }
     5202  return false;
     5203}
     5204
     5205
     5206BOOLEAN newDegree(leftv res, leftv args)
     5207{
     5208  leftv h=args;
     5209  if((h != NULL)&&(h->Typ() == POLY_CMD))
     5210  {
     5211    poly p= (poly)h->Data();
     5212    h   = h->next;
     5213    if((h != NULL)&&(h->Typ() == INT_CMD))
     5214    {
     5215       int num= (int)(long)h->Data();
     5216       res->rtyp =INT_CMD;
     5217       res->data =(void *)(long)redefinedeg( p, num);
     5218    }
     5219  }
     5220  return false;
     5221}
     5222
     5223
     5224
     5225BOOLEAN nonf2f(leftv res, leftv args)
     5226{
     5227  leftv h=args;
     5228  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5229  {
     5230     ideal h1= (ideal)h->Data();
     5231     res->rtyp =IDEAL_CMD;
     5232     res->data =complementsimplex(h1);
     5233  }
     5234  return false;
     5235}
     5236
     5237
     5238
     5239BOOLEAN dimsim(leftv res, leftv args)
     5240{
     5241  leftv h=args;
     5242  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5243  {
     5244     ideal h1= (ideal)h->Data();
     5245     res->rtyp =INT_CMD;
     5246     res->data =(void *)(long)dim_sim(h1);
     5247  }
     5248  return false;
     5249}
     5250
     5251
     5252
     5253BOOLEAN numdim(leftv res, leftv args)
     5254{
     5255  leftv h=args;
     5256  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
     5257  {
     5258    ideal h1= (ideal)h->Data();
     5259    h   = h->next;
     5260    if((h != NULL)&&(h->Typ() == INT_CMD))
     5261    {
     5262       int num= (int)(long)h->Data();
     5263       res->rtyp =INT_CMD;
     5264       res->data =(void *)(long)num4dim( h1, num);
     5265    }
     5266  }
     5267  return false;
     5268}
    31175269
    31185270/**************************************interface T2****************************************/
    … …  
    31275279  p->iiAddCproc("","findaset",FALSE,fa);
    31285280  p->iiAddCproc("","fgp",FALSE,fgp);
    3129   p->iiAddCproc("","idcomplement",FALSE,idcomplement);
     5281  p->iiAddCproc("","fgpl",FALSE,fgpl);
     5282  p->iiAddCproc("","idcomplements",FALSE,idcomplement);
    31305283  p->iiAddCproc("","genst",FALSE,genstt);
    31315284  p->iiAddCproc("","sgp",FALSE,sgp);
     5285  p->iiAddCproc("","sgpl",FALSE,sgpl);
    31325286  p->iiAddCproc("","Links",FALSE,Links);
    3133 }
    3134 
    3135 
    3136 
    3137 extern "C" int SI_MOD_INIT(cohomo)(SModulFunctions* p)
     5287  p->iiAddCproc("","eqsolve1",FALSE,eqsolve1);
     5288  p->iiAddCproc("","pb",FALSE,pb);
     5289  p->iiAddCproc("","pa",FALSE,pa);
     5290  p->iiAddCproc("","makeSimplex",FALSE,makeSimplex);
     5291  p->iiAddCproc("","isSim",FALSE,isSim);
     5292  p->iiAddCproc("","nfaces1",FALSE,nfaces1);
     5293  p->iiAddCproc("","nfaces2",FALSE,nfaces2);
     5294  p->iiAddCproc("","nfaces3",FALSE,nfaces3);
     5295  p->iiAddCproc("","comedg",FALSE,comedg);
     5296  p->iiAddCproc("","tsets",FALSE,tsets);
     5297  p->iiAddCproc("","valency",FALSE,Valency);
     5298  p->iiAddCproc("","nab",FALSE,nabvl);
     5299  p->iiAddCproc("","tnab",FALSE,tnabvl);
     5300  p->iiAddCproc("","mab",FALSE,mabvl);
     5301  p->iiAddCproc("","SRideal",FALSE,SRideal);
     5302  p->iiAddCproc("","Linkn",FALSE,linkn);
     5303  p->iiAddCproc("","Existb",FALSE,existsub);
     5304  p->iiAddCproc("","pConstant",FALSE,pConstant);
     5305  p->iiAddCproc("","pChange",FALSE,pChange);
     5306  p->iiAddCproc("","pNew",FALSE,p_New);
     5307  p->iiAddCproc("","pSupport",FALSE,support);
     5308  p->iiAddCproc("","psMinusp",FALSE,psMinusp);
     5309  p->iiAddCproc("","cNew",FALSE,cNew);
     5310  p->iiAddCproc("","isoNumber",FALSE,isoNumber);
     5311  p->iiAddCproc("","vsInsec",FALSE,vsIntersec);
     5312  p->iiAddCproc("","getnabt",FALSE,nabtvl);
     5313  p->iiAddCproc("","idmodulo",FALSE,idModulo);
     5314  p->iiAddCproc("","ndegree",FALSE,newDegree);
     5315  p->iiAddCproc("","nonf2f",FALSE,nonf2f);
     5316  p->iiAddCproc("","ifIsom",FALSE,ifIsomorphism);
     5317  p->iiAddCproc("","stellarsubdivision",FALSE,stellarsubdivision);
     5318  p->iiAddCproc("","star",FALSE,stars);
     5319  p->iiAddCproc("","numdim",FALSE,numdim);
     5320  p->iiAddCproc("","dimsim",FALSE,dimsim);
     5321  p->iiAddCproc("","bprime",FALSE,bprime);
     5322  p->iiAddCproc("","remainpart",FALSE,stellarremain);
     5323  p->iiAddCproc("","idminus",FALSE,idminus);
     5324  p->iiAddCproc("","time1",FALSE,t1h);
     5325
     5326}
     5327
     5328
     5329
     5330extern "C" int SI_MOD_INIT0(cohomo)(SModulFunctions* p)
    31385331{
    31395332  firstorderdef_setup(p);

  • gfanlib/Makefile.am

    r5c73a12 rff2859d  
    1717SOURCES = gfanlib_circuittableint.cpp gfanlib_mixedvolume.cpp gfanlib_paralleltraverser.cpp gfanlib_polyhedralfan.cpp gfanlib_polymakefile.cpp gfanlib_symmetriccomplex.cpp gfanlib_symmetry.cpp gfanlib_traversal.cpp gfanlib_zcone.cpp gfanlib_zfan.cpp
    1818libgfan_la_SOURCES = $(SOURCES)
    19 libgfan_la_LDFLAGS = $(SINGULAR_LDFLAGS) $(CDDGMPLDFLAGS)
     19libgfan_la_LDFLAGS = $(SINGULAR_LDFLAGS) $(CDDGMPLDFLAGS) -release ${PACKAGE_VERSION}
    2020
    2121libgfan_includedir =$(includedir)/gfanlib
Note: See TracChangeset for help on using the changeset viewer.