source: git/kernel/ideals.cc @ cbeafc2

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