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("internal 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;
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 static int isMonoIdBasesSame(ideal J, ideal Ob)
1431 {
1432  /*
1433  * polynomials of J and Ob are assumed to
1434  * be already sorted. J and Ob are
1435  * represented by the minimal generating set
1436  */
1437  int i, s;
1438  s = 1;
1439  int JCount = IDELEMS(J);
1440  int ObCount = IDELEMS(Ob);
1441 
1442  if(idIs0(J))
1443  {
1444  return(1);
1445  }
1446  if(JCount != ObCount)
1447  {
1448  return(0);
1449  }
1450 
1451  for(i = 0; i < JCount; i++)
1452  {
1453  if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1454  {
1455  return(0);
1456  }
1457  }
1458  return(s);
1459 }
1460 
1461 static int CountOnIdUptoTruncationIndex(ideal I, int tr)
1462 {
1463  /*
1464  * I must be sorted in ascending order
1465  * counts the number of polys in I upto
1466  * degree less or equal to tr
1467  */
1468 
1469  //case when I=1;
1470  if(p_Totaldegree(I->m[0], currRing) == 0)
1471  {
1472  return(1);
1473  }
1474 
1475  int count = 0;
1476  for(int i = 0; i < IDELEMS(I); i++)
1477  {
1478  if(p_Totaldegree(I->m[i], currRing) > tr)
1479  {
1480  return (count);
1481  }
1482  count = count + 1;
1483  }
1484 
1485  return(count);
1486 }
1487 
1488 static int isMonoIdBasesSame_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
1489 {
1490  /*
1491  * polynomials of J and obc are assumed to
1492  * be already sorted. J and Ob are
1493  * represented by the minimal generating set.
1494  * checks if J and Ob are same in polys upto deg <=tr
1495  */
1496 
1497  int i, s;
1498  s = 1;
1499  //when J is null
1500  if(JCount == 0)
1501  {
1502  return(1);
1503  }
1504 
1505  if(JCount != ObCount)
1506  {
1507  return(0);
1508  }
1509 
1510  for(i = 0; i< JCount; i++)
1511  {
1512  if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1513  {
1514  return(0);
1515  }
1516  }
1517 
1518  return(s);
1519 }
1520 
1521 static int positionInOrbit_IG_Case(ideal I, poly w, std::vector<ideal> idorb, std::vector<poly> polist, int trInd)
1522 {
1523  /*
1524  * compares the ideal I with ideals in the Orbit 'idorb'
1525  * upto degree trInd - max(deg of w, deg of word in polist) polynomials;
1526  * I and ideals in the Orbit are sorted,
1527  * Orbit is ordered,
1528  *
1529  * returns 0 if I is not equal to any of the ideals
1530  * in the Orbit else returns position of the matched ideal
1531  */
1532 
1533  int ps = 0;
1534  int i, s = 0;
1535  int orbCount = idorb.size();
1536 
1537  if(idIs0(I))
1538  {
1539  return(1);
1540  }
1541 
1542  int degw = p_Totaldegree(w, currRing);
1543  int degp;
1544  int dtr;
1545  int dtrp;
1546 
1547  dtr = trInd - degw;
1548  int IwCount;
1549 
1550  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1551 
1552  if(IwCount == 0)
1553  {
1554  return(1);
1555  }
1556 
1557  int ObCount;
1558 
1559  bool flag2 = FALSE;
1560 
1561  for(i = 1;i < orbCount; i++)
1562  {
1563  degp = p_Totaldegree(polist[i], currRing);
1564  if(degw > degp)
1565  {
1566  dtr = trInd - degw;
1567 
1568  ObCount = 0;
1569  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1570  if(ObCount == 0)
1571  {continue;}
1572  if(flag2)
1573  {
1574  IwCount = 0;
1575  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1576  flag2 = FALSE;
1577  }
1578  }
1579  else
1580  {
1581  flag2 = TRUE;
1582  dtrp = trInd - degp;
1583  ObCount = 0;
1584  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtrp);
1585  IwCount = 0;
1586  IwCount = CountOnIdUptoTruncationIndex(I, dtrp);
1587  }
1588 
1589  s = isMonoIdBasesSame_IG_Case(I, IwCount, idorb[i], ObCount);
1590 
1591  if(s)
1592  {
1593  ps = i + 1;
1594  break;
1595  }
1596  }
1597  return(ps);
1598 }
1599 
1600 static int positionInOrbit_FG_Case(ideal I, poly, std::vector<ideal> idorb, std::vector<poly>, int)
1601 {
1602  /*
1603  * compares the ideal I with ideals in the Orbit 'idorb'
1604  * I and ideals in the Orbit are sorted,
1605  * Orbit is ordered,
1606  *
1607  * returns 0 if I is not equal to any of the ideals
1608  * in the Orbit else returns position of the matched ideal
1609  */
1610  int ps = 0;
1611  int i, s = 0;
1612  int OrbCount = idorb.size();
1613 
1614  if(idIs0(I))
1615  {
1616  return(1);
1617  }
1618 
1619  for(i = 1; i < OrbCount; i++)
1620  {
1621  s = isMonoIdBasesSame(I, idorb[i]);
1622  if(s)
1623  {
1624  ps = i + 1;
1625  break;
1626  }
1627  }
1628 
1629  return(ps);
1630 }
1631 
1632 static int monCompare( const void *m, const void *n)
1633 {
1634  /* compares monomials */
1635 
1636  return(p_Compare(*(poly*) m, *(poly*)n, currRing));
1637 }
1638 
1640 {
1641  /*
1642  * sorts the monomial ideal in ascending order
1643  * order must be a total degree
1644  */
1645 
1646  qsort(I->m, IDELEMS(I), sizeof(poly), monCompare);
1647 
1648 }
1649 
1650 
1651 
1652 static ideal minimalMonomialsGenSet(ideal I)
1653 {
1654  /*
1655  * eliminates monomials which
1656  * can be generated by others in I
1657  */
1658  //first sort monomials of the ideal
1659 
1660  idSkipZeroes(I);
1661 
1663 
1664  int i, k;
1665  int ICount = IDELEMS(I);
1666 
1667  for(k = ICount - 1; k >=1; k--)
1668  {
1669  for(i = 0; i < k; i++)
1670  {
1671 
1672  if(p_LmDivisibleBy(I->m[i], I->m[k], currRing))
1673  {
1674  pDelete(&(I->m[k]));
1675  break;
1676  }
1677  }
1678  }
1679 
1680  idSkipZeroes(I);
1681  return(I);
1682 }
1683 
1684 static poly shiftInMon(poly p, int i, int lV, const ring r)
1685 {
1686  /*
1687  * shifts the varibles of monomial p in the i^th layer,
1688  * p remains unchanged,
1689  * creates new poly and returns it for the colon ideal
1690  */
1691  poly smon = p_One(r);
1692  int j, sh, cnt;
1693  cnt = r->N;
1694  sh = i*lV;
1695  int *e=(int *)omAlloc((r->N+1)*sizeof(int));
1696  int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
1697  p_GetExpV(p, e, r);
1698 
1699  for(j = 1; j <= cnt; j++)
1700  {
1701  if(e[j] == 1)
1702  {
1703  s[j+sh] = e[j];
1704  }
1705  }
1706 
1707  p_SetExpV(smon, s, currRing);
1708  omFree(e);
1709  omFree(s);
1710 
1711  p_SetComp(smon, p_GetComp(p, currRing), currRing);
1712  p_Setm(smon, currRing);
1713 
1714  return(smon);
1715 }
1716 
1717 static poly deleteInMon(poly w, int i, int lV, const ring r)
1718 {
1719  /*
1720  * deletes the variables upto i^th layer of monomial w
1721  * w remains unchanged
1722  * creates new poly and returns it for the colon ideal
1723  */
1724 
1725  poly dw = p_One(currRing);
1726  int *e = (int *)omAlloc((r->N+1)*sizeof(int));
1727  int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
1728  p_GetExpV(w, e, r);
1729  int j, cnt;
1730  cnt = i*lV;
1731  /*
1732  for(j=1;j<=cnt;j++)
1733  {
1734  e[j]=0;
1735  }*/
1736  for(j = (cnt+1); j < (r->N+1); j++)
1737  {
1738  s[j] = e[j];
1739  }
1740 
1741  p_SetExpV(dw, s, currRing);//new exponents
1742  omFree(e);
1743  omFree(s);
1744 
1746  p_Setm(dw, currRing);
1747 
1748  return(dw);
1749 }
1750 
1751 static void TwordMap(poly p, poly w, int lV, int d, ideal Jwi, bool &flag)
1752 {
1753  /*
1754  * computes T_w(p) in a new poly object and places it
1755  * in a colon ideal Jwi of I
1756  * p and w remain unchanged
1757  * the new polys for Jwi are constructed by sub-routines
1758  * deleteInMon, shiftInMon, p_Divide,
1759  * places the result in Jwi and deletes the new polys
1760  * coming in dw, smon, qmon
1761  */
1762  int i;
1763  poly smon, dw;
1764  poly qmonp = NULL;
1765  bool del;
1766 
1767  for(i = 0;i <= d - 1; i++)
1768  {
1769  dw = deleteInMon(w, i, lV, currRing);
1770  smon = shiftInMon(p, i, lV, currRing);
1771  del = TRUE;
1772 
1773  if(pLmDivisibleBy(smon, w))
1774  {
1775  flag = TRUE;
1776  del = FALSE;
1777 
1778  pDelete(&dw);
1779  pDelete(&smon);
1780 
1781  //delete all monomials of Jwi
1782  //and make Jwi =1
1783 
1784  for(int j = 0;j < IDELEMS(Jwi); j++)
1785  {
1786  pDelete(&Jwi->m[j]);
1787  }
1788 
1790  break;
1791  }
1792 
1793  if(pLmDivisibleBy(dw, smon))
1794  {
1795  del = FALSE;
1796  qmonp = p_Divide(smon, dw, currRing);
1797  idInsertMonomials(Jwi, shiftInMon(qmonp, -d, lV, currRing));
1798 
1799  //shiftInMon(qmonp, -d, lV, currRing):returns a new poly,
1800  //qmonp remains unchanged, delete it
1801  pDelete(&qmonp);
1802  pDelete(&dw);
1803  pDelete(&smon);
1804  }
1805  //in case both if are false, delete dw and smon
1806  if(del)
1807  {
1808  pDelete(&dw);
1809  pDelete(&smon);
1810  }
1811  }
1812 
1813 }
1814 
1815 static ideal colonIdeal(ideal S, poly w, int lV, ideal Jwi)
1816 {
1817  /*
1818  * computes the colon ideal of two-sided ideal S
1819  * w.r.t. word w and save it on Jwi
1820  * keeps S and w unchanged
1821  */
1822 
1823  if(idIs0(S))
1824  {
1825  return(S);
1826  }
1827 
1828  int i, d;
1829  d = p_Totaldegree(w, currRing);
1830  bool flag = FALSE;
1831  int SCount = IDELEMS(S);
1832  for(i = 0; i < SCount; i++)
1833  {
1834  TwordMap(S->m[i], w, lV, d, Jwi, flag);
1835  if(flag)
1836  {
1837  break;
1838  }
1839  }
1840 
1841  Jwi = minimalMonomialsGenSet(Jwi);
1842  return(Jwi);
1843 }
1844 
1845 void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp)
1846 {
1847  /*
1848  * It is based on iterative right colon operation to the
1849  * monomial ideals of the free associative algebras.
1850  * The algorithm terminates for the monomial right
1851  * ideals whose monomials define regular formal language,
1852  * that is, all the monomials of ideal can be obtained from
1853  * finite subsets by applying the finite number
1854  * of elementary operations.
1855  */
1856 
1857  int trInd;
1858  S = minimalMonomialsGenSet(S);
1859 
1860  int (*POS)(ideal, poly, std::vector<ideal>, std::vector<poly>, int);
1861  if(IG_CASE)
1862  {
1863  if(idIs0(S))
1864  {
1865  WerrorS("wrong input:not the infinitely gen. case");
1866  return;
1867  }
1868  trInd = p_Totaldegree(S->m[IDELEMS(S)-1], currRing);
1869  POS = &positionInOrbit_IG_Case;
1870  }
1871  else
1872  {
1873  POS = &positionInOrbit_FG_Case;
1874  }
1875 
1876  std::vector<ideal > idorb;
1877  std::vector< poly > polist;
1878 
1879  ideal orb_init = idInit(1, 1);
1880  idorb.push_back(orb_init);
1881 
1882  polist.push_back( p_One(currRing));
1883 
1884  std::vector< std::vector<int> > posMat;
1885  std::vector<int> posRow(lV,0);
1886  std::vector<int> C;
1887 
1888  int ds, is, ps;
1889  int lpcnt = 0;
1890 
1891  poly w, wi;
1892  ideal Jwi;
1893 
1894  while(lpcnt < idorb.size())
1895  {
1896  w = NULL;
1897  w = polist[lpcnt];
1898 
1899  if(lpcnt >= 1)
1900  {
1901  if(p_Totaldegree(idorb[lpcnt]->m[0], currRing) != 0)
1902  {
1903  C.push_back(1);
1904  }
1905  else
1906  C.push_back(0);
1907  }
1908  else
1909  C.push_back(1);
1910 
1911  ds = p_Totaldegree(w, currRing);
1912  lpcnt++;
1913 
1914  for(is = 1; is <= lV; is++)
1915  {
1916  wi = NULL;
1917  //make new copy of word w=polist[lpcnt];
1918  //in wi and update it (next colon word)
1919  //if corresponding to wi get a new ideal(colon of S),
1920  //keep it in the polist else delete it
1921 
1922  wi = pCopy(w);
1923  p_SetExp(wi, (ds*lV)+is, 1, currRing);
1924  p_Setm(wi, currRing);
1925 
1926  Jwi = NULL;
1927  //Jwi stores colon ideal of S w.r.t. wi
1928  //if get a new ideal place it in the idorb
1929  //otherwise delete it
1930  Jwi = idInit(1,1);
1931 
1932  Jwi = colonIdeal(S, wi, lV, Jwi);
1933  ps = (*POS)(Jwi, wi, idorb, polist, trInd);
1934 
1935  if(ps == 0) // finds a new ideal
1936  {
1937  posRow[is-1] = idorb.size();
1938 
1939  idorb.push_back(Jwi);
1940  polist.push_back(wi);
1941  }
1942  else // ideal is already there in the orbit
1943  {
1944  posRow[is-1]=ps-1;
1945  idDelete(&Jwi);
1946  pDelete(&wi);
1947  }
1948  }
1949  posMat.push_back(posRow);
1950  posRow.resize(lV,0);
1951  }
1952  int lO = C.size();//size of the orbit
1953  PrintLn();
1954  Print("Maximal length of words = %ld\n", p_Totaldegree(polist[lO-1], currRing));
1955  Print("\nOrbit length = %d\n", lO);
1956  PrintLn();
1957 
1958  if(odp)
1959  {
1960  Print("Words description of the Orbit: \n");
1961  for(is = 0; is < lO; is++)
1962  {
1963  pWrite0(polist[is]);
1964  PrintS(" ");
1965  }
1966  PrintLn();
1967  }
1968 
1969  for(is = idorb.size()-1; is >= 0; is--)
1970  {
1971  idDelete(&idorb[is]);
1972  }
1973  for(is = polist.size()-1; is >= 0; is--)
1974  {
1975  pDelete(&polist[is]);
1976  }
1977 
1978  idorb.resize(0);
1979  polist.resize(0);
1980 
1981  int adjMatrix[lO][lO];
1982  memset(adjMatrix, 0, lO*lO*sizeof(int));
1983  int rowCount, colCount;
1984  int tm = 0;
1985  if(!mgrad)
1986  {
1987  for(rowCount = 0; rowCount < lO; rowCount++)
1988  {
1989  for(colCount = 0; colCount < lV; colCount++)
1990  {
1991  tm = posMat[rowCount][colCount];
1992  adjMatrix[rowCount][tm] = adjMatrix[rowCount][tm] + 1;
1993  }
1994  }
1995  }
1996 
1997  ring r = currRing;
1998  int npar;
1999  char** tt;
2000  TransExtInfo p;
2001  if(!mgrad)
2002  {
2003  tt=(char**)omalloc(sizeof(char*));
2004  tt[0] = omStrDup("t");
2005  npar = 1;
2006  }
2007  else
2008  {
2009  tt=(char**)omalloc(lV*sizeof(char*));
2010  for(is = 0; is < lV; is++)
2011  {
2012  tt[is] = (char*)omalloc(7*sizeof(char)); //if required enlarge it later
2013  sprintf (tt[is], "t(%d)", is+1);
2014  }
2015  npar = lV;
2016  }
2017 
2018  p.r = rDefault(0, npar, tt);
2019  coeffs cf = nInitChar(n_transExt, &p);
2020  char** xx = (char**)omalloc(sizeof(char*));
2021  xx[0] = omStrDup("x");
2022  ring R = rDefault(cf, 1, xx);
2023  rChangeCurrRing(R);//rWrite(R);
2024  /*
2025  * matrix corresponding to the orbit of the ideal
2026  */
2027  matrix mR = mpNew(lO, lO);
2028  matrix cMat = mpNew(lO,1);
2029  poly rc;
2030 
2031  if(!mgrad)
2032  {
2033  for(rowCount = 0; rowCount < lO; rowCount++)
2034  {
2035  for(colCount = 0; colCount < lO; colCount++)
2036  {
2037  if(adjMatrix[rowCount][colCount] != 0)
2038  {
2039  MATELEM(mR, rowCount + 1, colCount + 1) = p_ISet(adjMatrix[rowCount][colCount], R);
2040  p_SetCoeff(MATELEM(mR, rowCount + 1, colCount + 1), n_Mult(pGetCoeff(mR->m[lO*rowCount+colCount]),n_Param(1, R->cf), R->cf), R);
2041  }
2042  }
2043  }
2044  }
2045  else
2046  {
2047  for(rowCount = 0; rowCount < lO; rowCount++)
2048  {
2049  for(colCount = 0; colCount < lV; colCount++)
2050  {
2051  rc=NULL;
2052  rc=p_One(R);
2053  p_SetCoeff(rc, n_Mult(pGetCoeff(rc), n_Param(colCount+1, R->cf),R->cf), R);
2054  MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1)=p_Add_q(rc,MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1), R);
2055  }
2056  }
2057  }
2058 
2059  for(rowCount = 0; rowCount < lO; rowCount++)
2060  {
2061  if(C[rowCount] != 0)
2062  {
2063  MATELEM(cMat, rowCount + 1, 1) = p_ISet(C[rowCount], R);
2064  }
2065  }
2066 
2067  matrix u;
2068  unitMatrix(lO, u); //unit matrix
2069  matrix gMat = mp_Sub(u, mR, R);
2070  char* s;
2071  if(odp)
2072  {
2073  PrintS("\nlinear system:\n");
2074  if(!mgrad)
2075  {
2076  for(rowCount = 0; rowCount < lO; rowCount++)
2077  {
2078  Print("H(%d) = ", rowCount+1);
2079  for(colCount = 0; colCount < lV; colCount++)
2080  {
2081  StringSetS(""); nWrite(n_Param(1, R->cf));
2082  s = StringEndS(); PrintS(s);
2083  Print("*"); omFree(s);
2084  Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2085  }
2086  Print(" %d\n", C[rowCount] );
2087  }
2088  PrintS("where H(1) represents the series corresp. to input ideal\n");
2089  PrintS("and i^th summand in the rhs of an eqn. is according\n");
2090  PrintS("to the right colon map corresp. to the i^th variable\n");
2091  }
2092  else
2093  {
2094  for(rowCount = 0; rowCount < lO; rowCount++)
2095  {
2096  Print("H(%d) = ", rowCount+1);
2097  for(colCount = 0; colCount < lV; colCount++)
2098  {
2099  StringSetS(""); nWrite(n_Param(colCount+1, R->cf));
2100  s = StringEndS(); PrintS(s);
2101  Print("*");omFree(s);
2102  Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2103  }
2104  Print(" %d\n", C[rowCount] );
2105  }
2106  PrintS("where H(1) represents the series corresp. to input ideal\n");
2107  }
2108  }
2109  posMat.resize(0);
2110  C.resize(0);
2111  matrix pMat;
2112  matrix lMat;
2113  matrix uMat;
2114  luDecomp(gMat, pMat, lMat, uMat, R);
2115  matrix H_serVec = mpNew(lO, 1);
2116  matrix Hnot;
2117  luSolveViaLUDecomp(pMat, lMat, uMat, cMat, H_serVec, Hnot);
2118 
2119  mp_Delete(&mR, R);
2120  mp_Delete(&u, R);
2121  mp_Delete(&pMat, R);
2122  mp_Delete(&lMat, R);
2123  mp_Delete(&uMat, R);
2124  mp_Delete(&cMat, R);
2125  mp_Delete(&gMat, R);
2126  mp_Delete(&Hnot, R);
2127  //print the Hilbert series and Orbit length
2128  PrintLn();
2129  Print("Hilbert series:");
2130  PrintLn();
2131  pWrite(H_serVec->m[0]);
2132  PrintLn();
2133  if(!mgrad)
2134  {
2135  omFree(tt[0]);
2136  }
2137  else
2138  {
2139  for(is = lV-1; is >= 0; is--)
2140 
2141  omFree( tt[is]);
2142  }
2143  omFree(tt);
2144  omFree(xx[0]);
2145  omFree(xx);
2146  rChangeCurrRing(r);
2147  rKill(R);
2148 }
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...
#define nWrite(n)
Definition: numbers.h:29
static int positionInOrbit_FG_Case(ideal I, poly, std::vector< ideal > idorb, std::vector< poly >, int)
Definition: hilb.cc:1600
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:1430
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: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:583
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
void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp)
Definition: hilb.cc:1845
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:98
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:115
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: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:1751
int k
Definition: cfEzgcd.cc:93
char * StringEndS()
Definition: reporter.cc:151
#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
void pWrite0(poly p)
Definition: polys.h:291
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:1521
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
void rKill(ring r)
Definition: ipshell.cc:6056
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: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 StringSetS(const char *st)
Definition: reporter.cc:128
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:1652
void hStaircase(scfmon stc, int *Nstc, varset var, int Nvar)
Definition: hutil.cc:319
void sortMonoIdeal_pCompare(ideal I)
Definition: hilb.cc:1639
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
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:4752
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
void PrintS(const char *s)
Definition: reporter.cc:284
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:614
static int isMonoIdBasesSame_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
Definition: hilb.cc:1488
ring rDefault(const coeffs cf, int N, char **n, int ord_size, rRingOrder_t *ord, int *block0, int *block1, int **wvhdl)
Definition: ring.cc:113
#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:1632
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:1461
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:1815
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:3594
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
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:1684
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
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:877
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:1717
#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:341
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