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