source: git/kernel/ideals.cc @ ad3aa7

fieker-DuValspielwiese
Last change on this file since ad3aa7 was ad3aa7, checked in by Hans Schoenemann <hannes@…>, 13 years ago
id_MinDegW -> simpleideals.h, simpleideals.cc
  • Property mode set to 100644
File size: 56.6 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: ideals.cc 14320 2011-07-04 14:48:27Z hannes $ */
5/*
6* ABSTRACT - all basic methods to manipulate ideals
7*/
8
9/* includes */
10#include "mod2.h"
11
12
13#ifndef NDEBUG
14# define MYTEST 0
15#else /* ifndef NDEBUG */
16# define MYTEST 0
17#endif /* ifndef NDEBUG */
18
19#include <omalloc/omalloc.h>
20
21#include <misc/options.h>
22#include <misc/intvec.h>
23
24#include <coeffs/coeffs.h>
25#include <coeffs/numbers.h>
26
27#include <polys/polys.h>
28#include <polys/monomials/ring.h>
29#include <polys/matpol.h>
30#include <polys/weight.h>
31#include <polys/sparsmat.h>
32#include <polys/prCopy.h>
33#include <polys/nc/nc.h>
34
35
36#include <kernel/ideals.h>
37
38#include <kernel/febase.h>
39#include <kernel/kstd1.h>
40#include <kernel/syz.h>
41
42#include <kernel/longrat.h>
43
44
45/* #define WITH_OLD_MINOR */
46#define pCopy_noCheck(p) pCopy(p)
47
48static poly * idpower;
49/*collects the monomials in makemonoms, must be allocated befor*/
50static int idpowerpoint;
51/*index of the actual monomial in idpower*/
52static poly * givenideal;
53/*the ideal from which a power is computed*/
54
55/*0 implementation*/
56
57/*2
58*returns a minimized set of generators of h1
59*/
60ideal idMinBase (ideal h1)
61{
62  ideal h2, h3,h4,e;
63  int j,k;
64  int i,l,ll;
65  intvec * wth;
66  BOOLEAN homog;
67
68  homog = idHomModule(h1,currQuotient,&wth);
69  if (rHasGlobalOrdering(currRing))
70  {
71    if(!homog)
72    {
73      WarnS("minbase applies only to the local or homogeneous case");
74      e=idCopy(h1);
75      return e;
76    }
77    else
78    {
79      ideal re=kMin_std(h1,currQuotient,(tHomog)homog,&wth,h2,NULL,0,3);
80      idDelete(&re);
81      return h2;
82    }
83  }
84  e=idInit(1,h1->rank);
85  if (idIs0(h1))
86  {
87    return e;
88  }
89  pEnlargeSet(&(e->m),IDELEMS(e),15);
90  IDELEMS(e) = 16;
91  h2 = kStd(h1,currQuotient,isNotHomog,NULL);
92  h3 = idMaxIdeal(1);
93  h4=idMult(h2,h3);
94  idDelete(&h3);
95  h3=kStd(h4,currQuotient,isNotHomog,NULL);
96  k = IDELEMS(h3);
97  while ((k > 0) && (h3->m[k-1] == NULL)) k--;
98  j = -1;
99  l = IDELEMS(h2);
100  while ((l > 0) && (h2->m[l-1] == NULL)) l--;
101  for (i=l-1; i>=0; i--)
102  {
103    if (h2->m[i] != NULL)
104    {
105      ll = 0;
106      while ((ll < k) && ((h3->m[ll] == NULL)
107      || !pDivisibleBy(h3->m[ll],h2->m[i])))
108        ll++;
109      if (ll >= k)
110      {
111        j++;
112        if (j > IDELEMS(e)-1)
113        {
114          pEnlargeSet(&(e->m),IDELEMS(e),16);
115          IDELEMS(e) += 16;
116        }
117        e->m[j] = pCopy(h2->m[i]);
118      }
119    }
120  }
121  idDelete(&h2);
122  idDelete(&h3);
123  idDelete(&h4);
124  if (currQuotient!=NULL)
125  {
126    h3=idInit(1,e->rank);
127    h2=kNF(h3,currQuotient,e);
128    idDelete(&h3);
129    idDelete(&e);
130    e=h2;
131  }
132  idSkipZeroes(e);
133  return e;
134}
135
136
137/*3
138*multiplies p with t (!cas) or  (t-1)
139*the index of t is:1, so we have to shift all variables
140*p is NOT in the actual ring, it has no t
141*/
142static poly pMultWithT (poly p,BOOLEAN cas)
143{
144  /*qp is the working pointer in p*/
145  /*result is the result, qresult is the working pointer*/
146  /*pp is p in the actual ring(shifted), qpp the working pointer*/
147  poly result,qp,pp;
148  poly qresult=NULL;
149  poly qpp=NULL;
150  int  i,j,lex;
151  number n;
152
153  pp = NULL;
154  result = NULL;
155  qp = p;
156  while (qp != NULL)
157  {
158    i = 0;
159    if (result == NULL)
160    {/*first monomial*/
161      result = pInit();
162      qresult = result;
163    }
164    else
165    {
166      qresult->next = pInit();
167      pIter(qresult);
168    }
169    for (j=(currRing->N)-1; j>0; j--)
170    {
171      lex = pGetExp(qp,j);
172      pSetExp(qresult,j+1,lex);/*copy all variables*/
173    }
174    lex = pGetComp(qp);
175    pSetComp(qresult,lex);
176    n=nCopy(pGetCoeff(qp));
177    pSetCoeff0(qresult,n);
178    qresult->next = NULL;
179    pSetm(qresult);
180    /*qresult is now qp brought into the actual ring*/
181    if (cas)
182    { /*case: mult with t-1*/
183      pSetExp(qresult,1,0);
184      pSetm(qresult);
185      if (pp == NULL)
186      { /*first monomial*/
187        pp = pCopy(qresult);
188        qpp = pp;
189      }
190      else
191      {
192        qpp->next = pCopy(qresult);
193        pIter(qpp);
194      }
195      pGetCoeff(qpp)=nNeg(pGetCoeff(qpp));
196      /*now qpp contains -1*qp*/
197    }
198    pSetExp(qresult,1,1);/*this is mult. by t*/
199    pSetm(qresult);
200    pIter(qp);
201  }
202  /*
203  *now p is processed:
204  *result contains t*p
205  * if cas: pp contains -1*p (in the new ring)
206  */
207  if (cas)  qresult->next = pp;
208  /*  else      qresult->next = NULL;*/
209  return result;
210}
211
212/*2
213*initialized a field with r numbers between beg and end for the
214*procedure idNextChoise
215*/
216ideal idSectWithElim (ideal h1,ideal h2)
217// does not destroy h1,h2
218{
219  if (TEST_OPT_PROT) PrintS("intersect by elimination method\n");
220  assume(!idIs0(h1));
221  assume(!idIs0(h2));
222  assume(IDELEMS(h1)<=IDELEMS(h2));
223  assume(id_RankFreeModule(h1,currRing)==0);
224  assume(id_RankFreeModule(h2,currRing)==0);
225  // add a new variable:
226  int j;
227  ring origRing=currRing;
228  ring r=rCopy0(origRing);
229  r->N++;
230  r->block0[0]=1;
231  r->block1[0]= r->N;
232  omFree(r->order);
233  r->order=(int*)omAlloc0(3*sizeof(int*));
234  r->order[0]=ringorder_dp;
235  r->order[1]=ringorder_C;
236  char **names=(char**)omAlloc0(rVar(r) * sizeof(char_ptr));
237  for (j=0;j<r->N-1;j++) names[j]=r->names[j];
238  names[r->N-1]=omStrDup("@");
239  omFree(r->names);
240  r->names=names;
241  rComplete(r,TRUE);
242  // fetch h1, h2
243  ideal h;
244  h1=idrCopyR(h1,origRing,r);
245  h2=idrCopyR(h2,origRing,r);
246  // switch to temp. ring r
247  rChangeCurrRing(r);
248  // create 1-t, t
249  poly omt=p_One(currRing);
250  p_SetExp(omt,r->N,1,currRing);
251  poly t=p_Copy(omt,currRing);
252  p_Setm(omt,currRing);
253  omt=p_Neg(omt,currRing);
254  omt=p_Add_q(omt,pOne(),currRing);
255  // compute (1-t)*h1
256  h1=(ideal)mp_MultP((matrix)h1,omt,currRing);
257  // compute t*h2
258  h2=(ideal)mp_MultP((matrix)h2,pCopy(t),currRing);
259  // (1-t)h1 + t*h2
260  h=idInit(IDELEMS(h1)+IDELEMS(h2),1);
261  int l;
262  for (l=IDELEMS(h1)-1; l>=0; l--)
263  {
264    h->m[l] = h1->m[l];  h1->m[l]=NULL;
265  }
266  j=IDELEMS(h1);
267  for (l=IDELEMS(h2)-1; l>=0; l--)
268  {
269    h->m[l+j] = h2->m[l];  h2->m[l]=NULL;
270  }
271  idDelete(&h1);
272  idDelete(&h2);
273  // eliminate t:
274
275  ideal res=idElimination(h,t);
276  // cleanup
277  idDelete(&h);
278  res=idrMoveR(res,r,origRing);
279  rChangeCurrRing(origRing);
280  rKill(r);
281  return res;
282}
283/*2
284* h3 := h1 intersect h2
285*/
286ideal idSect (ideal h1,ideal h2)
287{
288  int i,j,k,length;
289  int flength = id_RankFreeModule(h1,currRing);
290  int slength = id_RankFreeModule(h2,currRing);
291  int rank=si_min(flength,slength);
292  if ((idIs0(h1)) || (idIs0(h2)))  return idInit(1,rank);
293
294  ideal first,second,temp,temp1,result;
295  poly p,q;
296
297  if (IDELEMS(h1)<IDELEMS(h2))
298  {
299    first = h1;
300    second = h2;
301  }
302  else
303  {
304    first = h2;
305    second = h1;
306    int t=flength; flength=slength; slength=t;
307  }
308  length  = si_max(flength,slength);
309  if (length==0)
310  {
311    if ((currQuotient==NULL)
312    && (currRing->OrdSgn==1)
313    && (!rIsPluralRing(currRing))
314    && ((TEST_V_INTERSECT_ELIM) || (!TEST_V_INTERSECT_SYZ)))
315      return idSectWithElim(first,second);
316    else length = 1;
317  }
318  if (TEST_OPT_PROT) PrintS("intersect by syzygy methods\n");
319  j = IDELEMS(first);
320
321  ring orig_ring=currRing;
322  ring syz_ring=rAssure_SyzComp(orig_ring,TRUE);
323  rChangeCurrRing(syz_ring);
324  rSetSyzComp(length, syz_ring);
325
326  while ((j>0) && (first->m[j-1]==NULL)) j--;
327  temp = idInit(j /*IDELEMS(first)*/+IDELEMS(second),length+j);
328  k = 0;
329  for (i=0;i<j;i++)
330  {
331    if (first->m[i]!=NULL)
332    {
333      if (syz_ring==orig_ring)
334        temp->m[k] = pCopy(first->m[i]);
335      else
336        temp->m[k] = prCopyR(first->m[i], orig_ring, syz_ring);
337      q = pOne();
338      pSetComp(q,i+1+length);
339      pSetmComp(q);
340      if (flength==0) p_Shift(&(temp->m[k]),1,currRing);
341      p = temp->m[k];
342      while (pNext(p)!=NULL) pIter(p);
343      pNext(p) = q;
344      k++;
345    }
346  }
347  for (i=0;i<IDELEMS(second);i++)
348  {
349    if (second->m[i]!=NULL)
350    {
351      if (syz_ring==orig_ring)
352        temp->m[k] = pCopy(second->m[i]);
353      else
354        temp->m[k] = prCopyR(second->m[i], orig_ring,currRing);
355      if (slength==0) p_Shift(&(temp->m[k]),1,currRing);
356      k++;
357    }
358  }
359  intvec *w=NULL;
360  temp1 = kStd(temp,currQuotient,testHomog,&w,NULL,length);
361  if (w!=NULL) delete w;
362  idDelete(&temp);
363  if(syz_ring!=orig_ring)
364    rChangeCurrRing(orig_ring);
365
366  result = idInit(IDELEMS(temp1),rank);
367  j = 0;
368  for (i=0;i<IDELEMS(temp1);i++)
369  {
370    if ((temp1->m[i]!=NULL)
371    && (p_GetComp(temp1->m[i],syz_ring)>length))
372    {
373      if(syz_ring==orig_ring)
374      {
375        p = temp1->m[i];
376      }
377      else
378      {
379        p = prMoveR(temp1->m[i], syz_ring,orig_ring);
380      }
381      temp1->m[i]=NULL;
382      while (p!=NULL)
383      {
384        q = pNext(p);
385        pNext(p) = NULL;
386        k = pGetComp(p)-1-length;
387        pSetComp(p,0);
388        pSetmComp(p);
389        /* Warning! multiply only from the left! it's very important for Plural */
390        result->m[j] = pAdd(result->m[j],pMult(p,pCopy(first->m[k])));
391        p = q;
392      }
393      j++;
394    }
395  }
396  if(syz_ring!=orig_ring)
397  {
398    rChangeCurrRing(syz_ring);
399    idDelete(&temp1);
400    rChangeCurrRing(orig_ring);
401    rKill(syz_ring);
402  }
403  else
404  {
405    idDelete(&temp1);
406  }
407
408  idSkipZeroes(result);
409  if (TEST_OPT_RETURN_SB)
410  {
411     w=NULL;
412     temp1=kStd(result,currQuotient,testHomog,&w);
413     if (w!=NULL) delete w;
414     idDelete(&result);
415     idSkipZeroes(temp1);
416     return temp1;
417  }
418  else //temp1=kInterRed(result,currQuotient);
419    return result;
420}
421
422/*2
423* ideal/module intersection for a list of objects
424* given as 'resolvente'
425*/
426ideal idMultSect(resolvente arg, int length)
427{
428  int i,j=0,k=0,syzComp,l,maxrk=-1,realrki;
429  ideal bigmat,tempstd,result;
430  poly p;
431  int isIdeal=0;
432  intvec * w=NULL;
433
434  /* find 0-ideals and max rank -----------------------------------*/
435  for (i=0;i<length;i++)
436  {
437    if (!idIs0(arg[i]))
438    {
439      realrki=id_RankFreeModule(arg[i],currRing);
440      k++;
441      j += IDELEMS(arg[i]);
442      if (realrki>maxrk) maxrk = realrki;
443    }
444    else
445    {
446      if (arg[i]!=NULL)
447      {
448        return idInit(1,arg[i]->rank);
449      }
450    }
451  }
452  if (maxrk == 0)
453  {
454    isIdeal = 1;
455    maxrk = 1;
456  }
457  /* init -----------------------------------------------------------*/
458  j += maxrk;
459  syzComp = k*maxrk;
460
461  ring orig_ring=currRing;
462  ring syz_ring=rAssure_SyzComp(orig_ring,TRUE);
463  rChangeCurrRing(syz_ring);
464  rSetSyzComp(syzComp, syz_ring);
465
466  bigmat = idInit(j,(k+1)*maxrk);
467  /* create unit matrices ------------------------------------------*/
468  for (i=0;i<maxrk;i++)
469  {
470    for (j=0;j<=k;j++)
471    {
472      p = pOne();
473      pSetComp(p,i+1+j*maxrk);
474      pSetmComp(p);
475      bigmat->m[i] = pAdd(bigmat->m[i],p);
476    }
477  }
478  /* enter given ideals ------------------------------------------*/
479  i = maxrk;
480  k = 0;
481  for (j=0;j<length;j++)
482  {
483    if (arg[j]!=NULL)
484    {
485      for (l=0;l<IDELEMS(arg[j]);l++)
486      {
487        if (arg[j]->m[l]!=NULL)
488        {
489          if (syz_ring==orig_ring)
490            bigmat->m[i] = pCopy(arg[j]->m[l]);
491          else
492            bigmat->m[i] = prCopyR(arg[j]->m[l], orig_ring,currRing);
493          p_Shift(&(bigmat->m[i]),k*maxrk+isIdeal,currRing);
494          i++;
495        }
496      }
497      k++;
498    }
499  }
500  /* std computation --------------------------------------------*/
501  tempstd = kStd(bigmat,currQuotient,testHomog,&w,NULL,syzComp);
502  if (w!=NULL) delete w;
503  idDelete(&bigmat);
504
505  if(syz_ring!=orig_ring)
506    rChangeCurrRing(orig_ring);
507
508  /* interprete result ----------------------------------------*/
509  result = idInit(IDELEMS(tempstd),maxrk);
510  k = 0;
511  for (j=0;j<IDELEMS(tempstd);j++)
512  {
513    if ((tempstd->m[j]!=NULL) && (p_GetComp(tempstd->m[j],syz_ring)>syzComp))
514    {
515      if (syz_ring==orig_ring)
516        p = pCopy(tempstd->m[j]);
517      else
518        p = prCopyR(tempstd->m[j], syz_ring,currRing);
519      p_Shift(&p,-syzComp-isIdeal,currRing);
520      result->m[k] = p;
521      k++;
522    }
523  }
524  /* clean up ----------------------------------------------------*/
525  if(syz_ring!=orig_ring)
526    rChangeCurrRing(syz_ring);
527  idDelete(&tempstd);
528  if(syz_ring!=orig_ring)
529  {
530    rChangeCurrRing(orig_ring);
531    rKill(syz_ring);
532  }
533  idSkipZeroes(result);
534  return result;
535}
536
537/*2
538*computes syzygies of h1,
539*if quot != NULL it computes in the quotient ring modulo "quot"
540*works always in a ring with ringorder_s
541*/
542static ideal idPrepare (ideal  h1, tHomog hom, int syzcomp, intvec **w)
543{
544  ideal   h2, h3;
545  int     i;
546  int     j,jj=0,k;
547  poly    p,q;
548
549  if (idIs0(h1)) return NULL;
550  k = id_RankFreeModule(h1,currRing);
551  h2=idCopy(h1);
552  i = IDELEMS(h2)-1;
553  if (k == 0)
554  {
555    for (j=0; j<=i; j++) p_Shift(&(h2->m[j]),1,currRing);
556    k = 1;
557  }
558  if (syzcomp<k)
559  {
560    Warn("syzcomp too low, should be %d instead of %d",k,syzcomp);
561    syzcomp = k;
562    rSetSyzComp(k,currRing);
563  }
564  h2->rank = syzcomp+i+1;
565
566  //if (hom==testHomog)
567  //{
568  //  if(idHomIdeal(h1,currQuotient))
569  //  {
570  //    hom=TRUE;
571  //  }
572  //}
573
574#if MYTEST
575#ifdef RDEBUG
576  Print("Prepare::h2: ");
577  idPrint(h2);
578
579  for(j=0;j<IDELEMS(h2);j++) pTest(h2->m[j]);
580
581#endif
582#endif
583
584  for (j=0; j<=i; j++)
585  {
586    p = h2->m[j];
587    q = pOne();
588    pSetComp(q,syzcomp+1+j);
589    pSetmComp(q);
590    if (p!=NULL)
591    {
592      while (pNext(p)) pIter(p);
593      p->next = q;
594    }
595    else
596      h2->m[j]=q;
597  }
598
599#ifdef PDEBUG
600  for(j=0;j<IDELEMS(h2);j++) pTest(h2->m[j]);
601
602#if MYTEST
603#ifdef RDEBUG
604  Print("Prepare::Input: ");
605  idPrint(h2);
606
607  Print("Prepare::currQuotient: ");
608  idPrint(currQuotient);
609#endif
610#endif
611
612#endif
613
614  idTest(h2);
615
616  h3 = kStd(h2,currQuotient,hom,w,NULL,syzcomp);
617
618#if MYTEST
619#ifdef RDEBUG
620  Print("Prepare::Output: ");
621  idPrint(h3);
622  for(j=0;j<IDELEMS(h2);j++) pTest(h3->m[j]);
623#endif
624#endif
625
626
627  idDelete(&h2);
628  return h3;
629}
630
631/*2
632* compute the syzygies of h1 in R/quot,
633* weights of components are in w
634* if setRegularity, return the regularity in deg
635* do not change h1,  w
636*/
637ideal idSyzygies (ideal  h1, tHomog h,intvec **w, BOOLEAN setSyzComp,
638                  BOOLEAN setRegularity, int *deg)
639{
640  ideal s_h1;
641  poly  p;
642  int   j, k, length=0,reg;
643  BOOLEAN isMonomial=TRUE;
644  int ii, idElemens_h1;
645
646  assume(h1 != NULL);
647
648  idElemens_h1=IDELEMS(h1);
649#ifdef PDEBUG
650  for(ii=0;ii<idElemens_h1 /*IDELEMS(h1)*/;ii++) pTest(h1->m[ii]);
651#endif
652  if (idIs0(h1))
653  {
654    ideal result=idFreeModule(idElemens_h1/*IDELEMS(h1)*/);
655    int curr_syz_limit=rGetCurrSyzLimit(currRing);
656    if (curr_syz_limit>0)
657    for (ii=0;ii<idElemens_h1/*IDELEMS(h1)*/;ii++)
658    {
659      if (h1->m[ii]!=NULL)
660        p_Shift(&h1->m[ii],curr_syz_limit,currRing);
661    }
662    return result;
663  }
664  int slength=(int)id_RankFreeModule(h1,currRing);
665  k=si_max(1,slength /*id_RankFreeModule(h1)*/);
666
667  assume(currRing != NULL);
668  ring orig_ring=currRing;
669  ring syz_ring=rAssure_SyzComp(orig_ring,TRUE);
670  rChangeCurrRing(syz_ring);
671
672  if (setSyzComp)
673    rSetSyzComp(k,syz_ring);
674
675  if (orig_ring != syz_ring)
676  {
677    s_h1=idrCopyR_NoSort(h1,orig_ring,syz_ring);
678  }
679  else
680  {
681    s_h1 = h1;
682  }
683
684  idTest(s_h1);
685
686  ideal s_h3=idPrepare(s_h1,h,k,w); // main (syz) GB computation
687
688  if (s_h3==NULL)
689  {
690    return idFreeModule( idElemens_h1 /*IDELEMS(h1)*/);
691  }
692
693  if (orig_ring != syz_ring)
694  {
695    idDelete(&s_h1);
696    for (j=0; j<IDELEMS(s_h3); j++)
697    {
698      if (s_h3->m[j] != NULL)
699      {
700        if (p_MinComp(s_h3->m[j],syz_ring) > k)
701          p_Shift(&s_h3->m[j], -k,currRing);
702        else
703          pDelete(&s_h3->m[j]);
704      }
705    }
706    idSkipZeroes(s_h3);
707    s_h3->rank -= k;
708    rChangeCurrRing(orig_ring);
709    s_h3 = idrMoveR_NoSort(s_h3, syz_ring, orig_ring);
710    rKill(syz_ring);
711    #ifdef HAVE_PLURAL
712    if (rIsPluralRing(currRing))
713    {
714      idDelMultiples(s_h3);
715      idSkipZeroes(s_h3);
716    }
717    #endif
718    idTest(s_h3);
719    return s_h3;
720  }
721
722  ideal e = idInit(IDELEMS(s_h3), s_h3->rank);
723
724  for (j=IDELEMS(s_h3)-1; j>=0; j--)
725  {
726    if (s_h3->m[j] != NULL)
727    {
728      if (p_MinComp(s_h3->m[j],syz_ring) <= k)
729      {
730        e->m[j] = s_h3->m[j];
731        isMonomial=isMonomial && (pNext(s_h3->m[j])==NULL);
732        pDelete(&pNext(s_h3->m[j]));
733        s_h3->m[j] = NULL;
734      }
735    }
736  }
737
738  idSkipZeroes(s_h3);
739  idSkipZeroes(e);
740
741  if ((deg != NULL)
742  && (!isMonomial)
743  && (!TEST_OPT_NOTREGULARITY)
744  && (setRegularity)
745  && (h==isHomog)
746  && (!rIsPluralRing(currRing))
747  )
748  {
749    ring dp_C_ring = rAssure_dp_C(syz_ring);
750    if (dp_C_ring != syz_ring)
751    {
752      rChangeCurrRing(dp_C_ring);
753      e = idrMoveR_NoSort(e, syz_ring, dp_C_ring);
754    }
755    resolvente res = sySchreyerResolvente(e,-1,&length,TRUE, TRUE);
756    intvec * dummy = syBetti(res,length,&reg, *w);
757    *deg = reg+2;
758    delete dummy;
759    for (j=0;j<length;j++)
760    {
761      if (res[j]!=NULL) idDelete(&(res[j]));
762    }
763    omFreeSize((ADDRESS)res,length*sizeof(ideal));
764    idDelete(&e);
765    if (dp_C_ring != syz_ring)
766    {
767      rChangeCurrRing(syz_ring);
768      rKill(dp_C_ring);
769    }
770  }
771  else
772  {
773    idDelete(&e);
774  }
775  idTest(s_h3);
776  if (currQuotient != NULL)
777  {
778    ideal ts_h3=kStd(s_h3,currQuotient,h,w);
779    idDelete(&s_h3);
780    s_h3 = ts_h3;
781  }
782  return s_h3;
783}
784
785/*2
786*/
787ideal idXXX (ideal  h1, int k)
788{
789  ideal s_h1;
790  int j;
791  intvec *w=NULL;
792
793  assume(currRing != NULL);
794  ring orig_ring=currRing;
795  ring syz_ring=rAssure_SyzComp(orig_ring,TRUE);
796  rChangeCurrRing(syz_ring);
797
798  rSetSyzComp(k,syz_ring);
799
800  if (orig_ring != syz_ring)
801  {
802    s_h1=idrCopyR_NoSort(h1,orig_ring, syz_ring);
803  }
804  else
805  {
806    s_h1 = h1;
807  }
808
809  ideal s_h3=kStd(s_h1,NULL,testHomog,&w,NULL,k);
810
811  if (s_h3==NULL)
812  {
813    return idFreeModule(IDELEMS(h1));
814  }
815
816  if (orig_ring != syz_ring)
817  {
818    idDelete(&s_h1);
819    idSkipZeroes(s_h3);
820    rChangeCurrRing(orig_ring);
821    s_h3 = idrMoveR_NoSort(s_h3, syz_ring, orig_ring);
822    rKill(syz_ring);
823    idTest(s_h3);
824    return s_h3;
825  }
826
827  idSkipZeroes(s_h3);
828  idTest(s_h3);
829  return s_h3;
830}
831
832/*
833*computes a standard basis for h1 and stores the transformation matrix
834* in ma
835*/
836ideal idLiftStd (ideal  h1, matrix* ma, tHomog hi, ideal * syz)
837{
838  int   i, j, k, t, inputIsIdeal=id_RankFreeModule(h1,currRing);
839  poly  p=NULL, q, qq;
840  intvec *w=NULL;
841
842  idDelete((ideal*)ma);
843  BOOLEAN lift3=FALSE;
844  if (syz!=NULL) { lift3=TRUE; idDelete(syz); }
845  if (idIs0(h1))
846  {
847    *ma=mpNew(1,0);
848    if (lift3)
849    {
850      *syz=idFreeModule(IDELEMS(h1));
851      int curr_syz_limit=rGetCurrSyzLimit(currRing);
852      if (curr_syz_limit>0)
853      for (int ii=0;ii<IDELEMS(h1);ii++)
854      {
855        if (h1->m[ii]!=NULL)
856          p_Shift(&h1->m[ii],curr_syz_limit,currRing);
857      }
858    }
859    return idInit(1,h1->rank);
860  }
861
862  BITSET save_verbose=verbose;
863
864  k=si_max(1,(int)id_RankFreeModule(h1,currRing));
865
866  if ((k==1) && (!lift3)) verbose |=Sy_bit(V_IDLIFT);
867
868  ring orig_ring = currRing;
869  ring syz_ring = rAssure_SyzComp(orig_ring,TRUE);
870  rChangeCurrRing(syz_ring);
871  rSetSyzComp(k,syz_ring);
872
873  ideal s_h1=h1;
874
875  if (orig_ring != syz_ring)
876    s_h1 = idrCopyR_NoSort(h1,orig_ring,syz_ring);
877  else
878    s_h1 = h1;
879
880  ideal s_h3=idPrepare(s_h1,hi,k,&w); // main (syz) GB computation
881
882  ideal s_h2 = idInit(IDELEMS(s_h3), s_h3->rank);
883
884  if (lift3) (*syz)=idInit(IDELEMS(s_h3),IDELEMS(h1));
885
886  if (w!=NULL) delete w;
887  i = 0;
888
889  // now sort the result, SB : leave in s_h3
890  //                      T:  put in s_h2
891  //                      syz: put in *syz
892  for (j=0; j<IDELEMS(s_h3); j++)
893  {
894    if (s_h3->m[j] != NULL)
895    {
896      //if (p_MinComp(s_h3->m[j],syz_ring) <= k)
897      if (pGetComp(s_h3->m[j]) <= k) // syz_ring == currRing
898      {
899        i++;
900        q = s_h3->m[j];
901        while (pNext(q) != NULL)
902        {
903          if (pGetComp(pNext(q)) > k)
904          {
905            s_h2->m[j] = pNext(q);
906            pNext(q) = NULL;
907          }
908          else
909          {
910            pIter(q);
911          }
912        }
913        if (!inputIsIdeal) p_Shift(&(s_h3->m[j]), -1,currRing);
914      }
915      else
916      {
917        // we a syzygy here:
918        if (lift3)
919        {
920          p_Shift(&s_h3->m[j], -k,currRing);
921          (*syz)->m[j]=s_h3->m[j];
922          s_h3->m[j]=NULL;
923        }
924        else
925          pDelete(&(s_h3->m[j]));
926      }
927    }
928  }
929  idSkipZeroes(s_h3);
930  //extern char * iiStringMatrix(matrix im, int dim,char ch);
931  //PrintS("SB: ----------------------------------------\n");
932  //PrintS(iiStringMatrix((matrix)s_h3,k,'\n'));
933  //PrintLn();
934  //PrintS("T: ----------------------------------------\n");
935  //PrintS(iiStringMatrix((matrix)s_h2,h1->rank,'\n'));
936  //PrintLn();
937
938  if (lift3) idSkipZeroes(*syz);
939
940  j = IDELEMS(s_h1);
941
942
943  if (syz_ring!=orig_ring)
944  {
945    idDelete(&s_h1);
946    rChangeCurrRing(orig_ring);
947  }
948
949  *ma = mpNew(j,i);
950
951  i = 1;
952  for (j=0; j<IDELEMS(s_h2); j++)
953  {
954    if (s_h2->m[j] != NULL)
955    {
956      q = prMoveR( s_h2->m[j], syz_ring,orig_ring);
957      s_h2->m[j] = NULL;
958
959      while (q != NULL)
960      {
961        p = q;
962        pIter(q);
963        pNext(p) = NULL;
964        t=pGetComp(p);
965        pSetComp(p,0);
966        pSetmComp(p);
967        MATELEM(*ma,t-k,i) = pAdd(MATELEM(*ma,t-k,i),p);
968      }
969      i++;
970    }
971  }
972  idDelete(&s_h2);
973
974  for (i=0; i<IDELEMS(s_h3); i++)
975  {
976    s_h3->m[i] = prMoveR_NoSort(s_h3->m[i], syz_ring,orig_ring);
977  }
978  if (lift3)
979  {
980    for (i=0; i<IDELEMS(*syz); i++)
981    {
982      (*syz)->m[i] = prMoveR_NoSort((*syz)->m[i], syz_ring,orig_ring);
983    }
984  }
985
986  if (syz_ring!=orig_ring) rKill(syz_ring);
987  verbose = save_verbose;
988  return s_h3;
989}
990
991static void idPrepareStd(ideal s_temp, int k)
992{
993  int j,rk=id_RankFreeModule(s_temp,currRing);
994  poly p,q;
995
996  if (rk == 0)
997  {
998    for (j=0; j<IDELEMS(s_temp); j++)
999    {
1000      if (s_temp->m[j]!=NULL) pSetCompP(s_temp->m[j],1);
1001    }
1002    k = si_max(k,1);
1003  }
1004  for (j=0; j<IDELEMS(s_temp); j++)
1005  {
1006    if (s_temp->m[j]!=NULL)
1007    {
1008      p = s_temp->m[j];
1009      q = pOne();
1010      //pGetCoeff(q)=nNeg(pGetCoeff(q));   //set q to -1
1011      pSetComp(q,k+1+j);
1012      pSetmComp(q);
1013      while (pNext(p)) pIter(p);
1014      pNext(p) = q;
1015    }
1016  }
1017}
1018
1019/*2
1020*computes a representation of the generators of submod with respect to those
1021* of mod
1022*/
1023
1024ideal idLift(ideal mod, ideal submod,ideal *rest, BOOLEAN goodShape,
1025             BOOLEAN isSB, BOOLEAN divide, matrix *unit)
1026{
1027  int lsmod =id_RankFreeModule(submod,currRing), i, j, k;
1028  int comps_to_add=0;
1029  poly p;
1030
1031  if (idIs0(submod))
1032  {
1033    if (unit!=NULL)
1034    {
1035      *unit=mpNew(1,1);
1036      MATELEM(*unit,1,1)=pOne();
1037    }
1038    if (rest!=NULL)
1039    {
1040      *rest=idInit(1,mod->rank);
1041    }
1042    return idInit(1,mod->rank);
1043  }
1044  if (idIs0(mod)) /* and not idIs0(submod) */
1045  {
1046    WerrorS("2nd module does not lie in the first");
1047    #if 0
1048    if (unit!=NULL)
1049    {
1050      i=IDELEMS(submod);
1051      *unit=mpNew(i,i);
1052      for (j=i;j>0;j--)
1053      {
1054        MATELEM(*unit,j,j)=pOne();
1055      }
1056    }
1057    if (rest!=NULL)
1058    {
1059      *rest=idCopy(submod);
1060    }
1061    return idInit(1,mod->rank);
1062    #endif
1063    return idInit(IDELEMS(submod),submod->rank);
1064  }
1065  if (unit!=NULL)
1066  {
1067    comps_to_add = IDELEMS(submod);
1068    while ((comps_to_add>0) && (submod->m[comps_to_add-1]==NULL))
1069      comps_to_add--;
1070  }
1071  k=si_max(id_RankFreeModule(mod,currRing),id_RankFreeModule(submod,currRing));
1072  if  ((k!=0) && (lsmod==0)) lsmod=1;
1073  k=si_max(k,(int)mod->rank);
1074  if (k<submod->rank) { WarnS("rk(submod) > rk(mod) ?");k=submod->rank; }
1075
1076  ring orig_ring=currRing;
1077  ring syz_ring=rAssure_SyzComp(orig_ring,TRUE);
1078  rChangeCurrRing(syz_ring);
1079  rSetSyzComp(k,syz_ring);
1080
1081  ideal s_mod, s_temp;
1082  if (orig_ring != syz_ring)
1083  {
1084    s_mod = idrCopyR_NoSort(mod,orig_ring,syz_ring);
1085    s_temp = idrCopyR_NoSort(submod,orig_ring,syz_ring);
1086  }
1087  else
1088  {
1089    s_mod = mod;
1090    s_temp = idCopy(submod);
1091  }
1092  ideal s_h3;
1093  if (isSB)
1094  {
1095    s_h3 = idCopy(s_mod);
1096    idPrepareStd(s_h3, k+comps_to_add);
1097  }
1098  else
1099  {
1100    s_h3 = idPrepare(s_mod,(tHomog)FALSE,k+comps_to_add,NULL);
1101  }
1102  if (!goodShape)
1103  {
1104    for (j=0;j<IDELEMS(s_h3);j++)
1105    {
1106      if ((s_h3->m[j] != NULL) && (pMinComp(s_h3->m[j]) > k))
1107        pDelete(&(s_h3->m[j]));
1108    }
1109  }
1110  idSkipZeroes(s_h3);
1111  if (lsmod==0)
1112  {
1113    for (j=IDELEMS(s_temp);j>0;j--)
1114    {
1115      if (s_temp->m[j-1]!=NULL)
1116        p_Shift(&(s_temp->m[j-1]),1,currRing);
1117    }
1118  }
1119  if (unit!=NULL)
1120  {
1121    for(j = 0;j<comps_to_add;j++)
1122    {
1123      p = s_temp->m[j];
1124      if (p!=NULL)
1125      {
1126        while (pNext(p)!=NULL) pIter(p);
1127        pNext(p) = pOne();
1128        pIter(p);
1129        pSetComp(p,1+j+k);
1130        pSetmComp(p);
1131        p = pNeg(p);
1132      }
1133    }
1134  }
1135  ideal s_result = kNF(s_h3,currQuotient,s_temp,k);
1136  s_result->rank = s_h3->rank;
1137  ideal s_rest = idInit(IDELEMS(s_result),k);
1138  idDelete(&s_h3);
1139  idDelete(&s_temp);
1140
1141  for (j=0;j<IDELEMS(s_result);j++)
1142  {
1143    if (s_result->m[j]!=NULL)
1144    {
1145      if (pGetComp(s_result->m[j])<=k)
1146      {
1147        if (!divide)
1148        {
1149          if (isSB)
1150          {
1151            WarnS("first module not a standardbasis\n"
1152              "// ** or second not a proper submodule");
1153          }
1154          else
1155            WerrorS("2nd module does not lie in the first");
1156          idDelete(&s_result);
1157          idDelete(&s_rest);
1158          s_result=idInit(IDELEMS(submod),submod->rank);
1159          break;
1160        }
1161        else
1162        {
1163          p = s_rest->m[j] = s_result->m[j];
1164          while ((pNext(p)!=NULL) && (pGetComp(pNext(p))<=k)) pIter(p);
1165          s_result->m[j] = pNext(p);
1166          pNext(p) = NULL;
1167        }
1168      }
1169      p_Shift(&(s_result->m[j]),-k,currRing);
1170      pNeg(s_result->m[j]);
1171    }
1172  }
1173  if ((lsmod==0) && (!idIs0(s_rest)))
1174  {
1175    for (j=IDELEMS(s_rest);j>0;j--)
1176    {
1177      if (s_rest->m[j-1]!=NULL)
1178      {
1179        p_Shift(&(s_rest->m[j-1]),-1,currRing);
1180        s_rest->m[j-1] = s_rest->m[j-1];
1181      }
1182    }
1183  }
1184  if(syz_ring!=orig_ring)
1185  {
1186    idDelete(&s_mod);
1187    rChangeCurrRing(orig_ring);
1188    s_result = idrMoveR_NoSort(s_result, syz_ring, orig_ring);
1189    s_rest = idrMoveR_NoSort(s_rest, syz_ring, orig_ring);
1190    rKill(syz_ring);
1191  }
1192  if (rest!=NULL)
1193    *rest = s_rest;
1194  else
1195    idDelete(&s_rest);
1196//idPrint(s_result);
1197  if (unit!=NULL)
1198  {
1199    *unit=mpNew(comps_to_add,comps_to_add);
1200    int i;
1201    for(i=0;i<IDELEMS(s_result);i++)
1202    {
1203      poly p=s_result->m[i];
1204      poly q=NULL;
1205      while(p!=NULL)
1206      {
1207        if(pGetComp(p)<=comps_to_add)
1208        {
1209          pSetComp(p,0);
1210          if (q!=NULL)
1211          {
1212            pNext(q)=pNext(p);
1213          }
1214          else
1215          {
1216            pIter(s_result->m[i]);
1217          }
1218          pNext(p)=NULL;
1219          MATELEM(*unit,i+1,i+1)=pAdd(MATELEM(*unit,i+1,i+1),p);
1220          if(q!=NULL)   p=pNext(q);
1221          else          p=s_result->m[i];
1222        }
1223        else
1224        {
1225          q=p;
1226          pIter(p);
1227        }
1228      }
1229      p_Shift(&s_result->m[i],-comps_to_add,currRing);
1230    }
1231  }
1232  return s_result;
1233}
1234
1235/*2
1236*computes division of P by Q with remainder up to (w-weighted) degree n
1237*P, Q, and w are not changed
1238*/
1239void idLiftW(ideal P,ideal Q,int n,matrix &T, ideal &R,short *w)
1240{
1241  long N=0;
1242  int i;
1243  for(i=IDELEMS(Q)-1;i>=0;i--)
1244    if(w==NULL)
1245      N=si_max(N,pDeg(Q->m[i]));
1246    else
1247      N=si_max(N,pDegW(Q->m[i],w));
1248  N+=n;
1249
1250  T=mpNew(IDELEMS(Q),IDELEMS(P));
1251  R=idInit(IDELEMS(P),P->rank);
1252
1253  for(i=IDELEMS(P)-1;i>=0;i--)
1254  {
1255    poly p;
1256    if(w==NULL)
1257      p=ppJet(P->m[i],N);
1258    else
1259      p=ppJetW(P->m[i],N,w);
1260
1261    int j=IDELEMS(Q)-1;
1262    while(p!=NULL)
1263    {
1264      if(pDivisibleBy(Q->m[j],p))
1265      {
1266        poly p0=p_DivideM(pHead(p),pHead(Q->m[j]),currRing);
1267        if(w==NULL)
1268          p=pJet(pSub(p,ppMult_mm(Q->m[j],p0)),N);
1269        else
1270          p=pJetW(pSub(p,ppMult_mm(Q->m[j],p0)),N,w);
1271        pNormalize(p);
1272        if((w==NULL)&&(pDeg(p0)>n)||(w!=NULL)&&(pDegW(p0,w)>n))
1273          pDelete(&p0);
1274        else
1275          MATELEM(T,j+1,i+1)=pAdd(MATELEM(T,j+1,i+1),p0);
1276        j=IDELEMS(Q)-1;
1277      }
1278      else
1279      {
1280        if(j==0)
1281        {
1282          poly p0=p;
1283          pIter(p);
1284          pNext(p0)=NULL;
1285          if(((w==NULL)&&(pDeg(p0)>n))
1286          ||((w!=NULL)&&(pDegW(p0,w)>n)))
1287            pDelete(&p0);
1288          else
1289            R->m[i]=pAdd(R->m[i],p0);
1290          j=IDELEMS(Q)-1;
1291        }
1292        else
1293          j--;
1294      }
1295    }
1296  }
1297}
1298
1299/*2
1300*computes the quotient of h1,h2 : internal routine for idQuot
1301*BEWARE: the returned ideals may contain incorrectly ordered polys !
1302*
1303*/
1304static ideal idInitializeQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb,
1305                               BOOLEAN *addOnlyOne, int *kkmax)
1306{
1307  ideal temph1;
1308  poly     p,q = NULL;
1309  int i,l,ll,k,kkk,kmax;
1310  int j = 0;
1311  int k1 = id_RankFreeModule(h1,currRing);
1312  int k2 = id_RankFreeModule(h2,currRing);
1313  tHomog   hom=isNotHomog;
1314
1315  k=si_max(k1,k2);
1316  if (k==0)
1317    k = 1;
1318  if ((k2==0) && (k>1)) *addOnlyOne = FALSE;
1319
1320  intvec * weights;
1321  hom = (tHomog)idHomModule(h1,currQuotient,&weights);
1322  if (/**addOnlyOne &&*/ (!h1IsStb))
1323    temph1 = kStd(h1,currQuotient,hom,&weights,NULL);
1324  else
1325    temph1 = idCopy(h1);
1326  if (weights!=NULL) delete weights;
1327  idTest(temph1);
1328/*--- making a single vector from h2 ---------------------*/
1329  for (i=0; i<IDELEMS(h2); i++)
1330  {
1331    if (h2->m[i] != NULL)
1332    {
1333      p = pCopy(h2->m[i]);
1334      if (k2 == 0)
1335        p_Shift(&p,j*k+1,currRing);
1336      else
1337        p_Shift(&p,j*k,currRing);
1338      q = pAdd(q,p);
1339      j++;
1340    }
1341  }
1342  *kkmax = kmax = j*k+1;
1343/*--- adding a monomial for the result (syzygy) ----------*/
1344  p = q;
1345  while (pNext(p)!=NULL) pIter(p);
1346  pNext(p) = pOne();
1347  pIter(p);
1348  pSetComp(p,kmax);
1349  pSetmComp(p);
1350/*--- constructing the big matrix ------------------------*/
1351  ideal h4 = idInit(16,kmax+k-1);
1352  h4->m[0] = q;
1353  if (k2 == 0)
1354  {
1355    if (k > IDELEMS(h4))
1356    {
1357      pEnlargeSet(&(h4->m),IDELEMS(h4),k-IDELEMS(h4));
1358      IDELEMS(h4) = k;
1359    }
1360    for (i=1; i<k; i++)
1361    {
1362      if (h4->m[i-1]!=NULL)
1363      {
1364        p = pCopy_noCheck(h4->m[i-1]);
1365        p_Shift(&p,1,currRing);
1366        h4->m[i] = p;
1367      }
1368    }
1369  }
1370  idSkipZeroes(h4);
1371  kkk = IDELEMS(h4);
1372  i = IDELEMS(temph1);
1373  for (l=0; l<i; l++)
1374  {
1375    if(temph1->m[l]!=NULL)
1376    {
1377      for (ll=0; ll<j; ll++)
1378      {
1379        p = pCopy(temph1->m[l]);
1380        if (k1 == 0)
1381          p_Shift(&p,ll*k+1,currRing);
1382        else
1383          p_Shift(&p,ll*k,currRing);
1384        if (kkk >= IDELEMS(h4))
1385        {
1386          pEnlargeSet(&(h4->m),IDELEMS(h4),16);
1387          IDELEMS(h4) += 16;
1388        }
1389        h4->m[kkk] = p;
1390        kkk++;
1391      }
1392    }
1393  }
1394/*--- if h2 goes in as single vector - the h1-part is just SB ---*/
1395  if (*addOnlyOne)
1396  {
1397    idSkipZeroes(h4);
1398    p = h4->m[0];
1399    for (i=0;i<IDELEMS(h4)-1;i++)
1400    {
1401      h4->m[i] = h4->m[i+1];
1402    }
1403    h4->m[IDELEMS(h4)-1] = p;
1404    test |= Sy_bit(OPT_SB_1);
1405  }
1406  idDelete(&temph1);
1407  return h4;
1408}
1409/*2
1410*computes the quotient of h1,h2
1411*/
1412ideal idQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb, BOOLEAN resultIsIdeal)
1413{
1414  // first check for special case h1:(0)
1415  if (idIs0(h2))
1416  {
1417    ideal res;
1418    if (resultIsIdeal)
1419    {
1420      res = idInit(1,1);
1421      res->m[0] = pOne();
1422    }
1423    else
1424      res = idFreeModule(h1->rank);
1425    return res;
1426  }
1427  BITSET old_test=test;
1428  int i,l,ll,k,kkk,kmax;
1429  BOOLEAN  addOnlyOne=TRUE;
1430  tHomog   hom=isNotHomog;
1431  intvec * weights1;
1432
1433  ideal s_h4 = idInitializeQuot (h1,h2,h1IsStb,&addOnlyOne,&kmax);
1434
1435  hom = (tHomog)idHomModule(s_h4,currQuotient,&weights1);
1436
1437  ring orig_ring=currRing;
1438  ring syz_ring=rAssure_SyzComp(orig_ring,TRUE);
1439  rChangeCurrRing(syz_ring);
1440  rSetSyzComp(kmax-1,syz_ring);
1441  if (orig_ring!=syz_ring)
1442  //  s_h4 = idrMoveR_NoSort(s_h4,orig_ring, syz_ring);
1443    s_h4 = idrMoveR(s_h4,orig_ring, syz_ring);
1444  idTest(s_h4);
1445  #if 0
1446  void ipPrint_MA0(matrix m, const char *name);
1447  matrix m=idModule2Matrix(idCopy(s_h4));
1448  PrintS("start:\n");
1449  ipPrint_MA0(m,"Q");
1450  idDelete((ideal *)&m);
1451  PrintS("last elem:");wrp(s_h4->m[IDELEMS(s_h4)-1]);PrintLn();
1452  #endif
1453  ideal s_h3;
1454  if (addOnlyOne)
1455  {
1456    s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,0/*kmax-1*/,IDELEMS(s_h4)-1);
1457  }
1458  else
1459  {
1460    s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,kmax-1);
1461  }
1462  test = old_test;
1463  #if 0
1464  // only together with the above debug stuff
1465  idSkipZeroes(s_h3);
1466  m=idModule2Matrix(idCopy(s_h3));
1467  Print("result, kmax=%d:\n",kmax);
1468  ipPrint_MA0(m,"S");
1469  idDelete((ideal *)&m);
1470  #endif
1471  idTest(s_h3);
1472  if (weights1!=NULL) delete weights1;
1473  idDelete(&s_h4);
1474
1475  for (i=0;i<IDELEMS(s_h3);i++)
1476  {
1477    if ((s_h3->m[i]!=NULL) && (pGetComp(s_h3->m[i])>=kmax))
1478    {
1479      if (resultIsIdeal)
1480        p_Shift(&s_h3->m[i],-kmax,currRing);
1481      else
1482        p_Shift(&s_h3->m[i],-kmax+1,currRing);
1483    }
1484    else
1485      pDelete(&s_h3->m[i]);
1486  }
1487  if (resultIsIdeal)
1488    s_h3->rank = 1;
1489  else
1490    s_h3->rank = h1->rank;
1491  if(syz_ring!=orig_ring)
1492  {
1493    rChangeCurrRing(orig_ring);
1494    s_h3 = idrMoveR_NoSort(s_h3, syz_ring, orig_ring);
1495    rKill(syz_ring);
1496  }
1497  idSkipZeroes(s_h3);
1498  idTest(s_h3);
1499  return s_h3;
1500}
1501
1502/*2
1503* eliminate delVar (product of vars) in h1
1504*/
1505ideal idElimination (ideal h1,poly delVar,intvec *hilb)
1506{
1507  int    i,j=0,k,l;
1508  ideal  h,hh, h3;
1509  int    *ord,*block0,*block1;
1510  int    ordersize=2;
1511  int    **wv;
1512  tHomog hom;
1513  intvec * w;
1514  ring tmpR;
1515  ring origR = currRing;
1516
1517  if (delVar==NULL)
1518  {
1519    return idCopy(h1);
1520  }
1521  if ((currQuotient!=NULL) && rIsPluralRing(origR))
1522  {
1523    WerrorS("cannot eliminate in a qring");
1524    return idCopy(h1);
1525  }
1526  if (idIs0(h1)) return idInit(1,h1->rank);
1527#ifdef HAVE_PLURAL
1528  if (rIsPluralRing(origR))
1529    /* in the NC case, we have to check the admissibility of */
1530    /* the subalgebra to be intersected with */
1531  {
1532    if ((ncRingType(origR) != nc_skew) && (ncRingType(origR) != nc_exterior)) /* in (quasi)-commutative algebras every subalgebra is admissible */
1533    {
1534      if (nc_CheckSubalgebra(delVar,origR))
1535      {
1536        WerrorS("no elimination is possible: subalgebra is not admissible");
1537        return idCopy(h1);
1538      }
1539    }
1540  }
1541#endif
1542  hom=(tHomog)idHomModule(h1,NULL,&w); //sets w to weight vector or NULL
1543  h3=idInit(16,h1->rank);
1544  for (k=0;; k++)
1545  {
1546    if (origR->order[k]!=0) ordersize++;
1547    else break;
1548  }
1549#if 0
1550  if (rIsPluralRing(origR)) // we have too keep the odering: it may be needed
1551                            // for G-algebra
1552  {
1553    for (k=0;k<ordersize-1; k++)
1554    {
1555      block0[k+1] = origR->block0[k];
1556      block1[k+1] = origR->block1[k];
1557      ord[k+1] = origR->order[k];
1558      if (origR->wvhdl[k]!=NULL) wv[k+1] = (int*) omMemDup(origR->wvhdl[k]);
1559    }
1560  }
1561  else
1562  {
1563    block0[1] = 1;
1564    block1[1] = (currRing->N);
1565    if (origR->OrdSgn==1) ord[1] = ringorder_wp;
1566    else                  ord[1] = ringorder_ws;
1567    wv[1]=(int*)omAlloc0((currRing->N)*sizeof(int));
1568    double wNsqr = (double)2.0 / (double)(currRing->N);
1569    wFunctional = wFunctionalBuch;
1570    int  *x= (int * )omAlloc(2 * ((currRing->N) + 1) * sizeof(int));
1571    int sl=IDELEMS(h1) - 1;
1572    wCall(h1->m, sl, x, wNsqr);
1573    for (sl = (currRing->N); sl!=0; sl--)
1574      wv[1][sl-1] = x[sl + (currRing->N) + 1];
1575    omFreeSize((ADDRESS)x, 2 * ((currRing->N) + 1) * sizeof(int));
1576
1577    ord[2]=ringorder_C;
1578    ord[3]=0;
1579  }
1580#else
1581#endif
1582  if ((hom==TRUE) && (origR->OrdSgn==1) && (!rIsPluralRing(origR)))
1583  {
1584    #if 1
1585    // we change to an ordering:
1586    // aa(1,1,1,...,0,0,0),wp(...),C
1587    // this seems to be better than version 2 below,
1588    // according to Tst/../elimiate_[3568].tat (- 17 %)
1589    ord=(int*)omAlloc0(4*sizeof(int));
1590    block0=(int*)omAlloc0(4*sizeof(int));
1591    block1=(int*)omAlloc0(4*sizeof(int));
1592    wv=(int**) omAlloc0(4*sizeof(int**));
1593    block0[0] = block0[1] = 1;
1594    block1[0] = block1[1] = rVar(origR);
1595    wv[0]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
1596    // use this special ordering: like ringorder_a, except that pFDeg, pWeights
1597    // ignore it
1598    ord[0] = ringorder_aa;
1599    for (j=0;j<rVar(origR);j++)
1600      if (pGetExp(delVar,j+1)!=0) wv[0][j]=1;
1601    BOOLEAN wp=FALSE;
1602    for (j=0;j<rVar(origR);j++)
1603      if (pWeight(j+1,origR)!=1) { wp=TRUE;break; }
1604    if (wp)
1605    {
1606      wv[1]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
1607      for (j=0;j<rVar(origR);j++)
1608        wv[1][j]=pWeight(j+1,origR);
1609      ord[1] = ringorder_wp;
1610    }
1611    else
1612      ord[1] = ringorder_dp;
1613    #else
1614    // we change to an ordering:
1615    // a(w1,...wn),wp(1,...0.....),C
1616    ord=(int*)omAlloc0(4*sizeof(int));
1617    block0=(int*)omAlloc0(4*sizeof(int));
1618    block1=(int*)omAlloc0(4*sizeof(int));
1619    wv=(int**) omAlloc0(4*sizeof(int**));
1620    block0[0] = block0[1] = 1;
1621    block1[0] = block1[1] = rVar(origR);
1622    wv[0]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
1623    wv[1]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
1624    ord[0] = ringorder_a;
1625    for (j=0;j<rVar(origR);j++)
1626      wv[0][j]=pWeight(j+1,origR);
1627    ord[1] = ringorder_wp;
1628    for (j=0;j<rVar(origR);j++)
1629      if (pGetExp(delVar,j+1)!=0) wv[1][j]=1;
1630    #endif
1631    ord[2] = ringorder_C;
1632    ord[3] = 0;
1633  }
1634  else
1635  {
1636    // we change to an ordering:
1637    // aa(....),orig_ordering
1638    ord=(int*)omAlloc0(ordersize*sizeof(int));
1639    block0=(int*)omAlloc0(ordersize*sizeof(int));
1640    block1=(int*)omAlloc0(ordersize*sizeof(int));
1641    wv=(int**) omAlloc0(ordersize*sizeof(int**));
1642    for (k=0;k<ordersize-1; k++)
1643    {
1644      block0[k+1] = origR->block0[k];
1645      block1[k+1] = origR->block1[k];
1646      ord[k+1] = origR->order[k];
1647      if (origR->wvhdl[k]!=NULL) wv[k+1] = (int*) omMemDup(origR->wvhdl[k]);
1648    }
1649    block0[0] = 1;
1650    block1[0] = rVar(origR);
1651    wv[0]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
1652    for (j=0;j<rVar(origR);j++)
1653      if (pGetExp(delVar,j+1)!=0) wv[0][j]=1;
1654    // use this special ordering: like ringorder_a, except that pFDeg, pWeights
1655    // ignore it
1656    ord[0] = ringorder_aa;
1657  }
1658  // fill in tmp ring to get back the data later on
1659  tmpR  = rCopy0(origR,FALSE,FALSE); // qring==NULL
1660  //rUnComplete(tmpR);
1661  tmpR->p_Procs=NULL;
1662  tmpR->order = ord;
1663  tmpR->block0 = block0;
1664  tmpR->block1 = block1;
1665  tmpR->wvhdl = wv;
1666  rComplete(tmpR, 1);
1667
1668#ifdef HAVE_PLURAL
1669  /* update nc structure on tmpR */
1670  if (rIsPluralRing(origR))
1671  {
1672    if ( nc_rComplete(origR, tmpR, false) ) // no quotient ideal!
1673    {
1674      Werror("no elimination is possible: ordering condition is violated");
1675      // cleanup
1676      rDelete(tmpR);
1677      if (w!=NULL)
1678        delete w;
1679      return idCopy(h1);
1680    }
1681  }
1682#endif
1683  // change into the new ring
1684  //pChangeRing((currRing->N),currRing->OrdSgn,ord,block0,block1,wv);
1685  rChangeCurrRing(tmpR);
1686
1687  //h = idInit(IDELEMS(h1),h1->rank);
1688  // fetch data from the old ring
1689  //for (k=0;k<IDELEMS(h1);k++) h->m[k] = prCopyR( h1->m[k], origR);
1690  h=idrCopyR(h1,origR,currRing);
1691  if (origR->qideal!=NULL)
1692  {
1693    WarnS("eliminate in q-ring: experimental");
1694    ideal q=idrCopyR(origR->qideal,origR,currRing);
1695    ideal s=idSimpleAdd(h,q);
1696    idDelete(&h);
1697    idDelete(&q);
1698    h=s;
1699  }
1700  // compute kStd
1701#if 1
1702  //rWrite(tmpR);PrintLn();
1703  BITSET save=test;
1704  //test |=1;
1705  //Print("h: %d gen, rk=%d\n",IDELEMS(h),h->rank);
1706  //extern char * showOption();
1707  //Print("%s\n",showOption());
1708  hh = kStd(h,NULL,hom,&w,hilb);
1709  test=save;
1710  idDelete(&h);
1711#else
1712  extern ideal kGroebner(ideal F, ideal Q);
1713  hh=kGroebner(h,NULL);
1714#endif
1715  // go back to the original ring
1716  rChangeCurrRing(origR);
1717  i = IDELEMS(hh)-1;
1718  while ((i >= 0) && (hh->m[i] == NULL)) i--;
1719  j = -1;
1720  // fetch data from temp ring
1721  for (k=0; k<=i; k++)
1722  {
1723    l=(currRing->N);
1724    while ((l>0) && (p_GetExp( hh->m[k],l,tmpR)*pGetExp(delVar,l)==0)) l--;
1725    if (l==0)
1726    {
1727      j++;
1728      if (j >= IDELEMS(h3))
1729      {
1730        pEnlargeSet(&(h3->m),IDELEMS(h3),16);
1731        IDELEMS(h3) += 16;
1732      }
1733      h3->m[j] = prMoveR( hh->m[k], tmpR,origR);
1734      hh->m[k] = NULL;
1735    }
1736  }
1737  id_Delete(&hh, tmpR);
1738  idSkipZeroes(h3);
1739  rDelete(tmpR);
1740  if (w!=NULL)
1741    delete w;
1742  return h3;
1743}
1744
1745/*2
1746* compute the which-th ar-minor of the matrix a
1747*/
1748poly idMinor(matrix a, int ar, unsigned long which, ideal R)
1749{
1750  int     i,j,k,size;
1751  unsigned long curr;
1752  int *rowchoise,*colchoise;
1753  BOOLEAN rowch,colch;
1754  ideal result;
1755  matrix tmp;
1756  poly p,q;
1757
1758  i = binom(a->rows(),ar);
1759  j = binom(a->cols(),ar);
1760
1761  rowchoise=(int *)omAlloc(ar*sizeof(int));
1762  colchoise=(int *)omAlloc(ar*sizeof(int));
1763  if ((i>512) || (j>512) || (i*j >512)) size=512;
1764  else size=i*j;
1765  result=idInit(size,1);
1766  tmp=mpNew(ar,ar);
1767  k = 0; /* the index in result*/
1768  curr = 0; /* index of current minor */
1769  idInitChoise(ar,1,a->rows(),&rowch,rowchoise);
1770  while (!rowch)
1771  {
1772    idInitChoise(ar,1,a->cols(),&colch,colchoise);
1773    while (!colch)
1774    {
1775      if (curr == which)
1776      {
1777        for (i=1; i<=ar; i++)
1778        {
1779          for (j=1; j<=ar; j++)
1780          {
1781            MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
1782          }
1783        }
1784        p = mp_DetBareiss(tmp,currRing);
1785        if (p!=NULL)
1786        {
1787          if (R!=NULL)
1788          {
1789            q = p;
1790            p = kNF(R,currQuotient,q);
1791            pDelete(&q);
1792          }
1793          /*delete the matrix tmp*/
1794          for (i=1; i<=ar; i++)
1795          {
1796            for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
1797          }
1798          idDelete((ideal*)&tmp);
1799          omFreeSize((ADDRESS)rowchoise,ar*sizeof(int));
1800          omFreeSize((ADDRESS)colchoise,ar*sizeof(int));
1801          return (p);
1802        }
1803      }
1804      curr++;
1805      idGetNextChoise(ar,a->cols(),&colch,colchoise);
1806    }
1807    idGetNextChoise(ar,a->rows(),&rowch,rowchoise);
1808  }
1809  return (poly) 1;
1810}
1811
1812#ifdef WITH_OLD_MINOR
1813/*2
1814* compute all ar-minors of the matrix a
1815*/
1816ideal idMinors(matrix a, int ar, ideal R)
1817{
1818  int     i,j,k,size;
1819  int *rowchoise,*colchoise;
1820  BOOLEAN rowch,colch;
1821  ideal result;
1822  matrix tmp;
1823  poly p,q;
1824
1825  i = binom(a->rows(),ar);
1826  j = binom(a->cols(),ar);
1827
1828  rowchoise=(int *)omAlloc(ar*sizeof(int));
1829  colchoise=(int *)omAlloc(ar*sizeof(int));
1830  if ((i>512) || (j>512) || (i*j >512)) size=512;
1831  else size=i*j;
1832  result=idInit(size,1);
1833  tmp=mpNew(ar,ar);
1834  k = 0; /* the index in result*/
1835  idInitChoise(ar,1,a->rows(),&rowch,rowchoise);
1836  while (!rowch)
1837  {
1838    idInitChoise(ar,1,a->cols(),&colch,colchoise);
1839    while (!colch)
1840    {
1841      for (i=1; i<=ar; i++)
1842      {
1843        for (j=1; j<=ar; j++)
1844        {
1845          MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
1846        }
1847      }
1848      p = mp_DetBareiss(tmp,vcurrRing);
1849      if (p!=NULL)
1850      {
1851        if (R!=NULL)
1852        {
1853          q = p;
1854          p = kNF(R,currQuotient,q);
1855          pDelete(&q);
1856        }
1857        if (p!=NULL)
1858        {
1859          if (k>=size)
1860          {
1861            pEnlargeSet(&result->m,size,32);
1862            size += 32;
1863          }
1864          result->m[k] = p;
1865          k++;
1866        }
1867      }
1868      idGetNextChoise(ar,a->cols(),&colch,colchoise);
1869    }
1870    idGetNextChoise(ar,a->rows(),&rowch,rowchoise);
1871  }
1872  /*delete the matrix tmp*/
1873  for (i=1; i<=ar; i++)
1874  {
1875    for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
1876  }
1877  idDelete((ideal*)&tmp);
1878  if (k==0)
1879  {
1880    k=1;
1881    result->m[0]=NULL;
1882  }
1883  omFreeSize((ADDRESS)rowchoise,ar*sizeof(int));
1884  omFreeSize((ADDRESS)colchoise,ar*sizeof(int));
1885  pEnlargeSet(&result->m,size,k-size);
1886  IDELEMS(result) = k;
1887  return (result);
1888}
1889#else
1890/*2
1891* compute all ar-minors of the matrix a
1892* the caller of mpRecMin
1893* the elements of the result are not in R (if R!=NULL)
1894*/
1895ideal idMinors(matrix a, int ar, ideal R)
1896{
1897  int elems=0;
1898  int r=a->nrows,c=a->ncols;
1899  int i;
1900  matrix b;
1901  ideal result,h;
1902  ring origR;
1903  ring tmpR;
1904  long bound;
1905
1906  if((ar<=0) || (ar>r) || (ar>c))
1907  {
1908    Werror("%d-th minor, matrix is %dx%d",ar,r,c);
1909    return NULL;
1910  }
1911  h = idMatrix2Module(mp_Copy(a,currRing));
1912  bound = sm_ExpBound(h,c,r,ar,currRing);
1913  idDelete(&h);
1914  tmpR=sm_RingChange(origR,bound);
1915  b = mpNew(r,c);
1916  for (i=r*c-1;i>=0;i--)
1917  {
1918    if (a->m[i])
1919      b->m[i] = prCopyR(a->m[i],origR,currRing);
1920  }
1921  if (R!=NULL)
1922  {
1923    R = idrCopyR(R,origR,currRing);
1924    //if (ar>1) // otherwise done in mpMinorToResult
1925    //{
1926    //  matrix bb=(matrix)kNF(R,currQuotient,(ideal)b);
1927    //  bb->rank=b->rank; bb->nrows=b->nrows; bb->ncols=b->ncols;
1928    //  idDelete((ideal*)&b); b=bb;
1929    //}
1930  }
1931  result=idInit(32,1);
1932  if(ar>1) mpRecMin(ar-1,result,elems,b,r,c,NULL,R);
1933  else mpMinorToResult(result,elems,b,r,c,R);
1934  idDelete((ideal *)&b);
1935  if (R!=NULL) idDelete(&R);
1936  idSkipZeroes(result);
1937  rChangeCurrRing(origR);
1938  result = idrMoveR(result,tmpR,origR);
1939  smKillModifiedRing(tmpR);
1940  idTest(result);
1941  return result;
1942}
1943#endif
1944
1945/*2
1946*returns TRUE if id1 is a submodule of id2
1947*/
1948BOOLEAN idIsSubModule(ideal id1,ideal id2)
1949{
1950  int i;
1951  poly p;
1952
1953  if (idIs0(id1)) return TRUE;
1954  for (i=0;i<IDELEMS(id1);i++)
1955  {
1956    if (id1->m[i] != NULL)
1957    {
1958      p = kNF(id2,currQuotient,id1->m[i]);
1959      if (p != NULL)
1960      {
1961        pDelete(&p);
1962        return FALSE;
1963      }
1964    }
1965  }
1966  return TRUE;
1967}
1968
1969BOOLEAN idTestHomModule(ideal m, ideal Q, intvec *w)
1970{
1971  if ((Q!=NULL) && (!idHomIdeal(Q,NULL)))  { PrintS(" Q not hom\n"); return FALSE;}
1972  if (idIs0(m)) return TRUE;
1973
1974  int cmax=-1;
1975  int i;
1976  poly p=NULL;
1977  int length=IDELEMS(m);
1978  polyset P=m->m;
1979  for (i=length-1;i>=0;i--)
1980  {
1981    p=P[i];
1982    if (p!=NULL) cmax=si_max(cmax,(int)pMaxComp(p)+1);
1983  }
1984  if (w != NULL)
1985  if (w->length()+1 < cmax)
1986  {
1987    // Print("length: %d - %d \n", w->length(),cmax);
1988    return FALSE;
1989  }
1990
1991  if(w!=NULL)
1992    pSetModDeg(w);
1993
1994  for (i=length-1;i>=0;i--)
1995  {
1996    p=P[i];
1997    poly q=p;
1998    if (p!=NULL)
1999    {
2000      int d=currRing->pFDeg(p,currRing);
2001      loop
2002      {
2003        pIter(p);
2004        if (p==NULL) break;
2005        if (d!=currRing->pFDeg(p,currRing))
2006        {
2007          //pWrite(q); wrp(p); Print(" -> %d - %d\n",d,pFDeg(p,currRing));
2008          if(w!=NULL)
2009            pSetModDeg(NULL);
2010          return FALSE;
2011        }
2012      }
2013    }
2014  }
2015
2016  if(w!=NULL)
2017    pSetModDeg(NULL);
2018
2019  return TRUE;
2020}
2021
2022ideal idSeries(int n,ideal M,matrix U,intvec *w)
2023{
2024  for(int i=IDELEMS(M)-1;i>=0;i--)
2025  {
2026    if(U==NULL)
2027      M->m[i]=pSeries(n,M->m[i],NULL,w);
2028    else
2029    {
2030      M->m[i]=pSeries(n,M->m[i],MATELEM(U,i+1,i+1),w);
2031      MATELEM(U,i+1,i+1)=NULL;
2032    }
2033  }
2034  if(U!=NULL)
2035    idDelete((ideal*)&U);
2036  return M;
2037}
2038
2039matrix idDiff(matrix i, int k)
2040{
2041  int e=MATCOLS(i)*MATROWS(i);
2042  matrix r=mpNew(MATROWS(i),MATCOLS(i));
2043  r->rank=i->rank;
2044  int j;
2045  for(j=0; j<e; j++)
2046  {
2047    r->m[j]=pDiff(i->m[j],k);
2048  }
2049  return r;
2050}
2051
2052matrix idDiffOp(ideal I, ideal J,BOOLEAN multiply)
2053{
2054  matrix r=mpNew(IDELEMS(I),IDELEMS(J));
2055  int i,j;
2056  for(i=0; i<IDELEMS(I); i++)
2057  {
2058    for(j=0; j<IDELEMS(J); j++)
2059    {
2060      MATELEM(r,i+1,j+1)=pDiffOp(I->m[i],J->m[j],multiply);
2061    }
2062  }
2063  return r;
2064}
2065
2066/*3
2067*handles for some ideal operations the ring/syzcomp managment
2068*returns all syzygies (componentwise-)shifted by -syzcomp
2069*or -syzcomp-1 (in case of ideals as input)
2070static ideal idHandleIdealOp(ideal arg,int syzcomp,int isIdeal=FALSE)
2071{
2072  ring orig_ring=currRing;
2073  ring syz_ring=rCurrRingAssure_SyzComp();
2074  rSetSyzComp(length,syz_ring);
2075
2076  ideal s_temp;
2077  if (orig_ring!=syz_ring)
2078    s_temp=idrMoveR_NoSort(arg,orig_ring, syz_ring);
2079  else
2080    s_temp=arg;
2081
2082  ideal s_temp1 = kStd(s_temp,currQuotient,testHomog,&w,NULL,length);
2083  if (w!=NULL) delete w;
2084
2085  if (syz_ring!=orig_ring)
2086  {
2087    idDelete(&s_temp);
2088    rChangeCurrRing(orig_ring);
2089  }
2090
2091  idDelete(&temp);
2092  ideal temp1=idRingCopy(s_temp1,syz_ring);
2093
2094  if (syz_ring!=orig_ring)
2095  {
2096    rChangeCurrRing(syz_ring);
2097    idDelete(&s_temp1);
2098    rChangeCurrRing(orig_ring);
2099    rKill(syz_ring);
2100  }
2101
2102  for (i=0;i<IDELEMS(temp1);i++)
2103  {
2104    if ((temp1->m[i]!=NULL)
2105    && (pGetComp(temp1->m[i])<=length))
2106    {
2107      pDelete(&(temp1->m[i]));
2108    }
2109    else
2110    {
2111      p_Shift(&(temp1->m[i]),-length,currRing);
2112    }
2113  }
2114  temp1->rank = rk;
2115  idSkipZeroes(temp1);
2116
2117  return temp1;
2118}
2119*/
2120/*2
2121* represents (h1+h2)/h2=h1/(h1 intersect h2)
2122*/
2123//ideal idModulo (ideal h2,ideal h1)
2124ideal idModulo (ideal h2,ideal h1, tHomog hom, intvec ** w)
2125{
2126  intvec *wtmp=NULL;
2127
2128  int i,j,k,rk,flength=0,slength,length;
2129  poly p,q;
2130
2131  if (idIs0(h2))
2132    return idFreeModule(si_max(1,h2->ncols));
2133  if (!idIs0(h1))
2134    flength = id_RankFreeModule(h1,currRing);
2135  slength = id_RankFreeModule(h2,currRing);
2136  length  = si_max(flength,slength);
2137  if (length==0)
2138  {
2139    length = 1;
2140  }
2141  ideal temp = idInit(IDELEMS(h2),length+IDELEMS(h2));
2142  if ((w!=NULL)&&((*w)!=NULL))
2143  {
2144    //Print("input weights:");(*w)->show(1);PrintLn();
2145    int d;
2146    int k;
2147    wtmp=new intvec(length+IDELEMS(h2));
2148    for (i=0;i<length;i++)
2149      ((*wtmp)[i])=(**w)[i];
2150    for (i=0;i<IDELEMS(h2);i++)
2151    {
2152      poly p=h2->m[i];
2153      if (p!=NULL)
2154      {
2155        d = pDeg(p);
2156        k= pGetComp(p);
2157        if (slength>0) k--;
2158        d +=((**w)[k]);
2159        ((*wtmp)[i+length]) = d;
2160      }
2161    }
2162    //Print("weights:");wtmp->show(1);PrintLn();
2163  }
2164  for (i=0;i<IDELEMS(h2);i++)
2165  {
2166    temp->m[i] = pCopy(h2->m[i]);
2167    q = pOne();
2168    pSetComp(q,i+1+length);
2169    pSetmComp(q);
2170    if(temp->m[i]!=NULL)
2171    {
2172      if (slength==0) p_Shift(&(temp->m[i]),1,currRing);
2173      p = temp->m[i];
2174      while (pNext(p)!=NULL) pIter(p);
2175      pNext(p) = q;
2176    }
2177    else
2178      temp->m[i]=q;
2179  }
2180  rk = k = IDELEMS(h2);
2181  if (!idIs0(h1))
2182  {
2183    pEnlargeSet(&(temp->m),IDELEMS(temp),IDELEMS(h1));
2184    IDELEMS(temp) += IDELEMS(h1);
2185    for (i=0;i<IDELEMS(h1);i++)
2186    {
2187      if (h1->m[i]!=NULL)
2188      {
2189        temp->m[k] = pCopy(h1->m[i]);
2190        if (flength==0) p_Shift(&(temp->m[k]),1,currRing);
2191        k++;
2192      }
2193    }
2194  }
2195
2196  ring orig_ring=currRing;
2197  ring syz_ring=rAssure_SyzComp(orig_ring.TRUE);
2198  rChangeCurrRing(syz_ring);
2199  rSetSyzComp(length, syz_ring);
2200  ideal s_temp;
2201
2202  if (syz_ring != orig_ring)
2203  {
2204    s_temp = idrMoveR_NoSort(temp, orig_ring, syz_ring);
2205  }
2206  else
2207  {
2208    s_temp = temp;
2209  }
2210
2211  idTest(s_temp);
2212  ideal s_temp1 = kStd(s_temp,currQuotient,hom,&wtmp,NULL,length);
2213
2214  //if (wtmp!=NULL)  Print("output weights:");wtmp->show(1);PrintLn();
2215  if ((w!=NULL) && (*w !=NULL) && (wtmp!=NULL))
2216  {
2217    delete *w;
2218    *w=new intvec(IDELEMS(h2));
2219    for (i=0;i<IDELEMS(h2);i++)
2220      ((**w)[i])=(*wtmp)[i+length];
2221  }
2222  if (wtmp!=NULL) delete wtmp;
2223
2224  for (i=0;i<IDELEMS(s_temp1);i++)
2225  {
2226    if ((s_temp1->m[i]!=NULL)
2227    && (pGetComp(s_temp1->m[i])<=length))
2228    {
2229      pDelete(&(s_temp1->m[i]));
2230    }
2231    else
2232    {
2233      p_Shift(&(s_temp1->m[i]),-length,currRing);
2234    }
2235  }
2236  s_temp1->rank = rk;
2237  idSkipZeroes(s_temp1);
2238
2239  if (syz_ring!=orig_ring)
2240  {
2241    rChangeCurrRing(orig_ring);
2242    s_temp1 = idrMoveR_NoSort(s_temp1, syz_ring, orig_ring);
2243    rKill(syz_ring);
2244    // Hmm ... here seems to be a memory leak
2245    // However, simply deleting it causes memory trouble
2246    // idDelete(&s_temp);
2247  }
2248  else
2249  {
2250    idDelete(&temp);
2251  }
2252  idTest(s_temp1);
2253  return s_temp1;
2254}
2255
2256/*
2257*computes module-weights for liftings of homogeneous modules
2258*/
2259intvec * idMWLift(ideal mod,intvec * weights)
2260{
2261  if (idIs0(mod)) return new intvec(2);
2262  int i=IDELEMS(mod);
2263  while ((i>0) && (mod->m[i-1]==NULL)) i--;
2264  intvec *result = new intvec(i+1);
2265  while (i>0)
2266  {
2267    (*result)[i]=currRing->pFDeg(mod->m[i],currRing)+(*weights)[pGetComp(mod->m[i])];
2268  }
2269  return result;
2270}
2271
2272/*2
2273*sorts the kbase for idCoef* in a special way (lexicographically
2274*with x_max,...,x_1)
2275*/
2276ideal idCreateSpecialKbase(ideal kBase,intvec ** convert)
2277{
2278  int i;
2279  ideal result;
2280
2281  if (idIs0(kBase)) return NULL;
2282  result = idInit(IDELEMS(kBase),kBase->rank);
2283  *convert = idSort(kBase,FALSE);
2284  for (i=0;i<(*convert)->length();i++)
2285  {
2286    result->m[i] = pCopy(kBase->m[(**convert)[i]-1]);
2287  }
2288  return result;
2289}
2290
2291/*2
2292*returns the index of a given monom in the list of the special kbase
2293*/
2294int idIndexOfKBase(poly monom, ideal kbase)
2295{
2296  int j=IDELEMS(kbase);
2297
2298  while ((j>0) && (kbase->m[j-1]==NULL)) j--;
2299  if (j==0) return -1;
2300  int i=(currRing->N);
2301  while (i>0)
2302  {
2303    loop
2304    {
2305      if (pGetExp(monom,i)>pGetExp(kbase->m[j-1],i)) return -1;
2306      if (pGetExp(monom,i)==pGetExp(kbase->m[j-1],i)) break;
2307      j--;
2308      if (j==0) return -1;
2309    }
2310    if (i==1)
2311    {
2312      while(j>0)
2313      {
2314        if (pGetComp(monom)==pGetComp(kbase->m[j-1])) return j-1;
2315        if (pGetComp(monom)>pGetComp(kbase->m[j-1])) return -1;
2316        j--;
2317      }
2318    }
2319    i--;
2320  }
2321  return -1;
2322}
2323
2324/*2
2325*decomposes the monom in a part of coefficients described by the
2326*complement of how and a monom in variables occuring in how, the
2327*index of which in kbase is returned as integer pos (-1 if it don't
2328*exists)
2329*/
2330poly idDecompose(poly monom, poly how, ideal kbase, int * pos)
2331{
2332  int i;
2333  poly coeff=pOne(), base=pOne();
2334
2335  for (i=1;i<=(currRing->N);i++)
2336  {
2337    if (pGetExp(how,i)>0)
2338    {
2339      pSetExp(base,i,pGetExp(monom,i));
2340    }
2341    else
2342    {
2343      pSetExp(coeff,i,pGetExp(monom,i));
2344    }
2345  }
2346  pSetComp(base,pGetComp(monom));
2347  pSetm(base);
2348  pSetCoeff(coeff,nCopy(pGetCoeff(monom)));
2349  pSetm(coeff);
2350  *pos = idIndexOfKBase(base,kbase);
2351  if (*pos<0)
2352    pDelete(&coeff);
2353  pDelete(&base);
2354  return coeff;
2355}
2356
2357/*2
2358*returns a matrix A of coefficients with kbase*A=arg
2359*if all monomials in variables of how occur in kbase
2360*the other are deleted
2361*/
2362matrix idCoeffOfKBase(ideal arg, ideal kbase, poly how)
2363{
2364  matrix result;
2365  ideal tempKbase;
2366  poly p,q;
2367  intvec * convert;
2368  int i=IDELEMS(kbase),j=IDELEMS(arg),k,pos;
2369#if 0
2370  while ((i>0) && (kbase->m[i-1]==NULL)) i--;
2371  if (idIs0(arg))
2372    return mpNew(i,1);
2373  while ((j>0) && (arg->m[j-1]==NULL)) j--;
2374  result = mpNew(i,j);
2375#else
2376  result = mpNew(i, j);
2377  while ((j>0) && (arg->m[j-1]==NULL)) j--;
2378#endif
2379
2380  tempKbase = idCreateSpecialKbase(kbase,&convert);
2381  for (k=0;k<j;k++)
2382  {
2383    p = arg->m[k];
2384    while (p!=NULL)
2385    {
2386      q = idDecompose(p,how,tempKbase,&pos);
2387      if (pos>=0)
2388      {
2389        MATELEM(result,(*convert)[pos],k+1) =
2390            pAdd(MATELEM(result,(*convert)[pos],k+1),q);
2391      }
2392      else
2393        pDelete(&q);
2394      pIter(p);
2395    }
2396  }
2397  idDelete(&tempKbase);
2398  return result;
2399}
2400
2401static void idDeleteComps(ideal arg,int* red_comp,int del)
2402// red_comp is an array [0..args->rank]
2403{
2404  int i,j;
2405  poly p;
2406
2407  for (i=IDELEMS(arg)-1;i>=0;i--)
2408  {
2409    p = arg->m[i];
2410    while (p!=NULL)
2411    {
2412      j = pGetComp(p);
2413      if (red_comp[j]!=j)
2414      {
2415        pSetComp(p,red_comp[j]);
2416        pSetmComp(p);
2417      }
2418      pIter(p);
2419    }
2420  }
2421  (arg->rank) -= del;
2422}
2423
2424/*2
2425* returns the presentation of an isomorphic, minimally
2426* embedded  module (arg represents the quotient!)
2427*/
2428ideal idMinEmbedding(ideal arg,BOOLEAN inPlace, intvec **w)
2429{
2430  if (idIs0(arg)) return idInit(1,arg->rank);
2431  int i,next_gen,next_comp;
2432  ideal res=arg;
2433  if (!inPlace) res = idCopy(arg);
2434  res->rank=si_max(res->rank,id_RankFreeModule(res,currRing));
2435  int *red_comp=(int*)omAlloc((res->rank+1)*sizeof(int));
2436  for (i=res->rank;i>=0;i--) red_comp[i]=i;
2437
2438  int del=0;
2439  loop
2440  {
2441    next_gen = idReadOutPivot(res,&next_comp);
2442    if (next_gen<0) break;
2443    del++;
2444    syGaussForOne(res,next_gen,next_comp,0,IDELEMS(res));
2445    for(i=next_comp+1;i<=arg->rank;i++) red_comp[i]--;
2446    if ((w !=NULL)&&(*w!=NULL))
2447    {
2448      for(i=next_comp;i<(*w)->length();i++) (**w)[i-1]=(**w)[i];
2449    }
2450  }
2451
2452  idDeleteComps(res,red_comp,del);
2453  idSkipZeroes(res);
2454  omFree(red_comp);
2455
2456  if ((w !=NULL)&&(*w!=NULL) &&(del>0))
2457  {
2458    intvec *wtmp=new intvec((*w)->length()-del);
2459    for(i=0;i<res->rank;i++) (*wtmp)[i]=(**w)[i];
2460    delete *w;
2461    *w=wtmp;
2462  }
2463  return res;
2464}
2465
2466#include <polys/clapsing.h>
2467
2468#ifdef HAVE_FACTORY
2469poly id_GCD(poly f, poly g, const ring r)
2470{
2471  ring save_r=currRing;
2472  rChangeCurrRing(r);
2473  ideal I=idInit(2,1); I->m[0]=f; I->m[1]=g;
2474  intvec *w = NULL;
2475  ideal S=idSyzygies(I,testHomog,&w);
2476  if (w!=NULL) delete w;
2477  poly gg=pTakeOutComp(&(S->m[0]),2);
2478  idDelete(&S);
2479  poly gcd_p=singclap_pdivide(f,gg,r);
2480  pDelete(&gg);
2481  rChangeCurrRing(save_r);
2482  return gcd_p;
2483}
2484#endif
2485
2486/*2
2487* xx,q: arrays of length 0..rl-1
2488* xx[i]: SB mod q[i]
2489* assume: char=0
2490* assume: q[i]!=0
2491* destroys xx
2492*/
2493#ifdef HAVE_FACTORY
2494ideal idChineseRemainder(ideal *xx, number *q, int rl)
2495{
2496  int cnt=IDELEMS(xx[0])*xx[0]->nrows;
2497  ideal result=idInit(cnt,xx[0]->rank);
2498  result->nrows=xx[0]->nrows; // for lifting matrices
2499  result->ncols=xx[0]->ncols; // for lifting matrices
2500  int i,j;
2501  poly r,h,hh,res_p;
2502  number *x=(number *)omAlloc(rl*sizeof(number));
2503  for(i=cnt-1;i>=0;i--)
2504  {
2505    res_p=NULL;
2506    loop
2507    {
2508      r=NULL;
2509      for(j=rl-1;j>=0;j--)
2510      {
2511        h=xx[j]->m[i];
2512        if ((h!=NULL)
2513        &&((r==NULL)||(pLmCmp(r,h)==-1)))
2514          r=h;
2515      }
2516      if (r==NULL) break;
2517      h=pHead(r);
2518      for(j=rl-1;j>=0;j--)
2519      {
2520        hh=xx[j]->m[i];
2521        if ((hh!=NULL) && (pLmCmp(r,hh)==0))
2522        {
2523          x[j]=pGetCoeff(hh);
2524          hh=pLmFreeAndNext(hh);
2525          xx[j]->m[i]=hh;
2526        }
2527        else
2528          x[j]=nlInit(0, currRing);
2529      }
2530      number n=nlChineseRemainder(x,q,rl);
2531      for(j=rl-1;j>=0;j--)
2532      {
2533        x[j]=NULL; // nlInit(0...) takes no memory
2534      }
2535      if (nlIsZero(n)) pDelete(&h);
2536      else
2537      {
2538        pSetCoeff(h,n);
2539        //Print("new mon:");pWrite(h);
2540        res_p=pAdd(res_p,h);
2541      }
2542    }
2543    result->m[i]=res_p;
2544  }
2545  omFree(x);
2546  for(i=rl-1;i>=0;i--) idDelete(&(xx[i]));
2547  omFree(xx);
2548  return result;
2549}
2550#endif
2551/* currently unsed:
2552ideal idChineseRemainder(ideal *xx, intvec *iv)
2553{
2554  int rl=iv->length();
2555  number *q=(number *)omAlloc(rl*sizeof(number));
2556  int i;
2557  for(i=0; i<rl; i++)
2558  {
2559    q[i]=nInit((*iv)[i]);
2560  }
2561  return idChineseRemainder(xx,q,rl);
2562}
2563*/
2564/*
2565 * lift ideal with coeffs over Z (mod N) to Q via Farey
2566 */
2567ideal idFarey(ideal x, number N)
2568{
2569  int cnt=IDELEMS(x)*x->nrows;
2570  ideal result=idInit(cnt,x->rank);
2571  result->nrows=x->nrows; // for lifting matrices
2572  result->ncols=x->ncols; // for lifting matrices
2573
2574  int i;
2575  for(i=cnt-1;i>=0;i--)
2576  {
2577    poly h=pCopy(x->m[i]);
2578    result->m[i]=h;
2579    while(h!=NULL)
2580    {
2581      number c=pGetCoeff(h);
2582      pSetCoeff0(h,nlFarey(c,N));
2583      nDelete(&c);
2584      pIter(h);
2585    }
2586    while((result->m[i]!=NULL)&&(nIsZero(pGetCoeff(result->m[i]))))
2587    {
2588      pLmDelete(&(result->m[i]));
2589    }
2590    h=result->m[i];
2591    while((h!=NULL) && (pNext(h)!=NULL))
2592    {
2593      if(nIsZero(pGetCoeff(pNext(h))))
2594      {
2595        pLmDelete(&pNext(h));
2596      }
2597      else pIter(h);
2598    }
2599  }
2600  return result;
2601}
Note: See TracBrowser for help on using the repository browser.