source: git/kernel/ideals.cc @ ce89df

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