source: git/libpolys/polys/simpleideals.cc @ 6f3273

spielwiese
Last change on this file since 6f3273 was 6f3273, checked in by Hans Schoenemann <hannes@…>, 13 years ago
fix BOOLEAN id_IsConstant(ideal id, const ri ng r)
  • Property mode set to 100644
File size: 53.8 KB
RevLine 
[35aab3]1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
[341696]4/* $Id$ */
[35aab3]5/*
6* ABSTRACT - all basic methods to manipulate ideals
7*/
8
[2ad10e9]9
[35aab3]10/* includes */
[2ad10e9]11#include "config.h"
12#include <misc/auxiliary.h>
[f71e8c5]13#include <misc/options.h>
14#include <omalloc/omalloc.h>
15#include <polys/monomials/p_polys.h>
[f5c2d02]16#include <misc/intvec.h>
[2ad10e9]17
18#include "simpleideals.h"
[35aab3]19
20/*2
21* initialise an ideal
22*/
23ideal idInit(int idsize, int rank)
24{
25  /*- initialise an ideal -*/
26  ideal hh = (ideal )omAllocBin(sip_sideal_bin);
27  hh->nrows = 1;
28  hh->rank = rank;
29  IDELEMS(hh) = idsize;
30  if (idsize>0)
31  {
32    hh->m = (poly *)omAlloc0(idsize*sizeof(poly));
33  }
34  else
35    hh->m=NULL;
36  return hh;
37}
38
[e9c3b2]39#ifdef PDEBUG
[e070895]40// this is only for outputting an ideal within the debugger
[645a19]41void idShow(const ideal id, const ring lmRing, const ring tailRing, const int debugPrint)
[35aab3]42{
[645a19]43  assume( debugPrint >= 0 );
[bead81]44
[52e2f6]45  if( id == NULL )
[f44fb9]46    PrintS("(NULL)");
[52e2f6]47  else
[35aab3]48  {
[6867f5]49    Print("Module of rank %ld,real rank %ld and %d generators.\n",
[f71e8c5]50          id->rank,id_RankFreeModule(id, lmRing, tailRing),IDELEMS(id));
[645a19]51
52    int j = (id->ncols*id->nrows) - 1;
53    while ((j > 0) && (id->m[j]==NULL)) j--;
54    for (int i = 0; i <= j; i++)
[35aab3]55    {
[645a19]56      Print("generator %d: ",i); p_DebugPrint(id->m[i], lmRing, tailRing, debugPrint);
[35aab3]57    }
58  }
59}
[e070895]60#endif
[35aab3]61
[dd5534]62/* index of generator with leading term in ground ring (if any);
63   otherwise -1 */
[f71e8c5]64int id_PosConstant(ideal id, const ring r)
[dd5534]65{
66  int k;
67  for (k = IDELEMS(id)-1; k>=0; k--)
68  {
[f71e8c5]69    if (p_LmIsConstantComp(id->m[k], r) == TRUE)
[dd5534]70      return k;
71  }
72  return -1;
73}
74
[35aab3]75/*2
76* initialise the maximal ideal (at 0)
77*/
[f71e8c5]78ideal id_MaxIdeal (const ring r)
[35aab3]79{
80  int l;
81  ideal hh=NULL;
82
[f71e8c5]83  hh=idInit(rVar(r),1);
84  for (l=0; l<rVar(r); l++)
[35aab3]85  {
[f71e8c5]86    hh->m[l] = p_One(r);
87    p_SetExp(hh->m[l],l+1,1,r);
88    p_Setm(hh->m[l],r);
[35aab3]89  }
90  return hh;
91}
92
93/*2
94* deletes an ideal/matrix
95*/
96void id_Delete (ideal * h, ring r)
97{
98  int j,elems;
99  if (*h == NULL)
100    return;
101  elems=j=(*h)->nrows*(*h)->ncols;
102  if (j>0)
103  {
104    do
105    {
106      p_Delete(&((*h)->m[--j]), r);
107    }
108    while (j>0);
109    omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*elems);
110  }
111  omFreeBin((ADDRESS)*h, sip_sideal_bin);
112  *h=NULL;
113}
114
115
116/*2
117* Shallowdeletes an ideal/matrix
118*/
119void id_ShallowDelete (ideal *h, ring r)
120{
121  int j,elems;
122  if (*h == NULL)
123    return;
124  elems=j=(*h)->nrows*(*h)->ncols;
125  if (j>0)
126  {
127    do
128    {
129      p_ShallowDelete(&((*h)->m[--j]), r);
130    }
131    while (j>0);
132    omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*elems);
133  }
134  omFreeBin((ADDRESS)*h, sip_sideal_bin);
135  *h=NULL;
136}
137
138/*2
139*gives an ideal the minimal possible size
140*/
141void idSkipZeroes (ideal ide)
142{
143  int k;
144  int j = -1;
145  BOOLEAN change=FALSE;
146  for (k=0; k<IDELEMS(ide); k++)
147  {
148    if (ide->m[k] != NULL)
149    {
150      j++;
151      if (change)
152      {
153        ide->m[j] = ide->m[k];
154      }
155    }
156    else
157    {
158      change=TRUE;
159    }
160  }
161  if (change)
162  {
163    if (j == -1)
164      j = 0;
165    else
166    {
167      for (k=j+1; k<IDELEMS(ide); k++)
168        ide->m[k] = NULL;
169    }
170    pEnlargeSet(&(ide->m),IDELEMS(ide),j+1-IDELEMS(ide));
171    IDELEMS(ide) = j+1;
172  }
173}
174
[2b3caae]175/*2
176* copies the first k (>= 1) entries of the given ideal
177* and returns these as a new ideal
178* (Note that the copied polynomials may be zero.)
179*/
[f71e8c5]180ideal id_CopyFirstK (const ideal ide, const int k,const ring r)
[2b3caae]181{
182  ideal newI = idInit(k, 0);
183  for (int i = 0; i < k; i++)
[f71e8c5]184    newI->m[i] = p_Copy(ide->m[i],r);
[2b3caae]185  return newI;
186}
187
[35aab3]188/*2
189* ideal id = (id[i])
190* result is leadcoeff(id[i]) = 1
191*/
[9aa29b]192void id_Norm(ideal id, const ring r)
[35aab3]193{
[699567]194  for (int i=IDELEMS(id)-1; i>=0; i--)
[35aab3]195  {
196    if (id->m[i] != NULL)
197    {
[9aa29b]198      p_Norm(id->m[i],r);
[35aab3]199    }
200  }
201}
202
203/*2
[dd5534]204* ideal id = (id[i]), c any unit
[35aab3]205* if id[i] = c*id[j] then id[j] is deleted for j > i
206*/
[f5c2d02]207void id_DelMultiples(ideal id, const ring r)
[35aab3]208{
[699567]209  int i, j;
210  int k = IDELEMS(id)-1;
211  for (i=k; i>=0; i--)
[35aab3]212  {
213    if (id->m[i]!=NULL)
214    {
[699567]215      for (j=k; j>i; j--)
[35aab3]216      {
[dd5534]217        if (id->m[j]!=NULL)
[35aab3]218        {
[dd5534]219#ifdef HAVE_RINGS
[f5c2d02]220          if (rField_is_Ring(r))
[dd5534]221          {
222            /* if id[j] = c*id[i] then delete id[j].
223               In the below cases of a ground field, we
224               check whether id[i] = c*id[j] and, if so,
225               delete id[j] for historical reasons (so
226               that previous output does not change) */
[f5c2d02]227            if (p_ComparePolys(id->m[j], id->m[i],r)) p_Delete(&id->m[j],r);
[dd5534]228          }
229          else
230          {
[f5c2d02]231            if (p_ComparePolys(id->m[i], id->m[j],r)) p_Delete(&id->m[j],r);
[dd5534]232          }
233#else
[f5c2d02]234          if (p_ComparePolys(id->m[i], id->m[j],r)) p_Delete(&id->m[j],r);
[3d0808]235#endif   
[35aab3]236        }
237      }
238    }
239  }
240}
241
242/*2
243* ideal id = (id[i])
244* if id[i] = id[j] then id[j] is deleted for j > i
245*/
[4a08e7]246void id_DelEquals(ideal id, const ring r)
[35aab3]247{
[7ac29f]248  int i, j;
249  int k = IDELEMS(id)-1;
250  for (i=k; i>=0; i--)
[35aab3]251  {
[7ac29f]252    if (id->m[i]!=NULL)
[35aab3]253    {
[7ac29f]254      for (j=k; j>i; j--)
[35aab3]255      {
[7ac29f]256        if ((id->m[j]!=NULL)
[4a08e7]257        && (p_EqualPolys(id->m[i], id->m[j],r)))
[7ac29f]258        {
[4a08e7]259          p_Delete(&id->m[j],r);
[7ac29f]260        }
[35aab3]261      }
262    }
263  }
264}
265
266//
[a8b44d]267// Delete id[j], if Lm(j) == Lm(i) and both LC(j), LC(i) are units and j > i
[35aab3]268//
[119853]269void id_DelLmEquals(ideal id, const ring r)
[35aab3]270{
[7ac29f]271  int i, j;
272  int k = IDELEMS(id)-1;
273  for (i=k; i>=0; i--)
[35aab3]274  {
[73df93]275    if (id->m[i] != NULL)
[35aab3]276    {
[7ac29f]277      for (j=k; j>i; j--)
[35aab3]278      {
[7ac29f]279        if ((id->m[j] != NULL)
[119853]280        && p_LmEqual(id->m[i], id->m[j],r)
[a8b44d]281#ifdef HAVE_RINGS
[c9c118]282        && n_IsUnit(pGetCoeff(id->m[i]),r->cf) && n_IsUnit(pGetCoeff(id->m[j]),r->cf)
[a8b44d]283#endif
284        )
[35aab3]285        {
[119853]286          p_Delete(&id->m[j],r);
[35aab3]287        }
288      }
289    }
290  }
291}
292
[a8b44d]293//
294// delete id[j], if LT(j) == coeff*mon*LT(i) and vice versa, i.e.,
295// delete id[i], if LT(i) == coeff*mon*LT(j)
296//
[3d0808]297void id_DelDiv(ideal id, const ring r)
[35aab3]298{
[7ac29f]299  int i, j;
300  int k = IDELEMS(id)-1;
301  for (i=k; i>=0; i--)
[35aab3]302  {
[73df93]303    if (id->m[i] != NULL)
[35aab3]304    {
[7ac29f]305      for (j=k; j>i; j--)
[35aab3]306      {
[7ac29f]307        if (id->m[j]!=NULL)
[35aab3]308        {
[a8b44d]309#ifdef HAVE_RINGS
[3d0808]310          if (rField_is_Ring(r))
[a8b44d]311          {
[3d0808]312            if (p_DivisibleByRingCase(id->m[i], id->m[j],r))
[a8b44d]313            {
[3d0808]314              p_Delete(&id->m[j],r);
315            }
316            else if (p_DivisibleByRingCase(id->m[j], id->m[i],r))
317            {
318              p_Delete(&id->m[i],r);
319              break;
[a8b44d]320            }
321          }
322          else
323          {
324#endif
325          /* the case of a ground field: */
[3d0808]326          if (p_DivisibleBy(id->m[i], id->m[j],r))
[7ac29f]327          {
[3d0808]328            p_Delete(&id->m[j],r);
[7ac29f]329          }
[3d0808]330          else if (p_DivisibleBy(id->m[j], id->m[i],r))
[7ac29f]331          {
[3d0808]332            p_Delete(&id->m[i],r);
[7ac29f]333            break;
334          }
[a8b44d]335#ifdef HAVE_RINGS
336          }
[3d0808]337#endif   
[35aab3]338        }
339      }
340    }
341  }
342}
343
344/*2
345*test if the ideal has only constant polynomials
346*/
[6f3273]347BOOLEAN id_IsConstant(ideal id, const ri ng r)
[35aab3]348{
349  int k;
350  for (k = IDELEMS(id)-1; k>=0; k--)
351  {
[6f3273]352    if (!p_IsConstantPoly(id->m[k],r))
[35aab3]353      return FALSE;
354  }
355  return TRUE;
356}
357
358/*2
359* copy an ideal
360*/
[d523f3]361
362ideal id_Copy (ideal h1, const ring r)
363{
364  int i;
365  ideal h2;
366
367//#ifdef TEST
368  if (h1 == NULL)
369  {
370    h2=idInit(1,1);
371  }
372  else
373//#endif
374  {
375    h2=idInit(IDELEMS(h1),h1->rank);
376    for (i=IDELEMS(h1)-1; i>=0; i--)
377      h2->m[i] = p_Copy(h1->m[i],r);
378  }
379  return h2;
380}
[35aab3]381
382#ifdef PDEBUG
[85e68dd]383void idDBTest(ideal h1, int level, const char *f,const int l)
[35aab3]384{
385  int i;
386
387  if (h1 != NULL)
388  {
389    // assume(IDELEMS(h1) > 0); for ideal/module, does not apply to matrix
390    omCheckAddrSize(h1,sizeof(*h1));
391    omdebugAddrSize(h1->m,h1->ncols*h1->nrows*sizeof(poly));
392    /* to be able to test matrices: */
393    for (i=(h1->ncols*h1->nrows)-1; i>=0; i--)
394      _p_Test(h1->m[i], currRing, level);
395    int new_rk=idRankFreeModule(h1);
396    if(new_rk > h1->rank)
397    {
398      dReportError("wrong rank %d (should be %d) in %s:%d\n",
399                   h1->rank, new_rk, f,l);
400      omPrintAddrInfo(stderr, h1, " for ideal");
401      h1->rank=new_rk;
402    }
403  }
404}
405#endif
406
407/*3
408* for idSort: compare a and b revlex inclusive module comp.
409*/
[f71e8c5]410static int p_Comp_RevLex(poly a, poly b,BOOLEAN nolex, const ring r)
[35aab3]411{
412  if (b==NULL) return 1;
413  if (a==NULL) return -1;
414
[3d0808]415  if (nolex)
[2c872b]416  {
417    int r=pLmCmp(a,b);
418    if (r!=0) return r;
419    number h=nSub(pGetCoeff(a),pGetCoeff(b));
420    r = -1+nIsZero(h)+2*nGreaterZero(h); /* -1: <, 0:==, 1: > */
421    nDelete(&h);
422    return r;
423  }
[f71e8c5]424  int l=rVar(r);
[35aab3]425  while ((l>0) && (pGetExp(a,l)==pGetExp(b,l))) l--;
426  if (l==0)
427  {
[2c872b]428    if (pGetComp(a)==pGetComp(b))
429    {
430      number h=nSub(pGetCoeff(a),pGetCoeff(b));
[3b81b81]431      int r = -1+nIsZero(h)+2*nGreaterZero(h); /* -1: <, 0:==, 1: > */
[2c872b]432      nDelete(&h);
433      return r;
434    }
[35aab3]435    if (pGetComp(a)>pGetComp(b)) return 1;
436  }
437  else if (pGetExp(a,l)>pGetExp(b,l))
438    return 1;
439  return -1;
440}
441
442/*2
443*sorts the ideal w.r.t. the actual ringordering
444*uses lex-ordering when nolex = FALSE
445*/
446intvec *idSort(ideal id,BOOLEAN nolex)
447{
448  poly p,q;
449  intvec * result = new intvec(IDELEMS(id));
450  int i, j, actpos=0, newpos, l;
451  int diff, olddiff, lastcomp, newcomp;
452  BOOLEAN notFound;
453
454  for (i=0;i<IDELEMS(id);i++)
455  {
456    if (id->m[i]!=NULL)
457    {
458      notFound = TRUE;
459      newpos = actpos / 2;
460      diff = (actpos+1) / 2;
461      diff = (diff+1) / 2;
462      lastcomp = pComp_RevLex(id->m[i],id->m[(*result)[newpos]],nolex);
463      if (lastcomp<0)
464      {
465        newpos -= diff;
466      }
467      else if (lastcomp>0)
468      {
469        newpos += diff;
470      }
471      else
472      {
473        notFound = FALSE;
474      }
475      //while ((newpos>=0) && (newpos<actpos) && (notFound))
476      while (notFound && (newpos>=0) && (newpos<actpos))
477      {
478        newcomp = pComp_RevLex(id->m[i],id->m[(*result)[newpos]],nolex);
479        olddiff = diff;
480        if (diff>1)
481        {
482          diff = (diff+1) / 2;
483          if ((newcomp==1)
484          && (actpos-newpos>1)
485          && (diff>1)
486          && (newpos+diff>=actpos))
487          {
488            diff = actpos-newpos-1;
489          }
490          else if ((newcomp==-1)
491          && (diff>1)
492          && (newpos<diff))
493          {
494            diff = newpos;
495          }
496        }
497        if (newcomp<0)
498        {
499          if ((olddiff==1) && (lastcomp>0))
500            notFound = FALSE;
501          else
502            newpos -= diff;
503        }
504        else if (newcomp>0)
505        {
506          if ((olddiff==1) && (lastcomp<0))
507          {
508            notFound = FALSE;
509            newpos++;
510          }
511          else
512          {
513            newpos += diff;
514          }
515        }
516        else
517        {
518          notFound = FALSE;
519        }
520        lastcomp = newcomp;
521        if (diff==0) notFound=FALSE; /*hs*/
522      }
523      if (newpos<0) newpos = 0;
524      if (newpos>actpos) newpos = actpos;
525      while ((newpos<actpos) && (pComp_RevLex(id->m[i],id->m[(*result)[newpos]],nolex)==0))
526        newpos++;
527      for (j=actpos;j>newpos;j--)
528      {
529        (*result)[j] = (*result)[j-1];
530      }
531      (*result)[newpos] = i;
532      actpos++;
533    }
534  }
535  for (j=0;j<actpos;j++) (*result)[j]++;
536  return result;
537}
538
539/*2
540* concat the lists h1 and h2 without zeros
541*/
542ideal idSimpleAdd (ideal h1,ideal h2)
543{
544  int i,j,r,l;
545  ideal result;
546
547  if (h1==NULL) return idCopy(h2);
548  if (h2==NULL) return idCopy(h1);
549  j = IDELEMS(h1)-1;
550  while ((j >= 0) && (h1->m[j] == NULL)) j--;
551  i = IDELEMS(h2)-1;
552  while ((i >= 0) && (h2->m[i] == NULL)) i--;
553  r = si_max(h1->rank,h2->rank);
554  if (i+j==(-2))
555    return idInit(1,r);
556  else
557    result=idInit(i+j+2,r);
558  for (l=j; l>=0; l--)
559  {
560    result->m[l] = pCopy(h1->m[l]);
561  }
562  r = i+j+1;
563  for (l=i; l>=0; l--, r--)
564  {
565    result->m[r] = pCopy(h2->m[l]);
566  }
567  return result;
568}
569
[e070895]570/*2
[ded085]571* insert h2 into h1 (if h2 is not the zero polynomial)
572* return TRUE iff h2 was indeed inserted
[e070895]573*/
[ded085]574BOOLEAN idInsertPoly (ideal h1, poly h2)
[e070895]575{
[ded085]576  if (h2==NULL) return FALSE;
[e070895]577  int j = IDELEMS(h1)-1;
578  while ((j >= 0) && (h1->m[j] == NULL)) j--;
579  j++;
580  if (j==IDELEMS(h1))
581  {
582    pEnlargeSet(&(h1->m),IDELEMS(h1),16);
583    IDELEMS(h1)+=16;
584  }
585  h1->m[j]=h2;
[ded085]586  return TRUE;
[e070895]587}
588
[1a68d1d]589/*2
[2b3caae]590* insert h2 into h1 depending on the two boolean parameters:
591* - if zeroOk is true, then h2 will also be inserted when it is zero
592* - if duplicateOk is true, then h2 will also be inserted when it is
593*   already present in h1
[ded085]594* return TRUE iff h2 was indeed inserted
[1a68d1d]595*/
[2b3caae]596BOOLEAN idInsertPolyWithTests (ideal h1, const int validEntries,
597  const poly h2, const bool zeroOk, const bool duplicateOk)
[1a68d1d]598{
[2b3caae]599  if ((!zeroOk) && (h2 == NULL)) return FALSE;
600  if (!duplicateOk)
[1a68d1d]601  {
[2b3caae]602    bool h2FoundInH1 = false;
603    int i = 0;
604    while ((i < validEntries) && (!h2FoundInH1))
605    {
606      h2FoundInH1 = pEqualPolys(h1->m[i], h2);
607      i++;
608    }
609    if (h2FoundInH1) return FALSE;
[1a68d1d]610  }
[2b3caae]611  if (validEntries == IDELEMS(h1))
612  {
613    pEnlargeSet(&(h1->m), IDELEMS(h1), 16);
614    IDELEMS(h1) += 16;
615  }
616  h1->m[validEntries] = h2;
617  return TRUE;
[1a68d1d]618}
619
[35aab3]620/*2
621* h1 + h2
622*/
623ideal idAdd (ideal h1,ideal h2)
624{
625  ideal result = idSimpleAdd(h1,h2);
[10ea45f]626  idCompactify(result);
[35c62a9]627  return result;
[35aab3]628}
629
630/*2
631* h1 * h2
632*/
633ideal  idMult (ideal h1,ideal  h2)
634{
635  int i,j,k;
636  ideal  hh;
637
638  j = IDELEMS(h1);
639  while ((j > 0) && (h1->m[j-1] == NULL)) j--;
640  i = IDELEMS(h2);
641  while ((i > 0) && (h2->m[i-1] == NULL)) i--;
642  j = j * i;
643  if (j == 0)
644    hh = idInit(1,1);
645  else
646    hh=idInit(j,1);
647  if (h1->rank<h2->rank)
648    hh->rank = h2->rank;
649  else
650    hh->rank = h1->rank;
651  if (j==0) return hh;
652  k = 0;
653  for (i=0; i<IDELEMS(h1); i++)
654  {
655    if (h1->m[i] != NULL)
656    {
657      for (j=0; j<IDELEMS(h2); j++)
658      {
659        if (h2->m[j] != NULL)
660        {
661          hh->m[k] = ppMult_qq(h1->m[i],h2->m[j]);
662          k++;
663        }
664      }
665    }
666  }
667  {
[10ea45f]668    idCompactify(hh);
669    return hh;
[35aab3]670  }
671}
672
673/*2
674*returns true if h is the zero ideal
675*/
676BOOLEAN idIs0 (ideal h)
677{
678  int i;
679
680  if (h == NULL) return TRUE;
[9dd6270]681  i = IDELEMS(h)-1;
682  while ((i >= 0) && (h->m[i] == NULL))
[35aab3]683  {
684    i--;
685  }
[9dd6270]686  if (i < 0)
[35aab3]687    return TRUE;
688  else
689    return FALSE;
690}
691
692/*2
693* return the maximal component number found in any polynomial in s
694*/
695long idRankFreeModule (ideal s, ring lmRing, ring tailRing)
696{
697  if (s!=NULL)
698  {
699    int  j=0;
700
701    if (rRing_has_Comp(tailRing) && rRing_has_Comp(lmRing))
702    {
703      int  l=IDELEMS(s);
704      poly *p=s->m;
705      int  k;
706      for (; l != 0; l--)
707      {
708        if (*p!=NULL)
709        {
710          pp_Test(*p, lmRing, tailRing);
711          k = p_MaxComp(*p, lmRing, tailRing);
712          if (k>j) j = k;
713        }
714        p++;
715      }
716    }
717    return j;
718  }
719  return -1;
720}
721
722BOOLEAN idIsModule(ideal id, ring r)
723{
724  if (id != NULL && rRing_has_Comp(r))
725  {
726    int j, l = IDELEMS(id);
727    for (j=0; j<l; j++)
728    {
729      if (id->m[j] != NULL && p_GetComp(id->m[j], r) > 0) return TRUE;
730    }
731  }
732  return FALSE;
733}
734
735
736/*2
737*returns true if id is homogenous with respect to the aktual weights
738*/
739BOOLEAN idHomIdeal (ideal id, ideal Q)
740{
741  int i;
742  BOOLEAN     b;
743  if ((id == NULL) || (IDELEMS(id) == 0)) return TRUE;
744  i = 0;
745  b = TRUE;
746  while ((i < IDELEMS(id)) && b)
747  {
748    b = pIsHomogeneous(id->m[i]);
749    i++;
750  }
751  if ((b) && (Q!=NULL) && (IDELEMS(Q)>0))
752  {
753    i=0;
754    while ((i < IDELEMS(Q)) && b)
755    {
756      b = pIsHomogeneous(Q->m[i]);
757      i++;
758    }
759  }
760  return b;
761}
762
763/*2
764*the minimal index of used variables - 1
765*/
766int pLowVar (poly p)
767{
768  int k,l,lex;
769
770  if (p == NULL) return -1;
771
772  k = 32000;/*a very large dummy value*/
773  while (p != NULL)
774  {
775    l = 1;
776    lex = pGetExp(p,l);
777    while ((l < pVariables) && (lex == 0))
778    {
779      l++;
780      lex = pGetExp(p,l);
781    }
782    l--;
783    if (l < k) k = l;
784    pIter(p);
785  }
786  return k;
787}
788
789/*3
790*multiplies p with t (!cas) or  (t-1)
791*the index of t is:1, so we have to shift all variables
792*p is NOT in the actual ring, it has no t
793*/
794static poly pMultWithT (poly p,BOOLEAN cas)
795{
796  /*qp is the working pointer in p*/
797  /*result is the result, qresult is the working pointer*/
798  /*pp is p in the actual ring(shifted), qpp the working pointer*/
799  poly result,qp,pp;
800  poly qresult=NULL;
801  poly qpp=NULL;
802  int  i,j,lex;
803  number n;
804
805  pp = NULL;
806  result = NULL;
807  qp = p;
808  while (qp != NULL)
809  {
810    i = 0;
811    if (result == NULL)
812    {/*first monomial*/
813      result = pInit();
814      qresult = result;
815    }
816    else
817    {
818      qresult->next = pInit();
819      pIter(qresult);
820    }
821    for (j=pVariables-1; j>0; j--)
822    {
823      lex = pGetExp(qp,j);
824      pSetExp(qresult,j+1,lex);/*copy all variables*/
825    }
826    lex = pGetComp(qp);
827    pSetComp(qresult,lex);
828    n=nCopy(pGetCoeff(qp));
829    pSetCoeff0(qresult,n);
830    qresult->next = NULL;
831    pSetm(qresult);
832    /*qresult is now qp brought into the actual ring*/
833    if (cas)
834    { /*case: mult with t-1*/
835      pSetExp(qresult,1,0);
836      pSetm(qresult);
837      if (pp == NULL)
838      { /*first monomial*/
839        pp = pCopy(qresult);
840        qpp = pp;
841      }
842      else
843      {
844        qpp->next = pCopy(qresult);
845        pIter(qpp);
846      }
847      pGetCoeff(qpp)=nNeg(pGetCoeff(qpp));
848      /*now qpp contains -1*qp*/
849    }
850    pSetExp(qresult,1,1);/*this is mult. by t*/
851    pSetm(qresult);
852    pIter(qp);
853  }
854  /*
855  *now p is processed:
856  *result contains t*p
857  * if cas: pp contains -1*p (in the new ring)
858  */
859  if (cas)  qresult->next = pp;
860  /*  else      qresult->next = NULL;*/
861  return result;
862}
863
864/*2
865* verschiebt die Indizees der Modulerzeugenden um i
866*/
867void pShift (poly * p,int i)
868{
869  poly qp1 = *p,qp2 = *p;/*working pointers*/
870  int     j = pMaxComp(*p),k = pMinComp(*p);
871
872  if (j+i < 0) return ;
873  while (qp1 != NULL)
874  {
875    if ((pGetComp(qp1)+i > 0) || ((j == -i) && (j == k)))
876    {
[959263]877      pAddComp(qp1,i);
[35aab3]878      pSetmComp(qp1);
879      qp2 = qp1;
880      pIter(qp1);
881    }
882    else
883    {
884      if (qp2 == *p)
885      {
886        pIter(*p);
[fb82895]887        pLmDelete(&qp2);
[35aab3]888        qp2 = *p;
889        qp1 = *p;
890      }
891      else
892      {
893        qp2->next = qp1->next;
[fb82895]894        if (qp1!=NULL) pLmDelete(&qp1);
[35aab3]895        qp1 = qp2->next;
896      }
897    }
898  }
899}
900
901/*2
902*initialized a field with r numbers between beg and end for the
903*procedure idNextChoise
904*/
905void idInitChoise (int r,int beg,int end,BOOLEAN  *endch,int * choise)
906{
907  /*returns the first choise of r numbers between beg and end*/
908  int i;
909  for (i=0; i<r; i++)
910  {
911    choise[i] = 0;
912  }
913  if (r <= end-beg+1)
914    for (i=0; i<r; i++)
915    {
916      choise[i] = beg+i;
917    }
918  if (r > end-beg+1)
919    *endch = TRUE;
920  else
921    *endch = FALSE;
922}
923
924/*2
925*returns the next choise of r numbers between beg and end
926*/
927void idGetNextChoise (int r,int end,BOOLEAN  *endch,int * choise)
928{
929  int i = r-1,j;
930  while ((i >= 0) && (choise[i] == end))
931  {
932    i--;
933    end--;
934  }
935  if (i == -1)
936    *endch = TRUE;
937  else
938  {
939    choise[i]++;
940    for (j=i+1; j<r; j++)
941    {
942      choise[j] = choise[i]+j-i;
943    }
944    *endch = FALSE;
945  }
946}
947
948/*2
949*takes the field choise of d numbers between beg and end, cancels the t-th
950*entree and searches for the ordinal number of that d-1 dimensional field
951* w.r.t. the algorithm of construction
952*/
953int idGetNumberOfChoise(int t, int d, int begin, int end, int * choise)
954{
955  int * localchoise,i,result=0;
956  BOOLEAN b=FALSE;
957
958  if (d<=1) return 1;
959  localchoise=(int*)omAlloc((d-1)*sizeof(int));
960  idInitChoise(d-1,begin,end,&b,localchoise);
961  while (!b)
962  {
963    result++;
964    i = 0;
965    while ((i<t) && (localchoise[i]==choise[i])) i++;
966    if (i>=t)
967    {
968      i = t+1;
969      while ((i<d) && (localchoise[i-1]==choise[i])) i++;
970      if (i>=d)
[f71e8c5]971      {
972        omFreeSize((ADDRESS)localchoise,(d-1)*sizeof(int));
973        return result;
[35aab3]974      }
975    }
[f71e8c5]976    idGetNextChoise(d-1,end,&b,localchoise);
[35aab3]977  }
[f71e8c5]978  omFreeSize((ADDRESS)localchoise,(d-1)*sizeof(int));
979  return 0;
[35aab3]980}
981
982/*2
[f71e8c5]983*computes the binomial coefficient
[35aab3]984*/
[f71e8c5]985int binom (int n,int r)
986{
987  int i,result;
[35aab3]988
[f71e8c5]989  if (r==0) return 1;
990  if (n-r<r) return binom(n,n-r);
991  result = n-r+1;
992  for (i=2;i<=r;i++)
[35aab3]993  {
[f71e8c5]994    result *= n-r+i;
995    if (result<0)
[35aab3]996    {
[f71e8c5]997      WarnS("overflow in binomials");
998      return 0;
[35aab3]999    }
[f71e8c5]1000    result /= i;
[35aab3]1001  }
[f71e8c5]1002  return result;
[35aab3]1003}
[f71e8c5]1004
[35aab3]1005/*2
[f71e8c5]1006*the free module of rank i
[35aab3]1007*/
[f71e8c5]1008ideal idFreeModule (int i)
[35aab3]1009{
[f71e8c5]1010  int j;
1011  ideal h;
1012
1013  h=idInit(i,i);
1014  for (j=0; j<i; j++)
[35aab3]1015  {
[f71e8c5]1016    h->m[j] = pOne();
1017    pSetComp(h->m[j],j+1);
1018    pSetmComp(h->m[j]);
[35aab3]1019  }
[f71e8c5]1020  return h;
1021}
[35aab3]1022
[f71e8c5]1023ideal idSectWithElim (ideal h1,ideal h2)
1024// does not destroy h1,h2
1025{
1026  if (TEST_OPT_PROT) PrintS("intersect by elimination method\n");
1027  assume(!idIs0(h1));
1028  assume(!idIs0(h2));
1029  assume(IDELEMS(h1)<=IDELEMS(h2));
1030  assume(idRankFreeModule(h1)==0);
1031  assume(idRankFreeModule(h2)==0);
1032  // add a new variable:
1033  int j;
1034  ring origRing=currRing;
1035  ring r=rCopy0(origRing);
1036  r->N++;
1037  r->block0[0]=1;
1038  r->block1[0]= r->N;
1039  omFree(r->order);
1040  r->order=(int*)omAlloc0(3*sizeof(int*));
1041  r->order[0]=ringorder_dp;
1042  r->order[1]=ringorder_C;
1043  char **names=(char**)omAlloc0(rVar(r) * sizeof(char_ptr));
1044  for (j=0;j<r->N-1;j++) names[j]=r->names[j];
1045  names[r->N-1]=omStrDup("@");
1046  omFree(r->names);
1047  r->names=names;
1048  rComplete(r,TRUE);
1049  // fetch h1, h2
1050  ideal h;
1051  h1=idrCopyR(h1,origRing,r);
1052  h2=idrCopyR(h2,origRing,r);
1053  // switch to temp. ring r
1054  rChangeCurrRing(r);
1055  // create 1-t, t
1056  poly omt=pOne();
1057  pSetExp(omt,r->N,1);
1058  poly t=pCopy(omt);
1059  pSetm(omt);
1060  omt=pNeg(omt);
1061  omt=pAdd(omt,pOne());
1062  // compute (1-t)*h1
1063  h1=(ideal)mpMultP((matrix)h1,omt);
1064  // compute t*h2
1065  h2=(ideal)mpMultP((matrix)h2,pCopy(t));
1066  // (1-t)h1 + t*h2
1067  h=idInit(IDELEMS(h1)+IDELEMS(h2),1);
1068  int l;
1069  for (l=IDELEMS(h1)-1; l>=0; l--)
[35aab3]1070  {
[f71e8c5]1071    h->m[l] = h1->m[l];  h1->m[l]=NULL;
[35aab3]1072  }
[f71e8c5]1073  j=IDELEMS(h1);
1074  for (l=IDELEMS(h2)-1; l>=0; l--)
[35aab3]1075  {
[f71e8c5]1076    h->m[l+j] = h2->m[l];  h2->m[l]=NULL;
[35aab3]1077  }
[f71e8c5]1078  idDelete(&h1);
1079  idDelete(&h2);
1080  // eliminate t:
1081
1082  ideal res=idElimination(h,t);
[3d0808]1083  // cleanup
[f71e8c5]1084  idDelete(&h);
1085  res=idrMoveR(res,r,origRing);
1086  rChangeCurrRing(origRing);
1087  rKill(r);
1088  return res;
[35aab3]1089}
1090
1091/*2
1092*computes recursively all monomials of a certain degree
1093*in every step the actvar-th entry in the exponential
1094*vector is incremented and the other variables are
1095*computed by recursive calls of makemonoms
1096*if the last variable is reached, the difference to the
1097*degree is computed directly
1098*vars is the number variables
1099*actvar is the actual variable to handle
1100*deg is the degree of the monomials to compute
1101*monomdeg is the actual degree of the monomial in consideration
1102*/
1103static void makemonoms(int vars,int actvar,int deg,int monomdeg)
1104{
1105  poly p;
1106  int i=0;
1107
1108  if ((idpowerpoint == 0) && (actvar ==1))
1109  {
1110    idpower[idpowerpoint] = pOne();
1111    monomdeg = 0;
1112  }
1113  while (i<=deg)
1114  {
1115    if (deg == monomdeg)
1116    {
1117      pSetm(idpower[idpowerpoint]);
1118      idpowerpoint++;
1119      return;
1120    }
1121    if (actvar == vars)
1122    {
1123      pSetExp(idpower[idpowerpoint],actvar,deg-monomdeg);
1124      pSetm(idpower[idpowerpoint]);
1125      pTest(idpower[idpowerpoint]);
1126      idpowerpoint++;
1127      return;
1128    }
1129    else
1130    {
1131      p = pCopy(idpower[idpowerpoint]);
1132      makemonoms(vars,actvar+1,deg,monomdeg);
1133      idpower[idpowerpoint] = p;
1134    }
1135    monomdeg++;
1136    pSetExp(idpower[idpowerpoint],actvar,pGetExp(idpower[idpowerpoint],actvar)+1);
1137    pSetm(idpower[idpowerpoint]);
1138    pTest(idpower[idpowerpoint]);
1139    i++;
1140  }
1141}
1142
1143/*2
1144*returns the deg-th power of the maximal ideal of 0
1145*/
1146ideal idMaxIdeal(int deg)
1147{
1148  if (deg < 0)
1149  {
1150    WarnS("maxideal: power must be non-negative");
1151  }
1152  if (deg < 1)
1153  {
1154    ideal I=idInit(1,1);
1155    I->m[0]=pOne();
1156    return I;
1157  }
1158  if (deg == 1)
1159  {
1160    return idMaxIdeal();
1161  }
1162
1163  int vars = currRing->N;
1164  int i = binom(vars+deg-1,deg);
1165  if (i<=0) return idInit(1,1);
1166  ideal id=idInit(i,1);
1167  idpower = id->m;
1168  idpowerpoint = 0;
1169  makemonoms(vars,1,deg,0);
1170  idpower = NULL;
1171  idpowerpoint = 0;
1172  return id;
1173}
1174
1175/*2
1176*computes recursively all generators of a certain degree
1177*of the ideal "givenideal"
1178*elms is the number elements in the given ideal
1179*actelm is the actual element to handle
1180*deg is the degree of the power to compute
1181*gendeg is the actual degree of the generator in consideration
1182*/
1183static void makepotence(int elms,int actelm,int deg,int gendeg)
1184{
1185  poly p;
1186  int i=0;
1187
1188  if ((idpowerpoint == 0) && (actelm ==1))
1189  {
1190    idpower[idpowerpoint] = pOne();
1191    gendeg = 0;
1192  }
1193  while (i<=deg)
1194  {
1195    if (deg == gendeg)
1196    {
1197      idpowerpoint++;
1198      return;
1199    }
1200    if (actelm == elms)
1201    {
1202      p=pPower(pCopy(givenideal[actelm-1]),deg-gendeg);
1203      idpower[idpowerpoint]=pMult(idpower[idpowerpoint],p);
1204      idpowerpoint++;
1205      return;
1206    }
1207    else
1208    {
1209      p = pCopy(idpower[idpowerpoint]);
1210      makepotence(elms,actelm+1,deg,gendeg);
1211      idpower[idpowerpoint] = p;
1212    }
1213    gendeg++;
1214    idpower[idpowerpoint]=pMult(idpower[idpowerpoint],pCopy(givenideal[actelm-1]));
1215    i++;
1216  }
1217}
1218
1219/*2
1220*returns the deg-th power of the ideal gid
1221*/
1222//ideal idPower(ideal gid,int deg)
1223//{
1224//  int i;
1225//  ideal id;
1226//
1227//  if (deg < 1) deg = 1;
1228//  i = binom(IDELEMS(gid)+deg-1,deg);
1229//  id=idInit(i,1);
1230//  idpower = id->m;
1231//  givenideal = gid->m;
1232//  idpowerpoint = 0;
1233//  makepotence(IDELEMS(gid),1,deg,0);
1234//  idpower = NULL;
1235//  givenideal = NULL;
1236//  idpowerpoint = 0;
1237//  return id;
1238//}
1239static void idNextPotence(ideal given, ideal result,
1240  int begin, int end, int deg, int restdeg, poly ap)
1241{
1242  poly p;
1243  int i;
1244
1245  p = pPower(pCopy(given->m[begin]),restdeg);
1246  i = result->nrows;
1247  result->m[i] = pMult(pCopy(ap),p);
1248//PrintS(".");
1249  (result->nrows)++;
1250  if (result->nrows >= IDELEMS(result))
1251  {
1252    pEnlargeSet(&(result->m),IDELEMS(result),16);
1253    IDELEMS(result) += 16;
1254  }
1255  if (begin == end) return;
1256  for (i=restdeg-1;i>0;i--)
1257  {
1258    p = pPower(pCopy(given->m[begin]),i);
1259    p = pMult(pCopy(ap),p);
1260    idNextPotence(given, result, begin+1, end, deg, restdeg-i, p);
1261    pDelete(&p);
1262  }
1263  idNextPotence(given, result, begin+1, end, deg, restdeg, ap);
1264}
1265
1266ideal idPower(ideal given,int exp)
1267{
1268  ideal result,temp;
1269  poly p1;
1270  int i;
1271
1272  if (idIs0(given)) return idInit(1,1);
1273  temp = idCopy(given);
1274  idSkipZeroes(temp);
1275  i = binom(IDELEMS(temp)+exp-1,exp);
1276  result = idInit(i,1);
1277  result->nrows = 0;
1278//Print("ideal contains %d elements\n",i);
1279  p1=pOne();
1280  idNextPotence(temp,result,0,IDELEMS(temp)-1,exp,exp,p1);
1281  pDelete(&p1);
1282  idDelete(&temp);
1283  result->nrows = 1;
1284  idDelEquals(result);
[ff2fd1]1285  idSkipZeroes(result);
[35aab3]1286  return result;
1287}
1288
[0a64b14]1289/*2
1290* compute the which-th ar-minor of the matrix a
1291*/
1292poly idMinor(matrix a, int ar, unsigned long which, ideal R)
1293{
1294  int     i,j,k,size;
1295  unsigned long curr;
1296  int *rowchoise,*colchoise;
1297  BOOLEAN rowch,colch;
1298  ideal result;
1299  matrix tmp;
1300  poly p,q;
1301
1302  i = binom(a->rows(),ar);
1303  j = binom(a->cols(),ar);
1304
1305  rowchoise=(int *)omAlloc(ar*sizeof(int));
1306  colchoise=(int *)omAlloc(ar*sizeof(int));
1307  if ((i>512) || (j>512) || (i*j >512)) size=512;
1308  else size=i*j;
1309  result=idInit(size,1);
1310  tmp=mpNew(ar,ar);
1311  k = 0; /* the index in result*/
1312  curr = 0; /* index of current minor */
1313  idInitChoise(ar,1,a->rows(),&rowch,rowchoise);
1314  while (!rowch)
1315  {
1316    idInitChoise(ar,1,a->cols(),&colch,colchoise);
1317    while (!colch)
1318    {
1319      if (curr == which)
1320      {
1321        for (i=1; i<=ar; i++)
1322        {
1323          for (j=1; j<=ar; j++)
1324          {
1325            MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
1326          }
1327        }
1328        p = mpDetBareiss(tmp);
1329        if (p!=NULL)
1330        {
1331          if (R!=NULL)
1332          {
1333            q = p;
1334            p = kNF(R,currQuotient,q);
1335            pDelete(&q);
1336          }
1337          /*delete the matrix tmp*/
1338          for (i=1; i<=ar; i++)
1339          {
1340            for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
1341          }
1342          idDelete((ideal*)&tmp);
1343          omFreeSize((ADDRESS)rowchoise,ar*sizeof(int));
1344          omFreeSize((ADDRESS)colchoise,ar*sizeof(int));
1345          return (p);
1346        }
1347      }
1348      curr++;
1349      idGetNextChoise(ar,a->cols(),&colch,colchoise);
1350    }
1351    idGetNextChoise(ar,a->rows(),&rowch,rowchoise);
1352  }
1353  return (poly) 1;
1354}
1355
[35aab3]1356#ifdef WITH_OLD_MINOR
1357/*2
1358* compute all ar-minors of the matrix a
1359*/
1360ideal idMinors(matrix a, int ar, ideal R)
1361{
1362  int     i,j,k,size;
1363  int *rowchoise,*colchoise;
1364  BOOLEAN rowch,colch;
1365  ideal result;
1366  matrix tmp;
1367  poly p,q;
1368
1369  i = binom(a->rows(),ar);
1370  j = binom(a->cols(),ar);
1371
1372  rowchoise=(int *)omAlloc(ar*sizeof(int));
1373  colchoise=(int *)omAlloc(ar*sizeof(int));
1374  if ((i>512) || (j>512) || (i*j >512)) size=512;
1375  else size=i*j;
1376  result=idInit(size,1);
1377  tmp=mpNew(ar,ar);
1378  k = 0; /* the index in result*/
1379  idInitChoise(ar,1,a->rows(),&rowch,rowchoise);
1380  while (!rowch)
1381  {
1382    idInitChoise(ar,1,a->cols(),&colch,colchoise);
1383    while (!colch)
1384    {
1385      for (i=1; i<=ar; i++)
1386      {
1387        for (j=1; j<=ar; j++)
1388        {
1389          MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
1390        }
1391      }
1392      p = mpDetBareiss(tmp);
1393      if (p!=NULL)
1394      {
1395        if (R!=NULL)
1396        {
1397          q = p;
1398          p = kNF(R,currQuotient,q);
1399          pDelete(&q);
1400        }
1401        if (p!=NULL)
1402        {
1403          if (k>=size)
1404          {
1405            pEnlargeSet(&result->m,size,32);
1406            size += 32;
1407          }
1408          result->m[k] = p;
1409          k++;
1410        }
1411      }
1412      idGetNextChoise(ar,a->cols(),&colch,colchoise);
1413    }
1414    idGetNextChoise(ar,a->rows(),&rowch,rowchoise);
1415  }
1416  /*delete the matrix tmp*/
1417  for (i=1; i<=ar; i++)
1418  {
1419    for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
1420  }
1421  idDelete((ideal*)&tmp);
1422  if (k==0)
1423  {
1424    k=1;
1425    result->m[0]=NULL;
1426  }
1427  omFreeSize((ADDRESS)rowchoise,ar*sizeof(int));
1428  omFreeSize((ADDRESS)colchoise,ar*sizeof(int));
1429  pEnlargeSet(&result->m,size,k-size);
1430  IDELEMS(result) = k;
1431  return (result);
1432}
1433#else
1434/*2
1435* compute all ar-minors of the matrix a
1436* the caller of mpRecMin
1437* the elements of the result are not in R (if R!=NULL)
1438*/
1439ideal idMinors(matrix a, int ar, ideal R)
1440{
1441  int elems=0;
1442  int r=a->nrows,c=a->ncols;
1443  int i;
1444  matrix b;
1445  ideal result,h;
1446  ring origR;
[ab453da]1447  ring tmpR;
[0b5e3d]1448  long bound;
[35aab3]1449
1450  if((ar<=0) || (ar>r) || (ar>c))
1451  {
1452    Werror("%d-th minor, matrix is %dx%d",ar,r,c);
1453    return NULL;
1454  }
1455  h = idMatrix2Module(mpCopy(a));
1456  bound = smExpBound(h,c,r,ar);
1457  idDelete(&h);
[ab453da]1458  tmpR=smRingChange(&origR,bound);
[35aab3]1459  b = mpNew(r,c);
1460  for (i=r*c-1;i>=0;i--)
1461  {
1462    if (a->m[i])
1463      b->m[i] = prCopyR(a->m[i],origR);
1464  }
[7012d0]1465  if (R!=NULL)
1466  {
1467    R = idrCopyR(R,origR);
[3060a7]1468    //if (ar>1) // otherwise done in mpMinorToResult
1469    //{
1470    //  matrix bb=(matrix)kNF(R,currQuotient,(ideal)b);
1471    //  bb->rank=b->rank; bb->nrows=b->nrows; bb->ncols=b->ncols;
1472    //  idDelete((ideal*)&b); b=bb;
1473    //}
[7012d0]1474  }
[35aab3]1475  result=idInit(32,1);
1476  if(ar>1) mpRecMin(ar-1,result,elems,b,r,c,NULL,R);
1477  else mpMinorToResult(result,elems,b,r,c,R);
1478  idDelete((ideal *)&b);
[7012d0]1479  if (R!=NULL) idDelete(&R);
[35aab3]1480  idSkipZeroes(result);
1481  rChangeCurrRing(origR);
[ab453da]1482  result = idrMoveR(result,tmpR);
1483  smKillModifiedRing(tmpR);
[35aab3]1484  idTest(result);
1485  return result;
1486}
1487#endif
1488
1489/*2
1490*skips all zeroes and double elements, searches also for units
1491*/
[10ea45f]1492void idCompactify(ideal id)
[35aab3]1493{
1494  int i,j;
1495  BOOLEAN b=FALSE;
1496
1497  i = IDELEMS(id)-1;
1498  while ((! b) && (i>=0))
1499  {
1500    b=pIsUnit(id->m[i]);
1501    i--;
1502  }
1503  if (b)
1504  {
[10ea45f]1505    for(i=IDELEMS(id)-1;i>=0;i--) pDelete(&id->m[i]);
1506    id->m[0]=pOne();
[35aab3]1507  }
1508  else
1509  {
[10ea45f]1510    idDelMultiples(id);
[35aab3]1511  }
[962de7]1512  idSkipZeroes(id);
[35aab3]1513}
1514
1515/*2
1516*returns TRUE if id1 is a submodule of id2
1517*/
1518BOOLEAN idIsSubModule(ideal id1,ideal id2)
1519{
1520  int i;
1521  poly p;
1522
1523  if (idIs0(id1)) return TRUE;
1524  for (i=0;i<IDELEMS(id1);i++)
1525  {
1526    if (id1->m[i] != NULL)
1527    {
1528      p = kNF(id2,currQuotient,id1->m[i]);
1529      if (p != NULL)
1530      {
1531        pDelete(&p);
1532        return FALSE;
1533      }
1534    }
1535  }
1536  return TRUE;
1537}
1538
1539/*2
1540* returns the ideals of initial terms
1541*/
1542ideal idHead(ideal h)
1543{
1544  ideal m = idInit(IDELEMS(h),h->rank);
1545  int i;
1546
1547  for (i=IDELEMS(h)-1;i>=0; i--)
1548  {
1549    if (h->m[i]!=NULL) m->m[i]=pHead(h->m[i]);
1550  }
1551  return m;
1552}
1553
1554ideal idHomogen(ideal h, int varnum)
1555{
1556  ideal m = idInit(IDELEMS(h),h->rank);
1557  int i;
1558
1559  for (i=IDELEMS(h)-1;i>=0; i--)
1560  {
1561    m->m[i]=pHomogen(h->m[i],varnum);
1562  }
1563  return m;
1564}
1565
1566/*------------------type conversions----------------*/
1567ideal idVec2Ideal(poly vec)
1568{
1569   ideal result=idInit(1,1);
1570   omFree((ADDRESS)result->m);
1571   result->m=NULL; // remove later
1572   pVec2Polys(vec, &(result->m), &(IDELEMS(result)));
1573   return result;
1574}
1575
[ca3e7b]1576#define NEW_STUFF
[35aab3]1577#ifndef NEW_STUFF
1578// converts mat to module, destroys mat
1579ideal idMatrix2Module(matrix mat)
1580{
1581  int mc=MATCOLS(mat);
1582  int mr=MATROWS(mat);
1583  ideal result = idInit(si_max(mc,1),si_max(mr,1));
1584  int i,j;
1585  poly h;
1586
1587  for(j=0;j<mc /*MATCOLS(mat)*/;j++) /* j is also index in result->m */
1588  {
1589    for (i=1;i<=mr /*MATROWS(mat)*/;i++)
1590    {
1591      h = MATELEM(mat,i,j+1);
1592      if (h!=NULL)
1593      {
1594        MATELEM(mat,i,j+1)=NULL;
1595        pSetCompP(h,i);
1596        result->m[j] = pAdd(result->m[j],h);
1597      }
1598    }
1599  }
1600  // obachman: need to clean this up
1601  idDelete((ideal*) &mat);
1602  return result;
1603}
1604#else
1605
[2ad10e9]1606#include "sbuckets.h"
[35aab3]1607
1608// converts mat to module, destroys mat
1609ideal idMatrix2Module(matrix mat)
1610{
1611  int mc=MATCOLS(mat);
1612  int mr=MATROWS(mat);
1613  ideal result = idInit(si_max(mc,1),si_max(mr,1));
1614  int i,j, l;
1615  poly h;
1616  poly p;
[cbeafc2]1617  sBucket_pt bucket = sBucketCreate(currRing);
[35aab3]1618
1619  for(j=0;j<mc /*MATCOLS(mat)*/;j++) /* j is also index in result->m */
1620  {
1621    for (i=1;i<=mr /*MATROWS(mat)*/;i++)
1622    {
1623      h = MATELEM(mat,i,j+1);
1624      if (h!=NULL)
1625      {
[ca3e7b]1626        l=pLength(h);
[35aab3]1627        MATELEM(mat,i,j+1)=NULL;
[cbeafc2]1628        p_SetCompP(h,i, currRing);
[35aab3]1629        sBucket_Merge_p(bucket, h, l);
1630      }
1631    }
1632    sBucketClearMerge(bucket, &(result->m[j]), &l);
1633  }
[cbeafc2]1634  sBucketDestroy(&bucket);
[35aab3]1635
1636  // obachman: need to clean this up
1637  idDelete((ideal*) &mat);
1638  return result;
1639}
1640#endif
1641
1642/*2
1643* converts a module into a matrix, destroyes the input
1644*/
1645matrix idModule2Matrix(ideal mod)
1646{
1647  matrix result = mpNew(mod->rank,IDELEMS(mod));
1648  int i,cp;
1649  poly p,h;
1650
1651  for(i=0;i<IDELEMS(mod);i++)
1652  {
[d0164d9]1653    p=pReverse(mod->m[i]);
[35aab3]1654    mod->m[i]=NULL;
1655    while (p!=NULL)
1656    {
1657      h=p;
1658      pIter(p);
1659      pNext(h)=NULL;
1660//      cp = si_max(1,pGetComp(h));     // if used for ideals too
1661      cp = pGetComp(h);
1662      pSetComp(h,0);
1663      pSetmComp(h);
1664#ifdef TEST
1665      if (cp>mod->rank)
1666      {
[6867f5]1667        Print("## inv. rank %ld -> %d\n",mod->rank,cp);
[35aab3]1668        int k,l,o=mod->rank;
1669        mod->rank=cp;
1670        matrix d=mpNew(mod->rank,IDELEMS(mod));
1671        for (l=1; l<=o; l++)
1672        {
1673          for (k=1; k<=IDELEMS(mod); k++)
1674          {
1675            MATELEM(d,l,k)=MATELEM(result,l,k);
1676            MATELEM(result,l,k)=NULL;
1677          }
1678        }
1679        idDelete((ideal *)&result);
1680        result=d;
1681      }
1682#endif
1683      MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
1684    }
1685  }
1686  // obachman 10/99: added the following line, otherwise memory leack!
1687  idDelete(&mod);
1688  return result;
1689}
1690
1691matrix idModule2formatedMatrix(ideal mod,int rows, int cols)
1692{
1693  matrix result = mpNew(rows,cols);
1694  int i,cp,r=idRankFreeModule(mod),c=IDELEMS(mod);
1695  poly p,h;
1696
1697  if (r>rows) r = rows;
1698  if (c>cols) c = cols;
1699  for(i=0;i<c;i++)
1700  {
[bafaec0]1701    p=pReverse(mod->m[i]);
[35aab3]1702    mod->m[i]=NULL;
1703    while (p!=NULL)
1704    {
1705      h=p;
1706      pIter(p);
1707      pNext(h)=NULL;
1708      cp = pGetComp(h);
1709      if (cp<=r)
1710      {
1711        pSetComp(h,0);
1712        pSetmComp(h);
1713        MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
1714      }
1715      else
1716        pDelete(&h);
1717    }
1718  }
1719  idDelete(&mod);
1720  return result;
1721}
1722
1723/*2
1724* substitute the n-th variable by the monomial e in id
1725* destroy id
1726*/
1727ideal  idSubst(ideal id, int n, poly e)
1728{
1729  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
1730  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
1731
1732  res->rank = id->rank;
1733  for(k--;k>=0;k--)
1734  {
1735    res->m[k]=pSubst(id->m[k],n,e);
1736    id->m[k]=NULL;
1737  }
1738  idDelete(&id);
1739  return res;
1740}
1741
1742BOOLEAN idHomModule(ideal m, ideal Q, intvec **w)
1743{
1744  if (w!=NULL) *w=NULL;
1745  if ((Q!=NULL) && (!idHomIdeal(Q,NULL))) return FALSE;
[43ebb1]1746  if (idIs0(m))
1747  {
[a12776]1748    if (w!=NULL) (*w)=new intvec(m->rank);
[43ebb1]1749    return TRUE;
1750  }
[35aab3]1751
[4e63600]1752  long cmax=1,order=0,ord,* diff,diffmin=32000;
1753  int *iscom;
1754  int i,j;
[35aab3]1755  poly p=NULL;
[1f5db38]1756  pFDegProc d;
1757  if (pLexOrder && (currRing->order[0]==ringorder_lp))
[99bdcf]1758     d=p_Totaldegree;
[bead81]1759  else
[1f5db38]1760     d=pFDeg;
[35aab3]1761  int length=IDELEMS(m);
1762  polyset P=m->m;
1763  polyset F=(polyset)omAlloc(length*sizeof(poly));
1764  for (i=length-1;i>=0;i--)
1765  {
1766    p=F[i]=P[i];
[4e63600]1767    cmax=si_max(cmax,(long)pMaxComp(p));
[35aab3]1768  }
[4e63600]1769  cmax++;
1770  diff = (long *)omAlloc0(cmax*sizeof(long));
[35aab3]1771  if (w!=NULL) *w=new intvec(cmax-1);
1772  iscom = (int *)omAlloc0(cmax*sizeof(int));
1773  i=0;
1774  while (i<=length)
1775  {
1776    if (i<length)
1777    {
1778      p=F[i];
[4e63600]1779      while ((p!=NULL) && (iscom[pGetComp(p)]==0)) pIter(p);
[35aab3]1780    }
1781    if ((p==NULL) && (i<length))
1782    {
1783      i++;
1784    }
1785    else
1786    {
[4e63600]1787      if (p==NULL) /* && (i==length) */
[35aab3]1788      {
1789        i=0;
1790        while ((i<length) && (F[i]==NULL)) i++;
1791        if (i>=length) break;
1792        p = F[i];
1793      }
[1f5db38]1794      //if (pLexOrder && (currRing->order[0]==ringorder_lp))
1795      //  order=pTotaldegree(p);
1796      //else
[35aab3]1797      //  order = p->order;
[1f5db38]1798      //  order = pFDeg(p,currRing);
1799      order = d(p,currRing) +diff[pGetComp(p)];
1800      //order += diff[pGetComp(p)];
[35aab3]1801      p = F[i];
1802//Print("Actual p=F[%d]: ",i);pWrite(p);
1803      F[i] = NULL;
1804      i=0;
1805    }
1806    while (p!=NULL)
1807    {
[4e63600]1808      if (pLexOrder && (currRing->order[0]==ringorder_lp))
1809        ord=pTotaldegree(p);
1810      else
[35aab3]1811      //  ord = p->order;
[4e63600]1812        ord = pFDeg(p,currRing);
1813      if (iscom[pGetComp(p)]==0)
[35aab3]1814      {
1815        diff[pGetComp(p)] = order-ord;
1816        iscom[pGetComp(p)] = 1;
1817/*
1818*PrintS("new diff: ");
1819*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
1820*PrintLn();
1821*PrintS("new iscom: ");
1822*for (j=0;j<cmax;j++) Print("%d ",iscom[j]);
1823*PrintLn();
1824*Print("new set %d, order %d, ord %d, diff %d\n",pGetComp(p),order,ord,diff[pGetComp(p)]);
1825*/
1826      }
1827      else
1828      {
1829/*
1830*PrintS("new diff: ");
1831*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
1832*PrintLn();
1833*Print("order %d, ord %d, diff %d\n",order,ord,diff[pGetComp(p)]);
1834*/
[4e63600]1835        if (order != (ord+diff[pGetComp(p)]))
[35aab3]1836        {
1837          omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
[4e63600]1838          omFreeSize((ADDRESS) diff,cmax*sizeof(long));
[35aab3]1839          omFreeSize((ADDRESS) F,length*sizeof(poly));
1840          delete *w;*w=NULL;
1841          return FALSE;
1842        }
1843      }
1844      pIter(p);
1845    }
1846  }
1847  omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
1848  omFreeSize((ADDRESS) F,length*sizeof(poly));
[4e63600]1849  for (i=1;i<cmax;i++) (**w)[i-1]=(int)(diff[i]);
[35aab3]1850  for (i=1;i<cmax;i++)
1851  {
1852    if (diff[i]<diffmin) diffmin=diff[i];
1853  }
1854  if (w!=NULL)
1855  {
1856    for (i=1;i<cmax;i++)
1857    {
[4e63600]1858      (**w)[i-1]=(int)(diff[i]-diffmin);
[35aab3]1859    }
1860  }
[4e63600]1861  omFreeSize((ADDRESS) diff,cmax*sizeof(long));
[35aab3]1862  return TRUE;
1863}
1864
[30b8381]1865BOOLEAN idTestHomModule(ideal m, ideal Q, intvec *w)
1866{
1867  if ((Q!=NULL) && (!idHomIdeal(Q,NULL)))  { PrintS(" Q not hom\n"); return FALSE;}
1868  if (idIs0(m)) return TRUE;
1869
1870  int cmax=-1;
1871  int i;
1872  poly p=NULL;
1873  int length=IDELEMS(m);
1874  polyset P=m->m;
1875  for (i=length-1;i>=0;i--)
1876  {
1877    p=P[i];
1878    if (p!=NULL) cmax=si_max(cmax,(int)pMaxComp(p)+1);
1879  }
[8324cc]1880  if (w != NULL)
[30b8381]1881  if (w->length()+1 < cmax)
[4d13e7]1882  {
[30b8381]1883    // Print("length: %d - %d \n", w->length(),cmax);
1884    return FALSE;
1885  }
[8324cc]1886
1887  if(w!=NULL)
1888    pSetModDeg(w);
1889
[30b8381]1890  for (i=length-1;i>=0;i--)
1891  {
1892    p=P[i];
1893    poly q=p;
1894    if (p!=NULL)
1895    {
[b130fb]1896      int d=pFDeg(p,currRing);
[30b8381]1897      loop
1898      {
1899        pIter(p);
1900        if (p==NULL) break;
[4d13e7]1901        if (d!=pFDeg(p,currRing))
[30b8381]1902        {
[4d13e7]1903          //pWrite(q); wrp(p); Print(" -> %d - %d\n",d,pFDeg(p,currRing));
[8324cc]1904          if(w!=NULL)
1905            pSetModDeg(NULL);
[4d13e7]1906          return FALSE;
[30b8381]1907        }
1908      }
1909    }
1910  }
[cbc7e3]1911
[8324cc]1912  if(w!=NULL)
1913    pSetModDeg(NULL);
[cbc7e3]1914
[30b8381]1915  return TRUE;
1916}
1917
[35aab3]1918ideal idJet(ideal i,int d)
1919{
1920  ideal r=idInit((i->nrows)*(i->ncols),i->rank);
1921  r->nrows = i-> nrows;
1922  r->ncols = i-> ncols;
1923  //r->rank = i-> rank;
1924  int k;
1925  for(k=(i->nrows)*(i->ncols)-1;k>=0; k--)
1926  {
1927    r->m[k]=ppJet(i->m[k],d);
1928  }
1929  return r;
1930}
1931
1932ideal idJetW(ideal i,int d, intvec * iv)
1933{
1934  ideal r=idInit(IDELEMS(i),i->rank);
1935  if (ecartWeights!=NULL)
1936  {
1937    WerrorS("cannot compute weighted jets now");
1938  }
1939  else
1940  {
1941    short *w=iv2array(iv);
1942    int k;
1943    for(k=0; k<IDELEMS(i); k++)
1944    {
1945      r->m[k]=ppJetW(i->m[k],d,w);
1946    }
1947    omFreeSize((ADDRESS)w,(pVariables+1)*sizeof(short));
1948  }
1949  return r;
1950}
1951
1952int idMinDegW(ideal M,intvec *w)
1953{
1954  int d=-1;
1955  for(int i=0;i<IDELEMS(M);i++)
1956  {
1957    int d0=pMinDeg(M->m[i],w);
1958    if(-1<d0&&(d0<d||d==-1))
1959      d=d0;
1960  }
1961  return d;
1962}
1963
1964ideal idSeries(int n,ideal M,matrix U,intvec *w)
1965{
1966  for(int i=IDELEMS(M)-1;i>=0;i--)
1967  {
1968    if(U==NULL)
1969      M->m[i]=pSeries(n,M->m[i],NULL,w);
1970    else
1971    {
1972      M->m[i]=pSeries(n,M->m[i],MATELEM(U,i+1,i+1),w);
1973      MATELEM(U,i+1,i+1)=NULL;
1974    }
1975  }
1976  if(U!=NULL)
1977    idDelete((ideal*)&U);
1978  return M;
1979}
1980
1981matrix idDiff(matrix i, int k)
1982{
1983  int e=MATCOLS(i)*MATROWS(i);
1984  matrix r=mpNew(MATROWS(i),MATCOLS(i));
[360507]1985  r->rank=i->rank;
[35aab3]1986  int j;
1987  for(j=0; j<e; j++)
1988  {
1989    r->m[j]=pDiff(i->m[j],k);
1990  }
1991  return r;
1992}
1993
1994matrix idDiffOp(ideal I, ideal J,BOOLEAN multiply)
1995{
1996  matrix r=mpNew(IDELEMS(I),IDELEMS(J));
1997  int i,j;
1998  for(i=0; i<IDELEMS(I); i++)
1999  {
2000    for(j=0; j<IDELEMS(J); j++)
2001    {
2002      MATELEM(r,i+1,j+1)=pDiffOp(I->m[i],J->m[j],multiply);
2003    }
2004  }
2005  return r;
2006}
2007
[b3930d]2008int idElem(const ideal F)
[35aab3]2009{
[b3930d]2010  int i=0,j=IDELEMS(F)-1;
[35aab3]2011
[b3930d]2012  while(j>=0)
[35aab3]2013  {
[b3930d]2014    if ((F->m)[j]!=NULL) i++;
2015    j--;
[35aab3]2016  }
2017  return i;
2018}
2019
2020/*
2021*computes module-weights for liftings of homogeneous modules
2022*/
2023intvec * idMWLift(ideal mod,intvec * weights)
2024{
2025  if (idIs0(mod)) return new intvec(2);
2026  int i=IDELEMS(mod);
2027  while ((i>0) && (mod->m[i-1]==NULL)) i--;
2028  intvec *result = new intvec(i+1);
2029  while (i>0)
2030  {
[b130fb]2031    (*result)[i]=pFDeg(mod->m[i],currRing)+(*weights)[pGetComp(mod->m[i])];
[35aab3]2032  }
2033  return result;
2034}
2035
2036/*2
2037*sorts the kbase for idCoef* in a special way (lexicographically
2038*with x_max,...,x_1)
2039*/
2040ideal idCreateSpecialKbase(ideal kBase,intvec ** convert)
2041{
2042  int i;
2043  ideal result;
2044
2045  if (idIs0(kBase)) return NULL;
2046  result = idInit(IDELEMS(kBase),kBase->rank);
2047  *convert = idSort(kBase,FALSE);
2048  for (i=0;i<(*convert)->length();i++)
2049  {
2050    result->m[i] = pCopy(kBase->m[(**convert)[i]-1]);
2051  }
2052  return result;
2053}
2054
2055/*2
2056*returns the index of a given monom in the list of the special kbase
2057*/
2058int idIndexOfKBase(poly monom, ideal kbase)
2059{
2060  int j=IDELEMS(kbase);
2061
2062  while ((j>0) && (kbase->m[j-1]==NULL)) j--;
2063  if (j==0) return -1;
2064  int i=pVariables;
2065  while (i>0)
2066  {
2067    loop
2068    {
2069      if (pGetExp(monom,i)>pGetExp(kbase->m[j-1],i)) return -1;
2070      if (pGetExp(monom,i)==pGetExp(kbase->m[j-1],i)) break;
2071      j--;
2072      if (j==0) return -1;
2073    }
2074    if (i==1)
2075    {
2076      while(j>0)
2077      {
2078        if (pGetComp(monom)==pGetComp(kbase->m[j-1])) return j-1;
2079        if (pGetComp(monom)>pGetComp(kbase->m[j-1])) return -1;
2080        j--;
2081      }
2082    }
2083    i--;
2084  }
2085  return -1;
2086}
2087
2088/*2
2089*decomposes the monom in a part of coefficients described by the
2090*complement of how and a monom in variables occuring in how, the
2091*index of which in kbase is returned as integer pos (-1 if it don't
2092*exists)
2093*/
2094poly idDecompose(poly monom, poly how, ideal kbase, int * pos)
2095{
2096  int i;
2097  poly coeff=pOne(), base=pOne();
2098
2099  for (i=1;i<=pVariables;i++)
2100  {
2101    if (pGetExp(how,i)>0)
2102    {
2103      pSetExp(base,i,pGetExp(monom,i));
2104    }
2105    else
2106    {
2107      pSetExp(coeff,i,pGetExp(monom,i));
2108    }
2109  }
2110  pSetComp(base,pGetComp(monom));
2111  pSetm(base);
2112  pSetCoeff(coeff,nCopy(pGetCoeff(monom)));
2113  pSetm(coeff);
2114  *pos = idIndexOfKBase(base,kbase);
2115  if (*pos<0)
2116    pDelete(&coeff);
2117  pDelete(&base);
2118  return coeff;
2119}
2120
2121/*2
2122*returns a matrix A of coefficients with kbase*A=arg
2123*if all monomials in variables of how occur in kbase
2124*the other are deleted
2125*/
2126matrix idCoeffOfKBase(ideal arg, ideal kbase, poly how)
2127{
2128  matrix result;
2129  ideal tempKbase;
2130  poly p,q;
2131  intvec * convert;
2132  int i=IDELEMS(kbase),j=IDELEMS(arg),k,pos;
2133#if 0
2134  while ((i>0) && (kbase->m[i-1]==NULL)) i--;
2135  if (idIs0(arg))
2136    return mpNew(i,1);
2137  while ((j>0) && (arg->m[j-1]==NULL)) j--;
2138  result = mpNew(i,j);
2139#else
2140  result = mpNew(i, j);
2141  while ((j>0) && (arg->m[j-1]==NULL)) j--;
2142#endif
2143
2144  tempKbase = idCreateSpecialKbase(kbase,&convert);
2145  for (k=0;k<j;k++)
2146  {
2147    p = arg->m[k];
2148    while (p!=NULL)
2149    {
2150      q = idDecompose(p,how,tempKbase,&pos);
2151      if (pos>=0)
2152      {
2153        MATELEM(result,(*convert)[pos],k+1) =
2154            pAdd(MATELEM(result,(*convert)[pos],k+1),q);
2155      }
2156      else
2157        pDelete(&q);
2158      pIter(p);
2159    }
2160  }
2161  idDelete(&tempKbase);
2162  return result;
2163}
2164
2165/*3
[b8f199]2166* searches for the next unit in the components of the module arg and
2167* returns the first one;
[35aab3]2168*/
[b8f199]2169static int idReadOutPivot(ideal arg,int* comp)
[35aab3]2170{
[1d138c]2171  if (idIs0(arg)) return -1;
[8421b8]2172  int i=0,j, generator=-1;
2173  int rk_arg=arg->rank; //idRankFreeModule(arg);
2174  int * componentIsUsed =(int *)omAlloc((rk_arg+1)*sizeof(int));
[fc7902]2175  poly p;
[35aab3]2176
[8421b8]2177  while ((generator<0) && (i<IDELEMS(arg)))
[35aab3]2178  {
[8421b8]2179    memset(componentIsUsed,0,(rk_arg+1)*sizeof(int));
[35aab3]2180    p = arg->m[i];
2181    while (p!=NULL)
2182    {
2183      j = pGetComp(p);
[8421b8]2184      if (componentIsUsed[j]==0)
[35aab3]2185      {
[b8f199]2186#ifdef HAVE_RINGS
[4a7fc0]2187        if (pLmIsConstantComp(p) &&
2188            (!rField_is_Ring(currRing) || nIsUnit(pGetCoeff(p))))
[b8f199]2189        {
2190#else
[35aab3]2191        if (pLmIsConstantComp(p))
2192        {
[b8f199]2193#endif
[35aab3]2194          generator = i;
[8421b8]2195          componentIsUsed[j] = 1;
[35aab3]2196        }
2197        else
2198        {
[8421b8]2199          componentIsUsed[j] = -1;
[35aab3]2200        }
2201      }
[8421b8]2202      else if (componentIsUsed[j]>0)
[35aab3]2203      {
[8421b8]2204        (componentIsUsed[j])++;
[35aab3]2205      }
2206      pIter(p);
2207    }
2208    i++;
2209  }
2210  i = 0;
2211  *comp = -1;
2212  for (j=0;j<=rk_arg;j++)
2213  {
[8421b8]2214    if (componentIsUsed[j]>0)
[35aab3]2215    {
[8421b8]2216      if ((*comp==-1) || (componentIsUsed[j]<i))
[35aab3]2217      {
2218        *comp = j;
[8421b8]2219        i= componentIsUsed[j];
[35aab3]2220      }
2221    }
2222  }
[8421b8]2223  omFree(componentIsUsed);
[35aab3]2224  return generator;
2225}
2226
[955025]2227#if 0
[35aab3]2228static void idDeleteComp(ideal arg,int red_comp)
2229{
2230  int i,j;
2231  poly p;
2232
2233  for (i=IDELEMS(arg)-1;i>=0;i--)
2234  {
2235    p = arg->m[i];
2236    while (p!=NULL)
2237    {
2238      j = pGetComp(p);
2239      if (j>red_comp)
2240      {
2241        pSetComp(p,j-1);
2242        pSetm(p);
2243      }
2244      pIter(p);
2245    }
2246  }
2247  (arg->rank)--;
2248}
[955025]2249#endif
2250
2251static void idDeleteComps(ideal arg,int* red_comp,int del)
2252// red_comp is an array [0..args->rank]
2253{
2254  int i,j;
2255  poly p;
2256
2257  for (i=IDELEMS(arg)-1;i>=0;i--)
2258  {
2259    p = arg->m[i];
2260    while (p!=NULL)
2261    {
2262      j = pGetComp(p);
2263      if (red_comp[j]!=j)
2264      {
2265        pSetComp(p,red_comp[j]);
2266        pSetmComp(p);
2267      }
2268      pIter(p);
2269    }
2270  }
2271  (arg->rank) -= del;
2272}
[35aab3]2273
2274/*2
2275* returns the presentation of an isomorphic, minimally
2276* embedded  module (arg represents the quotient!)
2277*/
[cf108d]2278ideal idMinEmbedding(ideal arg,BOOLEAN inPlace, intvec **w)
[35aab3]2279{
2280  if (idIs0(arg)) return idInit(1,arg->rank);
[3504d7]2281  int i,next_gen,next_comp;
[35aab3]2282  ideal res=arg;
2283  if (!inPlace) res = idCopy(arg);
[8421b8]2284  res->rank=si_max(res->rank,idRankFreeModule(res));
[955025]2285  int *red_comp=(int*)omAlloc((res->rank+1)*sizeof(int));
2286  for (i=res->rank;i>=0;i--) red_comp[i]=i;
[8421b8]2287
[07b3e1]2288  int del=0;
[35aab3]2289  loop
2290  {
[b8f199]2291    next_gen = idReadOutPivot(res,&next_comp);
[35aab3]2292    if (next_gen<0) break;
[07b3e1]2293    del++;
[35aab3]2294    syGaussForOne(res,next_gen,next_comp,0,IDELEMS(res));
[955025]2295    for(i=next_comp+1;i<=arg->rank;i++) red_comp[i]--;
[cf108d]2296    if ((w !=NULL)&&(*w!=NULL))
2297    {
[07b3e1]2298      for(i=next_comp;i<(*w)->length();i++) (**w)[i-1]=(**w)[i];
[3504d7]2299    }
2300  }
[955025]2301
2302  idDeleteComps(res,red_comp,del);
2303  idSkipZeroes(res);
2304  omFree(red_comp);
2305
[07b3e1]2306  if ((w !=NULL)&&(*w!=NULL) &&(del>0))
[3504d7]2307  {
[07b3e1]2308    intvec *wtmp=new intvec((*w)->length()-del);
2309    for(i=0;i<res->rank;i++) (*wtmp)[i]=(**w)[i];
[3504d7]2310    delete *w;
2311    *w=wtmp;
[35aab3]2312  }
2313  return res;
2314}
2315
2316/*2
2317* transpose a module
2318*/
2319ideal idTransp(ideal a)
2320{
2321  int r = a->rank, c = IDELEMS(a);
2322  ideal b =  idInit(r,c);
2323
2324  for (int i=c; i>0; i--)
2325  {
2326    poly p=a->m[i-1];
2327    while(p!=NULL)
2328    {
2329      poly h=pHead(p);
2330      int co=pGetComp(h)-1;
2331      pSetComp(h,i);
2332      pSetmComp(h);
2333      b->m[co]=pAdd(b->m[co],h);
2334      pIter(p);
2335    }
2336  }
2337  return b;
2338}
2339
2340intvec * idQHomWeight(ideal id)
2341{
2342  poly head, tail;
2343  int k;
2344  int in=IDELEMS(id)-1, ready=0, all=0,
2345      coldim=pVariables, rowmax=2*coldim;
2346  if (in<0) return NULL;
2347  intvec *imat=new intvec(rowmax+1,coldim,0);
2348
2349  do
2350  {
2351    head = id->m[in--];
2352    if (head!=NULL)
2353    {
2354      tail = pNext(head);
2355      while (tail!=NULL)
2356      {
2357        all++;
2358        for (k=1;k<=coldim;k++)
2359          IMATELEM(*imat,all,k) = pGetExpDiff(head,tail,k);
2360        if (all==rowmax)
2361        {
2362          ivTriangIntern(imat, ready, all);
2363          if (ready==coldim)
2364          {
2365            delete imat;
2366            return NULL;
2367          }
2368        }
2369        pIter(tail);
2370      }
2371    }
2372  } while (in>=0);
2373  if (all>ready)
2374  {
2375    ivTriangIntern(imat, ready, all);
2376    if (ready==coldim)
2377    {
2378      delete imat;
2379      return NULL;
2380    }
2381  }
2382  intvec *result = ivSolveKern(imat, ready);
2383  delete imat;
2384  return result;
2385}
2386
2387BOOLEAN idIsZeroDim(ideal I)
2388{
2389  BOOLEAN *UsedAxis=(BOOLEAN *)omAlloc0(pVariables*sizeof(BOOLEAN));
2390  int i,n;
2391  poly po;
2392  BOOLEAN res=TRUE;
2393  for(i=IDELEMS(I)-1;i>=0;i--)
2394  {
2395    po=I->m[i];
2396    if ((po!=NULL) &&((n=pIsPurePower(po))!=0)) UsedAxis[n-1]=TRUE;
2397  }
2398  for(i=pVariables-1;i>=0;i--)
2399  {
2400    if(UsedAxis[i]==FALSE) {res=FALSE; break;} // not zero-dim.
2401  }
2402  omFreeSize(UsedAxis,pVariables*sizeof(BOOLEAN));
2403  return res;
2404}
2405
2406void idNormalize(ideal I)
2407{
2408  if (rField_has_simple_inverse()) return; /* Z/p, GF(p,n), R, long R/C */
2409  int i;
2410  for(i=IDELEMS(I)-1;i>=0;i--)
2411  {
[896561]2412    pNormalize(I->m[i]);
[35aab3]2413  }
2414}
[225d94]2415
[2ad10e9]2416// #include <kernel/clapsing.h>
[225d94]2417
[2f6fc61]2418#ifdef HAVE_FACTORY
[225d94]2419poly id_GCD(poly f, poly g, const ring r)
2420{
2421  ring save_r=currRing;
2422  rChangeCurrRing(r);
2423  ideal I=idInit(2,1); I->m[0]=f; I->m[1]=g;
2424  intvec *w = NULL;
2425  ideal S=idSyzygies(I,testHomog,&w);
2426  if (w!=NULL) delete w;
2427  poly gg=pTakeOutComp(&(S->m[0]),2);
2428  idDelete(&S);
2429  poly gcd_p=singclap_pdivide(f,gg);
2430  pDelete(&gg);
2431  rChangeCurrRing(save_r);
2432  return gcd_p;
2433}
[2f6fc61]2434#endif
[bba835]2435
2436/*2
2437* xx,q: arrays of length 0..rl-1
2438* xx[i]: SB mod q[i]
2439* assume: char=0
2440* assume: q[i]!=0
2441* destroys xx
2442*/
[2f6fc61]2443#ifdef HAVE_FACTORY
[bba835]2444ideal idChineseRemainder(ideal *xx, number *q, int rl)
2445{
[07969d]2446  int cnt=IDELEMS(xx[0])*xx[0]->nrows;
2447  ideal result=idInit(cnt,xx[0]->rank);
2448  result->nrows=xx[0]->nrows; // for lifting matrices
2449  result->ncols=xx[0]->ncols; // for lifting matrices
[bba835]2450  int i,j;
[cbc7e3]2451  poly r,h,hh,res_p;
[bba835]2452  number *x=(number *)omAlloc(rl*sizeof(number));
[07969d]2453  for(i=cnt-1;i>=0;i--)
[bba835]2454  {
2455    res_p=NULL;
2456    loop
2457    {
2458      r=NULL;
2459      for(j=rl-1;j>=0;j--)
2460      {
2461        h=xx[j]->m[i];
[4d8843]2462        if ((h!=NULL)
2463        &&((r==NULL)||(pLmCmp(r,h)==-1)))
2464          r=h;
[bba835]2465      }
2466      if (r==NULL) break;
[cbc7e3]2467      h=pHead(r);
[bba835]2468      for(j=rl-1;j>=0;j--)
2469      {
[cbc7e3]2470        hh=xx[j]->m[i];
2471        if ((hh!=NULL) && (pLmCmp(r,hh)==0))
2472        {
[38f763]2473          x[j]=pGetCoeff(hh);
[cbc7e3]2474          hh=pLmFreeAndNext(hh);
2475          xx[j]->m[i]=hh;
[bba835]2476        }
2477        else
[8391d8]2478          x[j]=nlInit(0, currRing);
[bba835]2479      }
2480      number n=nlChineseRemainder(x,q,rl);
2481      for(j=rl-1;j>=0;j--)
2482      {
[38f763]2483        x[j]=NULL; // nlInit(0...) takes no memory
[bba835]2484      }
[a8ef67]2485      if (nlIsZero(n)) pDelete(&h);
2486      else
2487      {
2488        pSetCoeff(h,n);
2489        //Print("new mon:");pWrite(h);
2490        res_p=pAdd(res_p,h);
2491      }
[bba835]2492    }
2493    result->m[i]=res_p;
2494  }
2495  omFree(x);
2496  for(i=rl-1;i>=0;i--) idDelete(&(xx[i]));
2497  omFree(xx);
2498  return result;
2499}
[2f6fc61]2500#endif
[3580b7]2501/* currently unsed:
[94eaef]2502ideal idChineseRemainder(ideal *xx, intvec *iv)
2503{
2504  int rl=iv->length();
2505  number *q=(number *)omAlloc(rl*sizeof(number));
2506  int i;
2507  for(i=0; i<rl; i++)
2508  {
2509    q[i]=nInit((*iv)[i]);
2510  }
2511  return idChineseRemainder(xx,q,rl);
[cbc7e3]2512}
[3580b7]2513*/
[3149a5]2514/*
2515 * lift ideal with coeffs over Z (mod N) to Q via Farey
2516 */
2517ideal idFarey(ideal x, number N)
2518{
[b86768]2519  int cnt=IDELEMS(x)*x->nrows;
2520  ideal result=idInit(cnt,x->rank);
2521  result->nrows=x->nrows; // for lifting matrices
2522  result->ncols=x->ncols; // for lifting matrices
2523
[3149a5]2524  int i;
[b86768]2525  for(i=cnt-1;i>=0;i--)
[3149a5]2526  {
2527    poly h=pCopy(x->m[i]);
2528    result->m[i]=h;
2529    while(h!=NULL)
2530    {
2531      number c=pGetCoeff(h);
2532      pSetCoeff0(h,nlFarey(c,N));
2533      nDelete(&c);
2534      pIter(h);
2535    }
[b86768]2536    while((result->m[i]!=NULL)&&(nIsZero(pGetCoeff(result->m[i]))))
2537    {
2538      pLmDelete(&(result->m[i]));
2539    }
2540    h=result->m[i];
2541    while((h!=NULL) && (pNext(h)!=NULL))
2542    {
2543      if(nIsZero(pGetCoeff(pNext(h))))
2544      {
2545        pLmDelete(&pNext(h));
2546      }
2547      else pIter(h);
2548    }
[3149a5]2549  }
2550  return result;
2551}
[90a60f]2552
2553/*2
2554* transpose a module
2555* NOTE: just a version of "ideal idTransp(ideal)" which works within specified ring.
2556*/
2557ideal id_Transp(ideal a, const ring rRing)
2558{
2559  int r = a->rank, c = IDELEMS(a);
2560  ideal b =  idInit(r,c);
2561
2562  for (int i=c; i>0; i--)
2563  {
2564    poly p=a->m[i-1];
2565    while(p!=NULL)
2566    {
2567      poly h=p_Head(p, rRing);
2568      int co=p_GetComp(h, rRing)-1;
2569      p_SetComp(h, i, rRing);
2570      p_Setm(h, rRing);
2571      b->m[co] = p_Add_q(b->m[co], h, rRing);
2572      pIter(p);
2573    }
2574  }
2575  return b;
2576}
2577
2578
2579
2580/*2
2581* The following is needed to compute the image of certain map used in
2582* the computation of cohomologies via BGG
2583* let M = { w_1, ..., w_k }, k = size(M) == ncols(M), n = nvars(currRing).
2584* assuming that nrows(M) <= m*n; the procedure computes:
2585* transpose(M) * transpose( var(1) I_m | ... | var(n) I_m ) :== transpose(module{f_1, ... f_k}),
2586* where f_i = \sum_{j=1}^{m} (w_i, v_j) gen(j),  (w_i, v_j) is a `scalar` multiplication.
2587* that is, if w_i = (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m) then
2588
2589  (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m)
2590*  var_1  ... var_1  |  var_2  ...  var_2  | ... |  var_n  ...  var(n)
2591*  gen_1  ... gen_m  |  gen_1  ...  gen_m  | ... |  gen_1  ...  gen_m
2592+ =>
2593  f_i =
2594
2595   a^1_1 * var(1) * gen(1) + ... + a^1_m * var(1) * gen(m) +
2596   a^2_1 * var(2) * gen(1) + ... + a^2_m * var(2) * gen(m) +
2597                             ...
2598   a^n_1 * var(n) * gen(1) + ... + a^n_m * var(n) * gen(m);
2599
2600   NOTE: for every f_i we run only ONCE along w_i saving partial sums into a temporary array of polys of size m
2601*/
[9c1b63]2602ideal id_TensorModuleMult(const int m, const ideal M, const ring rRing)
[90a60f]2603{
2604// #ifdef DEBU
2605//  WarnS("tensorModuleMult!!!!");
2606
2607  assume(m > 0);
2608  assume(M != NULL);
2609
2610  const int n = rRing->N;
2611
2612  assume(M->rank <= m * n);
2613
2614  const int k = IDELEMS(M);
2615
2616  ideal idTemp =  idInit(k,m); // = {f_1, ..., f_k }
2617
2618  for( int i = 0; i < k; i++ ) // for every w \in M
2619  {
2620    poly pTempSum = NULL;
2621
2622    poly w = M->m[i];
2623
2624    while(w != NULL) // for each term of w...
2625    {
2626      poly h = p_Head(w, rRing);
2627
2628      const int  gen = p_GetComp(h, rRing); // 1 ...
2629
2630      assume(gen > 0);
2631      assume(gen <= n*m);
2632
2633      // TODO: write a formula with %, / instead of while!
2634      /*
2635      int c = gen;
2636      int v = 1;
2637      while(c > m)
2638      {
2639        c -= m;
2640        v++;
2641      }
2642      */
2643
[592906]2644      int cc = gen % m;
[90a60f]2645      if( cc == 0) cc = m;
2646      int vv = 1 + (gen - cc) / m;
2647
2648//      assume( cc == c );
2649//      assume( vv == v );
2650
2651      //  1<= c <= m
2652      assume( cc > 0 );
2653      assume( cc <= m );
2654
2655      assume( vv > 0 );
2656      assume( vv <= n );
2657
2658      assume( (cc + (vv-1)*m) == gen );
2659
[9c1b63]2660      p_IncrExp(h, vv, rRing); // h *= var(j) && //      p_AddExp(h, vv, 1, rRing);
[592906]2661      p_SetComp(h, cc, rRing);
[90a60f]2662
2663      p_Setm(h, rRing);         // addjust degree after the previous steps!
2664
2665      pTempSum = p_Add_q(pTempSum, h, rRing); // it is slow since h will be usually put to the back of pTempSum!!!
2666
2667      pIter(w);
2668    }
2669
2670    idTemp->m[i] = pTempSum;
2671  }
2672
2673  // simplify idTemp???
2674
2675  ideal idResult = id_Transp(idTemp, rRing);
2676
2677  id_Delete(&idTemp, rRing);
2678
2679  return(idResult);
2680}
Note: See TracBrowser for help on using the repository browser.