My Project
Loading...
Searching...
No Matches
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
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
67static 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
167static 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 }
183 //idDegSortTest(res);
184 return(res);
185}
186
187//idQuot(I,p) for I monomial ideal, p a ideal with a single monomial.
188ideal 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);
226 {
227 p_Delete(&res->m[i],currRing);
228 }
229 else
230 {
231 p_Delete(&I->m[i],currRing);
232 }
233 }
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
250static 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)
258static 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 {
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)
290static 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);
305 flag = FALSE;
306 }
307 }
308 if(!p_IsOne(m, currRing))
309 {
310 return(m);
311 }
312 }
314 m = ChoosePVar(I);
315 return(m);
316}
317
318//chooses 1 \neq p \not\in S. This choice should be made optimal
319static 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)
327static 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)
353static 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
363static 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);
402 }
403}
404
405//tests if an ideal is Square Free, if no, returns the variable which appears at powers >1
406static 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
435static 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
474static 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 }
496 return(m);
497}
498
499//the Roune Slice Algorithm
500static 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 }
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 }
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);
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);
585 break;
586 }
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
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
656void 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);
675 I = Itmp;
676 //printf("\n-------------RouneSlice--------------\n");
677 rouneslice(I,S,q,X->m[0],prune, moreprune, steps, NNN, hilbertcoef, hilbpower);
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
732void 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
750static 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*/
776void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree)
777{
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))
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
808static 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
835static 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
866static 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
893static 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
928static 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
1006static 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
1037static 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
1111static int monCompare( const void *m, const void *n)
1112{
1113 /* compares monomials */
1114
1115 return(p_Compare(*(poly*) m, *(poly*)n, currRing));
1116}
1117
1118static 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
1131static 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
1163static 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
1196static 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
1230static 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
1291static 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
1326void 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);
1369 }
1370 else
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;
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
1673ideal 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"
1702#include "polys/simpleideals.h"
1703#include "coeffs/coeffs.h"
1704
1705#include "kernel/ideals.h"
1706
1707static BOOLEAN p_Div_hi(poly p, const int* exp_q, const ring src)
1708{
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__)
1729static int compare_rp(void *arg,const void *pp1, const void *pp2)
1730#else
1731static 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
1747static 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
1761static 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}
1830poly 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
1929static 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
1952intvec* 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
2009static 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
2020static 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
2036intvec* 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/* ------------------------------------------------------------------ */
2071static int hMinModulweight(intvec *modulweight)
2072{
2073 if(modulweight==NULL) return 0;
2074 return modulweight->min_in();
2075}
2076
2077static 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}
2094static 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
2142static 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
2198static 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
2239static 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
2304static 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);
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
2450intvec * 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
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
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
intvec * hSecondSeries(intvec *hseries1)
Definition: hilb.cc:697
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
static intvec * hSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q)
Definition: hilb.cc:2304
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
intvec * hFirstSeries(ideal A, intvec *module_w, ideal Q, intvec *wdegree)
Definition: hilb.cc:2036
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
intvec * hFirstSeries0(ideal A, ideal Q, intvec *wdegree, const ring src, const ring Qt)
Definition: hilb.cc:1952
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
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 int64 * hAddHilb(int Nv, int x, int64 *pol, int *lp)
Definition: hilb.cc:2094
static ideal minimalMonomialGenSet(ideal I)
Definition: hilb.cc:1131
intvec * hFirstSeries1(ideal S, intvec *modulweight, ideal Q, intvec *wdegree)
Definition: hilb.cc:2450
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
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
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
#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:873
static matrix mu(matrix A, const ring R)
Definition: matpol.cc:2025
matrix mp_Sub(matrix a, matrix b, const ring R)
Definition: matpol.cc:189
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:2197
poly p_MDivide(poly a, poly b, const ring r)
Definition: p_polys.cc:1492
int p_Compare(const poly a, const poly b, const ring R)
Definition: p_polys.cc:4849
unsigned long p_GetShortExpVector(const poly p, const ring r)
Definition: p_polys.cc:4780
poly p_One(const ring r)
Definition: p_polys.cc:1313
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3696
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
#define Q
Definition: sirandom.c:26
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
void prune(Variable &alpha)
Definition: variable.cc:261