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