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