source: git/Singular/ideals.cc @ 416465

spielwiese
Last change on this file since 416465 was 416465, checked in by Olaf Bachmann <obachman@…>, 24 years ago
* bug-fixes from work with Thomas git-svn-id: file:///usr/local/Singular/svn/trunk@3826 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 64.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: ideals.cc,v 1.73 1999-11-15 17:20:06 obachman Exp $ */
5/*
6* ABSTRACT - all basic methods to manipulate ideals
7*/
8
9/* includes */
10#include "mod2.h"
11#include "tok.h"
12#include "mmemory.h"
13#include "febase.h"
14#include "numbers.h"
15#include "polys.h"
16#include "ipid.h"
17#include "ring.h"
18#include "kstd1.h"
19#include "matpol.h"
20#include "weight.h"
21#include "intvec.h"
22#include "syz.h"
23#include "ideals.h"
24#include "lists.h"
25#include "prCopy.h"
26
27
28/* #define WITH_OLD_MINOR */
29#define pCopy_noCheck(p) pCopy(p)
30
31static poly * idpower;
32/*collects the monomials in makemonoms, must be allocated befor*/
33static int idpowerpoint;
34/*index of the actual monomial in idpower*/
35static poly * givenideal;
36/*the ideal from which a power is computed*/
37
38/*0 implementation*/
39
40/*2
41* initialise an ideal
42*/
43#ifdef PDEBUG
44ideal idDBInit(int idsize, int rank, char *f, int l)
45#else
46ideal idInit(int idsize, int rank)
47#endif
48{
49  /*- initialise an ideal -*/
50#if defined(MDEBUG) && defined(PDEBUG)
51  ideal hh = (ideal )mmDBAllocBlock(sizeof(sip_sideal),f,l);
52#else
53  ideal hh = (ideal )AllocSizeOf(sip_sideal);
54#endif
55  hh->nrows = 1;
56  hh->rank = rank;
57  IDELEMS(hh) = idsize;
58  if (idsize>0)
59  {
60#if defined(MDEBUG) && defined(PDEBUG)
61    hh->m = (poly *)mmDBAllocBlock0(idsize*sizeof(poly),f,l);
62#else
63    hh->m = (poly *)Alloc0(idsize*sizeof(poly));
64#endif
65  }
66  else
67    hh->m=NULL;
68  return hh;
69}
70
71#ifndef __OPTIMIZE__
72// this is mainly for outputting an ideal within the debugger
73void idPrint(ideal id)
74{
75  Print("Module of rank %d,real rank %d and %d generators.\n",
76         id->rank,idRankFreeModule(id),IDELEMS(id));
77  for (int i=0;i<IDELEMS(id);i++)
78  {
79    if (id->m[i]!=NULL)
80    {
81      Print("generator %d: ",i);pWrite(id->m[i]);
82    }
83  }
84}
85#endif
86
87/*2
88* initialise the maximal ideal (at 0)
89*/
90ideal idMaxIdeal (void)
91{
92  int l;
93  ideal hh=NULL;
94
95  hh=idInit(pVariables,1);
96  for (l=0; l<pVariables; l++)
97  {
98    hh->m[l] = pOne();
99    pSetExp(hh->m[l],l+1,1);
100    pSetm(hh->m[l]);
101  }
102  return hh;
103}
104
105/*2
106* deletes an ideal/matrix
107*/
108#ifdef PDEBUG
109void idDBDelete (ideal* h, char *f, int l)
110#else
111void idDelete (ideal * h)
112#endif
113{
114  int j,elems;
115  if (*h == NULL)
116    return;
117  elems=j=(*h)->nrows*(*h)->ncols;
118  if (j>0)
119  {
120    do
121    {
122      #if defined(MDEBUG) && defined(PDEBUG)
123      pDBDelete(&((*h)->m[--j]),mm_specHeap, f,l);
124      #else
125      pDelete(&((*h)->m[--j]));
126      #endif
127    }
128    while (j>0);
129    #if defined(MDEBUG) && defined(PDEBUG)
130    mmDBFreeBlock((ADDRESS)((*h)->m),sizeof(poly)*elems,f,l);
131    #else
132    Free((ADDRESS)((*h)->m),sizeof(poly)*elems);
133    #endif
134  }
135  #if defined(MDEBUG) && defined(PDEBUG)
136  mmDBFreeBlock((ADDRESS)(*h),sizeof(sip_sideal),f,l);
137  #else
138  FreeSizeOf((ADDRESS)*h,sip_sideal);
139  #endif
140  *h=NULL;
141}
142
143/*2
144*gives an ideal the minimal possible size
145*/
146void idSkipZeroes (ideal ide)
147{
148  int k;
149  int j = -1;
150  BOOLEAN change=FALSE;
151  for (k=0; k<IDELEMS(ide); k++)
152  {
153    if (ide->m[k] != NULL)
154    {
155      j++;
156      if (change)
157      {
158        ide->m[j] = ide->m[k];
159      }
160    }
161    else
162    {
163      change=TRUE;
164    }
165  }
166  if (change)
167  {
168    if (j == -1)
169      j = 0;
170    else
171    {
172      for (k=j+1; k<IDELEMS(ide); k++)
173        ide->m[k] = NULL;
174    }
175    pEnlargeSet(&(ide->m),IDELEMS(ide),j+1-IDELEMS(ide));
176    IDELEMS(ide) = j+1;
177  }
178}
179
180/*2
181* ideal id = (id[i])
182* result is leadcoeff(id[i]) = 1
183*/
184void idNorm(ideal id)
185{
186  for (int i=0; i<IDELEMS(id); i++)
187  {
188    if (id->m[i] != NULL)
189    {
190      pNorm(id->m[i]);
191    }
192  }
193}
194
195/*2
196* ideal id = (id[i]), c any number
197* if id[i] = c*id[j] then id[j] is deleted for j > i
198*/
199void idDelMultiples(ideal id)
200{
201  int i, j, t;
202  int k = IDELEMS(id), l = k;
203  for (i=k-2; i>=0; i--)
204  {
205    if (id->m[i]!=NULL)
206    {
207      for (j=l-1; j>i; j--)
208      {
209        if ((id->m[j]!=NULL)
210        && (pComparePolys(id->m[i], id->m[j])))
211        {
212          pDelete(&id->m[j]);
213          l--;
214          for(t=j; t<l; t++)
215          {
216            id->m[t] = id->m[t+1];
217          }
218        }
219      }
220    }
221  }
222  if (l != k)
223  {
224    pEnlargeSet(&id->m, k, l-k);
225    IDELEMS(id) = l;
226  }
227}
228
229/*2
230* ideal id = (id[i])
231* if id[i] = id[j] then id[j] is deleted for j > i
232*/
233void idDelEquals(ideal id)
234{
235  int i, j, t;
236  int k = IDELEMS(id), l = k;
237  for (i=k-2; i>=0; i--)
238  {
239    for (j=l-1; j>i; j--)
240    {
241      if (pEqualPolys(id->m[i], id->m[j]))
242      {
243        pDelete(&id->m[j]);
244        l--;
245        for(t=j; t<l; t++)
246        {
247          id->m[t] = id->m[t+1];
248        }
249      }
250    }
251  }
252  if (l != k)
253  {
254    pEnlargeSet(&id->m, k, l-k);
255    IDELEMS(id) = l;
256  }
257}
258
259//
260// Delete id[j], if Lm(j) == Lm(i) and j > i
261//
262void idDelLmEquals(ideal id)
263{
264  int i, j, t;
265  int k = IDELEMS(id), l = k;
266  for (i=k-2; i>=0; i--)
267  {
268    for (j=l-1; j>i; j--)
269    {
270      if (pLmEqual(id->m[i], id->m[j]))
271      {
272        pDelete(&id->m[j]);
273        l--;
274        for(t=j; t<l; t++)
275        {
276          id->m[t] = id->m[t+1];
277        }
278      }
279    }
280  }
281  if (l != k)
282  {
283    pEnlargeSet(&id->m, k, l-k);
284    IDELEMS(id) = l;
285  }
286}
287
288void idDelDiv(ideal id)
289{
290  int i, j, t;
291  int k = IDELEMS(id), l = k;
292  for (i=k-2; i>=0; i--)
293  {
294    for (j=l-1; j>i; j--)
295    {
296
297      if (((id->m[j] != NULL) && pDivisibleBy(id->m[i], id->m[j])) ||
298          (id->m[i] == NULL && id->m[j] == NULL))
299      {
300        pDelete(&id->m[j]);
301        l--;
302        for(t=j; t<l; t++)
303        {
304          id->m[t] = id->m[t+1];
305        }
306      }
307    }
308  }
309  if (l != k)
310  {
311    pEnlargeSet(&id->m, k, l-k);
312    IDELEMS(id) = l;
313  }
314}
315
316
317/*2
318* copy an ideal
319*/
320#ifdef PDEBUG
321ideal idDBCopy(ideal h1,char *f,int l)
322#else
323ideal idCopy (ideal h1)
324#endif
325{
326  int i;
327  ideal h2;
328
329#ifdef PDEBUG
330  idDBTest(h1,f,l);
331#endif
332//#ifdef TEST
333  if (h1 == NULL)
334  {
335#ifdef PDEBUG
336    h2=idDBInit(1,1,f,l);
337#else
338    h2=idInit(1,1);
339#endif
340  }
341  else
342//#endif
343  {
344#ifdef PDEBUG
345    h2=idDBInit(IDELEMS(h1),h1->rank,f,l);
346#else
347    h2=idInit(IDELEMS(h1),h1->rank);
348#endif
349    for (i=IDELEMS(h1)-1; i>=0; i--)
350#if defined(PDEBUG) && defined(MDEBUG)
351      h2->m[i] = pDBCopy(h1->m[i],f,l);
352#else
353      h2->m[i] = pCopy(h1->m[i]);
354#endif
355  }
356  return h2;
357}
358
359#ifdef PDEBUG
360void idDBTest(ideal h1,char *f,int l)
361{
362  int i;
363
364  if (h1 != NULL)
365  {
366    #ifdef MDEBUG
367    mmDBTestBlock(h1,sizeof(*h1),f,l);
368    /* to be able to test matrices: */
369    mmDBTestBlock(h1->m,h1->ncols*h1->nrows*sizeof(poly),f,l);
370    #endif
371    /* to be able to test matrices: */
372    for (i=(h1->ncols*h1->nrows)-1; i>=0; i--)
373      pDBTest(h1->m[i],f,l);
374    int new_rk=idRankFreeModule(h1);
375    if(new_rk > h1->rank)
376    {
377      Print("wrong rank %d (should be %d) in %s:%d\n",
378      h1->rank, new_rk, f,l);
379      h1->rank=new_rk;
380    }
381  }
382}
383#endif
384
385/*3
386* for idSort: compare a and b revlex inclusive module comp.
387*/
388static int pComp_RevLex(poly a, poly b)
389{
390 int l=pVariables;
391 while ((l>0) && (pGetExp(a,l)==pGetExp(b,l))) l--;
392 if (l==0)
393 {
394   if (pGetComp(a)==pGetComp(b)) return 0;
395   if (pGetComp(a)>pGetComp(b)) return 1;
396 }
397 else if (pGetExp(a,l)>pGetExp(b,l))
398   return 1;
399 return -1;
400}
401
402/*2
403*sorts the ideal w.r.t. the actual ringordering
404*uses lex-ordering when nolex = FALSE
405*/
406intvec *idSort(ideal id,BOOLEAN nolex)
407{
408  poly p,q;
409  intvec * result = NewIntvec1(IDELEMS(id));
410  int i, j, actpos=0, newpos, l;
411  int diff, olddiff, lastcomp, newcomp;
412  BOOLEAN notFound;
413
414  pCompProc oldComp=pComp0;
415
416  if (!nolex) pComp0=pComp_RevLex;
417
418  for (i=0;i<IDELEMS(id);i++)
419  {
420    if (id->m[i]!=NULL)
421    {
422      notFound = TRUE;
423      newpos = actpos / 2;
424      diff = (actpos+1) / 2;
425      diff = (diff+1) / 2;
426      lastcomp = pComp0(id->m[i],id->m[(*result)[newpos]]);
427      if (lastcomp<0)
428      {
429        newpos -= diff;
430      }
431      else if (lastcomp>0)
432      {
433        newpos += diff;
434      }
435      else
436      {
437        notFound = FALSE;
438      }
439      while ((newpos>=0) && (newpos<actpos) && (notFound))
440      {
441        newcomp = pComp0(id->m[i],id->m[(*result)[newpos]]);
442        olddiff = diff;
443        if (diff>1)
444        {
445          diff = (diff+1) / 2;
446          if ((newcomp==1)
447          && (actpos-newpos>1)
448          && (diff>1)
449          && (newpos+diff>=actpos))
450          {
451            diff = actpos-newpos-1;
452          }
453          else if ((newcomp==-1)
454          && (diff>1)
455          && (newpos<diff))
456          {
457            diff = newpos;
458          }
459        }
460        if (newcomp<0)
461        {
462          if ((olddiff==1) && (lastcomp>0))
463            notFound = FALSE;
464          else
465            newpos -= diff;
466        }
467        else if (newcomp>0)
468        {
469          if ((olddiff==1) && (lastcomp<0))
470          {
471            notFound = FALSE;
472            newpos++;
473          }
474          else
475          {
476            newpos += diff;
477          }
478        }
479        else
480        {
481          notFound = FALSE;
482        }
483        lastcomp = newcomp;
484        if (diff==0) notFound=FALSE; /*hs*/
485      }
486      if (newpos<0) newpos = 0;
487      if (newpos>actpos) newpos = actpos;
488      while ((newpos<actpos) && (pComp0(id->m[i],id->m[(*result)[newpos]])==0))
489        newpos++;
490      for (j=actpos;j>newpos;j--)
491      {
492        (*result)[j] = (*result)[j-1];
493      }
494      (*result)[newpos] = i;
495      actpos++;
496    }
497  }
498  for (j=0;j<actpos;j++) (*result)[j]++;
499  pComp0=oldComp;
500  return result;
501}
502
503/*2
504* concat the lists h1 and h2 without zeros
505*/
506ideal idSimpleAdd (ideal h1,ideal h2)
507{
508  int i,j,r,l;
509  ideal result;
510
511  if (h1==NULL) return idCopy(h2);
512  if (h2==NULL) return idCopy(h1);
513  j = IDELEMS(h1)-1;
514  while ((j >= 0) && (h1->m[j] == NULL)) j--;
515  i = IDELEMS(h2)-1;
516  while ((i >= 0) && (h2->m[i] == NULL)) i--;
517  r = max(h1->rank,h2->rank);
518  if (i+j==(-2))
519    return idInit(1,r);
520  else
521    result=idInit(i+j+2,r);
522  for (l=j; l>=0; l--)
523  {
524    result->m[l] = pCopy(h1->m[l]);
525  }
526  r = i+j+1;
527  for (l=i; l>=0; l--, r--)
528  {
529    result->m[r] = pCopy(h2->m[l]);
530  }
531  return result;
532}
533
534/*2
535* h1 + h2
536*/
537ideal idAdd (ideal h1,ideal h2)
538{
539  ideal result = idSimpleAdd(h1,h2);
540  ideal tmp = idCompactify(result);
541
542  idDelete(&result);
543  return tmp;
544}
545
546/*2
547* h1 * h2
548*/
549ideal  idMult (ideal h1,ideal  h2)
550{
551  int i,j,k;
552  ideal  hh;
553
554  j = IDELEMS(h1);
555  while ((j > 0) && (h1->m[j-1] == NULL)) j--;
556  i = IDELEMS(h2);
557  while ((i > 0) && (h2->m[i-1] == NULL)) i--;
558  j = j * i;
559  if (j == 0)
560    hh = idInit(1,1);
561  else
562    hh=idInit(j,1);
563  if (h1->rank<h2->rank)
564    hh->rank = h2->rank;
565  else
566    hh->rank = h1->rank;
567  if (j==0) return hh;
568  k = 0;
569  for (i=0; i<IDELEMS(h1); i++)
570  {
571    if (h1->m[i] != NULL)
572    {
573      for (j=0; j<IDELEMS(h2); j++)
574      {
575        if (h2->m[j] != NULL)
576        {
577          hh->m[k] = pMult(pCopy(h1->m[i]),pCopy(h2->m[j]));
578          k++;
579        }
580      }
581    }
582  }
583  {
584    ideal tmp = idCompactify(hh);
585    idDelete(&hh);
586    return tmp;
587    //return hh;
588  }
589}
590
591/*2
592*returns true if h is the zero ideal
593*/
594BOOLEAN idIs0 (ideal h)
595{
596  int i;
597
598  if (h == NULL) return TRUE;
599  i = IDELEMS(h);
600  while ((i > 0) && (h->m[i-1] == NULL))
601  {
602    i--;
603  }
604  if (i == 0)
605    return TRUE;
606  else
607    return FALSE;
608}
609
610/*2
611* return the maximal component number found in any polynomial in s
612*/
613int idRankFreeModule (ideal s)
614{
615  if (s!=NULL)
616  {
617    int  j=0;
618    int  l=IDELEMS(s);
619    poly *p=s->m;
620    int  k;
621
622    for (; l != 0; l--)
623    {
624      if (*p!=NULL)
625      {
626        k = pMaxComp(*p);
627        if (k>j) j = k;
628      }
629      p++;
630    }
631    return j;
632  }
633  return -1;
634}
635
636/*2
637*returns true if id is homogenous with respect to the aktual weights
638*/
639BOOLEAN idHomIdeal (ideal id, ideal Q)
640{
641  int i;
642  BOOLEAN     b;
643  if ((id == NULL) || (IDELEMS(id) == 0)) return TRUE;
644  i = 0;
645  b = TRUE;
646  while ((i < IDELEMS(id)) && b)
647  {
648    b = pIsHomogeneous(id->m[i]);
649    i++;
650  }
651  if ((b) && (Q!=NULL) && (IDELEMS(Q)>0))
652  {
653    i=0;
654    while ((i < IDELEMS(Q)) && b)
655    {
656      b = pIsHomogeneous(Q->m[i]);
657      i++;
658    }
659  }
660  return b;
661}
662
663/*2
664*returns a minimized set of generators of h1
665*/
666ideal idMinBase (ideal h1)
667{
668  ideal h2, h3,h4,e;
669  int j,k;
670  int i,l,ll;
671  intvec * wth;
672  BOOLEAN homog;
673
674  homog = idHomModule(h1,currQuotient,&wth);
675  if ((currRing->OrdSgn == 1) && (!homog))
676  {
677    Warn("minbase applies only to the local or homogeneous case");
678    e=idCopy(h1);
679    return e;
680  }
681  if ((currRing->OrdSgn == 1) && (homog))
682  {
683    lists re=min_std(h1,currQuotient,(tHomog)homog,&wth,NULL,0,3);
684    h2 = (ideal)re->m[1].data;
685    re->m[1].data = NULL;
686    re->m[1].rtyp = NONE;
687    re->Clean();
688    return h2;
689  }
690  e=idInit(1,h1->rank);
691  if (idIs0(h1))
692  {
693    return e;
694  }
695  pEnlargeSet(&(e->m),IDELEMS(e),15);
696  IDELEMS(e) = 16;
697  h2 = kStd(h1,currQuotient,isNotHomog,NULL);
698  h3 = idMaxIdeal();
699  h4=idMult(h2,h3);
700  idDelete(&h3);
701  h3=kStd(h4,currQuotient,isNotHomog,NULL);
702  k = IDELEMS(h3);
703  while ((k > 0) && (h3->m[k-1] == NULL)) k--;
704  j = -1;
705  l = IDELEMS(h2);
706  while ((l > 0) && (h2->m[l-1] == NULL)) l--;
707  for (i=l-1; i>=0; i--)
708  {
709    if (h2->m[i] != NULL)
710    {
711      ll = 0;
712      while ((ll < k) && ((h3->m[ll] == NULL)
713      || !pDivisibleBy(h3->m[ll],h2->m[i])))
714        ll++;
715      if (ll >= k)
716      {
717        j++;
718        if (j > IDELEMS(e)-1)
719        {
720          pEnlargeSet(&(e->m),IDELEMS(e),16);
721          IDELEMS(e) += 16;
722        }
723        e->m[j] = pCopy(h2->m[i]);
724      }
725    }
726  }
727  idDelete(&h2);
728  idDelete(&h3);
729  idDelete(&h4);
730  if (currQuotient!=NULL)
731  {
732    h3=idInit(1,e->rank);
733    h2=kNF(h3,currQuotient,e);
734    idDelete(&h3);
735    idDelete(&e);
736    e=h2;
737  }
738  idSkipZeroes(e);
739  return e;
740}
741
742/*2
743*the minimal index of used variables - 1
744*/
745int pLowVar (poly p)
746{
747  int k,l,lex;
748
749  if (p == NULL) return -1;
750
751  k = 32000;/*a very large dummy value*/
752  while (p != NULL)
753  {
754    l = 1;
755    lex = pGetExp(p,l);
756    while ((l < pVariables) && (lex == 0))
757    {
758      l++;
759      lex = pGetExp(p,l);
760    }
761    l--;
762    if (l < k) k = l;
763    pIter(p);
764  }
765  return k;
766}
767
768/*3
769*multiplies p with t (!cas) or  (t-1)
770*the index of t is:1, so we have to shift all variables
771*p is NOT in the actual ring, it has no t
772*/
773static poly pMultWithT (poly p,BOOLEAN cas)
774{
775  /*qp is the working pointer in p*/
776  /*result is the result, qresult is the working pointer*/
777  /*pp is p in the actual ring(shifted), qpp the working pointer*/
778  poly result,qp,pp;
779  poly qresult=NULL;
780  poly qpp=NULL;
781  int  i,j,lex;
782  number n;
783
784  pp = NULL;
785  result = NULL;
786  qp = p;
787  while (qp != NULL)
788  {
789    i = 0;
790    if (result == NULL)
791    {/*first monomial*/
792      result = pInit();
793      qresult = result;
794    }
795    else
796    {
797      qresult->next = pInit();
798      pIter(qresult);
799    }
800    for (j=pVariables-1; j>0; j--)
801    {
802      lex = pGetExp(qp,j);
803      pSetExp(qresult,j+1,lex);/*copy all variables*/
804    }
805    lex = pGetComp(qp);
806    pSetComp(qresult,lex);
807    n=nCopy(pGetCoeff(qp));
808    pSetCoeff0(qresult,n);
809    qresult->next = NULL;
810    pSetm(qresult);
811    /*qresult is now qp brought into the actual ring*/
812    if (cas)
813    { /*case: mult with t-1*/
814      pSetExp(qresult,1,0);
815      pSetm(qresult);
816      if (pp == NULL)
817      { /*first monomial*/
818        pp = pCopy(qresult);
819        qpp = pp;
820      }
821      else
822      {
823        qpp->next = pCopy(qresult);
824        pIter(qpp);
825      }
826      pGetCoeff(qpp)=nNeg(pGetCoeff(qpp));
827      /*now qpp contains -1*qp*/
828    }
829    pSetExp(qresult,1,1);/*this is mult. by t*/
830    pSetm(qresult);
831    pIter(qp);
832  }
833  /*
834  *now p is processed:
835  *result contains t*p
836  * if cas: pp contains -1*p (in the new ring)
837  */
838  if (cas)  qresult->next = pp;
839  /*  else      qresult->next = NULL;*/
840  return result;
841}
842
843/*3
844*deletes the place of t in p (t: variable with index 1)
845*p is NOT in the actual ring: it has pVariables+1 variables
846*/
847static poly pDivByT (poly * p,int size)
848{
849
850  poly result=NULL,
851       resultp=NULL , /** working pointer in result **/
852       pp;
853  int  i,j;
854
855  while (*p != NULL)
856  {
857    i = 0;
858    if (result == NULL)
859    {/*the first monomial*/
860      result = pInit();
861      resultp = result;
862      resultp->next = NULL;
863    }
864    else
865    {
866      resultp->next = pInit();
867      pIter(resultp);
868      resultp->next = NULL;
869    }
870    for (j=1; j<=pVariables; j++)
871    {
872      pSetExp(resultp,j,pGetExp(*p,j+1));
873    }
874    pSetComp(resultp,pGetComp(*p));
875    pSetCoeff0(resultp,pGetCoeff(*p));
876    pSetm(resultp);
877    pp = (*p)->next;
878    Free((ADDRESS)*p,size);
879    *p = pp;
880  }
881  return result;
882}
883
884/*2
885*dehomogenized the generators of the ideal id1 with the leading
886*monomial of p replaced by n
887*/
888ideal idDehomogen (ideal id1,poly p,number n)
889{
890  int i;
891  ideal result;
892
893  if (idIs0(id1))
894  {
895    return idInit(1,id1->rank);
896  }
897  result=idInit(IDELEMS(id1),id1->rank);
898  for (i=0; i<IDELEMS(id1); i++)
899  {
900    result->m[i] = pDehomogen(id1->m[i],p,n);
901  }
902  return result;
903}
904
905/*2
906* verschiebt die Indizees der Modulerzeugenden um i
907*/
908void pShift (poly * p,int i)
909{
910  poly qp1 = *p,qp2 = *p;/*working pointers*/
911  int     j = pMaxComp(*p),k = pMinComp(*p);
912
913  if (j+i < 0) return ;
914  while (qp1 != NULL)
915  {
916    if ((pGetComp(qp1)+i > 0) || ((j == -i) && (j == k)))
917    {
918      pSetComp(qp1,pGetComp(qp1)+i);
919      pSetmComp(qp1);
920      qp2 = qp1;
921      pIter(qp1);
922    }
923    else
924    {
925      if (qp2 == *p)
926      {
927        pIter(*p);
928        qp2->next = NULL;
929        pDelete(&qp2);
930        qp2 = *p;
931        qp1 = *p;
932      }
933      else
934      {
935        qp2->next = qp1->next;
936        qp1->next = NULL;
937        pDelete(&qp1);
938        qp1 = qp2->next;
939      }
940    }
941  }
942}
943
944/*2
945*initialized a field with r numbers between beg and end for the
946*procedure idNextChoise
947*/
948void idInitChoise (int r,int beg,int end,BOOLEAN  *endch,int * choise)
949{
950  /*returns the first choise of r numbers between beg and end*/
951  int i;
952  for (i=0; i<r; i++)
953  {
954    choise[i] = 0;
955  }
956  if (r <= end-beg+1)
957    for (i=0; i<r; i++)
958    {
959      choise[i] = beg+i;
960    }
961  if (r > end-beg+1)
962    *endch = TRUE;
963  else
964    *endch = FALSE;
965}
966
967/*2
968*returns the next choise of r numbers between beg and end
969*/
970void idGetNextChoise (int r,int end,BOOLEAN  *endch,int * choise)
971{
972  int i = r-1,j;
973  while ((i >= 0) && (choise[i] == end))
974  {
975    i--;
976    end--;
977  }
978  if (i == -1)
979    *endch = TRUE;
980  else
981  {
982    choise[i]++;
983    for (j=i+1; j<r; j++)
984    {
985      choise[j] = choise[i]+j-i;
986    }
987    *endch = FALSE;
988  }
989}
990
991/*2
992*takes the field choise of d numbers between beg and end, cancels the t-th
993*entree and searches for the ordinal number of that d-1 dimensional field
994* w.r.t. the algorithm of construction
995*/
996int idGetNumberOfChoise(int t, int d, int begin, int end, int * choise)
997{
998  int * localchoise,i,result=0;
999  BOOLEAN b=FALSE;
1000
1001  if (d<=1) return 1;
1002  localchoise=(int*)Alloc((d-1)*sizeof(int));
1003  idInitChoise(d-1,begin,end,&b,localchoise);
1004  while (!b)
1005  {
1006    result++;
1007    i = 0;
1008    while ((i<t) && (localchoise[i]==choise[i])) i++;
1009    if (i>=t)
1010    {
1011      i = t+1;
1012      while ((i<d) && (localchoise[i-1]==choise[i])) i++;
1013      if (i>=d)
1014      {
1015        Free((ADDRESS)localchoise,(d-1)*sizeof(int));
1016        return result;
1017      }
1018    }
1019    idGetNextChoise(d-1,end,&b,localchoise);
1020  }
1021  Free((ADDRESS)localchoise,(d-1)*sizeof(int));
1022  return 0;
1023}
1024
1025/*2
1026*computes the binomial coefficient
1027*/
1028int binom (int n,int r)
1029{
1030  int i,result;
1031
1032  if (r==0) return 1;
1033  if (n-r<r) return binom(n,n-r);
1034  result = n-r+1;
1035  for (i=2;i<=r;i++)
1036  {
1037    result *= n-r+i;
1038    result /= i;
1039  }
1040  return result;
1041}
1042
1043/*2
1044*the free module of rank i
1045*/
1046ideal idFreeModule (int i)
1047{
1048  int j;
1049  ideal h;
1050
1051  h=idInit(i,i);
1052  for (j=0; j<i; j++)
1053  {
1054    h->m[j] = pOne();
1055    pSetComp(h->m[j],j+1);
1056    pSetmComp(h->m[j]);
1057  }
1058  return h;
1059}
1060
1061/*2
1062* h3 := h1 intersect h2
1063*/
1064ideal idSect (ideal h1,ideal h2)
1065{
1066  ideal first=h2,second=h1,temp,temp1,result;
1067  int i,j,k,flength,slength,length,rank=min(h1->rank,h2->rank);
1068  intvec *w;
1069  poly p,q;
1070
1071  if ((idIs0(h1)) && (idIs0(h2)))  return idInit(1,rank);
1072  if (IDELEMS(h1)<IDELEMS(h2))
1073  {
1074    first = h1;
1075    second = h2;
1076  }
1077  flength = idRankFreeModule(first);
1078  slength = idRankFreeModule(second);
1079  length  = max(flength,slength);
1080  if (length==0)
1081  {
1082    length = 1;
1083  }
1084  j = IDELEMS(first);
1085  temp = idInit(j /*IDELEMS(first)*/,length+j);
1086
1087  ring orig_ring=currRing;
1088  ring syz_ring=rCurrRingAssureSyzComp();
1089  rSetSyzComp(length);
1090
1091  while ((j>0) && (first->m[j-1]==NULL)) j--;
1092  k = 0;
1093  for (i=0;i<j;i++)
1094  {
1095    if (first->m[i]!=NULL)
1096    {
1097      if (syz_ring==orig_ring)
1098        temp->m[k] = pCopy(first->m[i]);
1099      else
1100        temp->m[k] = prCopyR(first->m[i], orig_ring);
1101      q = pOne();
1102      pSetComp(q,i+1+length);
1103      pSetmComp(q);
1104      if (flength==0) pShift(&(temp->m[k]),1);
1105      p = temp->m[k];
1106      while (pNext(p)!=NULL) pIter(p);
1107      pNext(p) = q;
1108      k++;
1109    }
1110  }
1111  pEnlargeSet(&(temp->m),IDELEMS(temp),j+IDELEMS(second)-IDELEMS(temp));
1112  IDELEMS(temp) = j+IDELEMS(second);
1113  for (i=0;i<IDELEMS(second);i++)
1114  {
1115    if (second->m[i]!=NULL)
1116    {
1117      if (syz_ring==orig_ring)
1118        temp->m[k] = pCopy(second->m[i]);
1119      else
1120        temp->m[k] = prCopyR(second->m[i], orig_ring);
1121      if (slength==0) pShift(&(temp->m[k]),1);
1122      k++;
1123    }
1124  }
1125  temp1 = kStd(temp,currQuotient,testHomog,&w,NULL,length);
1126  if (w!=NULL) delete w;
1127  idDelete(&temp);
1128
1129  if(syz_ring!=orig_ring)
1130    rChangeCurrRing(orig_ring,TRUE);
1131
1132  result = idInit(IDELEMS(temp1),rank);
1133  j = 0;
1134  for (i=0;i<IDELEMS(temp1);i++)
1135  {
1136    if ((temp1->m[i]!=NULL)
1137    && (pRingGetComp(syz_ring,temp1->m[i])>length))
1138    {
1139      if(syz_ring==orig_ring)
1140        p = pCopy(temp1->m[i]);
1141      else
1142        p = prCopyR(temp1->m[i], syz_ring);
1143      while (p!=NULL)
1144      {
1145        q = pNext(p);
1146        pNext(p) = NULL;
1147        k = pGetComp(p)-1-length;
1148        pSetComp(p,0);
1149        pSetmComp(p);
1150        result->m[j] = pAdd(result->m[j],pMult(pCopy(first->m[k]),p));
1151        p = q;
1152      }
1153      j++;
1154    }
1155  }
1156  if(syz_ring!=orig_ring)
1157  {
1158    rChangeCurrRing(syz_ring,FALSE);
1159    idDelete(&temp1);
1160    rChangeCurrRing(orig_ring,TRUE);
1161    rKill(syz_ring);
1162  }
1163  else
1164  {
1165    idDelete(&temp1);
1166  }
1167
1168  idSkipZeroes(result);
1169  return result;
1170}
1171
1172/*2
1173* ideal/module intersection for a list of objects
1174* given as 'resolvente'
1175*/
1176ideal idMultSect(resolvente arg, int length)
1177{
1178  int i,j=0,k=0,syzComp,l,maxrk=-1,realrki;
1179  ideal bigmat,tempstd,result;
1180  poly p;
1181  int isIdeal=0;
1182  intvec * w=NULL;
1183
1184  /* find 0-ideals and max rank -----------------------------------*/
1185  for (i=0;i<length;i++)
1186  {
1187    if (!idIs0(arg[i]))
1188    {
1189      realrki=idRankFreeModule(arg[i]);
1190      k++;
1191      j += IDELEMS(arg[i]);
1192      if (realrki>maxrk) maxrk = realrki;
1193    }
1194    else
1195    {
1196      if (arg[i]!=NULL)
1197      {
1198        return idInit(1,arg[i]->rank);
1199      }
1200    }
1201  }
1202  if (maxrk == 0)
1203  {
1204    isIdeal = 1;
1205    maxrk = 1;
1206  }
1207  /* init -----------------------------------------------------------*/
1208  j += maxrk;
1209  syzComp = k*maxrk;
1210
1211  ring orig_ring=currRing;
1212  ring syz_ring=rCurrRingAssureSyzComp();
1213  rSetSyzComp(syzComp);
1214
1215  bigmat = idInit(j,(k+1)*maxrk);
1216  /* create unit matrices ------------------------------------------*/
1217  for (i=0;i<maxrk;i++)
1218  {
1219    for (j=0;j<=k;j++)
1220    {
1221      p = pOne();
1222      pSetComp(p,i+1+j*maxrk);
1223      pSetmComp(p);
1224      bigmat->m[i] = pAdd(bigmat->m[i],p);
1225    }
1226  }
1227  /* enter given ideals ------------------------------------------*/
1228  i = maxrk;
1229  k = 0;
1230  for (j=0;j<length;j++)
1231  {
1232    if (arg[j]!=NULL)
1233    {
1234      for (l=0;l<IDELEMS(arg[j]);l++)
1235      {
1236        if (arg[j]->m[l]!=NULL)
1237        {
1238          if (syz_ring==orig_ring)
1239            bigmat->m[i] = pCopy(arg[j]->m[l]);
1240          else
1241            bigmat->m[i] = prCopyR(arg[j]->m[l], orig_ring);
1242          pShift(&(bigmat->m[i]),k*maxrk+isIdeal);
1243          i++;
1244        }
1245      }
1246      k++;
1247    }
1248  }
1249  /* std computation --------------------------------------------*/
1250  tempstd = kStd(bigmat,currQuotient,testHomog,&w,NULL,syzComp);
1251  if (w!=NULL) delete w;
1252  idDelete(&bigmat);
1253
1254  if(syz_ring!=orig_ring)
1255    rChangeCurrRing(orig_ring,TRUE);
1256
1257  /* interprete result ----------------------------------------*/
1258  result = idInit(8,maxrk);
1259  k = 0;
1260  for (j=0;j<IDELEMS(tempstd);j++)
1261  {
1262    if ((tempstd->m[j]!=NULL) && (pRingGetComp(syz_ring,tempstd->m[j])>syzComp))
1263    {
1264      if (k>=IDELEMS(result))
1265      {
1266        pEnlargeSet(&(result->m),IDELEMS(result),8);
1267        IDELEMS(result) += 8;
1268      }
1269      if (syz_ring==orig_ring)
1270        p = pCopy(tempstd->m[j]);
1271      else
1272        p = prCopyR(tempstd->m[j], syz_ring);
1273      pShift(&p,-syzComp-isIdeal);
1274      result->m[k] = p;
1275      k++;
1276    }
1277  }
1278  /* clean up ----------------------------------------------------*/
1279  if(syz_ring!=orig_ring)
1280    rChangeCurrRing(syz_ring,FALSE);
1281  idDelete(&tempstd);
1282  if(syz_ring!=orig_ring)
1283  {
1284    rChangeCurrRing(orig_ring,TRUE);
1285    rKill(syz_ring);
1286  }
1287  idSkipZeroes(result);
1288  return result;
1289}
1290
1291/*2
1292*computes syzygies of h1,
1293*if quot != NULL it computes in the quotient ring modulo "quot"
1294*works always in a ring with ringorder_s
1295*/
1296static ideal idPrepare (ideal  h1, tHomog h, int syzcomp, intvec **w)
1297{
1298  ideal   h2, h3;
1299  int     i;
1300  int     j,jj=0,k;
1301  poly    p,q;
1302
1303  if (idIs0(h1)) return NULL;
1304  k = idRankFreeModule(h1);
1305  h2=idCopy(h1);
1306  i = IDELEMS(h2)-1;
1307  if (k == 0)
1308  {
1309    for (j=0; j<=i; j++) pShift(&(h2->m[j]),1);
1310    k = 1;
1311  }
1312  if (syzcomp<k)
1313  {
1314    Warn("syzcomp too low, should be %d instead of %d",k,syzcomp);
1315    syzcomp = k;
1316    rSetSyzComp(k);
1317  }
1318  h2->rank = syzcomp+i+1;
1319  for (j=0; j<=i; j++)
1320  {
1321    p = h2->m[j];
1322    q = pOne();
1323    pSetComp(q,syzcomp+1+j);
1324    pSetmComp(q);
1325    if (p!=NULL)
1326    {
1327      while (pNext(p)) pIter(p);
1328      p->next = q;
1329    }
1330    else
1331      h2->m[j]=q;
1332  }
1333
1334#ifdef PDEBUG
1335  for(j=0;j<IDELEMS(h2);j++) pTest(h2->m[j]);
1336#endif
1337  h3=kStd(h2,currQuotient,h,w,NULL,syzcomp);
1338  idDelete(&h2);
1339  return h3;
1340}
1341
1342/*2
1343* compute the syzygies of h1 in R/quot,
1344* weights of components are in w
1345* if setRegularity, return the regularity in deg
1346* do not change h1,  w
1347*/
1348ideal idSyzygies (ideal  h1, tHomog h,intvec **w, BOOLEAN setSyzComp,
1349                  BOOLEAN setRegularity, int &deg)
1350{
1351  ideal s_h1;
1352  poly  p;
1353  int   i, j, k, length=0,reg;
1354  BOOLEAN isMonomial=TRUE;
1355
1356#ifdef PDEBUG
1357  int ii;
1358  for(ii=0;ii<IDELEMS(h1);ii++) pTest(h1->m[ii]);
1359#endif
1360  if (idIs0(h1))
1361    return idFreeModule(IDELEMS(h1));
1362  k=max(1,idRankFreeModule(h1));
1363
1364  assume(currRing != NULL);
1365  ring orig_ring=currRing;
1366  ring syz_ring=rCurrRingAssureSyzComp();
1367
1368  if (setSyzComp)
1369    rSetSyzComp(k);
1370
1371  if (orig_ring != syz_ring)
1372  {
1373    s_h1=idrCopyR_NoSort(h1,orig_ring);
1374  }
1375  else
1376  {
1377    s_h1 = h1;
1378  }
1379
1380  ideal s_h3=idPrepare(s_h1,h,k,w);
1381
1382  if (s_h3==NULL)
1383  {
1384    return idFreeModule(IDELEMS(h1));
1385  }
1386
1387  if (orig_ring != syz_ring)
1388  {
1389    idDelete(&s_h1);
1390    for (j=0; j<IDELEMS(s_h3); j++)
1391    {
1392      if (s_h3->m[j] != NULL) 
1393      {
1394        if (pMinComp(s_h3->m[j],syz_ring) > k)
1395          pShift(&s_h3->m[j], -k);
1396        else
1397          pDelete(&s_h3->m[j]);
1398      }
1399    }
1400    idSkipZeroes(s_h3);
1401    s_h3->rank -= k;
1402    rChangeCurrRing(orig_ring, TRUE);
1403    s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
1404    rKill(syz_ring);
1405    idTest(s_h3);
1406    return s_h3;
1407  }
1408 
1409  ideal e = idInit(IDELEMS(s_h3), s_h3->rank);
1410
1411  for (j=0; j<IDELEMS(s_h3); j++)
1412  {
1413    if (s_h3->m[j] != NULL)
1414    {
1415      if (pMinComp(s_h3->m[j],syz_ring) <= k)
1416      {
1417        e->m[j] = s_h3->m[j];
1418        isMonomial=isMonomial && (pNext(s_h3->m[j])==NULL);
1419        pDelete(&pNext(s_h3->m[j]));
1420        s_h3->m[j] = NULL;
1421      }
1422    }
1423  }
1424 
1425  idSkipZeroes(s_h3);
1426  idSkipZeroes(e);
1427 
1428  if ((!isMonomial)
1429  && (!TEST_OPT_NOTREGULARITY)
1430  && (setRegularity)
1431  && (h==isHomog))
1432  {
1433    ring dp_C_ring = rCurrRingAssure_dp_C();
1434    if (dp_C_ring != syz_ring)
1435      e = idrMoveR_NoSort(e, syz_ring);
1436    resolvente res = sySchreyerResolvente(e,-1,&length,TRUE, TRUE);
1437    intvec * dummy = syBetti(res,length,&reg, *w);
1438    deg = reg+2;
1439    delete dummy;
1440    for (j=0;j<length;j++)
1441    {
1442      if (res[j]!=NULL) idDelete(&(res[j]));
1443    }
1444    Free((ADDRESS)res,length*sizeof(ideal));
1445    idDelete(&e);
1446    if (dp_C_ring != syz_ring)
1447    {
1448      rChangeCurrRing(syz_ring, TRUE);
1449      rKill(dp_C_ring);
1450    }
1451  }
1452  else
1453  {
1454    idDelete(&e);
1455  }
1456  idTest(s_h3);
1457  return s_h3;
1458}
1459
1460/*
1461*computes a standard basis for h1 and stores the transformation matrix
1462* in ma
1463*/
1464ideal idLiftStd (ideal  h1, matrix* ma, tHomog h)
1465{
1466  int   i, j, k, t, inputIsIdeal=idRankFreeModule(h1);
1467  poly  p=NULL, q, qq;
1468  intvec *w=NULL;
1469
1470  idDelete((ideal*)ma);
1471  *ma=mpNew(1,0);
1472  if (idIs0(h1))
1473    return idInit(1,h1->rank);
1474  k=max(1,idRankFreeModule(h1));
1475
1476  ring orig_ring=currRing;
1477  ring syz_ring=rCurrRingAssureSyzComp();
1478  rSetSyzComp(k);
1479
1480  ideal s_h1=h1;
1481
1482  if (orig_ring != syz_ring)
1483    s_h1 = idrCopyR_NoSort(h1,orig_ring);
1484  else
1485    s_h1 = h1;
1486
1487  ideal s_h3=idPrepare(s_h1,h,k,&w);
1488  ideal s_h2 = idInit(IDELEMS(s_h3), s_h3->rank);
1489
1490  if (w!=NULL) delete w;
1491  i = 0;
1492  for (j=0; j<IDELEMS(s_h3); j++)
1493  {
1494    if ((s_h3->m[j] != NULL) && (pMinComp(s_h3->m[j],syz_ring) <= k))
1495    {
1496      i++;
1497      q = s_h3->m[j];
1498      while (pNext(q) != NULL)
1499      {
1500        if (pGetComp(pNext(q)) > k)
1501        {
1502          s_h2->m[j] = pNext(q);
1503          pNext(q) = NULL;
1504        }
1505        else
1506        {
1507          pIter(q);
1508        }
1509      }
1510      if (!inputIsIdeal) pShift(&(s_h3->m[j]), -1);
1511    }
1512    else
1513    {
1514      pDelete(&(s_h3->m[j]));
1515    }
1516  }
1517
1518  idSkipZeroes(s_h3);
1519  j = IDELEMS(s_h1);
1520
1521  if (syz_ring!=orig_ring)
1522  {
1523    idDelete(&s_h1);
1524    rChangeCurrRing(orig_ring,TRUE);
1525  }
1526 
1527  idDelete((ideal*)ma);
1528  *ma = mpNew(j,i);
1529
1530  i = 1;
1531  for (j=0; j<IDELEMS(s_h2); j++)
1532  {
1533    if (s_h2->m[j] != NULL)
1534    {
1535      q = prMoveR( s_h2->m[j], syz_ring);
1536      s_h2->m[j] = NULL;
1537
1538      while (q != NULL)
1539      {
1540        p = q;
1541        pIter(q);
1542        pNext(p) = NULL;
1543        t=pGetComp(p);
1544        pSetComp(p,0);
1545        pSetmComp(p);
1546        MATELEM(*ma,t-k,i) = pAdd(MATELEM(*ma,t-k,i),p);
1547      }
1548      i++;
1549    }
1550  }
1551  idDelete(&s_h2);
1552
1553  for (i=0; i<IDELEMS(s_h3); i++)
1554  {
1555    s_h3->m[i] = prMoveR_NoSort(s_h3->m[i], syz_ring);
1556  }
1557
1558  if (syz_ring!=orig_ring) rKill(syz_ring);
1559  return s_h3;
1560}
1561
1562/*2
1563*computes a representation of the generators of submod with respect to those
1564* of mod
1565*/
1566ideal idLiftNonStB (ideal  mod, ideal submod,BOOLEAN goodShape)
1567{
1568  int   lsmod =idRankFreeModule(submod), i, j, k;
1569
1570  if (idIs0(mod))
1571    return idInit(1,mod->rank);
1572
1573  k=idRankFreeModule(mod);
1574  if  ((k!=0) && (lsmod==0)) lsmod=1;
1575  k=max(k,1);
1576
1577  ring orig_ring=currRing;
1578  ring syz_ring=rCurrRingAssureSyzComp();
1579  rSetSyzComp(k);
1580
1581  ideal s_mod, s_temp;
1582  if (orig_ring != syz_ring)
1583  {
1584    s_mod = idrCopyR_NoSort(mod,orig_ring);
1585    s_temp = idrCopyR_NoSort(submod,orig_ring);
1586  }
1587  else
1588  {
1589    s_mod = mod;
1590    s_temp = idCopy(submod);
1591  }
1592
1593  ideal s_h3=idPrepare(s_mod,(tHomog)FALSE,k,NULL);
1594
1595  if (!goodShape)
1596  {
1597    for (j=0;j<IDELEMS(s_h3);j++)
1598    {
1599      if ((s_h3->m[j] != NULL) && (pMinComp(s_h3->m[j]) > k))
1600        pDelete(&(s_h3->m[j]));
1601    }
1602  }
1603  idSkipZeroes(s_h3);
1604  if (lsmod==0)
1605  {
1606    for (j=IDELEMS(s_temp);j>0;j--)
1607    {
1608      if (s_temp->m[j-1]!=NULL)
1609        pShift(&(s_temp->m[j-1]),1);
1610    }
1611  }
1612  ideal s_result = kNF(s_h3,currQuotient,s_temp,k);
1613  s_result->rank = s_h3->rank;
1614  idDelete(&s_h3);
1615  idDelete(&s_temp);
1616
1617  for (j=0;j<IDELEMS(s_result);j++)
1618  {
1619    if (s_result->m[j]!=NULL)
1620    {
1621      if (pGetComp(s_result->m[j])<=k)
1622      {
1623        WerrorS("2nd module lies not in the first");
1624        idDelete(&s_result);
1625        s_result=idInit(1,1);
1626        break;
1627      }
1628      else
1629      {
1630        pShift(&(s_result->m[j]),-k);
1631        pNeg(s_result->m[j]);
1632      }
1633    }
1634  }
1635
1636  if(syz_ring!=orig_ring)
1637  {
1638    idDelete(&s_mod);
1639    rChangeCurrRing(orig_ring,TRUE);
1640    s_result = idrMoveR_NoSort(s_result, syz_ring);
1641    rKill(syz_ring);
1642  }
1643  return s_result;
1644}
1645
1646/*2
1647*computes a representation of the generators of submod with respect to those
1648* of mod which is given as standardbasis,
1649* uses currQuotient as the quotient ideal (if not NULL)
1650*/
1651ideal  idLift (ideal  mod,ideal submod)
1652{
1653  int   j,k;
1654  poly  p,q;
1655  BOOLEAN reported=FALSE;
1656
1657  if (idIs0(mod)) return idInit(1,mod->rank);
1658
1659  k = idRankFreeModule(mod);
1660
1661
1662  ring orig_ring=currRing;
1663  ring syz_ring=rCurrRingAssureSyzComp();
1664  rSetSyzComp(max(k,1));
1665
1666  ideal s_result=idInit(IDELEMS(submod),submod->rank);
1667  ideal s_temp;
1668 
1669  if (syz_ring != orig_ring)
1670  {
1671    s_temp = idrCopyR_NoSort(mod,orig_ring);
1672  }
1673  else
1674  {
1675    s_temp = idCopy(mod);
1676  }
1677 
1678  if (k == 0)
1679  {
1680    for (j=0; j<IDELEMS(s_temp); j++)
1681    {
1682      if (s_temp->m[j]!=NULL) pSetCompP(s_temp->m[j],1);
1683    }
1684    k = 1;
1685  }
1686  for (j=0; j<IDELEMS(s_temp); j++)
1687  {
1688    if (s_temp->m[j]!=NULL)
1689    {
1690      p = s_temp->m[j];
1691      q = pOne();
1692      pGetCoeff(q)=nNeg(pGetCoeff(q));   //set q to -1
1693      pSetComp(q,k+1+j);
1694      pSetmComp(q);
1695      while (pNext(p)) pIter(p);
1696      pNext(p) = q;
1697    }
1698  }
1699
1700  for (j=0; j<IDELEMS(submod); j++)
1701  {
1702    if (submod->m[j]!=NULL)
1703    {
1704      if (syz_ring==orig_ring)
1705        p = pCopy(submod->m[j]);
1706      else
1707        p=prCopyR(submod->m[j], orig_ring);
1708      if (pGetComp(p)==0) pSetCompP(p,1);
1709      q = kNF(s_temp,currQuotient,p,k);
1710      pDelete(&p);
1711      if (q!=NULL)
1712      {
1713        if (pMinComp(q)<=k)
1714        {
1715          if (!reported)
1716          {
1717            WarnS("first module not a standardbasis\n"
1718            "// ** or second not a proper submodule");
1719            reported=TRUE;
1720          }
1721          pDelete(&q);
1722        }
1723        else
1724        {
1725          pShift(&q,-k);
1726          s_result->m[j] = q;
1727        }
1728      }
1729    }
1730  }
1731  idDelete(&s_temp);
1732
1733  if(syz_ring!=orig_ring)
1734  {
1735    rChangeCurrRing(orig_ring,TRUE);
1736    s_result = idrMoveR_NoSort(s_result, syz_ring);
1737    rKill(syz_ring);
1738  }
1739
1740  return s_result;
1741}
1742
1743
1744/*2
1745*computes the quotient of h1,h2
1746*/
1747static ideal idInitializeQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb,
1748                               BOOLEAN *addOnlyOne, int *kkmax)
1749{
1750  ideal temph1;
1751  poly     p,q = NULL;
1752  int i,l,ll,k,kkk,kmax;
1753  int j = 0;
1754  int k1 = idRankFreeModule(h1);
1755  int k2 = idRankFreeModule(h2);
1756  tHomog   hom=isNotHomog;
1757
1758  k=max(k1,k2);
1759  if (k==0)
1760    k = 1;
1761  if ((k2==0) && (k>1)) *addOnlyOne = FALSE;
1762
1763  intvec * weights;
1764  hom = (tHomog)idHomModule(h1,currQuotient,&weights);
1765  if (addOnlyOne && (!h1IsStb))
1766    temph1 = kStd(h1,currQuotient,hom,&weights,NULL);
1767  else
1768    temph1 = idCopy(h1);
1769  if (weights!=NULL) delete weights;
1770  idTest(temph1);
1771/*--- making a single vector from h2 ---------------------*/
1772  for (i=0; i<IDELEMS(h2); i++)
1773  {
1774    if (h2->m[i] != NULL)
1775    {
1776      p = pCopy(h2->m[i]);
1777      if (k2 == 0)
1778        pShift(&p,j*k+1);
1779      else
1780        pShift(&p,j*k);
1781      q = pAdd(q,p);
1782      j++;
1783    }
1784  }
1785  *kkmax = kmax = j*k+1;
1786/*--- adding a monomial for the result (syzygy) ----------*/
1787  p = q;
1788  while (pNext(p)!=NULL) pIter(p);
1789  pNext(p) = pOne();
1790  pIter(p);
1791  pSetComp(p,kmax);
1792  pSetmComp(p);
1793/*--- constructing the big matrix ------------------------*/
1794  ideal h4 = idInit(16,kmax);
1795  h4->m[0] = q;
1796  if (k2 == 0)
1797  {
1798    if (k > IDELEMS(h4))
1799    {
1800      pEnlargeSet(&(h4->m),IDELEMS(h4),k-IDELEMS(h4));
1801      IDELEMS(h4) = k;
1802    }
1803    for (i=1; i<k; i++)
1804    {
1805      p = pCopy_noCheck(h4->m[i-1]);
1806      pShift(&p,1);
1807      h4->m[i] = p;
1808    }
1809  }
1810
1811  kkk = IDELEMS(h4);
1812  i = IDELEMS(temph1);
1813  while ((i>0) && (temph1->m[i-1]==NULL)) i--;
1814  for (l=0; l<i; l++)
1815  {
1816    if(temph1->m[l]!=NULL)
1817    {
1818      for (ll=0; ll<j; ll++)
1819      {
1820        p = pCopy(temph1->m[l]);
1821        if (k1 == 0)
1822          pShift(&p,ll*k+1);
1823        else
1824          pShift(&p,ll*k);
1825        if (kkk >= IDELEMS(h4))
1826        {
1827          pEnlargeSet(&(h4->m),IDELEMS(h4),16);
1828          IDELEMS(h4) += 16;
1829        }
1830        h4->m[kkk] = p;
1831        kkk++;
1832      }
1833    }
1834  }
1835/*--- if h2 goes in as single vector - the h1-part is just SB ---*/
1836  if (*addOnlyOne)
1837  {
1838    p = h4->m[0];
1839    for (i=0;i<IDELEMS(h4)-1;i++)
1840    {
1841      h4->m[i] = h4->m[i+1];
1842    }
1843    h4->m[IDELEMS(h4)-1] = p;
1844    idSkipZeroes(h4);
1845    test |= Sy_bit(OPT_SB_1);
1846  }
1847  idDelete(&temph1);
1848  idTest(h4);
1849  return h4;
1850}
1851/*2
1852*computes the quotient of h1,h2
1853*/
1854ideal idQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb, BOOLEAN resultIsIdeal)
1855{
1856  // first check for special case h1:(0)
1857  if (idIs0(h2))
1858  {
1859    ideal res;
1860    if (resultIsIdeal)
1861    {
1862      res = idInit(1,1);
1863      res->m[0] = pOne();
1864    }
1865    else
1866      res = idFreeModule(h1->rank);
1867    return res;
1868  }
1869  BITSET old_test=test;
1870  poly     p,q = NULL;
1871  int i,l,ll,k,kkk,kmax;
1872  BOOLEAN  addOnlyOne=TRUE;
1873  tHomog   hom=isNotHomog;
1874  intvec * weights1;
1875
1876  ideal s_h4 = idInitializeQuot (h1,h2,h1IsStb,&addOnlyOne,&kmax);
1877  hom = (tHomog)idHomModule(s_h4,currQuotient,&weights1);
1878  ring orig_ring=currRing;
1879  ring syz_ring=rCurrRingAssureSyzComp();
1880  rSetSyzComp(kmax-1);
1881  if (orig_ring!=syz_ring)
1882    s_h4 = idrMoveR_NoSort(s_h4,syz_ring);
1883
1884  ideal s_h3;
1885  if (addOnlyOne)
1886  {
1887    s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,kmax-1,IDELEMS(s_h4)-1);
1888  }
1889  else
1890  {
1891    s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,kmax-1);
1892  }
1893  idTest(s_h3);
1894  if (weights1!=NULL) delete weights1;
1895  idDelete(&s_h4);
1896
1897
1898  for (i=0;i<IDELEMS(s_h3);i++)
1899  {
1900    if ((s_h3->m[i]!=NULL) && (pGetComp(s_h3->m[i])>=kmax))
1901    {
1902      if (resultIsIdeal)
1903        pShift(&s_h3->m[i],-kmax);
1904      else
1905        pShift(&s_h3->m[i],-kmax+1);
1906    }
1907    else
1908      pDelete(&s_h3->m[i]);
1909  }
1910  if (resultIsIdeal)
1911    s_h3->rank = 1;
1912  else
1913    s_h3->rank = h1->rank;
1914  if(syz_ring!=orig_ring)
1915  {
1916//    pDelete(&q);
1917    rChangeCurrRing(orig_ring,TRUE);
1918    s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
1919    rKill(syz_ring);
1920  }
1921  idSkipZeroes(s_h3);
1922  test = old_test;
1923  idTest(s_h3);
1924  return s_h3;
1925}
1926
1927/*2
1928*computes recursively all monomials of a certain degree
1929*in every step the actvar-th entry in the exponential
1930*vector is incremented and the other variables are
1931*computed by recursive calls of makemonoms
1932*if the last variable is reached, the difference to the
1933*degree is computed directly
1934*vars is the number variables
1935*actvar is the actual variable to handle
1936*deg is the degree of the monomials to compute
1937*monomdeg is the actual degree of the monomial in consideration
1938*/
1939static void makemonoms(int vars,int actvar,int deg,int monomdeg)
1940{
1941  poly p;
1942  int i=0;
1943
1944  if ((idpowerpoint == 0) && (actvar ==1))
1945  {
1946    idpower[idpowerpoint] = pOne();
1947    monomdeg = 0;
1948  }
1949  while (i<=deg)
1950  {
1951    if (deg == monomdeg)
1952    {
1953      pSetm(idpower[idpowerpoint]);
1954      idpowerpoint++;
1955      return;
1956    }
1957    if (actvar == vars)
1958    {
1959      pSetExp(idpower[idpowerpoint],actvar,deg-monomdeg);
1960      pSetm(idpower[idpowerpoint]);
1961      pTest(idpower[idpowerpoint]);
1962      idpowerpoint++;
1963      return;
1964    }
1965    else
1966    {
1967      p = pCopy(idpower[idpowerpoint]);
1968      makemonoms(vars,actvar+1,deg,monomdeg);
1969      idpower[idpowerpoint] = p;
1970    }
1971    monomdeg++;
1972    pSetExp(idpower[idpowerpoint],actvar,pGetExp(idpower[idpowerpoint],actvar)+1);
1973    pSetm(idpower[idpowerpoint]);
1974    pTest(idpower[idpowerpoint]);
1975    i++;
1976  }
1977}
1978
1979/*2
1980*returns the deg-th power of the maximal ideal of 0
1981*/
1982ideal idMaxIdeal(int deg)
1983{
1984  if (deg < 0)
1985  {
1986    WarnS("maxideal: power must be non-negative");
1987  }
1988  if (deg < 1)
1989  {
1990    ideal I=idInit(1,1);
1991    I->m[0]=pOne();
1992    return I;
1993  }
1994  if (deg == 1)
1995  {
1996    return idMaxIdeal();
1997  }
1998
1999  int vars = currRing->N;
2000  int i = binom(vars+deg-1,deg);
2001  ideal id=idInit(i,1);
2002  idpower = id->m;
2003  idpowerpoint = 0;
2004  makemonoms(vars,1,deg,0);
2005  idpower = NULL;
2006  idpowerpoint = 0;
2007  return id;
2008}
2009
2010/*2
2011*computes recursively all generators of a certain degree
2012*of the ideal "givenideal"
2013*elms is the number elements in the given ideal
2014*actelm is the actual element to handle
2015*deg is the degree of the power to compute
2016*gendeg is the actual degree of the generator in consideration
2017*/
2018static void makepotence(int elms,int actelm,int deg,int gendeg)
2019{
2020  poly p;
2021  int i=0;
2022
2023  if ((idpowerpoint == 0) && (actelm ==1))
2024  {
2025    idpower[idpowerpoint] = pOne();
2026    gendeg = 0;
2027  }
2028  while (i<=deg)
2029  {
2030    if (deg == gendeg)
2031    {
2032      idpowerpoint++;
2033      return;
2034    }
2035    if (actelm == elms)
2036    {
2037      p=pPower(pCopy(givenideal[actelm-1]),deg-gendeg);
2038      idpower[idpowerpoint]=pMult(idpower[idpowerpoint],p);
2039      idpowerpoint++;
2040      return;
2041    }
2042    else
2043    {
2044      p = pCopy(idpower[idpowerpoint]);
2045      makepotence(elms,actelm+1,deg,gendeg);
2046      idpower[idpowerpoint] = p;
2047    }
2048    gendeg++;
2049    idpower[idpowerpoint]=pMult(idpower[idpowerpoint],pCopy(givenideal[actelm-1]));
2050    i++;
2051  }
2052}
2053
2054/*2
2055*returns the deg-th power of the ideal gid
2056*/
2057//ideal idPower(ideal gid,int deg)
2058//{
2059//  int i;
2060//  ideal id;
2061//
2062//  if (deg < 1) deg = 1;
2063//  i = binom(IDELEMS(gid)+deg-1,deg);
2064//  id=idInit(i,1);
2065//  idpower = id->m;
2066//  givenideal = gid->m;
2067//  idpowerpoint = 0;
2068//  makepotence(IDELEMS(gid),1,deg,0);
2069//  idpower = NULL;
2070//  givenideal = NULL;
2071//  idpowerpoint = 0;
2072//  return id;
2073//}
2074static void idNextPotence(ideal given, ideal result,
2075  int begin, int end, int deg, int restdeg, poly ap)
2076{
2077  poly p;
2078  int i;
2079
2080  p = pPower(pCopy(given->m[begin]),restdeg);
2081  i = result->nrows;
2082  result->m[i] = pMult(pCopy(ap),p);
2083//PrintS(".");
2084  (result->nrows)++;
2085  if (result->nrows >= IDELEMS(result))
2086  {
2087    pEnlargeSet(&(result->m),IDELEMS(result),16);
2088    IDELEMS(result) += 16;
2089  }
2090  if (begin == end) return;
2091  for (i=restdeg-1;i>0;i--)
2092  {
2093    p = pPower(pCopy(given->m[begin]),i);
2094    p = pMult(pCopy(ap),p);
2095    idNextPotence(given, result, begin+1, end, deg, restdeg-i, p);
2096    pDelete(&p);
2097  }
2098  idNextPotence(given, result, begin+1, end, deg, restdeg, ap);
2099}
2100
2101ideal idPower(ideal given,int exp)
2102{
2103  ideal result,temp;
2104  poly p1;
2105  int i;
2106
2107  if (idIs0(given)) return idInit(1,1);
2108  temp = idCopy(given);
2109  idSkipZeroes(temp);
2110  i = binom(IDELEMS(temp)+exp-1,exp);
2111  result = idInit(i,1);
2112  result->nrows = 0;
2113//Print("ideal contains %d elements\n",i);
2114  p1=pOne();
2115  idNextPotence(temp,result,0,IDELEMS(temp)-1,exp,exp,p1);
2116  pDelete(&p1);
2117  idDelete(&temp);
2118  result->nrows = 1;
2119  idSkipZeroes(result);
2120  idDelEquals(result);
2121  return result;
2122}
2123
2124/*2
2125* eliminate delVar (product of vars) in h1
2126*/
2127ideal idElimination (ideal h1,poly delVar,intvec *hilb)
2128{
2129  int    i,j=0,k,l;
2130  ideal  h,hh, h3;
2131  int    *ord,*block0,*block1;
2132  int    ordersize=2;
2133  int    **wv;
2134  tHomog hom;
2135  intvec * w;
2136  sip_sring tmpR;
2137  ring origR = currRing;
2138
2139  if (delVar==NULL)
2140  {
2141    return idCopy(h1);
2142  }
2143  if (currQuotient!=NULL)
2144  {
2145    WerrorS("cannot eliminate in a qring");
2146    return idCopy(h1);
2147  }
2148  if (idIs0(h1)) return idInit(1,h1->rank);
2149  hom=(tHomog)idHomModule(h1,NULL,&w); //sets w to weight vector or NULL
2150  h3=idInit(16,h1->rank);
2151  for (k=0;; k++)
2152  {
2153    if (currRing->order[k]!=0) ordersize++;
2154    else break;
2155  }
2156  ord=(int*)Alloc0(ordersize*sizeof(int));
2157  block0=(int*)Alloc(ordersize*sizeof(int));
2158  block1=(int*)Alloc(ordersize*sizeof(int));
2159  for (k=0;; k++)
2160  {
2161    if (currRing->order[k]!=0)
2162    {
2163      block0[k+1] = currRing->block0[k];
2164      block1[k+1] = currRing->block1[k];
2165      ord[k+1] = currRing->order[k];
2166    }
2167    else
2168      break;
2169  }
2170  block0[0] = 1;
2171  block1[0] = pVariables;
2172  wv=(int**) Alloc0(ordersize*sizeof(int**));
2173  memcpy4(wv+1,currRing->wvhdl,(ordersize-1)*sizeof(int**));
2174  wv[0]=(int*)AllocL((pVariables+1)*sizeof(int));
2175  memset(wv[0],0,(pVariables+1)*sizeof(int));
2176  for (j=0;j<pVariables;j++)
2177    if (pGetExp(delVar,j+1)!=0) wv[0][j]=1;
2178  ord[0] = ringorder_a;
2179
2180  // fill in tmp ring to get back the data later on
2181  tmpR  = *origR;
2182  tmpR.order = ord;
2183  tmpR.block0 = block0;
2184  tmpR.block1 = block1;
2185  tmpR.wvhdl = wv;
2186  rComplete(&tmpR, 1);
2187
2188  // change into the new ring
2189  //pChangeRing(pVariables,currRing->OrdSgn,ord,block0,block1,wv);
2190  rChangeCurrRing(&tmpR, TRUE);
2191  currRing = &tmpR;
2192  h = idInit(IDELEMS(h1),1);
2193  // fetch data from the old ring
2194  for (k=0;k<IDELEMS(h1);k++) h->m[k] = prCopyR( h1->m[k], origR);
2195  // compute kStd
2196  hh = kStd(h,NULL,hom,&w,hilb);
2197  idDelete(&h);
2198
2199  // go back to the original ring
2200  rChangeCurrRing(origR,TRUE);
2201  i = IDELEMS(hh)-1;
2202  while ((i >= 0) && (hh->m[i] == NULL)) i--;
2203  j = -1;
2204  // fetch data from temp ring
2205  for (k=0; k<=i; k++)
2206  {
2207    l=pVariables;
2208    while ((l>0) && (pRingGetExp(&tmpR, hh->m[k],l)*pGetExp(delVar,l)==0)) l--;
2209    if (l==0)
2210    {
2211      j++;
2212      if (j >= IDELEMS(h3))
2213      {
2214        pEnlargeSet(&(h3->m),IDELEMS(h3),16);
2215        IDELEMS(h3) += 16;
2216      }
2217      h3->m[j] = prCopyR( hh->m[k], &tmpR);
2218    }
2219  }
2220  rChangeCurrRing(&tmpR, FALSE);
2221  idDelete(&hh);
2222  rChangeCurrRing(origR, TRUE);
2223  idSkipZeroes(h3);
2224  FreeL((ADDRESS)wv[0]);
2225  Free((ADDRESS)wv,ordersize*sizeof(int**));
2226  Free((ADDRESS)ord,ordersize*sizeof(int));
2227  Free((ADDRESS)block0,ordersize*sizeof(int));
2228  Free((ADDRESS)block1,ordersize*sizeof(int));
2229  rUnComplete(&tmpR);
2230  if (w!=NULL)
2231    delete w;
2232  return h3;
2233}
2234
2235//void idEnterSet (poly p,ideal r, int * next)
2236//{
2237//
2238//  if ((*next) == IDELEMS(r)-1)
2239//  {
2240//    pEnlargeSet(&(r->m),IDELEMS(r),16);
2241//    IDELEMS(r)+=16;
2242//  }
2243//  int at;
2244//  int i;
2245//  if (*next==0) at=0;
2246//  else
2247//  {
2248//    int an = 0;
2249//    int en= *next-1;
2250//    int c;
2251//    if (pComp0(r->m[(*next)-1],p)!= 1)
2252//      at=*next;
2253//    else
2254//    {
2255//      loop
2256//      {
2257//        if (an >= en-1)
2258//        {
2259//          if (pComp0(r->m[an],p) == 1)
2260//          {
2261//            at=an; break;
2262//          }
2263//          else
2264//          {
2265//            at=en; break;
2266//          }
2267//        }
2268//        i=(an+en) / 2;
2269//        if (pComp0(r->m[i],p) == 1) en=i;
2270//        else                       an=i;
2271//      }
2272//    }
2273//  }
2274//  if (pComp(r->m[at],p)==0)
2275//  {
2276//    pDelete(&p);
2277//  }
2278//  else
2279//  {
2280//    (*next)++;
2281//    for (i=(*next); i>=at+1; i--)
2282//    {
2283//      r->m[i] = r->m[i-1];
2284//    }
2285//    /*- save result -*/
2286//    r->m[at] = p;
2287//  }
2288//}
2289
2290#ifdef WITH_OLD_MINOR
2291/*2
2292* compute all ar-minors of the matrix a
2293*/
2294ideal idMinors(matrix a, int ar)
2295{
2296  int     i,j,k,size;
2297  int *rowchoise,*colchoise;
2298  BOOLEAN rowch,colch;
2299  ideal result;
2300  matrix tmp;
2301  poly p;
2302
2303  i = binom(a->rows(),ar);
2304  j = binom(a->cols(),ar);
2305
2306  rowchoise=(int *)Alloc(ar*sizeof(int));
2307  colchoise=(int *)Alloc(ar*sizeof(int));
2308  if ((i>512) || (j>512) || (i*j >512)) size=512;
2309  else size=i*j;
2310  result=idInit(size,1);
2311  tmp=mpNew(ar,ar);
2312  k = 0; /* the index in result*/
2313  idInitChoise(ar,1,a->rows(),&rowch,rowchoise);
2314  while (!rowch)
2315  {
2316    idInitChoise(ar,1,a->cols(),&colch,colchoise);
2317    while (!colch)
2318    {
2319      for (i=1; i<=ar; i++)
2320      {
2321        for (j=1; j<=ar; j++)
2322        {
2323          MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
2324        }
2325      }
2326      p = mpDetBareiss(tmp);
2327      if (p!=NULL)
2328      {
2329        if (k>=size)
2330        {
2331          pEnlargeSet(&result->m,size,32);
2332          size += 32;
2333        }
2334        result->m[k] = p;
2335        k++;
2336      }
2337      idGetNextChoise(ar,a->cols(),&colch,colchoise);
2338    }
2339    idGetNextChoise(ar,a->rows(),&rowch,rowchoise);
2340  }
2341  /*delete the matrix tmp*/
2342  for (i=1; i<=ar; i++)
2343  {
2344    for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
2345  }
2346  idDelete((ideal*)&tmp);
2347  if (k==0)
2348  {
2349    k=1;
2350    result->m[0]=NULL;
2351  }
2352  Free((ADDRESS)rowchoise,ar*sizeof(int));
2353  Free((ADDRESS)colchoise,ar*sizeof(int));
2354  pEnlargeSet(&result->m,size,k-size);
2355  IDELEMS(result) = k;
2356  return (result);
2357}
2358#else
2359/*2
2360* compute all ar-minors of the matrix a
2361* the caller of mpRecMin
2362*/
2363ideal idMinors(matrix a, int ar)
2364{
2365  ideal result;
2366  int elems=0;
2367
2368  if((ar<=0) || (ar>min(a->ncols,a->nrows)))
2369  {
2370    Werror("%d-th minor, matrix is %dx%d",ar,a->ncols,a->nrows);
2371    return NULL;
2372  }
2373  a = mpCopy(a);
2374  result=idInit(32,1);
2375  if(ar>1) mpRecMin(ar-1,result,elems,a,a->nrows,a->ncols,NULL);
2376  else mpMinorToResult(result,elems,a,a->nrows,a->ncols);
2377  idDelete((ideal *)&a);
2378  idSkipZeroes(result);
2379  idTest(result);
2380  return result;
2381}
2382#endif
2383
2384/*2
2385*returns TRUE if p is a unit element in the current ring
2386*/
2387BOOLEAN pIsUnit(poly p)
2388{
2389  int i;
2390
2391  if (p == NULL) return FALSE;
2392  i = 1;
2393  while (i<=pVariables && pGetExp(p,i) == 0) i++;
2394  if (i > pVariables && (pGetComp(p) == 0))
2395  {
2396    if (currRing->OrdSgn == 1 && pNext(p) !=NULL) return FALSE;
2397    return TRUE;
2398  }
2399  return FALSE;
2400}
2401
2402/*2
2403*skips all zeroes and double elements, searches also for units
2404*/
2405ideal idCompactify(ideal id)
2406{
2407  int i,j;
2408  BOOLEAN b=FALSE;
2409
2410  i = IDELEMS(id)-1;
2411  while ((! b) && (i>=0))
2412  {
2413    b=pIsUnit(id->m[i]);
2414    i--;
2415  }
2416  if (b)
2417  {
2418    ideal result=idInit(1,id->rank);
2419    result->m[0]=pOne();
2420    return result;
2421  }
2422  else
2423  {
2424    ideal result=idCopy(id);
2425    for (i=1;i<IDELEMS(result);i++)
2426    {
2427      if (result->m[i]!=NULL)
2428      {
2429        for (j=0;j<i;j++)
2430        {
2431          if ((result->m[j]!=NULL)
2432          && (pComparePolys(result->m[i],result->m[j])))
2433          {
2434            pDelete(&(result->m[j]));
2435          }
2436        }
2437      }
2438    }
2439    idSkipZeroes(result);
2440    return result;
2441  }
2442}
2443
2444/*2
2445*returns TRUE if id1 is a submodule of id2
2446*/
2447BOOLEAN idIsSubModule(ideal id1,ideal id2)
2448{
2449  int i;
2450  poly p;
2451
2452  if (idIs0(id1)) return TRUE;
2453  for (i=0;i<IDELEMS(id1);i++)
2454  {
2455    if (id1->m[i] != NULL)
2456    {
2457      p = kNF(id2,currQuotient,id1->m[i]);
2458      if (p != NULL)
2459      {
2460        pDelete(&p);
2461        return FALSE;
2462      }
2463    }
2464  }
2465  return TRUE;
2466}
2467
2468/*2
2469* returns the ideals of initial terms
2470*/
2471ideal idHead(ideal h)
2472{
2473  ideal m = idInit(IDELEMS(h),h->rank);
2474  int i;
2475
2476  for (i=IDELEMS(h)-1;i>=0; i--)
2477  {
2478    if (h->m[i]!=NULL) m->m[i]=pHead(h->m[i]);
2479  }
2480  return m;
2481}
2482
2483ideal idHomogen(ideal h, int varnum)
2484{
2485  ideal m = idInit(IDELEMS(h),h->rank);
2486  int i;
2487
2488  for (i=IDELEMS(h)-1;i>=0; i--)
2489  {
2490    m->m[i]=pHomogen(h->m[i],varnum);
2491  }
2492  return m;
2493}
2494
2495/*------------------type conversions----------------*/
2496ideal idVec2Ideal(poly vec)
2497{
2498   ideal result=idInit(1,1);
2499   FreeSizeOf((ADDRESS)result->m,poly);
2500   result->m=NULL; // remove later
2501   pVec2Polys(vec, &(result->m), &(IDELEMS(result)));
2502   return result;
2503}
2504
2505// converts mat to module, destroys mat
2506ideal idMatrix2Module(matrix mat)
2507{
2508  ideal result = idInit(MATCOLS(mat),MATROWS(mat));
2509  int i,j;
2510  poly h;
2511#ifdef DRING
2512  poly p;
2513#endif
2514
2515  for(j=0;j<MATCOLS(mat);j++) /* j is also index in result->m */
2516  {
2517    for (i=1;i<=MATROWS(mat);i++)
2518    {
2519      h = MATELEM(mat,i,j+1);
2520      if (h!=NULL)
2521      {
2522        MATELEM(mat,i,j+1)=NULL;
2523        pSetCompP(h,i);
2524#ifdef DRING
2525        pdSetDFlagP(h,0);
2526#endif
2527        result->m[j] = pAdd(result->m[j],h);
2528      }
2529    }
2530  }
2531  // obachman: need to clean this up
2532  idDelete((ideal*) &mat);
2533  return result;
2534}
2535
2536/*2
2537* converts a module into a matrix, destroyes the input
2538*/
2539matrix idModule2Matrix(ideal mod)
2540{
2541  matrix result = mpNew(mod->rank,IDELEMS(mod));
2542  int i,cp;
2543  poly p,h;
2544
2545  for(i=0;i<IDELEMS(mod);i++)
2546  {
2547    p=mod->m[i];
2548    mod->m[i]=NULL;
2549    while (p!=NULL)
2550    {
2551      h=p;
2552      pIter(p);
2553      pNext(h)=NULL;
2554//      cp = max(1,pGetComp(h));     // if used for ideals too
2555      cp = pGetComp(h);
2556      pSetComp(h,0);
2557      pSetmComp(h);
2558#ifdef TEST
2559      if (cp>mod->rank)
2560      {
2561        Print("## inv. rank %d -> %d\n",mod->rank,cp);
2562        int k,l,o=mod->rank;
2563        mod->rank=cp;
2564        matrix d=mpNew(mod->rank,IDELEMS(mod));
2565        for (l=1; l<=o; l++)
2566        {
2567          for (k=1; k<=IDELEMS(mod); k++)
2568          {
2569            MATELEM(d,l,k)=MATELEM(result,l,k);
2570            MATELEM(result,l,k)=NULL;
2571          }
2572        }
2573        idDelete((ideal *)&result);
2574        result=d;
2575      }
2576#endif
2577      MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
2578    }
2579  }
2580// obachman 10/99: added the following line, otherwise memory lack!
2581  idDelete(&mod);
2582  return result;
2583}
2584
2585matrix idModule2formatedMatrix(ideal mod,int rows, int cols)
2586{
2587  matrix result = mpNew(rows,cols);
2588  int i,cp,r=idRankFreeModule(mod),c=IDELEMS(mod);
2589  poly p,h;
2590
2591  if (r>rows) r = rows;
2592  if (c>cols) c = cols;
2593  for(i=0;i<c;i++)
2594  {
2595    p=mod->m[i];
2596    mod->m[i]=NULL;
2597    while (p!=NULL)
2598    {
2599      h=p;
2600      pIter(p);
2601      pNext(h)=NULL;
2602      cp = pGetComp(h);
2603      if (cp<=r)
2604      {
2605        pSetComp(h,0);
2606        pSetmComp(h);
2607        MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
2608      }
2609      else
2610        pDelete(&h);
2611    }
2612  }
2613  idDelete(&mod);
2614  return result;
2615}
2616
2617/*2
2618* substitute the n-th variable by the monomial e in id
2619* destroy id
2620*/
2621ideal  idSubst(ideal id, int n, poly e)
2622{
2623  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
2624  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
2625
2626  res->rank = id->rank;
2627  for(k--;k>=0;k--)
2628  {
2629    res->m[k]=pSubst(id->m[k],n,e);
2630    id->m[k]=NULL;
2631  }
2632  idDelete(&id);
2633  return res;
2634}
2635
2636BOOLEAN idHomModule(ideal m, ideal Q, intvec **w)
2637{
2638  if (w!=NULL) *w=NULL;
2639  if ((Q!=NULL) && (!idHomIdeal(Q,NULL))) return FALSE;
2640  if (idIs0(m)) return TRUE;
2641
2642  int i,j,cmax=2,order=0,ord,* diff,* iscom,diffmin=32000;
2643  poly p=NULL;
2644  int length=IDELEMS(m);
2645  polyset P=m->m;
2646  polyset F=(polyset)Alloc(length*sizeof(poly));
2647  for (i=length-1;i>=0;i--)
2648  {
2649    p=F[i]=P[i];
2650    cmax=max(cmax,pMaxComp(p)+1);
2651  }
2652  diff = (int *)Alloc0(cmax*sizeof(int));
2653  if (w!=NULL) *w=NewIntvec1(cmax-1);
2654  iscom = (int *)Alloc0(cmax*sizeof(int));
2655  i=0;
2656  while (i<=length)
2657  {
2658    if (i<length)
2659    {
2660      p=F[i];
2661      while ((p!=NULL) && (!iscom[pGetComp(p)])) pIter(p);
2662    }
2663    if ((p==NULL) && (i<length))
2664    {
2665      i++;
2666    }
2667    else
2668    {
2669      if (p==NULL)
2670      {
2671        i=0;
2672        while ((i<length) && (F[i]==NULL)) i++;
2673        if (i>=length) break;
2674        p = F[i];
2675      }
2676      if (pLexOrder)
2677        order=pTotaldegree(p);
2678      else
2679      //  order = p->order;
2680        order = pFDeg(p);
2681      order += diff[pGetComp(p)];
2682      p = F[i];
2683//Print("Actual p=F[%d]: ",i);pWrite(p);
2684      F[i] = NULL;
2685      i=0;
2686    }
2687    while (p!=NULL)
2688    {
2689      //if (pLexOrder)
2690      //  ord=pTotaldegree(p);
2691      //else
2692      //  ord = p->order;
2693      ord = pFDeg(p);
2694      if (!iscom[pGetComp(p)])
2695      {
2696        diff[pGetComp(p)] = order-ord;
2697        iscom[pGetComp(p)] = 1;
2698/*
2699*PrintS("new diff: ");
2700*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
2701*PrintLn();
2702*PrintS("new iscom: ");
2703*for (j=0;j<cmax;j++) Print("%d ",iscom[j]);
2704*PrintLn();
2705*Print("new set %d, order %d, ord %d, diff %d\n",pGetComp(p),order,ord,diff[pGetComp(p)]);
2706*/
2707      }
2708      else
2709      {
2710/*
2711*PrintS("new diff: ");
2712*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
2713*PrintLn();
2714*Print("order %d, ord %d, diff %d\n",order,ord,diff[pGetComp(p)]);
2715*/
2716        if (order != ord+diff[pGetComp(p)])
2717        {
2718          Free((ADDRESS) iscom,cmax*sizeof(int));
2719          Free((ADDRESS) diff,cmax*sizeof(int));
2720          Free((ADDRESS) F,length*sizeof(poly));
2721          delete *w;*w=NULL;
2722          return FALSE;
2723        }
2724      }
2725      pIter(p);
2726    }
2727  }
2728  Free((ADDRESS) iscom,cmax*sizeof(int));
2729  Free((ADDRESS) F,length*sizeof(poly));
2730  for (i=1;i<cmax;i++) (**w)[i-1]=diff[i];
2731  for (i=1;i<cmax;i++)
2732  {
2733    if (diff[i]<diffmin) diffmin=diff[i];
2734  }
2735  if (w!=NULL)
2736  {
2737    for (i=1;i<cmax;i++)
2738    {
2739      (**w)[i-1]=diff[i]-diffmin;
2740    }
2741  }
2742  Free((ADDRESS) diff,cmax*sizeof(int));
2743  return TRUE;
2744}
2745
2746ideal idJet(ideal i,int d)
2747{
2748  ideal r=idInit((i->nrows)*(i->ncols),i->rank);
2749  r->nrows = i-> nrows;
2750  r->ncols = i-> ncols;
2751  //r->rank = i-> rank;
2752  int k;
2753  for(k=(i->nrows)*(i->ncols)-1;k>=0; k--)
2754  {
2755    r->m[k]=pJet(i->m[k],d);
2756  }
2757  return r;
2758}
2759
2760ideal idJetW(ideal i,int d, intvec * iv)
2761{
2762  ideal r=idInit(IDELEMS(i),i->rank);
2763  if (ecartWeights!=NULL)
2764  {
2765    WerrorS("cannot compute weighted jets now");
2766  }
2767  else
2768  {
2769    short *w=iv2array(iv);
2770    int k;
2771    for(k=0; k<IDELEMS(i); k++)
2772    {
2773      r->m[k]=pJetW(i->m[k],d,w);
2774    }
2775    Free((ADDRESS)w,(pVariables+1)*sizeof(short));
2776  }
2777  return r;
2778}
2779
2780matrix idDiff(matrix i, int k)
2781{
2782  int e=MATCOLS(i)*MATROWS(i);
2783  matrix r=mpNew(MATROWS(i),MATCOLS(i));
2784  int j;
2785  for(j=0; j<e; j++)
2786  {
2787    r->m[j]=pDiff(i->m[j],k);
2788  }
2789  return r;
2790}
2791
2792matrix idDiffOp(ideal I, ideal J,BOOLEAN multiply)
2793{
2794  matrix r=mpNew(IDELEMS(I),IDELEMS(J));
2795  int i,j;
2796  for(i=0; i<IDELEMS(I); i++)
2797  {
2798    for(j=0; j<IDELEMS(J); j++)
2799    {
2800      MATELEM(r,i+1,j+1)=pDiffOp(I->m[i],J->m[j],multiply);
2801    }
2802  }
2803  return r;
2804}
2805
2806/*3
2807*handles for some ideal operations the ring/syzcomp managment
2808*returns all syzygies (componentwise-)shifted by -syzcomp
2809*or -syzcomp-1 (in case of ideals as input)
2810static ideal idHandleIdealOp(ideal arg,int syzcomp,int isIdeal=FALSE)
2811{
2812  ring orig_ring=currRing;
2813  ring syz_ring=rCurrRingAssureSyzComp();
2814  rSetSyzComp(length);
2815
2816  ideal s_temp;
2817  if (orig_ring!=syz_ring)
2818    s_temp=idrMoveR_NoSort(arg,orig_ring);
2819  else
2820    s_temp=arg;
2821
2822  ideal s_temp1 = kStd(s_temp,currQuotient,testHomog,&w,NULL,length);
2823  if (w!=NULL) delete w;
2824
2825  if (syz_ring!=orig_ring)
2826  {
2827    idDelete(&s_temp);
2828    rChangeCurrRing(orig_ring,TRUE);
2829  }
2830
2831  idDelete(&temp);
2832  ideal temp1=idRingCopy(s_temp1,syz_ring);
2833
2834  if (syz_ring!=orig_ring)
2835  {
2836    rChangeCurrRing(syz_ring,FALSE);
2837    idDelete(&s_temp1);
2838    rChangeCurrRing(orig_ring,TRUE);
2839    rKill(syz_ring);
2840  }
2841
2842  for (i=0;i<IDELEMS(temp1);i++)
2843  {
2844    if ((temp1->m[i]!=NULL)
2845    && (pGetComp(temp1->m[i])<=length))
2846    {
2847      pDelete(&(temp1->m[i]));
2848    }
2849    else
2850    {
2851      pShift(&(temp1->m[i]),-length);
2852    }
2853  }
2854  temp1->rank = rk;
2855  idSkipZeroes(temp1);
2856
2857  return temp1;
2858}
2859*/
2860/*2
2861* represents (h1+h2)/h2=h1/(h1 intersect h2)
2862*/
2863ideal idModulo (ideal h2,ideal h1)
2864{
2865  int i,j,k,rk,flength=0,slength,length;
2866  intvec * w;
2867  poly p,q;
2868
2869  if (idIs0(h2))
2870    return idFreeModule(max(1,h2->ncols));
2871  if (!idIs0(h1))
2872    flength = idRankFreeModule(h1);
2873  slength = idRankFreeModule(h2);
2874  length  = max(flength,slength);
2875  if (length==0)
2876  {
2877    length = 1;
2878  }
2879  ideal temp = idInit(IDELEMS(h2),length+IDELEMS(h2));
2880  for (i=0;i<IDELEMS(h2);i++)
2881  {
2882    temp->m[i] = pCopy(h2->m[i]);
2883    q = pOne();
2884    pSetComp(q,i+1+length);
2885    pSetmComp(q);
2886    if(temp->m[i]!=NULL)
2887    {
2888      if (slength==0) pShift(&(temp->m[i]),1);
2889      p = temp->m[i];
2890      while (pNext(p)!=NULL) pIter(p);
2891      pNext(p) = q;
2892    }
2893    else
2894      temp->m[i]=q;
2895  }
2896  rk = k = IDELEMS(h2);
2897  if (!idIs0(h1))
2898  {
2899    pEnlargeSet(&(temp->m),IDELEMS(temp),IDELEMS(h1));
2900    IDELEMS(temp) += IDELEMS(h1);
2901    for (i=0;i<IDELEMS(h1);i++)
2902    {
2903      if (h1->m[i]!=NULL)
2904      {
2905        temp->m[k] = pCopy(h1->m[i]);
2906        if (flength==0) pShift(&(temp->m[k]),1);
2907        k++;
2908      }
2909    }
2910  }
2911
2912  ring orig_ring=currRing;
2913  ring syz_ring=rCurrRingAssureSyzComp();
2914  rSetSyzComp(length);
2915  ideal s_temp;
2916 
2917  if (syz_ring != orig_ring)
2918  {
2919    s_temp = idrCopyR_NoSort(temp, orig_ring);
2920  }
2921  else
2922  {
2923    s_temp = temp;
2924  }
2925 
2926  ideal s_temp1 = kStd(s_temp,currQuotient,testHomog,&w,NULL,length);
2927  if (w!=NULL) delete w;
2928
2929  for (i=0;i<IDELEMS(s_temp1);i++)
2930  {
2931    if ((s_temp1->m[i]!=NULL)
2932    && (pGetComp(s_temp1->m[i])<=length))
2933    {
2934      pDelete(&(s_temp1->m[i]));
2935    }
2936    else
2937    {
2938      pShift(&(s_temp1->m[i]),-length);
2939    }
2940  }
2941  s_temp1->rank = rk;
2942  idSkipZeroes(s_temp1);
2943
2944  if (syz_ring!=orig_ring)
2945  {
2946    rChangeCurrRing(orig_ring,TRUE);
2947    s_temp = idrMoveR_NoSort(s_temp, syz_ring);
2948    idDelete(&temp);
2949    rKill(syz_ring);
2950  }
2951  return s_temp1;
2952}
2953
2954int idElem(ideal F)
2955{
2956  int i=0,j=0;
2957
2958  while(j<IDELEMS(F))
2959  {
2960   if ((F->m)[j]!=NULL) i++;
2961   j++;
2962  }
2963  return i;
2964}
2965
2966/*
2967*computes module-weights for liftings of homogeneous modules
2968*/
2969intvec * idMWLift(ideal mod,intvec * weights)
2970{
2971  if (idIs0(mod)) return NewIntvec1(2);
2972  int i=IDELEMS(mod);
2973  while ((i>0) && (mod->m[i-1]==NULL)) i--;
2974  intvec *result = NewIntvec1(i+1);
2975  while (i>0)
2976  {
2977    (*result)[i]=pFDeg(mod->m[i])+(*weights)[pGetComp(mod->m[i])];
2978  }
2979  return result;
2980}
2981
2982/*2
2983*sorts the kbase for idCoef* in a special way (lexicographically
2984*with x_max,...,x_1)
2985*/
2986ideal idCreateSpecialKbase(ideal kBase,intvec ** convert)
2987{
2988  int i;
2989  ideal result;
2990
2991  if (idIs0(kBase)) return NULL;
2992  result = idInit(IDELEMS(kBase),kBase->rank);
2993  *convert = idSort(kBase,FALSE);
2994  for (i=0;i<(*convert)->length();i++)
2995  {
2996    result->m[i] = pCopy(kBase->m[(**convert)[i]-1]);
2997  }
2998  return result;
2999}
3000
3001/*2
3002*returns the index of a given monom in the list of the special kbase
3003*/
3004int idIndexOfKBase(poly monom, ideal kbase)
3005{
3006  int j=IDELEMS(kbase);
3007
3008  while ((j>0) && (kbase->m[j-1]==NULL)) j--;
3009  if (j==0) return -1;
3010  int i=pVariables;
3011  while (i>0)
3012  {
3013    loop
3014    {
3015      if (pGetExp(monom,i)>pGetExp(kbase->m[j-1],i)) return -1;
3016      if (pGetExp(monom,i)==pGetExp(kbase->m[j-1],i)) break;
3017      j--;
3018      if (j==0) return -1;
3019    }
3020    if (i==1)
3021    {
3022      while(j>0)
3023      {
3024        if (pGetComp(monom)==pGetComp(kbase->m[j-1])) return j-1;
3025        if (pGetComp(monom)>pGetComp(kbase->m[j-1])) return -1;
3026        j--;
3027      }
3028    }
3029    i--;
3030  }
3031  return -1;
3032}
3033
3034/*2
3035*decomposes the monom in a part of coefficients described by the
3036*complement of how and a monom in variables occuring in how, the
3037*index of which in kbase is returned as integer pos (-1 if it don't
3038*exists)
3039*/
3040poly idDecompose(poly monom, poly how, ideal kbase, int * pos)
3041{
3042  int i;
3043  poly coeff=pOne(), base=pOne();
3044
3045  for (i=1;i<=pVariables;i++)
3046  {
3047    if (pGetExp(how,i)>0)
3048    {
3049      pSetExp(base,i,pGetExp(monom,i));
3050    }
3051    else
3052    {
3053      pSetExp(coeff,i,pGetExp(monom,i));
3054    }
3055  }
3056  pSetComp(base,pGetComp(monom));
3057  pSetm(base);
3058  pSetCoeff(coeff,nCopy(pGetCoeff(monom)));
3059  pSetm(coeff);
3060  *pos = idIndexOfKBase(base,kbase);
3061  if (*pos<0)
3062    pDelete(&coeff);
3063  pDelete(&base);
3064  return coeff;
3065}
3066
3067/*2
3068*returns a matrix A of coefficients with kbase*A=arg
3069*if all monomials in variables of how occur in kbase
3070*the other are deleted
3071*/
3072matrix idCoeffOfKBase(ideal arg, ideal kbase, poly how)
3073{
3074  matrix result;
3075  ideal tempKbase;
3076  poly p,q;
3077  intvec * convert;
3078  int i=IDELEMS(kbase),j=IDELEMS(arg),k,pos;
3079#if 0
3080  while ((i>0) && (kbase->m[i-1]==NULL)) i--;
3081  if (idIs0(arg))
3082    return mpNew(i,1);
3083  while ((j>0) && (arg->m[j-1]==NULL)) j--;
3084  result = mpNew(i,j);
3085#else
3086  result = mpNew(i, j);
3087  while ((j>0) && (arg->m[j-1]==NULL)) j--;
3088#endif
3089
3090  tempKbase = idCreateSpecialKbase(kbase,&convert);
3091  for (k=0;k<j;k++)
3092  {
3093    p = arg->m[k];
3094    while (p!=NULL)
3095    {
3096      q = idDecompose(p,how,tempKbase,&pos);
3097      if (pos>=0)
3098      {
3099        MATELEM(result,(*convert)[pos],k+1) =
3100            pAdd(MATELEM(result,(*convert)[pos],k+1),q);
3101      }
3102      else
3103        pDelete(&q);
3104      pIter(p);
3105    }
3106  }
3107  idDelete(&tempKbase);
3108  return result;
3109}
3110
3111intvec * idQHomWeights(ideal id)
3112{
3113  intvec * imat=NewIntvec3(2*pVariables,pVariables,0);
3114  poly actHead=NULL,wPoint=NULL;
3115  int actIndex,i=-1,j=1,k;
3116  BOOLEAN notReady=TRUE;
3117
3118  while (notReady)
3119  {
3120    if (wPoint==NULL)
3121    {
3122      i++;
3123      while ((i<IDELEMS(id))
3124      && ((id->m[i]==NULL) || (pNext(id->m[i])==NULL)))
3125        i++;
3126      if (i<IDELEMS(id))
3127      {
3128        actHead = id->m[i];
3129        wPoint = pNext(actHead);
3130      }
3131    }
3132    while ((wPoint!=NULL) && (j<=2*pVariables))
3133    {
3134      for (k=1;k<=pVariables;k++)
3135        IMATELEM(*imat,j,k) += pGetExp(actHead,k)-pGetExp(wPoint,k);
3136      pIter(wPoint);
3137      j++;
3138    }
3139    if ((i>=IDELEMS(id)) || (j>2*pVariables))
3140    {
3141      ivTriangMat(imat,1,1);
3142      j = ivFirstEmptyRow(imat);
3143      if ((i>=IDELEMS(id)) || (j>pVariables)) notReady=FALSE;
3144    }
3145  }
3146  intvec *result=NULL;
3147  if (j<=pVariables)
3148  {
3149    result=ivSolveIntMat(imat);
3150  }
3151  //else
3152  //{
3153  //  WerrorS("not homogeneous");
3154  //}
3155  delete imat;
3156  return result;
3157}
3158
3159/*3
3160* searches for units in the components of the module arg and
3161* returns the first one
3162*/
3163static int idReadOutUnits(ideal arg,int* comp)
3164{
3165  if (idIs0(arg)) return -1;
3166  int i=0,j,rk_arg=idRankFreeModule(arg),generator=-1;
3167  intvec * componentIsUsed =new intvec(rk_arg+1);
3168  poly p,q;
3169
3170  while ((i<IDELEMS(arg)) && (generator<0))
3171  {
3172    for (j=rk_arg;j>=0;j--) 
3173      (*componentIsUsed)[j]=0;
3174    p = arg->m[i];
3175    while (p!=NULL)
3176    { 
3177      j = pGetComp(p);
3178      if ((*componentIsUsed)[j]==0)
3179      {
3180        if (pIsConstantComp(p))
3181        {
3182          generator = i;
3183          (*componentIsUsed)[j] = 1;
3184        }
3185        else
3186        {
3187          (*componentIsUsed)[j] = -1;
3188        }
3189      }
3190      else if ((*componentIsUsed)[j]>0)
3191      {
3192        ((*componentIsUsed)[j])++;
3193      }
3194      pIter(p);
3195    }
3196    i++;
3197  }
3198  i = 0;
3199  *comp = -1;
3200  for (j=0;j<=rk_arg;j++)
3201  {
3202    if ((*componentIsUsed)[j]>0)
3203    {
3204      if ((*comp==-1) || ((*componentIsUsed)[j]<i))
3205      {
3206        *comp = j;
3207        i= (*componentIsUsed)[j];
3208      }
3209    }
3210  }
3211  return generator;
3212}
3213
3214static void idDeleteComp(ideal arg,int red_comp)
3215{
3216  int i,j;
3217  poly p;
3218
3219  for (i=IDELEMS(arg)-1;i>=0;i--)
3220  {
3221    p = arg->m[i];
3222    while (p!=NULL)
3223    {
3224      j = pGetComp(p);
3225      if (j>red_comp)
3226      {
3227        pSetComp(p,j-1);
3228        pSetm(p);
3229      }
3230      pIter(p);
3231    }
3232  }
3233  (arg->rank)--;
3234}
3235
3236/*2
3237* returns the presentation of an isomorphic, minimally
3238* embedded  module
3239*/
3240ideal idMinEmbedding(ideal arg,BOOLEAN inPlace)
3241{
3242  if (idIs0(arg)) return idInit(1,arg->rank);
3243  int next_gen,next_comp;
3244  ideal res=arg;
3245
3246  if (!inPlace) res = idCopy(arg);
3247  loop
3248  {
3249    next_gen = idReadOutUnits(res,&next_comp);
3250    if (next_gen<0) break;
3251    syGaussForOne(res,next_gen,next_comp,0,IDELEMS(res));
3252    idDeleteComp(res,next_comp);
3253  }
3254  idSkipZeroes(res);
3255  return res;
3256}
3257
3258/*2
3259* transpose a module
3260*/
3261ideal idTransp(ideal a)
3262{
3263  int r = a->rank, c = IDELEMS(a);
3264  ideal b =  idInit(r,c);
3265
3266  for (int i=c; i>0; i--)
3267  {
3268    poly p=a->m[i-1];
3269    while(p!=NULL)
3270    {
3271      poly h=pHead(p);
3272      int co=pGetComp(h)-1;
3273      pSetComp(h,i);
3274      pSetmComp(h);
3275      b->m[co]=pAdd(b->m[co],h);
3276      pIter(p);
3277    }
3278  }
3279  return b;
3280}
3281
Note: See TracBrowser for help on using the repository browser.