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