source: git/kernel/ideals.cc @ f7d39b

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