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