My Project
hilb.cc
Go to the documentation of this file.
1/****************************************
2* Computer Algebra System SINGULAR *
3****************************************/
4/*
5* ABSTRACT - Hilbert series
6*/
7
8#include "kernel/mod2.h"
9
10#include "misc/mylimits.h"
11#include "misc/intvec.h"
12
16
19#include "polys/simpleideals.h"
20
21
22// #include "kernel/structs.h"
23// #include "kernel/polys.h"
24//ADICHANGES:
25#include "kernel/ideals.h"
26// #include "kernel/GBEngine/kstd1.h"
27
28# define omsai 1
29#if omsai
30
32#include "coeffs/coeffs.h"
34#include "coeffs/numbers.h"
35#include <vector>
36#include "Singular/ipshell.h"
37
38
39#include <Singular/ipshell.h>
40#include <ctime>
41#include <iostream>
42#endif
43
47
48
49static int hMinModulweight(intvec *modulweight)
50{
51 int i,j,k;
52
53 if(modulweight==NULL) return 0;
54 j=(*modulweight)[0];
55 for(i=modulweight->rows()-1;i!=0;i--)
56 {
57 k=(*modulweight)[i];
58 if(k<j) j=k;
59 }
60 return j;
61}
62
63static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
64{
65 int i, j;
66 int x, y, z = 1;
67 int *p;
68 for (i = Nvar; i>0; i--)
69 {
70 x = 0;
71 for (j = 0; j < Nstc; j++)
72 {
73 y = stc[j][var[i]];
74 if (y > x)
75 x = y;
76 }
77 z += x;
78 j = i - 1;
79 if (z > Ql[j])
80 {
81 if (z>(MAX_INT_VAL)/2)
82 {
83 WerrorS("internal arrays too big");
84 return;
85 }
86 p = (int *)omAlloc((unsigned long)z * sizeof(int));
87 if (Ql[j]!=0)
88 {
89 if (j==0)
90 memcpy(p, Qpol[j], Ql[j] * sizeof(int));
91 omFreeSize((ADDRESS)Qpol[j], Ql[j] * sizeof(int));
92 }
93 if (j==0)
94 {
95 for (x = Ql[j]; x < z; x++)
96 p[x] = 0;
97 }
98 Ql[j] = z;
99 Qpol[j] = p;
100 }
101 }
102}
103
104static int *hAddHilb(int Nv, int x, int *pol, int *lp)
105{
106 int l = *lp, ln, i;
107 int *pon;
108 *lp = ln = l + x;
109 pon = Qpol[Nv];
110 memcpy(pon, pol, l * sizeof(int));
111 if (l > x)
112 {/*pon[i] -= pol[i - x];*/
113 for (i = x; i < l; i++)
114 { int64 t=pon[i];
115 int64 t2=pol[i - x];
116 t-=t2;
117 if ((t>=INT_MIN)&&(t<=INT_MAX)) pon[i]=t;
118 else if (!errorreported) WerrorS("int overflow in hilb 1");
119 }
120 for (i = l; i < ln; i++)
121 { /*pon[i] = -pol[i - x];*/
122 int64 t= -pol[i - x];
123 if ((t>=INT_MIN)&&(t<=INT_MAX)) pon[i]=t;
124 else if (!errorreported) WerrorS("int overflow in hilb 2");
125 }
126 }
127 else
128 {
129 for (i = l; i < x; i++)
130 pon[i] = 0;
131 for (i = x; i < ln; i++)
132 pon[i] = -pol[i - x];
133 }
134 return pon;
135}
136
137static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
138{
139 int l = lp, x, i, j;
140 int *pl;
141 int *p;
142 p = pol;
143 for (i = Nv; i>0; i--)
144 {
145 x = pure[var[i + 1]];
146 if (x!=0)
147 p = hAddHilb(i, x, p, &l);
148 }
149 pl = *Qpol;
150 j = Q0[Nv + 1];
151 for (i = 0; i < l; i++)
152 { /* pl[i + j] += p[i];*/
153 int64 t=pl[i+j];
154 int64 t2=p[i];
155 t+=t2;
156 if ((t>=INT_MIN)&&(t<=INT_MAX)) pl[i+j]=t;
157 else if (!errorreported) WerrorS("int overflow in hilb 3");
158 }
159 x = pure[var[1]];
160 if (x!=0)
161 {
162 j += x;
163 for (i = 0; i < l; i++)
164 { /* pl[i + j] -= p[i];*/
165 int64 t=pl[i+j];
166 int64 t2=p[i];
167 t-=t2;
168 if ((t>=INT_MIN)&&(t<=INT_MAX)) pl[i+j]=t;
169 else if (!errorreported) WerrorS("int overflow in hilb 4");
170 }
171 }
172 j += l;
173 if (j > hLength)
174 hLength = j;
175}
176
177static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var,
178 int Nvar, int *pol, int Lpol)
179{
180 int iv = Nvar -1, ln, a, a0, a1, b, i;
181 int x, x0;
182 scmon pn;
183 scfmon sn;
184 int *pon;
185 if (Nstc==0)
186 {
187 hLastHilb(pure, iv, var, pol, Lpol);
188 return;
189 }
190 x = a = 0;
191 pn = hGetpure(pure);
192 sn = hGetmem(Nstc, stc, stcmem[iv]);
193 hStepS(sn, Nstc, var, Nvar, &a, &x);
194 Q0[iv] = Q0[Nvar];
195 ln = Lpol;
196 pon = pol;
197 if (a == Nstc)
198 {
199 x = pure[var[Nvar]];
200 if (x!=0)
201 pon = hAddHilb(iv, x, pon, &ln);
202 hHilbStep(pn, sn, a, var, iv, pon, ln);
203 return;
204 }
205 else
206 {
207 pon = hAddHilb(iv, x, pon, &ln);
208 hHilbStep(pn, sn, a, var, iv, pon, ln);
209 }
210 b = a;
211 x0 = 0;
212 loop
213 {
214 Q0[iv] += (x - x0);
215 a0 = a;
216 x0 = x;
217 hStepS(sn, Nstc, var, Nvar, &a, &x);
218 hElimS(sn, &b, a0, a, var, iv);
219 a1 = a;
220 hPure(sn, a0, &a1, var, iv, pn, &i);
221 hLex2S(sn, b, a0, a1, var, iv, hwork);
222 b += (a1 - a0);
223 ln = Lpol;
224 if (a < Nstc)
225 {
226 pon = hAddHilb(iv, x - x0, pol, &ln);
227 hHilbStep(pn, sn, b, var, iv, pon, ln);
228 }
229 else
230 {
231 x = pure[var[Nvar]];
232 if (x!=0)
233 pon = hAddHilb(iv, x - x0, pol, &ln);
234 else
235 pon = pol;
236 hHilbStep(pn, sn, b, var, iv, pon, ln);
237 return;
238 }
239 }
240}
241
242/*
243*basic routines
244*/
245static void hWDegree(intvec *wdegree)
246{
247 int i, k;
248 int x;
249
250 for (i=(currRing->N); i; i--)
251 {
252 x = (*wdegree)[i-1];
253 if (x != 1)
254 {
255 for (k=hNexist-1; k>=0; k--)
256 {
257 hexist[k][i] *= x;
258 }
259 }
260 }
261}
262// ---------------------------------- ADICHANGES ---------------------------------------------
263//!!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
264
265#if 0 // unused
266//Tests if the ideal is sorted by degree
267static bool idDegSortTest(ideal I)
268{
269 if((I == NULL)||(idIs0(I)))
270 {
271 return(TRUE);
272 }
273 for(int i = 0; i<IDELEMS(I)-1; i++)
274 {
275 if(p_Totaldegree(I->m[i],currRing)>p_Totaldegree(I->m[i+1],currRing))
276 {
277 idPrint(I);
278 WerrorS("Ideal is not deg sorted!!");
279 return(FALSE);
280 }
281 }
282 return(TRUE);
283}
284#endif
285
286//adds the new polynomial at the coresponding position
287//and simplifies the ideal, destroys p
288static void SortByDeg_p(ideal I, poly p)
289{
290 int i,j;
291 if(idIs0(I))
292 {
293 I->m[0] = p;
294 return;
295 }
296 idSkipZeroes(I);
297 #if 1
298 for(i = 0; (i<IDELEMS(I)) && (p_Totaldegree(I->m[i],currRing)<=p_Totaldegree(p,currRing)); i++)
299 {
300 if(p_DivisibleBy( I->m[i],p, currRing))
301 {
303 return;
304 }
305 }
306 for(i = IDELEMS(I)-1; (i>=0) && (p_Totaldegree(I->m[i],currRing)>=p_Totaldegree(p,currRing)); i--)
307 {
308 if(p_DivisibleBy(p,I->m[i], currRing))
309 {
310 p_Delete(&I->m[i],currRing);
311 }
312 }
313 if(idIs0(I))
314 {
315 idSkipZeroes(I);
316 I->m[0] = p;
317 return;
318 }
319 #endif
320 idSkipZeroes(I);
321 //First I take the case when all generators have the same degree
322 if(p_Totaldegree(I->m[0],currRing) == p_Totaldegree(I->m[IDELEMS(I)-1],currRing))
323 {
325 {
326 idInsertPoly(I,p);
327 idSkipZeroes(I);
328 for(i=IDELEMS(I)-1;i>=1; i--)
329 {
330 I->m[i] = I->m[i-1];
331 }
332 I->m[0] = p;
333 return;
334 }
336 {
337 idInsertPoly(I,p);
338 idSkipZeroes(I);
339 return;
340 }
341 }
343 {
344 idInsertPoly(I,p);
345 idSkipZeroes(I);
346 for(i=IDELEMS(I)-1;i>=1; i--)
347 {
348 I->m[i] = I->m[i-1];
349 }
350 I->m[0] = p;
351 return;
352 }
354 {
355 idInsertPoly(I,p);
356 idSkipZeroes(I);
357 return;
358 }
359 for(i = IDELEMS(I)-2; ;)
360 {
362 {
363 idInsertPoly(I,p);
364 idSkipZeroes(I);
365 for(j = IDELEMS(I)-1; j>=i+1;j--)
366 {
367 I->m[j] = I->m[j-1];
368 }
369 I->m[i] = p;
370 return;
371 }
373 {
374 idInsertPoly(I,p);
375 idSkipZeroes(I);
376 for(j = IDELEMS(I)-1; j>=i+2;j--)
377 {
378 I->m[j] = I->m[j-1];
379 }
380 I->m[i+1] = p;
381 return;
382 }
383 i--;
384 }
385}
386
387//it sorts the ideal by the degrees
388static ideal SortByDeg(ideal I)
389{
390 if(idIs0(I))
391 {
392 return id_Copy(I,currRing);
393 }
394 int i;
395 ideal res;
396 idSkipZeroes(I);
397 res = idInit(1,1);
398 for(i = 0; i<=IDELEMS(I)-1;i++)
399 {
400 SortByDeg_p(res, I->m[i]);
401 I->m[i]=NULL; // I->m[i] is now in res
402 }
404 //idDegSortTest(res);
405 return(res);
406}
407
408//idQuot(I,p) for I monomial ideal, p a ideal with a single monomial.
409ideal idQuotMon(ideal Iorig, ideal p)
410{
411 if(idIs0(Iorig))
412 {
413 ideal res = idInit(1,1);
414 res->m[0] = poly(0);
415 return(res);
416 }
417 if(idIs0(p))
418 {
419 ideal res = idInit(1,1);
420 res->m[0] = pOne();
421 return(res);
422 }
423 ideal I = id_Head(Iorig,currRing);
424 ideal res = idInit(IDELEMS(I),1);
425 int i,j;
426 int dummy;
427 for(i = 0; i<IDELEMS(I); i++)
428 {
429 res->m[i] = p_Head(I->m[i], currRing);
430 for(j = 1; (j<=currRing->N) ; j++)
431 {
432 dummy = p_GetExp(p->m[0], j, currRing);
433 if(dummy > 0)
434 {
435 if(p_GetExp(I->m[i], j, currRing) < dummy)
436 {
437 p_SetExp(res->m[i], j, 0, currRing);
438 }
439 else
440 {
441 p_SetExp(res->m[i], j, p_GetExp(I->m[i], j, currRing) - dummy, currRing);
442 }
443 }
444 }
445 p_Setm(res->m[i], currRing);
447 {
448 p_Delete(&res->m[i],currRing);
449 }
450 else
451 {
452 p_Delete(&I->m[i],currRing);
453 }
454 }
456 idSkipZeroes(I);
457 if(!idIs0(res))
458 {
459 for(i = 0; i<=IDELEMS(res)-1; i++)
460 {
461 SortByDeg_p(I,res->m[i]);
462 res->m[i]=NULL; // is now in I
463 }
464 }
466 //idDegSortTest(I);
467 return(I);
468}
469
470//id_Add for monomials
471static void idAddMon(ideal I, ideal p)
472{
473 SortByDeg_p(I,p->m[0]);
474 p->m[0]=NULL; // is now in I
475 //idSkipZeroes(I);
476}
477
478//searches for a variable that is not yet used (assumes that the ideal is sqrfree)
479static poly ChoosePVar (ideal I)
480{
481 bool flag=TRUE;
482 int i,j;
483 poly res;
484 for(i=1;i<=currRing->N;i++)
485 {
486 flag=TRUE;
487 for(j=IDELEMS(I)-1;(j>=0)&&(flag);j--)
488 {
489 if(p_GetExp(I->m[j], i, currRing)>0)
490 {
491 flag=FALSE;
492 }
493 }
494
495 if(flag == TRUE)
496 {
497 res = p_ISet(1, currRing);
498 p_SetExp(res, i, 1, currRing);
500 return(res);
501 }
502 else
503 {
505 }
506 }
507 return(NULL); //i.e. it is the maximal ideal
508}
509
510#if 0 // unused
511//choice XL: last entry divided by x (xy10z15 -> y9z14)
512static poly ChoosePXL(ideal I)
513{
514 int i,j,dummy=0;
515 poly m;
516 for(i = IDELEMS(I)-1; (i>=0) && (dummy == 0); i--)
517 {
518 for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
519 {
520 if(p_GetExp(I->m[i],j, currRing)>1)
521 {
522 dummy = 1;
523 }
524 }
525 }
526 m = p_Copy(I->m[i+1],currRing);
527 for(j = 1; j<=currRing->N; j++)
528 {
529 dummy = p_GetExp(m,j,currRing);
530 if(dummy >= 1)
531 {
532 p_SetExp(m, j, dummy-1, currRing);
533 }
534 }
535 if(!p_IsOne(m, currRing))
536 {
537 p_Setm(m, currRing);
538 return(m);
539 }
540 m = ChoosePVar(I);
541 return(m);
542}
543#endif
544
545#if 0 // unused
546//choice XF: first entry divided by x (xy10z15 -> y9z14)
547static poly ChoosePXF(ideal I)
548{
549 int i,j,dummy=0;
550 poly m;
551 for(i =0 ; (i<=IDELEMS(I)-1) && (dummy == 0); i++)
552 {
553 for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
554 {
555 if(p_GetExp(I->m[i],j, currRing)>1)
556 {
557 dummy = 1;
558 }
559 }
560 }
561 m = p_Copy(I->m[i-1],currRing);
562 for(j = 1; j<=currRing->N; j++)
563 {
564 dummy = p_GetExp(m,j,currRing);
565 if(dummy >= 1)
566 {
567 p_SetExp(m, j, dummy-1, currRing);
568 }
569 }
570 if(!p_IsOne(m, currRing))
571 {
572 p_Setm(m, currRing);
573 return(m);
574 }
575 m = ChoosePVar(I);
576 return(m);
577}
578#endif
579
580#if 0 // unused
581//choice OL: last entry the first power (xy10z15 -> xy9z15)
582static poly ChoosePOL(ideal I)
583{
584 int i,j,dummy;
585 poly m;
586 for(i = IDELEMS(I)-1;i>=0;i--)
587 {
588 m = p_Copy(I->m[i],currRing);
589 for(j=1;j<=currRing->N;j++)
590 {
591 dummy = p_GetExp(m,j,currRing);
592 if(dummy > 0)
593 {
594 p_SetExp(m,j,dummy-1,currRing);
596 }
597 }
598 if(!p_IsOne(m, currRing))
599 {
600 return(m);
601 }
602 else
603 {
605 }
606 }
607 m = ChoosePVar(I);
608 return(m);
609}
610#endif
611
612#if 0 // unused
613//choice OF: first entry the first power (xy10z15 -> xy9z15)
614static poly ChoosePOF(ideal I)
615{
616 int i,j,dummy;
617 poly m;
618 for(i = 0 ;i<=IDELEMS(I)-1;i++)
619 {
620 m = p_Copy(I->m[i],currRing);
621 for(j=1;j<=currRing->N;j++)
622 {
623 dummy = p_GetExp(m,j,currRing);
624 if(dummy > 0)
625 {
626 p_SetExp(m,j,dummy-1,currRing);
628 }
629 }
630 if(!p_IsOne(m, currRing))
631 {
632 return(m);
633 }
634 else
635 {
637 }
638 }
639 m = ChoosePVar(I);
640 return(m);
641}
642#endif
643
644#if 0 // unused
645//choice VL: last entry the first variable with power (xy10z15 -> y)
646static poly ChoosePVL(ideal I)
647{
648 int i,j,dummy;
649 bool flag = TRUE;
650 poly m = p_ISet(1,currRing);
651 for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
652 {
653 flag = TRUE;
654 for(j=1;(j<=currRing->N) && (flag);j++)
655 {
656 dummy = p_GetExp(I->m[i],j,currRing);
657 if(dummy >= 2)
658 {
659 p_SetExp(m,j,1,currRing);
661 flag = FALSE;
662 }
663 }
664 if(!p_IsOne(m, currRing))
665 {
666 return(m);
667 }
668 }
669 m = ChoosePVar(I);
670 return(m);
671}
672#endif
673
674#if 0 // unused
675//choice VF: first entry the first variable with power (xy10z15 -> y)
676static poly ChoosePVF(ideal I)
677{
678 int i,j,dummy;
679 bool flag = TRUE;
680 poly m = p_ISet(1,currRing);
681 for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
682 {
683 flag = TRUE;
684 for(j=1;(j<=currRing->N) && (flag);j++)
685 {
686 dummy = p_GetExp(I->m[i],j,currRing);
687 if(dummy >= 2)
688 {
689 p_SetExp(m,j,1,currRing);
691 flag = FALSE;
692 }
693 }
694 if(!p_IsOne(m, currRing))
695 {
696 return(m);
697 }
698 }
699 m = ChoosePVar(I);
700 return(m);
701}
702#endif
703
704//choice JL: last entry just variable with power (xy10z15 -> y10)
705static poly ChoosePJL(ideal I)
706{
707 int i,j,dummy;
708 bool flag = TRUE;
709 poly m = p_ISet(1,currRing);
710 for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
711 {
712 flag = TRUE;
713 for(j=1;(j<=currRing->N) && (flag);j++)
714 {
715 dummy = p_GetExp(I->m[i],j,currRing);
716 if(dummy >= 2)
717 {
718 p_SetExp(m,j,dummy-1,currRing);
720 flag = FALSE;
721 }
722 }
723 if(!p_IsOne(m, currRing))
724 {
725 return(m);
726 }
727 }
729 m = ChoosePVar(I);
730 return(m);
731}
732
733#if 0
734//choice JF: last entry just variable with power -1 (xy10z15 -> y9)
735static poly ChoosePJF(ideal I)
736{
737 int i,j,dummy;
738 bool flag = TRUE;
739 poly m = p_ISet(1,currRing);
740 for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
741 {
742 flag = TRUE;
743 for(j=1;(j<=currRing->N) && (flag);j++)
744 {
745 dummy = p_GetExp(I->m[i],j,currRing);
746 if(dummy >= 2)
747 {
748 p_SetExp(m,j,dummy-1,currRing);
750 flag = FALSE;
751 }
752 }
753 if(!p_IsOne(m, currRing))
754 {
755 return(m);
756 }
757 }
758 m = ChoosePVar(I);
759 return(m);
760}
761#endif
762
763//chooses 1 \neq p \not\in S. This choice should be made optimal
764static poly ChooseP(ideal I)
765{
766 poly m;
767 // TEST TO SEE WHICH ONE IS BETTER
768 //m = ChoosePXL(I);
769 //m = ChoosePXF(I);
770 //m = ChoosePOL(I);
771 //m = ChoosePOF(I);
772 //m = ChoosePVL(I);
773 //m = ChoosePVF(I);
774 m = ChoosePJL(I);
775 //m = ChoosePJF(I);
776 return(m);
777}
778
779///searches for a monomial of degree d>=2 and divides it by a variable (result monomial of deg d-1)
780static poly SearchP(ideal I)
781{
782 int i,j,exp;
783 poly res;
784 if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)<=1)
785 {
786 res = ChoosePVar(I);
787 return(res);
788 }
789 i = IDELEMS(I)-1;
790 res = p_Copy(I->m[i], currRing);
791 for(j=1;j<=currRing->N;j++)
792 {
793 exp = p_GetExp(I->m[i], j, currRing);
794 if(exp > 0)
795 {
796 p_SetExp(res, j, exp - 1, currRing);
798 break;
799 }
800 }
801 assume( j <= currRing->N );
802 return(res);
803}
804
805//test if the ideal is of the form (x1, ..., xr)
806static bool JustVar(ideal I)
807{
808 #if 0
809 int i,j;
810 bool foundone;
811 for(i=0;i<=IDELEMS(I)-1;i++)
812 {
813 foundone = FALSE;
814 for(j = 1;j<=currRing->N;j++)
815 {
816 if(p_GetExp(I->m[i], j, currRing)>0)
817 {
818 if(foundone == TRUE)
819 {
820 return(FALSE);
821 }
822 foundone = TRUE;
823 }
824 }
825 }
826 return(TRUE);
827 #else
828 if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)>1)
829 {
830 return(FALSE);
831 }
832 return(TRUE);
833 #endif
834}
835
836//computes the Euler Characteristic of the ideal
837static void eulerchar (ideal I, int variables, mpz_ptr ec)
838{
839 loop
840 {
841 mpz_t dummy;
842 if(JustVar(I) == TRUE)
843 {
844 if(IDELEMS(I) == variables)
845 {
846 mpz_init(dummy);
847 if((variables % 2) == 0)
848 mpz_set_ui(dummy, 1);
849 else
850 mpz_set_si(dummy, -1);
851 mpz_add(ec, ec, dummy);
852 mpz_clear(dummy);
853 }
854 return;
855 }
856 ideal p = idInit(1,1);
857 p->m[0] = SearchP(I);
858 //idPrint(I);
859 //idPrint(p);
860 //printf("\nNow get in idQuotMon\n");
861 ideal Ip = idQuotMon(I,p);
862 //idPrint(Ip);
863 //Ip = SortByDeg(Ip);
864 int i,howmanyvarinp = 0;
865 for(i = 1;i<=currRing->N;i++)
866 {
867 if(p_GetExp(p->m[0],i,currRing)>0)
868 {
869 howmanyvarinp++;
870 }
871 }
872 eulerchar(Ip, variables-howmanyvarinp, ec);
873 id_Delete(&Ip, currRing);
874 idAddMon(I,p);
876 }
877}
878
879//tests if an ideal is Square Free, if no, returns the variable which appears at powers >1
880static poly SqFree (ideal I)
881{
882 int i,j;
883 bool flag=TRUE;
884 poly notsqrfree = NULL;
885 if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)<=1)
886 {
887 return(notsqrfree);
888 }
889 for(i=IDELEMS(I)-1;(i>=0)&&(flag);i--)
890 {
891 for(j=1;(j<=currRing->N)&&(flag);j++)
892 {
893 if(p_GetExp(I->m[i],j,currRing)>1)
894 {
895 flag=FALSE;
896 notsqrfree = p_ISet(1,currRing);
897 p_SetExp(notsqrfree,j,1,currRing);
898 }
899 }
900 }
901 if(notsqrfree != NULL)
902 {
903 p_Setm(notsqrfree,currRing);
904 }
905 return(notsqrfree);
906}
907
908//checks if a polynomial is in I
909static bool IsIn(poly p, ideal I)
910{
911 //assumes that I is ordered by degree
912 if(idIs0(I))
913 {
914 if(p==poly(0))
915 {
916 return(TRUE);
917 }
918 else
919 {
920 return(FALSE);
921 }
922 }
923 if(p==poly(0))
924 {
925 return(FALSE);
926 }
927 int i,j;
928 bool flag;
929 for(i = 0;i<IDELEMS(I);i++)
930 {
931 flag = TRUE;
932 for(j = 1;(j<=currRing->N) &&(flag);j++)
933 {
934 if(p_GetExp(p, j, currRing)<p_GetExp(I->m[i], j, currRing))
935 {
936 flag = FALSE;
937 }
938 }
939 if(flag)
940 {
941 return(TRUE);
942 }
943 }
944 return(FALSE);
945}
946
947//computes the lcm of min I, I monomial ideal
948static poly LCMmon(ideal I)
949{
950 if(idIs0(I))
951 {
952 return(NULL);
953 }
954 poly m;
955 int dummy,i,j;
956 m = p_ISet(1,currRing);
957 for(i=1;i<=currRing->N;i++)
958 {
959 dummy=0;
960 for(j=IDELEMS(I)-1;j>=0;j--)
961 {
962 if(p_GetExp(I->m[j],i,currRing) > dummy)
963 {
964 dummy = p_GetExp(I->m[j],i,currRing);
965 }
966 }
967 p_SetExp(m,i,dummy,currRing);
968 }
970 return(m);
971}
972
973//the Roune Slice Algorithm
974void rouneslice(ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int* &hilbpower)
975{
976 loop
977 {
978 (steps)++;
979 int i,j;
980 int dummy;
981 poly m;
982 ideal p;
983 //----------- PRUNING OF S ---------------
984 //S SHOULD IN THIS POINT BE ORDERED BY DEGREE
985 for(i=IDELEMS(S)-1;i>=0;i--)
986 {
987 if(IsIn(S->m[i],I))
988 {
989 p_Delete(&S->m[i],currRing);
990 prune++;
991 }
992 }
993 idSkipZeroes(S);
994 //----------------------------------------
995 for(i=IDELEMS(I)-1;i>=0;i--)
996 {
997 m = p_Head(I->m[i],currRing);
998 for(j=1;j<=currRing->N;j++)
999 {
1000 dummy = p_GetExp(m,j,currRing);
1001 if(dummy > 0)
1002 {
1003 p_SetExp(m,j,dummy-1,currRing);
1004 }
1005 }
1006 p_Setm(m, currRing);
1007 if(IsIn(m,S))
1008 {
1009 p_Delete(&I->m[i],currRing);
1010 //printf("\n Deleted, since pi(m) is in S\n");pWrite(m);
1011 }
1013 }
1014 idSkipZeroes(I);
1015 //----------- MORE PRUNING OF S ------------
1016 m = LCMmon(I);
1017 if(m != NULL)
1018 {
1019 for(i=0;i<IDELEMS(S);i++)
1020 {
1021 if(!(p_DivisibleBy(S->m[i], m, currRing)))
1022 {
1023 S->m[i] = NULL;
1024 j++;
1025 moreprune++;
1026 }
1027 else
1028 {
1029 if(pLmEqual(S->m[i],m))
1030 {
1031 S->m[i] = NULL;
1032 moreprune++;
1033 }
1034 }
1035 }
1036 idSkipZeroes(S);
1037 }
1039 /*printf("\n---------------------------\n");
1040 printf("\n I\n");idPrint(I);
1041 printf("\n S\n");idPrint(S);
1042 printf("\n q\n");pWrite(q);
1043 getchar();*/
1044
1045 if(idIs0(I))
1046 {
1047 id_Delete(&I, currRing);
1048 id_Delete(&S, currRing);
1049 break;
1050 }
1051 m = LCMmon(I);
1052 if(!p_DivisibleBy(x,m, currRing))
1053 {
1054 //printf("\nx does not divide lcm(I)");
1055 //printf("\nEmpty set");pWrite(q);
1056 id_Delete(&I, currRing);
1057 id_Delete(&S, currRing);
1058 p_Delete(&m, currRing);
1059 break;
1060 }
1061 p_Delete(&m, currRing);
1062 m = SqFree(I);
1063 if(m==NULL)
1064 {
1065 //printf("\n Corner: ");
1066 //pWrite(q);
1067 //printf("\n With the facets of the dual simplex:\n");
1068 //idPrint(I);
1069 mpz_t ec;
1070 mpz_init(ec);
1071 mpz_ptr ec_ptr = ec;
1072 eulerchar(I, currRing->N, ec_ptr);
1073 bool flag = FALSE;
1074 if(NNN==0)
1075 {
1076 hilbertcoef = (mpz_ptr)omAlloc((NNN+1)*sizeof(mpz_t));
1077 hilbpower = (int*)omAlloc((NNN+1)*sizeof(int));
1078 mpz_init_set( &hilbertcoef[NNN], ec);
1079 hilbpower[NNN] = p_Totaldegree(q,currRing);
1080 NNN++;
1081 }
1082 else
1083 {
1084 //I look if the power appears already
1085 for(i = 0;(i<NNN)&&(flag == FALSE)&&(p_Totaldegree(q,currRing)>=hilbpower[i]);i++)
1086 {
1087 if((hilbpower[i]) == (p_Totaldegree(q,currRing)))
1088 {
1089 flag = TRUE;
1090 mpz_add(&hilbertcoef[i],&hilbertcoef[i],ec_ptr);
1091 }
1092 }
1093 if(flag == FALSE)
1094 {
1095 hilbertcoef = (mpz_ptr)omRealloc(hilbertcoef, (NNN+1)*sizeof(mpz_t));
1096 hilbpower = (int*)omRealloc(hilbpower, (NNN+1)*sizeof(int));
1097 mpz_init(&hilbertcoef[NNN]);
1098 for(j = NNN; j>i; j--)
1099 {
1100 mpz_set(&hilbertcoef[j],&hilbertcoef[j-1]);
1101 hilbpower[j] = hilbpower[j-1];
1102 }
1103 mpz_set( &hilbertcoef[i], ec);
1104 hilbpower[i] = p_Totaldegree(q,currRing);
1105 NNN++;
1106 }
1107 }
1108 mpz_clear(ec);
1109 id_Delete(&I, currRing);
1110 id_Delete(&S, currRing);
1111 break;
1112 }
1113 else
1114 p_Delete(&m, currRing);
1115 m = ChooseP(I);
1116 p = idInit(1,1);
1117 p->m[0] = m;
1118 ideal Ip = idQuotMon(I,p);
1119 ideal Sp = idQuotMon(S,p);
1120 poly pq = pp_Mult_mm(q,m,currRing);
1121 rouneslice(Ip, Sp, pq, x, prune, moreprune, steps, NNN, hilbertcoef,hilbpower);
1122 idAddMon(S,p);
1123 p->m[0]=NULL;
1124 id_Delete(&p, currRing); // p->m[0] was also in S
1125 p_Delete(&pq,currRing);
1126 }
1127}
1128
1129//it computes the first hilbert series by means of Roune Slice Algorithm
1130void slicehilb(ideal I)
1131{
1132 //printf("Adi changes are here: \n");
1133 int i, NNN = 0;
1134 int steps = 0, prune = 0, moreprune = 0;
1135 mpz_ptr hilbertcoef;
1136 int *hilbpower;
1137 ideal S = idInit(1,1);
1138 poly q = p_One(currRing);
1139 ideal X = idInit(1,1);
1140 X->m[0]=p_One(currRing);
1141 for(i=1;i<=currRing->N;i++)
1142 {
1143 p_SetExp(X->m[0],i,1,currRing);
1144 }
1145 p_Setm(X->m[0],currRing);
1146 I = id_Mult(I,X,currRing);
1147 ideal Itmp = SortByDeg(I);
1148 id_Delete(&I,currRing);
1149 I = Itmp;
1150 //printf("\n-------------RouneSlice--------------\n");
1151 rouneslice(I,S,q,X->m[0],prune, moreprune, steps, NNN, hilbertcoef, hilbpower);
1152 id_Delete(&X,currRing);
1153 p_Delete(&q,currRing);
1154 //printf("\nIn total Prune got rid of %i elements\n",prune);
1155 //printf("\nIn total More Prune got rid of %i elements\n",moreprune);
1156 //printf("\nSteps of rouneslice: %i\n\n", steps);
1157 printf("\n// %8d t^0",1);
1158 for(i = 0; i<NNN; i++)
1159 {
1160 if(mpz_sgn(&hilbertcoef[i])!=0)
1161 {
1162 gmp_printf("\n// %8Zd t^%d",&hilbertcoef[i],hilbpower[i]);
1163 }
1164 }
1165 PrintLn();
1166 omFreeSize(hilbertcoef, (NNN)*sizeof(mpz_t));
1167 omFreeSize(hilbpower, (NNN)*sizeof(int));
1168 //printf("\n-------------------------------------\n");
1169}
1170
1171// -------------------------------- END OF CHANGES -------------------------------------------
1172static intvec * hSeries(ideal S, intvec *modulweight,
1173 int /*notstc*/, intvec *wdegree, ideal Q, ring tailRing)
1174{
1175// id_TestTail(S, currRing, tailRing);
1176
1177 intvec *work, *hseries1=NULL;
1178 int mc;
1179 int p0;
1180 int i, j, k, l, ii, mw;
1181 hexist = hInit(S, Q, &hNexist, tailRing);
1182 if (hNexist==0)
1183 {
1184 hseries1=new intvec(2);
1185 (*hseries1)[0]=1;
1186 (*hseries1)[1]=0;
1187 return hseries1;
1188 }
1189
1190 #if 0
1191 if (wdegree == NULL)
1192 hWeight();
1193 else
1194 hWDegree(wdegree);
1195 #else
1196 if (wdegree != NULL) hWDegree(wdegree);
1197 #endif
1198
1199 p0 = 1;
1200 hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
1201 hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
1202 hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
1203 stcmem = hCreate((currRing->N) - 1);
1204 Qpol = (int **)omAlloc(((currRing->N) + 1) * sizeof(int *));
1205 Ql = (int *)omAlloc0(((currRing->N) + 1) * sizeof(int));
1206 Q0 = (int *)omAlloc(((currRing->N) + 1) * sizeof(int));
1207 *Qpol = NULL;
1208 hLength = k = j = 0;
1209 mc = hisModule;
1210 if (mc!=0)
1211 {
1212 mw = hMinModulweight(modulweight);
1213 hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1214 }
1215 else
1216 {
1217 mw = 0;
1218 hstc = hexist;
1219 hNstc = hNexist;
1220 }
1221 loop
1222 {
1223 if (mc!=0)
1224 {
1225 hComp(hexist, hNexist, mc, hstc, &hNstc);
1226 if (modulweight != NULL)
1227 j = (*modulweight)[mc-1]-mw;
1228 }
1229 if (hNstc!=0)
1230 {
1231 hNvar = (currRing->N);
1232 for (i = hNvar; i>=0; i--)
1233 hvar[i] = i;
1234 //if (notstc) // TODO: no mon divides another
1236 hSupp(hstc, hNstc, hvar, &hNvar);
1237 if (hNvar!=0)
1238 {
1239 if ((hNvar > 2) && (hNstc > 10))
1242 memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
1243 hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
1245 Q0[hNvar] = 0;
1246 hHilbStep(hpure, hstc, hNstc, hvar, hNvar, &p0, 1);
1247 }
1248 }
1249 else
1250 {
1251 if(*Qpol!=NULL)
1252 (**Qpol)++;
1253 else
1254 {
1255 *Qpol = (int *)omAlloc(sizeof(int));
1256 hLength = *Ql = **Qpol = 1;
1257 }
1258 }
1259 if (*Qpol!=NULL)
1260 {
1261 i = hLength;
1262 while ((i > 0) && ((*Qpol)[i - 1] == 0))
1263 i--;
1264 if (i > 0)
1265 {
1266 l = i + j;
1267 if (l > k)
1268 {
1269 work = new intvec(l);
1270 for (ii=0; ii<k; ii++)
1271 (*work)[ii] = (*hseries1)[ii];
1272 if (hseries1 != NULL)
1273 delete hseries1;
1274 hseries1 = work;
1275 k = l;
1276 }
1277 while (i > 0)
1278 {
1279 (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
1280 (*Qpol)[i - 1] = 0;
1281 i--;
1282 }
1283 }
1284 }
1285 mc--;
1286 if (mc <= 0)
1287 break;
1288 }
1289 if (k==0)
1290 {
1291 hseries1=new intvec(2);
1292 (*hseries1)[0]=0;
1293 (*hseries1)[1]=0;
1294 }
1295 else
1296 {
1297 l = k+1;
1298 while ((*hseries1)[l-2]==0) l--;
1299 if (l!=k)
1300 {
1301 work = new intvec(l);
1302 for (ii=l-2; ii>=0; ii--)
1303 (*work)[ii] = (*hseries1)[ii];
1304 delete hseries1;
1305 hseries1 = work;
1306 }
1307 (*hseries1)[l-1] = mw;
1308 }
1309 for (i = 0; i <= (currRing->N); i++)
1310 {
1311 if (Ql[i]!=0)
1312 omFreeSize((ADDRESS)Qpol[i], Ql[i] * sizeof(int));
1313 }
1314 omFreeSize((ADDRESS)Q0, ((currRing->N) + 1) * sizeof(int));
1315 omFreeSize((ADDRESS)Ql, ((currRing->N) + 1) * sizeof(int));
1316 omFreeSize((ADDRESS)Qpol, ((currRing->N) + 1) * sizeof(int *));
1317 hKill(stcmem, (currRing->N) - 1);
1318 omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
1319 omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
1320 omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1322 if (hisModule!=0)
1323 omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1324 return hseries1;
1325}
1326
1327
1328intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
1329{
1330 id_TestTail(S, currRing, tailRing);
1331 if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
1332 return hSeries(S, modulweight, 0, wdegree, Q, tailRing);
1333}
1334
1335intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
1336{
1337 id_TestTail(S, currRing, tailRing);
1338 if (Q!= NULL) id_TestTail(Q, currRing, tailRing);
1339
1340 intvec *hseries1= hSeries(S, modulweight, 1, wdegree, Q, tailRing);
1341 if (errorreported) { delete hseries1; hseries1=NULL; }
1342 return hseries1;
1343}
1344
1346{
1347 intvec *work, *hseries2;
1348 int i, j, k, t, l;
1349 int s;
1350 if (hseries1 == NULL)
1351 return NULL;
1352 work = new intvec(hseries1);
1353 k = l = work->length()-1;
1354 s = 0;
1355 for (i = k-1; i >= 0; i--)
1356 s += (*work)[i];
1357 loop
1358 {
1359 if ((s != 0) || (k == 1))
1360 break;
1361 s = 0;
1362 t = (*work)[k-1];
1363 k--;
1364 for (i = k-1; i >= 0; i--)
1365 {
1366 j = (*work)[i];
1367 (*work)[i] = -t;
1368 s += t;
1369 t += j;
1370 }
1371 }
1372 hseries2 = new intvec(k+1);
1373 for (i = k-1; i >= 0; i--)
1374 (*hseries2)[i] = (*work)[i];
1375 (*hseries2)[k] = (*work)[l];
1376 delete work;
1377 return hseries2;
1378}
1379
1380void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
1381{
1382 int i, j, k;
1383 int m;
1384 *co = *mu = 0;
1385 if ((s1 == NULL) || (s2 == NULL))
1386 return;
1387 i = s1->length();
1388 j = s2->length();
1389 if (j > i)
1390 return;
1391 m = 0;
1392 for(k=j-2; k>=0; k--)
1393 m += (*s2)[k];
1394 *mu = m;
1395 *co = i - j;
1396}
1397
1398static void hPrintHilb(intvec *hseries,intvec *modul_weight)
1399{
1400 int i, j, l, k;
1401 if (hseries == NULL)
1402 return;
1403 l = hseries->length()-1;
1404 k = (*hseries)[l];
1405 if ((modul_weight!=NULL)&&(modul_weight->compare(0)!=0))
1406 {
1407 char *s=modul_weight->ivString(1,0,1);
1408 Print("module weights:%s\n",s);
1409 omFree(s);
1410 }
1411 for (i = 0; i < l; i++)
1412 {
1413 j = (*hseries)[i];
1414 if (j != 0)
1415 {
1416 Print("// %8d t^%d\n", j, i+k);
1417 }
1418 }
1419}
1420
1421/*
1422*caller
1423*/
1424void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
1425{
1426 id_TestTail(S, currRing, tailRing);
1427
1428 intvec *hseries1 = hFirstSeries(S, modulweight, Q, wdegree, tailRing);
1429 if (errorreported) return;
1430
1431 hPrintHilb(hseries1,modulweight);
1432
1433 const int l = hseries1->length()-1;
1434
1435 intvec *hseries2 = (l > 1) ? hSecondSeries(hseries1) : hseries1;
1436
1437 int co, mu;
1438 hDegreeSeries(hseries1, hseries2, &co, &mu);
1439
1440 PrintLn();
1441 hPrintHilb(hseries2,modulweight);
1442 if ((l == 1) &&(mu == 0))
1444 else
1445 scPrintDegree(co, mu);
1446 if (l>1)
1447 delete hseries1;
1448 delete hseries2;
1449}
1450
1451/***********************************************************************
1452 Computation of Hilbert series of non-commutative monomial algebras
1453************************************************************************/
1454
1455
1456static void idInsertMonomial(ideal I, poly p)
1457{
1458 /*
1459 * It adds monomial in I and if required,
1460 * enlarge the size of poly-set by 16.
1461 * It does not make copy of p.
1462 */
1463
1464 if(I == NULL)
1465 {
1466 return;
1467 }
1468
1469 int j = IDELEMS(I) - 1;
1470 while ((j >= 0) && (I->m[j] == NULL))
1471 {
1472 j--;
1473 }
1474 j++;
1475 if (j == IDELEMS(I))
1476 {
1477 pEnlargeSet(&(I->m), IDELEMS(I), 16);
1478 IDELEMS(I) +=16;
1479 }
1480 I->m[j] = p;
1481}
1482
1483static int comapreMonoIdBases(ideal J, ideal Ob)
1484{
1485 /*
1486 * Monomials of J and Ob are assumed to
1487 * be already sorted. J and Ob are
1488 * represented by the minimal generating set.
1489 */
1490 int i, s;
1491 s = 1;
1492 int JCount = IDELEMS(J);
1493 int ObCount = IDELEMS(Ob);
1494
1495 if(idIs0(J))
1496 {
1497 return(1);
1498 }
1499 if(JCount != ObCount)
1500 {
1501 return(0);
1502 }
1503
1504 for(i = 0; i < JCount; i++)
1505 {
1506 if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1507 {
1508 return(0);
1509 }
1510 }
1511 return(s);
1512}
1513
1514static int CountOnIdUptoTruncationIndex(ideal I, int tr)
1515{
1516 /*
1517 * The ideal I must be sorted in increasing total degree.
1518 * It counts the number of monomials in I up to
1519 * degree less than or equal to tr.
1520 */
1521
1522 //case when I=1;
1523 if(p_Totaldegree(I->m[0], currRing) == 0)
1524 {
1525 return(1);
1526 }
1527
1528 int count = 0;
1529 for(int i = 0; i < IDELEMS(I); i++)
1530 {
1531 if(p_Totaldegree(I->m[i], currRing) > tr)
1532 {
1533 return (count);
1534 }
1535 count = count + 1;
1536 }
1537
1538 return(count);
1539}
1540
1541static int comapreMonoIdBases_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
1542{
1543 /*
1544 * Monomials of J and Ob are assumed to
1545 * be already sorted in increasing degrees.
1546 * J and Ob are represented by the minimal
1547 * generating set. It checks if J and Ob have
1548 * same monomials up to deg <=tr.
1549 */
1550
1551 int i, s;
1552 s = 1;
1553 //when J is null
1554 //
1555 if(JCount != ObCount)
1556 {
1557 return(0);
1558 }
1559
1560 if(JCount == 0)
1561 {
1562 return(1);
1563 }
1564
1565 for(i = 0; i< JCount; i++)
1566 {
1567 if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1568 {
1569 return(0);
1570 }
1571 }
1572
1573 return(s);
1574}
1575
1576static int positionInOrbit_IG_Case(ideal I, poly w, std::vector<ideal> idorb, std::vector<poly> polist, int trInd, int trunDegHs)
1577{
1578 /*
1579 * It compares the ideal I with ideals in the set 'idorb'
1580 * up to total degree =
1581 * trInd - max(deg of w, deg of word in polist) polynomials.
1582 *
1583 * It returns 0 if I is not equal to any ideal in the
1584 * 'idorb' else returns position of the matched ideal.
1585 */
1586
1587 int ps = 0;
1588 int i, s = 0;
1589 int orbCount = idorb.size();
1590
1591 if(idIs0(I))
1592 {
1593 return(1);
1594 }
1595
1596 int degw = p_Totaldegree(w, currRing);
1597 int degp;
1598 int dtr;
1599 int dtrp;
1600
1601 dtr = trInd - degw;
1602 int IwCount;
1603
1604 IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1605
1606 if(IwCount == 0)
1607 {
1608 return(1);
1609 }
1610
1611 int ObCount;
1612
1613 bool flag2 = FALSE;
1614
1615 for(i = 1;i < orbCount; i++)
1616 {
1617 degp = p_Totaldegree(polist[i], currRing);
1618 if(degw > degp)
1619 {
1620 dtr = trInd - degw;
1621
1622 ObCount = 0;
1623 ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1624 if(ObCount == 0)
1625 {continue;}
1626 if(flag2)
1627 {
1628 IwCount = 0;
1629 IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1630 flag2 = FALSE;
1631 }
1632 }
1633 else
1634 {
1635 flag2 = TRUE;
1636 dtrp = trInd - degp;
1637 ObCount = 0;
1638 ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtrp);
1639 IwCount = 0;
1640 IwCount = CountOnIdUptoTruncationIndex(I, dtrp);
1641 }
1642
1643 s = comapreMonoIdBases_IG_Case(I, IwCount, idorb[i], ObCount);
1644
1645 if(s)
1646 {
1647 ps = i + 1;
1648 break;
1649 }
1650 }
1651 return(ps);
1652}
1653
1654static int positionInOrbit_FG_Case(ideal I, poly, std::vector<ideal> idorb, std::vector<poly>, int, int)
1655{
1656 /*
1657 * It compares the ideal I with ideals in the set 'idorb'.
1658 * I and ideals of 'idorb' are sorted.
1659 *
1660 * It returns 0 if I is not equal to any ideal of 'idorb'
1661 * else returns position of the matched ideal.
1662 */
1663 int ps = 0;
1664 int i, s = 0;
1665 int OrbCount = idorb.size();
1666
1667 if(idIs0(I))
1668 {
1669 return(1);
1670 }
1671
1672 for(i = 1; i < OrbCount; i++)
1673 {
1674 s = comapreMonoIdBases(I, idorb[i]);
1675 if(s)
1676 {
1677 ps = i + 1;
1678 break;
1679 }
1680 }
1681
1682 return(ps);
1683}
1684
1685static int positionInOrbitTruncationCase(ideal I, poly w, std::vector<ideal> idorb, std::vector<poly> polist, int , int trunDegHs)
1686{
1687 /*
1688 * It compares the ideal I with ideals in the set 'idorb'.
1689 * I and ideals in 'idorb' are sorted.
1690
1691 * returns 0 if I is not equal to any ideal of 'idorb'
1692 * else returns position of the matched ideal.
1693 */
1694
1695 int ps = 0;
1696 int i, s = 0;
1697 int OrbCount = idorb.size();
1698 int dtr=0; int IwCount, ObCount;
1699 dtr = trunDegHs - 1 - p_Totaldegree(w, currRing);
1700
1701 if(idIs0(I))
1702 {
1703 for(i = 1; i < OrbCount; i++)
1704 {
1705 if(p_Totaldegree(w, currRing) == p_Totaldegree(polist[i], currRing))
1706 {
1707 if(idIs0(idorb[i]))
1708 return(i+1);
1709 ObCount=0;
1710 ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1711 if(ObCount==0)
1712 {
1713 ps = i + 1;
1714 break;
1715 }
1716 }
1717 }
1718
1719 return(ps);
1720 }
1721
1722 IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1723
1724 if(p_Totaldegree(I->m[0], currRing)==0)
1725 {
1726 for(i = 1; i < OrbCount; i++)
1727 {
1728 if(idIs0(idorb[i]))
1729 continue;
1730 if(p_Totaldegree(idorb[i]->m[0], currRing)==0)
1731 {
1732 ps = i + 1;
1733 break;
1734 }
1735 }
1736 return(ps);
1737 }
1738
1739 for(i = 1; i < OrbCount; i++)
1740 {
1741 if(p_Totaldegree(w, currRing) == p_Totaldegree(polist[i], currRing))
1742 {
1743 if(idIs0(idorb[i]))
1744 continue;
1745 ObCount=0;
1746 ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1747 s = comapreMonoIdBases_IG_Case(I, IwCount, idorb[i], ObCount);
1748 if(s)
1749 {
1750 ps = i + 1;
1751 break;
1752 }
1753 }
1754 }
1755
1756 return(ps);
1757}
1758
1759static int monCompare( const void *m, const void *n)
1760{
1761 /* compares monomials */
1762
1763 return(p_Compare(*(poly*) m, *(poly*)n, currRing));
1764}
1765
1767{
1768 /*
1769 * sorts monomial ideal in ascending order
1770 * order must be a total degree
1771 */
1772
1773 qsort(I->m, IDELEMS(I), sizeof(poly), monCompare);
1774
1775}
1776
1777
1778
1779static ideal minimalMonomialGenSet(ideal I)
1780{
1781 /*
1782 * eliminates monomials which
1783 * can be generated by others in I
1784 */
1785 //first sort monomials of the ideal
1786
1787 idSkipZeroes(I);
1788
1790
1791 int i, k;
1792 int ICount = IDELEMS(I);
1793
1794 for(k = ICount - 1; k >=1; k--)
1795 {
1796 for(i = 0; i < k; i++)
1797 {
1798
1799 if(p_LmDivisibleBy(I->m[i], I->m[k], currRing))
1800 {
1801 pDelete(&(I->m[k]));
1802 break;
1803 }
1804 }
1805 }
1806
1807 idSkipZeroes(I);
1808 return(I);
1809}
1810
1811static poly shiftInMon(poly p, int i, int lV, const ring r)
1812{
1813 /*
1814 * shifts the varibles of monomial p in the i^th layer,
1815 * p remains unchanged,
1816 * creates new poly and returns it for the colon ideal
1817 */
1818 poly smon = p_One(r);
1819 int j, sh, cnt;
1820 cnt = r->N;
1821 sh = i*lV;
1822 int *e=(int *)omAlloc((r->N+1)*sizeof(int));
1823 int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
1824 p_GetExpV(p, e, r);
1825
1826 for(j = 1; j <= cnt; j++)
1827 {
1828 if(e[j] == 1)
1829 {
1830 s[j+sh] = e[j];
1831 }
1832 }
1833
1834 p_SetExpV(smon, s, currRing);
1835 omFree(e);
1836 omFree(s);
1837
1839 p_Setm(smon, currRing);
1840
1841 return(smon);
1842}
1843
1844static poly deleteInMon(poly w, int i, int lV, const ring r)
1845{
1846 /*
1847 * deletes the variables upto i^th layer of monomial w
1848 * w remains unchanged
1849 * creates new poly and returns it for the colon ideal
1850 */
1851
1852 poly dw = p_One(currRing);
1853 int *e = (int *)omAlloc((r->N+1)*sizeof(int));
1854 int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
1855 p_GetExpV(w, e, r);
1856 int j, cnt;
1857 cnt = i*lV;
1858 /*
1859 for(j=1;j<=cnt;j++)
1860 {
1861 e[j]=0;
1862 }*/
1863 for(j = (cnt+1); j < (r->N+1); j++)
1864 {
1865 s[j] = e[j];
1866 }
1867
1868 p_SetExpV(dw, s, currRing);//new exponents
1869 omFree(e);
1870 omFree(s);
1871
1873 p_Setm(dw, currRing);
1874
1875 return(dw);
1876}
1877
1878static void TwordMap(poly p, poly w, int lV, int d, ideal Jwi, bool &flag)
1879{
1880 /*
1881 * computes T_w(p) in a new poly object and places it
1882 * in Jwi which stores elements of colon ideal of I,
1883 * p and w remain unchanged,
1884 * the new polys for Jwi are constructed by sub-routines
1885 * deleteInMon, shiftInMon, p_MDivide,
1886 * places the result in Jwi and deletes the new polys
1887 * coming in dw, smon, qmon
1888 */
1889 int i;
1890 poly smon, dw;
1891 poly qmonp = NULL;
1892 bool del;
1893
1894 for(i = 0;i <= d - 1; i++)
1895 {
1896 dw = deleteInMon(w, i, lV, currRing);
1897 smon = shiftInMon(p, i, lV, currRing);
1898 del = TRUE;
1899
1900 if(pLmDivisibleBy(smon, w))
1901 {
1902 flag = TRUE;
1903 del = FALSE;
1904
1905 pDelete(&dw);
1906 pDelete(&smon);
1907
1908 //delete all monomials of Jwi
1909 //and make Jwi =1
1910
1911 for(int j = 0;j < IDELEMS(Jwi); j++)
1912 {
1913 pDelete(&Jwi->m[j]);
1914 }
1915
1917 break;
1918 }
1919
1920 if(pLmDivisibleBy(dw, smon))
1921 {
1922 del = FALSE;
1923 qmonp = p_MDivide(smon, dw, currRing);
1924 idInsertMonomial(Jwi, shiftInMon(qmonp, -d, lV, currRing));
1925 pLmFree(&qmonp);
1926 pDelete(&dw);
1927 pDelete(&smon);
1928 }
1929 //in case both if are false, delete dw and smon
1930 if(del)
1931 {
1932 pDelete(&dw);
1933 pDelete(&smon);
1934 }
1935 }
1936
1937}
1938
1939static ideal colonIdeal(ideal S, poly w, int lV, ideal Jwi, int trunDegHs)
1940{
1941 /*
1942 * It computes the right colon ideal of a two-sided ideal S
1943 * w.r.t. word w and save it in a new object Jwi.
1944 * It keeps S and w unchanged.
1945 */
1946
1947 if(idIs0(S))
1948 {
1949 return(S);
1950 }
1951
1952 int i, d;
1953 d = p_Totaldegree(w, currRing);
1954 if(trunDegHs !=0 && d >= trunDegHs)
1955 {
1957 return(Jwi);
1958 }
1959 bool flag = FALSE;
1960 int SCount = IDELEMS(S);
1961 for(i = 0; i < SCount; i++)
1962 {
1963 TwordMap(S->m[i], w, lV, d, Jwi, flag);
1964 if(flag)
1965 {
1966 break;
1967 }
1968 }
1969
1970 Jwi = minimalMonomialGenSet(Jwi);
1971 return(Jwi);
1972}
1973
1974void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp, int trunDegHs)
1975{
1976
1977 /* new story:
1978 no lV is needed, i.e. it is to be determined
1979 the rest is extracted from the interface input list in extra.cc and makes the input of this proc
1980 called from extra.cc
1981 */
1982
1983 /*
1984 * This is based on iterative right colon operations on a
1985 * two-sided monomial ideal of the free associative algebra.
1986 * The algorithm terminates for those monomial ideals
1987 * whose monomials define "regular formal languages",
1988 * that is, all monomials of the input ideal can be obtained
1989 * from finite languages by applying finite number of
1990 * rational operations.
1991 */
1992
1993 int trInd;
1994 S = minimalMonomialGenSet(S);
1995 if( !idIs0(S) && p_Totaldegree(S->m[0], currRing)==0)
1996 {
1997 PrintS("Hilbert Series:\n 0\n");
1998 return;
1999 }
2000 int (*POS)(ideal, poly, std::vector<ideal>, std::vector<poly>, int, int);
2001 if(trunDegHs != 0)
2002 {
2003 Print("\nTruncation degree = %d\n",trunDegHs);
2005 }
2006 else
2007 {
2008 if(IG_CASE)
2009 {
2010 if(idIs0(S))
2011 {
2012 WerrorS("wrong input: it is not an infinitely gen. case");
2013 return;
2014 }
2015 trInd = p_Totaldegree(S->m[IDELEMS(S)-1], currRing);
2017 }
2018 else
2020 }
2021 std::vector<ideal > idorb;
2022 std::vector< poly > polist;
2023
2024 ideal orb_init = idInit(1, 1);
2025 idorb.push_back(orb_init);
2026
2027 polist.push_back( p_One(currRing));
2028
2029 std::vector< std::vector<int> > posMat;
2030 std::vector<int> posRow(lV,0);
2031 std::vector<int> C;
2032
2033 int ds, is, ps;
2034 unsigned long lpcnt = 0;
2035
2036 poly w, wi;
2037 ideal Jwi;
2038
2039 while(lpcnt < idorb.size())
2040 {
2041 w = NULL;
2042 w = polist[lpcnt];
2043 if(lpcnt >= 1 && idIs0(idorb[lpcnt]) == FALSE)
2044 {
2045 if(p_Totaldegree(idorb[lpcnt]->m[0], currRing) != 0)
2046 {
2047 C.push_back(1);
2048 }
2049 else
2050 C.push_back(0);
2051 }
2052 else
2053 {
2054 C.push_back(1);
2055 }
2056
2057 ds = p_Totaldegree(w, currRing);
2058 lpcnt++;
2059
2060 for(is = 1; is <= lV; is++)
2061 {
2062 wi = NULL;
2063 //make new copy 'wi' of word w=polist[lpcnt]
2064 //and update it (for the colon operation).
2065 //if corresponding to wi, right colon operation gives
2066 //a new (right colon) ideal of S,
2067 //keep 'wi' in the polist else delete it
2068
2069 wi = pCopy(w);
2070 p_SetExp(wi, (ds*lV)+is, 1, currRing);
2071 p_Setm(wi, currRing);
2072 Jwi = NULL;
2073 //Jwi stores (right) colon ideal of S w.r.t. word
2074 //wi if colon operation gives a new ideal place it
2075 //in the vector of ideals 'idorb'
2076 //otherwise delete it
2077
2078 Jwi = idInit(1,1);
2079
2080 Jwi = colonIdeal(S, wi, lV, Jwi, trunDegHs);
2081 ps = (*POS)(Jwi, wi, idorb, polist, trInd, trunDegHs);
2082
2083 if(ps == 0) // finds a new ideal
2084 {
2085 posRow[is-1] = idorb.size();
2086
2087 idorb.push_back(Jwi);
2088 polist.push_back(wi);
2089 }
2090 else // ideal is already there in the set
2091 {
2092 posRow[is-1]=ps-1;
2093 idDelete(&Jwi);
2094 pDelete(&wi);
2095 }
2096 }
2097 posMat.push_back(posRow);
2098 posRow.resize(lV,0);
2099 }
2100 int lO = C.size();//size of the orbit
2101 PrintLn();
2102 Print("maximal length of words = %ld\n", p_Totaldegree(polist[lO-1], currRing));
2103 Print("\nlength of the Orbit = %d", lO);
2104 PrintLn();
2105
2106 if(odp)
2107 {
2108 Print("words description of the Orbit: \n");
2109 for(is = 0; is < lO; is++)
2110 {
2111 pWrite0(polist[is]);
2112 PrintS(" ");
2113 }
2114 PrintLn();
2115 PrintS("\nmaximal degree, #(sum_j R(w,w_j))");
2116 PrintLn();
2117 for(is = 0; is < lO; is++)
2118 {
2119 if(idIs0(idorb[is]))
2120 {
2121 PrintS("NULL\n");
2122 }
2123 else
2124 {
2125 Print("%ld, %d \n",p_Totaldegree(idorb[is]->m[IDELEMS(idorb[is])-1], currRing),IDELEMS(idorb[is]));
2126 }
2127 }
2128 }
2129
2130 for(is = idorb.size()-1; is >= 0; is--)
2131 {
2132 idDelete(&idorb[is]);
2133 }
2134 for(is = polist.size()-1; is >= 0; is--)
2135 {
2136 pDelete(&polist[is]);
2137 }
2138
2139 idorb.resize(0);
2140 polist.resize(0);
2141
2142 int adjMatrix[lO][lO];
2143 memset(adjMatrix, 0, lO*lO*sizeof(int));
2144 int rowCount, colCount;
2145 int tm = 0;
2146 if(!mgrad)
2147 {
2148 for(rowCount = 0; rowCount < lO; rowCount++)
2149 {
2150 for(colCount = 0; colCount < lV; colCount++)
2151 {
2152 tm = posMat[rowCount][colCount];
2153 adjMatrix[rowCount][tm] = adjMatrix[rowCount][tm] + 1;
2154 }
2155 }
2156 }
2157
2158 ring r = currRing;
2159 int npar;
2160 char** tt;
2162 if(!mgrad)
2163 {
2164 tt=(char**)omAlloc(sizeof(char*));
2165 tt[0] = omStrDup("t");
2166 npar = 1;
2167 }
2168 else
2169 {
2170 tt=(char**)omalloc(lV*sizeof(char*));
2171 for(is = 0; is < lV; is++)
2172 {
2173 tt[is] = (char*)omAlloc(7*sizeof(char)); //if required enlarge it later
2174 sprintf (tt[is], "t%d", is+1);
2175 }
2176 npar = lV;
2177 }
2178
2179 p.r = rDefault(0, npar, tt);
2181 char** xx = (char**)omAlloc(sizeof(char*));
2182 xx[0] = omStrDup("x");
2183 ring R = rDefault(cf, 1, xx);
2184 rChangeCurrRing(R);//rWrite(R);
2185 /*
2186 * matrix corresponding to the orbit of the ideal
2187 */
2188 matrix mR = mpNew(lO, lO);
2189 matrix cMat = mpNew(lO,1);
2190 poly rc;
2191
2192 if(!mgrad)
2193 {
2194 for(rowCount = 0; rowCount < lO; rowCount++)
2195 {
2196 for(colCount = 0; colCount < lO; colCount++)
2197 {
2198 if(adjMatrix[rowCount][colCount] != 0)
2199 {
2200 MATELEM(mR, rowCount + 1, colCount + 1) = p_ISet(adjMatrix[rowCount][colCount], R);
2201 p_SetCoeff(MATELEM(mR, rowCount + 1, colCount + 1), n_Mult(pGetCoeff(mR->m[lO*rowCount+colCount]),n_Param(1, R->cf), R->cf), R);
2202 }
2203 }
2204 }
2205 }
2206 else
2207 {
2208 for(rowCount = 0; rowCount < lO; rowCount++)
2209 {
2210 for(colCount = 0; colCount < lV; colCount++)
2211 {
2212 rc=NULL;
2213 rc=p_One(R);
2214 p_SetCoeff(rc, n_Mult(pGetCoeff(rc), n_Param(colCount+1, R->cf),R->cf), R);
2215 MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1)=p_Add_q(rc,MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1), R);
2216 }
2217 }
2218 }
2219
2220 for(rowCount = 0; rowCount < lO; rowCount++)
2221 {
2222 if(C[rowCount] != 0)
2223 {
2224 MATELEM(cMat, rowCount + 1, 1) = p_ISet(C[rowCount], R);
2225 }
2226 }
2227
2228 matrix u;
2229 unitMatrix(lO, u); //unit matrix
2230 matrix gMat = mp_Sub(u, mR, R);
2231
2232 char* s;
2233
2234 if(odp)
2235 {
2236 PrintS("\nlinear system:\n");
2237 if(!mgrad)
2238 {
2239 for(rowCount = 0; rowCount < lO; rowCount++)
2240 {
2241 Print("H(%d) = ", rowCount+1);
2242 for(colCount = 0; colCount < lV; colCount++)
2243 {
2244 StringSetS(""); nWrite(n_Param(1, R->cf));
2245 s = StringEndS(); PrintS(s);
2246 Print("*"); omFree(s);
2247 Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2248 }
2249 Print(" %d\n", C[rowCount] );
2250 }
2251 PrintS("where H(1) represents the series corresp. to input ideal\n");
2252 PrintS("and i^th summand in the rhs of an eqn. is according\n");
2253 PrintS("to the right colon map corresp. to the i^th variable\n");
2254 }
2255 else
2256 {
2257 for(rowCount = 0; rowCount < lO; rowCount++)
2258 {
2259 Print("H(%d) = ", rowCount+1);
2260 for(colCount = 0; colCount < lV; colCount++)
2261 {
2262 StringSetS(""); nWrite(n_Param(colCount+1, R->cf));
2263 s = StringEndS(); PrintS(s);
2264 Print("*");omFree(s);
2265 Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2266 }
2267 Print(" %d\n", C[rowCount] );
2268 }
2269 PrintS("where H(1) represents the series corresp. to input ideal\n");
2270 }
2271 }
2272 PrintLn();
2273 posMat.resize(0);
2274 C.resize(0);
2275 matrix pMat;
2276 matrix lMat;
2277 matrix uMat;
2278 matrix H_serVec = mpNew(lO, 1);
2279 matrix Hnot;
2280
2281 //std::clock_t start;
2282 //start = std::clock();
2283
2284 luDecomp(gMat, pMat, lMat, uMat, R);
2285 luSolveViaLUDecomp(pMat, lMat, uMat, cMat, H_serVec, Hnot);
2286
2287 //to print system solving time
2288 //if(odp){
2289 //std::cout<<"solving time of the system = "<<(std::clock()-start)/(double)(CLOCKS_PER_SEC / 1000)<<" ms"<<std::endl;}
2290
2291 mp_Delete(&mR, R);
2292 mp_Delete(&u, R);
2293 mp_Delete(&pMat, R);
2294 mp_Delete(&lMat, R);
2295 mp_Delete(&uMat, R);
2296 mp_Delete(&cMat, R);
2297 mp_Delete(&gMat, R);
2298 mp_Delete(&Hnot, R);
2299 //print the Hilbert series and length of the Orbit
2300 PrintLn();
2301 Print("Hilbert series:");
2302 PrintLn();
2303 pWrite(H_serVec->m[0]);
2304 if(!mgrad)
2305 {
2306 omFree(tt[0]);
2307 }
2308 else
2309 {
2310 for(is = lV-1; is >= 0; is--)
2311
2312 omFree( tt[is]);
2313 }
2314 omFree(tt);
2315 omFree(xx[0]);
2316 omFree(xx);
2317 rChangeCurrRing(r);
2318 rKill(R);
2319}
2320
2321ideal RightColonOperation(ideal S, poly w, int lV)
2322{
2323 /*
2324 * This returns right colon ideal of a monomial two-sided ideal of
2325 * the free associative algebra with respect to a monomial 'w'
2326 * (S:_R w).
2327 */
2328 S = minimalMonomialGenSet(S);
2329 ideal Iw = idInit(1,1);
2330 Iw = colonIdeal(S, w, lV, Iw, 0);
2331 return (Iw);
2332}
2333
#define NULL
Definition: auxiliary.h:104
long int64
Definition: auxiliary.h:68
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
void * ADDRESS
Definition: auxiliary.h:119
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 k
Definition: cfEzgcd.cc:99
Variable x
Definition: cfModGcd.cc:4084
int p
Definition: cfModGcd.cc:4080
CanonicalForm cf
Definition: cfModGcd.cc:4085
CanonicalForm b
Definition: cfModGcd.cc:4105
void mu(int **points, int sizePoints)
Definition: intvec.h:23
int length() const
Definition: intvec.h:94
int compare(const intvec *o) const
Definition: intvec.cc:206
char * ivString(int not_mat=1, int spaces=0, int dim=2) const
Definition: intvec.cc:58
int rows() const
Definition: intvec.h:96
poly * m
Definition: matpol.h:18
Coefficient rings, fields and other domains suitable for Singular polynomials.
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of 'a' and 'b', i.e., a*b
Definition: coeffs.h:637
static FORCE_INLINE number n_Param(const int iParameter, const coeffs r)
return the (iParameter^th) parameter as a NEW number NOTE: parameter numbering: 1....
Definition: coeffs.h:807
@ n_transExt
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition: coeffs.h:39
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:358
#define Print
Definition: emacs.cc:80
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
CanonicalForm res
Definition: facAbsFact.cc:60
const CanonicalForm & w
Definition: facAbsFact.cc:51
int j
Definition: facHensel.cc:110
void FACTORY_PUBLIC prune(Variable &alpha)
Definition: variable.cc:261
VAR short errorreported
Definition: feFopen.cc:23
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define STATIC_VAR
Definition: globaldefs.h:7
void scPrintDegree(int co, int mu)
Definition: hdegree.cc:881
static void idInsertMonomial(ideal I, poly p)
Definition: hilb.cc:1456
static int comapreMonoIdBases_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
Definition: hilb.cc:1541
STATIC_VAR int * Q0
Definition: hilb.cc:45
intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1328
static poly SqFree(ideal I)
Definition: hilb.cc:880
static void idAddMon(ideal I, ideal p)
Definition: hilb.cc:471
static int comapreMonoIdBases(ideal J, ideal Ob)
Definition: hilb.cc:1483
static void TwordMap(poly p, poly w, int lV, int d, ideal Jwi, bool &flag)
Definition: hilb.cc:1878
static poly ChooseP(ideal I)
Definition: hilb.cc:764
static poly deleteInMon(poly w, int i, int lV, const ring r)
Definition: hilb.cc:1844
intvec * hSecondSeries(intvec *hseries1)
Definition: hilb.cc:1345
STATIC_VAR int hLength
Definition: hilb.cc:46
static int CountOnIdUptoTruncationIndex(ideal I, int tr)
Definition: hilb.cc:1514
static poly ChoosePJL(ideal I)
Definition: hilb.cc:705
static int monCompare(const void *m, const void *n)
Definition: hilb.cc:1759
static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hilb.cc:63
static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
Definition: hilb.cc:137
static void hPrintHilb(intvec *hseries, intvec *modul_weight)
Definition: hilb.cc:1398
static int positionInOrbitTruncationCase(ideal I, poly w, std::vector< ideal > idorb, std::vector< poly > polist, int, int trunDegHs)
Definition: hilb.cc:1685
static poly LCMmon(ideal I)
Definition: hilb.cc:948
void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp, int trunDegHs)
Definition: hilb.cc:1974
void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1424
static ideal colonIdeal(ideal S, poly w, int lV, ideal Jwi, int trunDegHs)
Definition: hilb.cc:1939
static int hMinModulweight(intvec *modulweight)
Definition: hilb.cc:49
static poly shiftInMon(poly p, int i, int lV, const ring r)
Definition: hilb.cc:1811
intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1335
static intvec * hSeries(ideal S, intvec *modulweight, int, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1172
static poly ChoosePVar(ideal I)
Definition: hilb.cc:479
static int positionInOrbit_FG_Case(ideal I, poly, std::vector< ideal > idorb, std::vector< poly >, int, int)
Definition: hilb.cc:1654
static ideal SortByDeg(ideal I)
Definition: hilb.cc:388
static int * hAddHilb(int Nv, int x, int *pol, int *lp)
Definition: hilb.cc:104
static bool IsIn(poly p, ideal I)
Definition: hilb.cc:909
static void eulerchar(ideal I, int variables, mpz_ptr ec)
Definition: hilb.cc:837
ideal RightColonOperation(ideal S, poly w, int lV)
Definition: hilb.cc:2321
static void hWDegree(intvec *wdegree)
Definition: hilb.cc:245
void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
Definition: hilb.cc:1380
static poly SearchP(ideal I)
searches for a monomial of degree d>=2 and divides it by a variable (result monomial of deg d-1)
Definition: hilb.cc:780
static ideal minimalMonomialGenSet(ideal I)
Definition: hilb.cc:1779
STATIC_VAR int * Ql
Definition: hilb.cc:45
void sortMonoIdeal_pCompare(ideal I)
Definition: hilb.cc:1766
ideal idQuotMon(ideal Iorig, ideal p)
Definition: hilb.cc:409
static int positionInOrbit_IG_Case(ideal I, poly w, std::vector< ideal > idorb, std::vector< poly > polist, int trInd, int trunDegHs)
Definition: hilb.cc:1576
static void SortByDeg_p(ideal I, poly p)
!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
Definition: hilb.cc:288
void slicehilb(ideal I)
Definition: hilb.cc:1130
static bool JustVar(ideal I)
Definition: hilb.cc:806
void rouneslice(ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int *&hilbpower)
Definition: hilb.cc:974
STATIC_VAR int ** Qpol
Definition: hilb.cc:44
static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var, int Nvar, int *pol, int Lpol)
Definition: hilb.cc:177
monf hCreate(int Nvar)
Definition: hutil.cc:999
void hComp(scfmon exist, int Nexist, int ak, scfmon stc, int *Nstc)
Definition: hutil.cc:157
scfmon hInit(ideal S, ideal Q, int *Nexist, ring tailRing)
Definition: hutil.cc:31
void hLex2S(scfmon rad, int e1, int a2, int e2, varset var, int Nvar, scfmon w)
Definition: hutil.cc:815
VAR scfmon hstc
Definition: hutil.cc:16
VAR varset hvar
Definition: hutil.cc:18
void hKill(monf xmem, int Nvar)
Definition: hutil.cc:1013
VAR int hNexist
Definition: hutil.cc:19
void hElimS(scfmon stc, int *e1, int a2, int e2, varset var, int Nvar)
Definition: hutil.cc:675
void hLexS(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:509
void hDelete(scfmon ev, int ev_length)
Definition: hutil.cc:143
VAR monf stcmem
Definition: hutil.cc:21
scfmon hGetmem(int lm, scfmon old, monp monmem)
Definition: hutil.cc:1026
void hPure(scfmon stc, int a, int *Nstc, varset var, int Nvar, scmon pure, int *Npure)
Definition: hutil.cc:624
VAR scfmon hwork
Definition: hutil.cc:16
void hSupp(scfmon stc, int Nstc, varset var, int *Nvar)
Definition: hutil.cc:177
VAR scmon hpure
Definition: hutil.cc:17
VAR int hisModule
Definition: hutil.cc:20
void hStepS(scfmon stc, int Nstc, varset var, int Nvar, int *a, int *x)
Definition: hutil.cc:952
void hStaircase(scfmon stc, int *Nstc, varset var, int Nvar)
Definition: hutil.cc:316
void hOrdSupp(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:205
VAR int hNpure
Definition: hutil.cc:19
VAR scfmon hexist
Definition: hutil.cc:16
scmon hGetpure(scmon p)
Definition: hutil.cc:1055
VAR int hNstc
Definition: hutil.cc:19
VAR int hNvar
Definition: hutil.cc:19
scmon * scfmon
Definition: hutil.h:15
int * varset
Definition: hutil.h:16
int * scmon
Definition: hutil.h:14
#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
ideal id_Copy(ideal h1, const ring r)
copy an ideal
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
#define idPrint(id)
Definition: ideals.h:46
STATIC_VAR int variables
void rKill(ring r)
Definition: ipshell.cc:6255
STATIC_VAR jList * Q
Definition: janet.cc:30
bool unitMatrix(const int n, matrix &unitMat, const ring R)
Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring R)
LU-decomposition of a given (m x n)-matrix.
bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, const matrix bVec, matrix &xVec, matrix &H)
Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decompositi...
void mp_Delete(matrix *a, const ring r)
Definition: matpol.cc:880
matrix mp_Sub(matrix a, matrix b, const ring R)
Definition: matpol.cc:196
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define assume(x)
Definition: mod2.h:387
#define p_GetComp(p, r)
Definition: monomials.h:64
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
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:357
const int MAX_INT_VAL
Definition: mylimits.h:12
The main handler for Singular numbers which are suitable for Singular polynomials.
#define nWrite(n)
Definition: numbers.h:29
#define omStrDup(s)
Definition: omAllocDecl.h:263
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omalloc(size)
Definition: omAllocDecl.h:228
#define omRealloc(addr, size)
Definition: omAllocDecl.h:225
#define omFree(addr)
Definition: omAllocDecl.h:261
#define omAlloc0(size)
Definition: omAllocDecl.h:211
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1292
poly p_MDivide(poly a, poly b, const ring r)
Definition: p_polys.cc:1479
int p_Compare(const poly a, const poly b, const ring R)
Definition: p_polys.cc:4934
poly p_One(const ring r)
Definition: p_polys.cc:1308
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3766
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:896
static poly p_Head(poly p, const ring r)
copy the i(leading) term of p
Definition: p_polys.h:826
#define p_LmEqual(p1, p2, r)
Definition: p_polys.h:1691
static void p_SetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1504
static poly pp_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:991
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent @Note: VarOffset encodes the position in p->exp
Definition: p_polys.h:488
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:247
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:233
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:412
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent @Note: the integer VarOffset encodes:
Definition: p_polys.h:469
static BOOLEAN p_IsOne(const poly p, const ring R)
either poly(1) or gen(k)?!
Definition: p_polys.h:1979
static BOOLEAN p_LmDivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1863
static BOOLEAN p_DivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1872
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:861
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1480
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:812
static long p_Totaldegree(poly p, const ring r)
Definition: p_polys.h:1467
void rChangeCurrRing(ring r)
Definition: polys.cc:15
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define pDelete(p_ptr)
Definition: polys.h:186
#define pLmEqual(p1, p2)
Definition: polys.h:111
void pWrite0(poly p)
Definition: polys.h:309
#define pLmDivisibleBy(a, b)
like pDivisibleBy, except that it is assumed that a!=NULL, b!=NULL
Definition: polys.h:140
static void pLmFree(poly p)
frees the space of the monomial m, assumes m != NULL coef is not freed, m is not advanced
Definition: polys.h:70
void pWrite(poly p)
Definition: polys.h:308
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
void StringSetS(const char *st)
Definition: reporter.cc:128
void PrintS(const char *s)
Definition: reporter.cc:284
char * StringEndS()
Definition: reporter.cc:151
void PrintLn()
Definition: reporter.cc:310
ring rDefault(const coeffs cf, int N, char **n, int ord_size, rRingOrder_t *ord, int *block0, int *block1, int **wvhdl, unsigned long bitmask)
Definition: ring.cc:102
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:594
int status int void size_t count
Definition: si_signals.h:59
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
ideal id_Head(ideal h, const ring r)
returns the ideals of initial terms
ideal id_Mult(ideal h1, ideal h2, const ring R)
h1 * h2 one h_i must be an ideal (with at least one column) the other h_i may be a module (with no co...
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
#define IDELEMS(i)
Definition: simpleideals.h:23
#define id_TestTail(A, lR, tR)
Definition: simpleideals.h:77
#define R
Definition: sirandom.c:27
#define loop
Definition: structs.h:80
struct for passing initialization parameters to naInitChar
Definition: transext.h:88