My Project
cohomo.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /*
5 * ABSTRACT - Stainly
6 */
7 
8 #include "kernel/mod2.h"
9 
10 #if !defined(__CYGWIN__) || defined(STATIC_VERSION)
11 // acces from a module to routines from the main program
12 // does not work on windows (restrict of the dynamic linker),
13 // a static version is required:
14 // ./configure --with-builtinmodules=cohomo,...
15 
16 
17 #include "omalloc/omalloc.h"
18 #include "misc/mylimits.h"
19 #include "libpolys/misc/intvec.h"
20 #include <assert.h>
21 #include <unistd.h>
22 
26 #include "cohomo.h"//for my thing
27 #include "kernel/GBEngine/tgb.h"//
28 #include "Singular/ipid.h"//for ggetid
29 #include "polys/monomials/ring.h"
31 #include "polys/simpleideals.h"
32 #include "Singular/lists.h"
33 #include "kernel/linear_algebra/linearAlgebra.h"//for printNumber
34 #include "kernel/GBEngine/kstd1.h"//for gb
35 #include <kernel/ideals.h>
36 #if 1
37 
41 #include <coeffs/numbers.h>
42 #include <vector>
43 #include <Singular/ipshell.h>
44 #include <Singular/libsingular.h>
45 #include <time.h>
46 
47 /***************************print(only for debugging)***********************************************/
48 //print vector of integers.
49 static void listprint(std::vector<int> vec)
50 {
51  unsigned i;
52  for(i=0;i<vec.size();i++)
53  {
54  Print(" _[%d]=%d\n",i+1,vec[i]);
55  PrintLn();
56  }
57  if(vec.size()==0)
58  {
59  PrintS(" _[1]= \n");
60  PrintLn();
61  }
62 }
63 
64 //print vector of vectors of integers.
65 static void listsprint(std::vector<std::vector<int> > posMat)
66 {
67  unsigned i;
68  for(i=0;i<posMat.size();i++)
69  {
70  Print("[%d]:\n",i+1);
71  listprint(posMat[i]);
72  Print("\n");
73  PrintLn();
74  }
75  if(posMat.size()==0)
76  {
77  PrintS("[1]:\n");
78  PrintLn();
79  }
80 }
81 
82 
83 //print ideal.
84 static void id_print(ideal h)
85 {
86  int i;
87  for(i=0;i<IDELEMS(h);i++)
88  {
89  Print(" [%d]\n",i+1);
90  pWrite(h->m[i]);
91  PrintLn();
92  }
93 }
94 
95 //only for T^2,
96 //print vector of polynomials.
97 static void lpprint( std::vector<poly> pv)
98 {
99  for(unsigned i=0;i<pv.size();i++)
100  {
101  Print(" _[%d]=",i+1);
102  pWrite(pv[i]);
103  }
104  if(pv.size()==0)
105  {
106  PrintS(" _[1]= \n");
107  PrintLn();
108  }
109 }
110 
111 //print vector of vectors of polynomials.
112 static void lpsprint(std::vector<std::vector<poly> > pvs)
113 {
114  for(unsigned i=0;i<pvs.size();i++)
115  {
116  Print("[%d]:\n",i+1);
117  lpprint(pvs[i]);
118  Print("\n");
119  PrintLn();
120  }
121  if(pvs.size()==0)
122  {
123  PrintS("[1]:\n");
124  PrintLn();
125  }
126 }
127 
128 /*************operations for vectors (regard vectors as sets)*********/
129 
130 //returns true if integer n is in vector vec,
131 //otherwise, returns false
132 static bool IsinL(int a, std::vector<int> vec)
133 {
134  unsigned i;
135  for(i=0;i<vec.size();i++)
136  {
137  if(a==vec[i])
138  {
139  return true;
140  }
141  }
142  return false;
143 }
144 
145 //returns intersection of vectors p and q,
146 //returns empty if they are disjoint
147 static std::vector<int> vecIntersection(std::vector<int> p, std::vector<int> q)
148 {
149  unsigned i;
150  std::vector<int> inte;
151  for(i=0;i<p.size();i++)
152  {
153  if(IsinL(p[i],q))
154  inte.push_back(p[i]);
155  }
156  return inte;
157 }
158 
159 //returns true if vec1 is contained in vec2
160 static bool vsubset(std::vector<int> vec1, std::vector<int> vec2)
161 {
162  int i;
163  if(vec1.size()>vec2.size())
164  return false;
165  for(i=0;i<vec1.size();i++)
166  {
167  if(!IsinL(vec1[i],vec2))
168  return false;
169  }
170  return true;
171 }
172 
173 //not strictly equal(order doesn't matter)
174 static bool vEvl(std::vector<int> vec1, std::vector<int> vec2)
175 {
176  if(vec1.size()==0 && vec2.size()==0)
177  return true;
178  if(vsubset(vec1,vec2)&&vsubset(vec2,vec1))
179  return true;
180  return false;
181 }
182 
183 //the length of vec must be same to it of the elements of vecs
184 //returns true if vec is as same as some element of vecs ((not strictly same))
185 //returns false if vec is not in vecs
186 static bool vInvsl(std::vector<int> vec, std::vector<std::vector<int> > vecs)
187 {
188  int i;
189  for(i=0;i<vecs.size();i++)
190  {
191  if(vEvl(vec,vecs[i]))
192  {
193  return true;
194  }
195  }
196  return false;
197 }
198 
199 //returns the union of two vectors(as the union of sets)
200 static std::vector<int> vecUnion(std::vector<int> vec1, std::vector<int> vec2)
201 {
202  std::vector<int> vec=vec1;
203  unsigned i;
204  for(i=0;i<vec2.size();i++)
205  {
206  if(!IsinL(vec2[i],vec))
207  vec.push_back(vec2[i]);
208  }
209  return vec;
210 }
211 
212 static std::vector<int> vecMinus(std::vector<int> vec1,std::vector<int> vec2)
213 {
214  std::vector<int> vec;
215  for(unsigned i=0;i<vec1.size();i++)
216  {
217  if(!IsinL(vec1[i],vec2))
218  {
219  vec.push_back(vec1[i]);
220  }
221  }
222  return vec;
223 }
224 
225 static std::vector<std::vector<int> > vsMinusv(std::vector<std::vector<int> > vecs, std::vector<int> vec)
226 {
227  int i;
228  std::vector<std::vector<int> > rem;
229  for(i=0;i<vecs.size();i++)
230  {
231  if(!vEvl(vecs[i],vec))
232  {
233  rem.push_back(vecs[i]);
234  }
235  }
236  return (rem);
237 }
238 
239 static std::vector<std::vector<int> > vsUnion(std::vector<std::vector<int> > vs1, std::vector<std::vector<int> > vs2)
240 {
241  int i;
242  std::vector<std::vector<int> > vs=vs1;
243  for(i=0;i<vs2.size();i++)
244  {
245  if(!vInvsl(vs2[i],vs))
246  {
247  vs.push_back(vs2[i]);
248  }
249  }
250  return vs;
251 }
252 
253 static std::vector<std::vector<int> > vsIntersection(std::vector<std::vector<int> > vs1, std::vector<std::vector<int> > vs2)
254 {
255  int i;
256  std::vector<std::vector<int> > vs;
257  for(i=0;i<vs2.size();i++)
258  {
259  if(vInvsl(vs2[i],vs1))
260  {
261  vs.push_back(vs2[i]);
262  }
263  }
264  return vs;
265 }
266 
267 /*************************************for transition between ideal and vectors******************************************/
268 
269 //P should be monomial,
270 // vector version of poly support(poly p)
271 static std::vector<int> support1(poly p)
272 {
273  int j;
274  std::vector<int> supset;
275  if(p==0) return supset;
276  for(j=1;j<=rVar(currRing);j++)
277  {
278  if(pGetExp(p,j)>0)
279  {
280  supset.push_back(j);
281  }
282  }
283  return (supset);
284 }
285 
286 //simplicial complex(the faces set is ideal h)
287 static std::vector<std::vector<int> > supports(ideal h)
288 {
289  std::vector<std::vector<int> > vecs;
290  std::vector<int> vec;
291  if(!idIs0(h))
292  {
293  for(int s=0;s<IDELEMS(h);s++)
294  {
295  vec=support1(h->m[s]);
296  vecs.push_back(vec);
297  }
298  }
299  return vecs;
300 }
301 
302 // only for eqsolve1
303 // p could be any polynomial
304 static std::vector<int> support2(poly p)
305 {
306  int j;
307  poly q;
308  std::vector<int> supset;
309  for(j=1;j<=rVar(currRing);j++)
310  {
311  q=pCopy(p);
312  while (q!=NULL)
313  {
314  if(p_GetExp(q,j,currRing)!=0)
315  {
316  supset.push_back(j);
317  break;
318  }
319  q=pNext(q);
320  }
321  }
322  return (supset);
323 }
324 
325 //the supports of ideal
326 static std::vector<std::vector<int> > supports2(ideal h)
327 {
328  std::vector<std::vector<int> > vecs;
329  std::vector<int> vec;
330  if(!idIs0(h))
331  {
332  for(int s=0;s<IDELEMS(h);s++)
333  {
334  vec=support2(h->m[s]);
335  vecs.push_back(vec);
336  }
337  }
338  return vecs;
339 }
340 
341 //convert the vector(vbase[i] are the coefficients of x_{i+1}) to a polynomial w.r.t. current ring
342 //vector vbase has length of currRing->N.
343 static poly pMake(std::vector<int> vbase)
344 {
345  int n=vbase.size(); poly p,q=0;
346  for(int i=0;i<n;i++)
347  {
348  if(vbase[i]!=0)
349  {
350  p = pOne();pSetExp(p, i+1, 1);pSetm(p);pSetCoeff(p, nInit(vbase[i]));
351  q = pAdd(q, p);
352  }
353  }
354  return q;
355 }
356 
357 //convert the vectors to a ideal(for T^1)
358 static ideal idMake(std::vector<std::vector<int> > vecs)
359 {
360  int lv=vecs.size(), i;
361  poly p;
362  ideal id_re=idInit(1,1);
363  for(i=0;i<lv;i++)
364  {
365  p=pMake(vecs[i]);
366  idInsertPoly(id_re, p);
367  }
368  idSkipZeroes(id_re);
369  return id_re;
370 }
371 
372 /*****************************quotient ring of two ideals*********************/
373 
374 //the quotient ring of h1 respect to h2
375 static ideal idmodulo(ideal h1,ideal h2)
376 {
377  ideal gb=kStd(h2,NULL,testHomog,NULL,NULL,0,0,NULL);
378  idSkipZeroes(gb);
379  ideal idq=kNF(gb,NULL,h1);
380  idSkipZeroes(idq);
381  return idq;
382 }
383 
384 //returns the coeff of the monomial of polynomial p which involves the mth varialbe
385 //assume the polynomial p has form of y1+y2+...
386 static int pcoef(poly p, int m)
387 {
388  int i,co; poly q=pCopy(p);
389  for(i=1;i<=currRing->N;i++)
390  {
391  if(p_GetExp(q,m,currRing)!=0)
392  {
393  co=n_Int(pGetCoeff(q),currRing->cf);
394  return co;
395  }
396  else
397  q=pNext(q);
398  }
399  if(q!=NULL)
400  co=0;
401  return co;
402 }
403 
404 //returns true if p involves the mth variable of the current ring
405 static bool vInp(int m,poly p)
406 {
407  poly q=pCopy(p);
408  while (q!=NULL)
409  {
410  if(p_GetExp(q,m,currRing)!=0)
411  {
412  return true;
413  }
414  q=pNext(q);
415  }
416  return false;
417 }
418 
419 //returns the vector w.r.t. polynomial p
420 static std::vector<int> vMake(poly p)
421 {
422  int i;
423  std::vector<int> vbase;
424  for(i=1;i<=currRing->N;i++)
425  {
426  if(vInp(i,p))
427  {
428  vbase.push_back(pcoef(p,i));
429  }
430  else
431  {
432  vbase.push_back(0);
433  }
434  }
435  return (vbase);
436 }
437 
438 //returns the vectors w.r.t. ideal h
439 static std::vector<std::vector<int> > vsMake(ideal h)
440 {
441  std::vector<int> vec;
442  std::vector<std::vector<int> > vecs;
443  int i;
444  for(i=0;i<IDELEMS(h);i++)
445  {
446  vec=vMake(h->m[i]);
447  vecs.push_back(vec);
448  }
449  return vecs;
450 }
451 
452 //the quotient ring of two ideals which are represented by vectors,
453 //the result is also represented by vector.
454 static std::vector<std::vector<int> > vecqring(std::vector<std::vector<int> > vec1, std::vector<std::vector<int> > vec2)
455 {
456  ideal h1=idMake(vec1), h2=idMake(vec2);
457  ideal h=idmodulo(h1,h2);
458  std::vector<std::vector<int> > vecs= vsMake(h);
459  return vecs;
460 }
461 
462 /****************************************************************/
463 //construct a monomial based on the support of it
464 //returns a squarefree monomial
465 static poly pMaken(std::vector<int> vbase)
466 {
467  int n=vbase.size();
468  poly p,q=pOne();
469  for(int i=0;i<n;i++)
470  {
471  p = pOne();pSetExp(p, vbase[i], 1);pSetm(p);pSetCoeff(p, nInit(1));
472  //pWrite(p);
473  q=pp_Mult_mm(q,p,currRing);
474  }
475  return q;
476 }
477 
478 // returns a ideal according to a set of supports
479 static ideal idMaken(std::vector<std::vector<int> > vecs)
480 {
481  ideal id_re=idInit(1,1);
482  poly p;
483  int i,lv=vecs.size();
484  for(i=0;i<lv;i++)
485  {
486  p=pMaken(vecs[i]);
487  idInsertPoly(id_re, p);
488  }
489  idSkipZeroes(id_re);
490  //id_print(id_re);
491  return id_re;
492 }
493 
494 /********************************new version for stanley reisner ideal ***********************************************/
495 
496 static std::vector<std::vector<int> > b_subsets(std::vector<int> vec)
497 {
498  int i,j;
499  std::vector<int> bv;
500  std::vector<std::vector<int> > vecs;
501  for(i=0;i<vec.size();i++)
502  {
503  bv.push_back(vec[i]);
504  vecs.push_back(bv);
505  bv.clear();
506  }
507  //listsprint(vecs);
508  for(i=0;i<vecs.size();i++)
509  {
510  for(j=i+1;j<vecs.size();j++)
511  {
512  bv=vecUnion(vecs[i], vecs[j]);
513  if(!vInvsl(bv,vecs))
514  vecs.push_back(bv);
515  }
516  }
517  //listsprint(vecs);
518  return(vecs);
519 }
520 
521 //the number of the variables
522 static int idvert(ideal h)
523 {
524  int i, j, vert=0;
525  if(idIs0(h))
526  return vert;
527  for(i=currRing->N;i>0;i--)
528  {
529  for(j=0;j<IDELEMS(h);j++)
530  {
531  if(pGetExp(h->m[j],i)>0)
532  {
533  vert=i;
534  return vert;
535  }
536  }
537  }
538  return vert;
539 }
540 
541 static int pvert(poly p)
542 {
543  int i, vert=0;
544  for(i=currRing->N;i>0;i--)
545  {
546  if(pGetExp(p,i)>0)
547  {
548  vert=i;
549  return vert;
550  }
551  }
552  return vert;
553 }
554 
555 /*
556 //full complex
557 static std::vector<std::vector<int> > fullcomplex(ideal h)
558 {
559  int vert=vertnum(h), i, j;
560  //Print("there are %d vertices\n", vert);
561  std::vector<std::vector<int> > fmons;
562  std::vector<int> pre;
563  for(i=1;i<=vert;i++)
564  {
565  pre.push_back(i);
566  }
567  fmons=b_subsets(pre);
568  return fmons;
569 }*/
570 
571 /*
572 //all the squarefree monomials whose degree is less or equal to n
573 static std::vector<std::vector<int> > sfrmons(ideal h, int n)
574 {
575  int vert=vertnum(h), i, j, time=0;
576  std::vector<std::vector<int> > fmons, pres, pres0, pres1;
577  std::vector<int> pre;
578  for(i=1;i<=vert;i++)
579  {
580  pre.push_back(i);
581  pres0.push_back(pre);
582  }
583  pres=pres0;
584  for(i=0;i<size(pres),time<=n;i++)
585  {
586  time++;
587  pre=pres[i];
588  for(j=0;j<size(pres0);j++)
589  {
590  pre=vecUnion(pre, pres0[j]);
591  if(pre.)
592  }
593  }
594  return fmons;
595 }
596 */
597 
598 /*
599 static ideal id_complement(ideal h)
600 {
601  int i,j;
602  std::vector<std::vector<int> > full=fullcomplex(h), hvs=supports(h), res;
603  for(i=0;i<full.size();i++)
604  {
605  if(!vInvsl(full[i], hvs))
606  {
607  res.push_back(full[i]);
608  }
609  }
610  return idMaken(res);
611 }*/
612 
613 
614 /*****************About simplicial complex and stanley-reisner ideal and ring **************************/
615 
616 //h1 minus h2
617 static ideal idMinus(ideal h1,ideal h2)
618 {
619  ideal h=idInit(1,1);
620  int i,j,eq=0;
621  for(i=0;i<IDELEMS(h1);i++)
622  {
623  eq=0;
624  for(j=0;j<IDELEMS(h2);j++)
625  {
626  if(p_EqualPolys(pCopy(h1->m[i]),pCopy(h2->m[j]), currRing))
627  {
628  eq=1;
629  break;
630  }
631  }
632  if(eq==0)
633  {
634  idInsertPoly(h, pCopy(h1->m[i]));
635  }
636  }
637  idSkipZeroes(h);
638  return h;
639 }
640 
641 //If poly P is squarefree, returns 1
642 //returns 0 otherwise,
643 static bool p_Ifsfree(poly P)
644 {
645  int i,sf=1;
646  for(i=1;i<=rVar(currRing);i++)
647  {
648  if (pGetExp(P,i)>1)
649  {
650  sf=0;
651  break;
652  }
653  }
654  return sf;
655 }
656 
657 //returns the set of all squarefree monomials of degree deg in ideal h
658 static ideal sfreemon(ideal h,int deg)
659 {
660  int j;
661  ideal temp;
662  temp=idInit(1,1);
663  if(!idIs0(h))
664  {
665  for(j=0;j<IDELEMS(h);j++)
666  {
667  if((p_Ifsfree(h->m[j]))&&(pTotaldegree(h->m[j])==deg))
668  {
669  idInsertPoly(temp, h->m[j]);
670  }
671  }
672  idSkipZeroes(temp);
673  }
674  return temp;
675 }
676 
677 //full simplex represented by ideal.
678 //(all the squarefree monomials over the polynomial ring)
679 static ideal id_sfmon(ideal h)
680 {
681  ideal asfmons,sfmons,mons;
682  int j, vert=idvert(h);
683  mons=id_MaxIdeal(1, currRing);
684  asfmons=sfreemon(mons,1);
685  for(j=2;j<=vert;j++)
686  {
687  mons=id_MaxIdeal(j, currRing);
688  sfmons=sfreemon(mons,j);
689  asfmons=id_Add(asfmons,sfmons,currRing);
690  }
691  return asfmons;
692 }
693 
694 //if the input ideal is simplicial complex, returns the stanley-reisner ideal,
695 //if the input ideal is stanley-reisner ideal, returns the monomial ideal according to simplicial complex.
696 //(nonfaces and faces).
697 //returns the complement of the ideal h (consisting of only squarefree polynomials)
698 static ideal id_complement(ideal h)
699 {
700  int j, vert=idvert(h);
701  ideal i1=id_sfmon(h);
702  ideal i3=idInit(1,1);
703  poly p;
704  for(j=0;j<IDELEMS(i1);j++)
705  {
706  p=pCopy(i1->m[j]);
707  if(pvert(p)<=vert)
708  {
709  idInsertPoly(i3, p);
710  }
711  }
712  ideal i2=idMinus(i3,h);
713  idSkipZeroes(i2);
714  return (i2);
715 }
716 
717 //Returns true if p is one of the generators of ideal X
718 //returns false otherwise
719 static bool IsInX(poly p,ideal X)
720 {
721  int i;
722  for(i=0;i<IDELEMS(X);i++)
723  {
724  if(pEqualPolys(p,X->m[i]))
725  {
726  //PrintS("yes\n");
727  return(true);
728  }
729  }
730  //PrintS("no\n");
731  return(false);
732 }
733 
734 //returns the monomials in the quotient ring R/(h1+h2) which have degree deg.
735 static ideal qringadd(ideal h1, ideal h2, int deg)
736 {
737  ideal h,qrh;
738  h=idAdd(h1,h2);
739  qrh=scKBase(deg,h);
740  return qrh;
741 }
742 
743 //returns the maximal degree of the monomials in ideal h
744 static int id_maxdeg(ideal h)
745 {
746  int i,max;
747  max=pTotaldegree(h->m[0]);
748  for(i=1;i<IDELEMS(h);i++)
749  {
750  if(pTotaldegree(h->m[i]) > max)
751  max=pTotaldegree(h->m[i]);
752  }
753  return (max);
754 }
755 
756 //input ideal h (a squarefree monomial ideal) is the ideal associated to simplicial complex,
757 //and returns the Stanley-Reisner ideal(minimal generators)
758 static ideal idsrRing(ideal h)
759 {
760  int i,n;
761  ideal pp,qq,rsr,ppp,hc=idCopy(h);
762  for(i=1;i<=rVar(currRing);i++)
763  {
764  pp=sfreemon(hc,i);
765  pp=scKBase(i,pp);//quotient ring (R/I_i)_i
766  if(!idIs0(pp))
767  {
768  pp=sfreemon(pp,i);
769  rsr=pp;
770  //Print("This is the first quotient generators %d:\n",i);
771  //id_print(rsr);
772  break;
773  }
774  }
775  for(n=i+1;n<=rVar(currRing);n++)
776  {
777  qq=sfreemon(hc,n);
778  pp=qringadd(qq,rsr,n);
779  ppp=sfreemon(pp,n);
780  rsr=idAdd(rsr,ppp);
781  }
782  idSkipZeroes(rsr);
783  return rsr;
784 }
785 
786 //returns the set of all the polynomials could divide p
787 static ideal SimFacset(poly p)
788 {
789  int i,j,max=pTotaldegree(p);
790  ideal h1,mons,id_re=idInit(1,1);
791  for(i=1;i<max;i++)
792  {
793  mons=id_MaxIdeal(i, currRing);
794  h1=sfreemon(mons,i);
795 
796  for(j=0;j<IDELEMS(h1);j++)
797  {
798  if(p_DivisibleBy(h1->m[j],p,currRing))
799  {
800  idInsertPoly(id_re, h1->m[j]);
801  }
802  }
803  }
804  idSkipZeroes(id_re);
805  return id_re;
806 }
807 
808 static ideal idadda(ideal h1, ideal h2)
809 {
810  ideal h=idInit(1,1);
811  for(int i=0;i<IDELEMS(h1);i++)
812  {
813  if(!IsInX(h1->m[i],h))
814  {
815  idInsertPoly(h, h1->m[i]);
816  }
817  }
818  for(int i=0;i<IDELEMS(h2);i++)
819  {
820  if(!IsInX(h2->m[i],h))
821  {
822  idInsertPoly(h, h2->m[i]);
823  }
824  }
825  idSkipZeroes(h);
826  return h;
827 }
828 
829 //complicated version
830 //(returns false if it is not a simplicial complex and print the simplex)
831 //input h is need to be at least part of faces
832 static ideal IsSimplex(ideal h)
833 {
834  int i,max=id_maxdeg(h);
835  poly e=pOne();
836  ideal id_re, id_so=idCopy(h);
837  for(i=0;i<IDELEMS(h);i++)
838  {
839  id_re=SimFacset(h->m[i]);
840  if(!idIs0(id_re))
841  {
842  id_so=idadda(id_so, id_re);//idAdd(id_so,id_re);
843  }
844  }
845  idInsertPoly(id_so,e);
846  idSkipZeroes(id_so);
847  return (idMinus(id_so,h));
848 }
849 
850 //input is the subset of the Stainley-Reisner ideal
851 //returns the faces
852 //is not used
853 static ideal complementsimplex(ideal h)
854 {
855  int i,j;poly p,e=pOne();
856  ideal h1=idInit(1,1), pp, h3;
857  for(i=1;i<=rVar(currRing);i++)
858  {
859  p = pOne(); pSetExp(p, i, 2); pSetm(p); pSetCoeff(p, nInit(1));
860  idInsertPoly(h1, p);
861  }
862  idSkipZeroes(h1);
863  ideal h2=idAdd(h,h1);
864  pp=scKBase(1,h2);
865  h3=idCopy(pp);
866  for(j=2;j<=rVar(currRing);j++)
867  {
868  pp=scKBase(j,h2);
869  h3=idAdd(h3,pp);
870  }
871  idInsertPoly(h3, e);
872  idSkipZeroes(h3);
873  return (h3);
874 }
875 
876 static int dim_sim(ideal h)
877 {
878  int dim=pTotaldegree(h->m[0]), i;
879  for(i=1; i<IDELEMS(h);i++)
880  {
881  if(dim<pTotaldegree(h->m[i]))
882  {
883  dim=pTotaldegree(h->m[i]);
884  }
885  }
886  return dim;
887 }
888 
889 static int num4dim(ideal h, int n)
890 {
891  int num=0;
892  for(int i=0; i<IDELEMS(h); i++)
893  {
894  if(pTotaldegree(h->m[i])==n)
895  {
896  num++;
897  }
898  }
899  return num;
900 }
901 
902 /********************Procedures for T1(M method and N method) ***********/
903 
904 //h is ideal( monomial ideal) associated to simplicial complex
905 //returns the all the monomials x^b (x^b must be able to divide
906 //at least one monomial in Stanley-Reisner ring)
907 //not so efficient
908 static ideal findb(ideal h)
909 {
910  ideal ib=id_sfmon(h), nonf=id_complement(h), bset=idInit(1,1);
911  poly e=pOne();
912  int i,j;
913  for(i=0;i<IDELEMS(ib);i++)
914  {
915  for(j=0;j<IDELEMS(nonf);j++)
916  {
917  if(p_DivisibleBy(ib->m[i],nonf->m[j],currRing))
918  {
919  idInsertPoly(bset, ib->m[i]);
920  break;
921  }
922  }
923  }
924  idInsertPoly(bset,e);
925  idSkipZeroes(bset);
926  return bset;
927 }
928 
929 //h is ideal(monomial ideal associated to simplicial complex
930 //1.poly S is x^b
931 //2.and deg(x^a)=deg(x^b)
932 //3.x^a and x^a have disjoint supports
933 //returns all the possible x^a according conditions 1. 2. 3.
934 static ideal finda(ideal h,poly S,int ddeg)
935 {
936  poly e=pOne();
937  ideal h2=id_complement(h), aset=idInit(1,1);
938  int i,deg1=pTotaldegree(S);
939  int tdeg=deg1+ddeg;
940  if(tdeg!=0)
941  {
942  std::vector<int> v,bv=support1(S),in;
943  std::vector<std::vector<int> > hvs=supports(h);
944  ideal ia=id_MaxIdeal(tdeg, currRing);
945  for(i=0;i<IDELEMS(ia);i++)
946  {
947  v=support1(ia->m[i]);
948  in=vecIntersection(v,bv);
949  if(vInvsl(v,hvs)&&in.size()==0)
950  {
951  idInsertPoly(aset, ia->m[i]);
952  }
953  }
954  idSkipZeroes(aset);
955  }
956  else idInsertPoly(aset,e);
957  return(aset);
958 }
959 
960 //returns true if support(p) union support(a) minus support(b) is face,
961 //otherwise returns false
962 //(the vector version of mabcondition)
963 static bool mabconditionv(std::vector<std::vector<int> > hvs,std::vector<int> pv,std::vector<int> av,std::vector<int> bv)
964 {
965  std::vector<int> uv=vecUnion(pv,av);
966  uv=vecMinus(uv,bv);
967  if(vInvsl(uv,hvs))
968  {
969  return(true);
970  }
971  return(false);
972 }
973 
974 // returns the set of nonfaces p where mabconditionv(h, p, a, b) is true
975 static std::vector<std::vector<int> > Mabv(ideal h,poly a,poly b)
976 {
977  std::vector<int> av=support1(a), bv=support1(b), pv, vec;
978  ideal h2=id_complement(h);
979  std::vector<std::vector<int> > hvs=supports(h), h2v=supports(h2), vecs;
980  for(unsigned i=0;i<h2v.size();i++)
981  {
982  pv=h2v[i];
983  if(mabconditionv(hvs,pv,av,bv))
984  {
985  vecs.push_back(pv);
986  }
987  }
988  return vecs;
989 }
990 
991 /***************************************************************************/
992 //For solving the equations which has form of x_i-x_j.(equations got from T_1)
993 /***************************************************************************/
994 
995 //subroutine for soleli1
996 static std::vector<int> eli1(std::vector<int> eq1,std::vector<int> eq2)
997 {
998  int i,j;
999  std::vector<int> eq;
1000  if(eq1[0]==eq2[0])
1001  {
1002  i=eq1[1];j=eq2[1];
1003  eq.push_back(i);
1004  eq.push_back(j);
1005  }
1006  else
1007  {
1008  eq=eq2;
1009  }
1010  return(eq);
1011 }
1012 
1013 /*
1014 //get triangular form(eqs.size()>0)
1015 static std::vector<std::vector<int> > soleli1( std::vector<std::vector<int> > eqs)
1016 {
1017  int i,j;
1018  std::vector<int> yaya;
1019  std::vector<std::vector<int> > pre=eqs, ppre, re;
1020  if(eqs.size()>0)
1021  {
1022  re.push_back(eqs[0]);
1023  pre.erase(pre.begin());
1024  }
1025  for(i=0;i<re.size(),pre.size()>0;i++)
1026  {
1027  yaya=eli1(re[i],pre[0]);
1028  re.push_back(yaya);
1029  for(j=1;j<pre.size();j++)
1030  {
1031  ppre.push_back(eli1(re[i],pre[j]));
1032  }
1033  pre=ppre;
1034  ppre.resize(0);
1035  }
1036  return re;
1037 }*/
1038 //make sure the first element is smaller that the second one
1039 static std::vector<int> keeporder( std::vector<int> vec)
1040 {
1041  std::vector<int> yaya;
1042  int n;
1043  if(vec[0]>vec[1])
1044  {
1045  n=vec[0];
1046  vec[0]=vec[1];
1047  vec[1]=n;
1048  }
1049  return vec;
1050 }
1051 
1052 static std::vector<std::vector<int> > soleli1( std::vector<std::vector<int> > eqs)
1053 {
1054  int i;
1055  std::vector<int> yaya;
1056  std::vector<std::vector<int> > pre=eqs, ppre, re;
1057  if(eqs.size()>0)
1058  {
1059  re.push_back(eqs[0]);
1060  pre.erase(pre.begin());
1061  }
1062  while(pre.size()>0)
1063  {
1064  yaya=keeporder(eli1(re[0],pre[0]));
1065  for(i=1;i<re.size();i++)
1066  {
1067  if(!vInvsl(yaya, re))
1068  {
1069  yaya=eli1(re[i],yaya);
1070  yaya=keeporder(yaya);
1071  }
1072  }
1073  if(!vInvsl(yaya, re))
1074  {
1075  re.push_back(yaya);
1076  }
1077  pre.erase(pre.begin());
1078  }
1079  return re;
1080 }
1081 
1082 // input is a set of equations who is of triangular form(every equations has a form of x_i-x_j)
1083 // n is the number of variables
1084 //get the free variables and the dimension
1085 static std::vector<int> freevars(int n, std::vector<int> bset, std::vector<std::vector<int> > gset)
1086 {
1087  int ql=gset.size(), bl=bset.size(), i;
1088  std::vector<int> mvar, fvar;
1089  for(i=0;i<bl;i++)
1090  {
1091  mvar.push_back(bset[i]);
1092  }
1093  for(i=0;i<ql;i++)
1094  {
1095  mvar.push_back(gset[i][0]);
1096  }
1097  for(i=1;i<=n;i++)
1098  {
1099  if(!IsinL(i,mvar))
1100  {
1101  fvar.push_back(i);
1102  }
1103  }
1104  return fvar;
1105 }
1106 
1107 //return the set of free variables except the vnum one
1108 static std::vector<int> fvarsvalue(int vnum, std::vector<int> fvars)
1109 {
1110  int i;
1111  std::vector<int> fset=fvars;
1112  for(i=0;i<fset.size();i++)
1113  {
1114  if(fset[i]==vnum)
1115  {
1116  fset.erase(fset.begin()+i);
1117  break;
1118  }
1119  }
1120  return fset;
1121 }
1122 
1123 //returns the simplified bset and gset
1124 //enlarge bset, simplify gset
1125 static std::vector<std::vector<int> > vAbsorb( std::vector<int> bset,std::vector<std::vector<int> > gset)
1126 {
1127  std::vector<int> badset=bset;
1128  int i,j,m, bl=bset.size(), gl=gset.size();
1129  for(i=0;i<bl;i++)
1130  {
1131  m=badset[i];
1132  for(j=0;j<gl;j++)
1133  {
1134  if(gset[j][0]==m && !IsinL(gset[j][1],badset))
1135  {
1136  badset.push_back(gset[j][1]);
1137  gset.erase(gset.begin()+j);
1138  j--;
1139  gl--;
1140  bl++;
1141  }
1142  else if(!IsinL(gset[j][0],badset) && gset[j][1]==m)
1143  {
1144  badset.push_back(gset[j][0]);
1145  gset.erase(gset.begin()+j);
1146  j--;
1147  gl--;
1148  bl++;
1149  }
1150  else if(IsinL(gset[j][0],badset) && IsinL(gset[j][1],badset))
1151  {
1152  gset.erase(gset.begin()+j);
1153  j--;
1154  gl--;
1155  }
1156  else
1157  {
1158  ;
1159  }
1160  }
1161  }
1162  if(badset.size()==0) badset.push_back(0);
1163  gset.push_back(badset);
1164  return gset;
1165 }
1166 
1167 //the labels of new variables are started with 1
1168 //returns a vector of solution space according to index
1169 static std::vector<int> vecbase1(int num, std::vector<int> oset)
1170 {
1171  int i;
1172  std::vector<int> base;
1173  for(i=0;i<num;i++)
1174  {
1175  if(IsinL(i+1,oset))
1176  base.push_back(1);
1177  else
1178  base.push_back(0);
1179  }
1180  return base;
1181 }
1182 
1183 //returns a vector which has length of n,
1184 //and all the entries are 0.
1185 static std::vector<int> make0(int n)
1186 {
1187  int i;
1188  std::vector<int> vec;
1189  for(i=0;i<n;i++)
1190  {
1191  vec.push_back(0);
1192  }
1193  return vec;
1194 }
1195 
1196 //returns a vector which has length of n,
1197 //and all the entries are 1.
1198 static std::vector<int> make1(int n)
1199 {
1200  int i;
1201  std::vector<int> vec;
1202  for(i=0;i<n;i++)
1203  {
1204  vec.push_back(1);
1205  }
1206  return vec;
1207 }
1208 
1209 //input gset must be the triangular form after zero absorbing according to the badset,
1210 //bset must be the zero set after absorbing.
1211 static std::vector<int> ofindbases1(int num, int vnum, std::vector<int> bset,std::vector<std::vector<int> > gset)
1212 {
1213  std::vector<std::vector<int> > goodset;
1214  std::vector<int> fvars=freevars(num, bset, gset), oset, base;
1215  std::vector<int> zset=fvarsvalue(vnum, fvars);
1216  zset=vecUnion(zset,bset);
1217  oset.push_back(vnum);
1218  goodset=vAbsorb(oset, gset);
1219  oset=goodset[goodset.size()-1];
1220  goodset.erase(goodset.end());
1221  base= vecbase1(num, oset);
1222  return base;
1223 }
1224 
1225 //input gset must be the triangular form after zero absorbing according to the badset
1226 //bset must be the zero set after absorbing
1227 static std::vector<std::vector<int> > ofindbases(int num, std::vector<int> bset,std::vector<std::vector<int> > gset)
1228 {
1229  int i,m;
1230  std::vector<std::vector<int> > bases;
1231  std::vector<int> fvars=freevars(num, bset, gset), base1;
1232  if (fvars.size()==0)
1233  {
1234  base1=make0(num);
1235  bases.push_back(base1);
1236  }
1237  else
1238  {
1239  for(i=0;i<fvars.size();i++)
1240  {
1241  m=fvars[i];
1242  base1=ofindbases1(num, m, bset, gset);
1243  bases.push_back(base1);
1244  }
1245  }
1246  //PrintS("They are the bases for the solution space:\n");
1247  //listsprint(bases);
1248  return bases;
1249 }
1250 
1251 //gset is a set of equations which have forms of x_i-x_j
1252 //num is the number of varialbes also the length of the set which we need to consider
1253 //output is trigular form of gset and badset where x_i=0
1254 static std::vector<std::vector<int> > eli2(int num, std::vector<int> bset,std::vector<std::vector<int> > gset)
1255 {
1256  std::vector<int> badset;
1257  std::vector<std::vector<int> > goodset, solve;
1258 //PrintS("This is the input bset\n");listprint(bset);
1259 //PrintS("This is the input gset\n");listsprint(gset);
1260  if(gset.size()!=0)//gset is not empty
1261  {
1262  //find all the variables which are zeroes
1263 
1264  if(bset.size()!=0)//bset is not empty
1265  {
1266  goodset=vAbsorb(bset, gset);//e.g. x_1=0, put x_i into the badset if x_i-x_1=0 or x_1-x_i=0
1267  int m=goodset.size();
1268  badset=goodset[m-1];
1269  goodset.erase(goodset.end());
1270  }
1271  else //bset is empty
1272  {
1273  goodset=gset;//badset is empty
1274  }//goodset is already the set which doesn't contain zero variables
1275 //PrintS("This is the badset after absorb \n");listprint(badset);
1276 //PrintS("This is the goodset after absorb \n");listsprint(goodset);
1277  goodset=soleli1(goodset);//get the triangular form of goodset
1278 //PrintS("This is the goodset after triangulization \n");listsprint(goodset);
1279  solve=ofindbases(num,badset,goodset);
1280  }
1281  else
1282  {
1283  solve=ofindbases(num,bset,gset);
1284  }
1285 //PrintS("This is the solution\n");listsprint(solve);
1286  return solve;
1287 }
1288 
1289 /********************************************************************/
1290 /************************links***********************************/
1291 
1292 //returns the links of face a in simplicial complex X
1293 static std::vector<std::vector<int> > links(poly a, ideal h)
1294 {
1295  int i;
1296  std::vector<std::vector<int> > lk,X=supports(h);
1297  std::vector<int> U,In,av=support1(a);
1298  for(i=0;i<X.size();i++)
1299  {
1300  U=vecUnion(av,X[i]);
1301  In=vecIntersection(av,X[i]);
1302  if( In.size()==0 && vInvsl(U,X))
1303  {
1304  //PrintS("The union of them is FACE and intersection is EMPTY!\n");
1305  lk.push_back(X[i]);
1306  }
1307  else
1308  {
1309  ;
1310  }
1311  }
1312  return lk;
1313 }
1314 
1315 static int redefinedeg(poly p, int num)
1316 {
1317  int deg=0, deg0;
1318  for(int i=1;i<=currRing->N;i++)
1319  {
1320  deg0=pGetExp(p, i);
1321  if(i>num)
1322  {
1323  deg= deg+2*deg0;
1324  }
1325  else
1326  {
1327  deg=deg+deg0;
1328  }
1329  }
1330  //Print("the new degree is: %d\n", deg);
1331  return (deg);
1332 }
1333 
1334 // the degree of variables should be same
1335 static ideal p_a(ideal h)
1336 {
1337  poly p;
1338  int i,j,deg=0,deg0;
1339  ideal aset=idCopy(h),ia,h1=idsrRing(h);
1340 //PrintS("idsrRing is:\n");id_print(h1);
1341  std::vector<int> as;
1342  std::vector<std::vector<int> > hvs=supports(h);
1343  for(i=0;i<IDELEMS(h1);i++)
1344  {
1345  deg0=pTotaldegree(h1->m[i]);
1346  if(deg < deg0)
1347  deg=deg0;
1348  }
1349  for(i=2;i<=deg;i++)
1350  {
1351  ia=id_MaxIdeal(i, currRing);
1352  for(j=0;j<IDELEMS(ia);j++)
1353  {
1354  p=pCopy(ia->m[j]);
1355  if(!IsInX(p,h))
1356  {
1357  as=support1(p);
1358  if(vInvsl(as,hvs))
1359  {
1360  idInsertPoly(aset, p);
1361  }
1362  }
1363  }
1364  }
1365  idSkipZeroes(aset);
1366  return(aset);
1367 }
1368 
1369 /*only for the exampels whose variables has degree more than 1*/
1370 /*ideal p_a(ideal h)
1371 {
1372  poly e=pOne(), p;
1373  int i,j,deg=0,deg0, ord=4;
1374  ideal aset=idCopy(h),ia,h1=idsrRing(h);
1375 //PrintS("idsrRing is:\n");id_print(h1);
1376  std::vector<int> as;
1377  std::vector<std::vector<int> > hvs=supports(h);
1378  for(i=0;i<IDELEMS(h1);i++)
1379  {
1380  deg0=redefinedeg(h1->m[i],ord);
1381  if(deg < deg0)
1382  deg=deg0;
1383  }
1384  for(i=2;i<=deg;i++)
1385  {
1386  ia=id_MaxIdeal(i, currRing);
1387  for(j=0;j<IDELEMS(ia);j++)
1388  {
1389  p=pCopy(ia->m[j]);
1390  if(!IsInX(p,h))
1391  {
1392  as=support1(p);
1393  if(vInvsl(as,hvs))
1394  {
1395  idInsertPoly(aset, p);
1396  }
1397  }
1398  }
1399  }
1400  idSkipZeroes(aset);
1401  return(aset);
1402 }*/
1403 
1404 static std::vector<int> vertset(std::vector<std::vector<int> > vecs)
1405 {
1406  int i,j;
1407  std::vector<int> vert;
1408  std::vector<std::vector<int> > vvs;
1409  for(i=1;i<=currRing->N;i++)
1410  {
1411  for(j=0;j<vecs.size();j++)
1412  {
1413  if(IsinL(i, vecs[j]))
1414  {
1415  if(!IsinL(i , vert))
1416  {
1417  vert.push_back(i);
1418  }
1419  break;
1420  }
1421  }
1422  }
1423  return (vert);
1424 }
1425 
1426 //smarter way
1427 static ideal p_b(ideal h, poly a)
1428 {
1429  std::vector<std::vector<int> > pbv,lk=links(a,h), res;
1430  std::vector<int> vert=vertset(lk), bv;
1431  res=b_subsets(vert);
1432  int i, adg=pTotaldegree(a);
1433  poly e=pOne();
1434  ideal idd=idInit(1,1);
1435  for(i=0;i<res.size();i++)
1436  {
1437  if(res[i].size()==adg)
1438  pbv.push_back(res[i]);
1439  }
1440  if(pEqualPolys(a,e))
1441  {
1442  idInsertPoly(idd, e);
1443  idSkipZeroes(idd);
1444  return (idd);
1445  }
1446  idd=idMaken(pbv);
1447  return(idd);
1448 }
1449 
1450 /*//dump way to get pb
1451 // the degree of variables should be same
1452 static ideal p_b(ideal h, poly a)
1453 {
1454  std::vector<std::vector<int> > pbv,lk=links(a,h),res;
1455 // PrintS("Its links are :\n");id_print(idMaken(lk));
1456  res=id_subsets(lk);
1457  //PrintS("res is :\n");listsprint(res);
1458  std::vector<int> bv;
1459  ideal bset=findb(h);
1460  int i,j,nu=res.size(),adg=pTotaldegree(a);
1461  poly e=pOne();ideal idd=idInit(1,1);
1462  for(i=0;i<res.size();i++)
1463  {
1464  if(res[i].size()==adg)
1465  pbv.push_back(res[i]);
1466  }
1467  if(pEqualPolys(a,e)){idInsertPoly(idd, e); idSkipZeroes(idd); return (idd);}
1468  for(i=0;i<nu;i++)
1469  {
1470  for(j=i+1;j<nu;j++)
1471  {
1472  if(res[i].size()!=0 && res[j].size()!=0)
1473  {
1474  bv = vecUnion(res[i], res[j]);
1475  if(IsInX(pMaken(bv),bset) && bv.size()==adg && !vInvsl(bv,pbv))
1476  {pbv.push_back(bv);}
1477  }
1478  }
1479  }
1480  idd=idMaken(pbv);
1481  //id_print(idd);
1482  return(idd);
1483 }*/
1484 
1485 // also only for the examples whose variables have degree more than 1(ndegreeb and p_b)
1486 /*int ndegreeb(std::vector<int> vec, int num)
1487 {
1488  int deg, deg0=0;
1489  for(int i=0;i<vec.size();i++)
1490  {
1491  if(vec[i]>num)
1492  {
1493  deg0++;
1494  }
1495  }
1496  deg=vec.size()+deg0;
1497  return(deg);
1498 }
1499 
1500 static ideal p_b(ideal h, poly a)
1501 {
1502  std::vector<std::vector<int> > pbv,lk=links(a,h),res;
1503 // PrintS("Its links are :\n");id_print(idMaken(lk));
1504  res=id_subsets(lk);
1505  //PrintS("res is :\n");listsprint(res);
1506  std::vector<int> bv;
1507  ideal bset=findb(h);
1508  int i,j,nu=res.size(),ord=4,adg=redefinedeg(a, ord);
1509  poly e=pOne();ideal idd=idInit(1,1);
1510  for(i=0;i<res.size();i++)
1511  {
1512  if(ndegreeb(res[i],ord)==adg)
1513  pbv.push_back(res[i]);
1514  }
1515  if(pEqualPolys(a,e)){idInsertPoly(idd, e); idSkipZeroes(idd); return (idd);}
1516  for(i=0;i<nu;i++)
1517  {
1518  for(j=i+1;j<nu;j++)
1519  {
1520  if(res[i].size()!=0 && res[j].size()!=0)
1521  {
1522  bv = vecUnion(res[i], res[j]);
1523  //PrintS("bv is :\n");listprint(bv);
1524  //Print("bv's degree is : %d\n", ndegreeb(bv,ord));
1525  if(IsInX(pMaken(bv),bset) && ndegreeb(bv,ord)==adg && !vInvsl(bv,pbv))
1526  {
1527  pbv.push_back(bv);
1528  }
1529  }
1530  }
1531  }
1532  idd=idMaken(pbv);
1533  //id_print(idd);
1534  return(idd);
1535 }*/
1536 
1537 //input is a squarefree monomial p
1538 //output is all the squarefree monomials which could divid p(including p itself?)
1539 static ideal psubset(poly p)
1540 {
1541  int i,j,max=pTotaldegree(p);
1542  ideal h1,mons, id_re=idInit(1,1);
1543  for(i=1;i<max;i++)
1544  {
1545  mons=id_MaxIdeal(i, currRing);
1546  h1=sfreemon(mons,i);
1547  for(j=0;j<IDELEMS(h1);j++)
1548  {
1549  if(p_DivisibleBy(h1->m[j],p,currRing))
1550  idInsertPoly(id_re, h1->m[j]);
1551  }
1552  }
1553  idSkipZeroes(id_re);
1554  //PrintS("This is the facset\n");
1555  //id_print(id_re);
1556  return id_re;
1557 }
1558 
1559 //inserts a new vector which has two elements a and b into vector gset (which is a vector of vectors)
1560 //(especially for gradedpiece1 and gradedpiece1n)
1561 static std::vector<std::vector<int> > listsinsertlist(std::vector<std::vector<int> > gset, int a, int b)
1562 {
1563  std::vector<int> eq;
1564  eq.push_back(a);
1565  eq.push_back(b);
1566  gset.push_back(eq);
1567  return gset;
1568 }
1569 
1570 static std::vector<int> makeequation(int i,int j, int t)
1571 {
1572  std::vector<int> equation;
1573  equation.push_back(i);
1574  equation.push_back(j);
1575  equation.push_back(t);
1576  //listprint(equation);
1577  return equation;
1578 }
1579 
1580 /****************************************************************/
1581 //only for solving the equations obtained from T^2
1582 //input should be a vector which has only 3 entries
1583 static poly pMake3(std::vector<int> vbase)
1584 {
1585  int co=1;
1586  poly p,q=0;
1587  for(int i=0;i<3;i++)
1588  {
1589  if(vbase[i]!=0)
1590  {
1591  if(i==1) co=-1;
1592  p = pOne();pSetExp(p, vbase[i], 1);pSetm(p);pSetCoeff(p, nInit(co));
1593  }
1594  else p=0;
1595  q = pAdd(q, p);
1596  co=1;
1597  }
1598  return q;
1599 }
1600 
1601 static ideal idMake3(std::vector<std::vector<int> > vecs)
1602 {
1603  ideal id_re=idInit(1,1);
1604  poly p;
1605  int i,lv=vecs.size();
1606  for(i=0;i<lv;i++)
1607  {
1608  p=pMake3(vecs[i]);
1609  idInsertPoly(id_re, p);
1610  }
1611  idSkipZeroes(id_re);
1612  return id_re;
1613 }
1614 
1615 /****************************************************************/
1616 
1617 //change the current ring to a new ring which is in num new variables
1618 static void equmab(int num)
1619 {
1620  int i;
1621  //Print("There are %d new variables for equations solving.\n",num);
1622  ring r=currRing;
1623  char** tt;
1624  coeffs cf=nCopyCoeff(r->cf);
1625  tt=(char**)omAlloc(num*sizeof(char *));
1626  for(i=0; i <num; i++)
1627  {
1628  tt[i] = (char*)omalloc(10); //if required enlarge it later
1629  sprintf (tt[i], "t(%d)", i+1);
1630  tt[i]=omStrDup(tt[i]);
1631  }
1632  ring R=rDefault(cf,num,tt,ringorder_lp);
1634  IDRING(h)=rCopy(R);
1635  rSetHdl(h);
1636 }
1637 
1638 //returns the trivial case of T^1
1639 //b must only contain one variable
1640 static std::vector<int> subspace1(std::vector<std::vector<int> > mv, std::vector<int> bv)
1641 {
1642  int i, num=mv.size();
1643  std::vector<int> base;
1644  for(i=0;i<num;i++)
1645  {
1646  if(IsinL(bv[0],mv[i]))
1647  base.push_back(1);
1648  else
1649  base.push_back(0);
1650  }
1651  return base;
1652 }
1653 
1654 /***************************only for T^2*************************************/
1655 //vbase only has two elements which records the position of the monomials in mv
1656 
1657 static std::vector<poly> pMakei(std::vector<std::vector<int> > mv,std::vector<int> vbase)
1658 {
1659  poly p;
1660  std::vector<poly> h1;
1661  int n=vbase.size();
1662  for(int i=0;i<n;i++)
1663  {
1664  p=pMaken(mv[vbase[i]]);
1665  h1.push_back(p);
1666  }
1667  return h1;
1668 }
1669 
1670 // returns a ideal according to a set of supports
1671 static std::vector<std::vector<poly> > idMakei(std::vector<std::vector<int> > mv,std::vector<std::vector<int> > vecs)
1672 {
1673  int i,lv=vecs.size();
1674  std::vector<std::vector<poly> > re;
1675  std::vector<poly> h;
1676  for(i=0;i<lv;i++)
1677  {
1678  h=pMakei(mv,vecs[i]);
1679  re.push_back(h);
1680  }
1681  //PrintS("This is the metrix M:\n");
1682  //listsprint(vecs);
1683  //PrintS("the ideal according to metrix M is:\n");
1684  return re;
1685 }
1686 
1687 /****************************************************************/
1688 
1689 //return the graded pieces of cohomology T^1 according to a,b
1690 //original method (only for debugging)
1691 static void gradedpiece1(ideal h,poly a,poly b)
1692 {
1693  int i,j,m;
1694  ideal sub=psubset(b);
1695  std::vector<int> av=support1(a), bv=support1(b), bad, vv;
1696  std::vector<std::vector<int> > hvs=supports(h), sbv=supports(sub), mv=Mabv(h,a,b),good;
1697  m=mv.size();
1698  ring r=currRing;
1699  if( m > 0 )
1700  {
1701  for(i=0;i<m;i++)
1702  {
1703  if(!vsubset(bv,mv[i]))
1704  {
1705  bad.push_back(i+1);
1706  }
1707  }
1708  for(i=0;i<m;i++)
1709  {
1710  for(j=i+1;j<m;j++)
1711  {
1712  vv=vecUnion(mv[i],mv[j]);
1713  if(mabconditionv(hvs,vv,av,bv))
1714  {
1715  good=listsinsertlist(good,i+1,j+1);
1716  }
1717  else
1718  {
1719  //PrintS("They are not in Mabt!\n");
1720  ;
1721  }
1722  }
1723  }
1724  std::vector<std::vector<int> > solve=eli2(m,bad,good);
1725  if(bv.size()!=1)
1726  {
1727  //PrintS("This is the solution of coefficients:\n");
1728  listsprint(solve);
1729  }
1730  else
1731  {
1732  std::vector<int> su=subspace1(mv,bv);
1733  //PrintS("This is the solution of subspace:\n");
1734  //listprint(su);
1735  std::vector<std::vector<int> > suu;
1736  suu.push_back(su);
1737  equmab(solve[0].size());
1738  std::vector<std::vector<int> > solves=vecqring(solve,suu);
1739  //PrintS("This is the solution of coefficients:\n");
1740  listsprint(solves);
1741  rChangeCurrRing(r);
1742  }
1743  }
1744  else
1745  {
1746  PrintS("No element considered!\n");
1747  }
1748 }
1749 
1750 //Returns true if b can divide p*q
1751 static bool condition1for2(std::vector<int > pv,std::vector<int > qv,std::vector<int > bv)
1752 {
1753  std::vector<int > vec=vecUnion(pv,qv);
1754  if(vsubset(bv,vec))
1755  {
1756  //PrintS("condition1for2 yes\n");
1757  return true;
1758  }
1759  //PrintS("condition1for2 no\n");
1760  return false;
1761 }
1762 
1763 //Returns true if support(p) union support(q) union support(s) union support(a) minus support(b) is face
1764 static bool condition2for2(std::vector<std::vector<int> > hvs, std::vector<int> pv, std::vector<int> qv, std::vector<int> sv, std::vector<int> av, std::vector<int> bv)
1765 {
1766  std::vector<int> vec=vecUnion(pv,qv);
1767  vec=vecUnion(vec,sv);
1768  if(mabconditionv(hvs,vec,av,bv))
1769  {
1770  //PrintS("condition2for2 yes\n");
1771  return (true);
1772  }
1773  //PrintS("condition2for2 no\n");
1774  return (false);
1775 }
1776 
1777 static bool condition3for2(std::vector<std::vector<int> > hvs, std::vector<int> pv, std::vector<int> qv, std::vector<int> av, std::vector<int> bv)
1778 {
1779  std::vector<int> v1,v2,v3;
1780  v1=vecIntersection(pv,qv);//intersection
1781  v2=vecUnion(pv,qv);
1782  v2=vecUnion(v2,av);
1783  v2=vecMinus(v2,bv);
1784  v3=vecUnion(v1,v2);
1785  if(vInvsl(v3,hvs))
1786  {
1787  //PrintS("condition3for2 yes\n");
1788  return(true);
1789  }
1790  //PrintS("condition3for2 no\n");
1791  return(false);
1792 }
1793 
1794 /****************solve the equations got from T^2*********************/
1795 
1796 static ideal getpresolve(ideal h)
1797 {
1798  //ring r=currRing;
1799  //assume (LIB "presolve.lib");
1800  sleftv a;a.Init();
1801  a.rtyp=IDEAL_CMD;a.data=(void*)h;
1802  idhdl solve=ggetid("elimlinearpart");
1803  if(solve==NULL)
1804  {
1805  WerrorS("presolve.lib are not loaded!");
1806  return NULL;
1807  }
1808  BOOLEAN sl=iiMake_proc(solve,NULL,&a);
1809  //PrintS("no errors here\n");
1810  if(sl)
1811  {
1812  WerrorS("error in solve!");
1813  }
1814  lists L=(lists) iiRETURNEXPR.Data();
1815  ideal re=(ideal)L->m[4].CopyD();
1816  //iiRETURNEXPR.CleanUp();
1817  iiRETURNEXPR.Init();
1818  //PrintS("no errors here\n");
1819  //idSkipZeroes(re);
1820  //id_print(re);
1821  return re;
1822 }
1823 
1824 static std::vector<int> numfree(ideal h)
1825 {
1826  int i,j;
1827  std::vector<int> fvar;
1828  for(j=1;j<=currRing->N;j++)
1829  {
1830  for(i=0;i<IDELEMS(h);i++)
1831  {
1832  if(vInp(j,h->m[i]))
1833  {
1834  fvar.push_back(j);
1835  break;
1836  }
1837  }
1838  }
1839  //Print("There are %d free variables in total\n",num);
1840  return fvar;
1841 }
1842 
1843 static std::vector<std::vector<int> > canonicalbase(int n)
1844 {
1845  std::vector<std::vector<int> > vecs;
1846  std::vector<int> vec;
1847  int i,j;
1848  for(i=0;i<n;i++)
1849  {
1850  for(j=0;j<n;j++)
1851  {
1852  if(i==j)
1853  vec.push_back(1);
1854  else
1855  vec.push_back(0);
1856  }
1857  vecs.push_back(vec);
1858  vec.clear();
1859  }
1860  return vecs;
1861 }
1862 
1863 static std::vector<std::vector<int> > getvector(ideal h,int n)
1864 {
1865  std::vector<int> vec;
1866  std::vector<std::vector<int> > vecs;
1867  ideal h2=idCopy(h);
1868  if(!idIs0(h))
1869  {
1870  ideal h1=getpresolve(h2);
1871  poly q,e=pOne();
1872  int lg=IDELEMS(h1),n,i,j,t;
1873  std::vector<int> fvar=numfree(h1);
1874  n=fvar.size();
1875  if(n==0)
1876  {
1877  vec=make0(IDELEMS(h1));vecs.push_back(vec);//listsprint(vecs);
1878  }
1879  else
1880  {
1881  for(t=0;t<n;t++)
1882  {
1883  vec.clear();
1884  for(i=0;i<lg;i++)
1885  {
1886  q=pCopy(h1->m[i]);
1887  //pWrite(q);
1888  if(q==0)
1889  {
1890  vec.push_back(0);
1891  }
1892  else
1893  {
1894  q=p_Subst(q, fvar[t], e,currRing);
1895  //Print("the %dth variable was substituted by 1:\n",fvar[t]);
1896  //pWrite(q);
1897  for(j=0;j<n;j++)
1898  {
1899  //Print("the %dth variable was substituted by 0:\n",fvar[j]);
1900  q=p_Subst(q, fvar[j],0,currRing);
1901  //pWrite(q);
1902  }
1903  if(q==0)
1904  {
1905  vec.push_back(0);
1906  }
1907  else
1908  {
1909  vec.push_back(n_Int(pGetCoeff(q),currRing->cf));
1910  }
1911  }
1912  }
1913  //listprint(vec);
1914  vecs.push_back(vec);
1915  }
1916  }
1917  }
1918  else
1919  {vecs=canonicalbase(n);}
1920  //listsprint(vecs);
1921  return vecs;
1922 }
1923 
1924 /**************************************************************************/
1925 
1926 //subspace of T2(find all the possible values of alpha)
1927 static std::vector<int> findalpha(std::vector<std::vector<int> > mv, std::vector<int> bv)
1928 {
1929  std::vector<int> alset;
1930  for(unsigned i=0;i<mv.size();i++)
1931  {
1932  if(vsubset(bv,mv[i]))
1933  {
1934  alset.push_back(i);
1935  }
1936  }
1937  //Print("This is the alpha set, and the subspace is dim-%ld\n",alset.size());
1938  //listprint(alset);
1939  return alset;
1940 }
1941 
1942 static std::vector<int> subspacet1(int num, std::vector<std::vector<int> > ntvs)
1943 {
1944  int i, j, t, n=ntvs.size();
1945  std::vector<int> subase;
1946  for(t=0;t<n;t++)
1947  {
1948  i=ntvs[t][0];
1949  j=ntvs[t][1];
1950  if(i==(num))
1951  {
1952  subase.push_back(1);
1953  }
1954  else if(j==num)
1955  {
1956  subase.push_back(-1);
1957  }
1958  else
1959  {
1960  subase.push_back(0);
1961  }
1962  }
1963  //Print("This is the basis w.r.t. %dth polynomial in alpha set\n",num);
1964  //listprint(subase);
1965  return subase;
1966 }
1967 
1968 //subspace for T^2(mab method)
1969 static std::vector<std::vector<int> > subspacet(std::vector<std::vector<int> > mv, std::vector<int> bv,std::vector<std::vector<int> > ntvs)
1970 {
1971  std::vector<int> alset=findalpha(mv,bv), subase;
1972  std::vector<std::vector<int> > subases;
1973  for(unsigned i=0;i<alset.size();i++)
1974  {
1975  subase=subspacet1(alset[i],ntvs);
1976  subases.push_back(subase);
1977  }
1978  //PrintS("These are the bases for the subspace:\n");
1979  //listsprint(subases);
1980  return subases;
1981 }
1982 
1983 static std::vector<std::vector<int> > mabtv(std::vector<std::vector<int> > hvs, std::vector<std::vector<int> > Mv, std::vector<int> av, std::vector<int> bv)
1984 {
1985  std::vector<int> v1,var;
1986  std::vector<std::vector<int> > vars;
1987  for(unsigned i=0;i<Mv.size();i++)
1988  {
1989  for(unsigned j=i+1;j<Mv.size();j++)
1990  {
1991  var.clear();
1992  v1=vecUnion(Mv[i],Mv[j]);
1993  if(mabconditionv(hvs, v1, av, bv))
1994  {
1995  var.push_back(i);
1996  var.push_back(j);
1997  vars.push_back(var);
1998  }
1999  }
2000  }
2001  return vars;
2002 }
2003 
2004 //fix the problem of the number of the new variables
2005 //original method for T^2(only for debugging)
2006 static void gradedpiece2(ideal h,poly a,poly b)
2007 {
2008  int t0,t1,t2,i,j,t,m;
2009  ideal sub=psubset(b);
2010  ring r=rCopy(currRing);
2011  std::vector<std::vector<int> > hvs=supports(h), mv=Mabv(h,a,b), mts, vecs,vars;
2012  std::vector<int> av=support1(a), bv=support1(b), vec,var;
2013  mts=mabtv(hvs,mv,av,bv);
2014  PrintS("The homomorphism should map onto:\n");
2015  lpsprint(idMakei(mv,mts));
2016  m=mv.size();
2017  if(m > 0)
2018  {
2019  vars=mabtv(hvs,mv,av,bv);
2020  int vn=vars.size();
2021  for(t0=0;t0<vars.size();t0++)
2022  {
2023  i=vars[t0][0];
2024  j=vars[t0][1];
2025  if(!condition1for2(mv[i],mv[j],bv))//condition 1
2026  {
2027  //PrintS("And they satisfy the condition 1.\n");
2028  vec=makeequation(t0+1,0,0);
2029  //PrintS("So the equation:\n");
2030  //pWrite(p);
2031  //PrintS("holds.\n");
2032  vecs.push_back(vec);
2033  vec.clear();
2034  }
2035  if(condition3for2(hvs,mv[i],mv[j],av,bv))//condition 3
2036  {
2037  //PrintS("And they satisfy the condition 3.\n");
2038  vec=makeequation(t0+1,0,0);
2039  //PrintS("So the equation: \n");
2040  //pWrite(p);
2041  //PrintS("holds.\n");
2042  vecs.push_back(vec);
2043  vec.clear();
2044  }
2045  for(t1=t0+1;t1<vars.size();t1++)
2046  {
2047  for(t2=t1+1;t2<vars.size();t2++)
2048  {
2049  if(vars[t0][0]==vars[t1][0]&&vars[t1][1]==vars[t2][1]&&vars[t0][1]==vars[t2][0])
2050  {
2051  i=vars[t0][0];
2052  j=vars[t0][1];
2053  t=vars[t1][1];
2054  if(condition2for2(hvs,mv[i],mv[j],mv[t],av,bv))//condition 2
2055  {
2056  vec=makeequation(t0+1,t1+1,t2+1);
2057  vecs.push_back(vec);
2058  vec.clear();
2059  }
2060  }
2061  }
2062  }
2063  }
2064  //PrintS("this is EQUATIONS:\n");
2065  //listsprint(vecs);
2066  equmab(vn);
2067  ideal id_re=idMake3(vecs);
2068  //id_print(id_re);
2069  std::vector<std::vector<int> > re=getvector(id_re,vn);
2070  PrintS("this is the solution for ideal :\n");
2071  listsprint(re);
2072  rChangeCurrRing(r);
2073  std::vector<std::vector<int> > sub=subspacet(mv, bv,vars);
2074  PrintS("this is the solution for subspace:\n");
2075  listsprint(sub);
2076  equmab(vn);
2077  std::vector<std::vector<int> > solve=vecqring(re, sub);
2078  PrintS("This is the solution of coefficients:\n");
2079  listsprint(solve);
2080  rChangeCurrRing(r);
2081  }
2082  else
2083  {
2084  PrintS("No element considered!");
2085  }
2086 }
2087 
2088 /**********************************************************************/
2089 //For the method of N_{a-b}
2090 
2091 //returns true if pv(support of monomial) satisfies pv union av minus bv is in hvs
2092 static bool nabconditionv(std::vector<std::vector<int> > hvs, std::vector<int> pv, std::vector<int> av, std::vector<int> bv)
2093 {
2094  std::vector<int> vec1=vecIntersection(pv,bv), vec2=vecUnion(pv,bv);
2095  int s1=vec1.size();
2096  if(!vInvsl(vec2,hvs) && s1==0 && vsubset(av,pv))
2097  {
2098  //PrintS("nab condition satisfied\n");
2099  return(true);
2100  }
2101  //PrintS("nab condition not satisfied\n");
2102  return(false);
2103 }
2104 
2105 //returns N_{a-b}
2106 static std::vector<std::vector<int> > Nabv(std::vector<std::vector<int> > hvs, std::vector<int> av, std::vector<int> bv)
2107 {
2108  std::vector<std::vector<int> > vecs;
2109  int num=hvs.size();
2110  for(int i=0;i<num;i++)
2111  {
2112  if(nabconditionv(hvs,hvs[i],av,bv))
2113  {
2114  //PrintS("satisfy:\n");
2115  vecs.push_back(hvs[i]);
2116  }
2117  }
2118  return vecs;
2119 }
2120 
2121 //returns true if pv union qv union av minus bv is in hvs
2122 //hvs is simplicial complex
2123 static bool nabtconditionv(std::vector<std::vector<int> > hvs, std::vector<int> pv, std::vector<int> qv)
2124 {
2125  std::vector<int> v1;
2126  v1=vecUnion(pv,qv);
2127  if(vInvsl(v1,hvs))
2128  {
2129  return (true);
2130  }
2131  return (false);
2132 }
2133 
2134 //returns N_{a-b}^(2)
2135 static std::vector<std::vector<int> > nabtv(std::vector<std::vector<int> > hvs, std::vector<std::vector<int> > Nv, std::vector<int> av, std::vector<int> bv)
2136 {
2137  std::vector<int> v1,var;
2138  std::vector<std::vector<int> > vars;
2139  for(unsigned i=0;i<Nv.size();i++)
2140  {
2141  for(unsigned j=i+1;j<Nv.size();j++)
2142  {
2143  var.clear();
2144  if(nabtconditionv(hvs, Nv[i], Nv[j]))
2145  {
2146  var.push_back(i);
2147  var.push_back(j);
2148  vars.push_back(var);
2149  }
2150  }
2151  }
2152  return vars;
2153 }
2154 
2155 //p must be the monomial which is a face
2156 // ideal sub=psubset(b); bvs=supports(sub);
2157 static bool tNab(std::vector<std::vector<int> > hvs, std::vector<int> pv, std::vector<std::vector<int> > bvs)
2158 {
2159  std::vector<int> sv;
2160  if(bvs.size()<=1) return false;
2161  for(unsigned i=0;i<bvs.size();i++)
2162  {
2163  sv=vecUnion(pv,bvs[i]);
2164  if(!vInvsl(sv,hvs))
2165  {
2166  return true;
2167  }
2168  }
2169  return false;
2170 }
2171 
2172 static std::vector<int> tnab(std::vector<std::vector<int> > hvs,std::vector<std::vector<int> > nvs,std::vector<std::vector<int> > bvs)
2173 {
2174  std::vector<int> pv, vec;
2175  for(unsigned j=0;j<nvs.size();j++)
2176  {
2177  pv=nvs[j];
2178  if(tNab(hvs, pv, bvs))
2179  {
2180  vec.push_back(j);
2181  }
2182  }
2183  return vec;
2184 }
2185 
2186 //the image phi(pv)=pv union av minus bv
2187 static std::vector<int> phimage(std::vector<int> pv, std::vector<int> av, std::vector<int> bv)
2188 {
2189  std::vector<int> qv=vecUnion(pv,av);
2190  qv=vecMinus(qv,bv);
2191  return qv;
2192 }
2193 
2194 //mvs and nvs are the supports of ideal Mab and Nab
2195 //vecs is the solution of nab
2196 static std::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)
2197 {
2198  int j;
2199  std::vector<int> pv, base;
2200  std::vector<std::vector<int> > bases;
2201  for(unsigned t=0;t<vecs.size();t++)
2202  {
2203  for(unsigned i=0;i<mvs.size();i++)
2204  {
2205  pv=phimage(mvs[i],av,bv);
2206  for( j=0;j<nvs.size();j++)
2207  {
2208  if(vEvl(pv,nvs[j]))
2209  {
2210  base.push_back(vecs[t][j]);
2211  break;
2212  }
2213  }
2214  if(j==nvs.size())
2215  {
2216  base.push_back(0);
2217  }
2218  }
2219  if(base.size()!=mvs.size())
2220  {
2221  //WerrorS("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1");
2222  WerrorS("Errors in Equations solving (Values Finding)!");
2223  usleep(1000000);
2224  assert(false);
2225 
2226  }
2227  bases.push_back(base);
2228  base.clear();
2229  }
2230  return bases;
2231 }
2232 
2233 static intvec *Tmat(std::vector<std::vector<int> > vecs)
2234 {
2235  //std::vector<std::vector<int> > solve=gradedpiece1n(h,a,b);
2236  //Print("the size of solve is: %ld\n",solve.size());
2237  //vtm(solve);
2238  intvec *m;
2239  int i,j, a=vecs.size();
2240  if(a==0)
2241  {
2242  m=new intvec(1,1,10);
2243  }
2244  else
2245  {
2246  int b=vecs[0].size();
2247  m=new intvec(a,b,0);
2248  for(i=1;i<=a;i++)
2249  {
2250  for(j=1;j<=b;j++)
2251  {
2252  IMATELEM(*m,i,j)=vecs[i-1][j-1];
2253  }
2254  }
2255  }
2256  return (m);
2257 }
2258 
2259 //returns the set of position number of minimal gens in M
2260 static std::vector<int> gensindex(ideal M, ideal ids)
2261 {
2262  int i;
2263  std::vector<int> vec,index;
2264  if(!idIs0(M))
2265  {
2266  std::vector<std::vector<int> > vecs=supports(ids);
2267  for(i=0;i<IDELEMS(M);i++)
2268  {
2269  vec=support1(M->m[i]);
2270  if(vInvsl(vec,vecs))
2271  index.push_back(i);
2272  }
2273  }
2274  return (index);
2275 }
2276 
2277 static ideal mingens(ideal h, poly a, poly b)
2278 {
2279  int i;
2280  std::vector<std::vector<int> > mv=Mabv(h,a,b);
2281  ideal M=idMaken(mv), hi=idInit(1,1);
2282  std::vector<int> index = gensindex(M, idsrRing(h));
2283  for(i=0;i<index.size();i++)
2284  {
2285  idInsertPoly(hi,M->m[index[i]]);
2286  }
2287  idSkipZeroes(hi);
2288  return (hi);
2289 }
2290 
2291 static std::vector<std::vector<int> > minisolve(std::vector<std::vector<int> > solve, std::vector<int> index)
2292 {
2293  int i,j;
2294  std::vector<int> vec,solm;
2295  std::vector<std::vector<int> > solsm;
2296  for(i=0;i<solve.size();i++)
2297  {
2298  vec=solve[i];
2299  for(j=0;j<vec.size();j++)
2300  {
2301  if(IsinL(j,index))
2302  solm.push_back(vec[j]);
2303  }
2304  solsm.push_back(solm);
2305  solm.clear();
2306  }
2307  return (solsm);
2308 }
2309 
2310 //T_1 graded piece(N method)
2311 //frame of the most efficient version
2312 //regardless of links
2313 static intvec * gradedpiece1n(ideal h,poly a,poly b)
2314 {
2315  int i,j,co,n;
2316  std::vector<std::vector<int> > hvs=supports(h),mv=Mabv(h,a,b),sbv,nv,good,solve;
2317  std::vector<int> av=support1(a), bv=support1(b), bad, tnv, index;
2318  ideal sub=psubset(b),M;
2319  sbv=supports(sub);
2320  nv=Nabv(hvs,av,bv);
2321  M=idMaken(mv);
2322  index = gensindex(M, idsrRing(h));
2323  n=nv.size();
2324  ring r=currRing;
2325  if(n > 0)
2326  {
2327  tnv=tnab(hvs,nv,sbv);
2328  for(i=0;i<tnv.size();i++)
2329  {
2330  co=tnv[i];
2331  bad.push_back(co+1);
2332  }
2333  for(i=0;i<n;i++)
2334  {
2335  for(j=i+1;j<n;j++)
2336  {
2337  if(nabtconditionv(hvs,nv[i],nv[j]))
2338  {
2339  good=listsinsertlist(good,i+1,j+1);
2340  }
2341  else
2342  {
2343  ;
2344  }
2345  }
2346  }
2347  solve=eli2(n,bad,good);
2348  if(bv.size()!=1)
2349  {;
2350  //PrintS("This is the solution of coefficients:\n");
2351  //listsprint(solve);
2352  }
2353  else
2354  {
2355  std::vector<int> su=make1(n);
2356  std::vector<std::vector<int> > suu;
2357  suu.push_back(su);
2358  equmab(n);
2359  solve=vecqring(solve,suu);
2360  //PrintS("This is the solution of coefficients:\n");
2361  //listsprint(solve);
2362  rChangeCurrRing(r);
2363  }
2364  solve=value1(mv,nv,solve,av,bv);
2365  }
2366  else
2367  {
2368  //PrintS("No element considered here!\n");
2369  solve.clear();
2370  }
2371  //PrintS("This is the solution of final coefficients:\n");
2372  //listsprint(solve);
2374  intvec *sl=Tmat(solve);
2375  //sl->show(0,0);
2376  return sl;
2377 }
2378 
2379 //for debugging
2380 static void T1(ideal h)
2381 {
2382  ideal bi=findb(h),ai;
2383  int mm=0;
2384  id_print(bi);
2385  poly a,b;
2386  std::vector<std::vector<int> > solve;
2387  for(int i=0;i<IDELEMS(bi);i++)
2388  {
2389  //PrintS("This is aset according to:");
2390  b=pCopy(bi->m[i]);
2391  pWrite(b);
2392  ai=finda(h,b,0);
2393  if(!idIs0(ai))
2394  {
2395  id_print(ai);
2396  for(int j=0;j<IDELEMS(ai);j++)
2397  {
2398  //PrintS("This is a:");
2399  a=pCopy(ai->m[j]);
2400  //pWrite(a);
2401  intvec * solve=gradedpiece1n(h, a, b);
2402  if (IMATELEM(*solve,1,1)!=10)
2403  mm++;
2404  }
2405  }
2406  }
2407  Print("Finished %d!\n",mm);
2408 }
2409 
2410 static bool condition2for2nv(std::vector<std::vector<int> > hvs, std::vector<int> pv, std::vector<int> qv, std::vector<int> fv)
2411 {
2412  std::vector<int> vec=vecUnion(pv,qv);
2413  vec=vecUnion(vec,fv);
2414  if(vInvsl(vec,hvs))
2415  {
2416  //PrintS("condition2for2 yes\n");
2417  return (true);
2418  }
2419  //PrintS("condition2for2 no\n");
2420  return (false);
2421 }
2422 
2423 //for subspace of T2(find all the possible values of alpha)
2424 static std::vector<int> findalphan(std::vector<std::vector<int> > N, std::vector<int> tN)
2425 {
2426  int i;std::vector<int> alset,vec;
2427  for(i=0;i<N.size();i++)
2428  {
2429  // vec=N[i];
2430  if(!IsinL(i,tN))
2431  {
2432  alset.push_back(i);
2433  }
2434  }
2435  //listprint(alset);
2436  return alset;
2437 }
2438 
2439 //subspace of T^2 (nab method)
2440 static std::vector<std::vector<int> > subspacetn(std::vector<std::vector<int> > N, std::vector<int> tN, std::vector<std::vector<int> > ntvs)
2441 {
2442  int i;
2443  std::vector<int> alset=findalphan(N,tN), subase;
2444  std::vector<std::vector<int> > subases;
2445  for(i=0;i<alset.size();i++)
2446  {
2447  subase=subspacet1(alset[i],ntvs);
2448  subases.push_back(subase);
2449  }
2450  //PrintS("These are the bases for the subspace:\n");
2451  //listsprint(subases);
2452  return subases;
2453 }
2454 
2455 //mts Mabt
2456 //nts Nabt
2457 //mvs Mab
2458 //nvs Nab
2459 static std::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)
2460 {
2461  int row,col,j;
2462  std::vector<int> pv,qv, base;
2463  std::vector<std::vector<int> > bases;
2464  //PrintS("This is the nabt:\n");
2465  //listsprint(nts);
2466  //PrintS("nabt ends:\n");
2467  //PrintS("This is the mabt:\n");
2468  //listsprint(mts);
2469  //PrintS("mabt ends:\n");
2470  for(unsigned t=0;t<vecs.size();t++)
2471  {
2472  for(unsigned i=0;i<mts.size();i++)
2473  {
2474  row=mts[i][0];
2475  col=mts[i][1];
2476  pv=phimage(mvs[row],av,bv);
2477  qv=phimage(mvs[col],av,bv);
2478  if(vEvl(pv,qv))
2479  base.push_back(0);
2480  else
2481  {
2482  for(j=0;j<nts.size();j++)
2483  {
2484  row=nts[j][0];
2485  col=nts[j][1];
2486  if(vEvl(pv,nvs[row])&&vEvl(qv,nvs[col]))
2487  {
2488  base.push_back(vecs[t][j]);break;
2489  }
2490  else if(vEvl(pv,nvs[col])&&vEvl(qv,nvs[row]))
2491  {
2492  base.push_back(-vecs[t][j]);break;
2493  }
2494  }
2495  if(j==nts.size()) {base.push_back(0);}
2496  }
2497  }
2498  if(base.size()!=mts.size())
2499  {
2500  WerrorS("Errors in Values Finding(value2)!");
2501  //WerrorS("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1");
2502  usleep(1000000);
2503  assert(false);
2504  }
2505  bases.push_back(base);
2506  base.clear();
2507  }
2508  return bases;
2509 }
2510 
2511 static ideal genst(ideal h, poly a, poly b)
2512 {
2513  std::vector<std::vector<int> > hvs=supports(h),mv,mts;
2514  std::vector<int> av=support1(a), bv=support1(b);
2515  mv=Mabv(h,a,b);
2516  mts=mabtv(hvs,mv,av,bv);
2517  std::vector<std::vector<poly> > pvs=idMakei(mv,mts);
2518  ideal gens=idInit(1,1);
2519  for(unsigned i=0;i<pvs.size();i++)
2520  {
2521  idInsertPoly(gens,pvs[i][0]);
2522  idInsertPoly(gens,pvs[i][1]);
2523  }
2524  idSkipZeroes(gens);
2525  return (gens);
2526 }
2527 
2528 static intvec * gradedpiece2n(ideal h,poly a,poly b)
2529 {
2530  int i,j,t,n;
2531  std::vector<std::vector<int> > hvs=supports(h),nv,mv,mts,sbv,vecs,vars,ntvs,solve;
2532  std::vector<int> av=support1(a), bv=support1(b),tnv,vec,var;
2533  ideal sub=psubset(b);
2534  sbv=supports(sub);
2535  nv=Nabv(hvs,av,bv);
2536  n=nv.size();
2537  tnv=tnab(hvs,nv,sbv);
2538  ring r=currRing;
2539  mv=Mabv(h,a,b);
2540  mts=mabtv(hvs,mv,av,bv);
2541  //PrintS("The relations are:\n");
2542  //listsprint(mts);
2543  //PrintS("The homomorphism should map onto:\n");
2544  //lpsprint(idMakei(mv,mts));
2545  if(n>0)
2546  {
2547  ntvs=nabtv( hvs, nv, av, bv);
2548  //PrintS("The current homomorphism map onto###:\n");
2549  //lpsprint(idMakei(nv,ntvs));
2550  int l=ntvs.size();
2551  for(int t0=0;t0<l;t0++)
2552  {
2553  i=ntvs[t0][0];
2554  j=ntvs[t0][1];
2555  if(tNab(hvs,nv[i],sbv)&&tNab(hvs,nv[j],sbv))//condition 1
2556  {
2557  vec=makeequation(t0+1,0,0);
2558  vecs.push_back(vec);
2559  vec.clear();
2560  }
2561  for(int t1=t0+1;t1<ntvs.size();t1++)
2562  {
2563  for(int t2=t1+1;t2<ntvs.size();t2++)
2564  {
2565  if(ntvs[t0][0]==ntvs[t1][0]&&ntvs[t1][1]==ntvs[t2][1]&&ntvs[t0][1]==ntvs[t2][0])
2566  {
2567  i=ntvs[t0][0];
2568  j=ntvs[t0][1];
2569  t=ntvs[t1][1];
2570  if(condition2for2nv(hvs,nv[i],nv[j],nv[t]))
2571  {
2572  vec=makeequation(t0+1,t1+1,t2+1);
2573  vecs.push_back(vec);
2574  vec.clear();
2575  }
2576  }
2577  }
2578  }
2579  }
2580  //PrintS("this is EQUATIONS:\n");
2581  //listsprint(vecs);
2582  if(n==1) l=1;
2583  equmab(l);
2584  ideal id_re=idMake3(vecs);
2585  //id_print(id_re);
2586  std::vector<std::vector<int> > re=getvector(id_re,l);
2587  //PrintS("this is the solution for ideal :\n");
2588  //listsprint(re);
2589  rChangeCurrRing(r);
2590  std::vector<std::vector<int> > sub=subspacetn(nv, tnv,ntvs);
2591  //PrintS("this is the solution for subspace:\n");
2592  //listsprint(sub);
2593  equmab(l);
2594  solve=vecqring(re, sub);
2595  //PrintS("This is the solution of coefficients:\n");
2596  //listsprint(solve);
2597  rChangeCurrRing(r);
2598  solve=value2(mv,nv,mts,ntvs,solve,av,bv);
2599  }
2600  else
2601  solve.clear();
2602  intvec *sl=Tmat(solve);
2603  return sl;
2604 }
2605 
2606 //for debugging
2607 static void T2(ideal h)
2608 {
2609  ideal bi=findb(h),ai;
2610  id_print(bi);
2611  poly a,b;
2612  int mm=0,gp=0;
2613  std::vector<int> bv,av;
2614  std::vector<std::vector<int> > solve;
2615  for(int i=0;i<IDELEMS(bi);i++)
2616  {
2617  b=pCopy(bi->m[i]);
2618  //bv=support1(b);
2619  //PrintS("This is aset according to:");
2620  pWrite(b);
2621 //if(bv.size()==2)
2622  //{
2623  ai=finda(h,b,0);
2624  if(!idIs0(ai))
2625  {
2626  PrintS("This is a set according to current b:\n");
2627  id_print(ai);
2628  for(int j=0;j<IDELEMS(ai);j++)
2629  {
2630  PrintS("This is a:");
2631  a=pCopy(ai->m[j]);
2632  pWrite(a);
2633  PrintS("This is b:");
2634  pWrite(b);
2635  intvec *solve=gradedpiece2n(h, a, b);
2636  delete solve;
2637  gp++;
2638  }
2639  }
2640  mm=mm+1;
2641  }
2642  if(mm==IDELEMS(bi))
2643  PrintS("Finished!\n");
2644  Print("There are %d graded pieces in total.\n",gp);
2645 }
2646 
2647 /*****************************for links*******************************************/
2648 //the image phi(pv)=pv minus av minus bv
2649 static std::vector<int> phimagel(std::vector<int> fv, std::vector<int> av, std::vector<int> bv)
2650 {
2651  std::vector<int> nv;
2652  nv=vecMinus(fv,bv);
2653  nv=vecMinus(nv,av);
2654  return nv;
2655 }
2656 
2657 //mvs and nvs are the supports of ideal Mab and Nab
2658 //vecs is the solution of nab
2659 static std::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)
2660 {
2661  int j;
2662  std::vector<int> pv;
2663  std::vector<int> base;
2664  std::vector<std::vector<int> > bases;
2665  for(unsigned t=0;t<vecs.size();t++)
2666  {
2667  for(unsigned i=0;i<mvs.size();i++)
2668  {
2669  pv=phimagel(mvs[i], av, bv);
2670  for(j=0;j<lks.size();j++)
2671  {
2672  if(vEvl(pv,lks[j]))
2673  {
2674  base.push_back(vecs[t][j]);break;
2675  }
2676  }
2677  //if(j==lks.size()) {base.push_back(0);}
2678  }
2679  if(base.size()!=mvs.size())
2680  {
2681  WerrorS("Errors in Values Finding(value1l)!");
2682  usleep(1000000);
2683  assert(false);
2684  }
2685  bases.push_back(base);
2686  base.clear();
2687  }
2688  return bases;
2689 }
2690 
2691 /***************************************************/
2693 /**************************************************/
2694 
2695 static void TimeShow(clock_t t_construct, clock_t t_solve, clock_t t_value ,clock_t t_total)
2696 {
2697  Print("The time of value matching for first order deformation: %.2f sec ;\n", ((double) t_value)/CLOCKS_PER_SEC);
2698  Print("The total time of fpiece: %.2f sec ;\n", ((double) t_total)/CLOCKS_PER_SEC);
2699  Print("The time of equations construction for fpiece: %.2f sec ;\n", ((double) t_construct)/CLOCKS_PER_SEC);
2700  Print("The total time of equations solving for fpiece: %.2f sec ;\n", ((double) t_solve)/CLOCKS_PER_SEC);
2701  PrintS("__________________________________________________________\n");
2702 }
2703 
2704 static std::vector<std::vector<int> > gpl(ideal h,poly a,poly b)
2705 {
2706  int i,j,co;
2707  std::vector<std::vector<int> > hvs=supports(h),sbv,nv,mv,good,solve;
2708  std::vector<int> av=support1(a), bv=support1(b),index,bad,tnv;
2709  ideal sub=psubset(b);
2710  sbv=supports(sub);
2711  nv=Nabv(hvs,av,bv);
2712  mv=Mabv(h,a,b);
2713  ideal M=idMaken(mv);
2714  index = gensindex(M, idsrRing(h));
2715  int n=nv.size();
2716  ring r=currRing;
2717  t_begin=clock();
2718  if(n > 0)
2719  {
2720  tnv=tnab(hvs,nv,sbv);
2721  for(i=0;i<tnv.size();i++)
2722  {
2723  co=tnv[i];
2724  bad.push_back(co+1);
2725  }
2726  for(i=0;i<n;i++)
2727  {
2728  for(j=i+1;j<n;j++)
2729  {
2730  if(nabtconditionv(hvs,nv[i],nv[j]))
2731  {
2732  good=listsinsertlist(good,i+1,j+1);
2733  }
2734  else
2735  {
2736  ;
2737  }
2738  }
2739  }
2741  t_begin=clock();
2742  solve=eli2(n,bad,good);
2743  t_solve=t_solve+clock()-t_begin;
2744  if(bv.size()!=1)
2745  {;
2746  }
2747  else
2748  {
2749  std::vector<int> su=make1(n);
2750  std::vector<std::vector<int> > suu;
2751  suu.push_back(su);
2752  equmab(n);
2753  solve=vecqring(solve,suu);
2754  rChangeCurrRing(r);
2755  }
2756  }
2757  else
2758  {
2759  solve.clear();
2760  }
2761  //listsprint(solve);
2762  //sl->show(0,0);
2763  return solve;
2764 }
2765 
2766 //T^1
2767 //only need to consider the links of a, and reduce a to empty set
2768 static intvec * gradedpiece1nl(ideal h,poly a,poly b, int set)
2769 {
2770  t_start=clock();
2771  poly e=pOne();
2772  std::vector<int> av=support1(a),bv=support1(b),index, em;
2773  std::vector<std::vector<int> > solve, hvs=supports(h), lks=links(a,h), mv=Mabv(h,a,b), nvl;
2774  ideal id_links=idMaken(lks);
2775  ideal M=idMaken(mv);
2776  index = gensindex(M, idsrRing(h));
2777  solve=gpl(id_links,e,b);
2778  t_mark=clock();
2779  nvl=Nabv(lks,em,bv);
2780  solve=value1l(mv, nvl , solve, av, bv);
2781  if(set==1)
2782  {
2784  }
2785  intvec *sl=Tmat(solve);
2786  t_value=t_value+clock()-t_mark;
2787  t_total=t_total+clock()-t_start;
2788  return sl;
2789 }
2790 
2791 //for finding values of T^2
2792 static std::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)
2793 {
2794  std::vector<int> pv,qv,base;
2795  int row,col,j;
2796  std::vector<std::vector<int> > bases;
2797  if(vecs.size()==0)
2798  {
2799 
2800  }
2801  for(unsigned t=0;t<vecs.size();t++)
2802  {
2803  for(unsigned i=0;i<mts.size();i++)
2804  {
2805  row=mts[i][0];
2806  col=mts[i][1];
2807  pv=phimagel(mvs[row],av,bv);
2808  qv=phimagel(mvs[col],av,bv);
2809  if(vEvl(pv,qv))
2810  base.push_back(0);
2811  else
2812  {
2813  for(j=0;j<lkts.size();j++)
2814  {
2815  row=lkts[j][0];
2816  col=lkts[j][1];
2817  if(vEvl(pv,lks[row])&&vEvl(qv,lks[col]))
2818  {
2819  base.push_back(vecs[t][j]);break;
2820  }
2821  else if(vEvl(qv,lks[row])&&vEvl(pv,lks[col]))
2822  {
2823  base.push_back(-vecs[t][j]);break;
2824  }
2825  }
2826  //if(j==lkts.size())
2827  //{
2828  //base.push_back(0);
2829  //}
2830  }
2831  }
2832  if(base.size()!=mts.size())
2833  {
2834  WerrorS("Errors in Values Finding!");
2835  //WerrorS("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1");
2836  usleep(1000000);
2837  assert(false);
2838  }
2839  bases.push_back(base);
2840  base.clear();
2841  }
2842  return bases;
2843 }
2844 
2845 static std::vector<std::vector<int> > gpl2(ideal h,poly a,poly b)
2846 {
2847  int i,j,t,n;
2848  std::vector<std::vector<int> > hvs=supports(h),sbv,nv,mv,mts,vecs,vars,ntvs,solve;
2849  std::vector<int> av=support1(a), bv=support1(b),vec,var,tnv;
2850  ideal sub=psubset(b);
2851  sbv=supports(sub);
2852  nv=Nabv(hvs,av,bv);
2853  n=nv.size();
2854  tnv=tnab(hvs,nv,sbv);
2855  ring r=currRing;
2856  mv=Mabv(h,a,b);
2857  mts=mabtv(hvs,mv,av,bv);
2858  if(n>0)
2859  {
2860  ntvs=nabtv( hvs, nv, av, bv);
2861  int l=ntvs.size();
2862  if(l>0)
2863  {
2864  for(int t0=0;t0<l;t0++)
2865  {
2866  i=ntvs[t0][0];
2867  j=ntvs[t0][1];
2868  if(tNab(hvs,nv[i],sbv)&&tNab(hvs,nv[j],sbv))//condition 1
2869  {
2870  vec=makeequation(t0+1,0,0);
2871  vecs.push_back(vec);
2872  vec.clear();
2873  }
2874  for(int t1=t0+1;t1<ntvs.size();t1++)
2875  {
2876  for(int t2=t1+1;t2<ntvs.size();t2++)
2877  {
2878  if(ntvs[t0][0]==ntvs[t1][0]&&ntvs[t1][1]==ntvs[t2][1]&&ntvs[t0][1]==ntvs[t2][0])
2879  {
2880  i=ntvs[t0][0];
2881  j=ntvs[t0][1];
2882  t=ntvs[t1][1];
2883  if(condition2for2nv(hvs,nv[i],nv[j],nv[t]))
2884  {
2885  vec=makeequation(t0+1,t1+1,t2+1);
2886  vecs.push_back(vec);
2887  vec.clear();
2888  }
2889  }
2890  }
2891  }
2892  }
2893  if(n==1) {l=1;}
2894  equmab(l);
2895  ideal id_re=idMake3(vecs);
2896  std::vector<std::vector<int> > re=getvector(id_re,l);
2897  rChangeCurrRing(r);
2898  std::vector<std::vector<int> > sub=subspacetn(nv, tnv,ntvs);
2899  equmab(l);
2900  solve=vecqring(re, sub);
2901  rChangeCurrRing(r);
2902  }
2903  else
2904  {
2905  solve.clear();
2906  }
2907  }
2908  else
2909  solve.clear();
2910  return solve;
2911 }
2912 
2913 static intvec * gradedpiece2nl(ideal h,poly a,poly b)
2914 {
2915  poly e=pOne();
2916  std::vector<int> av=support1(a), bv=support1(b), em;
2917  std::vector<std::vector<int> > hvs=supports(h), mv=Mabv(h,a,b),mts,solve,lks,nvl,ntsl;
2918  mts=mabtv(hvs,mv,av,bv);
2919  lks=links(a,h);
2920  ideal id_links=idMaken(lks);
2921 //PrintS("This is the links of a:\n"); id_print(id_links);
2922  nvl=Nabv(lks,em,bv);
2923 //PrintS("This is the N set:\n"); id_print(idMaken(nvl));
2924  ntsl=nabtv(lks,nvl,em,bv);
2925 //PrintS("This is N^2:\n"); listsprint(ntsl);
2926  solve=gpl2(id_links,e,b);
2927 //PrintS("This is pre solution of N:\n"); listsprint(solve);
2928  if(solve.size() > 0)
2929  {
2930  solve=value2l(mv, nvl, mts, ntsl, solve, av, bv);
2931  }
2932 //PrintS("This is solution of N:\n"); listsprint(solve);
2933  intvec *sl=Tmat(solve);
2934  return sl;
2935 }
2936 
2937 //for debugging
2938 /*
2939 void Tlink(ideal h,poly a,poly b,int n)
2940 {
2941  std::vector<std::vector<int> > hvs=supports(h);
2942  std::vector<int> av=support1(a);
2943  std::vector<int> bv=support1(b);
2944  std::vector<std::vector<int> > vec=links(a, h);
2945  PrintS("This is the links of a:\n");
2946  listsprint(vec);
2947  ideal li=idMaken(vec);
2948  PrintS("This is the links of a(ideal version):\n");
2949  id_print(li);
2950  poly p=pOne();
2951  PrintS("1************************************************\n");
2952  PrintS("This is T_1 (m):\n");
2953  gradedpiece1(li,p,b);
2954  PrintS("2************************************************\n");
2955  PrintS("This is T_2 (m):\n");
2956  gradedpiece2(li,p,b);
2957  PrintS("3************************************************\n");
2958  PrintS("This is T_1 (n):\n");
2959  gradedpiece1n(li,p,b);
2960  PrintS("4************************************************\n");
2961  PrintS("This is T_2 (n):\n");
2962  gradedpiece2n(li,p,b);
2963 }
2964 */
2965 
2966 /******************************for triangulation***********************************/
2967 
2968 //returns all the faces which are triangles
2969 static ideal trisets(ideal h)
2970 {
2971  int i;
2972  ideal ids=idInit(1,1);
2973  std::vector<int> pv;
2974  for(i=0;i<IDELEMS(h);i++)
2975  {
2976  pv= support1(h->m[i]);
2977  if(pv.size()==3)
2978  idInsertPoly(ids, pCopy(h->m[i]));
2979  }
2980  idSkipZeroes(ids);
2981  return ids;
2982 }
2983 
2984 // case 1 new faces
2985 static std::vector<std::vector<int> > triface(poly p, int vert)
2986 {
2987  std::vector<int> vec, fv=support1(p);
2988  std::vector<std::vector<int> > fvs0, fvs;
2989  vec.push_back(vert);
2990  fvs.push_back(vec);
2991  fvs0=b_subsets(fv);
2992  fvs0=vsMinusv(fvs0,fv);
2993  for(unsigned i=0;i<fvs0.size();i++)
2994  {
2995  vec=fvs0[i];
2996  vec.push_back(vert);
2997  fvs.push_back(vec);
2998  }
2999  return (fvs);
3000 }
3001 
3002 // the size of p's support must be 3
3003 //returns the new complex which is a triangulation based on the face p
3004 static ideal triangulations1(ideal h, poly p, int vert)
3005 {
3006  std::vector<int> vec, pv=support1(p);
3007  std::vector<std::vector<int> > vecs=supports(h),vs,vs0;
3008  vs0=triface(p,vert);
3009  vecs=vsMinusv(vecs, pv);
3010  vecs=vsUnion(vecs,vs0);
3011  //PrintS("This is the new simplicial complex according to the face \n"); pWrite(p);
3012  //PrintS("is:\n");
3013  //listsprint(vecs);
3014  ideal re=idMaken(vecs);
3015  return re;
3016 }
3017 
3018 /*
3019 static ideal triangulations1(ideal h)
3020 {
3021  int i,vert=currRing->N+1;
3022  std::vector<int> vec;
3023  std::vector<std::vector<int> > vecs=supports(h),vs,vs0;
3024  for (i=0;i<vecs.size();i++)
3025  {
3026  if((vecs[i]).size()==3)
3027  {
3028  vs0=triface(vecs[i],vert);
3029  vs=vsMinusv(vecs,vecs[i]);
3030  vs=vsUnion(vs,vs0);
3031  PrintS("This is the new simplicial complex according to the face \n");listprint(vecs[i]);
3032  PrintS("is:\n");
3033  listsprint(vs);
3034  }
3035  //else if((vecs[i]).size()==4)
3036  //tetraface(vecs[i]);
3037  }
3038  //ideal hh=idMaken(vs);
3039  return h;
3040 }*/
3041 
3042 static std::vector<int> commonedge(poly p, poly q)
3043 {
3044  std::vector<int> ev, fv1= support1(p), fv2= support2(q);
3045  for(unsigned i=0;i<fv1.size();i++)
3046  {
3047  if(IsinL(fv1[i], fv2))
3048  ev.push_back(fv1[i]);
3049  }
3050  return ev;
3051 }
3052 
3053 static intvec *edgemat(poly p, poly q)
3054 {
3055  intvec *m;
3056  int i;
3057  std::vector<int> dg=commonedge(p, q);
3058  int lg=dg.size();
3059  m=new intvec(lg);
3060  if(lg!=0)
3061  {
3062  m=new intvec(lg);
3063  for(i=0;i<lg;i++)
3064  {
3065  (*m)[i]=dg[i];
3066  }
3067  }
3068  return (m);
3069 }
3070 
3071 // case 2 the new face
3072 static std::vector<std::vector<int> > tetraface(poly p, poly q, int vert)
3073 {
3074  std::vector<int> ev=commonedge(p, q), vec, fv1=support1(p), fv2=support1(q);
3075  std::vector<std::vector<int> > fvs1, fvs2, fvs;
3076  vec.push_back(vert);
3077  fvs.push_back(vec);
3078  fvs1=b_subsets(fv1);
3079  fvs2=b_subsets(fv2);
3080  fvs1=vsMinusv(fvs1, fv1);
3081  fvs2=vsMinusv(fvs2, fv2);
3082  fvs2=vsUnion(fvs1, fvs2);
3083  fvs2=vsMinusv(fvs2, ev);
3084  for(unsigned i=0;i<fvs2.size();i++)
3085  {
3086  vec=fvs2[i];
3087  vec.push_back(vert);
3088  fvs.push_back(vec);
3089  }
3090  return (fvs);
3091 }
3092 
3093 //if p and q have a common edge
3094 static ideal triangulations2(ideal h, poly p, poly q, int vert)
3095 {
3096  std::vector<int> ev, fv1=support1(p), fv2=support1(q);
3097  std::vector<std::vector<int> > vecs=supports(h), vs1;
3098  ev=commonedge(p, q);
3099  vecs=vsMinusv(vecs, ev);
3100  vecs=vsMinusv(vecs,fv1);
3101  vecs=vsMinusv(vecs,fv2);
3102  vs1=tetraface(p, q, vert);
3103  vecs=vsUnion(vecs,vs1);
3104  ideal hh=idMaken(vecs);
3105  return hh;
3106 }
3107 
3108 // case 2 the new face
3109 static std::vector<std::vector<int> > penface(poly p, poly q, poly g, int vert)
3110 {
3111  int en=0;
3112  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);
3113  std::vector<std::vector<int> > fvs1, fvs2, fvs3, fvs, evec;
3114  evec.push_back(ev1);
3115  evec.push_back(ev2);
3116  evec.push_back(ev3);
3117  for(unsigned i=0;i<evec.size();i++)
3118  {
3119  if(evec[i].size()==2)
3120  {
3121  en++;
3122  }
3123  }
3124  if(en==2)
3125  {
3126  vec.push_back(vert);
3127  fvs.push_back(vec);
3128  fvs1=b_subsets(fv1);
3129  fvs2=b_subsets(fv2);
3130  fvs3=b_subsets(fv3);
3131  fvs1=vsMinusv(fvs1, fv1);
3132  fvs2=vsMinusv(fvs2, fv2);
3133  fvs3=vsMinusv(fvs3, fv3);
3134  fvs3=vsUnion(fvs3, fvs2);
3135  fvs3=vsUnion(fvs3, fvs1);
3136  for(unsigned i=0;i<evec.size();i++)
3137  {
3138  if(evec[i].size()==2)
3139  {
3140  fvs3=vsMinusv(fvs3, evec[i]);
3141  }
3142  }
3143  for(unsigned i=0;i<fvs3.size();i++)
3144  {
3145  vec=fvs3[i];
3146  vec.push_back(vert);
3147  fvs.push_back(vec);
3148  }
3149  }
3150  return (fvs);
3151 }
3152 
3153 static ideal triangulations3(ideal h, poly p, poly q, poly g, int vert)
3154 {
3155  std::vector<int> ev1=commonedge(p, q), ev2=commonedge(p, g), ev3=commonedge(q, g), fv1=support1(p), fv2=support1(q), fv3=support1(g);
3156  std::vector<std::vector<int> > vecs=supports(h), vs1, evec;
3157  evec.push_back(ev1);
3158  evec.push_back(ev2);
3159  evec.push_back(ev3);
3160  for(unsigned i=0;i<evec.size();i++)
3161  {
3162  if(evec[i].size()==2)
3163  {
3164  vecs=vsMinusv(vecs, evec[i]);
3165  }
3166  }
3167  vecs=vsMinusv(vecs,fv1);
3168  vecs=vsMinusv(vecs,fv2);
3169  vecs=vsMinusv(vecs,fv3);
3170  vs1=penface(p, q, g, vert);
3171  vecs=vsUnion(vecs,vs1);
3172  ideal hh=idMaken(vecs);
3173  return hh;
3174 }
3175 
3176 //returns p's valency in h
3177 //p must be a vertex
3178 static int valency(ideal h, poly p)
3179 {
3180  int val=0;
3181  std::vector<int> ev=support1(pCopy(p));
3182  int ver=ev[0];
3183 //PrintS("the vertex is :\n"); listprint(p);
3184  std::vector<std::vector<int> > vecs=supports(idCopy(h));
3185  for(unsigned i=0;i<vecs.size();i++)
3186  {
3187  if(vecs[i].size()==2 && IsinL(ver, vecs[i]))
3188  val++;
3189  }
3190  return (val);
3191 }
3192 
3193 /*ideal triangulations2(ideal h)
3194 {
3195  int i,j,vert=currRing->N+1;
3196  std::vector<int> ev;
3197  std::vector<std::vector<int> > vecs=supports(h),vs,vs0,vs1;
3198  vs0=tetrasets(h);
3199  for (i=0;i<vs0.size();i++)
3200  {
3201  for(j=i;j<vs0.size();j++)
3202  {
3203  ev=commonedge(vs0[i],vs0[j]);
3204  if(ev.size()==2)
3205  {
3206  vecs=vsMinusv(vecs, ev);
3207  vs=vsMinusv(vecs,vs0[i]);
3208  vs=vsMinusv(vecs,vs0[j]);
3209  vs1=tetraface(vs0[i],vs0[j],vert);
3210  vs=vsUnion(vs,vs1);
3211  PrintS("This is the new simplicial complex according to the face 1 \n");listprint(vecs[i]);
3212 PrintS("face 2: \n");
3213  PrintS("is:\n");
3214  listsprint(vs);
3215  }
3216  }
3217 
3218  //else if((vecs[i]).size()==4)
3219  //tetraface(vecs[i]);
3220  }
3221  //ideal hh=idMaken(vs);
3222  return h;
3223 }*/
3224 
3225 /*********************************For computation of X_n***********************************/
3226 static std::vector<std::vector<int> > vsMinusvs(std::vector<std::vector<int> > vs1, std::vector<std::vector<int> > vs2)
3227 {
3228  std::vector<std::vector<int> > vs=vs1;
3229  for(unsigned i=0;i<vs2.size();i++)
3230  {
3231  vs=vsMinusv(vs, vs2[i]);
3232  }
3233  return vs;
3234 }
3235 
3236 static std::vector<std::vector<int> > vs_subsets(std::vector<std::vector<int> > vs)
3237 {
3238  std::vector<std::vector<int> > sset, bv;
3239  for(unsigned i=0;i<vs.size();i++)
3240  {
3241  bv=b_subsets(vs[i]);
3242  sset=vsUnion(sset, bv);
3243  }
3244  return sset;
3245 }
3246 
3247 static std::vector<std::vector<int> > p_constant(ideal Xo, ideal Sigma)
3248 {
3249  std::vector<std::vector<int> > xs=supports(idCopy(Xo)), ss=supports(idCopy(Sigma)), fvs1;
3250  fvs1=vs_subsets(ss);
3251  fvs1=vsMinusvs(xs, fvs1);
3252  return fvs1;
3253 }
3254 
3255 static std::vector<std::vector<int> > p_change(ideal Sigma)
3256 {
3257  std::vector<std::vector<int> > ss=supports(idCopy(Sigma)), fvs;
3258  fvs=vs_subsets(ss);
3259  return (fvs);
3260 }
3261 
3262 static std::vector<std::vector<int> > p_new(ideal Xo, ideal Sigma)
3263 {
3264  int vert=0;
3265  std::vector<std::vector<int> > ss=supports(idCopy(Sigma)), fvs;
3266  for(int i=1;i<=currRing->N;i++)
3267  {
3268  for(int j=0;j<IDELEMS(Xo);j++)
3269  {
3270  if(pGetExp(Xo->m[j],i)>0)
3271  {
3272  vert=i+1;
3273  break;
3274  }
3275  }
3276  }
3277  int typ=ss.size();
3278  if(typ==1)
3279  {
3280  fvs=triface(Sigma->m[0], vert);
3281  }
3282  else if(typ==2)
3283  {
3284  fvs=tetraface(Sigma->m[0], Sigma->m[1], vert);
3285  }
3286  else
3287  {
3288  fvs=penface(Sigma->m[0], Sigma->m[1], Sigma->m[2], vert);
3289  }
3290  return (fvs);
3291 }
3292 
3293 static ideal c_New(ideal Io, ideal sig)
3294 {
3295  std::vector<std::vector<int> > vs1=p_constant(Io, sig), vs2=p_change(sig), vs3=p_new(Io, sig), vsig=supports(sig), vs;
3296  std::vector<int> ev;
3297  int ednum=vsig.size();
3298  if(ednum==2)
3299  {
3300  vsig.push_back(commonedge(sig->m[0], sig->m[1]));
3301  }
3302  else if(ednum==3)
3303  {
3304  for(int i=0;i<IDELEMS(sig);i++)
3305  {
3306  for(int j=i+1;j<IDELEMS(sig);j++)
3307  {
3308  ev=commonedge(sig->m[i], sig->m[j]);
3309  if(ev.size()==2)
3310  {
3311  vsig.push_back(ev);
3312  }
3313  }
3314  }
3315  }
3316 //PrintS("the first part is:\n");id_print(idMaken(vs1));
3317 //PrintS("the second part is:\n");id_print(idMaken(vsig));
3318 //PrintS("the third part is:\n");id_print(idMaken(vs3));
3319  vs2=vsMinusvs(vs2, vsig);
3320 //PrintS("the constant part2 is:\n");id_print(idMaken(vs2));
3321  vs=vsUnion(vs2, vs1);
3322 //PrintS("the constant part is:\n");id_print(idMaken(vs));
3323  vs=vsUnion(vs, vs3);
3324 //PrintS("the whole part is:\n");id_print(idMaken(vs));
3325  return(idMaken(vs));
3326 }
3327 
3328 static std::vector<std::vector<int> > phi1(poly a, ideal Sigma)
3329 {
3330  std::vector<std::vector<int> > ss=supports(idCopy(Sigma)), fvs;
3331  std::vector<int> av=support1(a), intvec, vv;
3332  for(unsigned i=0;i<ss.size();i++)
3333  {
3334  intvec=vecIntersection(ss[i], av);
3335  if(intvec.size()==av.size())
3336  {
3337  vv=vecMinus(ss[i], av);
3338  fvs.push_back(vv);
3339  }
3340  }
3341  return fvs;
3342 }
3343 
3344 static std::vector<std::vector<int> > phi2(poly a, ideal Xo, ideal Sigma)
3345 {
3346 
3347  std::vector<std::vector<int> > ss=p_new(Sigma, Xo), fvs;
3348  std::vector<int> av=support1(a), intvec, vv;
3349  for(unsigned i=0;i<ss.size();i++)
3350  {
3351  intvec=vecIntersection(ss[i], av);
3352  if(intvec.size()==av.size())
3353  {
3354  vv=vecMinus(ss[i], av);
3355  fvs.push_back(vv);
3356  }
3357  }
3358  return fvs;
3359 }
3360 
3361 static std::vector<std::vector<int> > links_new(poly a, ideal Xo, ideal Sigma, int vert, int ord)
3362 {
3363  std::vector<int> av=support1(a);
3364  std::vector<std::vector<int> > lko, lkn, lk1, lk2;
3365  lko=links(a, Xo);
3366  if(ord==1)
3367  return lko;
3368  if(ord==2)
3369  {
3370  lk1=phi1(a, Sigma);
3371  lk2=phi2(a, Xo, Sigma);
3372  lkn=vsMinusvs(lko, lk1);
3373  lkn=vsUnion(lkn, lk2);
3374  return lkn;
3375  }
3376  if(ord==3)
3377  {
3378  lkn=phi2(a, Xo, Sigma);
3379  return lkn;
3380  }
3381  WerrorS("Cannot find the links smartly!");
3382  return lko;
3383 }
3384 
3385 //returns 1 if there is a real divisor of b not in Xs
3386 static int existIn(poly b, ideal Xs)
3387 {
3388  std::vector<int> bv=support1(pCopy(b));
3389  std::vector<std::vector<int> > xvs=supports(idCopy(Xs)), bs=b_subsets(bv);
3390  bs=vsMinusv(bs, bv);
3391  for(unsigned i=0;i<bs.size();i++)
3392  {
3393  if(!vInvsl(bs[i], xvs))
3394  {
3395  return 1;
3396  }
3397  }
3398  return 0;
3399 }
3400 
3401 static int isoNum(poly p, ideal I, poly a, poly b)
3402 {
3403  int i;
3404  std::vector<std::vector<int> > vs=supports(idCopy(I));
3405  std::vector<int> v1=support1(a), v2=support1(b), v=support1(p);
3406  std::vector<int> vp, iv=phimagel(v, v1, v2);
3407  for(i=0;i<IDELEMS(I);i++)
3408  {
3409  vp=support1(pCopy(I->m[i]));
3410  if(vEvl(iv, phimagel(vp, v1, v2)))
3411  {
3412  return (i+1);
3413  }
3414  }
3415  return (0);
3416 }
3417 
3418 static int ifIso(poly p, poly q, poly f, poly g, poly a, poly b)
3419 {
3420  std::vector<int> va=support1(a), vb=support1(b), vp=support1(p), vq=support1(q), vf=support1(f), vg=support1(g);
3421  std::vector<int> v1=phimagel(vp, va, vb), v2=phimagel(vq, va, vb), v3=phimagel(vf, va, vb), v4=phimagel(vg, va, vb);
3422  if((vEvl(v1, v3)&& vEvl(v2,v4))||(vEvl(v1, v4)&& vEvl(v2,v3)) )
3423  {
3424  return (1);
3425  }
3426  return (0);
3427 }
3428 
3429 static ideal idMinusp(ideal I, poly p)
3430 {
3431  ideal h=idInit(1,1);
3432  int i;
3433  for(i=0;i<IDELEMS(I);i++)
3434  {
3435  if(!p_EqualPolys(I->m[i], p, currRing))
3436  {
3437  idInsertPoly(h, pCopy(I->m[i]));
3438  }
3439  }
3440  idSkipZeroes(h);
3441  return h;
3442 }
3443 
3444 /****************************for the interface of .lib*********************************/
3445 
3446 static std::vector<int> v_minus(std::vector<int> v1, std::vector<int> v2)
3447 {
3448  std::vector<int> vec;
3449  for(unsigned i=0;i<v1.size();i++)
3450  {
3451  vec.push_back(v1[i]-v2[i]);
3452  }
3453  return vec;
3454 }
3455 
3456 static std::vector<int> gdegree(poly a, poly b)
3457 {
3458  int i;
3459  std::vector<int> av,bv;
3460  for(i=1;i<=currRing->N;i++)
3461  {
3462  av.push_back(pGetExp(a,i));
3463  bv.push_back(pGetExp(b,i));
3464  }
3465  std::vector<int> vec=v_minus(av,bv);
3466  //PrintS("The degree is:\n");
3467  //listprint(vec);
3468  return vec;
3469 }
3470 
3471 /********************************for stellar subdivision******************************/
3472 
3473 static std::vector<std::vector<int> > star(poly a, ideal h)
3474 {
3475  int i;
3476  std::vector<std::vector<int> > st,X=supports(h);
3477  std::vector<int> U,av=support1(a);
3478  for(i=0;i<X.size();i++)
3479  {
3480  U=vecUnion(av,X[i]);
3481  if(vInvsl(U,X))
3482  {
3483  st.push_back(X[i]);
3484  }
3485  }
3486  return st;
3487 }
3488 
3489 static std::vector<std::vector<int> > boundary(poly a)
3490 {
3491  std::vector<int> av=support1(a), vec;
3492  std::vector<std::vector<int> > vecs;
3493  vecs=b_subsets(av);
3494  vecs.push_back(vec);
3495  vecs=vsMinusv(vecs, av);
3496  return vecs;
3497 }
3498 
3499 static std::vector<std::vector<int> > stellarsub(poly a, ideal h)
3500 {
3501  std::vector<std::vector<int> > vecs_minus, vecs_plus, lk=links(a,h), hvs=supports(h), sub, bys=boundary(a);
3502  std::vector<int> av=support1(a), vec, vec_n;
3503  int i,j,vert=0;
3504  for(i=1;i<=currRing->N;i++)
3505  {
3506  for(j=0;j<IDELEMS(h);j++)
3507  {
3508  if(pGetExp(h->m[j],i)>0)
3509  {
3510  vert=i+1;
3511  break;
3512  }
3513  }
3514  }
3515  vec_n.push_back(vert);
3516  for(i=0;i<lk.size();i++)
3517  {
3518  vec=vecUnion(av, lk[i]);
3519  vecs_minus.push_back(vec);
3520  for(j=0;j<bys.size();j++)
3521  {
3522  vec=vecUnion(lk[i], vec_n);
3523  vec=vecUnion(vec, bys[j]);
3524  vecs_plus.push_back(vec);
3525  }
3526  }
3527  sub=vsMinusvs(hvs, vecs_minus);
3528  sub=vsUnion(sub, vecs_plus);
3529  return(sub);
3530 }
3531 
3532 static std::vector<std::vector<int> > bsubsets_1(poly b)
3533 {
3534  std::vector<int> bvs=support1(b), vs;
3535  std::vector<std::vector<int> > bset;
3536  for(unsigned i=0;i<bvs.size();i++)
3537  {
3538  for(int j=0;j!=i; j++)
3539  {
3540  vs.push_back(bvs[j]);
3541  }
3542  bset.push_back(vs);
3543  vs.resize(0);
3544  }
3545  return bset;
3546 }
3547 
3548 /***************************for time testing******************************/
3549 static ideal T_1h(ideal h)
3550 {
3551  int i, j;
3552  //std::vector < intvec > T1;
3553  ideal ai=p_a(h), bi;
3554  //intvec *L;
3555  for(i=0;i<IDELEMS(ai);i++)
3556  {
3557  bi=p_b(h,ai->m[i]);
3558  if(!idIs0(bi))
3559  {
3560  for(j=0;j<IDELEMS(bi);j++)
3561  {
3562  //PrintS("This is for:\n");pWrite(ai->m[i]); pWrite(bi->m[j]);
3563  gradedpiece1nl(h,ai->m[i],bi->m[j], 0);
3564  //PrintS("Succeed!\n");
3565  //T1.push_back(L);
3566  }
3567  }
3568  }
3570  return h;
3571 }
3572 
3573 /**************************************interface T1****************************************/
3574 /*
3575 static BOOLEAN makeqring(leftv res, leftv args)
3576 {
3577  leftv h=args;
3578  ideal h2= id_complement( hh);
3579  if((h != NULL)&&(h->Typ() == POLY_CMD))
3580  {
3581  poly p= (poly)h->Data();
3582  h = h->next;
3583  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3584  {
3585  ideal hh=(ideal)h->Data();
3586  ideal h2=id_complement(hh);
3587  ideal h1=id_Init(1,1);
3588  idInsertPoly(h1,p);
3589  ideal gb=kStd(h2,NULL,testHomog,NULL,NULL,0,0,NULL);
3590  ideal idq=kNF(gb,NULL,h1);
3591  idSkipZeroes(h1);
3592  res->rtyp =POLY_CMD;
3593  res->data =h1->m[0];
3594  }
3595  }
3596  }
3597  return false;
3598 }*/
3599 
3601 {
3602  leftv h=args;
3603  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3604  {
3605  ideal hh=(ideal)h->Data();
3606  res->rtyp =IDEAL_CMD;
3607  res->data =idsrRing(hh);
3608  }
3609  return false;
3610 }
3611 
3613 {
3614  leftv h=args;
3615  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3616  {
3617  ideal hh=(ideal)h->Data();
3618  ideal h2= id_complement(hh);
3619  res->rtyp =IDEAL_CMD;
3620  res->data =h2;
3621  }
3622  return false;
3623 }
3624 
3625 static BOOLEAN t1h(leftv res, leftv args)
3626 {
3627  leftv h=args;
3628  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3629  {
3630  ideal hh=(ideal)h->Data();
3631  res->rtyp =IDEAL_CMD;
3632  res->data =T_1h(hh);
3633  }
3634  return false;
3635 }
3636 
3637 static BOOLEAN idsr(leftv res, leftv args)
3638 {
3639  leftv h=args;
3640  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3641  {
3642  ideal h1= (ideal)h->Data();
3643  h = h->next;
3644  if((h != NULL)&&(h->Typ() == POLY_CMD))
3645  {
3646  poly p= (poly)h->Data();
3647  h = h->next;
3648  if((h != NULL)&&(h->Typ() == POLY_CMD))
3649  {
3650  poly q= (poly)h->Data();
3651  res->rtyp =IDEAL_CMD;
3652  res->data =mingens(h1,p,q);
3653  }
3654  }
3655  }
3656  return false;
3657 }
3658 
3659 static intvec *dmat(poly a, poly b)
3660 {
3661  intvec *m;
3662  int i;
3663  std::vector<int> dg=gdegree(a,b);
3664  int lg=dg.size();
3665  m=new intvec(lg);
3666  if(lg!=0)
3667  {
3668  m=new intvec(lg);
3669  for(i=0;i<lg;i++)
3670  {
3671  (*m)[i]=dg[i];
3672  }
3673  }
3674  return (m);
3675 }
3676 
3677 static BOOLEAN gd(leftv res, leftv args)
3678 {
3679  leftv h=args;
3680  if((h != NULL)&&(h->Typ() == POLY_CMD))
3681  {
3682  poly p= (poly)h->Data();
3683  h = h->next;
3684  if((h != NULL)&&(h->Typ() == POLY_CMD))
3685  {
3686  poly q= (poly)h->Data();
3687  res->rtyp =INTVEC_CMD;
3688  res->data =dmat(p,q);
3689  }
3690  }
3691  return false;
3692 }
3693 
3695 {
3696  leftv h=args;
3697  if((h != NULL)&&(h->Typ() == POLY_CMD))
3698  {
3699  poly p= (poly)h->Data();
3700  h = h->next;
3701  if((h != NULL)&&(h->Typ() == POLY_CMD))
3702  {
3703  poly q= (poly)h->Data();
3704  res->rtyp =INTVEC_CMD;
3705  res->data =edgemat(p,q);
3706  }
3707  }
3708  return false;
3709 }
3710 
3711 static BOOLEAN fb(leftv res, leftv args)
3712 {
3713  leftv h=args;
3714  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3715  {
3716  ideal h1= (ideal)h->Data();
3717  res->rtyp =IDEAL_CMD;
3718  res->data =findb(h1);
3719  }
3720  return false;
3721 }
3722 
3723 static BOOLEAN pa(leftv res, leftv args)
3724 {
3725  leftv h=args;
3726  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3727  {
3728  ideal h1= (ideal)h->Data();
3729  res->rtyp =IDEAL_CMD;
3730  res->data =p_a(h1);
3731  }
3732  return false;
3733 }
3734 
3736 {
3737  leftv h=args;
3738  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3739  {
3740  ideal h1= (ideal)h->Data();
3741  res->rtyp =IDEAL_CMD;
3742  res->data =complementsimplex(h1);
3743  }
3744  return false;
3745 }
3746 
3747 static BOOLEAN pb(leftv res, leftv args)
3748 {
3749  leftv h=args;
3750  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3751  {
3752  ideal h1= (ideal)h->Data();
3753  h = h->next;
3754  if((h != NULL)&&(h->Typ() == POLY_CMD))
3755  {
3756  poly p= (poly)h->Data();
3757  res->rtyp =IDEAL_CMD;
3758  res->data =p_b(h1,p);
3759  }
3760  }
3761  return false;
3762 }
3763 
3764 static BOOLEAN fa(leftv res, leftv args)
3765 {
3766  leftv h=args;
3767  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3768  {
3769  ideal h1= (ideal)h->Data();
3770  h = h->next;
3771  if((h != NULL)&&(h->Typ() == POLY_CMD))
3772  {
3773  poly q= (poly)h->Data();
3774  h = h->next;
3775  if((h != NULL)&&(h->Typ() == INT_CMD))
3776  {
3777  int d= (int)(long)h->Data();
3778  res->rtyp =IDEAL_CMD;
3779  res->data =finda(h1,q,d);
3780  }
3781  }
3782  }
3783  return false;
3784 }
3785 
3786 static BOOLEAN fgp(leftv res, leftv args)
3787 {
3788  leftv h=args;
3789  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3790  {
3791  ideal h1= (ideal)h->Data();
3792  h = h->next;
3793  if((h != NULL)&&(h->Typ() == POLY_CMD))
3794  {
3795  poly p= (poly)h->Data();
3796  h = h->next;
3797  if((h != NULL)&&(h->Typ() == POLY_CMD))
3798  {
3799  poly q= (poly)h->Data();
3800  res->rtyp =INTVEC_CMD;
3801  res->data =gradedpiece1n(h1,p,q);
3802  }
3803  }
3804  }
3805  return false;
3806 }
3807 
3808 static BOOLEAN fgpl(leftv res, leftv args)
3809 {
3810  leftv h=args;
3811  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3812  {
3813  ideal h1= (ideal)h->Data();
3814  h = h->next;
3815  if((h != NULL)&&(h->Typ() == POLY_CMD))
3816  {
3817  poly p= (poly)h->Data();
3818  h = h->next;
3819  if((h != NULL)&&(h->Typ() == POLY_CMD))
3820  {
3821  poly q= (poly)h->Data();
3822  h = h->next;
3823  if((h != NULL)&&(h->Typ() == INT_CMD))
3824  {
3825  int d= (int)(long)h->Data();
3826  res->rtyp =INTVEC_CMD;
3827  res->data =gradedpiece1nl(h1,p,q,d);
3828  }
3829  }
3830  }
3831  }
3832  return false;
3833 }
3834 
3836 {
3837  leftv h=args;
3838  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3839  {
3840  ideal h1= (ideal)h->Data();
3841  h = h->next;
3842  if((h != NULL)&&(h->Typ() == POLY_CMD))
3843  {
3844  poly p= (poly)h->Data();
3845  h = h->next;
3846  if((h != NULL)&&(h->Typ() == POLY_CMD))
3847  {
3848  poly q= (poly)h->Data();
3849  res->rtyp =IDEAL_CMD;
3850  res->data =genst(h1,p,q);
3851  }
3852  }
3853  }
3854  return false;
3855 }
3856 
3857 static BOOLEAN sgp(leftv res, leftv args)
3858 {
3859  leftv h=args;
3860  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3861  {
3862  ideal h1= (ideal)h->Data();
3863  h = h->next;
3864  if((h != NULL)&&(h->Typ() == POLY_CMD))
3865  {
3866  poly p= (poly)h->Data();
3867  h = h->next;
3868  if((h != NULL)&&(h->Typ() == POLY_CMD))
3869  {
3870  poly q= (poly)h->Data();
3871  res->rtyp =INTVEC_CMD;
3872  res->data =gradedpiece2n(h1,p,q);
3873  }
3874  }
3875  }
3876  return false;
3877 }
3878 
3879 static BOOLEAN sgpl(leftv res, leftv args)
3880 {
3881  leftv h=args;
3882  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3883  {
3884  ideal h1= (ideal)h->Data();
3885  h = h->next;
3886  if((h != NULL)&&(h->Typ() == POLY_CMD))
3887  {
3888  poly p= (poly)h->Data();
3889  h = h->next;
3890  if((h != NULL)&&(h->Typ() == POLY_CMD))
3891  {
3892  poly q= (poly)h->Data();
3893  res->rtyp =INTVEC_CMD;
3894  res->data =gradedpiece2nl(h1,p,q);
3895  }
3896  }
3897  }
3898  return false;
3899 }
3900 
3901 static BOOLEAN Links(leftv res, leftv args)
3902 {
3903  leftv h=args;
3904  if((h != NULL)&&(h->Typ() == POLY_CMD))
3905  {
3906  poly p= (poly)h->Data();
3907  h = h->next;
3908  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3909  {
3910  ideal h1= (ideal)h->Data();
3911  res->rtyp =IDEAL_CMD;
3912  std::vector<std::vector<int> > vecs=links(p,h1);
3913  res->data =idMaken(vecs);
3914  }
3915  }
3916  return false;
3917 }
3918 
3919 static BOOLEAN isSim(leftv res, leftv args)
3920 {
3921  leftv h=args;
3922  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3923  {
3924  ideal h1= (ideal)h->Data();
3925  res->rtyp =IDEAL_CMD;
3926  res->data =IsSimplex(h1);
3927  }
3928  return false;
3929 }
3930 
3932 {
3933  leftv h=args;
3934  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3935  {
3936  ideal h1= (ideal)h->Data();
3937  h = h->next;
3938  if((h != NULL)&&(h->Typ() == POLY_CMD))
3939  {
3940  poly p= (poly)h->Data();
3941  h = h->next;
3942  if((h != NULL)&&(h->Typ() == INT_CMD))
3943  {
3944  int d= (int)(long)h->Data();
3945  res->rtyp =IDEAL_CMD;
3946  res->data =triangulations1(h1, p, d);
3947  }
3948  }
3949  }
3950  return false;
3951 }
3952 
3954 {
3955  leftv h=args;
3956  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3957  {
3958  ideal h1= (ideal)h->Data();
3959  h = h->next;
3960  if((h != NULL)&&(h->Typ() == POLY_CMD))
3961  {
3962  poly p= (poly)h->Data();
3963  h = h->next;
3964  if((h != NULL)&&(h->Typ() == POLY_CMD))
3965  {
3966  poly q= (poly)h->Data();
3967  h = h->next;
3968  if((h != NULL)&&(h->Typ() == INT_CMD))
3969  {
3970  int d= (int)(long)h->Data();
3971  res->rtyp =IDEAL_CMD;
3972  res->data =triangulations2(h1,p,q,d);
3973  }
3974  }
3975  }
3976  }
3977  return false;
3978 }
3979 
3981 {
3982  leftv h=args;
3983  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
3984  {
3985  ideal h1= (ideal)h->Data();
3986  h = h->next;
3987  if((h != NULL)&&(h->Typ() == POLY_CMD))
3988  {
3989  poly p= (poly)h->Data();
3990  h = h->next;
3991  if((h != NULL)&&(h->Typ() == POLY_CMD))
3992  {
3993  poly q= (poly)h->Data();
3994  h = h->next;
3995  if((h != NULL)&&(h->Typ() == POLY_CMD))
3996  {
3997  poly g= (poly)h->Data();
3998  h = h->next;
3999  if((h != NULL)&&(h->Typ() == INT_CMD))
4000  {
4001  int d= (int)(long)h->Data();
4002  res->rtyp =IDEAL_CMD;
4003  res->data =triangulations3(h1,p,q,g,d);
4004  }
4005  }
4006  }
4007  }
4008  }
4009  return false;
4010 }
4011 
4013 {
4014  leftv h=args;int i;
4015  std::vector<int> bset,bs;
4016  std::vector<std::vector<int> > gset;
4017  if((h != NULL)&&(h->Typ() == INT_CMD))
4018  {
4019  int n= (int)(long)h->Data();
4020  h = h->next;
4021  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4022  {
4023  ideal bi= (ideal)h->Data();
4024  h = h->next;
4025  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4026  {
4027  ideal gi= (ideal)h->Data();
4028  for(i=0;i<IDELEMS(bi);i++)
4029  {
4030  bs=support1(bi->m[i]);
4031  if(bs.size()==1)
4032  bset.push_back(bs[0]);
4033  else if(bs.size()==0)
4034  ;
4035  else
4036  {
4037  WerrorS("Errors in T^1 Equations Solving!");
4038  usleep(1000000);
4039  assert(false);
4040  }
4041 
4042  }
4043  gset=supports2(gi);
4044  res->rtyp =INTVEC_CMD;
4045  std::vector<std::vector<int> > vecs=eli2(n,bset,gset);
4046  res->data =Tmat(vecs);
4047  }
4048  }
4049  }
4050  return false;
4051 }
4052 
4053 static BOOLEAN tsets(leftv res, leftv args)
4054 {
4055  leftv h=args;
4056  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4057  {
4058  ideal h1= (ideal)h->Data();
4059  res->rtyp =IDEAL_CMD;
4060  res->data =trisets(h1);
4061  }
4062  return false;
4063 }
4064 
4066 {
4067  leftv h=args;
4068  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4069  {
4070  ideal h1= (ideal)h->Data();
4071  h = h->next;
4072  if((h != NULL)&&(h->Typ() == POLY_CMD))
4073  {
4074  poly p= (poly)h->Data();
4075  res->rtyp =INT_CMD;
4076  res->data =(void *)(long)valency(h1,p);
4077  }
4078  }
4079  return false;
4080 }
4081 
4082 static BOOLEAN nabvl(leftv res, leftv args)
4083 {
4084  leftv h=args;
4085  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4086  {
4087  ideal h1= (ideal)h->Data();
4088  h = h->next;
4089  if((h != NULL)&&(h->Typ() == POLY_CMD))
4090  {
4091  poly p= (poly)h->Data();
4092  h = h->next;
4093  if((h != NULL)&&(h->Typ() == POLY_CMD))
4094  {
4095  poly q= (poly)h->Data();
4096  res->rtyp =IDEAL_CMD;
4097  std::vector<std::vector<int> > vecs=supports(h1);
4098  std::vector<int> pv=support1(p), qv=support1(q);
4099  res->data =idMaken(Nabv(vecs,pv,qv));
4100  }
4101  }
4102  }
4103  return false;
4104 }
4105 
4107 {
4108  leftv h=args;
4109  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4110  {
4111  ideal h1= (ideal)h->Data();
4112  h = h->next;
4113  if((h != NULL)&&(h->Typ() == POLY_CMD))
4114  {
4115  poly p= (poly)h->Data();
4116  h = h->next;
4117  if((h != NULL)&&(h->Typ() == POLY_CMD))
4118  {
4119  poly q= (poly)h->Data();
4120  res->rtyp =IDEAL_CMD;
4121  std::vector<std::vector<int> > vecs=supports(h1), sbv,tnbr;
4122  std::vector<int> pv=support1(p), qv=support1(q);
4123  std::vector<std::vector<int> > nvs=Nabv(vecs, pv, qv);
4124  ideal sub=psubset(q);
4125  sbv=supports(sub);
4126  std::vector<int> tnv =tnab(vecs,nvs,sbv);
4127  for(unsigned i=0;i<tnv.size();i++)
4128  {
4129  tnbr.push_back(nvs[tnv[i]]);
4130  }
4131  res->data =idMaken(tnbr);
4132  }
4133  }
4134  }
4135  return false;
4136 }
4137 
4139 {
4140  leftv h=args;
4141  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4142  {
4143  ideal h1= (ideal)h->Data();
4144  h = h->next;
4145  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4146  {
4147  ideal h2= (ideal)h->Data();
4148  res->rtyp =INT_CMD;
4149  std::vector<std::vector<int> > vs1=supports(h1), vs2=supports(h2);
4150  res->data =(void *)(long)(vsIntersection(vs1, vs2).size());
4151  }
4152  }
4153  return false;
4154 }
4155 
4156 static BOOLEAN mabvl(leftv res, leftv args)
4157 {
4158  leftv h=args;
4159  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4160  {
4161  ideal h1= (ideal)h->Data();
4162  h = h->next;
4163  if((h != NULL)&&(h->Typ() == POLY_CMD))
4164  {
4165  poly p= (poly)h->Data();
4166  h = h->next;
4167  if((h != NULL)&&(h->Typ() == POLY_CMD))
4168  {
4169  poly q= (poly)h->Data();
4170  res->rtyp =IDEAL_CMD;
4171  res->data =idMaken(Mabv(h1,p,q));
4172  }
4173  }
4174  }
4175  return false;
4176 }
4177 
4179 {
4180  leftv h=args;
4181  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4182  {
4183  ideal h1= (ideal)h->Data();
4184  h = h->next;
4185  if((h != NULL)&&(h->Typ() == POLY_CMD))
4186  {
4187  poly p= (poly)h->Data();
4188  h = h->next;
4189  if((h != NULL)&&(h->Typ() == POLY_CMD))
4190  {
4191  poly q= (poly)h->Data();
4192  std::vector<std::vector<int> > hvs=supports(h1), nv, ntvs;
4193  std::vector<int> av=support1(p), bv=support1(q);
4194  nv=Nabv(hvs,av,bv);
4195  ntvs=nabtv( hvs, nv, av, bv);
4196  std::vector<std::vector<poly> > pvs=idMakei(nv,ntvs);
4197  ideal gens=idInit(1,1);
4198  for(unsigned i=0;i<pvs.size();i++)
4199  {
4200  idInsertPoly(gens,pvs[i][0]);
4201  idInsertPoly(gens,pvs[i][1]);
4202  }
4203  idSkipZeroes(gens);
4204  res->rtyp =IDEAL_CMD;
4205  res->data =gens;
4206  }
4207  }
4208  }
4209  return false;
4210 }
4211 
4212 static BOOLEAN linkn(leftv res, leftv args)
4213 {
4214  leftv h=args;
4215  if((h != NULL)&&(h->Typ() == POLY_CMD))
4216  {
4217  poly a= (poly)h->Data();
4218  h = h->next;
4219  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4220  {
4221  ideal Xo= (ideal)h->Data();
4222  h = h->next;
4223  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4224  {
4225  ideal Sigma= (ideal)h->Data();
4226  h = h->next;
4227  if((h != NULL)&&(h->Typ() == INT_CMD))
4228  {
4229  int vert= (int)(long)h->Data();
4230  h = h->next;
4231  if((h != NULL)&&(h->Typ() == INT_CMD))
4232  {
4233  int ord= (int)(long)h->Data();
4234  res->rtyp =IDEAL_CMD;
4235  res->data =idMaken(links_new(a, Xo, Sigma, vert, ord));
4236  }
4237  }
4238  }
4239  }
4240  }
4241  return false;
4242 }
4243 
4245 {
4246  leftv h=args;
4247  if((h != NULL)&&(h->Typ() == POLY_CMD))
4248  {
4249  poly p= (poly)h->Data();
4250  h = h->next;
4251  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4252  {
4253  ideal h1= (ideal)h->Data();
4254  res->rtyp =INT_CMD;
4255  res->data =(void *)(long)existIn(p, h1);
4256  }
4257  }
4258  return false;
4259 }
4260 
4262 {
4263  leftv h=args;
4264  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4265  {
4266  ideal h1= (ideal)h->Data();
4267  h = h->next;
4268  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4269  {
4270  ideal h2= (ideal)h->Data();
4271  res->rtyp =IDEAL_CMD;
4272  res->data =idMaken(p_constant(h1,h2));
4273  }
4274  }
4275  return false;
4276 }
4277 
4279 {
4280  leftv h=args;
4281  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4282  {
4283  ideal h1= (ideal)h->Data();
4284  res->rtyp =IDEAL_CMD;
4285  res->data =idMaken(p_change(h1));
4286  }
4287  return false;
4288 }
4289 
4290 static BOOLEAN p_New(leftv res, leftv args)
4291 {
4292  leftv h=args;
4293  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4294  {
4295  ideal h1= (ideal)h->Data();
4296  h = h->next;
4297  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4298  {
4299  ideal h2= (ideal)h->Data();
4300  res->rtyp =IDEAL_CMD;
4301  res->data =idMaken(p_new(h1,h2));
4302  }
4303  }
4304  return false;
4305 }
4306 
4308 {
4309  leftv h=args;
4310  if((h != NULL)&&(h->Typ() == POLY_CMD))
4311  {
4312  poly p= (poly)h->Data();
4313  res->rtyp =INT_CMD;
4314  res->data =(void *)(long)(support1(p).size());
4315  }
4316  return false;
4317 }
4318 
4320 {
4321  leftv h=args;
4322  if((h != NULL)&&(h->Typ() == POLY_CMD))
4323  {
4324  poly p= (poly)h->Data();
4325  res->rtyp =IDEAL_CMD;
4326  res->data =idMaken(bsubsets_1(p));
4327  }
4328  return false;
4329 }
4330 
4332 {
4333  leftv h=args;
4334  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4335  {
4336  ideal h1= (ideal)h->Data();
4337  h = h->next;
4338  if((h != NULL)&&(h->Typ() == POLY_CMD))
4339  {
4340  poly p= (poly)h->Data();
4341  res->rtyp =IDEAL_CMD;
4342  res->data =idMinusp(h1, p);
4343  }
4344  }
4345  return false;
4346 }
4347 
4349 {
4350  leftv h=args;
4351  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4352  {
4353  ideal h1= (ideal)h->Data();
4354  h = h->next;
4355  if((h != NULL)&&(h->Typ() == POLY_CMD))
4356  {
4357  poly p= (poly)h->Data();
4358  std::vector<std::vector<int> > st=star(p, h1);
4359  std::vector<std::vector<int> > hvs=supports(h1);
4360  std::vector<std::vector<int> > re= vsMinusvs(hvs, st);
4361  res->rtyp =IDEAL_CMD;
4362  res->data =idMaken(re);
4363  }
4364  }
4365  return false;
4366 }
4367 
4368 static BOOLEAN cNew(leftv res, leftv args)
4369 {
4370  leftv h=args;
4371  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4372  {
4373  ideal h1= (ideal)h->Data();
4374  h = h->next;
4375  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4376  {
4377  ideal h2= (ideal)h->Data();
4378  res->rtyp =IDEAL_CMD;
4379  res->data =c_New(h1, h2);
4380  }
4381  }
4382  return false;
4383 }
4384 
4385 static BOOLEAN stars(leftv res, leftv args)
4386 {
4387  leftv h=args;
4388  if((h != NULL)&&(h->Typ() == POLY_CMD))
4389  {
4390  poly p= (poly)h->Data();
4391  h = h->next;
4392  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4393  {
4394  ideal h1= (ideal)h->Data();
4395  res->rtyp =IDEAL_CMD;
4396  res->data =idMaken(star(p, h1));
4397  }
4398  }
4399  return false;
4400 }
4401 
4403 {
4404  leftv h=args;
4405  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4406  {
4407  ideal h2= (ideal)h->Data();
4408  h = h->next;
4409  if((h != NULL)&&(h->Typ() == POLY_CMD))
4410  {
4411  poly p= (poly)h->Data();
4412  res->rtyp =IDEAL_CMD;
4413  res->data =idMaken(stellarsub(p, h2));
4414  }
4415  }
4416  return false;
4417 }
4418 
4420 {
4421  leftv h=args;
4422  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4423  {
4424  ideal h1= (ideal)h->Data();
4425  h = h->next;
4426  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4427  {
4428  ideal h2= (ideal)h->Data();
4429  res->rtyp =IDEAL_CMD;
4430  res->data =idmodulo(h1, h2);
4431  }
4432  }
4433  return false;
4434 }
4435 
4437 {
4438  leftv h=args;
4439  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4440  {
4441  ideal h1= (ideal)h->Data();
4442  h = h->next;
4443  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4444  {
4445  ideal h2= (ideal)h->Data();
4446  res->rtyp =IDEAL_CMD;
4447  res->data =idMinus(h1, h2);
4448  }
4449  }
4450  return false;
4451 }
4452 
4454 {
4455  leftv h=args;
4456  if((h != NULL)&&(h->Typ() == POLY_CMD))
4457  {
4458  poly p= (poly)h->Data();
4459  h = h->next;
4460  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4461  {
4462  ideal h1= (ideal)h->Data();
4463  h = h->next;
4464  if((h != NULL)&&(h->Typ() == POLY_CMD))
4465  {
4466  poly a= (poly)h->Data();
4467  h = h->next;
4468  if((h != NULL)&&(h->Typ() == POLY_CMD))
4469  {
4470  poly b= (poly)h->Data();
4471  res->rtyp =INT_CMD;
4472  res->data =(void *)(long)isoNum(p, h1, a, b);
4473  }
4474  }
4475  }
4476  }
4477  return false;
4478 }
4479 
4481 {
4482  leftv h=args;
4483  if((h != NULL)&&(h->Typ() == POLY_CMD))
4484  {
4485  poly p= (poly)h->Data();
4486  h = h->next;
4487  if((h != NULL)&&(h->Typ() == POLY_CMD))
4488  {
4489  poly q= (poly)h->Data();
4490  h = h->next;
4491  if((h != NULL)&&(h->Typ() == POLY_CMD))
4492  {
4493  poly f= (poly)h->Data();
4494  h = h->next;
4495  if((h != NULL)&&(h->Typ() == POLY_CMD))
4496  {
4497  poly g= (poly)h->Data();
4498  h = h->next;
4499  if((h != NULL)&&(h->Typ() == POLY_CMD))
4500  {
4501  poly a= (poly)h->Data();
4502  h = h->next;
4503  if((h != NULL)&&(h->Typ() == POLY_CMD))
4504  {
4505  poly b= (poly)h->Data();
4506  res->rtyp =INT_CMD;
4507  res->data =(void *)(long)ifIso(p,q,f,g, a, b);
4508  }
4509  }
4510  }
4511  }
4512  }
4513  }
4514  return false;
4515 }
4516 
4518 {
4519  leftv h=args;
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() == INT_CMD))
4525  {
4526  int num= (int)(long)h->Data();
4527  res->rtyp =INT_CMD;
4528  res->data =(void *)(long)redefinedeg( p, num);
4529  }
4530  }
4531  return false;
4532 }
4533 
4535 {
4536  leftv h=args;
4537  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4538  {
4539  ideal h1= (ideal)h->Data();
4540  res->rtyp =IDEAL_CMD;
4541  res->data =complementsimplex(h1);
4542  }
4543  return false;
4544 }
4545 
4547 {
4548  leftv h=args;
4549  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4550  {
4551  ideal h1= (ideal)h->Data();
4552  res->rtyp =INT_CMD;
4553  res->data =(void *)(long)dim_sim(h1);
4554  }
4555  return false;
4556 }
4557 
4559 {
4560  leftv h=args;
4561  if((h != NULL)&&(h->Typ() == IDEAL_CMD))
4562  {
4563  ideal h1= (ideal)h->Data();
4564  h = h->next;
4565  if((h != NULL)&&(h->Typ() == INT_CMD))
4566  {
4567  int num= (int)(long)h->Data();
4568  res->rtyp =INT_CMD;
4569  res->data =(void *)(long)num4dim( h1, num);
4570  }
4571  }
4572  return false;
4573 }
4574 
4575 /**************************************interface T2****************************************/
4576 
4578 {
4579  p->iiAddCproc("","mg",FALSE,idsr);
4580  p->iiAddCproc("","gd",FALSE,gd);
4581  p->iiAddCproc("","findbset",FALSE,fb);
4582  p->iiAddCproc("","findaset",FALSE,fa);
4583  p->iiAddCproc("","fgp",FALSE,fgp);
4584  p->iiAddCproc("","fgpl",FALSE,fgpl);
4585  p->iiAddCproc("","idcomplement",FALSE,idcomplement);
4586  p->iiAddCproc("","genst",FALSE,genstt);
4587  p->iiAddCproc("","sgp",FALSE,sgp);
4588  p->iiAddCproc("","sgpl",FALSE,sgpl);
4589  p->iiAddCproc("","Links",FALSE,Links);
4590  p->iiAddCproc("","eqsolve1",FALSE,eqsolve1);
4591  p->iiAddCproc("","pb",FALSE,pb);
4592  p->iiAddCproc("","pa",FALSE,pa);
4593  p->iiAddCproc("","makeSimplex",FALSE,makeSimplex);
4594  p->iiAddCproc("","isSim",FALSE,isSim);
4595  p->iiAddCproc("","nfaces1",FALSE,nfaces1);
4596  p->iiAddCproc("","nfaces2",FALSE,nfaces2);
4597  p->iiAddCproc("","nfaces3",FALSE,nfaces3);
4598  p->iiAddCproc("","comedg",FALSE,comedg);
4599  p->iiAddCproc("","tsets",FALSE,tsets);
4600  p->iiAddCproc("","valency",FALSE,Valency);
4601  p->iiAddCproc("","nab",FALSE,nabvl);
4602  p->iiAddCproc("","tnab",FALSE,tnabvl);
4603  p->iiAddCproc("","mab",FALSE,mabvl);
4604  p->iiAddCproc("","SRideal",FALSE,SRideal);
4605  p->iiAddCproc("","Linkn",FALSE,linkn);
4606  p->iiAddCproc("","Existb",FALSE,existsub);
4607  p->iiAddCproc("","pConstant",FALSE,pConstant);
4608  p->iiAddCproc("","pChange",FALSE,pChange);
4609  p->iiAddCproc("","pNew",FALSE,p_New);
4610  p->iiAddCproc("","pSupport",FALSE,support);
4611  p->iiAddCproc("","psMinusp",FALSE,psMinusp);
4612  p->iiAddCproc("","cNew",FALSE,cNew);
4613  p->iiAddCproc("","isoNumber",FALSE,isoNumber);
4614  p->iiAddCproc("","vsInsec",FALSE,vsIntersec);
4615  p->iiAddCproc("","getnabt",FALSE,nabtvl);
4616  p->iiAddCproc("","idmodulo",FALSE,idModulo);
4617  p->iiAddCproc("","ndegree",FALSE,newDegree);
4618  p->iiAddCproc("","nonf2f",FALSE,nonf2f);
4619  p->iiAddCproc("","ifIsom",FALSE,ifIsomorphism);
4620  p->iiAddCproc("","stellarsubdivision",FALSE,stellarsubdivision);
4621  p->iiAddCproc("","star",FALSE,stars);
4622  p->iiAddCproc("","numdim",FALSE,numdim);
4623  p->iiAddCproc("","dimsim",FALSE,dimsim);
4624  p->iiAddCproc("","bprime",FALSE,bprime);
4625  p->iiAddCproc("","remainpart",FALSE,stellarremain);
4626  p->iiAddCproc("","idminus",FALSE,idminus);
4627  p->iiAddCproc("","time1",FALSE,t1h);
4628 }
4629 
4630 extern "C" int SI_MOD_INIT(cohomo)(SModulFunctions* p)
4631 {
4633  return MAX_TOK;
4634 }
4635 #endif
4636 #endif
int BOOLEAN
Definition: auxiliary.h:87
#define FALSE
Definition: auxiliary.h:96
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
CanonicalForm num(const CanonicalForm &f)
Variable mvar(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int p
Definition: cfModGcd.cc:4078
g
Definition: cfModGcd.cc:4090
CanonicalForm cf
Definition: cfModGcd.cc:4083
CanonicalForm gp
Definition: cfModGcd.cc:4102
CanonicalForm b
Definition: cfModGcd.cc:4103
bool solve(int **extmat, int nrows, int ncols)
Definition: cf_linsys.cc:504
FILE * f
Definition: checklibs.c:9
Definition: idrec.h:35
Definition: intvec.h:23
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
void * CopyD(int t)
Definition: subexpr.cc:710
int rtyp
Definition: subexpr.h:91
void * Data()
Definition: subexpr.cc:1154
void Init()
Definition: subexpr.h:107
void * data
Definition: subexpr.h:88
Definition: lists.h:24
sleftv * m
Definition: lists.h:46
Coefficient rings, fields and other domains suitable for Singular polynomials.
static FORCE_INLINE long n_Int(number &n, const coeffs r)
conversion of n to an int; 0 if not possible in Z/pZ: the representing int lying in (-p/2 ....
Definition: coeffs.h:544
static FORCE_INLINE coeffs nCopyCoeff(const coeffs r)
"copy" coeffs, i.e. increment ref
Definition: coeffs.h:430
static BOOLEAN fa(leftv res, leftv args)
Definition: cohomo.cc:3764
static BOOLEAN pa(leftv res, leftv args)
Definition: cohomo.cc:3723
static BOOLEAN tsets(leftv res, leftv args)
Definition: cohomo.cc:4053
static ideal T_1h(ideal h)
Definition: cohomo.cc:3549
static std::vector< int > eli1(std::vector< int > eq1, std::vector< int > eq2)
Definition: cohomo.cc:996
static std::vector< std::vector< int > > b_subsets(std::vector< int > vec)
Definition: cohomo.cc:496
static bool IsinL(int a, std::vector< int > vec)
Definition: cohomo.cc:132
static std::vector< int > vertset(std::vector< std::vector< int > > vecs)
Definition: cohomo.cc:1404
static BOOLEAN tnabvl(leftv res, leftv args)
Definition: cohomo.cc:4106
BOOLEAN nfaces1(leftv res, leftv args)
Definition: cohomo.cc:3931
static std::vector< int > subspacet1(int num, std::vector< std::vector< int > > ntvs)
Definition: cohomo.cc:1942
VAR clock_t t_construct
Definition: cohomo.cc:2692
static BOOLEAN cNew(leftv res, leftv args)
Definition: cohomo.cc:4368
static BOOLEAN comedg(leftv res, leftv args)
Definition: cohomo.cc:3694
static poly pMaken(std::vector< int > vbase)
Definition: cohomo.cc:465
VAR clock_t t_start
Definition: cohomo.cc:2692
static std::vector< int > make1(int n)
Definition: cohomo.cc:1198
static std::vector< std::vector< int > > vecqring(std::vector< std::vector< int > > vec1, std::vector< std::vector< int > > vec2)
Definition: cohomo.cc:454
static BOOLEAN idModulo(leftv res, leftv args)
Definition: cohomo.cc:4419
static ideal idMaken(std::vector< std::vector< int > > vecs)
Definition: cohomo.cc:479
static std::vector< int > makeequation(int i, int j, int t)
Definition: cohomo.cc:1570
static intvec * dmat(poly a, poly b)
Definition: cohomo.cc:3659
static ideal idMake(std::vector< std::vector< int > > vecs)
Definition: cohomo.cc:358
static std::vector< std::vector< int > > triface(poly p, int vert)
Definition: cohomo.cc:2985
static std::vector< std::vector< int > > boundary(poly a)
Definition: cohomo.cc:3489
static std::vector< std::vector< int > > vsMinusv(std::vector< std::vector< int > > vecs, std::vector< int > vec)
Definition: cohomo.cc:225
static std::vector< int > vecbase1(int num, std::vector< int > oset)
Definition: cohomo.cc:1169
static BOOLEAN isoNumber(leftv res, leftv args)
Definition: cohomo.cc:4453
static std::vector< std::vector< int > > subspacet(std::vector< std::vector< int > > mv, std::vector< int > bv, std::vector< std::vector< int > > ntvs)
Definition: cohomo.cc:1969
static std::vector< std::vector< int > > nabtv(std::vector< std::vector< int > > hvs, std::vector< std::vector< int > > Nv, std::vector< int > av, std::vector< int > bv)
Definition: cohomo.cc:2135
static std::vector< int > gdegree(poly a, poly b)
Definition: cohomo.cc:3456
static BOOLEAN makeSimplex(leftv res, leftv args)
Definition: cohomo.cc:3735
static bool vInvsl(std::vector< int > vec, std::vector< std::vector< int > > vecs)
Definition: cohomo.cc:186
static std::vector< std::vector< int > > supports2(ideal h)
Definition: cohomo.cc:326
static BOOLEAN fgp(leftv res, leftv args)
Definition: cohomo.cc:3786
static std::vector< std::vector< int > > getvector(ideal h, int n)
Definition: cohomo.cc:1863
static bool tNab(std::vector< std::vector< int > > hvs, std::vector< int > pv, std::vector< std::vector< int > > bvs)
Definition: cohomo.cc:2157
static int valency(ideal h, poly p)
Definition: cohomo.cc:3178
static bool condition3for2(std::vector< std::vector< int > > hvs, std::vector< int > pv, std::vector< int > qv, std::vector< int > av, std::vector< int > bv)
Definition: cohomo.cc:1777
static std::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)
Definition: cohomo.cc:2659
static std::vector< int > vMake(poly p)
Definition: cohomo.cc:420
static void lpsprint(std::vector< std::vector< poly > > pvs)
Definition: cohomo.cc:112
static BOOLEAN stars(leftv res, leftv args)
Definition: cohomo.cc:4385
static BOOLEAN newDegree(leftv res, leftv args)
Definition: cohomo.cc:4517
static std::vector< int > vecUnion(std::vector< int > vec1, std::vector< int > vec2)
Definition: cohomo.cc:200
static std::vector< int > keeporder(std::vector< int > vec)
Definition: cohomo.cc:1039
static void T1(ideal h)
Definition: cohomo.cc:2380
static bool condition2for2(std::vector< std::vector< int > > hvs, std::vector< int > pv, std::vector< int > qv, std::vector< int > sv, std::vector< int > av, std::vector< int > bv)
Definition: cohomo.cc:1764
static std::vector< std::vector< int > > penface(poly p, poly q, poly g, int vert)
Definition: cohomo.cc:3109
static std::vector< std::vector< int > > vsUnion(std::vector< std::vector< int > > vs1, std::vector< std::vector< int > > vs2)
Definition: cohomo.cc:239
static ideal trisets(ideal h)
Definition: cohomo.cc:2969
static int pvert(poly p)
Definition: cohomo.cc:541
static int idvert(ideal h)
Definition: cohomo.cc:522
static std::vector< int > phimagel(std::vector< int > fv, std::vector< int > av, std::vector< int > bv)
Definition: cohomo.cc:2649
static std::vector< int > subspace1(std::vector< std::vector< int > > mv, std::vector< int > bv)
Definition: cohomo.cc:1640
static std::vector< int > fvarsvalue(int vnum, std::vector< int > fvars)
Definition: cohomo.cc:1108
static std::vector< std::vector< int > > stellarsub(poly a, ideal h)
Definition: cohomo.cc:3499
static void firstorderdef_setup(SModulFunctions *p)
Definition: cohomo.cc:4577
static std::vector< std::vector< int > > gpl2(ideal h, poly a, poly b)
Definition: cohomo.cc:2845
static ideal SimFacset(poly p)
Definition: cohomo.cc:787
static ideal psubset(poly p)
Definition: cohomo.cc:1539
static intvec * Tmat(std::vector< std::vector< int > > vecs)
Definition: cohomo.cc:2233
static BOOLEAN bprime(leftv res, leftv args)
Definition: cohomo.cc:4319
static void listprint(std::vector< int > vec)
Definition: cohomo.cc:49
static BOOLEAN pConstant(leftv res, leftv args)
Definition: cohomo.cc:4261
static bool vInp(int m, poly p)
Definition: cohomo.cc:405
static BOOLEAN dimsim(leftv res, leftv args)
Definition: cohomo.cc:4546
static std::vector< int > tnab(std::vector< std::vector< int > > hvs, std::vector< std::vector< int > > nvs, std::vector< std::vector< int > > bvs)
Definition: cohomo.cc:2172
static BOOLEAN fgpl(leftv res, leftv args)
Definition: cohomo.cc:3808
static bool p_Ifsfree(poly P)
Definition: cohomo.cc:643
static int redefinedeg(poly p, int num)
Definition: cohomo.cc:1315
static std::vector< int > make0(int n)
Definition: cohomo.cc:1185
static BOOLEAN psMinusp(leftv res, leftv args)
Definition: cohomo.cc:4331
static bool nabtconditionv(std::vector< std::vector< int > > hvs, std::vector< int > pv, std::vector< int > qv)
Definition: cohomo.cc:2123
static BOOLEAN SRideal(leftv res, leftv args)
Definition: cohomo.cc:3600
static std::vector< std::vector< int > > p_constant(ideal Xo, ideal Sigma)
Definition: cohomo.cc:3247
static std::vector< std::vector< int > > soleli1(std::vector< std::vector< int > > eqs)
Definition: cohomo.cc:1052
static std::vector< std::vector< int > > vs_subsets(std::vector< std::vector< int > > vs)
Definition: cohomo.cc:3236
static std::vector< std::vector< int > > minisolve(std::vector< std::vector< int > > solve, std::vector< int > index)
Definition: cohomo.cc:2291
static bool condition2for2nv(std::vector< std::vector< int > > hvs, std::vector< int > pv, std::vector< int > qv, std::vector< int > fv)
Definition: cohomo.cc:2410
static BOOLEAN stellarsubdivision(leftv res, leftv args)
Definition: cohomo.cc:4402
static intvec * edgemat(poly p, poly q)
Definition: cohomo.cc:3053
static std::vector< std::vector< int > > links_new(poly a, ideal Xo, ideal Sigma, int vert, int ord)
Definition: cohomo.cc:3361
static BOOLEAN idcomplement(leftv res, leftv args)
Definition: cohomo.cc:3612
static std::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)
Definition: cohomo.cc:2196
VAR clock_t t_begin
Definition: cohomo.cc:2692
static int isoNum(poly p, ideal I, poly a, poly b)
Definition: cohomo.cc:3401
static bool condition1for2(std::vector< int > pv, std::vector< int > qv, std::vector< int > bv)
Definition: cohomo.cc:1751
static BOOLEAN isSim(leftv res, leftv args)
Definition: cohomo.cc:3919
static ideal IsSimplex(ideal h)
Definition: cohomo.cc:832
static BOOLEAN gd(leftv res, leftv args)
Definition: cohomo.cc:3677
static BOOLEAN nabvl(leftv res, leftv args)
Definition: cohomo.cc:4082
static BOOLEAN fb(leftv res, leftv args)
Definition: cohomo.cc:3711
static ideal idmodulo(ideal h1, ideal h2)
Definition: cohomo.cc:375
VAR clock_t t_mark
Definition: cohomo.cc:2692
static ideal idadda(ideal h1, ideal h2)
Definition: cohomo.cc:808
static ideal idMake3(std::vector< std::vector< int > > vecs)
Definition: cohomo.cc:1601
static ideal complementsimplex(ideal h)
Definition: cohomo.cc:853
static ideal c_New(ideal Io, ideal sig)
Definition: cohomo.cc:3293
static std::vector< std::vector< int > > canonicalbase(int n)
Definition: cohomo.cc:1843
static std::vector< std::vector< int > > eli2(int num, std::vector< int > bset, std::vector< std::vector< int > > gset)
Definition: cohomo.cc:1254
static void T2(ideal h)
Definition: cohomo.cc:2607
static ideal qringadd(ideal h1, ideal h2, int deg)
Definition: cohomo.cc:735
static std::vector< std::vector< int > > tetraface(poly p, poly q, int vert)
Definition: cohomo.cc:3072
static BOOLEAN eqsolve1(leftv res, leftv args)
Definition: cohomo.cc:4012
static ideal finda(ideal h, poly S, int ddeg)
Definition: cohomo.cc:934
static BOOLEAN numdim(leftv res, leftv args)
Definition: cohomo.cc:4558
static std::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)
Definition: cohomo.cc:2792
static bool vEvl(std::vector< int > vec1, std::vector< int > vec2)
Definition: cohomo.cc:174
static BOOLEAN vsIntersec(leftv res, leftv args)
Definition: cohomo.cc:4138
static std::vector< std::vector< int > > gpl(ideal h, poly a, poly b)
Definition: cohomo.cc:2704
static BOOLEAN support(leftv res, leftv args)
Definition: cohomo.cc:4307
static void id_print(ideal h)
Definition: cohomo.cc:84
static std::vector< int > gensindex(ideal M, ideal ids)
Definition: cohomo.cc:2260
static BOOLEAN genstt(leftv res, leftv args)
Definition: cohomo.cc:3835
static std::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)
Definition: cohomo.cc:2459
static std::vector< std::vector< int > > listsinsertlist(std::vector< std::vector< int > > gset, int a, int b)
Definition: cohomo.cc:1561
static void TimeShow(clock_t t_construct, clock_t t_solve, clock_t t_value, clock_t t_total)
Definition: cohomo.cc:2695
static BOOLEAN nonf2f(leftv res, leftv args)
Definition: cohomo.cc:4534
static void equmab(int num)
Definition: cohomo.cc:1618
static BOOLEAN sgpl(leftv res, leftv args)
Definition: cohomo.cc:3879
static bool nabconditionv(std::vector< std::vector< int > > hvs, std::vector< int > pv, std::vector< int > av, std::vector< int > bv)
Definition: cohomo.cc:2092
static std::vector< int > support1(poly p)
Definition: cohomo.cc:271
static std::vector< int > findalpha(std::vector< std::vector< int > > mv, std::vector< int > bv)
Definition: cohomo.cc:1927
static ideal id_sfmon(ideal h)
Definition: cohomo.cc:679
static ideal triangulations1(ideal h, poly p, int vert)
Definition: cohomo.cc:3004
VAR clock_t t_value
Definition: cohomo.cc:2692
static std::vector< int > vecMinus(std::vector< int > vec1, std::vector< int > vec2)
Definition: cohomo.cc:212
static int id_maxdeg(ideal h)
Definition: cohomo.cc:744
static std::vector< int > numfree(ideal h)
Definition: cohomo.cc:1824
static std::vector< std::vector< poly > > idMakei(std::vector< std::vector< int > > mv, std::vector< std::vector< int > > vecs)
Definition: cohomo.cc:1671
static std::vector< std::vector< int > > vAbsorb(std::vector< int > bset, std::vector< std::vector< int > > gset)
Definition: cohomo.cc:1125
static ideal id_complement(ideal h)
Definition: cohomo.cc:698
static ideal p_a(ideal h)
Definition: cohomo.cc:1335
static ideal triangulations3(ideal h, poly p, poly q, poly g, int vert)
Definition: cohomo.cc:3153
static BOOLEAN sgp(leftv res, leftv args)
Definition: cohomo.cc:3857
static std::vector< int > support2(poly p)
Definition: cohomo.cc:304
static BOOLEAN idsr(leftv res, leftv args)
Definition: cohomo.cc:3637
static ideal genst(ideal h, poly a, poly b)
Definition: cohomo.cc:2511
static BOOLEAN linkn(leftv res, leftv args)
Definition: cohomo.cc:4212
static std::vector< std::vector< int > > p_change(ideal Sigma)
Definition: cohomo.cc:3255
static std::vector< std::vector< int > > bsubsets_1(poly b)
Definition: cohomo.cc:3532
static intvec * gradedpiece1nl(ideal h, poly a, poly b, int set)
Definition: cohomo.cc:2768
static std::vector< int > v_minus(std::vector< int > v1, std::vector< int > v2)
Definition: cohomo.cc:3446
static std::vector< std::vector< int > > phi2(poly a, ideal Xo, ideal Sigma)
Definition: cohomo.cc:3344
static BOOLEAN Valency(leftv res, leftv args)
Definition: cohomo.cc:4065
static bool mabconditionv(std::vector< std::vector< int > > hvs, std::vector< int > pv, std::vector< int > av, std::vector< int > bv)
Definition: cohomo.cc:963
static ideal idMinus(ideal h1, ideal h2)
Definition: cohomo.cc:617
int SI_MOD_INIT() cohomo(SModulFunctions *p)
Definition: cohomo.cc:4630
static int dim_sim(ideal h)
Definition: cohomo.cc:876
static std::vector< int > ofindbases1(int num, int vnum, std::vector< int > bset, std::vector< std::vector< int > > gset)
Definition: cohomo.cc:1211
static intvec * gradedpiece1n(ideal h, poly a, poly b)
Definition: cohomo.cc:2313
static BOOLEAN stellarremain(leftv res, leftv args)
Definition: cohomo.cc:4348
static void lpprint(std::vector< poly > pv)
Definition: cohomo.cc:97
static ideal getpresolve(ideal h)
Definition: cohomo.cc:1796
static BOOLEAN nabtvl(leftv res, leftv args)
Definition: cohomo.cc:4178
static std::vector< std::vector< int > > vsIntersection(std::vector< std::vector< int > > vs1, std::vector< std::vector< int > > vs2)
Definition: cohomo.cc:253
static std::vector< std::vector< int > > subspacetn(std::vector< std::vector< int > > N, std::vector< int > tN, std::vector< std::vector< int > > ntvs)
Definition: cohomo.cc:2440
static int pcoef(poly p, int m)
Definition: cohomo.cc:386
static ideal triangulations2(ideal h, poly p, poly q, int vert)
Definition: cohomo.cc:3094
VAR clock_t t_total
Definition: cohomo.cc:2692
static std::vector< int > findalphan(std::vector< std::vector< int > > N, std::vector< int > tN)
Definition: cohomo.cc:2424
static BOOLEAN nfaces3(leftv res, leftv args)
Definition: cohomo.cc:3980
static std::vector< int > commonedge(poly p, poly q)
Definition: cohomo.cc:3042
static std::vector< std::vector< int > > vsMake(ideal h)
Definition: cohomo.cc:439
static void gradedpiece1(ideal h, poly a, poly b)
Definition: cohomo.cc:1691
static BOOLEAN p_New(leftv res, leftv args)
Definition: cohomo.cc:4290
static BOOLEAN Links(leftv res, leftv args)
Definition: cohomo.cc:3901
static std::vector< std::vector< int > > Mabv(ideal h, poly a, poly b)
Definition: cohomo.cc:975
static poly pMake(std::vector< int > vbase)
Definition: cohomo.cc:343
static BOOLEAN nfaces2(leftv res, leftv args)
Definition: cohomo.cc:3953
static bool IsInX(poly p, ideal X)
Definition: cohomo.cc:719
static std::vector< std::vector< int > > mabtv(std::vector< std::vector< int > > hvs, std::vector< std::vector< int > > Mv, std::vector< int > av, std::vector< int > bv)
Definition: cohomo.cc:1983
static intvec * gradedpiece2nl(ideal h, poly a, poly b)
Definition: cohomo.cc:2913
static std::vector< std::vector< int > > supports(ideal h)
Definition: cohomo.cc:287
static BOOLEAN pb(leftv res, leftv args)
Definition: cohomo.cc:3747
VAR clock_t t_solve
Definition: cohomo.cc:2692
static std::vector< std::vector< int > > ofindbases(int num, std::vector< int > bset, std::vector< std::vector< int > > gset)
Definition: cohomo.cc:1227
static std::vector< int > vecIntersection(std::vector< int > p, std::vector< int > q)
Definition: cohomo.cc:147
static ideal sfreemon(ideal h, int deg)
Definition: cohomo.cc:658
static int existIn(poly b, ideal Xs)
Definition: cohomo.cc:3386
static std::vector< poly > pMakei(std::vector< std::vector< int > > mv, std::vector< int > vbase)
Definition: cohomo.cc:1657
static BOOLEAN pChange(leftv res, leftv args)
Definition: cohomo.cc:4278
static std::vector< std::vector< int > > Nabv(std::vector< std::vector< int > > hvs, std::vector< int > av, std::vector< int > bv)
Definition: cohomo.cc:2106
static std::vector< int > phimage(std::vector< int > pv, std::vector< int > av, std::vector< int > bv)
Definition: cohomo.cc:2187
static ideal idsrRing(ideal h)
Definition: cohomo.cc:758
static ideal idMinusp(ideal I, poly p)
Definition: cohomo.cc:3429
static BOOLEAN idminus(leftv res, leftv args)
Definition: cohomo.cc:4436
static std::vector< std::vector< int > > phi1(poly a, ideal Sigma)
Definition: cohomo.cc:3328
static ideal p_b(ideal h, poly a)
Definition: cohomo.cc:1427
static void listsprint(std::vector< std::vector< int > > posMat)
Definition: cohomo.cc:65
static BOOLEAN mabvl(leftv res, leftv args)
Definition: cohomo.cc:4156
static int num4dim(ideal h, int n)
Definition: cohomo.cc:889
static std::vector< int > freevars(int n, std::vector< int > bset, std::vector< std::vector< int > > gset)
Definition: cohomo.cc:1085
static std::vector< std::vector< int > > links(poly a, ideal h)
Definition: cohomo.cc:1293
static BOOLEAN ifIsomorphism(leftv res, leftv args)
Definition: cohomo.cc:4480
static int ifIso(poly p, poly q, poly f, poly g, poly a, poly b)
Definition: cohomo.cc:3418
static void gradedpiece2(ideal h, poly a, poly b)
Definition: cohomo.cc:2006
static std::vector< std::vector< int > > p_new(ideal Xo, ideal Sigma)
Definition: cohomo.cc:3262
static BOOLEAN existsub(leftv res, leftv args)
Definition: cohomo.cc:4244
static ideal mingens(ideal h, poly a, poly b)
Definition: cohomo.cc:2277
static bool vsubset(std::vector< int > vec1, std::vector< int > vec2)
Definition: cohomo.cc:160
static ideal findb(ideal h)
Definition: cohomo.cc:908
static BOOLEAN t1h(leftv res, leftv args)
Definition: cohomo.cc:3625
static std::vector< std::vector< int > > star(poly a, ideal h)
Definition: cohomo.cc:3473
static std::vector< std::vector< int > > vsMinusvs(std::vector< std::vector< int > > vs1, std::vector< std::vector< int > > vs2)
Definition: cohomo.cc:3226
static poly pMake3(std::vector< int > vbase)
Definition: cohomo.cc:1583
static intvec * gradedpiece2n(ideal h, poly a, poly b)
Definition: cohomo.cc:2528
#define Print
Definition: emacs.cc:80
const CanonicalForm int s
Definition: facAbsFact.cc:51
CanonicalForm res
Definition: facAbsFact.cc:60
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
bool bad
Definition: facFactorize.cc:64
fq_nmod_poly_t * vec
Definition: facHensel.cc:108
int j
Definition: facHensel.cc:110
static int max(int a, int b)
Definition: fast_mult.cc:264
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define VAR
Definition: globaldefs.h:5
@ IDEAL_CMD
Definition: grammar.cc:284
@ POLY_CMD
Definition: grammar.cc:289
@ RING_CMD
Definition: grammar.cc:281
ideal scKBase(int deg, ideal s, ideal Q, intvec *mv)
Definition: hdegree.cc:1449
BOOLEAN idInsertPoly(ideal h1, poly h2)
insert h2 into h1 (if h2 is not the zero polynomial) return TRUE iff h2 was indeed inserted
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
ideal idCopy(ideal A)
Definition: ideals.h:60
ideal idAdd(ideal h1, ideal h2)
h1 + h2
Definition: ideals.h:68
#define IMATELEM(M, I, J)
Definition: intvec.h:85
idhdl ggetid(const char *n)
Definition: ipid.cc:581
idhdl enterid(const char *s, int lev, int t, idhdl *root, BOOLEAN init, BOOLEAN search)
Definition: ipid.cc:279
#define IDROOT
Definition: ipid.h:19
#define IDRING(a)
Definition: ipid.h:127
BOOLEAN iiMake_proc(idhdl pn, package pack, leftv args)
Definition: iplib.cc:504
INST_VAR sleftv iiRETURNEXPR
Definition: iplib.cc:474
void rSetHdl(idhdl h)
Definition: ipshell.cc:5126
STATIC_VAR Poly * h
Definition: janet.cc:971
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:3182
ideal kStd(ideal F, ideal Q, tHomog h, intvec **w, intvec *hilb, int syzComp, int newIdeal, intvec *vw, s_poly_proc_t sp)
Definition: kstd1.cc:2447
void rem(unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
Definition: minpoly.cc:572
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
slists * lists
Definition: mpr_numeric.h:146
char N base
Definition: ValueTraits.h:144
The main handler for Singular numbers which are suitable for Singular polynomials.
#define nInit(i)
Definition: numbers.h:24
#define omStrDup(s)
Definition: omAllocDecl.h:263
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omalloc(size)
Definition: omAllocDecl.h:228
#define NULL
Definition: omList.c:12
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
poly p_Subst(poly p, int n, poly e, const ring r)
Definition: p_polys.cc:3954
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition: p_polys.cc:4508
static poly pp_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:1029
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent @Note: the integer VarOffset encodes:
Definition: p_polys.h:467
static BOOLEAN p_DivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1898
void rChangeCurrRing(ring r)
Definition: polys.cc:15
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define pAdd(p, q)
Definition: polys.h:203
static long pTotaldegree(poly p)
Definition: polys.h:282
#define pSetm(p)
Definition: polys.h:271
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
void pWrite(poly p)
Definition: polys.h:308
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pEqualPolys(p1, p2)
Definition: polys.h:399
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
ring rDefault(const coeffs cf, int N, char **n, int ord_size, rRingOrder_t *ord, int *block0, int *block1, int **wvhdl, unsigned long bitmask)
Definition: ring.cc:102
ring rCopy(ring r)
Definition: ring.cc:1731
@ ringorder_lp
Definition: ring.h:77
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:592
ideal id_Add(ideal h1, ideal h2, const ring r)
h1 + h2
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
ideal id_MaxIdeal(const ring r)
initialise the maximal ideal (at 0)
Definition: simpleideals.cc:98
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
#define IDELEMS(i)
Definition: simpleideals.h:23
#define R
Definition: sirandom.c:27
#define M
Definition: sirandom.c:25
@ testHomog
Definition: structs.h:38
#define assert(A)
Definition: svd_si.h:3
@ INTVEC_CMD
Definition: tok.h:101
@ INT_CMD
Definition: tok.h:96
@ MAX_TOK
Definition: tok.h:218
int dim(ideal I, ring r)
int tdeg(poly p)
Definition: walkSupport.cc:35