source: git/kernel/ideals.cc @ 62dd9b

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