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