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  }
466  id_Delete(&res,currRing);
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);
500  p_Setm(res,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);
798  p_Setm(res,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 
1833  p_SetComp(smon, p_GetComp(p, currRing), currRing);
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);
2168  coeffs cf = nInitChar(n_transExt, &p);
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 
int status int void size_t count
Definition: si_signals.h:59
ideal idQuotMon(ideal Iorig, ideal p)
Definition: hilb.cc:410
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
#define id_TestTail(A, lR, tR)
Definition: simpleideals.h:79
static int variables
int hNstc
Definition: hutil.cc:22
const CanonicalForm int s
Definition: facAbsFact.cc:55
void hLexS(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:512
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, const matrix bVec, matrix &xVec, matrix &H)
Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decompositi...
int j
Definition: facHensel.cc:105
#define nWrite(n)
Definition: numbers.h:30
void PrintLn()
Definition: reporter.cc:310
static int hLength
Definition: hilb.cc:47
static poly ChoosePVar(ideal I)
Definition: hilb.cc:480
#define Print
Definition: emacs.cc:80
scfmon hwork
Definition: hutil.cc:19
void mu(int **points, int sizePoints)
scfmon hGetmem(int lm, scfmon old, monp monmem)
Definition: hutil.cc:1029
int hNexist
Definition: hutil.cc:22
int * varset
Definition: hutil.h:16
void hElimS(scfmon stc, int *e1, int a2, int e2, varset var, int Nvar)
Definition: hutil.cc:678
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:358
void scPrintDegree(int co, int mu)
Definition: hdegree.cc:808
static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var, int Nvar, int *pol, int Lpol)
Definition: hilb.cc:178
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:94
void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp, int trunDegHs)
Definition: hilb.cc:1969
ideal id_Copy(ideal h1, const ring r)
copy an ideal
scmon hGetpure(scmon p)
Definition: hutil.cc:1058
static int positionInOrbitTruncationCase(ideal I, poly w, std::vector< ideal > idorb, std::vector< poly > polist, int, int trunDegHs)
Definition: hilb.cc:1680
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:247
void slicehilb(ideal I)
Definition: hilb.cc:1131
static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
Definition: hilb.cc:138
scmon * scfmon
Definition: hutil.h:15
#define p_GetComp(p, r)
Definition: monomials.h:71
static ideal minimalMonomialGenSet(ideal I)
Definition: hilb.cc:1774
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1466
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
int rows() const
Definition: intvec.h:94
scfmon hexist
Definition: hutil.cc:19
static void idAddMon(ideal I, ideal p)
Definition: hilb.cc:472
static int * Q0
Definition: hilb.cc:46
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:583
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
monf hCreate(int Nvar)
Definition: hutil.cc:1002
static bool JustVar(ideal I)
Definition: hilb.cc:807
long int64
Definition: auxiliary.h:66
int hNvar
Definition: hutil.cc:22
bool unitMatrix(const int n, matrix &unitMat, const ring R)
Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.
static bool IsIn(poly p, ideal I)
Definition: hilb.cc:910
static poly pp_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:987
#define TRUE
Definition: auxiliary.h:98
static long p_Totaldegree(poly p, const ring r)
Definition: p_polys.h:1453
void * ADDRESS
Definition: auxiliary.h:133
int hNpure
Definition: hutil.cc:22
void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
Definition: hilb.cc:1381
void pWrite(poly p)
Definition: polys.h:294
const int MAX_INT_VAL
Definition: mylimits.h:12
scmon hpure
Definition: hutil.cc:20
void WerrorS(const char *s)
Definition: feFopen.cc:24
static void TwordMap(poly p, poly w, int lV, int d, ideal Jwi, bool &flag)
Definition: hilb.cc:1873
int k
Definition: cfEzgcd.cc:92
char * StringEndS()
Definition: reporter.cc:151
#define Q
Definition: sirandom.c:25
#define pLmDivisibleBy(a, b)
like pDivisibleBy, except that it is assumed that a!=NULL, b!=NULL
Definition: polys.h:140
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:51
#define loop
Definition: structs.h:78
#define omAlloc(size)
Definition: omAllocDecl.h:210
void pWrite0(poly p)
Definition: polys.h:295
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:412
static poly LCMmon(ideal I)
Definition: hilb.cc:949
intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1329
static void p_SetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1481
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:812
void prune(Variable &alpha)
Definition: variable.cc:261
static FORCE_INLINE number n_Param(const int iParameter, const coeffs r)
return the (iParameter^th) parameter as a NEW number NOTE: parameter numbering: 1..n_NumberOfParameters(...)
Definition: coeffs.h:814
void hDelete(scfmon ev, int ev_length)
Definition: hutil.cc:146
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:637
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)
Definition: p_polys.h:825
#define idPrint(id)
Definition: ideals.h:46
Coefficient rings, fields and other domains suitable for Singular polynomials.
void hOrdSupp(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:208
Definition: intvec.h:17
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:1305
void hKill(monf xmem, int Nvar)
Definition: hutil.cc:1016
Variable var
Definition: int_poly.h:74
void rKill(ring r)
Definition: ipshell.cc:6076
varset hvar
Definition: hutil.cc:21
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent : the integer VarOffset encodes:
Definition: p_polys.h:469
#define omFree(addr)
Definition: omAllocDecl.h:261
#define assume(x)
Definition: mod2.h:390
static poly ChoosePJL(ideal I)
Definition: hilb.cc:706
static int * hAddHilb(int Nv, int x, int *pol, int *lp)
Definition: hilb.cc:105
The main handler for Singular numbers which are suitable for Singular polynomials.
static void hWDegree(intvec *wdegree)
Definition: hilb.cc:246
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:1649
void hPure(scfmon stc, int a, int *Nstc, varset var, int Nvar, scmon pure, int *Npure)
Definition: hutil.cc:627
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
void hStaircase(scfmon stc, int *Nstc, varset var, int Nvar)
Definition: hutil.cc:319
void sortMonoIdeal_pCompare(ideal I)
Definition: hilb.cc:1761
static BOOLEAN p_DivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1824
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:1478
int m
Definition: cfEzgcd.cc:121
static intvec * hSeries(ideal S, intvec *modulweight, int, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1173
int * scmon
Definition: hutil.h:14
struct for passing initialization parameters to naInitChar
Definition: transext.h:88
int p_Compare(const poly a, const poly b, const ring R)
Definition: p_polys.cc:4790
static BOOLEAN p_IsOne(const poly p, const ring R)
either poly(1) or gen(k)?!
Definition: p_polys.h:1931
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:818
void PrintS(const char *s)
Definition: reporter.cc:284
static void SortByDeg_p(ideal I, poly p)
!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
Definition: hilb.cc:289
#define pOne()
Definition: polys.h:301
#define IDELEMS(i)
Definition: simpleideals.h:24
static BOOLEAN p_LmDivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1815
void mp_Delete(matrix *a, const ring r)
Definition: matpol.cc:780
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:103
static int monCompare(const void *m, const void *n)
Definition: hilb.cc:1754
short errorreported
Definition: feFopen.cc:23
void rChangeCurrRing(ring r)
Definition: polys.cc:15
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:36
static int CountOnIdUptoTruncationIndex(ideal I, int tr)
Definition: hilb.cc:1509
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:857
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:37
#define p_LmEqual(p1, p2, r)
Definition: p_polys.h:1658
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:488
ideal RightColonOperation(ideal S, poly w, int lV)
Definition: hilb.cc:2309
static ideal SortByDeg(ideal I)
Definition: hilb.cc:389
CanonicalForm cf
Definition: cfModGcd.cc:4024
#define omalloc(size)
Definition: omAllocDecl.h:228
#define NULL
Definition: omList.c:10
static int * Ql
Definition: hilb.cc:46
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3633
monf stcmem
Definition: hutil.cc:24
int length() const
Definition: intvec.h:92
void hStepS(scfmon stc, int Nstc, varset var, int Nvar, int *a, int *x)
Definition: hutil.cc:955
intvec * hSecondSeries(intvec *hseries1)
Definition: hilb.cc:1346
#define R
Definition: sirandom.c:26
static poly shiftInMon(poly p, int i, int lV, const ring r)
Definition: hilb.cc:1806
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:173
Variable x
Definition: cfModGcd.cc:4023
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
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
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:233
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
poly p_MDivide(poly a, poly b, const ring r)
Definition: p_polys.cc:1454
static poly ChooseP(ideal I)
Definition: hilb.cc:765
int hisModule
Definition: hutil.cc:23
int p
Definition: cfModGcd.cc:4019
intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1336
matrix mp_Sub(matrix a, matrix b, const ring R)
Definition: matpol.cc:195
void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1419
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:892
scfmon hstc
Definition: hutil.cc:19
static int hMinModulweight(intvec *modulweight)
Definition: hilb.cc:50
static int comapreMonoIdBases_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
Definition: hilb.cc:1536
static ideal colonIdeal(ideal S, poly w, int lV, ideal Jwi, int trunDegHs)
Definition: hilb.cc:1934
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:160
static poly SqFree(ideal I)
Definition: hilb.cc:881
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1289
scfmon hInit(ideal S, ideal Q, int *Nexist, ring tailRing)
Definition: hutil.cc:34
#define pLmEqual(p1, p2)
Definition: polys.h:111
static poly deleteInMon(poly w, int i, int lV, const ring r)
Definition: hilb.cc:1839
#define omAlloc0(size)
Definition: omAllocDecl.h:211
int l
Definition: cfEzgcd.cc:93
static void idInsertMonomial(ideal I, poly p)
Definition: hilb.cc:1451
static int ** Qpol
Definition: hilb.cc:45
static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hilb.cc:64
#define pCopy(p)
return a copy of the poly
Definition: polys.h:172
#define MATELEM(mat, i, j)
Definition: matpol.h:28
static void hPrintHilb(intvec *hseries)
Definition: hilb.cc:1399
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:180
#define omStrDup(s)
Definition: omAllocDecl.h:263
static void eulerchar(ideal I, int variables, mpz_ptr ec)
Definition: hilb.cc:838