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