1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /*
5 * ABSTRACT - all basic methods to manipulate ideals
6 */
7
8
9 /* includes */
10
11
12
13 #include "misc/auxiliary.h"
14
15 #include "misc/options.h"
16 #include "misc/intvec.h"
17
18 #include "matpol.h"
19
20 #include "monomials/p_polys.h"
21 #include "weight.h"
22 #include "sbuckets.h"
23 #include "clapsing.h"
24
25 #include "simpleideals.h"
26
27 VAR omBin sip_sideal_bin = omGetSpecBin(sizeof(sip_sideal));
28
29 STATIC_VAR poly * idpower;
30 /*collects the monomials in makemonoms, must be allocated befor*/
31 STATIC_VAR int idpowerpoint;
32 /*index of the actual monomial in idpower*/
33
34 /// initialise an ideal / module
idInit(int idsize,int rank)35 ideal idInit(int idsize, int rank)
36 {
37 assume( idsize >= 0 && rank >= 0 );
38
39 ideal hh = (ideal)omAllocBin(sip_sideal_bin);
40
41 IDELEMS(hh) = idsize; // ncols
42 hh->nrows = 1; // ideal/module!
43
44 hh->rank = rank; // ideal: 1, module: >= 0!
45
46 if (idsize>0)
47 hh->m = (poly *)omAlloc0(idsize*sizeof(poly));
48 else
49 hh->m = NULL;
50
51 return hh;
52 }
53
54 #ifdef PDEBUG
55 // this is only for outputting an ideal within the debugger
56 // therefor it accept the otherwise illegal id==NULL
idShow(const ideal id,const ring lmRing,const ring tailRing,const int debugPrint)57 void idShow(const ideal id, const ring lmRing, const ring tailRing, const int debugPrint)
58 {
59 assume( debugPrint >= 0 );
60
61 if( id == NULL )
62 PrintS("(NULL)");
63 else
64 {
65 Print("Module of rank %ld,real rank %ld and %d generators.\n",
66 id->rank,id_RankFreeModule(id, lmRing, tailRing),IDELEMS(id));
67
68 int j = (id->ncols*id->nrows) - 1;
69 while ((j > 0) && (id->m[j]==NULL)) j--;
70 for (int i = 0; i <= j; i++)
71 {
72 Print("generator %d: ",i); p_wrp(id->m[i], lmRing, tailRing);PrintLn();
73 }
74 }
75 }
76 #endif
77
78 /// index of generator with leading term in ground ring (if any);
79 /// otherwise -1
id_PosConstant(ideal id,const ring r)80 int id_PosConstant(ideal id, const ring r)
81 {
82 id_Test(id, r);
83 const int N = IDELEMS(id) - 1;
84 const poly * m = id->m + N;
85
86 for (int k = N; k >= 0; --k, --m)
87 {
88 const poly p = *m;
89 if (p!=NULL)
90 if (p_LmIsConstantComp(p, r) == TRUE)
91 return k;
92 }
93
94 return -1;
95 }
96
97 /// initialise the maximal ideal (at 0)
id_MaxIdeal(const ring r)98 ideal id_MaxIdeal (const ring r)
99 {
100 int nvars;
101 #ifdef HAVE_SHIFTBBA
102 if (r->isLPring)
103 {
104 nvars = r->isLPring;
105 }
106 else
107 #endif
108 {
109 nvars = rVar(r);
110 }
111 ideal hh = idInit(nvars, 1);
112 for (int l=nvars-1; l>=0; l--)
113 {
114 hh->m[l] = p_One(r);
115 p_SetExp(hh->m[l],l+1,1,r);
116 p_Setm(hh->m[l],r);
117 }
118 id_Test(hh, r);
119 return hh;
120 }
121
122 /// deletes an ideal/module/matrix
id_Delete(ideal * h,ring r)123 void id_Delete (ideal * h, ring r)
124 {
125 if (*h == NULL)
126 return;
127
128 id_Test(*h, r);
129
130 const long elems = (long)(*h)->nrows * (long)(*h)->ncols;
131
132 if ( elems > 0 )
133 {
134 assume( (*h)->m != NULL );
135
136 if (r!=NULL)
137 {
138 long j = elems;
139 do
140 {
141 j--;
142 poly pp=((*h)->m[j]);
143 if (pp!=NULL) p_Delete(&pp, r);
144 }
145 while (j>0);
146 }
147
148 omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*elems);
149 }
150
151 omFreeBin((ADDRESS)*h, sip_sideal_bin);
152 *h=NULL;
153 }
154
155
156 /// Shallowdeletes an ideal/matrix
id_ShallowDelete(ideal * h,ring r)157 void id_ShallowDelete (ideal *h, ring r)
158 {
159 id_Test(*h, r);
160
161 if (*h == NULL)
162 return;
163
164 int j,elems;
165 elems=j=(*h)->nrows*(*h)->ncols;
166 if (j>0)
167 {
168 assume( (*h)->m != NULL );
169 do
170 {
171 p_ShallowDelete(&((*h)->m[--j]), r);
172 }
173 while (j>0);
174 omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*elems);
175 }
176 omFreeBin((ADDRESS)*h, sip_sideal_bin);
177 *h=NULL;
178 }
179
180 /// gives an ideal/module the minimal possible size
idSkipZeroes(ideal ide)181 void idSkipZeroes (ideal ide)
182 {
183 assume (ide != NULL);
184
185 int k;
186 int j = -1;
187 BOOLEAN change=FALSE;
188
189 for (k=0; k<IDELEMS(ide); k++)
190 {
191 if (ide->m[k] != NULL)
192 {
193 j++;
194 if (change)
195 {
196 ide->m[j] = ide->m[k];
197 }
198 }
199 else
200 {
201 change=TRUE;
202 }
203 }
204 if (change)
205 {
206 if (j == -1)
207 j = 0;
208 else
209 {
210 for (k=j+1; k<IDELEMS(ide); k++)
211 ide->m[k] = NULL;
212 }
213 pEnlargeSet(&(ide->m),IDELEMS(ide),j+1-IDELEMS(ide));
214 IDELEMS(ide) = j+1;
215 }
216 }
217
218 /// count non-zero elements
idElem(const ideal F)219 int idElem(const ideal F)
220 {
221 assume (F != NULL);
222
223 int i=0;
224
225 for(int j=IDELEMS(F)-1;j>=0;j--)
226 {
227 if ((F->m)[j]!=NULL) i++;
228 }
229 return i;
230 }
231
232 /// copies the first k (>= 1) entries of the given ideal/module
233 /// and returns these as a new ideal/module
234 /// (Note that the copied entries may be zero.)
id_CopyFirstK(const ideal ide,const int k,const ring r)235 ideal id_CopyFirstK (const ideal ide, const int k,const ring r)
236 {
237 id_Test(ide, r);
238
239 assume( ide != NULL );
240 assume( k <= IDELEMS(ide) );
241
242 ideal newI = idInit(k, ide->rank);
243
244 for (int i = 0; i < k; i++)
245 newI->m[i] = p_Copy(ide->m[i],r);
246
247 return newI;
248 }
249
250 /// ideal id = (id[i]), result is leadcoeff(id[i]) = 1
id_Norm(ideal id,const ring r)251 void id_Norm(ideal id, const ring r)
252 {
253 id_Test(id, r);
254 for (int i=IDELEMS(id)-1; i>=0; i--)
255 {
256 if (id->m[i] != NULL)
257 {
258 p_Norm(id->m[i],r);
259 }
260 }
261 }
262
263 /// ideal id = (id[i]), c any unit
264 /// if id[i] = c*id[j] then id[j] is deleted for j > i
id_DelMultiples(ideal id,const ring r)265 void id_DelMultiples(ideal id, const ring r)
266 {
267 id_Test(id, r);
268
269 int i, j;
270 int k = IDELEMS(id)-1;
271 for (i=k; i>=0; i--)
272 {
273 if (id->m[i]!=NULL)
274 {
275 for (j=k; j>i; j--)
276 {
277 if (id->m[j]!=NULL)
278 {
279 if (rField_is_Ring(r))
280 {
281 /* if id[j] = c*id[i] then delete id[j].
282 In the below cases of a ground field, we
283 check whether id[i] = c*id[j] and, if so,
284 delete id[j] for historical reasons (so
285 that previous output does not change) */
286 if (p_ComparePolys(id->m[j], id->m[i],r)) p_Delete(&id->m[j],r);
287 }
288 else
289 {
290 if (p_ComparePolys(id->m[i], id->m[j],r)) p_Delete(&id->m[j],r);
291 }
292 }
293 }
294 }
295 }
296 }
297
298 /// ideal id = (id[i])
299 /// if id[i] = id[j] then id[j] is deleted for j > i
id_DelEquals(ideal id,const ring r)300 void id_DelEquals(ideal id, const ring r)
301 {
302 id_Test(id, r);
303
304 int i, j;
305 int k = IDELEMS(id)-1;
306 for (i=k; i>=0; i--)
307 {
308 if (id->m[i]!=NULL)
309 {
310 for (j=k; j>i; j--)
311 {
312 if ((id->m[j]!=NULL)
313 && (p_EqualPolys(id->m[i], id->m[j],r)))
314 {
315 p_Delete(&id->m[j],r);
316 }
317 }
318 }
319 }
320 }
321
322 /// Delete id[j], if Lm(j) == Lm(i) and both LC(j), LC(i) are units and j > i
id_DelLmEquals(ideal id,const ring r)323 void id_DelLmEquals(ideal id, const ring r)
324 {
325 id_Test(id, r);
326
327 int i, j;
328 int k = IDELEMS(id)-1;
329 for (i=k; i>=0; i--)
330 {
331 if (id->m[i] != NULL)
332 {
333 for (j=k; j>i; j--)
334 {
335 if ((id->m[j] != NULL)
336 && p_LmEqual(id->m[i], id->m[j],r)
337 #ifdef HAVE_RINGS
338 && n_IsUnit(pGetCoeff(id->m[i]),r->cf) && n_IsUnit(pGetCoeff(id->m[j]),r->cf)
339 #endif
340 )
341 {
342 p_Delete(&id->m[j],r);
343 }
344 }
345 }
346 }
347 }
348
349 /// delete id[j], if LT(j) == coeff*mon*LT(i) and vice versa, i.e.,
350 /// delete id[i], if LT(i) == coeff*mon*LT(j)
id_DelDiv(ideal id,const ring r)351 void id_DelDiv(ideal id, const ring r)
352 {
353 id_Test(id, r);
354
355 int i, j;
356 int k = IDELEMS(id)-1;
357 for (i=k; i>=0; i--)
358 {
359 if (id->m[i] != NULL)
360 {
361 for (j=k; j>i; j--)
362 {
363 if (id->m[j]!=NULL)
364 {
365 #ifdef HAVE_RINGS
366 if (rField_is_Ring(r))
367 {
368 if (p_DivisibleByRingCase(id->m[i], id->m[j],r))
369 {
370 p_Delete(&id->m[j],r);
371 }
372 else if (p_DivisibleByRingCase(id->m[j], id->m[i],r))
373 {
374 p_Delete(&id->m[i],r);
375 break;
376 }
377 }
378 else
379 #endif
380 {
381 /* the case of a coefficient field: */
382 if (p_DivisibleBy(id->m[i], id->m[j],r))
383 {
384 p_Delete(&id->m[j],r);
385 }
386 else if (p_DivisibleBy(id->m[j], id->m[i],r))
387 {
388 p_Delete(&id->m[i],r);
389 break;
390 }
391 }
392 }
393 }
394 }
395 }
396 }
397
398 /// test if the ideal has only constant polynomials
399 /// NOTE: zero ideal/module is also constant
id_IsConstant(ideal id,const ring r)400 BOOLEAN id_IsConstant(ideal id, const ring r)
401 {
402 id_Test(id, r);
403
404 for (int k = IDELEMS(id)-1; k>=0; k--)
405 {
406 if (!p_IsConstantPoly(id->m[k],r))
407 return FALSE;
408 }
409 return TRUE;
410 }
411
412 /// copy an ideal
id_Copy(ideal h1,const ring r)413 ideal id_Copy(ideal h1, const ring r)
414 {
415 id_Test(h1, r);
416
417 ideal h2 = idInit(IDELEMS(h1), h1->rank);
418 for (int i=IDELEMS(h1)-1; i>=0; i--)
419 h2->m[i] = p_Copy(h1->m[i],r);
420 return h2;
421 }
422
423 #ifdef PDEBUG
424 /// Internal verification for ideals/modules and dense matrices!
id_DBTest(ideal h1,int level,const char * f,const int l,const ring r,const ring tailRing)425 void id_DBTest(ideal h1, int level, const char *f,const int l, const ring r, const ring tailRing)
426 {
427 if (h1 != NULL)
428 {
429 // assume(IDELEMS(h1) > 0); for ideal/module, does not apply to matrix
430 omCheckAddrSize(h1,sizeof(*h1));
431
432 assume( h1->ncols >= 0 );
433 assume( h1->nrows >= 0 ); // matrix case!
434
435 assume( h1->rank >= 0 );
436
437 const long n = ((long)h1->ncols * (long)h1->nrows);
438
439 assume( !( n > 0 && h1->m == NULL) );
440
441 if( h1->m != NULL && n > 0 )
442 omdebugAddrSize(h1->m, n * sizeof(poly));
443
444 long new_rk = 0; // inlining id_RankFreeModule(h1, r, tailRing);
445
446 /* to be able to test matrices: */
447 for (long i=n - 1; i >= 0; i--)
448 {
449 _pp_Test(h1->m[i], r, tailRing, level);
450 const long k = p_MaxComp(h1->m[i], r, tailRing);
451 if (k > new_rk) new_rk = k;
452 }
453
454 // dense matrices only contain polynomials:
455 // h1->nrows == h1->rank > 1 && new_rk == 0!
456 assume( !( h1->nrows == h1->rank && h1->nrows > 1 && new_rk > 0 ) ); //
457
458 if(new_rk > h1->rank)
459 {
460 dReportError("wrong rank %d (should be %d) in %s:%d\n",
461 h1->rank, new_rk, f,l);
462 omPrintAddrInfo(stderr, h1, " for ideal");
463 h1->rank = new_rk;
464 }
465 }
466 else
467 {
468 Print("error: ideal==NULL in %s:%d\n",f,l);
469 assume( h1 != NULL );
470 }
471 }
472 #endif
473
474 /// for idSort: compare a and b revlex inclusive module comp.
p_Comp_RevLex(poly a,poly b,BOOLEAN nolex,const ring R)475 static int p_Comp_RevLex(poly a, poly b,BOOLEAN nolex, const ring R)
476 {
477 if (b==NULL) return 1;
478 if (a==NULL) return -1;
479
480 if (nolex)
481 {
482 int r=p_LtCmp(a,b,R);
483 return r;
484 #if 0
485 if (r!=0) return r;
486 number h=n_Sub(pGetCoeff(a),pGetCoeff(b),R->cf);
487 r = -1+n_IsZero(h,R->cf)+2*n_GreaterZero(h,R->cf); /* -1: <, 0:==, 1: > */
488 n_Delete(&h, R->cf);
489 return r;
490 #endif
491 }
492 int l=rVar(R);
493 while ((l>0) && (p_GetExp(a,l,R)==p_GetExp(b,l,R))) l--;
494 if (l==0)
495 {
496 if (p_GetComp(a,R)==p_GetComp(b,R))
497 {
498 number h=n_Sub(pGetCoeff(a),pGetCoeff(b),R->cf);
499 int r = -1+n_IsZero(h,R->cf)+2*n_GreaterZero(h,R->cf); /* -1: <, 0:==, 1: > */
500 n_Delete(&h,R->cf);
501 return r;
502 }
503 if (p_GetComp(a,R)>p_GetComp(b,R)) return 1;
504 }
505 else if (p_GetExp(a,l,R)>p_GetExp(b,l,R))
506 return 1;
507 return -1;
508 }
509
510 // sorts the ideal w.r.t. the actual ringordering
511 // uses lex-ordering when nolex = FALSE
id_Sort(const ideal id,const BOOLEAN nolex,const ring r)512 intvec *id_Sort(const ideal id, const BOOLEAN nolex, const ring r)
513 {
514 id_Test(id, r);
515
516 intvec * result = new intvec(IDELEMS(id));
517 int i, j, actpos=0, newpos;
518 int diff, olddiff, lastcomp, newcomp;
519 BOOLEAN notFound;
520
521 for (i=0;i<IDELEMS(id);i++)
522 {
523 if (id->m[i]!=NULL)
524 {
525 notFound = TRUE;
526 newpos = actpos / 2;
527 diff = (actpos+1) / 2;
528 diff = (diff+1) / 2;
529 lastcomp = p_Comp_RevLex(id->m[i],id->m[(*result)[newpos]],nolex,r);
530 if (lastcomp<0)
531 {
532 newpos -= diff;
533 }
534 else if (lastcomp>0)
535 {
536 newpos += diff;
537 }
538 else
539 {
540 notFound = FALSE;
541 }
542 //while ((newpos>=0) && (newpos<actpos) && (notFound))
543 while (notFound && (newpos>=0) && (newpos<actpos))
544 {
545 newcomp = p_Comp_RevLex(id->m[i],id->m[(*result)[newpos]],nolex,r);
546 olddiff = diff;
547 if (diff>1)
548 {
549 diff = (diff+1) / 2;
550 if ((newcomp==1)
551 && (actpos-newpos>1)
552 && (diff>1)
553 && (newpos+diff>=actpos))
554 {
555 diff = actpos-newpos-1;
556 }
557 else if ((newcomp==-1)
558 && (diff>1)
559 && (newpos<diff))
560 {
561 diff = newpos;
562 }
563 }
564 if (newcomp<0)
565 {
566 if ((olddiff==1) && (lastcomp>0))
567 notFound = FALSE;
568 else
569 newpos -= diff;
570 }
571 else if (newcomp>0)
572 {
573 if ((olddiff==1) && (lastcomp<0))
574 {
575 notFound = FALSE;
576 newpos++;
577 }
578 else
579 {
580 newpos += diff;
581 }
582 }
583 else
584 {
585 notFound = FALSE;
586 }
587 lastcomp = newcomp;
588 if (diff==0) notFound=FALSE; /*hs*/
589 }
590 if (newpos<0) newpos = 0;
591 if (newpos>actpos) newpos = actpos;
592 while ((newpos<actpos) && (p_Comp_RevLex(id->m[i],id->m[(*result)[newpos]],nolex,r)==0))
593 newpos++;
594 for (j=actpos;j>newpos;j--)
595 {
596 (*result)[j] = (*result)[j-1];
597 }
598 (*result)[newpos] = i;
599 actpos++;
600 }
601 }
602 for (j=0;j<actpos;j++) (*result)[j]++;
603 return result;
604 }
605
606 /// concat the lists h1 and h2 without zeros
id_SimpleAdd(ideal h1,ideal h2,const ring R)607 ideal id_SimpleAdd (ideal h1,ideal h2, const ring R)
608 {
609 id_Test(h1, R);
610 id_Test(h2, R);
611
612 if ( idIs0(h1) )
613 {
614 ideal res=id_Copy(h2,R);
615 if (res->rank<h1->rank) res->rank=h1->rank;
616 return res;
617 }
618 if ( idIs0(h2) )
619 {
620 ideal res=id_Copy(h1,R);
621 if (res->rank<h2->rank) res->rank=h2->rank;
622 return res;
623 }
624
625 int j = IDELEMS(h1)-1;
626 while ((j >= 0) && (h1->m[j] == NULL)) j--;
627
628 int i = IDELEMS(h2)-1;
629 while ((i >= 0) && (h2->m[i] == NULL)) i--;
630
631 const int r = si_max(h1->rank, h2->rank);
632
633 ideal result = idInit(i+j+2,r);
634
635 int l;
636
637 for (l=j; l>=0; l--)
638 result->m[l] = p_Copy(h1->m[l],R);
639
640 j = i+j+1;
641 for (l=i; l>=0; l--, j--)
642 result->m[j] = p_Copy(h2->m[l],R);
643
644 return result;
645 }
646
647 /// insert h2 into h1 (if h2 is not the zero polynomial)
648 /// return TRUE iff h2 was indeed inserted
idInsertPoly(ideal h1,poly h2)649 BOOLEAN idInsertPoly (ideal h1, poly h2)
650 {
651 if (h2==NULL) return FALSE;
652 assume (h1 != NULL);
653
654 int j = IDELEMS(h1) - 1;
655
656 while ((j >= 0) && (h1->m[j] == NULL)) j--;
657 j++;
658 if (j==IDELEMS(h1))
659 {
660 pEnlargeSet(&(h1->m),IDELEMS(h1),16);
661 IDELEMS(h1)+=16;
662 }
663 h1->m[j]=h2;
664 return TRUE;
665 }
666
667 /// insert p into I on position pos
idInsertPolyOnPos(ideal I,poly p,int pos)668 BOOLEAN idInsertPolyOnPos (ideal I, poly p,int pos)
669 {
670 if (p==NULL) return FALSE;
671 assume (I != NULL);
672
673 int j = IDELEMS(I) - 1;
674
675 while ((j >= 0) && (I->m[j] == NULL)) j--;
676 j++;
677 if (j==IDELEMS(I))
678 {
679 pEnlargeSet(&(I->m),IDELEMS(I),IDELEMS(I)+1);
680 IDELEMS(I)+=1;
681 }
682 for(j = IDELEMS(I)-1;j>pos;j--)
683 I->m[j] = I->m[j-1];
684 I->m[pos]=p;
685 return TRUE;
686 }
687
688
689 /*! insert h2 into h1 depending on the two boolean parameters:
690 * - if zeroOk is true, then h2 will also be inserted when it is zero
691 * - if duplicateOk is true, then h2 will also be inserted when it is
692 * already present in h1
693 * return TRUE iff h2 was indeed inserted
694 */
id_InsertPolyWithTests(ideal h1,const int validEntries,const poly h2,const bool zeroOk,const bool duplicateOk,const ring r)695 BOOLEAN id_InsertPolyWithTests (ideal h1, const int validEntries,
696 const poly h2, const bool zeroOk, const bool duplicateOk, const ring r)
697 {
698 id_Test(h1, r);
699 p_Test(h2, r);
700
701 if ((!zeroOk) && (h2 == NULL)) return FALSE;
702 if (!duplicateOk)
703 {
704 bool h2FoundInH1 = false;
705 int i = 0;
706 while ((i < validEntries) && (!h2FoundInH1))
707 {
708 h2FoundInH1 = p_EqualPolys(h1->m[i], h2,r);
709 i++;
710 }
711 if (h2FoundInH1) return FALSE;
712 }
713 if (validEntries == IDELEMS(h1))
714 {
715 pEnlargeSet(&(h1->m), IDELEMS(h1), 16);
716 IDELEMS(h1) += 16;
717 }
718 h1->m[validEntries] = h2;
719 return TRUE;
720 }
721
722 /// h1 + h2
id_Add(ideal h1,ideal h2,const ring r)723 ideal id_Add (ideal h1,ideal h2, const ring r)
724 {
725 id_Test(h1, r);
726 id_Test(h2, r);
727
728 ideal result = id_SimpleAdd(h1,h2,r);
729 id_Compactify(result,r);
730 return result;
731 }
732
733 /// h1 * h2
734 /// one h_i must be an ideal (with at least one column)
735 /// the other h_i may be a module (with no columns at all)
id_Mult(ideal h1,ideal h2,const ring R)736 ideal id_Mult (ideal h1,ideal h2, const ring R)
737 {
738 id_Test(h1, R);
739 id_Test(h2, R);
740
741 int j = IDELEMS(h1);
742 while ((j > 0) && (h1->m[j-1] == NULL)) j--;
743
744 int i = IDELEMS(h2);
745 while ((i > 0) && (h2->m[i-1] == NULL)) i--;
746
747 j *= i;
748 int r = si_max( h2->rank, h1->rank );
749 if (j==0)
750 {
751 if ((IDELEMS(h1)>0) && (IDELEMS(h2)>0)) j=1;
752 return idInit(j, r);
753 }
754 ideal hh = idInit(j, r);
755
756 int k = 0;
757 for (i=0; i<IDELEMS(h1); i++)
758 {
759 if (h1->m[i] != NULL)
760 {
761 for (j=0; j<IDELEMS(h2); j++)
762 {
763 if (h2->m[j] != NULL)
764 {
765 hh->m[k] = pp_Mult_qq(h1->m[i],h2->m[j],R);
766 k++;
767 }
768 }
769 }
770 }
771
772 id_Compactify(hh,R);
773 return hh;
774 }
775
776 /// returns true if h is the zero ideal
idIs0(ideal h)777 BOOLEAN idIs0 (ideal h)
778 {
779 assume (h != NULL); // will fail :(
780 // if (h == NULL) return TRUE;
781
782 for( int i = IDELEMS(h)-1; i >= 0; i-- )
783 if(h->m[i] != NULL)
784 return FALSE;
785
786 return TRUE;
787
788 }
789
790 /// return the maximal component number found in any polynomial in s
id_RankFreeModule(ideal s,ring lmRing,ring tailRing)791 long id_RankFreeModule (ideal s, ring lmRing, ring tailRing)
792 {
793 id_TestTail(s, lmRing, tailRing);
794
795 long j = 0;
796
797 if (rRing_has_Comp(tailRing) && rRing_has_Comp(lmRing))
798 {
799 poly *p=s->m;
800 for (unsigned int l=IDELEMS(s); l > 0; --l, ++p)
801 if (*p != NULL)
802 {
803 pp_Test(*p, lmRing, tailRing);
804 const long k = p_MaxComp(*p, lmRing, tailRing);
805 if (k>j) j = k;
806 }
807 }
808
809 return j; // return -1;
810 }
811
812 /*2
813 *returns true if id is homogenous with respect to the aktual weights
814 */
id_HomIdeal(ideal id,ideal Q,const ring r)815 BOOLEAN id_HomIdeal (ideal id, ideal Q, const ring r)
816 {
817 int i;
818 BOOLEAN b;
819 i = 0;
820 b = TRUE;
821 while ((i < IDELEMS(id)) && b)
822 {
823 b = p_IsHomogeneous(id->m[i],r);
824 i++;
825 }
826 if ((b) && (Q!=NULL) && (IDELEMS(Q)>0))
827 {
828 i=0;
829 while ((i < IDELEMS(Q)) && b)
830 {
831 b = p_IsHomogeneous(Q->m[i],r);
832 i++;
833 }
834 }
835 return b;
836 }
837
838 /*2
839 *initialized a field with r numbers between beg and end for the
840 *procedure idNextChoise
841 */
idInitChoise(int r,int beg,int end,BOOLEAN * endch,int * choise)842 void idInitChoise (int r,int beg,int end,BOOLEAN *endch,int * choise)
843 {
844 /*returns the first choise of r numbers between beg and end*/
845 int i;
846 for (i=0; i<r; i++)
847 {
848 choise[i] = 0;
849 }
850 if (r <= end-beg+1)
851 for (i=0; i<r; i++)
852 {
853 choise[i] = beg+i;
854 }
855 if (r > end-beg+1)
856 *endch = TRUE;
857 else
858 *endch = FALSE;
859 }
860
861 /*2
862 *returns the next choise of r numbers between beg and end
863 */
idGetNextChoise(int r,int end,BOOLEAN * endch,int * choise)864 void idGetNextChoise (int r,int end,BOOLEAN *endch,int * choise)
865 {
866 int i = r-1,j;
867 while ((i >= 0) && (choise[i] == end))
868 {
869 i--;
870 end--;
871 }
872 if (i == -1)
873 *endch = TRUE;
874 else
875 {
876 choise[i]++;
877 for (j=i+1; j<r; j++)
878 {
879 choise[j] = choise[i]+j-i;
880 }
881 *endch = FALSE;
882 }
883 }
884
885 /*2
886 *takes the field choise of d numbers between beg and end, cancels the t-th
887 *entree and searches for the ordinal number of that d-1 dimensional field
888 * w.r.t. the algorithm of construction
889 */
idGetNumberOfChoise(int t,int d,int begin,int end,int * choise)890 int idGetNumberOfChoise(int t, int d, int begin, int end, int * choise)
891 {
892 int * localchoise,i,result=0;
893 BOOLEAN b=FALSE;
894
895 if (d<=1) return 1;
896 localchoise=(int*)omAlloc((d-1)*sizeof(int));
897 idInitChoise(d-1,begin,end,&b,localchoise);
898 while (!b)
899 {
900 result++;
901 i = 0;
902 while ((i<t) && (localchoise[i]==choise[i])) i++;
903 if (i>=t)
904 {
905 i = t+1;
906 while ((i<d) && (localchoise[i-1]==choise[i])) i++;
907 if (i>=d)
908 {
909 omFreeSize((ADDRESS)localchoise,(d-1)*sizeof(int));
910 return result;
911 }
912 }
913 idGetNextChoise(d-1,end,&b,localchoise);
914 }
915 omFreeSize((ADDRESS)localchoise,(d-1)*sizeof(int));
916 return 0;
917 }
918
919 /*2
920 *computes the binomial coefficient
921 */
binom(int n,int r)922 int binom (int n,int r)
923 {
924 int i;
925 int64 result;
926
927 if (r==0) return 1;
928 if (n-r<r) return binom(n,n-r);
929 result = n-r+1;
930 for (i=2;i<=r;i++)
931 {
932 result *= n-r+i;
933 result /= i;
934 }
935 if (result>MAX_INT_VAL)
936 {
937 WarnS("overflow in binomials");
938 result=0;
939 }
940 return (int)result;
941 }
942
943
944 /// the free module of rank i
id_FreeModule(int i,const ring r)945 ideal id_FreeModule (int i, const ring r)
946 {
947 assume(i >= 0);
948 if (r->isLPring)
949 {
950 PrintS("In order to address bimodules, the command freeAlgebra should be used.");
951 }
952 ideal h = idInit(i, i);
953
954 for (int j=0; j<i; j++)
955 {
956 h->m[j] = p_One(r);
957 p_SetComp(h->m[j],j+1,r);
958 p_SetmComp(h->m[j],r);
959 }
960
961 return h;
962 }
963
964 /*2
965 *computes recursively all monomials of a certain degree
966 *in every step the actvar-th entry in the exponential
967 *vector is incremented and the other variables are
968 *computed by recursive calls of makemonoms
969 *if the last variable is reached, the difference to the
970 *degree is computed directly
971 *vars is the number variables
972 *actvar is the actual variable to handle
973 *deg is the degree of the monomials to compute
974 *monomdeg is the actual degree of the monomial in consideration
975 */
makemonoms(int vars,int actvar,int deg,int monomdeg,const ring r)976 static void makemonoms(int vars,int actvar,int deg,int monomdeg, const ring r)
977 {
978 poly p;
979 int i=0;
980
981 if ((idpowerpoint == 0) && (actvar ==1))
982 {
983 idpower[idpowerpoint] = p_One(r);
984 monomdeg = 0;
985 }
986 while (i<=deg)
987 {
988 if (deg == monomdeg)
989 {
990 p_Setm(idpower[idpowerpoint],r);
991 idpowerpoint++;
992 return;
993 }
994 if (actvar == vars)
995 {
996 p_SetExp(idpower[idpowerpoint],actvar,deg-monomdeg,r);
997 p_Setm(idpower[idpowerpoint],r);
998 p_Test(idpower[idpowerpoint],r);
999 idpowerpoint++;
1000 return;
1001 }
1002 else
1003 {
1004 p = p_Copy(idpower[idpowerpoint],r);
1005 makemonoms(vars,actvar+1,deg,monomdeg,r);
1006 idpower[idpowerpoint] = p;
1007 }
1008 monomdeg++;
1009 p_SetExp(idpower[idpowerpoint],actvar,p_GetExp(idpower[idpowerpoint],actvar,r)+1,r);
1010 p_Setm(idpower[idpowerpoint],r);
1011 p_Test(idpower[idpowerpoint],r);
1012 i++;
1013 }
1014 }
1015
1016 #ifdef HAVE_SHIFTBBA
1017 /*2
1018 *computes recursively all letterplace monomials of a certain degree
1019 *vars is the number of original variables (lV)
1020 *deg is the degree of the monomials to compute
1021 *
1022 *NOTE: We use idpowerpoint as the last index of the previous call
1023 */
lpmakemonoms(int vars,int deg,const ring r)1024 static void lpmakemonoms(int vars, int deg, const ring r)
1025 {
1026 assume(deg <= r->N/r->isLPring);
1027 if (deg == 0)
1028 {
1029 idpower[0] = p_One(r);
1030 return;
1031 }
1032 else
1033 {
1034 lpmakemonoms(vars, deg - 1, r);
1035 }
1036
1037 int size = idpowerpoint + 1;
1038 for (int j = 2; j <= vars; j++)
1039 {
1040 for (int i = 0; i < size; i++)
1041 {
1042 idpowerpoint = (j-1)*size + i;
1043 idpower[idpowerpoint] = p_Copy(idpower[i], r);
1044 }
1045 }
1046 for (int j = 1; j <= vars; j++)
1047 {
1048 for (int i = 0; i < size; i++)
1049 {
1050 idpowerpoint = (j-1)*size + i;
1051 p_SetExp(idpower[idpowerpoint], ((deg - 1) * r->isLPring) + j, 1, r);
1052 p_Setm(idpower[idpowerpoint],r);
1053 p_Test(idpower[idpowerpoint],r);
1054 }
1055 }
1056 }
1057 #endif
1058
1059 /*2
1060 *returns the deg-th power of the maximal ideal of 0
1061 */
id_MaxIdeal(int deg,const ring r)1062 ideal id_MaxIdeal(int deg, const ring r)
1063 {
1064 if (deg < 1)
1065 {
1066 ideal I=idInit(1,1);
1067 I->m[0]=p_One(r);
1068 return I;
1069 }
1070 if (deg == 1
1071 #ifdef HAVE_SHIFTBBA
1072 && !r->isLPring
1073 #endif
1074 )
1075 {
1076 return id_MaxIdeal(r);
1077 }
1078
1079 int vars, i;
1080 #ifdef HAVE_SHIFTBBA
1081 if (r->isLPring)
1082 {
1083 vars = r->isLPring - r->LPncGenCount;
1084 i = 1;
1085 // i = vars^deg
1086 for (int j = 0; j < deg; j++)
1087 {
1088 i *= vars;
1089 }
1090 }
1091 else
1092 #endif
1093 {
1094 vars = rVar(r);
1095 i = binom(vars+deg-1,deg);
1096 }
1097 if (i<=0) return idInit(1,1);
1098 ideal id=idInit(i,1);
1099 idpower = id->m;
1100 idpowerpoint = 0;
1101 #ifdef HAVE_SHIFTBBA
1102 if (r->isLPring)
1103 {
1104 lpmakemonoms(vars, deg, r);
1105 }
1106 else
1107 #endif
1108 {
1109 makemonoms(vars,1,deg,0,r);
1110 }
1111 idpower = NULL;
1112 idpowerpoint = 0;
1113 return id;
1114 }
1115
id_NextPotence(ideal given,ideal result,int begin,int end,int deg,int restdeg,poly ap,const ring r)1116 static void id_NextPotence(ideal given, ideal result,
1117 int begin, int end, int deg, int restdeg, poly ap, const ring r)
1118 {
1119 poly p;
1120 int i;
1121
1122 p = p_Power(p_Copy(given->m[begin],r),restdeg,r);
1123 i = result->nrows;
1124 result->m[i] = p_Mult_q(p_Copy(ap,r),p,r);
1125 //PrintS(".");
1126 (result->nrows)++;
1127 if (result->nrows >= IDELEMS(result))
1128 {
1129 pEnlargeSet(&(result->m),IDELEMS(result),16);
1130 IDELEMS(result) += 16;
1131 }
1132 if (begin == end) return;
1133 for (i=restdeg-1;i>0;i--)
1134 {
1135 p = p_Power(p_Copy(given->m[begin],r),i,r);
1136 p = p_Mult_q(p_Copy(ap,r),p,r);
1137 id_NextPotence(given, result, begin+1, end, deg, restdeg-i, p,r);
1138 p_Delete(&p,r);
1139 }
1140 id_NextPotence(given, result, begin+1, end, deg, restdeg, ap,r);
1141 }
1142
id_Power(ideal given,int exp,const ring r)1143 ideal id_Power(ideal given,int exp, const ring r)
1144 {
1145 ideal result,temp;
1146 poly p1;
1147 int i;
1148
1149 if (idIs0(given)) return idInit(1,1);
1150 temp = id_Copy(given,r);
1151 idSkipZeroes(temp);
1152 i = binom(IDELEMS(temp)+exp-1,exp);
1153 result = idInit(i,1);
1154 result->nrows = 0;
1155 //Print("ideal contains %d elements\n",i);
1156 p1=p_One(r);
1157 id_NextPotence(temp,result,0,IDELEMS(temp)-1,exp,exp,p1,r);
1158 p_Delete(&p1,r);
1159 id_Delete(&temp,r);
1160 result->nrows = 1;
1161 id_DelEquals(result,r);
1162 idSkipZeroes(result);
1163 return result;
1164 }
1165
1166 /*2
1167 *skips all zeroes and double elements, searches also for units
1168 */
id_Compactify(ideal id,const ring r)1169 void id_Compactify(ideal id, const ring r)
1170 {
1171 int i;
1172 BOOLEAN b=FALSE;
1173
1174 i = IDELEMS(id)-1;
1175 while ((! b) && (i>=0))
1176 {
1177 b=p_IsUnit(id->m[i],r);
1178 i--;
1179 }
1180 if (b)
1181 {
1182 for(i=IDELEMS(id)-1;i>=0;i--) p_Delete(&id->m[i],r);
1183 id->m[0]=p_One(r);
1184 }
1185 else
1186 {
1187 id_DelMultiples(id,r);
1188 }
1189 idSkipZeroes(id);
1190 }
1191
1192 /// returns the ideals of initial terms
id_Head(ideal h,const ring r)1193 ideal id_Head(ideal h,const ring r)
1194 {
1195 ideal m = idInit(IDELEMS(h),h->rank);
1196
1197 for (int i=IDELEMS(h)-1;i>=0; i--)
1198 if (h->m[i]!=NULL)
1199 m->m[i]=p_Head(h->m[i],r);
1200
1201 return m;
1202 }
1203
id_Homogen(ideal h,int varnum,const ring r)1204 ideal id_Homogen(ideal h, int varnum,const ring r)
1205 {
1206 ideal m = idInit(IDELEMS(h),h->rank);
1207 int i;
1208
1209 for (i=IDELEMS(h)-1;i>=0; i--)
1210 {
1211 m->m[i]=p_Homogen(h->m[i],varnum,r);
1212 }
1213 return m;
1214 }
1215
1216 /*------------------type conversions----------------*/
id_Vec2Ideal(poly vec,const ring R)1217 ideal id_Vec2Ideal(poly vec, const ring R)
1218 {
1219 ideal result=idInit(1,1);
1220 omFree((ADDRESS)result->m);
1221 p_Vec2Polys(vec, &(result->m), &(IDELEMS(result)),R);
1222 return result;
1223 }
1224
1225 /// for julia: convert an array of poly to vector
id_Array2Vector(poly * m,unsigned n,const ring R)1226 poly id_Array2Vector(poly *m, unsigned n, const ring R)
1227 {
1228 poly h;
1229 int l;
1230 sBucket_pt bucket = sBucketCreate(R);
1231
1232 for(unsigned j=0;j<n ;j++)
1233 {
1234 h = m[j];
1235 if (h!=NULL)
1236 {
1237 h=p_Copy(h, R);
1238 l=pLength(h);
1239 p_SetCompP(h,j+1, R);
1240 sBucket_Merge_p(bucket, h, l);
1241 }
1242 }
1243 sBucketClearMerge(bucket, &h, &l);
1244 sBucketDestroy(&bucket);
1245 return h;
1246 }
1247
1248 /// converts mat to module, destroys mat
id_Matrix2Module(matrix mat,const ring R)1249 ideal id_Matrix2Module(matrix mat, const ring R)
1250 {
1251 int mc=MATCOLS(mat);
1252 int mr=MATROWS(mat);
1253 ideal result = idInit(mc,mr);
1254 int i,j,l;
1255 poly h;
1256 sBucket_pt bucket = sBucketCreate(R);
1257
1258 for(j=0;j<mc /*MATCOLS(mat)*/;j++) /* j is also index in result->m */
1259 {
1260 for (i=0;i<mr /*MATROWS(mat)*/;i++)
1261 {
1262 h = MATELEM0(mat,i,j);
1263 if (h!=NULL)
1264 {
1265 l=pLength(h);
1266 MATELEM0(mat,i,j)=NULL;
1267 p_SetCompP(h,i+1, R);
1268 sBucket_Merge_p(bucket, h, l);
1269 }
1270 }
1271 sBucketClearMerge(bucket, &(result->m[j]), &l);
1272 }
1273 sBucketDestroy(&bucket);
1274
1275 // obachman: need to clean this up
1276 id_Delete((ideal*) &mat,R);
1277 return result;
1278 }
1279
1280 /*2
1281 * converts a module into a matrix, destroyes the input
1282 */
id_Module2Matrix(ideal mod,const ring R)1283 matrix id_Module2Matrix(ideal mod, const ring R)
1284 {
1285 matrix result = mpNew(mod->rank,IDELEMS(mod));
1286 long i; long cp;
1287 poly p,h;
1288
1289 for(i=0;i<IDELEMS(mod);i++)
1290 {
1291 p=pReverse(mod->m[i]);
1292 mod->m[i]=NULL;
1293 while (p!=NULL)
1294 {
1295 h=p;
1296 pIter(p);
1297 pNext(h)=NULL;
1298 cp = si_max(1L,p_GetComp(h, R)); // if used for ideals too
1299 //cp = p_GetComp(h,R);
1300 p_SetComp(h,0,R);
1301 p_SetmComp(h,R);
1302 #ifdef TEST
1303 if (cp>mod->rank)
1304 {
1305 Print("## inv. rank %ld -> %ld\n",mod->rank,cp);
1306 int k,l,o=mod->rank;
1307 mod->rank=cp;
1308 matrix d=mpNew(mod->rank,IDELEMS(mod));
1309 for (l=0; l<o; l++)
1310 {
1311 for (k=0; k<IDELEMS(mod); k++)
1312 {
1313 MATELEM0(d,l,k)=MATELEM0(result,l,k);
1314 MATELEM0(result,l,k)=NULL;
1315 }
1316 }
1317 id_Delete((ideal *)&result,R);
1318 result=d;
1319 }
1320 #endif
1321 MATELEM0(result,cp-1,i) = p_Add_q(MATELEM0(result,cp-1,i),h,R);
1322 }
1323 }
1324 // obachman 10/99: added the following line, otherwise memory leack!
1325 id_Delete(&mod,R);
1326 return result;
1327 }
1328
id_Module2formatedMatrix(ideal mod,int rows,int cols,const ring R)1329 matrix id_Module2formatedMatrix(ideal mod,int rows, int cols, const ring R)
1330 {
1331 matrix result = mpNew(rows,cols);
1332 int i,cp,r=id_RankFreeModule(mod,R),c=IDELEMS(mod);
1333 poly p,h;
1334
1335 if (r>rows) r = rows;
1336 if (c>cols) c = cols;
1337 for(i=0;i<c;i++)
1338 {
1339 p=pReverse(mod->m[i]);
1340 mod->m[i]=NULL;
1341 while (p!=NULL)
1342 {
1343 h=p;
1344 pIter(p);
1345 pNext(h)=NULL;
1346 cp = p_GetComp(h,R);
1347 if (cp<=r)
1348 {
1349 p_SetComp(h,0,R);
1350 p_SetmComp(h,R);
1351 MATELEM0(result,cp-1,i) = p_Add_q(MATELEM0(result,cp-1,i),h,R);
1352 }
1353 else
1354 p_Delete(&h,R);
1355 }
1356 }
1357 id_Delete(&mod,R);
1358 return result;
1359 }
1360
id_ResizeModule(ideal mod,int rows,int cols,const ring R)1361 ideal id_ResizeModule(ideal mod,int rows, int cols, const ring R)
1362 {
1363 // columns?
1364 if (cols!=IDELEMS(mod))
1365 {
1366 for(int i=IDELEMS(mod)-1;i>=cols;i--) p_Delete(&mod->m[i],R);
1367 pEnlargeSet(&(mod->m),IDELEMS(mod),cols-IDELEMS(mod));
1368 IDELEMS(mod)=cols;
1369 }
1370 // rows?
1371 if (rows<mod->rank)
1372 {
1373 for(int i=IDELEMS(mod)-1;i>=0;i--)
1374 {
1375 if (mod->m[i]!=NULL)
1376 {
1377 while((mod->m[i]!=NULL) && (p_GetComp(mod->m[i],R)>rows))
1378 mod->m[i]=p_LmDeleteAndNext(mod->m[i],R);
1379 poly p=mod->m[i];
1380 while(pNext(p)!=NULL)
1381 {
1382 if (p_GetComp(pNext(p),R)>rows)
1383 pNext(p)=p_LmDeleteAndNext(pNext(p),R);
1384 else
1385 pIter(p);
1386 }
1387 }
1388 }
1389 }
1390 mod->rank=rows;
1391 return mod;
1392 }
1393
1394 /*2
1395 * substitute the n-th variable by the monomial e in id
1396 * destroy id
1397 */
id_Subst(ideal id,int n,poly e,const ring r)1398 ideal id_Subst(ideal id, int n, poly e, const ring r)
1399 {
1400 int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
1401 ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
1402
1403 res->rank = id->rank;
1404 for(k--;k>=0;k--)
1405 {
1406 res->m[k]=p_Subst(id->m[k],n,e,r);
1407 id->m[k]=NULL;
1408 }
1409 id_Delete(&id,r);
1410 return res;
1411 }
1412
id_HomModule(ideal m,ideal Q,intvec ** w,const ring R)1413 BOOLEAN id_HomModule(ideal m, ideal Q, intvec **w, const ring R)
1414 {
1415 if (w!=NULL) *w=NULL;
1416 if ((Q!=NULL) && (!id_HomIdeal(Q,NULL,R))) return FALSE;
1417 if (idIs0(m))
1418 {
1419 if (w!=NULL) (*w)=new intvec(m->rank);
1420 return TRUE;
1421 }
1422
1423 long cmax=1,order=0,ord,* diff,diffmin=32000;
1424 int *iscom;
1425 int i;
1426 poly p=NULL;
1427 pFDegProc d;
1428 if (R->pLexOrder && (R->order[0]==ringorder_lp))
1429 d=p_Totaldegree;
1430 else
1431 d=R->pFDeg;
1432 int length=IDELEMS(m);
1433 poly* P=m->m;
1434 poly* F=(poly*)omAlloc(length*sizeof(poly));
1435 for (i=length-1;i>=0;i--)
1436 {
1437 p=F[i]=P[i];
1438 cmax=si_max(cmax,p_MaxComp(p,R));
1439 }
1440 cmax++;
1441 diff = (long *)omAlloc0(cmax*sizeof(long));
1442 if (w!=NULL) *w=new intvec(cmax-1);
1443 iscom = (int *)omAlloc0(cmax*sizeof(int));
1444 i=0;
1445 while (i<=length)
1446 {
1447 if (i<length)
1448 {
1449 p=F[i];
1450 while ((p!=NULL) && (iscom[__p_GetComp(p,R)]==0)) pIter(p);
1451 }
1452 if ((p==NULL) && (i<length))
1453 {
1454 i++;
1455 }
1456 else
1457 {
1458 if (p==NULL) /* && (i==length) */
1459 {
1460 i=0;
1461 while ((i<length) && (F[i]==NULL)) i++;
1462 if (i>=length) break;
1463 p = F[i];
1464 }
1465 //if (pLexOrder && (currRing->order[0]==ringorder_lp))
1466 // order=pTotaldegree(p);
1467 //else
1468 // order = p->order;
1469 // order = pFDeg(p,currRing);
1470 order = d(p,R) +diff[__p_GetComp(p,R)];
1471 //order += diff[pGetComp(p)];
1472 p = F[i];
1473 //Print("Actual p=F[%d]: ",i);pWrite(p);
1474 F[i] = NULL;
1475 i=0;
1476 }
1477 while (p!=NULL)
1478 {
1479 if (R->pLexOrder && (R->order[0]==ringorder_lp))
1480 ord=p_Totaldegree(p,R);
1481 else
1482 // ord = p->order;
1483 ord = R->pFDeg(p,R);
1484 if (iscom[__p_GetComp(p,R)]==0)
1485 {
1486 diff[__p_GetComp(p,R)] = order-ord;
1487 iscom[__p_GetComp(p,R)] = 1;
1488 /*
1489 *PrintS("new diff: ");
1490 *for (j=0;j<cmax;j++) Print("%d ",diff[j]);
1491 *PrintLn();
1492 *PrintS("new iscom: ");
1493 *for (j=0;j<cmax;j++) Print("%d ",iscom[j]);
1494 *PrintLn();
1495 *Print("new set %d, order %d, ord %d, diff %d\n",pGetComp(p),order,ord,diff[pGetComp(p)]);
1496 */
1497 }
1498 else
1499 {
1500 /*
1501 *PrintS("new diff: ");
1502 *for (j=0;j<cmax;j++) Print("%d ",diff[j]);
1503 *PrintLn();
1504 *Print("order %d, ord %d, diff %d\n",order,ord,diff[pGetComp(p)]);
1505 */
1506 if (order != (ord+diff[__p_GetComp(p,R)]))
1507 {
1508 omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
1509 omFreeSize((ADDRESS) diff,cmax*sizeof(long));
1510 omFreeSize((ADDRESS) F,length*sizeof(poly));
1511 delete *w;*w=NULL;
1512 return FALSE;
1513 }
1514 }
1515 pIter(p);
1516 }
1517 }
1518 omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
1519 omFreeSize((ADDRESS) F,length*sizeof(poly));
1520 for (i=1;i<cmax;i++) (**w)[i-1]=(int)(diff[i]);
1521 for (i=1;i<cmax;i++)
1522 {
1523 if (diff[i]<diffmin) diffmin=diff[i];
1524 }
1525 if (w!=NULL)
1526 {
1527 for (i=1;i<cmax;i++)
1528 {
1529 (**w)[i-1]=(int)(diff[i]-diffmin);
1530 }
1531 }
1532 omFreeSize((ADDRESS) diff,cmax*sizeof(long));
1533 return TRUE;
1534 }
1535
id_Jet(const ideal i,int d,const ring R)1536 ideal id_Jet(const ideal i,int d, const ring R)
1537 {
1538 ideal r=idInit((i->nrows)*(i->ncols),i->rank);
1539 r->nrows = i-> nrows;
1540 r->ncols = i-> ncols;
1541 //r->rank = i-> rank;
1542
1543 for(long k=((long)(i->nrows))*((long)(i->ncols))-1;k>=0; k--)
1544 r->m[k]=pp_Jet(i->m[k],d,R);
1545
1546 return r;
1547 }
1548
id_JetW(const ideal i,int d,intvec * iv,const ring R)1549 ideal id_JetW(const ideal i,int d, intvec * iv, const ring R)
1550 {
1551 ideal r=idInit(IDELEMS(i),i->rank);
1552 if (ecartWeights!=NULL)
1553 {
1554 WerrorS("cannot compute weighted jets now");
1555 }
1556 else
1557 {
1558 int *w=iv2array(iv,R);
1559 int k;
1560 for(k=0; k<IDELEMS(i); k++)
1561 {
1562 r->m[k]=pp_JetW(i->m[k],d,w,R);
1563 }
1564 omFreeSize((ADDRESS)w,(rVar(R)+1)*sizeof(int));
1565 }
1566 return r;
1567 }
1568
1569 /*3
1570 * searches for the next unit in the components of the module arg and
1571 * returns the first one;
1572 */
id_ReadOutPivot(ideal arg,int * comp,const ring r)1573 int id_ReadOutPivot(ideal arg,int* comp, const ring r)
1574 {
1575 if (idIs0(arg)) return -1;
1576 int i=0,j, generator=-1;
1577 int rk_arg=arg->rank; //idRankFreeModule(arg);
1578 int * componentIsUsed =(int *)omAlloc((rk_arg+1)*sizeof(int));
1579 poly p;
1580
1581 while ((generator<0) && (i<IDELEMS(arg)))
1582 {
1583 memset(componentIsUsed,0,(rk_arg+1)*sizeof(int));
1584 p = arg->m[i];
1585 while (p!=NULL)
1586 {
1587 j = __p_GetComp(p,r);
1588 if (componentIsUsed[j]==0)
1589 {
1590 if (p_LmIsConstantComp(p,r) &&
1591 (!rField_is_Ring(r) || n_IsUnit(pGetCoeff(p),r->cf)))
1592 {
1593 generator = i;
1594 componentIsUsed[j] = 1;
1595 }
1596 else
1597 {
1598 componentIsUsed[j] = -1;
1599 }
1600 }
1601 else if (componentIsUsed[j]>0)
1602 {
1603 (componentIsUsed[j])++;
1604 }
1605 pIter(p);
1606 }
1607 i++;
1608 }
1609 i = 0;
1610 *comp = -1;
1611 for (j=0;j<=rk_arg;j++)
1612 {
1613 if (componentIsUsed[j]>0)
1614 {
1615 if ((*comp==-1) || (componentIsUsed[j]<i))
1616 {
1617 *comp = j;
1618 i= componentIsUsed[j];
1619 }
1620 }
1621 }
1622 omFree(componentIsUsed);
1623 return generator;
1624 }
1625
1626 #if 0
1627 static void idDeleteComp(ideal arg,int red_comp)
1628 {
1629 int i,j;
1630 poly p;
1631
1632 for (i=IDELEMS(arg)-1;i>=0;i--)
1633 {
1634 p = arg->m[i];
1635 while (p!=NULL)
1636 {
1637 j = pGetComp(p);
1638 if (j>red_comp)
1639 {
1640 pSetComp(p,j-1);
1641 pSetm(p);
1642 }
1643 pIter(p);
1644 }
1645 }
1646 (arg->rank)--;
1647 }
1648 #endif
1649
id_QHomWeight(ideal id,const ring r)1650 intvec * id_QHomWeight(ideal id, const ring r)
1651 {
1652 poly head, tail;
1653 int k;
1654 int in=IDELEMS(id)-1, ready=0, all=0,
1655 coldim=rVar(r), rowmax=2*coldim;
1656 if (in<0) return NULL;
1657 intvec *imat=new intvec(rowmax+1,coldim,0);
1658
1659 do
1660 {
1661 head = id->m[in--];
1662 if (head!=NULL)
1663 {
1664 tail = pNext(head);
1665 while (tail!=NULL)
1666 {
1667 all++;
1668 for (k=1;k<=coldim;k++)
1669 IMATELEM(*imat,all,k) = p_GetExpDiff(head,tail,k,r);
1670 if (all==rowmax)
1671 {
1672 ivTriangIntern(imat, ready, all);
1673 if (ready==coldim)
1674 {
1675 delete imat;
1676 return NULL;
1677 }
1678 }
1679 pIter(tail);
1680 }
1681 }
1682 } while (in>=0);
1683 if (all>ready)
1684 {
1685 ivTriangIntern(imat, ready, all);
1686 if (ready==coldim)
1687 {
1688 delete imat;
1689 return NULL;
1690 }
1691 }
1692 intvec *result = ivSolveKern(imat, ready);
1693 delete imat;
1694 return result;
1695 }
1696
id_IsZeroDim(ideal I,const ring r)1697 BOOLEAN id_IsZeroDim(ideal I, const ring r)
1698 {
1699 BOOLEAN *UsedAxis=(BOOLEAN *)omAlloc0(rVar(r)*sizeof(BOOLEAN));
1700 int i,n;
1701 poly po;
1702 BOOLEAN res=TRUE;
1703 for(i=IDELEMS(I)-1;i>=0;i--)
1704 {
1705 po=I->m[i];
1706 if ((po!=NULL) &&((n=p_IsPurePower(po,r))!=0)) UsedAxis[n-1]=TRUE;
1707 }
1708 for(i=rVar(r)-1;i>=0;i--)
1709 {
1710 if(UsedAxis[i]==FALSE) {res=FALSE; break;} // not zero-dim.
1711 }
1712 omFreeSize(UsedAxis,rVar(r)*sizeof(BOOLEAN));
1713 return res;
1714 }
1715
id_Normalize(ideal I,const ring r)1716 void id_Normalize(ideal I,const ring r) /* for ideal/matrix */
1717 {
1718 if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
1719 int i;
1720 for(i=I->nrows*I->ncols-1;i>=0;i--)
1721 {
1722 p_Normalize(I->m[i],r);
1723 }
1724 }
1725
id_MinDegW(ideal M,intvec * w,const ring r)1726 int id_MinDegW(ideal M,intvec *w, const ring r)
1727 {
1728 int d=-1;
1729 for(int i=0;i<IDELEMS(M);i++)
1730 {
1731 if (M->m[i]!=NULL)
1732 {
1733 int d0=p_MinDeg(M->m[i],w,r);
1734 if(-1<d0&&((d0<d)||(d==-1)))
1735 d=d0;
1736 }
1737 }
1738 return d;
1739 }
1740
1741 // #include "kernel/clapsing.h"
1742
1743 /*2
1744 * transpose a module
1745 */
id_Transp(ideal a,const ring rRing)1746 ideal id_Transp(ideal a, const ring rRing)
1747 {
1748 int r = a->rank, c = IDELEMS(a);
1749 ideal b = idInit(r,c);
1750
1751 int i;
1752 for (i=c; i>0; i--)
1753 {
1754 poly p=a->m[i-1];
1755 while(p!=NULL)
1756 {
1757 poly h=p_Head(p, rRing);
1758 int co=__p_GetComp(h, rRing)-1;
1759 p_SetComp(h, i, rRing);
1760 p_Setm(h, rRing);
1761 h->next=b->m[co];
1762 b->m[co]=h;
1763 pIter(p);
1764 }
1765 }
1766 for (i=IDELEMS(b)-1; i>=0; i--)
1767 {
1768 poly p=b->m[i];
1769 if(p!=NULL)
1770 {
1771 b->m[i]=p_SortMerge(p,rRing,TRUE);
1772 }
1773 }
1774 return b;
1775 }
1776
1777 /*2
1778 * The following is needed to compute the image of certain map used in
1779 * the computation of cohomologies via BGG
1780 * let M = { w_1, ..., w_k }, k = size(M) == ncols(M), n = nvars(currRing).
1781 * assuming that nrows(M) <= m*n; the procedure computes:
1782 * transpose(M) * transpose( var(1) I_m | ... | var(n) I_m ) :== transpose(module{f_1, ... f_k}),
1783 * where f_i = \sum_{j=1}^{m} (w_i, v_j) gen(j), (w_i, v_j) is a `scalar` multiplication.
1784 * that is, if w_i = (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m) then
1785
1786 (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m)
1787 * var_1 ... var_1 | var_2 ... var_2 | ... | var_n ... var(n)
1788 * gen_1 ... gen_m | gen_1 ... gen_m | ... | gen_1 ... gen_m
1789 + =>
1790 f_i =
1791
1792 a^1_1 * var(1) * gen(1) + ... + a^1_m * var(1) * gen(m) +
1793 a^2_1 * var(2) * gen(1) + ... + a^2_m * var(2) * gen(m) +
1794 ...
1795 a^n_1 * var(n) * gen(1) + ... + a^n_m * var(n) * gen(m);
1796
1797 NOTE: for every f_i we run only ONCE along w_i saving partial sums into a temporary array of polys of size m
1798 */
id_TensorModuleMult(const int m,const ideal M,const ring rRing)1799 ideal id_TensorModuleMult(const int m, const ideal M, const ring rRing)
1800 {
1801 // #ifdef DEBU
1802 // WarnS("tensorModuleMult!!!!");
1803
1804 assume(m > 0);
1805 assume(M != NULL);
1806
1807 const int n = rRing->N;
1808
1809 assume(M->rank <= m * n);
1810
1811 const int k = IDELEMS(M);
1812
1813 ideal idTemp = idInit(k,m); // = {f_1, ..., f_k }
1814
1815 for( int i = 0; i < k; i++ ) // for every w \in M
1816 {
1817 poly pTempSum = NULL;
1818
1819 poly w = M->m[i];
1820
1821 while(w != NULL) // for each term of w...
1822 {
1823 poly h = p_Head(w, rRing);
1824
1825 const int gen = __p_GetComp(h, rRing); // 1 ...
1826
1827 assume(gen > 0);
1828 assume(gen <= n*m);
1829
1830 // TODO: write a formula with %, / instead of while!
1831 /*
1832 int c = gen;
1833 int v = 1;
1834 while(c > m)
1835 {
1836 c -= m;
1837 v++;
1838 }
1839 */
1840
1841 int cc = gen % m;
1842 if( cc == 0) cc = m;
1843 int vv = 1 + (gen - cc) / m;
1844
1845 // assume( cc == c );
1846 // assume( vv == v );
1847
1848 // 1<= c <= m
1849 assume( cc > 0 );
1850 assume( cc <= m );
1851
1852 assume( vv > 0 );
1853 assume( vv <= n );
1854
1855 assume( (cc + (vv-1)*m) == gen );
1856
1857 p_IncrExp(h, vv, rRing); // h *= var(j) && // p_AddExp(h, vv, 1, rRing);
1858 p_SetComp(h, cc, rRing);
1859
1860 p_Setm(h, rRing); // addjust degree after the previous steps!
1861
1862 pTempSum = p_Add_q(pTempSum, h, rRing); // it is slow since h will be usually put to the back of pTempSum!!!
1863
1864 pIter(w);
1865 }
1866
1867 idTemp->m[i] = pTempSum;
1868 }
1869
1870 // simplify idTemp???
1871
1872 ideal idResult = id_Transp(idTemp, rRing);
1873
1874 id_Delete(&idTemp, rRing);
1875
1876 return(idResult);
1877 }
1878
id_ChineseRemainder(ideal * xx,number * q,int rl,const ring r)1879 ideal id_ChineseRemainder(ideal *xx, number *q, int rl, const ring r)
1880 {
1881 int cnt=0;int rw=0; int cl=0;
1882 int i,j;
1883 // find max. size of xx[.]:
1884 for(j=rl-1;j>=0;j--)
1885 {
1886 i=IDELEMS(xx[j])*xx[j]->nrows;
1887 if (i>cnt) cnt=i;
1888 if (xx[j]->nrows >rw) rw=xx[j]->nrows; // for lifting matrices
1889 if (xx[j]->ncols >cl) cl=xx[j]->ncols; // for lifting matrices
1890 }
1891 if (rw*cl !=cnt)
1892 {
1893 WerrorS("format mismatch in CRT");
1894 return NULL;
1895 }
1896 ideal result=idInit(cnt,xx[0]->rank);
1897 result->nrows=rw; // for lifting matrices
1898 result->ncols=cl; // for lifting matrices
1899 number *x=(number *)omAlloc(rl*sizeof(number));
1900 poly *p=(poly *)omAlloc(rl*sizeof(poly));
1901 CFArray inv_cache(rl);
1902 EXTERN_VAR int n_SwitchChinRem; //TEST
1903 int save_n_SwitchChinRem=n_SwitchChinRem;
1904 n_SwitchChinRem=1;
1905 for(i=cnt-1;i>=0;i--)
1906 {
1907 for(j=rl-1;j>=0;j--)
1908 {
1909 if(i>=IDELEMS(xx[j])*xx[j]->nrows) // out of range of this ideal
1910 p[j]=NULL;
1911 else
1912 p[j]=xx[j]->m[i];
1913 }
1914 result->m[i]=p_ChineseRemainder(p,x,q,rl,inv_cache,r);
1915 for(j=rl-1;j>=0;j--)
1916 {
1917 if(i<IDELEMS(xx[j])*xx[j]->nrows) xx[j]->m[i]=p[j];
1918 }
1919 }
1920 n_SwitchChinRem=save_n_SwitchChinRem;
1921 omFreeSize(p,rl*sizeof(poly));
1922 omFreeSize(x,rl*sizeof(number));
1923 for(i=rl-1;i>=0;i--) id_Delete(&(xx[i]),r);
1924 omFreeSize(xx,rl*sizeof(ideal));
1925 return result;
1926 }
1927
id_Shift(ideal M,int s,const ring r)1928 void id_Shift(ideal M, int s, const ring r)
1929 {
1930 // id_Test( M, r );
1931
1932 // assume( s >= 0 ); // negative is also possible // TODO: verify input ideal in such a case!?
1933
1934 for(int i=IDELEMS(M)-1; i>=0;i--)
1935 p_Shift(&(M->m[i]),s,r);
1936
1937 M->rank += s;
1938
1939 // id_Test( M, r );
1940 }
1941
id_Delete_Pos(const ideal I,const int p,const ring r)1942 ideal id_Delete_Pos(const ideal I, const int p, const ring r)
1943 {
1944 if ((p<0)||(p>=IDELEMS(I))) return NULL;
1945 ideal ret=idInit(IDELEMS(I)-1,I->rank);
1946 for(int i=0;i<p;i++) ret->m[i]=p_Copy(I->m[i],r);
1947 for(int i=p+1;i<IDELEMS(I);i++) ret->m[i-1]=p_Copy(I->m[i],r);
1948 return ret;
1949 }
1950