source: git/kernel/ideals.cc @ 42d5aa

spielwiese
Last change on this file since 42d5aa was 42d5aa, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: pFDeg git-svn-id: file:///usr/local/Singular/svn/trunk@7785 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 71.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: ideals.cc,v 1.8 2005-03-17 14:13:02 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   i, j, k, length=0,reg;
1387  BOOLEAN isMonomial=TRUE;
1388  int ii;
1389
1390#ifdef PDEBUG
1391  for(ii=0;ii<IDELEMS(h1);ii++) pTest(h1->m[ii]);
1392#endif
1393  if (idIs0(h1))
1394  {
1395    ideal result=idFreeModule(IDELEMS(h1));
1396    int curr_syz_limit=rGetCurrSyzLimit();
1397    if (curr_syz_limit>0)
1398    for (ii=0;ii<IDELEMS(h1);ii++)
1399    {
1400      if (h1->m[ii]!=NULL)
1401        pShift(&h1->m[ii],curr_syz_limit);
1402    }
1403    return result;
1404  }
1405  k=si_max(1,(int)idRankFreeModule(h1));
1406
1407  assume(currRing != NULL);
1408  ring orig_ring=currRing;
1409  ring syz_ring=rCurrRingAssure_SyzComp();
1410
1411  if (setSyzComp)
1412    rSetSyzComp(k);
1413
1414  if (orig_ring != syz_ring)
1415  {
1416    s_h1=idrCopyR_NoSort(h1,orig_ring);
1417  }
1418  else
1419  {
1420    s_h1 = h1;
1421  }
1422
1423  ideal s_h3=idPrepare(s_h1,h,k,w);
1424
1425  if (s_h3==NULL)
1426  {
1427    return idFreeModule(IDELEMS(h1));
1428  }
1429
1430  if (orig_ring != syz_ring)
1431  {
1432    idDelete(&s_h1);
1433    for (j=0; j<IDELEMS(s_h3); j++)
1434    {
1435      if (s_h3->m[j] != NULL)
1436      {
1437        if (p_MinComp(s_h3->m[j],syz_ring) > k)
1438          pShift(&s_h3->m[j], -k);
1439        else
1440          pDelete(&s_h3->m[j]);
1441      }
1442    }
1443    idSkipZeroes(s_h3);
1444    s_h3->rank -= k;
1445    rChangeCurrRing(orig_ring);
1446    s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
1447    rKill(syz_ring);
1448    idTest(s_h3);
1449    return s_h3;
1450  }
1451
1452  ideal e = idInit(IDELEMS(s_h3), s_h3->rank);
1453
1454  for (j=0; j<IDELEMS(s_h3); j++)
1455  {
1456    if (s_h3->m[j] != NULL)
1457    {
1458      if (p_MinComp(s_h3->m[j],syz_ring) <= k)
1459      {
1460        e->m[j] = s_h3->m[j];
1461        isMonomial=isMonomial && (pNext(s_h3->m[j])==NULL);
1462        pDelete(&pNext(s_h3->m[j]));
1463        s_h3->m[j] = NULL;
1464      }
1465    }
1466  }
1467
1468  idSkipZeroes(s_h3);
1469  idSkipZeroes(e);
1470
1471  if ((deg != NULL)
1472  && (!isMonomial)
1473  && (!TEST_OPT_NOTREGULARITY)
1474  && (setRegularity)
1475  && (h==isHomog)
1476  && (!rIsPluralRing(currRing))
1477  )
1478  {
1479    ring dp_C_ring = rCurrRingAssure_dp_C();
1480    if (dp_C_ring != syz_ring)
1481      e = idrMoveR_NoSort(e, syz_ring);
1482    resolvente res = sySchreyerResolvente(e,-1,&length,TRUE, TRUE);
1483    intvec * dummy = syBetti(res,length,&reg, *w);
1484    *deg = reg+2;
1485    delete dummy;
1486    for (j=0;j<length;j++)
1487    {
1488      if (res[j]!=NULL) idDelete(&(res[j]));
1489    }
1490    omFreeSize((ADDRESS)res,length*sizeof(ideal));
1491    idDelete(&e);
1492    if (dp_C_ring != syz_ring)
1493    {
1494      rChangeCurrRing(syz_ring);
1495      rKill(dp_C_ring);
1496    }
1497  }
1498  else
1499  {
1500    idDelete(&e);
1501  }
1502  idTest(s_h3);
1503  if (currQuotient != NULL)
1504  {
1505    ideal ts_h3=kStd(s_h3,currQuotient,h,w);
1506    idDelete(&s_h3);
1507    s_h3 = ts_h3;
1508  }
1509  return s_h3;
1510}
1511
1512/*2
1513*/
1514ideal idXXX (ideal  h1, int k)
1515{
1516  ideal s_h1;
1517  int j;
1518  intvec *w=NULL;
1519
1520  assume(currRing != NULL);
1521  ring orig_ring=currRing;
1522  ring syz_ring=rCurrRingAssure_SyzComp();
1523
1524  rSetSyzComp(k);
1525
1526  if (orig_ring != syz_ring)
1527  {
1528    s_h1=idrCopyR_NoSort(h1,orig_ring);
1529  }
1530  else
1531  {
1532    s_h1 = h1;
1533  }
1534
1535  ideal s_h3=kStd(s_h1,NULL,testHomog,&w,NULL,k);
1536
1537  if (s_h3==NULL)
1538  {
1539    return idFreeModule(IDELEMS(h1));
1540  }
1541
1542  if (orig_ring != syz_ring)
1543  {
1544    idDelete(&s_h1);
1545    idSkipZeroes(s_h3);
1546    rChangeCurrRing(orig_ring);
1547    s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
1548    rKill(syz_ring);
1549    idTest(s_h3);
1550    return s_h3;
1551  }
1552
1553  idSkipZeroes(s_h3);
1554  idTest(s_h3);
1555  return s_h3;
1556}
1557
1558/*
1559*computes a standard basis for h1 and stores the transformation matrix
1560* in ma
1561*/
1562ideal idLiftStd (ideal  h1, matrix* ma, tHomog h)
1563{
1564  int   i, j, k, t, inputIsIdeal=idRankFreeModule(h1);
1565  poly  p=NULL, q, qq;
1566  intvec *w=NULL;
1567
1568  idDelete((ideal*)ma);
1569  *ma=mpNew(1,0);
1570  if (idIs0(h1))
1571    return idInit(1,h1->rank);
1572  k=si_max(1,(int)idRankFreeModule(h1));
1573
1574  ring orig_ring=currRing;
1575  ring syz_ring=rCurrRingAssure_SyzComp();
1576  rSetSyzComp(k);
1577
1578  ideal s_h1=h1;
1579
1580  if (orig_ring != syz_ring)
1581    s_h1 = idrCopyR_NoSort(h1,orig_ring);
1582  else
1583    s_h1 = h1;
1584
1585  ideal s_h3=idPrepare(s_h1,h,k,&w);
1586  ideal s_h2 = idInit(IDELEMS(s_h3), s_h3->rank);
1587
1588  if (w!=NULL) delete w;
1589  i = 0;
1590  for (j=0; j<IDELEMS(s_h3); j++)
1591  {
1592    if ((s_h3->m[j] != NULL) && (p_MinComp(s_h3->m[j],syz_ring) <= k))
1593    {
1594      i++;
1595      q = s_h3->m[j];
1596      while (pNext(q) != NULL)
1597      {
1598        if (pGetComp(pNext(q)) > k)
1599        {
1600          s_h2->m[j] = pNext(q);
1601          pNext(q) = NULL;
1602        }
1603        else
1604        {
1605          pIter(q);
1606        }
1607      }
1608      if (!inputIsIdeal) pShift(&(s_h3->m[j]), -1);
1609    }
1610    else
1611    {
1612      pDelete(&(s_h3->m[j]));
1613    }
1614  }
1615
1616  idSkipZeroes(s_h3);
1617  j = IDELEMS(s_h1);
1618
1619  if (syz_ring!=orig_ring)
1620  {
1621    idDelete(&s_h1);
1622    rChangeCurrRing(orig_ring);
1623  }
1624
1625  idDelete((ideal*)ma);
1626  *ma = mpNew(j,i);
1627
1628  i = 1;
1629  for (j=0; j<IDELEMS(s_h2); j++)
1630  {
1631    if (s_h2->m[j] != NULL)
1632    {
1633      q = prMoveR( s_h2->m[j], syz_ring);
1634      s_h2->m[j] = NULL;
1635
1636      while (q != NULL)
1637      {
1638        p = q;
1639        pIter(q);
1640        pNext(p) = NULL;
1641        t=pGetComp(p);
1642        pSetComp(p,0);
1643        pSetmComp(p);
1644        MATELEM(*ma,t-k,i) = pAdd(MATELEM(*ma,t-k,i),p);
1645      }
1646      i++;
1647    }
1648  }
1649  idDelete(&s_h2);
1650
1651  for (i=0; i<IDELEMS(s_h3); i++)
1652  {
1653    s_h3->m[i] = prMoveR_NoSort(s_h3->m[i], syz_ring);
1654  }
1655
1656  if (syz_ring!=orig_ring) rKill(syz_ring);
1657  return s_h3;
1658}
1659
1660static void idPrepareStd(ideal s_temp, int k)
1661{
1662  int j,rk=idRankFreeModule(s_temp);
1663  poly p,q;
1664
1665  if (rk == 0)
1666  {
1667    for (j=0; j<IDELEMS(s_temp); j++)
1668    {
1669      if (s_temp->m[j]!=NULL) pSetCompP(s_temp->m[j],1);
1670    }
1671    k = si_max(k,1);
1672  }
1673  for (j=0; j<IDELEMS(s_temp); j++)
1674  {
1675    if (s_temp->m[j]!=NULL)
1676    {
1677      p = s_temp->m[j];
1678      q = pOne();
1679      //pGetCoeff(q)=nNeg(pGetCoeff(q));   //set q to -1
1680      pSetComp(q,k+1+j);
1681      pSetmComp(q);
1682      while (pNext(p)) pIter(p);
1683      pNext(p) = q;
1684    }
1685  }
1686}
1687
1688/*2
1689*computes a representation of the generators of submod with respect to those
1690* of mod
1691*/
1692
1693ideal idLift(ideal mod, ideal submod,ideal *rest, BOOLEAN goodShape,
1694             BOOLEAN isSB, BOOLEAN divide, matrix *unit)
1695{
1696  int lsmod =idRankFreeModule(submod), i, j, k;
1697  int comps_to_add=0;
1698  poly p;
1699
1700  if (idIs0(submod))
1701  {
1702    if (unit!=NULL)
1703    {
1704      *unit=mpNew(1,1);
1705      MATELEM(*unit,1,1)=pOne();
1706    }
1707    if (rest!=NULL)
1708    {
1709      *rest=idInit(1,mod->rank);
1710    }
1711    return idInit(1,mod->rank);
1712  }
1713  if (idIs0(mod))
1714  {
1715    if (unit!=NULL)
1716    {
1717      i=IDELEMS(submod);
1718      *unit=mpNew(i,i);
1719      for (j=i;j>0;j--)
1720      {
1721        MATELEM(*unit,j,j)=pOne();
1722      }
1723    }
1724    if (rest!=NULL)
1725    {
1726      *rest=idCopy(submod);
1727    }
1728    return idInit(1,mod->rank);
1729  }
1730  if (unit!=NULL)
1731  {
1732    comps_to_add = IDELEMS(submod);
1733    while ((comps_to_add>0) && (submod->m[comps_to_add-1]==NULL))
1734      comps_to_add--;
1735  }
1736  k=idRankFreeModule(mod);
1737  if  ((k!=0) && (lsmod==0)) lsmod=1;
1738  k=si_max(k,1);
1739
1740  ring orig_ring=currRing;
1741  ring syz_ring=rCurrRingAssure_SyzComp();
1742  rSetSyzComp(k);
1743
1744  ideal s_mod, s_temp;
1745  if (orig_ring != syz_ring)
1746  {
1747    s_mod = idrCopyR_NoSort(mod,orig_ring);
1748    s_temp = idrCopyR_NoSort(submod,orig_ring);
1749  }
1750  else
1751  {
1752    s_mod = mod;
1753    s_temp = idCopy(submod);
1754  }
1755  ideal s_h3;
1756  if (isSB)
1757  {
1758    s_h3 = idCopy(s_mod);
1759    idPrepareStd(s_h3, k+comps_to_add);
1760  }
1761  else
1762  {
1763    s_h3 = idPrepare(s_mod,(tHomog)FALSE,k+comps_to_add,NULL);
1764  }
1765  if (!goodShape)
1766  {
1767    for (j=0;j<IDELEMS(s_h3);j++)
1768    {
1769      if ((s_h3->m[j] != NULL) && (pMinComp(s_h3->m[j]) > k))
1770        pDelete(&(s_h3->m[j]));
1771    }
1772  }
1773  idSkipZeroes(s_h3);
1774  if (lsmod==0)
1775  {
1776    for (j=IDELEMS(s_temp);j>0;j--)
1777    {
1778      if (s_temp->m[j-1]!=NULL)
1779        pShift(&(s_temp->m[j-1]),1);
1780    }
1781  }
1782  if (unit!=NULL)
1783  {
1784    for(j = 0;j<comps_to_add;j++)
1785    {
1786      p = s_temp->m[j];
1787      if (p!=NULL)
1788      {
1789        while (pNext(p)!=NULL) pIter(p);
1790        pNext(p) = pOne();
1791        pIter(p);
1792        pSetComp(p,1+j+k);
1793        pSetmComp(p);
1794        p = pNeg(p);
1795      }
1796    }
1797  }
1798  ideal s_result = kNF(s_h3,currQuotient,s_temp,k);
1799  s_result->rank = s_h3->rank;
1800  ideal s_rest = idInit(IDELEMS(s_result),k);
1801  idDelete(&s_h3);
1802  idDelete(&s_temp);
1803
1804  for (j=0;j<IDELEMS(s_result);j++)
1805  {
1806    if (s_result->m[j]!=NULL)
1807    {
1808      if (pGetComp(s_result->m[j])<=k)
1809      {
1810        if (!divide)
1811        {
1812          if (isSB)
1813          {
1814            WarnS("first module not a standardbasis\n"
1815              "// ** or second not a proper submodule");
1816          }
1817          else
1818            WerrorS("2nd module lies not in the first");
1819          idDelete(&s_result);
1820          idDelete(&s_rest);
1821          s_result=idInit(IDELEMS(submod),submod->rank);
1822          break;
1823        }
1824        else
1825        {
1826          p = s_rest->m[j] = s_result->m[j];
1827          while ((pNext(p)!=NULL) && (pGetComp(pNext(p))<=k)) pIter(p);
1828          s_result->m[j] = pNext(p);
1829          pNext(p) = NULL;
1830        }
1831      }
1832      pShift(&(s_result->m[j]),-k);
1833      pNeg(s_result->m[j]);
1834    }
1835  }
1836  if ((lsmod==0) && (!idIs0(s_rest)))
1837  {
1838    for (j=IDELEMS(s_rest);j>0;j--)
1839    {
1840      if (s_rest->m[j-1]!=NULL)
1841      {
1842        pShift(&(s_rest->m[j-1]),-1);
1843        s_rest->m[j-1] = s_rest->m[j-1];
1844      }
1845    }
1846  }
1847  if(syz_ring!=orig_ring)
1848  {
1849    idDelete(&s_mod);
1850    rChangeCurrRing(orig_ring);
1851    s_result = idrMoveR_NoSort(s_result, syz_ring);
1852    s_rest = idrMoveR_NoSort(s_rest, syz_ring);
1853    rKill(syz_ring);
1854  }
1855  if (rest!=NULL)
1856    *rest = s_rest;
1857  else
1858    idDelete(&s_rest);
1859//idPrint(s_result);
1860  if (unit!=NULL)
1861  {
1862    *unit=mpNew(comps_to_add,comps_to_add);
1863    int i;
1864    for(i=0;i<IDELEMS(s_result);i++)
1865    {
1866      poly p=s_result->m[i];
1867      poly q=NULL;
1868      while(p!=NULL)
1869      {
1870        if(pGetComp(p)<=comps_to_add)
1871        {
1872          pSetComp(p,0);
1873          if (q!=NULL)
1874          {
1875            pNext(q)=pNext(p);
1876          }
1877          else
1878          {
1879            pIter(s_result->m[i]);
1880          }
1881          pNext(p)=NULL;
1882          MATELEM(*unit,i+1,i+1)=pAdd(MATELEM(*unit,i+1,i+1),p);
1883          if(q!=NULL)   p=pNext(q);
1884          else          p=s_result->m[i];
1885        }
1886        else
1887        {
1888          q=p;
1889          pIter(p);
1890        }
1891      }
1892      pShift(&s_result->m[i],-comps_to_add);
1893    }
1894  }
1895  return s_result;
1896}
1897
1898/*2
1899*computes division of P by Q with remainder up to (w-weighted) degree n
1900*P, Q, and w are not changed
1901*/
1902void idLiftW(ideal P,ideal Q,int n,matrix &T, ideal &R,short *w)
1903{
1904  long N=0;
1905  int i;
1906  for(i=IDELEMS(Q)-1;i>=0;i--)
1907    if(w==NULL)
1908      N=si_max(N,pDeg(Q->m[i]));
1909    else
1910      N=si_max(N,pDegW(Q->m[i],w));
1911  N+=n;
1912
1913  T=mpNew(IDELEMS(Q),IDELEMS(P));
1914  R=idInit(IDELEMS(P),P->rank);
1915
1916  for(i=IDELEMS(P)-1;i>=0;i--)
1917  {
1918    poly p;
1919    if(w==NULL)
1920      p=ppJet(P->m[i],N);
1921    else
1922      p=ppJetW(P->m[i],N,w);
1923
1924    int j=IDELEMS(Q)-1;
1925    while(p!=NULL)
1926    {
1927      if(pDivisibleBy(Q->m[j],p))
1928      {
1929        poly p0=pDivideM(pHead(p),pHead(Q->m[j]));
1930        if(w==NULL)
1931          p=pJet(pSub(p,ppMult_mm(Q->m[j],p0)),N);
1932        else
1933          p=pJetW(pSub(p,ppMult_mm(Q->m[j],p0)),N,w);
1934        pNormalize(p);
1935        if(w==NULL&&pDeg(p0)>n||w!=NULL&&pDegW(p0,w)>n)
1936          pDelete(&p0);
1937        else
1938          MATELEM(T,j+1,i+1)=pAdd(MATELEM(T,j+1,i+1),p0);
1939        j=IDELEMS(Q)-1;
1940      }
1941      else
1942      {
1943        if(j==0)
1944        {
1945          poly p0=p;
1946          pIter(p);
1947          pNext(p0)=NULL;
1948          if(w==NULL&&pDeg(p0)>n||w!=NULL&&pDegW(p0,w)>n)
1949            pDelete(&p0);
1950          else
1951            R->m[i]=pAdd(R->m[i],p0);
1952          j=IDELEMS(Q)-1;
1953        }
1954        else
1955          j--;
1956      }
1957    }
1958  }
1959}
1960
1961/*2
1962*computes the quotient of h1,h2 : internal routine for idQuot
1963*BEWARE: the returned ideals may contain incorrectly ordered polys !
1964*
1965*/
1966static ideal idInitializeQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb,
1967                               BOOLEAN *addOnlyOne, int *kkmax)
1968{
1969  ideal temph1;
1970  poly     p,q = NULL;
1971  int i,l,ll,k,kkk,kmax;
1972  int j = 0;
1973  int k1 = idRankFreeModule(h1);
1974  int k2 = idRankFreeModule(h2);
1975  tHomog   hom=isNotHomog;
1976
1977  k=si_max(k1,k2);
1978  if (k==0)
1979    k = 1;
1980  if ((k2==0) && (k>1)) *addOnlyOne = FALSE;
1981
1982  intvec * weights;
1983  hom = (tHomog)idHomModule(h1,currQuotient,&weights);
1984  if (addOnlyOne && (!h1IsStb))
1985    temph1 = kStd(h1,currQuotient,hom,&weights,NULL);
1986  else
1987    temph1 = idCopy(h1);
1988  if (weights!=NULL) delete weights;
1989  idTest(temph1);
1990/*--- making a single vector from h2 ---------------------*/
1991  for (i=0; i<IDELEMS(h2); i++)
1992  {
1993    if (h2->m[i] != NULL)
1994    {
1995      p = pCopy(h2->m[i]);
1996      if (k2 == 0)
1997        pShift(&p,j*k+1);
1998      else
1999        pShift(&p,j*k);
2000      q = pAdd(q,p);
2001      j++;
2002    }
2003  }
2004  *kkmax = kmax = j*k+1;
2005/*--- adding a monomial for the result (syzygy) ----------*/
2006  p = q;
2007  while (pNext(p)!=NULL) pIter(p);
2008  pNext(p) = pOne();
2009  pIter(p);
2010  pSetComp(p,kmax);
2011  pSetmComp(p);
2012/*--- constructing the big matrix ------------------------*/
2013  ideal h4 = idInit(16,kmax+k-1);
2014  h4->m[0] = q;
2015  if (k2 == 0)
2016  {
2017    if (k > IDELEMS(h4))
2018    {
2019      pEnlargeSet(&(h4->m),IDELEMS(h4),k-IDELEMS(h4));
2020      IDELEMS(h4) = k;
2021    }
2022    for (i=1; i<k; i++)
2023    {
2024      p = pCopy_noCheck(h4->m[i-1]);
2025      pShift(&p,1);
2026      h4->m[i] = p;
2027    }
2028  }
2029
2030  kkk = IDELEMS(h4);
2031  i = IDELEMS(temph1);
2032  while ((i>0) && (temph1->m[i-1]==NULL)) i--;
2033  for (l=0; l<i; l++)
2034  {
2035    if(temph1->m[l]!=NULL)
2036    {
2037      for (ll=0; ll<j; ll++)
2038      {
2039        p = pCopy(temph1->m[l]);
2040        if (k1 == 0)
2041          pShift(&p,ll*k+1);
2042        else
2043          pShift(&p,ll*k);
2044        if (kkk >= IDELEMS(h4))
2045        {
2046          pEnlargeSet(&(h4->m),IDELEMS(h4),16);
2047          IDELEMS(h4) += 16;
2048        }
2049        h4->m[kkk] = p;
2050        kkk++;
2051      }
2052    }
2053  }
2054/*--- if h2 goes in as single vector - the h1-part is just SB ---*/
2055  if (*addOnlyOne)
2056  {
2057    p = h4->m[0];
2058    for (i=0;i<IDELEMS(h4)-1;i++)
2059    {
2060      h4->m[i] = h4->m[i+1];
2061    }
2062    h4->m[IDELEMS(h4)-1] = p;
2063    idSkipZeroes(h4);
2064    test |= Sy_bit(OPT_SB_1);
2065  }
2066  idDelete(&temph1);
2067  return h4;
2068}
2069/*2
2070*computes the quotient of h1,h2
2071*/
2072ideal idQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb, BOOLEAN resultIsIdeal)
2073{
2074  // first check for special case h1:(0)
2075  if (idIs0(h2))
2076  {
2077    ideal res;
2078    if (resultIsIdeal)
2079    {
2080      res = idInit(1,1);
2081      res->m[0] = pOne();
2082    }
2083    else
2084      res = idFreeModule(h1->rank);
2085    return res;
2086  }
2087  BITSET old_test=test;
2088  poly     p,q = NULL;
2089  int i,l,ll,k,kkk,kmax;
2090  BOOLEAN  addOnlyOne=TRUE;
2091  tHomog   hom=isNotHomog;
2092  intvec * weights1;
2093
2094  ideal s_h4 = idInitializeQuot (h1,h2,h1IsStb,&addOnlyOne,&kmax);
2095
2096  hom = (tHomog)idHomModule(s_h4,currQuotient,&weights1);
2097
2098  ring orig_ring=currRing;
2099  ring syz_ring=rCurrRingAssure_SyzComp();
2100  rSetSyzComp(kmax-1);
2101  if (orig_ring!=syz_ring)
2102  //  s_h4 = idrMoveR_NoSort(s_h4,orig_ring);
2103    s_h4 = idrMoveR(s_h4,orig_ring);
2104  idTest(s_h4);
2105  ideal s_h3;
2106  if (addOnlyOne)
2107  {
2108    s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,kmax-1,IDELEMS(s_h4)-1);
2109  }
2110  else
2111  {
2112    s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,kmax-1);
2113  }
2114  idTest(s_h3);
2115  if (weights1!=NULL) delete weights1;
2116  idDelete(&s_h4);
2117
2118
2119  for (i=0;i<IDELEMS(s_h3);i++)
2120  {
2121    if ((s_h3->m[i]!=NULL) && (pGetComp(s_h3->m[i])>=kmax))
2122    {
2123      if (resultIsIdeal)
2124        pShift(&s_h3->m[i],-kmax);
2125      else
2126        pShift(&s_h3->m[i],-kmax+1);
2127    }
2128    else
2129      pDelete(&s_h3->m[i]);
2130  }
2131  if (resultIsIdeal)
2132    s_h3->rank = 1;
2133  else
2134    s_h3->rank = h1->rank;
2135  if(syz_ring!=orig_ring)
2136  {
2137//    pDelete(&q);
2138    rChangeCurrRing(orig_ring);
2139    s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
2140    rKill(syz_ring);
2141  }
2142  idSkipZeroes(s_h3);
2143  test = old_test;
2144  idTest(s_h3);
2145  return s_h3;
2146}
2147
2148/*2
2149*computes recursively all monomials of a certain degree
2150*in every step the actvar-th entry in the exponential
2151*vector is incremented and the other variables are
2152*computed by recursive calls of makemonoms
2153*if the last variable is reached, the difference to the
2154*degree is computed directly
2155*vars is the number variables
2156*actvar is the actual variable to handle
2157*deg is the degree of the monomials to compute
2158*monomdeg is the actual degree of the monomial in consideration
2159*/
2160static void makemonoms(int vars,int actvar,int deg,int monomdeg)
2161{
2162  poly p;
2163  int i=0;
2164
2165  if ((idpowerpoint == 0) && (actvar ==1))
2166  {
2167    idpower[idpowerpoint] = pOne();
2168    monomdeg = 0;
2169  }
2170  while (i<=deg)
2171  {
2172    if (deg == monomdeg)
2173    {
2174      pSetm(idpower[idpowerpoint]);
2175      idpowerpoint++;
2176      return;
2177    }
2178    if (actvar == vars)
2179    {
2180      pSetExp(idpower[idpowerpoint],actvar,deg-monomdeg);
2181      pSetm(idpower[idpowerpoint]);
2182      pTest(idpower[idpowerpoint]);
2183      idpowerpoint++;
2184      return;
2185    }
2186    else
2187    {
2188      p = pCopy(idpower[idpowerpoint]);
2189      makemonoms(vars,actvar+1,deg,monomdeg);
2190      idpower[idpowerpoint] = p;
2191    }
2192    monomdeg++;
2193    pSetExp(idpower[idpowerpoint],actvar,pGetExp(idpower[idpowerpoint],actvar)+1);
2194    pSetm(idpower[idpowerpoint]);
2195    pTest(idpower[idpowerpoint]);
2196    i++;
2197  }
2198}
2199
2200/*2
2201*returns the deg-th power of the maximal ideal of 0
2202*/
2203ideal idMaxIdeal(int deg)
2204{
2205  if (deg < 0)
2206  {
2207    WarnS("maxideal: power must be non-negative");
2208  }
2209  if (deg < 1)
2210  {
2211    ideal I=idInit(1,1);
2212    I->m[0]=pOne();
2213    return I;
2214  }
2215  if (deg == 1)
2216  {
2217    return idMaxIdeal();
2218  }
2219
2220  int vars = currRing->N;
2221  int i = binom(vars+deg-1,deg);
2222  if (i<=0) return idInit(1,1);
2223  ideal id=idInit(i,1);
2224  idpower = id->m;
2225  idpowerpoint = 0;
2226  makemonoms(vars,1,deg,0);
2227  idpower = NULL;
2228  idpowerpoint = 0;
2229  return id;
2230}
2231
2232/*2
2233*computes recursively all generators of a certain degree
2234*of the ideal "givenideal"
2235*elms is the number elements in the given ideal
2236*actelm is the actual element to handle
2237*deg is the degree of the power to compute
2238*gendeg is the actual degree of the generator in consideration
2239*/
2240static void makepotence(int elms,int actelm,int deg,int gendeg)
2241{
2242  poly p;
2243  int i=0;
2244
2245  if ((idpowerpoint == 0) && (actelm ==1))
2246  {
2247    idpower[idpowerpoint] = pOne();
2248    gendeg = 0;
2249  }
2250  while (i<=deg)
2251  {
2252    if (deg == gendeg)
2253    {
2254      idpowerpoint++;
2255      return;
2256    }
2257    if (actelm == elms)
2258    {
2259      p=pPower(pCopy(givenideal[actelm-1]),deg-gendeg);
2260      idpower[idpowerpoint]=pMult(idpower[idpowerpoint],p);
2261      idpowerpoint++;
2262      return;
2263    }
2264    else
2265    {
2266      p = pCopy(idpower[idpowerpoint]);
2267      makepotence(elms,actelm+1,deg,gendeg);
2268      idpower[idpowerpoint] = p;
2269    }
2270    gendeg++;
2271    idpower[idpowerpoint]=pMult(idpower[idpowerpoint],pCopy(givenideal[actelm-1]));
2272    i++;
2273  }
2274}
2275
2276/*2
2277*returns the deg-th power of the ideal gid
2278*/
2279//ideal idPower(ideal gid,int deg)
2280//{
2281//  int i;
2282//  ideal id;
2283//
2284//  if (deg < 1) deg = 1;
2285//  i = binom(IDELEMS(gid)+deg-1,deg);
2286//  id=idInit(i,1);
2287//  idpower = id->m;
2288//  givenideal = gid->m;
2289//  idpowerpoint = 0;
2290//  makepotence(IDELEMS(gid),1,deg,0);
2291//  idpower = NULL;
2292//  givenideal = NULL;
2293//  idpowerpoint = 0;
2294//  return id;
2295//}
2296static void idNextPotence(ideal given, ideal result,
2297  int begin, int end, int deg, int restdeg, poly ap)
2298{
2299  poly p;
2300  int i;
2301
2302  p = pPower(pCopy(given->m[begin]),restdeg);
2303  i = result->nrows;
2304  result->m[i] = pMult(pCopy(ap),p);
2305//PrintS(".");
2306  (result->nrows)++;
2307  if (result->nrows >= IDELEMS(result))
2308  {
2309    pEnlargeSet(&(result->m),IDELEMS(result),16);
2310    IDELEMS(result) += 16;
2311  }
2312  if (begin == end) return;
2313  for (i=restdeg-1;i>0;i--)
2314  {
2315    p = pPower(pCopy(given->m[begin]),i);
2316    p = pMult(pCopy(ap),p);
2317    idNextPotence(given, result, begin+1, end, deg, restdeg-i, p);
2318    pDelete(&p);
2319  }
2320  idNextPotence(given, result, begin+1, end, deg, restdeg, ap);
2321}
2322
2323ideal idPower(ideal given,int exp)
2324{
2325  ideal result,temp;
2326  poly p1;
2327  int i;
2328
2329  if (idIs0(given)) return idInit(1,1);
2330  temp = idCopy(given);
2331  idSkipZeroes(temp);
2332  i = binom(IDELEMS(temp)+exp-1,exp);
2333  result = idInit(i,1);
2334  result->nrows = 0;
2335//Print("ideal contains %d elements\n",i);
2336  p1=pOne();
2337  idNextPotence(temp,result,0,IDELEMS(temp)-1,exp,exp,p1);
2338  pDelete(&p1);
2339  idDelete(&temp);
2340  result->nrows = 1;
2341  idSkipZeroes(result);
2342  idDelEquals(result);
2343  return result;
2344}
2345
2346/*2
2347* eliminate delVar (product of vars) in h1
2348*/
2349ideal idElimination (ideal h1,poly delVar,intvec *hilb)
2350{
2351  int    i,j=0,k,l;
2352  ideal  h,hh, h3;
2353  int    *ord,*block0,*block1;
2354  int    ordersize=2;
2355  int    **wv;
2356  tHomog hom;
2357  intvec * w;
2358  sip_sring tmpR;
2359  ring origR = currRing;
2360
2361  if (delVar==NULL)
2362  {
2363    return idCopy(h1);
2364  }
2365  if (currQuotient!=NULL)
2366  {
2367    WerrorS("cannot eliminate in a qring");
2368    return idCopy(h1);
2369  }
2370  if (idIs0(h1)) return idInit(1,h1->rank);
2371#ifdef HAVE_PLURAL
2372  if (rIsPluralRing(currRing))
2373    /* in the NC case, we have to check the admissibility of */
2374    /* the subalgebra to be intersected with */
2375  {
2376    if (currRing->nc->type!=nc_skew)
2377    {
2378      if (!nc_CheckSubalgebra(delVar,currRing))
2379      {
2380        WerrorS("no elimination is possible: subalgebra is not admissible");
2381        return idCopy(h1);
2382      }
2383    }
2384  }
2385#endif
2386  hom=(tHomog)idHomModule(h1,NULL,&w); //sets w to weight vector or NULL
2387  h3=idInit(16,h1->rank);
2388  for (k=0;; k++)
2389  {
2390    if (currRing->order[k]!=0) ordersize++;
2391    else break;
2392  }
2393  ord=(int*)omAlloc0(ordersize*sizeof(int));
2394  block0=(int*)omAlloc(ordersize*sizeof(int));
2395  block1=(int*)omAlloc(ordersize*sizeof(int));
2396  for (k=0;; k++)
2397  {
2398    if (currRing->order[k]!=0)
2399    {
2400      block0[k+1] = currRing->block0[k];
2401      block1[k+1] = currRing->block1[k];
2402      ord[k+1] = currRing->order[k];
2403    }
2404    else
2405      break;
2406  }
2407  block0[0] = 1;
2408  block1[0] = pVariables;
2409  wv=(int**) omAlloc0(ordersize*sizeof(int**));
2410  memcpy4(wv+1,currRing->wvhdl,(ordersize-1)*sizeof(int**));
2411  wv[0]=(int*)omAlloc((pVariables+1)*sizeof(int));
2412  memset(wv[0],0,(pVariables+1)*sizeof(int));
2413  for (j=0;j<pVariables;j++)
2414    if (pGetExp(delVar,j+1)!=0) wv[0][j]=1;
2415  // use this special ordering: like ringorder_a, except that pFDeg, pWeights
2416  // ignore it
2417  ord[0] = ringorder_aa;
2418
2419  // fill in tmp ring to get back the data later on
2420  tmpR  = *origR;
2421  tmpR.order = ord;
2422  tmpR.block0 = block0;
2423  tmpR.block1 = block1;
2424  tmpR.wvhdl = wv;
2425  rComplete(&tmpR, 1);
2426
2427  // change into the new ring
2428  //pChangeRing(pVariables,currRing->OrdSgn,ord,block0,block1,wv);
2429  rChangeCurrRing(&tmpR);
2430  h = idInit(IDELEMS(h1),h1->rank);
2431  // fetch data from the old ring
2432  for (k=0;k<IDELEMS(h1);k++) h->m[k] = prCopyR( h1->m[k], origR);
2433  // compute kStd
2434  hh = kStd(h,NULL,hom,&w,hilb);
2435  idDelete(&h);
2436
2437  // go back to the original ring
2438  rChangeCurrRing(origR);
2439  i = IDELEMS(hh)-1;
2440  while ((i >= 0) && (hh->m[i] == NULL)) i--;
2441  j = -1;
2442  // fetch data from temp ring
2443  for (k=0; k<=i; k++)
2444  {
2445    l=pVariables;
2446    while ((l>0) && (p_GetExp( hh->m[k],l,&tmpR)*pGetExp(delVar,l)==0)) l--;
2447    if (l==0)
2448    {
2449      j++;
2450      if (j >= IDELEMS(h3))
2451      {
2452        pEnlargeSet(&(h3->m),IDELEMS(h3),16);
2453        IDELEMS(h3) += 16;
2454      }
2455      h3->m[j] = prCopyR( hh->m[k], &tmpR);
2456    }
2457  }
2458  id_Delete(&hh, &tmpR);
2459  idSkipZeroes(h3);
2460  omFree((ADDRESS)wv[0]);
2461  omFreeSize((ADDRESS)wv,ordersize*sizeof(int**));
2462  omFreeSize((ADDRESS)ord,ordersize*sizeof(int));
2463  omFreeSize((ADDRESS)block0,ordersize*sizeof(int));
2464  omFreeSize((ADDRESS)block1,ordersize*sizeof(int));
2465  rUnComplete(&tmpR);
2466  if (w!=NULL)
2467    delete w;
2468  return h3;
2469}
2470
2471#ifdef WITH_OLD_MINOR
2472/*2
2473* compute all ar-minors of the matrix a
2474*/
2475ideal idMinors(matrix a, int ar, ideal R)
2476{
2477  int     i,j,k,size;
2478  int *rowchoise,*colchoise;
2479  BOOLEAN rowch,colch;
2480  ideal result;
2481  matrix tmp;
2482  poly p,q;
2483
2484  i = binom(a->rows(),ar);
2485  j = binom(a->cols(),ar);
2486
2487  rowchoise=(int *)omAlloc(ar*sizeof(int));
2488  colchoise=(int *)omAlloc(ar*sizeof(int));
2489  if ((i>512) || (j>512) || (i*j >512)) size=512;
2490  else size=i*j;
2491  result=idInit(size,1);
2492  tmp=mpNew(ar,ar);
2493  k = 0; /* the index in result*/
2494  idInitChoise(ar,1,a->rows(),&rowch,rowchoise);
2495  while (!rowch)
2496  {
2497    idInitChoise(ar,1,a->cols(),&colch,colchoise);
2498    while (!colch)
2499    {
2500      for (i=1; i<=ar; i++)
2501      {
2502        for (j=1; j<=ar; j++)
2503        {
2504          MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
2505        }
2506      }
2507      p = mpDetBareiss(tmp);
2508      if (p!=NULL)
2509      {
2510        if (R!=NULL)
2511        {
2512          q = p;
2513          p = kNF(R,currQuotient,q);
2514          pDelete(&q);
2515        }
2516        if (p!=NULL)
2517        {
2518          if (k>=size)
2519          {
2520            pEnlargeSet(&result->m,size,32);
2521            size += 32;
2522          }
2523          result->m[k] = p;
2524          k++;
2525        }
2526      }
2527      idGetNextChoise(ar,a->cols(),&colch,colchoise);
2528    }
2529    idGetNextChoise(ar,a->rows(),&rowch,rowchoise);
2530  }
2531  /*delete the matrix tmp*/
2532  for (i=1; i<=ar; i++)
2533  {
2534    for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
2535  }
2536  idDelete((ideal*)&tmp);
2537  if (k==0)
2538  {
2539    k=1;
2540    result->m[0]=NULL;
2541  }
2542  omFreeSize((ADDRESS)rowchoise,ar*sizeof(int));
2543  omFreeSize((ADDRESS)colchoise,ar*sizeof(int));
2544  pEnlargeSet(&result->m,size,k-size);
2545  IDELEMS(result) = k;
2546  return (result);
2547}
2548#else
2549/*2
2550* compute all ar-minors of the matrix a
2551* the caller of mpRecMin
2552* the elements of the result are not in R (if R!=NULL)
2553*/
2554ideal idMinors(matrix a, int ar, ideal R)
2555{
2556  int elems=0;
2557  int r=a->nrows,c=a->ncols;
2558  int i;
2559  matrix b;
2560  ideal result,h;
2561  ring origR;
2562  sip_sring tmpR;
2563  Exponent_t bound;
2564
2565  if((ar<=0) || (ar>r) || (ar>c))
2566  {
2567    Werror("%d-th minor, matrix is %dx%d",ar,r,c);
2568    return NULL;
2569  }
2570  h = idMatrix2Module(mpCopy(a));
2571  bound = smExpBound(h,c,r,ar);
2572  idDelete(&h);
2573  smRingChange(&origR,tmpR,bound);
2574  b = mpNew(r,c);
2575  for (i=r*c-1;i>=0;i--)
2576  {
2577    if (a->m[i])
2578      b->m[i] = prCopyR(a->m[i],origR);
2579  }
2580  if (R) R = idrCopyR(R,origR);
2581  result=idInit(32,1);
2582  if(ar>1) mpRecMin(ar-1,result,elems,b,r,c,NULL,R);
2583  else mpMinorToResult(result,elems,b,r,c,R);
2584  idDelete((ideal *)&b);
2585  if (R) idDelete(&R);
2586  idSkipZeroes(result);
2587  rChangeCurrRing(origR);
2588  result = idrMoveR(result,&tmpR);
2589  smRingClean(origR,tmpR);
2590  idTest(result);
2591  return result;
2592}
2593#endif
2594
2595/*2
2596*skips all zeroes and double elements, searches also for units
2597*/
2598ideal idCompactify(ideal id)
2599{
2600  int i,j;
2601  BOOLEAN b=FALSE;
2602
2603  i = IDELEMS(id)-1;
2604  while ((! b) && (i>=0))
2605  {
2606    b=pIsUnit(id->m[i]);
2607    i--;
2608  }
2609  if (b)
2610  {
2611    ideal result=idInit(1,id->rank);
2612    result->m[0]=pOne();
2613    return result;
2614  }
2615  else
2616  {
2617    ideal result=idCopy(id);
2618    //if (IDELEMS(result) < 1)
2619    {
2620      for (i=1;i<IDELEMS(result);i++)
2621      {
2622        if (result->m[i]!=NULL)
2623        {
2624          for (j=0;j<i;j++)
2625          {
2626            if ((result->m[j]!=NULL)
2627            && (pComparePolys(result->m[i],result->m[j])))
2628            {
2629              pDelete(&(result->m[j]));
2630            }
2631          }
2632        }
2633      }
2634    }
2635//    else
2636//    {
2637//      intvec *v=idSort(result,TRUE);
2638//      for(i=0;i<IDELEMS(result);i++)
2639//        (*v)[i]--;
2640//      for (i=0;i<IDELEMS(result)-1;i++)
2641//      {
2642//        if (((*v)[i]>=0)
2643//        && (result->m[(*v)[i]]!=NULL))
2644//        {
2645//          j=i+1;
2646//          while ((j<IDELEMS(result))
2647//          && ((*v)[j]>=0)
2648//          && (result->m[(*v)[j]]!=NULL)
2649//          && (pComparePolys(result->m[(*v)[i]],result->m[(*v)[j]])))
2650//          {
2651//            pDelete(&(result->m[(*v)[j]]));
2652//            j++;
2653//          }
2654//        }
2655//      }
2656//      delete v;
2657//    }
2658    idSkipZeroes(result);
2659    return result;
2660  }
2661}
2662
2663/*2
2664*returns TRUE if id1 is a submodule of id2
2665*/
2666BOOLEAN idIsSubModule(ideal id1,ideal id2)
2667{
2668  int i;
2669  poly p;
2670
2671  if (idIs0(id1)) return TRUE;
2672  for (i=0;i<IDELEMS(id1);i++)
2673  {
2674    if (id1->m[i] != NULL)
2675    {
2676      p = kNF(id2,currQuotient,id1->m[i]);
2677      if (p != NULL)
2678      {
2679        pDelete(&p);
2680        return FALSE;
2681      }
2682    }
2683  }
2684  return TRUE;
2685}
2686
2687/*2
2688* returns the ideals of initial terms
2689*/
2690ideal idHead(ideal h)
2691{
2692  ideal m = idInit(IDELEMS(h),h->rank);
2693  int i;
2694
2695  for (i=IDELEMS(h)-1;i>=0; i--)
2696  {
2697    if (h->m[i]!=NULL) m->m[i]=pHead(h->m[i]);
2698  }
2699  return m;
2700}
2701
2702ideal idHomogen(ideal h, int varnum)
2703{
2704  ideal m = idInit(IDELEMS(h),h->rank);
2705  int i;
2706
2707  for (i=IDELEMS(h)-1;i>=0; i--)
2708  {
2709    m->m[i]=pHomogen(h->m[i],varnum);
2710  }
2711  return m;
2712}
2713
2714/*------------------type conversions----------------*/
2715ideal idVec2Ideal(poly vec)
2716{
2717   ideal result=idInit(1,1);
2718   omFree((ADDRESS)result->m);
2719   result->m=NULL; // remove later
2720   pVec2Polys(vec, &(result->m), &(IDELEMS(result)));
2721   return result;
2722}
2723
2724// #define NEW_STUFF
2725#ifndef NEW_STUFF
2726// converts mat to module, destroys mat
2727ideal idMatrix2Module(matrix mat)
2728{
2729  int mc=MATCOLS(mat);
2730  int mr=MATROWS(mat);
2731  ideal result = idInit(si_max(mc,1),si_max(mr,1));
2732  int i,j;
2733  poly h;
2734
2735  for(j=0;j<mc /*MATCOLS(mat)*/;j++) /* j is also index in result->m */
2736  {
2737    for (i=1;i<=mr /*MATROWS(mat)*/;i++)
2738    {
2739      h = MATELEM(mat,i,j+1);
2740      if (h!=NULL)
2741      {
2742        MATELEM(mat,i,j+1)=NULL;
2743        pSetCompP(h,i);
2744        result->m[j] = pAdd(result->m[j],h);
2745      }
2746    }
2747  }
2748  // obachman: need to clean this up
2749  idDelete((ideal*) &mat);
2750  return result;
2751}
2752#else
2753
2754#include "sbuckets.h"
2755
2756// converts mat to module, destroys mat
2757ideal idMatrix2Module(matrix mat)
2758{
2759  int mc=MATCOLS(mat);
2760  int mr=MATROWS(mat);
2761  ideal result = idInit(si_max(mc,1),si_max(mr,1));
2762  int i,j, l;
2763  poly h;
2764  poly p;
2765  sBucket_pt bucket = sBucketInit(currRing);
2766
2767  for(j=0;j<mc /*MATCOLS(mat)*/;j++) /* j is also index in result->m */
2768  {
2769    for (i=1;i<=mr /*MATROWS(mat)*/;i++)
2770    {
2771      h = MATELEM(mat,i,j+1);
2772      if (h!=NULL)
2773      {
2774        MATELEM(mat,i,j+1)=NULL;
2775        p_SetCompP(h,i, currRing, &l);
2776        sBucket_Merge_p(bucket, h, l);
2777      }
2778    }
2779    sBucketClearMerge(bucket, &(result->m[j]), &l);
2780  }
2781
2782  // obachman: need to clean this up
2783  idDelete((ideal*) &mat);
2784  return result;
2785}
2786#endif
2787
2788/*2
2789* converts a module into a matrix, destroyes the input
2790*/
2791matrix idModule2Matrix(ideal mod)
2792{
2793  matrix result = mpNew(mod->rank,IDELEMS(mod));
2794  int i,cp;
2795  poly p,h;
2796
2797  for(i=0;i<IDELEMS(mod);i++)
2798  {
2799    p=mod->m[i];
2800    mod->m[i]=NULL;
2801    while (p!=NULL)
2802    {
2803      h=p;
2804      pIter(p);
2805      pNext(h)=NULL;
2806//      cp = si_max(1,pGetComp(h));     // if used for ideals too
2807      cp = pGetComp(h);
2808      pSetComp(h,0);
2809      pSetmComp(h);
2810#ifdef TEST
2811      if (cp>mod->rank)
2812      {
2813        Print("## inv. rank %d -> %d\n",mod->rank,cp);
2814        int k,l,o=mod->rank;
2815        mod->rank=cp;
2816        matrix d=mpNew(mod->rank,IDELEMS(mod));
2817        for (l=1; l<=o; l++)
2818        {
2819          for (k=1; k<=IDELEMS(mod); k++)
2820          {
2821            MATELEM(d,l,k)=MATELEM(result,l,k);
2822            MATELEM(result,l,k)=NULL;
2823          }
2824        }
2825        idDelete((ideal *)&result);
2826        result=d;
2827      }
2828#endif
2829      MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
2830    }
2831  }
2832  // obachman 10/99: added the following line, otherwise memory leack!
2833  idDelete(&mod);
2834  return result;
2835}
2836
2837matrix idModule2formatedMatrix(ideal mod,int rows, int cols)
2838{
2839  matrix result = mpNew(rows,cols);
2840  int i,cp,r=idRankFreeModule(mod),c=IDELEMS(mod);
2841  poly p,h;
2842
2843  if (r>rows) r = rows;
2844  if (c>cols) c = cols;
2845  for(i=0;i<c;i++)
2846  {
2847    p=mod->m[i];
2848    mod->m[i]=NULL;
2849    while (p!=NULL)
2850    {
2851      h=p;
2852      pIter(p);
2853      pNext(h)=NULL;
2854      cp = pGetComp(h);
2855      if (cp<=r)
2856      {
2857        pSetComp(h,0);
2858        pSetmComp(h);
2859        MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
2860      }
2861      else
2862        pDelete(&h);
2863    }
2864  }
2865  idDelete(&mod);
2866  return result;
2867}
2868
2869/*2
2870* substitute the n-th variable by the monomial e in id
2871* destroy id
2872*/
2873ideal  idSubst(ideal id, int n, poly e)
2874{
2875  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
2876  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
2877
2878  res->rank = id->rank;
2879  for(k--;k>=0;k--)
2880  {
2881    res->m[k]=pSubst(id->m[k],n,e);
2882    id->m[k]=NULL;
2883  }
2884  idDelete(&id);
2885  return res;
2886}
2887
2888BOOLEAN idHomModule(ideal m, ideal Q, intvec **w)
2889{
2890  if (w!=NULL) *w=NULL;
2891  if ((Q!=NULL) && (!idHomIdeal(Q,NULL))) return FALSE;
2892  if (idIs0(m)) return TRUE;
2893
2894  int i,j,cmax=2,order=0,ord,* diff,* iscom,diffmin=32000;
2895  poly p=NULL;
2896  int length=IDELEMS(m);
2897  polyset P=m->m;
2898  polyset F=(polyset)omAlloc(length*sizeof(poly));
2899  for (i=length-1;i>=0;i--)
2900  {
2901    p=F[i]=P[i];
2902    cmax=si_max(cmax,(int)pMaxComp(p)+1);
2903  }
2904  diff = (int *)omAlloc0(cmax*sizeof(int));
2905  if (w!=NULL) *w=new intvec(cmax-1);
2906  iscom = (int *)omAlloc0(cmax*sizeof(int));
2907  i=0;
2908  while (i<=length)
2909  {
2910    if (i<length)
2911    {
2912      p=F[i];
2913      while ((p!=NULL) && (!iscom[pGetComp(p)])) pIter(p);
2914    }
2915    if ((p==NULL) && (i<length))
2916    {
2917      i++;
2918    }
2919    else
2920    {
2921      if (p==NULL)
2922      {
2923        i=0;
2924        while ((i<length) && (F[i]==NULL)) i++;
2925        if (i>=length) break;
2926        p = F[i];
2927      }
2928      if (pLexOrder)
2929        order=pTotaldegree(p);
2930      else
2931      //  order = p->order;
2932        order = pFDeg(p,currRing);
2933      order += diff[pGetComp(p)];
2934      p = F[i];
2935//Print("Actual p=F[%d]: ",i);pWrite(p);
2936      F[i] = NULL;
2937      i=0;
2938    }
2939    while (p!=NULL)
2940    {
2941      //if (pLexOrder)
2942      //  ord=pTotaldegree(p);
2943      //else
2944      //  ord = p->order;
2945      ord = pFDeg(p,currRing);
2946      if (!iscom[pGetComp(p)])
2947      {
2948        diff[pGetComp(p)] = order-ord;
2949        iscom[pGetComp(p)] = 1;
2950/*
2951*PrintS("new diff: ");
2952*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
2953*PrintLn();
2954*PrintS("new iscom: ");
2955*for (j=0;j<cmax;j++) Print("%d ",iscom[j]);
2956*PrintLn();
2957*Print("new set %d, order %d, ord %d, diff %d\n",pGetComp(p),order,ord,diff[pGetComp(p)]);
2958*/
2959      }
2960      else
2961      {
2962/*
2963*PrintS("new diff: ");
2964*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
2965*PrintLn();
2966*Print("order %d, ord %d, diff %d\n",order,ord,diff[pGetComp(p)]);
2967*/
2968        if (order != ord+diff[pGetComp(p)])
2969        {
2970          omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
2971          omFreeSize((ADDRESS) diff,cmax*sizeof(int));
2972          omFreeSize((ADDRESS) F,length*sizeof(poly));
2973          delete *w;*w=NULL;
2974          return FALSE;
2975        }
2976      }
2977      pIter(p);
2978    }
2979  }
2980  omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
2981  omFreeSize((ADDRESS) F,length*sizeof(poly));
2982  for (i=1;i<cmax;i++) (**w)[i-1]=diff[i];
2983  for (i=1;i<cmax;i++)
2984  {
2985    if (diff[i]<diffmin) diffmin=diff[i];
2986  }
2987  if (w!=NULL)
2988  {
2989    for (i=1;i<cmax;i++)
2990    {
2991      (**w)[i-1]=diff[i]-diffmin;
2992    }
2993  }
2994  omFreeSize((ADDRESS) diff,cmax*sizeof(int));
2995  return TRUE;
2996}
2997
2998BOOLEAN idTestHomModule(ideal m, ideal Q, intvec *w)
2999{
3000  if ((Q!=NULL) && (!idHomIdeal(Q,NULL)))  { PrintS(" Q not hom\n"); return FALSE;}
3001  if (idIs0(m)) return TRUE;
3002
3003  int cmax=-1;
3004  int i;
3005  poly p=NULL;
3006  int length=IDELEMS(m);
3007  polyset P=m->m;
3008  for (i=length-1;i>=0;i--)
3009  {
3010    p=P[i];
3011    if (p!=NULL) cmax=si_max(cmax,(int)pMaxComp(p)+1);
3012  }
3013  if (w->length()+1 < cmax)
3014  { 
3015    // Print("length: %d - %d \n", w->length(),cmax);
3016    return FALSE;
3017  }
3018  pSetModDeg(w);
3019  for (i=length-1;i>=0;i--)
3020  {
3021    p=P[i];
3022    poly q=p;
3023    if (p!=NULL)
3024    {
3025      int d=pFDeg(p,currRing);
3026      loop
3027      {
3028        pIter(p);
3029        if (p==NULL) break;
3030        if (d!=pFDeg(p,currRing)) 
3031        {
3032          //pWrite(q); wrp(p); Print(" -> %d - %d\n",d,pFDeg(p,currRing));
3033          return FALSE; 
3034        }
3035      }
3036    }
3037  }
3038  pSetModDeg(NULL);
3039  return TRUE;
3040}
3041
3042ideal idJet(ideal i,int d)
3043{
3044  ideal r=idInit((i->nrows)*(i->ncols),i->rank);
3045  r->nrows = i-> nrows;
3046  r->ncols = i-> ncols;
3047  //r->rank = i-> rank;
3048  int k;
3049  for(k=(i->nrows)*(i->ncols)-1;k>=0; k--)
3050  {
3051    r->m[k]=ppJet(i->m[k],d);
3052  }
3053  return r;
3054}
3055
3056ideal idJetW(ideal i,int d, intvec * iv)
3057{
3058  ideal r=idInit(IDELEMS(i),i->rank);
3059  if (ecartWeights!=NULL)
3060  {
3061    WerrorS("cannot compute weighted jets now");
3062  }
3063  else
3064  {
3065    short *w=iv2array(iv);
3066    int k;
3067    for(k=0; k<IDELEMS(i); k++)
3068    {
3069      r->m[k]=ppJetW(i->m[k],d,w);
3070    }
3071    omFreeSize((ADDRESS)w,(pVariables+1)*sizeof(short));
3072  }
3073  return r;
3074}
3075
3076int idMinDegW(ideal M,intvec *w)
3077{
3078  int d=-1;
3079  for(int i=0;i<IDELEMS(M);i++)
3080  {
3081    int d0=pMinDeg(M->m[i],w);
3082    if(-1<d0&&(d0<d||d==-1))
3083      d=d0;
3084  }
3085  return d;
3086}
3087
3088ideal idSeries(int n,ideal M,matrix U,intvec *w)
3089{
3090  for(int i=IDELEMS(M)-1;i>=0;i--)
3091  {
3092    if(U==NULL)
3093      M->m[i]=pSeries(n,M->m[i],NULL,w);
3094    else
3095    {
3096      M->m[i]=pSeries(n,M->m[i],MATELEM(U,i+1,i+1),w);
3097      MATELEM(U,i+1,i+1)=NULL;
3098    }
3099  }
3100  if(U!=NULL)
3101    idDelete((ideal*)&U);
3102  return M;
3103}
3104
3105matrix idDiff(matrix i, int k)
3106{
3107  int e=MATCOLS(i)*MATROWS(i);
3108  matrix r=mpNew(MATROWS(i),MATCOLS(i));
3109  int j;
3110  for(j=0; j<e; j++)
3111  {
3112    r->m[j]=pDiff(i->m[j],k);
3113  }
3114  return r;
3115}
3116
3117matrix idDiffOp(ideal I, ideal J,BOOLEAN multiply)
3118{
3119  matrix r=mpNew(IDELEMS(I),IDELEMS(J));
3120  int i,j;
3121  for(i=0; i<IDELEMS(I); i++)
3122  {
3123    for(j=0; j<IDELEMS(J); j++)
3124    {
3125      MATELEM(r,i+1,j+1)=pDiffOp(I->m[i],J->m[j],multiply);
3126    }
3127  }
3128  return r;
3129}
3130
3131/*3
3132*handles for some ideal operations the ring/syzcomp managment
3133*returns all syzygies (componentwise-)shifted by -syzcomp
3134*or -syzcomp-1 (in case of ideals as input)
3135static ideal idHandleIdealOp(ideal arg,int syzcomp,int isIdeal=FALSE)
3136{
3137  ring orig_ring=currRing;
3138  ring syz_ring=rCurrRingAssure_SyzComp();
3139  rSetSyzComp(length);
3140
3141  ideal s_temp;
3142  if (orig_ring!=syz_ring)
3143    s_temp=idrMoveR_NoSort(arg,orig_ring);
3144  else
3145    s_temp=arg;
3146
3147  ideal s_temp1 = kStd(s_temp,currQuotient,testHomog,&w,NULL,length);
3148  if (w!=NULL) delete w;
3149
3150  if (syz_ring!=orig_ring)
3151  {
3152    idDelete(&s_temp);
3153    rChangeCurrRing(orig_ring);
3154  }
3155
3156  idDelete(&temp);
3157  ideal temp1=idRingCopy(s_temp1,syz_ring);
3158
3159  if (syz_ring!=orig_ring)
3160  {
3161    rChangeCurrRing(syz_ring);
3162    idDelete(&s_temp1);
3163    rChangeCurrRing(orig_ring);
3164    rKill(syz_ring);
3165  }
3166
3167  for (i=0;i<IDELEMS(temp1);i++)
3168  {
3169    if ((temp1->m[i]!=NULL)
3170    && (pGetComp(temp1->m[i])<=length))
3171    {
3172      pDelete(&(temp1->m[i]));
3173    }
3174    else
3175    {
3176      pShift(&(temp1->m[i]),-length);
3177    }
3178  }
3179  temp1->rank = rk;
3180  idSkipZeroes(temp1);
3181
3182  return temp1;
3183}
3184*/
3185/*2
3186* represents (h1+h2)/h2=h1/(h1 intersect h2)
3187*/
3188ideal idModulo (ideal h2,ideal h1)
3189{
3190  int i,j,k,rk,flength=0,slength,length;
3191  intvec * w;
3192  poly p,q;
3193
3194  if (idIs0(h2))
3195    return idFreeModule(si_max(1,h2->ncols));
3196  if (!idIs0(h1))
3197    flength = idRankFreeModule(h1);
3198  slength = idRankFreeModule(h2);
3199  length  = si_max(flength,slength);
3200  if (length==0)
3201  {
3202    length = 1;
3203  }
3204  ideal temp = idInit(IDELEMS(h2),length+IDELEMS(h2));
3205  for (i=0;i<IDELEMS(h2);i++)
3206  {
3207    temp->m[i] = pCopy(h2->m[i]);
3208    q = pOne();
3209    pSetComp(q,i+1+length);
3210    pSetmComp(q);
3211    if(temp->m[i]!=NULL)
3212    {
3213      if (slength==0) pShift(&(temp->m[i]),1);
3214      p = temp->m[i];
3215      while (pNext(p)!=NULL) pIter(p);
3216      pNext(p) = q;
3217    }
3218    else
3219      temp->m[i]=q;
3220  }
3221  rk = k = IDELEMS(h2);
3222  if (!idIs0(h1))
3223  {
3224    pEnlargeSet(&(temp->m),IDELEMS(temp),IDELEMS(h1));
3225    IDELEMS(temp) += IDELEMS(h1);
3226    for (i=0;i<IDELEMS(h1);i++)
3227    {
3228      if (h1->m[i]!=NULL)
3229      {
3230        temp->m[k] = pCopy(h1->m[i]);
3231        if (flength==0) pShift(&(temp->m[k]),1);
3232        k++;
3233      }
3234    }
3235  }
3236
3237  ring orig_ring=currRing;
3238  ring syz_ring=rCurrRingAssure_SyzComp();
3239  rSetSyzComp(length);
3240  ideal s_temp;
3241
3242  if (syz_ring != orig_ring)
3243  {
3244    s_temp = idrMoveR_NoSort(temp, orig_ring);
3245  }
3246  else
3247  {
3248    s_temp = temp;
3249  }
3250
3251  idTest(s_temp);
3252  ideal s_temp1 = kStd(s_temp,currQuotient,testHomog,&w,NULL,length);
3253  if (w!=NULL) delete w;
3254
3255  for (i=0;i<IDELEMS(s_temp1);i++)
3256  {
3257    if ((s_temp1->m[i]!=NULL)
3258    && (pGetComp(s_temp1->m[i])<=length))
3259    {
3260      pDelete(&(s_temp1->m[i]));
3261    }
3262    else
3263    {
3264      pShift(&(s_temp1->m[i]),-length);
3265    }
3266  }
3267  s_temp1->rank = rk;
3268  idSkipZeroes(s_temp1);
3269
3270  if (syz_ring!=orig_ring)
3271  {
3272    rChangeCurrRing(orig_ring);
3273    s_temp1 = idrMoveR_NoSort(s_temp1, syz_ring);
3274    rKill(syz_ring);
3275    // Hmm ... here seems to be a memory leak
3276    // However, simply deleting it causes memory trouble
3277    // idDelete(&s_temp);
3278  }
3279  else
3280  {
3281    idDelete(&temp);
3282  }
3283  idTest(s_temp1);
3284  return s_temp1;
3285}
3286
3287int idElem(ideal F)
3288{
3289  int i=0,j=0;
3290
3291  while(j<IDELEMS(F))
3292  {
3293   if ((F->m)[j]!=NULL) i++;
3294   j++;
3295  }
3296  return i;
3297}
3298
3299/*
3300*computes module-weights for liftings of homogeneous modules
3301*/
3302intvec * idMWLift(ideal mod,intvec * weights)
3303{
3304  if (idIs0(mod)) return new intvec(2);
3305  int i=IDELEMS(mod);
3306  while ((i>0) && (mod->m[i-1]==NULL)) i--;
3307  intvec *result = new intvec(i+1);
3308  while (i>0)
3309  {
3310    (*result)[i]=pFDeg(mod->m[i],currRing)+(*weights)[pGetComp(mod->m[i])];
3311  }
3312  return result;
3313}
3314
3315/*2
3316*sorts the kbase for idCoef* in a special way (lexicographically
3317*with x_max,...,x_1)
3318*/
3319ideal idCreateSpecialKbase(ideal kBase,intvec ** convert)
3320{
3321  int i;
3322  ideal result;
3323
3324  if (idIs0(kBase)) return NULL;
3325  result = idInit(IDELEMS(kBase),kBase->rank);
3326  *convert = idSort(kBase,FALSE);
3327  for (i=0;i<(*convert)->length();i++)
3328  {
3329    result->m[i] = pCopy(kBase->m[(**convert)[i]-1]);
3330  }
3331  return result;
3332}
3333
3334/*2
3335*returns the index of a given monom in the list of the special kbase
3336*/
3337int idIndexOfKBase(poly monom, ideal kbase)
3338{
3339  int j=IDELEMS(kbase);
3340
3341  while ((j>0) && (kbase->m[j-1]==NULL)) j--;
3342  if (j==0) return -1;
3343  int i=pVariables;
3344  while (i>0)
3345  {
3346    loop
3347    {
3348      if (pGetExp(monom,i)>pGetExp(kbase->m[j-1],i)) return -1;
3349      if (pGetExp(monom,i)==pGetExp(kbase->m[j-1],i)) break;
3350      j--;
3351      if (j==0) return -1;
3352    }
3353    if (i==1)
3354    {
3355      while(j>0)
3356      {
3357        if (pGetComp(monom)==pGetComp(kbase->m[j-1])) return j-1;
3358        if (pGetComp(monom)>pGetComp(kbase->m[j-1])) return -1;
3359        j--;
3360      }
3361    }
3362    i--;
3363  }
3364  return -1;
3365}
3366
3367/*2
3368*decomposes the monom in a part of coefficients described by the
3369*complement of how and a monom in variables occuring in how, the
3370*index of which in kbase is returned as integer pos (-1 if it don't
3371*exists)
3372*/
3373poly idDecompose(poly monom, poly how, ideal kbase, int * pos)
3374{
3375  int i;
3376  poly coeff=pOne(), base=pOne();
3377
3378  for (i=1;i<=pVariables;i++)
3379  {
3380    if (pGetExp(how,i)>0)
3381    {
3382      pSetExp(base,i,pGetExp(monom,i));
3383    }
3384    else
3385    {
3386      pSetExp(coeff,i,pGetExp(monom,i));
3387    }
3388  }
3389  pSetComp(base,pGetComp(monom));
3390  pSetm(base);
3391  pSetCoeff(coeff,nCopy(pGetCoeff(monom)));
3392  pSetm(coeff);
3393  *pos = idIndexOfKBase(base,kbase);
3394  if (*pos<0)
3395    pDelete(&coeff);
3396  pDelete(&base);
3397  return coeff;
3398}
3399
3400/*2
3401*returns a matrix A of coefficients with kbase*A=arg
3402*if all monomials in variables of how occur in kbase
3403*the other are deleted
3404*/
3405matrix idCoeffOfKBase(ideal arg, ideal kbase, poly how)
3406{
3407  matrix result;
3408  ideal tempKbase;
3409  poly p,q;
3410  intvec * convert;
3411  int i=IDELEMS(kbase),j=IDELEMS(arg),k,pos;
3412#if 0
3413  while ((i>0) && (kbase->m[i-1]==NULL)) i--;
3414  if (idIs0(arg))
3415    return mpNew(i,1);
3416  while ((j>0) && (arg->m[j-1]==NULL)) j--;
3417  result = mpNew(i,j);
3418#else
3419  result = mpNew(i, j);
3420  while ((j>0) && (arg->m[j-1]==NULL)) j--;
3421#endif
3422
3423  tempKbase = idCreateSpecialKbase(kbase,&convert);
3424  for (k=0;k<j;k++)
3425  {
3426    p = arg->m[k];
3427    while (p!=NULL)
3428    {
3429      q = idDecompose(p,how,tempKbase,&pos);
3430      if (pos>=0)
3431      {
3432        MATELEM(result,(*convert)[pos],k+1) =
3433            pAdd(MATELEM(result,(*convert)[pos],k+1),q);
3434      }
3435      else
3436        pDelete(&q);
3437      pIter(p);
3438    }
3439  }
3440  idDelete(&tempKbase);
3441  return result;
3442}
3443
3444/*3
3445* searches for units in the components of the module arg and
3446* returns the first one
3447*/
3448static int idReadOutUnits(ideal arg,int* comp)
3449{
3450  if (idIs0(arg)) return -1;
3451  int i=0,j,rk_arg=idRankFreeModule(arg),generator=-1;
3452  intvec * componentIsUsed =new intvec(rk_arg+1);
3453  poly p,q;
3454
3455  while ((i<IDELEMS(arg)) && (generator<0))
3456  {
3457    for (j=rk_arg;j>=0;j--)
3458      (*componentIsUsed)[j]=0;
3459    p = arg->m[i];
3460    while (p!=NULL)
3461    {
3462      j = pGetComp(p);
3463      if ((*componentIsUsed)[j]==0)
3464      {
3465        if (pLmIsConstantComp(p))
3466        {
3467          generator = i;
3468          (*componentIsUsed)[j] = 1;
3469        }
3470        else
3471        {
3472          (*componentIsUsed)[j] = -1;
3473        }
3474      }
3475      else if ((*componentIsUsed)[j]>0)
3476      {
3477        ((*componentIsUsed)[j])++;
3478      }
3479      pIter(p);
3480    }
3481    i++;
3482  }
3483  i = 0;
3484  *comp = -1;
3485  for (j=0;j<=rk_arg;j++)
3486  {
3487    if ((*componentIsUsed)[j]>0)
3488    {
3489      if ((*comp==-1) || ((*componentIsUsed)[j]<i))
3490      {
3491        *comp = j;
3492        i= (*componentIsUsed)[j];
3493      }
3494    }
3495  }
3496  return generator;
3497}
3498
3499static void idDeleteComp(ideal arg,int red_comp)
3500{
3501  int i,j;
3502  poly p;
3503
3504  for (i=IDELEMS(arg)-1;i>=0;i--)
3505  {
3506    p = arg->m[i];
3507    while (p!=NULL)
3508    {
3509      j = pGetComp(p);
3510      if (j>red_comp)
3511      {
3512        pSetComp(p,j-1);
3513        pSetm(p);
3514      }
3515      pIter(p);
3516    }
3517  }
3518  (arg->rank)--;
3519}
3520
3521/*2
3522* returns the presentation of an isomorphic, minimally
3523* embedded  module (arg represents the quotient!)
3524*/
3525ideal idMinEmbedding(ideal arg,BOOLEAN inPlace)
3526{
3527  if (idIs0(arg)) return idInit(1,arg->rank);
3528  int next_gen,next_comp;
3529  ideal res=arg;
3530
3531  if (!inPlace) res = idCopy(arg);
3532  loop
3533  {
3534    next_gen = idReadOutUnits(res,&next_comp);
3535    if (next_gen<0) break;
3536    syGaussForOne(res,next_gen,next_comp,0,IDELEMS(res));
3537    idDeleteComp(res,next_comp);
3538  }
3539  idSkipZeroes(res);
3540  return res;
3541}
3542
3543/*2
3544* transpose a module
3545*/
3546ideal idTransp(ideal a)
3547{
3548  int r = a->rank, c = IDELEMS(a);
3549  ideal b =  idInit(r,c);
3550
3551  for (int i=c; i>0; i--)
3552  {
3553    poly p=a->m[i-1];
3554    while(p!=NULL)
3555    {
3556      poly h=pHead(p);
3557      int co=pGetComp(h)-1;
3558      pSetComp(h,i);
3559      pSetmComp(h);
3560      b->m[co]=pAdd(b->m[co],h);
3561      pIter(p);
3562    }
3563  }
3564  return b;
3565}
3566
3567intvec * idQHomWeight(ideal id)
3568{
3569  poly head, tail;
3570  int k;
3571  int in=IDELEMS(id)-1, ready=0, all=0,
3572      coldim=pVariables, rowmax=2*coldim;
3573  if (in<0) return NULL;
3574  intvec *imat=new intvec(rowmax+1,coldim,0);
3575
3576  do
3577  {
3578    head = id->m[in--];
3579    if (head!=NULL)
3580    {
3581      tail = pNext(head);
3582      while (tail!=NULL)
3583      {
3584        all++;
3585        for (k=1;k<=coldim;k++)
3586          IMATELEM(*imat,all,k) = pGetExpDiff(head,tail,k);
3587        if (all==rowmax)
3588        {
3589          ivTriangIntern(imat, ready, all);
3590          if (ready==coldim)
3591          {
3592            delete imat;
3593            return NULL;
3594          }
3595        }
3596        pIter(tail);
3597      }
3598    }
3599  } while (in>=0);
3600  if (all>ready)
3601  {
3602    ivTriangIntern(imat, ready, all);
3603    if (ready==coldim)
3604    {
3605      delete imat;
3606      return NULL;
3607    }
3608  }
3609  intvec *result = ivSolveKern(imat, ready);
3610  delete imat;
3611  return result;
3612}
3613
3614BOOLEAN idIsZeroDim(ideal I)
3615{
3616  BOOLEAN *UsedAxis=(BOOLEAN *)omAlloc0(pVariables*sizeof(BOOLEAN));
3617  int i,n;
3618  poly po;
3619  BOOLEAN res=TRUE;
3620  for(i=IDELEMS(I)-1;i>=0;i--)
3621  {
3622    po=I->m[i];
3623    if ((po!=NULL) &&((n=pIsPurePower(po))!=0)) UsedAxis[n-1]=TRUE;
3624  }
3625  for(i=pVariables-1;i>=0;i--)
3626  {
3627    if(UsedAxis[i]==FALSE) {res=FALSE; break;} // not zero-dim.
3628  }
3629  omFreeSize(UsedAxis,pVariables*sizeof(BOOLEAN));
3630  return res;
3631}
3632
3633void idNormalize(ideal I)
3634{
3635  if (rField_has_simple_inverse()) return; /* Z/p, GF(p,n), R, long R/C */
3636  int i;
3637  poly p;
3638  for(i=IDELEMS(I)-1;i>=0;i--)
3639  {
3640    p=I->m[i] ;
3641    while(p!=NULL)
3642    {
3643      nNormalize(pGetCoeff(p));
3644      pIter(p);
3645    }
3646  }
3647}
Note: See TracBrowser for help on using the repository browser.