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