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