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