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 #if 0 // unused
243 //Tests if the ideal is sorted by degree
244 static bool idDegSortTest(ideal I)
245 {
246  if((I == NULL)||(idIs0(I)))
247  {
248  return(TRUE);
249  }
250  for(int i = 0; i<IDELEMS(I)-1; i++)
251  {
252  if(p_Totaldegree(I->m[i],currRing)>p_Totaldegree(I->m[i+1],currRing))
253  {
254  idPrint(I);
255  WerrorS("Ideal is not deg sorted!!");
256  return(FALSE);
257  }
258  }
259  return(TRUE);
260 }
261 #endif
262 
263 //adds the new polynomial at the coresponding position
264 //and simplifies the ideal
265 static ideal SortByDeg_p(ideal I, poly p)
266 {
267  int i,j;
268  if((I == NULL) || (idIs0(I)))
269  {
270  ideal res = idInit(1,1);
271  res->m[0] = p;
272  return(res);
273  }
274  idSkipZeroes(I);
275  #if 1
276  for(i = 0; (i<IDELEMS(I)) && (p_Totaldegree(I->m[i],currRing)<=p_Totaldegree(p,currRing)); i++)
277  {
278  if(p_DivisibleBy( I->m[i],p, currRing))
279  {
280  return(I);
281  }
282  }
283  for(i = IDELEMS(I)-1; (i>=0) && (p_Totaldegree(I->m[i],currRing)>=p_Totaldegree(p,currRing)); i--)
284  {
285  if(p_DivisibleBy(p,I->m[i], currRing))
286  {
287  I->m[i] = NULL;
288  }
289  }
290  if(idIs0(I))
291  {
292  idSkipZeroes(I);
293  I->m[0] = p;
294  return(I);
295  }
296  #endif
297  idSkipZeroes(I);
298  //First I take the case when all generators have the same degree
299  if(p_Totaldegree(I->m[0],currRing) == p_Totaldegree(I->m[IDELEMS(I)-1],currRing))
300  {
302  {
303  idInsertPoly(I,p);
304  idSkipZeroes(I);
305  for(i=IDELEMS(I)-1;i>=1; i--)
306  {
307  I->m[i] = I->m[i-1];
308  }
309  I->m[0] = p;
310  return(I);
311  }
313  {
314  idInsertPoly(I,p);
315  idSkipZeroes(I);
316  return(I);
317  }
318  }
320  {
321  idInsertPoly(I,p);
322  idSkipZeroes(I);
323  for(i=IDELEMS(I)-1;i>=1; i--)
324  {
325  I->m[i] = I->m[i-1];
326  }
327  I->m[0] = p;
328  return(I);
329  }
331  {
332  idInsertPoly(I,p);
333  idSkipZeroes(I);
334  return(I);
335  }
336  for(i = IDELEMS(I)-2; ;)
337  {
339  {
340  idInsertPoly(I,p);
341  idSkipZeroes(I);
342  for(j = IDELEMS(I)-1; j>=i+1;j--)
343  {
344  I->m[j] = I->m[j-1];
345  }
346  I->m[i] = p;
347  return(I);
348  }
350  {
351  idInsertPoly(I,p);
352  idSkipZeroes(I);
353  for(j = IDELEMS(I)-1; j>=i+2;j--)
354  {
355  I->m[j] = I->m[j-1];
356  }
357  I->m[i+1] = p;
358  return(I);
359  }
360  i--;
361  }
362 }
363 
364 //it sorts the ideal by the degrees
365 static ideal SortByDeg(ideal I)
366 {
367  if(idIs0(I))
368  {
369  return(I);
370  }
371  int i;
372  ideal res;
373  idSkipZeroes(I);
374  res = idInit(1,1);
375  res->m[0] = poly(0);
376  for(i = 0; i<=IDELEMS(I)-1;i++)
377  {
378  res = SortByDeg_p(res, I->m[i]);
379  }
380  idSkipZeroes(res);
381  //idDegSortTest(res);
382  return(res);
383 }
384 
385 //idQuot(I,p) for I monomial ideal, p a ideal with a single monomial.
386 ideal idQuotMon(ideal Iorig, ideal p)
387 {
388  if(idIs0(Iorig))
389  {
390  ideal res = idInit(1,1);
391  res->m[0] = poly(0);
392  return(res);
393  }
394  if(idIs0(p))
395  {
396  ideal res = idInit(1,1);
397  res->m[0] = pOne();
398  return(res);
399  }
400  ideal I = idCopy(Iorig);
401  ideal res = idInit(IDELEMS(I),1);
402  int i,j;
403  int dummy;
404  for(i = 0; i<IDELEMS(I); i++)
405  {
406  res->m[i] = p_Copy(I->m[i], currRing);
407  for(j = 1; (j<=currRing->N) ; j++)
408  {
409  dummy = p_GetExp(p->m[0], j, currRing);
410  if(dummy > 0)
411  {
412  if(p_GetExp(I->m[i], j, currRing) < dummy)
413  {
414  p_SetExp(res->m[i], j, 0, currRing);
415  }
416  else
417  {
418  p_SetExp(res->m[i], j, p_GetExp(I->m[i], j, currRing) - dummy, currRing);
419  }
420  }
421  }
422  p_Setm(res->m[i], currRing);
423  if(p_Totaldegree(res->m[i],currRing) == p_Totaldegree(I->m[i],currRing))
424  {
425  res->m[i] = NULL; // pDelete
426  }
427  else
428  {
429  I->m[i] = NULL; // pDelete
430  }
431  }
432  idSkipZeroes(res);
433  idSkipZeroes(I);
434  if(!idIs0(res))
435  {
436  for(i = 0; i<=IDELEMS(res)-1; i++)
437  {
438  I = SortByDeg_p(I,res->m[i]);
439  }
440  }
441  //idDegSortTest(I);
442  return(I);
443 }
444 
445 //id_Add for monomials
446 static ideal idAddMon(ideal I, ideal p)
447 {
448  #if 1
449  I = SortByDeg_p(I,p->m[0]);
450  #else
451  I = id_Add(I,p,currRing);
452  #endif
453  //idSkipZeroes(I);
454  return(I);
455 }
456 
457 //searches for a variable that is not yet used (assumes that the ideal is sqrfree)
458 static poly ChoosePVar (ideal I)
459 {
460  bool flag=TRUE;
461  int i,j;
462  poly res;
463  for(i=1;i<=currRing->N;i++)
464  {
465  flag=TRUE;
466  for(j=IDELEMS(I)-1;(j>=0)&&(flag);j--)
467  {
468  if(p_GetExp(I->m[j], i, currRing)>0)
469  {
470  flag=FALSE;
471  }
472  }
473 
474  if(flag == TRUE)
475  {
476  res = p_ISet(1, currRing);
477  p_SetExp(res, i, 1, currRing);
479  return(res);
480  }
481  else
482  {
483  p_Delete(&res, currRing);
484  }
485  }
486  return(NULL); //i.e. it is the maximal ideal
487 }
488 
489 #if 0 // unused
490 //choice XL: last entry divided by x (xy10z15 -> y9z14)
491 static poly ChoosePXL(ideal I)
492 {
493  int i,j,dummy=0;
494  poly m;
495  for(i = IDELEMS(I)-1; (i>=0) && (dummy == 0); i--)
496  {
497  for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
498  {
499  if(p_GetExp(I->m[i],j, currRing)>1)
500  {
501  dummy = 1;
502  }
503  }
504  }
505  m = p_Copy(I->m[i+1],currRing);
506  for(j = 1; j<=currRing->N; j++)
507  {
508  dummy = p_GetExp(m,j,currRing);
509  if(dummy >= 1)
510  {
511  p_SetExp(m, j, dummy-1, currRing);
512  }
513  }
514  if(!p_IsOne(m, currRing))
515  {
516  p_Setm(m, currRing);
517  return(m);
518  }
519  m = ChoosePVar(I);
520  return(m);
521 }
522 #endif
523 
524 #if 0 // unused
525 //choice XF: first entry divided by x (xy10z15 -> y9z14)
526 static poly ChoosePXF(ideal I)
527 {
528  int i,j,dummy=0;
529  poly m;
530  for(i =0 ; (i<=IDELEMS(I)-1) && (dummy == 0); i++)
531  {
532  for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
533  {
534  if(p_GetExp(I->m[i],j, currRing)>1)
535  {
536  dummy = 1;
537  }
538  }
539  }
540  m = p_Copy(I->m[i-1],currRing);
541  for(j = 1; j<=currRing->N; j++)
542  {
543  dummy = p_GetExp(m,j,currRing);
544  if(dummy >= 1)
545  {
546  p_SetExp(m, j, dummy-1, currRing);
547  }
548  }
549  if(!p_IsOne(m, currRing))
550  {
551  p_Setm(m, currRing);
552  return(m);
553  }
554  m = ChoosePVar(I);
555  return(m);
556 }
557 #endif
558 
559 #if 0 // unused
560 //choice OL: last entry the first power (xy10z15 -> xy9z15)
561 static poly ChoosePOL(ideal I)
562 {
563  int i,j,dummy;
564  poly m;
565  for(i = IDELEMS(I)-1;i>=0;i--)
566  {
567  m = p_Copy(I->m[i],currRing);
568  for(j=1;j<=currRing->N;j++)
569  {
570  dummy = p_GetExp(m,j,currRing);
571  if(dummy > 0)
572  {
573  p_SetExp(m,j,dummy-1,currRing);
574  p_Setm(m,currRing);
575  }
576  }
577  if(!p_IsOne(m, currRing))
578  {
579  return(m);
580  }
581  else
582  {
583  p_Delete(&m,currRing);
584  }
585  }
586  m = ChoosePVar(I);
587  return(m);
588 }
589 #endif
590 
591 #if 0 // unused
592 //choice OF: first entry the first power (xy10z15 -> xy9z15)
593 static poly ChoosePOF(ideal I)
594 {
595  int i,j,dummy;
596  poly m;
597  for(i = 0 ;i<=IDELEMS(I)-1;i++)
598  {
599  m = p_Copy(I->m[i],currRing);
600  for(j=1;j<=currRing->N;j++)
601  {
602  dummy = p_GetExp(m,j,currRing);
603  if(dummy > 0)
604  {
605  p_SetExp(m,j,dummy-1,currRing);
606  p_Setm(m,currRing);
607  }
608  }
609  if(!p_IsOne(m, currRing))
610  {
611  return(m);
612  }
613  else
614  {
615  p_Delete(&m,currRing);
616  }
617  }
618  m = ChoosePVar(I);
619  return(m);
620 }
621 #endif
622 
623 #if 0 // unused
624 //choice VL: last entry the first variable with power (xy10z15 -> y)
625 static poly ChoosePVL(ideal I)
626 {
627  int i,j,dummy;
628  bool flag = TRUE;
629  poly m = p_ISet(1,currRing);
630  for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
631  {
632  flag = TRUE;
633  for(j=1;(j<=currRing->N) && (flag);j++)
634  {
635  dummy = p_GetExp(I->m[i],j,currRing);
636  if(dummy >= 2)
637  {
638  p_SetExp(m,j,1,currRing);
639  p_Setm(m,currRing);
640  flag = FALSE;
641  }
642  }
643  if(!p_IsOne(m, currRing))
644  {
645  return(m);
646  }
647  }
648  m = ChoosePVar(I);
649  return(m);
650 }
651 #endif
652 
653 #if 0 // unused
654 //choice VF: first entry the first variable with power (xy10z15 -> y)
655 static poly ChoosePVF(ideal I)
656 {
657  int i,j,dummy;
658  bool flag = TRUE;
659  poly m = p_ISet(1,currRing);
660  for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
661  {
662  flag = TRUE;
663  for(j=1;(j<=currRing->N) && (flag);j++)
664  {
665  dummy = p_GetExp(I->m[i],j,currRing);
666  if(dummy >= 2)
667  {
668  p_SetExp(m,j,1,currRing);
669  p_Setm(m,currRing);
670  flag = FALSE;
671  }
672  }
673  if(!p_IsOne(m, currRing))
674  {
675  return(m);
676  }
677  }
678  m = ChoosePVar(I);
679  return(m);
680 }
681 #endif
682 
683 //choice JL: last entry just variable with power (xy10z15 -> y10)
684 static poly ChoosePJL(ideal I)
685 {
686  int i,j,dummy;
687  bool flag = TRUE;
688  poly m = p_ISet(1,currRing);
689  for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
690  {
691  flag = TRUE;
692  for(j=1;(j<=currRing->N) && (flag);j++)
693  {
694  dummy = p_GetExp(I->m[i],j,currRing);
695  if(dummy >= 2)
696  {
697  p_SetExp(m,j,dummy-1,currRing);
698  p_Setm(m,currRing);
699  flag = FALSE;
700  }
701  }
702  if(!p_IsOne(m, currRing))
703  {
704  return(m);
705  }
706  }
707  m = ChoosePVar(I);
708  return(m);
709 }
710 
711 #if 0
712 //choice JF: last entry just variable with power -1 (xy10z15 -> y9)
713 static poly ChoosePJF(ideal I)
714 {
715  int i,j,dummy;
716  bool flag = TRUE;
717  poly m = p_ISet(1,currRing);
718  for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
719  {
720  flag = TRUE;
721  for(j=1;(j<=currRing->N) && (flag);j++)
722  {
723  dummy = p_GetExp(I->m[i],j,currRing);
724  if(dummy >= 2)
725  {
726  p_SetExp(m,j,dummy-1,currRing);
727  p_Setm(m,currRing);
728  flag = FALSE;
729  }
730  }
731  if(!p_IsOne(m, currRing))
732  {
733  return(m);
734  }
735  }
736  m = ChoosePVar(I);
737  return(m);
738 }
739 #endif
740 
741 //chooses 1 \neq p \not\in S. This choice should be made optimal
742 static poly ChooseP(ideal I)
743 {
744  poly m;
745  // TEST TO SEE WHICH ONE IS BETTER
746  //m = ChoosePXL(I);
747  //m = ChoosePXF(I);
748  //m = ChoosePOL(I);
749  //m = ChoosePOF(I);
750  //m = ChoosePVL(I);
751  //m = ChoosePVF(I);
752  m = ChoosePJL(I);
753  //m = ChoosePJF(I);
754  return(m);
755 }
756 
757 ///searches for a monomial of degree d>=2 and divides it by a variable (result monomial of deg d-1)
758 static poly SearchP(ideal I)
759 {
760  int i,j,exp;
761  poly res;
762  if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)<=1)
763  {
764  res = ChoosePVar(I);
765  return(res);
766  }
767  i = IDELEMS(I)-1;
768  res = p_Copy(I->m[i], currRing);
769  for(j=1;j<=currRing->N;j++)
770  {
771  exp = p_GetExp(I->m[i], j, currRing);
772  if(exp > 0)
773  {
774  p_SetExp(res, j, exp - 1, currRing);
776  break;
777  }
778  }
779  assume( j <= currRing->N );
780  return(res);
781 }
782 
783 //test if the ideal is of the form (x1, ..., xr)
784 static bool JustVar(ideal I)
785 {
786  #if 0
787  int i,j;
788  bool foundone;
789  for(i=0;i<=IDELEMS(I)-1;i++)
790  {
791  foundone = FALSE;
792  for(j = 1;j<=currRing->N;j++)
793  {
794  if(p_GetExp(I->m[i], j, currRing)>0)
795  {
796  if(foundone == TRUE)
797  {
798  return(FALSE);
799  }
800  foundone = TRUE;
801  }
802  }
803  }
804  return(TRUE);
805  #else
806  if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)>1)
807  {
808  return(FALSE);
809  }
810  return(TRUE);
811  #endif
812 }
813 
814 //computes the Euler Characteristic of the ideal
815 static void eulerchar (ideal I, int variables, mpz_ptr ec)
816 {
817  loop
818  {
819  mpz_t dummy;
820  if(JustVar(I) == TRUE)
821  {
822  if(IDELEMS(I) == variables)
823  {
824  mpz_init(dummy);
825  if((variables % 2) == 0)
826  {mpz_set_si(dummy, 1);}
827  else
828  {mpz_set_si(dummy, -1);}
829  mpz_add(ec, ec, dummy);
830  }
831  //mpz_clear(dummy);
832  return;
833  }
834  ideal p = idInit(1,1);
835  p->m[0] = SearchP(I);
836  //idPrint(I);
837  //idPrint(p);
838  //printf("\nNow get in idQuotMon\n");
839  ideal Ip = idQuotMon(I,p);
840  //idPrint(Ip);
841  //Ip = SortByDeg(Ip);
842  int i,howmanyvarinp = 0;
843  for(i = 1;i<=currRing->N;i++)
844  {
845  if(p_GetExp(p->m[0],i,currRing)>0)
846  {
847  howmanyvarinp++;
848  }
849  }
850  eulerchar(Ip, variables-howmanyvarinp, ec);
851  id_Delete(&Ip, currRing);
852  I = idAddMon(I,p);
853  }
854 }
855 
856 //tests if an ideal is Square Free, if no, returns the variable which appears at powers >1
857 static poly SqFree (ideal I)
858 {
859  int i,j;
860  bool flag=TRUE;
861  poly notsqrfree = NULL;
862  if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)<=1)
863  {
864  return(notsqrfree);
865  }
866  for(i=IDELEMS(I)-1;(i>=0)&&(flag);i--)
867  {
868  for(j=1;(j<=currRing->N)&&(flag);j++)
869  {
870  if(p_GetExp(I->m[i],j,currRing)>1)
871  {
872  flag=FALSE;
873  notsqrfree = p_ISet(1,currRing);
874  p_SetExp(notsqrfree,j,1,currRing);
875  }
876  }
877  }
878  if(notsqrfree != NULL)
879  {
880  p_Setm(notsqrfree,currRing);
881  }
882  return(notsqrfree);
883 }
884 
885 //checks if a polynomial is in I
886 static bool IsIn(poly p, ideal I)
887 {
888  //assumes that I is ordered by degree
889  if(idIs0(I))
890  {
891  if(p==poly(0))
892  {
893  return(TRUE);
894  }
895  else
896  {
897  return(FALSE);
898  }
899  }
900  if(p==poly(0))
901  {
902  return(FALSE);
903  }
904  int i,j;
905  bool flag;
906  for(i = 0;i<IDELEMS(I);i++)
907  {
908  flag = TRUE;
909  for(j = 1;(j<=currRing->N) &&(flag);j++)
910  {
911  if(p_GetExp(p, j, currRing)<p_GetExp(I->m[i], j, currRing))
912  {
913  flag = FALSE;
914  }
915  }
916  if(flag)
917  {
918  return(TRUE);
919  }
920  }
921  return(FALSE);
922 }
923 
924 //computes the lcm of min I, I monomial ideal
925 static poly LCMmon(ideal I)
926 {
927  if(idIs0(I))
928  {
929  return(NULL);
930  }
931  poly m;
932  int dummy,i,j;
933  m = p_ISet(1,currRing);
934  for(i=1;i<=currRing->N;i++)
935  {
936  dummy=0;
937  for(j=IDELEMS(I)-1;j>=0;j--)
938  {
939  if(p_GetExp(I->m[j],i,currRing) > dummy)
940  {
941  dummy = p_GetExp(I->m[j],i,currRing);
942  }
943  }
944  p_SetExp(m,i,dummy,currRing);
945  }
946  p_Setm(m,currRing);
947  return(m);
948 }
949 
950 //the Roune Slice Algorithm
951 void rouneslice(ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int* &hilbpower)
952 {
953  loop
954  {
955  (steps)++;
956  int i,j;
957  int dummy;
958  poly m;
959  ideal p;
960  //----------- PRUNING OF S ---------------
961  //S SHOULD IN THIS POINT BE ORDERED BY DEGREE
962  for(i=IDELEMS(S)-1;i>=0;i--)
963  {
964  if(IsIn(S->m[i],I))
965  {
966  S->m[i]=NULL;
967  prune++;
968  }
969  }
970  idSkipZeroes(S);
971  //----------------------------------------
972  for(i=IDELEMS(I)-1;i>=0;i--)
973  {
974  m = p_Copy(I->m[i],currRing);
975  for(j=1;j<=currRing->N;j++)
976  {
977  dummy = p_GetExp(m,j,currRing);
978  if(dummy > 0)
979  {
980  p_SetExp(m,j,dummy-1,currRing);
981  }
982  }
983  p_Setm(m, currRing);
984  if(IsIn(m,S))
985  {
986  I->m[i]=NULL;
987  //printf("\n Deleted, since pi(m) is in S\n");pWrite(m);
988  }
989  }
990  idSkipZeroes(I);
991  //----------- MORE PRUNING OF S ------------
992  m = LCMmon(I);
993  if(m != NULL)
994  {
995  for(i=0;i<IDELEMS(S);i++)
996  {
997  if(!(p_DivisibleBy(S->m[i], m, currRing)))
998  {
999  S->m[i] = NULL;
1000  j++;
1001  moreprune++;
1002  }
1003  else
1004  {
1005  if(pLmEqual(S->m[i],m))
1006  {
1007  S->m[i] = NULL;
1008  moreprune++;
1009  }
1010  }
1011  }
1012  idSkipZeroes(S);
1013  }
1014  /*printf("\n---------------------------\n");
1015  printf("\n I\n");idPrint(I);
1016  printf("\n S\n");idPrint(S);
1017  printf("\n q\n");pWrite(q);
1018  getchar();*/
1019 
1020  if(idIs0(I))
1021  {
1022  id_Delete(&I, currRing);
1023  id_Delete(&S, currRing);
1024  p_Delete(&m, currRing);
1025  break;
1026  }
1027  m = LCMmon(I);
1028  if(!p_DivisibleBy(x,m, currRing))
1029  {
1030  //printf("\nx does not divide lcm(I)");
1031  //printf("\nEmpty set");pWrite(q);
1032  id_Delete(&I, currRing);
1033  id_Delete(&S, currRing);
1034  p_Delete(&m, currRing);
1035  break;
1036  }
1037  m = SqFree(I);
1038  if(m==NULL)
1039  {
1040  //printf("\n Corner: ");
1041  //pWrite(q);
1042  //printf("\n With the facets of the dual simplex:\n");
1043  //idPrint(I);
1044  mpz_t ec;
1045  mpz_init(ec);
1046  mpz_ptr ec_ptr = ec;
1047  eulerchar(I, currRing->N, ec_ptr);
1048  bool flag = FALSE;
1049  if(NNN==0)
1050  {
1051  hilbertcoef = (mpz_ptr)omAlloc((NNN+1)*sizeof(mpz_t));
1052  hilbpower = (int*)omAlloc((NNN+1)*sizeof(int));
1053  mpz_init( &hilbertcoef[NNN]);
1054  mpz_set( &hilbertcoef[NNN], ec);
1055  mpz_clear(ec);
1056  hilbpower[NNN] = p_Totaldegree(q,currRing);
1057  NNN++;
1058  }
1059  else
1060  {
1061  //I look if the power appears already
1062  for(i = 0;(i<NNN)&&(flag == FALSE)&&(p_Totaldegree(q,currRing)>=hilbpower[i]);i++)
1063  {
1064  if((hilbpower[i]) == (p_Totaldegree(q,currRing)))
1065  {
1066  flag = TRUE;
1067  mpz_add(&hilbertcoef[i],&hilbertcoef[i],ec_ptr);
1068  }
1069  }
1070  if(flag == FALSE)
1071  {
1072  hilbertcoef = (mpz_ptr)omRealloc(hilbertcoef, (NNN+1)*sizeof(mpz_t));
1073  hilbpower = (int*)omRealloc(hilbpower, (NNN+1)*sizeof(int));
1074  mpz_init(&hilbertcoef[NNN]);
1075  for(j = NNN; j>i; j--)
1076  {
1077  mpz_set(&hilbertcoef[j],&hilbertcoef[j-1]);
1078  hilbpower[j] = hilbpower[j-1];
1079  }
1080  mpz_set( &hilbertcoef[i], ec);
1081  mpz_clear(ec);
1082  hilbpower[i] = p_Totaldegree(q,currRing);
1083  NNN++;
1084  }
1085  }
1086  break;
1087  }
1088  m = ChooseP(I);
1089  p = idInit(1,1);
1090  p->m[0] = m;
1091  ideal Ip = idQuotMon(I,p);
1092  ideal Sp = idQuotMon(S,p);
1093  poly pq = pp_Mult_mm(q,m,currRing);
1094  rouneslice(Ip, Sp, pq, x, prune, moreprune, steps, NNN, hilbertcoef,hilbpower);
1095  //id_Delete(&Ip, currRing);
1096  //id_Delete(&Sp, currRing);
1097  S = idAddMon(S,p);
1098  p->m[0]=NULL;
1099  id_Delete(&p, currRing); // p->m[0] was also in S
1100  p_Delete(&pq,currRing);
1101  }
1102 }
1103 
1104 //it computes the first hilbert series by means of Roune Slice Algorithm
1105 void slicehilb(ideal I)
1106 {
1107  //printf("Adi changes are here: \n");
1108  int i, NNN = 0;
1109  int steps = 0, prune = 0, moreprune = 0;
1110  mpz_ptr hilbertcoef;
1111  int *hilbpower;
1112  ideal S = idInit(1,1);
1113  poly q = p_ISet(1,currRing);
1114  ideal X = idInit(1,1);
1115  X->m[0]=p_One(currRing);
1116  for(i=1;i<=currRing->N;i++)
1117  {
1118  p_SetExp(X->m[0],i,1,currRing);
1119  }
1120  p_Setm(X->m[0],currRing);
1121  I = id_Mult(I,X,currRing);
1122  I = SortByDeg(I);
1123  //printf("\n-------------RouneSlice--------------\n");
1124  rouneslice(I,S,q,X->m[0],prune, moreprune, steps, NNN, hilbertcoef, hilbpower);
1125  //printf("\nIn total Prune got rid of %i elements\n",prune);
1126  //printf("\nIn total More Prune got rid of %i elements\n",moreprune);
1127  //printf("\nSteps of rouneslice: %i\n\n", steps);
1128  mpz_t coefhilb;
1129  mpz_t dummy;
1130  mpz_init(coefhilb);
1131  mpz_init(dummy);
1132  printf("\n// %8d t^0",1);
1133  for(i = 0; i<NNN; i++)
1134  {
1135  if(mpz_sgn(&hilbertcoef[i])!=0)
1136  {
1137  gmp_printf("\n// %8Zd t^%d",&hilbertcoef[i],hilbpower[i]);
1138  }
1139  }
1140  omFreeSize(hilbertcoef, (NNN)*sizeof(mpz_t));
1141  omFreeSize(hilbpower, (NNN)*sizeof(int));
1142  //printf("\n-------------------------------------\n");
1143 }
1144 
1145 // -------------------------------- END OF CHANGES -------------------------------------------
1146 static intvec * hSeries(ideal S, intvec *modulweight,
1147  int /*notstc*/, intvec *wdegree, ideal Q, ring tailRing)
1148 {
1149 // id_TestTail(S, currRing, tailRing);
1150 
1151  intvec *work, *hseries1=NULL;
1152  int mc;
1153  int p0;
1154  int i, j, k, l, ii, mw;
1155  hexist = hInit(S, Q, &hNexist, tailRing);
1156  if (hNexist==0)
1157  {
1158  hseries1=new intvec(2);
1159  (*hseries1)[0]=1;
1160  (*hseries1)[1]=0;
1161  return hseries1;
1162  }
1163 
1164  #if 0
1165  if (wdegree == NULL)
1166  hWeight();
1167  else
1168  hWDegree(wdegree);
1169  #else
1170  if (wdegree != NULL) hWDegree(wdegree);
1171  #endif
1172 
1173  p0 = 1;
1174  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
1175  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
1176  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
1177  stcmem = hCreate((currRing->N) - 1);
1178  Qpol = (int **)omAlloc(((currRing->N) + 1) * sizeof(int *));
1179  Ql = (int *)omAlloc0(((currRing->N) + 1) * sizeof(int));
1180  Q0 = (int *)omAlloc(((currRing->N) + 1) * sizeof(int));
1181  *Qpol = NULL;
1182  hLength = k = j = 0;
1183  mc = hisModule;
1184  if (mc!=0)
1185  {
1186  mw = hMinModulweight(modulweight);
1187  hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1188  }
1189  else
1190  {
1191  mw = 0;
1192  hstc = hexist;
1193  hNstc = hNexist;
1194  }
1195  loop
1196  {
1197  if (mc!=0)
1198  {
1199  hComp(hexist, hNexist, mc, hstc, &hNstc);
1200  if (modulweight != NULL)
1201  j = (*modulweight)[mc-1]-mw;
1202  }
1203  if (hNstc!=0)
1204  {
1205  hNvar = (currRing->N);
1206  for (i = hNvar; i>=0; i--)
1207  hvar[i] = i;
1208  //if (notstc) // TODO: no mon divides another
1210  hSupp(hstc, hNstc, hvar, &hNvar);
1211  if (hNvar!=0)
1212  {
1213  if ((hNvar > 2) && (hNstc > 10))
1216  memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
1217  hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
1218  hLexS(hstc, hNstc, hvar, hNvar);
1219  Q0[hNvar] = 0;
1220  hHilbStep(hpure, hstc, hNstc, hvar, hNvar, &p0, 1);
1221  }
1222  }
1223  else
1224  {
1225  if(*Qpol!=NULL)
1226  (**Qpol)++;
1227  else
1228  {
1229  *Qpol = (int *)omAlloc(sizeof(int));
1230  hLength = *Ql = **Qpol = 1;
1231  }
1232  }
1233  if (*Qpol!=NULL)
1234  {
1235  i = hLength;
1236  while ((i > 0) && ((*Qpol)[i - 1] == 0))
1237  i--;
1238  if (i > 0)
1239  {
1240  l = i + j;
1241  if (l > k)
1242  {
1243  work = new intvec(l);
1244  for (ii=0; ii<k; ii++)
1245  (*work)[ii] = (*hseries1)[ii];
1246  if (hseries1 != NULL)
1247  delete hseries1;
1248  hseries1 = work;
1249  k = l;
1250  }
1251  while (i > 0)
1252  {
1253  (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
1254  (*Qpol)[i - 1] = 0;
1255  i--;
1256  }
1257  }
1258  }
1259  mc--;
1260  if (mc <= 0)
1261  break;
1262  }
1263  if (k==0)
1264  {
1265  hseries1=new intvec(2);
1266  (*hseries1)[0]=0;
1267  (*hseries1)[1]=0;
1268  }
1269  else
1270  {
1271  l = k+1;
1272  while ((*hseries1)[l-2]==0) l--;
1273  if (l!=k)
1274  {
1275  work = new intvec(l);
1276  for (ii=l-2; ii>=0; ii--)
1277  (*work)[ii] = (*hseries1)[ii];
1278  delete hseries1;
1279  hseries1 = work;
1280  }
1281  (*hseries1)[l-1] = mw;
1282  }
1283  for (i = 0; i <= (currRing->N); i++)
1284  {
1285  if (Ql[i]!=0)
1286  omFreeSize((ADDRESS)Qpol[i], Ql[i] * sizeof(int));
1287  }
1288  omFreeSize((ADDRESS)Q0, ((currRing->N) + 1) * sizeof(int));
1289  omFreeSize((ADDRESS)Ql, ((currRing->N) + 1) * sizeof(int));
1290  omFreeSize((ADDRESS)Qpol, ((currRing->N) + 1) * sizeof(int *));
1291  hKill(stcmem, (currRing->N) - 1);
1292  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
1293  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
1294  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1296  if (hisModule!=0)
1297  omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1298  return hseries1;
1299 }
1300 
1301 
1302 intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
1303 {
1305  if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
1306  return hSeries(S, modulweight, 0, wdegree, Q, tailRing);
1307 }
1308 
1309 intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
1310 {
1312  if (Q!= NULL) id_TestTail(Q, currRing, tailRing);
1313 
1314  return hSeries(S, modulweight, 1, wdegree, Q, tailRing);
1315 }
1316 
1318 {
1319  intvec *work, *hseries2;
1320  int i, j, k, s, t, l;
1321  if (hseries1 == NULL)
1322  return NULL;
1323  work = new intvec(hseries1);
1324  k = l = work->length()-1;
1325  s = 0;
1326  for (i = k-1; i >= 0; i--)
1327  s += (*work)[i];
1328  loop
1329  {
1330  if ((s != 0) || (k == 1))
1331  break;
1332  s = 0;
1333  t = (*work)[k-1];
1334  k--;
1335  for (i = k-1; i >= 0; i--)
1336  {
1337  j = (*work)[i];
1338  (*work)[i] = -t;
1339  s += t;
1340  t += j;
1341  }
1342  }
1343  hseries2 = new intvec(k+1);
1344  for (i = k-1; i >= 0; i--)
1345  (*hseries2)[i] = (*work)[i];
1346  (*hseries2)[k] = (*work)[l];
1347  delete work;
1348  return hseries2;
1349 }
1350 
1351 void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
1352 {
1353  int m, i, j, k;
1354  *co = *mu = 0;
1355  if ((s1 == NULL) || (s2 == NULL))
1356  return;
1357  i = s1->length();
1358  j = s2->length();
1359  if (j > i)
1360  return;
1361  m = 0;
1362  for(k=j-2; k>=0; k--)
1363  m += (*s2)[k];
1364  *mu = m;
1365  *co = i - j;
1366 }
1367 
1368 static void hPrintHilb(intvec *hseries)
1369 {
1370  int i, j, l, k;
1371  if (hseries == NULL)
1372  return;
1373  l = hseries->length()-1;
1374  k = (*hseries)[l];
1375  for (i = 0; i < l; i++)
1376  {
1377  j = (*hseries)[i];
1378  if (j != 0)
1379  {
1380  Print("// %8d t^%d\n", j, i+k);
1381  }
1382  }
1383 }
1384 
1385 /*
1386 *caller
1387 */
1388 void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
1389 {
1391 
1392  intvec *hseries1 = hFirstSeries(S, modulweight, Q, wdegree, tailRing);
1393 
1394  hPrintHilb(hseries1);
1395 
1396  const int l = hseries1->length()-1;
1397 
1398  intvec *hseries2 = (l > 1) ? hSecondSeries(hseries1) : hseries1;
1399 
1400  int co, mu;
1401  hDegreeSeries(hseries1, hseries2, &co, &mu);
1402 
1403  PrintLn();
1404  hPrintHilb(hseries2);
1405  if ((l == 1) &&(mu == 0))
1406  scPrintDegree(rVar(currRing)+1, 0);
1407  else
1408  scPrintDegree(co, mu);
1409  if (l>1)
1410  delete hseries1;
1411  delete hseries2;
1412 }
1413 
1414 /***********************************************************************
1415  Computation of Hilbert series of non-commutative monomial algebras
1416 ************************************************************************/
1417 
1418 
1419 static void idInsertMonomials(ideal I, poly p)
1420 {
1421  /*
1422  * adds monomial in I and if required,
1423  * enlarges the size of poly-set by 16
1424  * does not make copy of p
1425  */
1426 
1427  if(I == NULL)
1428  {
1429  return;
1430  }
1431 
1432  int j = IDELEMS(I) - 1;
1433  while ((j >= 0) && (I->m[j] == NULL))
1434  {
1435  j--;
1436  }
1437  j++;
1438  if (j == IDELEMS(I))
1439  {
1440  pEnlargeSet(&(I->m), IDELEMS(I), 16);
1441  IDELEMS(I) +=16;
1442  }
1443  I->m[j] = p;
1444 }
1445 
1446 static int isMonoIdBasesSame(ideal J, ideal Ob)
1447 {
1448  /*
1449  * polynomials of J and Ob are assumed to
1450  * be already sorted. J and Ob are
1451  * represented by the minimal generating set
1452  */
1453  int i, s;
1454  s = 1;
1455  int JCount = IDELEMS(J);
1456  int ObCount = IDELEMS(Ob);
1457 
1458  if(idIs0(J))
1459  {
1460  return(1);
1461  }
1462  if(JCount != ObCount)
1463  {
1464  return(0);
1465  }
1466 
1467  for(i = 0; i < JCount; i++)
1468  {
1469  if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1470  {
1471  return(0);
1472  }
1473  }
1474  return(s);
1475 }
1476 
1477 static int CountOnIdUptoTruncationIndex(ideal I, int tr)
1478 {
1479  /*
1480  * I must be sorted in ascending order
1481  * counts the number of polys in I upto
1482  * degree less or equal to tr
1483  */
1484 
1485  //case when I=1;
1486  if(p_Totaldegree(I->m[0], currRing) == 0)
1487  {
1488  return(1);
1489  }
1490 
1491  int count = 0;
1492  for(int i = 0; i < IDELEMS(I); i++)
1493  {
1494  if(p_Totaldegree(I->m[i], currRing) > tr)
1495  {
1496  return (count);
1497  }
1498  count = count + 1;
1499  }
1500 
1501  return(count);
1502 }
1503 
1504 static int isMonoIdBasesSame_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
1505 {
1506  /*
1507  * polynomials of J and obc are assumed to
1508  * be already sorted. J and Ob are
1509  * represented by the minimal generating set.
1510  * checks if J and Ob are same in polys upto deg <=tr
1511  */
1512 
1513  int i, s;
1514  s = 1;
1515  //when J is null
1516  if(JCount == 0)
1517  {
1518  return(1);
1519  }
1520 
1521  if(JCount != ObCount)
1522  {
1523  return(0);
1524  }
1525 
1526  for(i = 0; i< JCount; i++)
1527  {
1528  if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1529  {
1530  return(0);
1531  }
1532  }
1533 
1534  return(s);
1535 }
1536 
1537 static int positionInOrbit_IG_Case(ideal I, poly w, std::vector<ideal> idorb, std::vector<poly> polist, int trInd)
1538 {
1539  /*
1540  * compares the ideal I with ideals in the Orbit 'idorb'
1541  * upto degree trInd - max(deg of w, deg of word in polist) polynomials;
1542  * I and ideals in the Orbit are sorted,
1543  * Orbit is ordered,
1544  *
1545  * returns 0 if I is not equal to any of the ideals
1546  * in the Orbit else returns position of the matched ideal
1547  */
1548 
1549  int ps = 0;
1550  int i, s = 0;
1551  int orbCount = idorb.size();
1552 
1553  if(idIs0(I))
1554  {
1555  return(1);
1556  }
1557 
1558  int degw = p_Totaldegree(w, currRing);
1559  int degp;
1560  int dtr;
1561  int dtrp;
1562 
1563  dtr = trInd - degw;
1564  int IwCount;
1565 
1566  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1567 
1568  if(IwCount == 0)
1569  {
1570  return(1);
1571  }
1572 
1573  int ObCount;
1574 
1575  bool flag2 = FALSE;
1576 
1577  for(i = 1;i < orbCount; i++)
1578  {
1579  degp = p_Totaldegree(polist[i], currRing);
1580  if(degw > degp)
1581  {
1582  dtr = trInd - degw;
1583 
1584  ObCount = 0;
1585  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1586  if(ObCount == 0)
1587  {continue;}
1588  if(flag2)
1589  {
1590  IwCount = 0;
1591  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1592  flag2 = FALSE;
1593  }
1594  }
1595  else
1596  {
1597  flag2 = TRUE;
1598  dtrp = trInd - degp;
1599  ObCount = 0;
1600  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtrp);
1601  IwCount = 0;
1602  IwCount = CountOnIdUptoTruncationIndex(I, dtrp);
1603  }
1604 
1605  s = isMonoIdBasesSame_IG_Case(I, IwCount, idorb[i], ObCount);
1606 
1607  if(s)
1608  {
1609  ps = i + 1;
1610  break;
1611  }
1612  }
1613  return(ps);
1614 }
1615 
1616 static int positionInOrbit_FG_Case(ideal I, poly, std::vector<ideal> idorb, std::vector<poly>, int)
1617 {
1618  /*
1619  * compares the ideal I with ideals in the Orbit 'idorb'
1620  * I and ideals in the Orbit are sorted,
1621  * Orbit is ordered,
1622  *
1623  * returns 0 if I is not equal to any of the ideals
1624  * in the Orbit else returns position of the matched ideal
1625  */
1626  int ps = 0;
1627  int i, s = 0;
1628  int OrbCount = idorb.size();
1629 
1630  if(idIs0(I))
1631  {
1632  return(1);
1633  }
1634 
1635  for(i = 1; i < OrbCount; i++)
1636  {
1637  s = isMonoIdBasesSame(I, idorb[i]);
1638  if(s)
1639  {
1640  ps = i + 1;
1641  break;
1642  }
1643  }
1644 
1645  return(ps);
1646 }
1647 
1648 static int monCompare( const void *m, const void *n)
1649 {
1650  /* compares monomials */
1651 
1652  return(p_Compare(*(poly*) m, *(poly*)n, currRing));
1653 }
1654 
1656 {
1657  /*
1658  * sorts the monomial ideal in ascending order
1659  * order must be a total degree
1660  */
1661 
1662  qsort(I->m, IDELEMS(I), sizeof(poly), monCompare);
1663 
1664 }
1665 
1666 
1667 
1668 static ideal minimalMonomialsGenSet(ideal I)
1669 {
1670  /*
1671  * eliminates monomials which
1672  * can be generated by others in I
1673  */
1674  //first sort monomials of the ideal
1675 
1676  idSkipZeroes(I);
1677 
1679 
1680  int i, k;
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]));
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 
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, d;
1845  d = p_Totaldegree(w, currRing);
1846  bool flag = FALSE;
1847  int SCount = IDELEMS(S);
1848  for(i = 0; i < SCount; i++)
1849  {
1850  TwordMap(S->m[i], w, lV, d, Jwi, flag);
1851  if(flag)
1852  {
1853  break;
1854  }
1855  }
1856 
1857  Jwi = minimalMonomialsGenSet(Jwi);
1858  return(Jwi);
1859 }
1860 
1861 void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp)
1862 {
1863  /*
1864  * It is based on iterative right colon operation to the
1865  * monomial ideals of the free associative algebras.
1866  * The algorithm terminates for the monomial right
1867  * ideals whose monomials define regular formal language,
1868  * that is, all the monomials of ideal can be obtained from
1869  * finite subsets by applying the finite number
1870  * of elementary operations.
1871  */
1872 
1873  int trInd;
1874  S = minimalMonomialsGenSet(S);
1875 
1876  int (*POS)(ideal, poly, std::vector<ideal>, std::vector<poly>, int);
1877  if(IG_CASE)
1878  {
1879  if(idIs0(S))
1880  {
1881  WerrorS("wrong input:not the infinitely gen. case");
1882  return;
1883  }
1884  trInd = p_Totaldegree(S->m[IDELEMS(S)-1], currRing);
1885  POS = &positionInOrbit_IG_Case;
1886  }
1887  else
1888  {
1889  POS = &positionInOrbit_FG_Case;
1890  }
1891 
1892  std::vector<ideal > idorb;
1893  std::vector< poly > polist;
1894 
1895  ideal orb_init = idInit(1, 1);
1896  idorb.push_back(orb_init);
1897 
1898  polist.push_back( p_One(currRing));
1899 
1900  std::vector< std::vector<int> > posMat;
1901  std::vector<int> posRow(lV,0);
1902  std::vector<int> C;
1903 
1904  int ds, is, ps;
1905  int lpcnt = 0;
1906 
1907  poly w, wi;
1908  ideal Jwi;
1909 
1910  while(lpcnt < idorb.size())
1911  {
1912  w = NULL;
1913  w = polist[lpcnt];
1914 
1915  if(lpcnt >= 1)
1916  {
1917  if(p_Totaldegree(idorb[lpcnt]->m[0], currRing) != 0)
1918  {
1919  C.push_back(1);
1920  }
1921  else
1922  C.push_back(0);
1923  }
1924  else
1925  C.push_back(1);
1926 
1927  ds = p_Totaldegree(w, currRing);
1928  lpcnt++;
1929 
1930  for(is = 1; is <= lV; is++)
1931  {
1932  wi = NULL;
1933  //make new copy of word w=polist[lpcnt];
1934  //in wi and update it (next colon word)
1935  //if corresponding to wi get a new ideal(colon of S),
1936  //keep it in the polist else delete it
1937 
1938  wi = pCopy(w);
1939  p_SetExp(wi, (ds*lV)+is, 1, currRing);
1940  p_Setm(wi, currRing);
1941 
1942  Jwi = NULL;
1943  //Jwi stores colon ideal of S w.r.t. wi
1944  //if get a new ideal place it in the idorb
1945  //otherwise delete it
1946  Jwi = idInit(1,1);
1947 
1948  Jwi = colonIdeal(S, wi, lV, Jwi);
1949  ps = (*POS)(Jwi, wi, idorb, polist, trInd);
1950 
1951  if(ps == 0) // finds a new ideal
1952  {
1953  posRow[is-1] = idorb.size();
1954 
1955  idorb.push_back(Jwi);
1956  polist.push_back(wi);
1957  }
1958  else // ideal is already there in the orbit
1959  {
1960  posRow[is-1]=ps-1;
1961  idDelete(&Jwi);
1962  pDelete(&wi);
1963  }
1964  }
1965  posMat.push_back(posRow);
1966  posRow.resize(lV,0);
1967  }
1968  int lO = C.size();//size of the orbit
1969  PrintLn();
1970  Print("Maximal length of words = %ld\n", p_Totaldegree(polist[lO-1], currRing));
1971  Print("\nOrbit length = %d\n", lO);
1972  PrintLn();
1973 
1974  if(odp)
1975  {
1976  Print("Words description of the Orbit: \n");
1977  for(is = 0; is < lO; is++)
1978  {
1979  pWrite0(polist[is]);
1980  PrintS(" ");
1981  }
1982  PrintLn();
1983  }
1984 
1985  for(is = idorb.size()-1; is >= 0; is--)
1986  {
1987  idDelete(&idorb[is]);
1988  }
1989  for(is = polist.size()-1; is >= 0; is--)
1990  {
1991  pDelete(&polist[is]);
1992  }
1993 
1994  idorb.resize(0);
1995  polist.resize(0);
1996 
1997  int adjMatrix[lO][lO];
1998  memset(adjMatrix, 0, lO*lO*sizeof(int));
1999  int rowCount, colCount;
2000  int tm = 0;
2001  if(!mgrad)
2002  {
2003  for(rowCount = 0; rowCount < lO; rowCount++)
2004  {
2005  for(colCount = 0; colCount < lV; colCount++)
2006  {
2007  tm = posMat[rowCount][colCount];
2008  adjMatrix[rowCount][tm] = adjMatrix[rowCount][tm] + 1;
2009  }
2010  }
2011  }
2012 
2013  ring r = currRing;
2014  int npar;
2015  char** tt;
2016  TransExtInfo p;
2017  if(!mgrad)
2018  {
2019  tt=(char**)omalloc(sizeof(char*));
2020  tt[0] = omStrDup("t");
2021  npar = 1;
2022  }
2023  else
2024  {
2025  tt=(char**)omalloc(lV*sizeof(char*));
2026  for(is = 0; is < lV; is++)
2027  {
2028  tt[is] = (char*)omalloc(7*sizeof(char)); //if required enlarge it later
2029  sprintf (tt[is], "t(%d)", is+1);
2030  }
2031  npar = lV;
2032  }
2033 
2034  p.r = rDefault(0, npar, tt);
2036  char** xx = (char**)omalloc(sizeof(char*));
2037  xx[0] = omStrDup("x");
2038  ring R = rDefault(cf, 1, xx);
2039  rChangeCurrRing(R);//rWrite(R);
2040  /*
2041  * matrix corresponding to the orbit of the ideal
2042  */
2043  matrix mR = mpNew(lO, lO);
2044  matrix cMat = mpNew(lO,1);
2045  poly rc;
2046 
2047  if(!mgrad)
2048  {
2049  for(rowCount = 0; rowCount < lO; rowCount++)
2050  {
2051  for(colCount = 0; colCount < lO; colCount++)
2052  {
2053  if(adjMatrix[rowCount][colCount] != 0)
2054  {
2055  MATELEM(mR, rowCount + 1, colCount + 1) = p_ISet(adjMatrix[rowCount][colCount], R);
2056  p_SetCoeff(MATELEM(mR, rowCount + 1, colCount + 1), n_Mult(pGetCoeff(mR->m[lO*rowCount+colCount]),n_Param(1, R->cf), R->cf), R);
2057  }
2058  }
2059  }
2060  }
2061  else
2062  {
2063  for(rowCount = 0; rowCount < lO; rowCount++)
2064  {
2065  for(colCount = 0; colCount < lV; colCount++)
2066  {
2067  rc=NULL;
2068  rc=p_One(R);
2069  p_SetCoeff(rc, n_Mult(pGetCoeff(rc), n_Param(colCount+1, R->cf),R->cf), R);
2070  MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1)=p_Add_q(rc,MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1), R);
2071  }
2072  }
2073  }
2074 
2075  for(rowCount = 0; rowCount < lO; rowCount++)
2076  {
2077  if(C[rowCount] != 0)
2078  {
2079  MATELEM(cMat, rowCount + 1, 1) = p_ISet(C[rowCount], R);
2080  }
2081  }
2082 
2083  matrix u;
2084  unitMatrix(lO, u); //unit matrix
2085  matrix gMat = mp_Sub(u, mR, R);
2086  char* s;
2087  if(odp)
2088  {
2089  PrintS("\nlinear system:\n");
2090  if(!mgrad)
2091  {
2092  for(rowCount = 0; rowCount < lO; rowCount++)
2093  {
2094  Print("H(%d) = ", rowCount+1);
2095  for(colCount = 0; colCount < lV; colCount++)
2096  {
2097  StringSetS(""); nWrite(n_Param(1, R->cf));
2098  s = StringEndS(); PrintS(s);
2099  Print("*"); omFree(s);
2100  Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2101  }
2102  Print(" %d\n", C[rowCount] );
2103  }
2104  PrintS("where H(1) represents the series corresp. to input ideal\n");
2105  PrintS("and i^th summand in the rhs of an eqn. is according\n");
2106  PrintS("to the right colon map corresp. to the i^th variable\n");
2107  }
2108  else
2109  {
2110  for(rowCount = 0; rowCount < lO; rowCount++)
2111  {
2112  Print("H(%d) = ", rowCount+1);
2113  for(colCount = 0; colCount < lV; colCount++)
2114  {
2115  StringSetS(""); nWrite(n_Param(colCount+1, R->cf));
2116  s = StringEndS(); PrintS(s);
2117  Print("*");omFree(s);
2118  Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2119  }
2120  Print(" %d\n", C[rowCount] );
2121  }
2122  PrintS("where H(1) represents the series corresp. to input ideal\n");
2123  }
2124  }
2125  posMat.resize(0);
2126  C.resize(0);
2127  matrix pMat;
2128  matrix lMat;
2129  matrix uMat;
2130  luDecomp(gMat, pMat, lMat, uMat, R);
2131  matrix H_serVec = mpNew(lO, 1);
2132  matrix Hnot;
2133  luSolveViaLUDecomp(pMat, lMat, uMat, cMat, H_serVec, Hnot);
2134 
2135  mp_Delete(&mR, R);
2136  mp_Delete(&u, R);
2137  mp_Delete(&pMat, R);
2138  mp_Delete(&lMat, R);
2139  mp_Delete(&uMat, R);
2140  mp_Delete(&cMat, R);
2141  mp_Delete(&gMat, R);
2142  mp_Delete(&Hnot, R);
2143  //print the Hilbert series and Orbit length
2144  PrintLn();
2145  Print("Hilbert series:");
2146  PrintLn();
2147  pWrite(H_serVec->m[0]);
2148  PrintLn();
2149  if(!mgrad)
2150  {
2151  omFree(tt[0]);
2152  }
2153  else
2154  {
2155  for(is = lV-1; is >= 0; is--)
2156 
2157  omFree( tt[is]);
2158  }
2159  omFree(tt);
2160  omFree(xx[0]);
2161  omFree(xx);
2162  rChangeCurrRing(r);
2163  rKill(R);
2164 }
int status int void size_t count
Definition: si_signals.h:59
ideal idQuotMon(ideal Iorig, ideal p)
Definition: hilb.cc:386
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:951
#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:1616
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:458
#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:1446
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:1105
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
monf hCreate(int Nvar)
Definition: hutil.cc:1002
static bool JustVar(ideal I)
Definition: hilb.cc:784
int hNvar
Definition: hutil.cc:22
void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp)
Definition: hilb.cc:1861
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:886
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:446
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:1351
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:1767
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
#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:1537
static poly LCMmon(ideal I)
Definition: hilb.cc:925
intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1302
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:1314
void hKill(monf xmem, int Nvar)
Definition: hutil.cc:1016
Variable var
Definition: int_poly.h:74
void rKill(ring r)
Definition: ipshell.cc:6057
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:684
static int * hAddHilb(int Nv, int x, int *pol, int *lp)
Definition: hilb.cc:103
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:1668
void hStaircase(scfmon stc, int *Nstc, varset var, int Nvar)
Definition: hutil.cc:319
void sortMonoIdeal_pCompare(ideal I)
Definition: hilb.cc:1655
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:1146
int * scmon
Definition: hutil.h:17
struct for passing initialization parameters to naInitChar
Definition: transext.h:93
int p_Compare(const poly a, const poly b, const ring R)
Definition: p_polys.cc:4760
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:1463
#define pOne()
Definition: polys.h:297
static int isMonoIdBasesSame_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
Definition: hilb.cc:1504
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:789
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
static int monCompare(const void *m, const void *n)
Definition: hilb.cc:1648
ideal idCopy(ideal A)
Definition: ideals.h:60
void rChangeCurrRing(ring r)
Definition: polys.cc:12
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:44
static int CountOnIdUptoTruncationIndex(ideal I, int tr)
Definition: hilb.cc:1477
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:1419
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 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:365
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:3602
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:1317
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:758
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:228
static ideal SortByDeg_p(ideal I, poly p)
!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
Definition: hilb.cc:265
static poly ChooseP(ideal I)
Definition: hilb.cc:742
p exp[i]
Definition: DebugPrint.cc:39
int hisModule
Definition: hutil.cc:23
intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1309
matrix mp_Sub(matrix a, matrix b, const ring R)
Definition: matpol.cc:203
void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1388
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:857
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1298
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:1368
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:815