source: git/kernel/ideals.cc @ 99f2bc

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