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