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