source: git/kernel/ideals.cc @ 0e8683

spielwiese
Last change on this file since 0e8683 was 85e68dd, checked in by Hans Schönemann <hannes@…>, 16 years ago
*hannes: gcc 4.2 git-svn-id: file:///usr/local/Singular/svn/trunk@10634 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.50 2008-03-19 17:44:08 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||w!=NULL&&pDegW(p0,w)>n)
1901            pDelete(&p0);
1902          else
1903            R->m[i]=pAdd(R->m[i],p0);
1904          j=IDELEMS(Q)-1;
1905        }
1906        else
1907          j--;
1908      }
1909    }
1910  }
1911}
1912
1913/*2
1914*computes the quotient of h1,h2 : internal routine for idQuot
1915*BEWARE: the returned ideals may contain incorrectly ordered polys !
1916*
1917*/
1918static ideal idInitializeQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb,
1919                               BOOLEAN *addOnlyOne, int *kkmax)
1920{
1921  ideal temph1;
1922  poly     p,q = NULL;
1923  int i,l,ll,k,kkk,kmax;
1924  int j = 0;
1925  int k1 = idRankFreeModule(h1);
1926  int k2 = idRankFreeModule(h2);
1927  tHomog   hom=isNotHomog;
1928
1929  k=si_max(k1,k2);
1930  if (k==0)
1931    k = 1;
1932  if ((k2==0) && (k>1)) *addOnlyOne = FALSE;
1933
1934  intvec * weights;
1935  hom = (tHomog)idHomModule(h1,currQuotient,&weights);
1936  if (addOnlyOne && (!h1IsStb))
1937    temph1 = kStd(h1,currQuotient,hom,&weights,NULL);
1938  else
1939    temph1 = idCopy(h1);
1940  if (weights!=NULL) delete weights;
1941  idTest(temph1);
1942/*--- making a single vector from h2 ---------------------*/
1943  for (i=0; i<IDELEMS(h2); i++)
1944  {
1945    if (h2->m[i] != NULL)
1946    {
1947      p = pCopy(h2->m[i]);
1948      if (k2 == 0)
1949        pShift(&p,j*k+1);
1950      else
1951        pShift(&p,j*k);
1952      q = pAdd(q,p);
1953      j++;
1954    }
1955  }
1956  *kkmax = kmax = j*k+1;
1957/*--- adding a monomial for the result (syzygy) ----------*/
1958  p = q;
1959  while (pNext(p)!=NULL) pIter(p);
1960  pNext(p) = pOne();
1961  pIter(p);
1962  pSetComp(p,kmax);
1963  pSetmComp(p);
1964/*--- constructing the big matrix ------------------------*/
1965  ideal h4 = idInit(16,kmax+k-1);
1966  h4->m[0] = q;
1967  if (k2 == 0)
1968  {
1969    if (k > IDELEMS(h4))
1970    {
1971      pEnlargeSet(&(h4->m),IDELEMS(h4),k-IDELEMS(h4));
1972      IDELEMS(h4) = k;
1973    }
1974    for (i=1; i<k; i++)
1975    {
1976      p = pCopy_noCheck(h4->m[i-1]);
1977      pShift(&p,1);
1978      h4->m[i] = p;
1979    }
1980  }
1981
1982  kkk = IDELEMS(h4);
1983  i = IDELEMS(temph1);
1984  while ((i>0) && (temph1->m[i-1]==NULL)) i--;
1985  for (l=0; l<i; l++)
1986  {
1987    if(temph1->m[l]!=NULL)
1988    {
1989      for (ll=0; ll<j; ll++)
1990      {
1991        p = pCopy(temph1->m[l]);
1992        if (k1 == 0)
1993          pShift(&p,ll*k+1);
1994        else
1995          pShift(&p,ll*k);
1996        if (kkk >= IDELEMS(h4))
1997        {
1998          pEnlargeSet(&(h4->m),IDELEMS(h4),16);
1999          IDELEMS(h4) += 16;
2000        }
2001        h4->m[kkk] = p;
2002        kkk++;
2003      }
2004    }
2005  }
2006/*--- if h2 goes in as single vector - the h1-part is just SB ---*/
2007  if (*addOnlyOne)
2008  {
2009    p = h4->m[0];
2010    for (i=0;i<IDELEMS(h4)-1;i++)
2011    {
2012      h4->m[i] = h4->m[i+1];
2013    }
2014    h4->m[IDELEMS(h4)-1] = p;
2015    idSkipZeroes(h4);
2016    test |= Sy_bit(OPT_SB_1);
2017  }
2018  idDelete(&temph1);
2019  return h4;
2020}
2021/*2
2022*computes the quotient of h1,h2
2023*/
2024ideal idQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb, BOOLEAN resultIsIdeal)
2025{
2026  // first check for special case h1:(0)
2027  if (idIs0(h2))
2028  {
2029    ideal res;
2030    if (resultIsIdeal)
2031    {
2032      res = idInit(1,1);
2033      res->m[0] = pOne();
2034    }
2035    else
2036      res = idFreeModule(h1->rank);
2037    return res;
2038  }
2039  BITSET old_test=test;
2040  poly     p,q = NULL;
2041  int i,l,ll,k,kkk,kmax;
2042  BOOLEAN  addOnlyOne=TRUE;
2043  tHomog   hom=isNotHomog;
2044  intvec * weights1;
2045
2046  ideal s_h4 = idInitializeQuot (h1,h2,h1IsStb,&addOnlyOne,&kmax);
2047
2048  hom = (tHomog)idHomModule(s_h4,currQuotient,&weights1);
2049
2050  ring orig_ring=currRing;
2051  ring syz_ring=rCurrRingAssure_SyzComp();
2052  rSetSyzComp(kmax-1);
2053  if (orig_ring!=syz_ring)
2054  //  s_h4 = idrMoveR_NoSort(s_h4,orig_ring);
2055    s_h4 = idrMoveR(s_h4,orig_ring);
2056  idTest(s_h4);
2057  ideal s_h3;
2058  if (addOnlyOne)
2059  {
2060    s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,kmax-1,IDELEMS(s_h4)-1);
2061  }
2062  else
2063  {
2064    s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,kmax-1);
2065  }
2066  idTest(s_h3);
2067  if (weights1!=NULL) delete weights1;
2068  idDelete(&s_h4);
2069
2070
2071  for (i=0;i<IDELEMS(s_h3);i++)
2072  {
2073    if ((s_h3->m[i]!=NULL) && (pGetComp(s_h3->m[i])>=kmax))
2074    {
2075      if (resultIsIdeal)
2076        pShift(&s_h3->m[i],-kmax);
2077      else
2078        pShift(&s_h3->m[i],-kmax+1);
2079    }
2080    else
2081      pDelete(&s_h3->m[i]);
2082  }
2083  if (resultIsIdeal)
2084    s_h3->rank = 1;
2085  else
2086    s_h3->rank = h1->rank;
2087  if(syz_ring!=orig_ring)
2088  {
2089//    pDelete(&q);
2090    rChangeCurrRing(orig_ring);
2091    s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
2092    rKill(syz_ring);
2093  }
2094  idSkipZeroes(s_h3);
2095  test = old_test;
2096  idTest(s_h3);
2097  return s_h3;
2098}
2099
2100/*2
2101*computes recursively all monomials of a certain degree
2102*in every step the actvar-th entry in the exponential
2103*vector is incremented and the other variables are
2104*computed by recursive calls of makemonoms
2105*if the last variable is reached, the difference to the
2106*degree is computed directly
2107*vars is the number variables
2108*actvar is the actual variable to handle
2109*deg is the degree of the monomials to compute
2110*monomdeg is the actual degree of the monomial in consideration
2111*/
2112static void makemonoms(int vars,int actvar,int deg,int monomdeg)
2113{
2114  poly p;
2115  int i=0;
2116
2117  if ((idpowerpoint == 0) && (actvar ==1))
2118  {
2119    idpower[idpowerpoint] = pOne();
2120    monomdeg = 0;
2121  }
2122  while (i<=deg)
2123  {
2124    if (deg == monomdeg)
2125    {
2126      pSetm(idpower[idpowerpoint]);
2127      idpowerpoint++;
2128      return;
2129    }
2130    if (actvar == vars)
2131    {
2132      pSetExp(idpower[idpowerpoint],actvar,deg-monomdeg);
2133      pSetm(idpower[idpowerpoint]);
2134      pTest(idpower[idpowerpoint]);
2135      idpowerpoint++;
2136      return;
2137    }
2138    else
2139    {
2140      p = pCopy(idpower[idpowerpoint]);
2141      makemonoms(vars,actvar+1,deg,monomdeg);
2142      idpower[idpowerpoint] = p;
2143    }
2144    monomdeg++;
2145    pSetExp(idpower[idpowerpoint],actvar,pGetExp(idpower[idpowerpoint],actvar)+1);
2146    pSetm(idpower[idpowerpoint]);
2147    pTest(idpower[idpowerpoint]);
2148    i++;
2149  }
2150}
2151
2152/*2
2153*returns the deg-th power of the maximal ideal of 0
2154*/
2155ideal idMaxIdeal(int deg)
2156{
2157  if (deg < 0)
2158  {
2159    WarnS("maxideal: power must be non-negative");
2160  }
2161  if (deg < 1)
2162  {
2163    ideal I=idInit(1,1);
2164    I->m[0]=pOne();
2165    return I;
2166  }
2167  if (deg == 1)
2168  {
2169    return idMaxIdeal();
2170  }
2171
2172  int vars = currRing->N;
2173  int i = binom(vars+deg-1,deg);
2174  if (i<=0) return idInit(1,1);
2175  ideal id=idInit(i,1);
2176  idpower = id->m;
2177  idpowerpoint = 0;
2178  makemonoms(vars,1,deg,0);
2179  idpower = NULL;
2180  idpowerpoint = 0;
2181  return id;
2182}
2183
2184/*2
2185*computes recursively all generators of a certain degree
2186*of the ideal "givenideal"
2187*elms is the number elements in the given ideal
2188*actelm is the actual element to handle
2189*deg is the degree of the power to compute
2190*gendeg is the actual degree of the generator in consideration
2191*/
2192static void makepotence(int elms,int actelm,int deg,int gendeg)
2193{
2194  poly p;
2195  int i=0;
2196
2197  if ((idpowerpoint == 0) && (actelm ==1))
2198  {
2199    idpower[idpowerpoint] = pOne();
2200    gendeg = 0;
2201  }
2202  while (i<=deg)
2203  {
2204    if (deg == gendeg)
2205    {
2206      idpowerpoint++;
2207      return;
2208    }
2209    if (actelm == elms)
2210    {
2211      p=pPower(pCopy(givenideal[actelm-1]),deg-gendeg);
2212      idpower[idpowerpoint]=pMult(idpower[idpowerpoint],p);
2213      idpowerpoint++;
2214      return;
2215    }
2216    else
2217    {
2218      p = pCopy(idpower[idpowerpoint]);
2219      makepotence(elms,actelm+1,deg,gendeg);
2220      idpower[idpowerpoint] = p;
2221    }
2222    gendeg++;
2223    idpower[idpowerpoint]=pMult(idpower[idpowerpoint],pCopy(givenideal[actelm-1]));
2224    i++;
2225  }
2226}
2227
2228/*2
2229*returns the deg-th power of the ideal gid
2230*/
2231//ideal idPower(ideal gid,int deg)
2232//{
2233//  int i;
2234//  ideal id;
2235//
2236//  if (deg < 1) deg = 1;
2237//  i = binom(IDELEMS(gid)+deg-1,deg);
2238//  id=idInit(i,1);
2239//  idpower = id->m;
2240//  givenideal = gid->m;
2241//  idpowerpoint = 0;
2242//  makepotence(IDELEMS(gid),1,deg,0);
2243//  idpower = NULL;
2244//  givenideal = NULL;
2245//  idpowerpoint = 0;
2246//  return id;
2247//}
2248static void idNextPotence(ideal given, ideal result,
2249  int begin, int end, int deg, int restdeg, poly ap)
2250{
2251  poly p;
2252  int i;
2253
2254  p = pPower(pCopy(given->m[begin]),restdeg);
2255  i = result->nrows;
2256  result->m[i] = pMult(pCopy(ap),p);
2257//PrintS(".");
2258  (result->nrows)++;
2259  if (result->nrows >= IDELEMS(result))
2260  {
2261    pEnlargeSet(&(result->m),IDELEMS(result),16);
2262    IDELEMS(result) += 16;
2263  }
2264  if (begin == end) return;
2265  for (i=restdeg-1;i>0;i--)
2266  {
2267    p = pPower(pCopy(given->m[begin]),i);
2268    p = pMult(pCopy(ap),p);
2269    idNextPotence(given, result, begin+1, end, deg, restdeg-i, p);
2270    pDelete(&p);
2271  }
2272  idNextPotence(given, result, begin+1, end, deg, restdeg, ap);
2273}
2274
2275ideal idPower(ideal given,int exp)
2276{
2277  ideal result,temp;
2278  poly p1;
2279  int i;
2280
2281  if (idIs0(given)) return idInit(1,1);
2282  temp = idCopy(given);
2283  idSkipZeroes(temp);
2284  i = binom(IDELEMS(temp)+exp-1,exp);
2285  result = idInit(i,1);
2286  result->nrows = 0;
2287//Print("ideal contains %d elements\n",i);
2288  p1=pOne();
2289  idNextPotence(temp,result,0,IDELEMS(temp)-1,exp,exp,p1);
2290  pDelete(&p1);
2291  idDelete(&temp);
2292  result->nrows = 1;
2293  idDelEquals(result);
2294  idSkipZeroes(result);
2295  return result;
2296}
2297
2298/*2
2299* eliminate delVar (product of vars) in h1
2300*/
2301ideal idElimination (ideal h1,poly delVar,intvec *hilb)
2302{
2303  int    i,j=0,k,l;
2304  ideal  h,hh, h3;
2305  int    *ord,*block0,*block1;
2306  int    ordersize=2;
2307  int    **wv;
2308  tHomog hom;
2309  intvec * w;
2310  ring tmpR;
2311  ring origR = currRing;
2312
2313  if (delVar==NULL)
2314  {
2315    return idCopy(h1);
2316  }
2317  if (currQuotient!=NULL)
2318  {
2319    WerrorS("cannot eliminate in a qring");
2320    return idCopy(h1);
2321  }
2322  if (idIs0(h1)) return idInit(1,h1->rank);
2323#ifdef HAVE_PLURAL
2324  if (rIsPluralRing(currRing))
2325    /* in the NC case, we have to check the admissibility of */
2326    /* the subalgebra to be intersected with */
2327  {
2328    if (ncRingType(currRing)!=nc_skew) /* in (quasi)-commutative algebras every subalgebra is admissible */
2329    {
2330      if (nc_CheckSubalgebra(delVar,currRing))
2331      {
2332        WerrorS("no elimination is possible: subalgebra is not admissible");
2333        return idCopy(h1);
2334      }
2335    }
2336  }
2337#endif
2338  hom=(tHomog)idHomModule(h1,NULL,&w); //sets w to weight vector or NULL
2339  h3=idInit(16,h1->rank);
2340  for (k=0;; k++)
2341  {
2342    if (origR->order[k]!=0) ordersize++;
2343    else break;
2344  }
2345  ord=(int*)omAlloc0(ordersize*sizeof(int));
2346  block0=(int*)omAlloc(ordersize*sizeof(int));
2347  block1=(int*)omAlloc(ordersize*sizeof(int));
2348  wv=(int**) omAlloc0(ordersize*sizeof(int**));
2349  for (k=0;k<ordersize-1; k++)
2350  {
2351    block0[k+1] = origR->block0[k];
2352    block1[k+1] = origR->block1[k];
2353    ord[k+1] = origR->order[k];
2354    if (origR->wvhdl[k]!=NULL) wv[k+1] = (int*) omMemDup(origR->wvhdl[k]);
2355  }
2356  block0[0] = 1;
2357  block1[0] = rVar(origR);
2358  wv[0]=(int*)omAlloc((pVariables+1)*sizeof(int));
2359  memset(wv[0],0,(pVariables+1)*sizeof(int));
2360  for (j=0;j<rVar(origR);j++)
2361    if (pGetExp(delVar,j+1)!=0) wv[0][j]=1;
2362  // use this special ordering: like ringorder_a, except that pFDeg, pWeights
2363  // ignore it
2364  ord[0] = ringorder_aa;
2365
2366  // fill in tmp ring to get back the data later on
2367  tmpR  = rCopy0(origR,FALSE,FALSE); // qring==NULL
2368  tmpR->order = ord;
2369  tmpR->block0 = block0;
2370  tmpR->block1 = block1;
2371  tmpR->wvhdl = wv;
2372  rComplete(tmpR, 1);
2373
2374#ifdef HAVE_PLURAL
2375  /* update nc structure on tmpR */
2376  if (rIsPluralRing(currRing))
2377  {
2378    BOOLEAN bBAD = FALSE;
2379    if ( nc_rComplete(origR, tmpR) )
2380    {
2381      Werror("error in nc_rComplete");
2382      bBAD = TRUE;
2383    }
2384    if (!bBAD)
2385    {
2386      /* tests the admissibility of the new elim. ordering */
2387      if ( nc_CheckOrdCondition( tmpR->nc->D, tmpR) )
2388      {
2389        Werror("no elimination is possible: ordering condition is violated");
2390        bBAD = TRUE;
2391      }
2392    }
2393    if (bBAD)
2394    {
2395      // cleanup
2396      rDelete(tmpR);
2397      if (w!=NULL)
2398      {
2399        delete w;
2400      }
2401      return idCopy(h1);
2402    }
2403  }
2404#endif
2405  // change into the new ring
2406  //pChangeRing(pVariables,currRing->OrdSgn,ord,block0,block1,wv);
2407  rChangeCurrRing(tmpR);
2408
2409  h = idInit(IDELEMS(h1),h1->rank);
2410  // fetch data from the old ring
2411  for (k=0;k<IDELEMS(h1);k++) h->m[k] = prCopyR( h1->m[k], origR);
2412  // compute kStd
2413#if 1
2414  hh = kStd(h,NULL,hom,&w,hilb);
2415  idDelete(&h);
2416#else
2417  extern ideal kGroebner(ideal F, ideal Q);
2418  hh=kGroebner(h,NULL);
2419#endif
2420  // go back to the original ring
2421  rChangeCurrRing(origR);
2422  i = IDELEMS(hh)-1;
2423  while ((i >= 0) && (hh->m[i] == NULL)) i--;
2424  j = -1;
2425  // fetch data from temp ring
2426  for (k=0; k<=i; k++)
2427  {
2428    l=pVariables;
2429    while ((l>0) && (p_GetExp( hh->m[k],l,tmpR)*pGetExp(delVar,l)==0)) l--;
2430    if (l==0)
2431    {
2432      j++;
2433      if (j >= IDELEMS(h3))
2434      {
2435        pEnlargeSet(&(h3->m),IDELEMS(h3),16);
2436        IDELEMS(h3) += 16;
2437      }
2438      h3->m[j] = prCopyR( hh->m[k], tmpR);
2439    }
2440  }
2441  id_Delete(&hh, tmpR);
2442  idSkipZeroes(h3);
2443  rDelete(tmpR);
2444  if (w!=NULL)
2445    delete w;
2446  return h3;
2447}
2448
2449#ifdef WITH_OLD_MINOR
2450/*2
2451* compute all ar-minors of the matrix a
2452*/
2453ideal idMinors(matrix a, int ar, ideal R)
2454{
2455  int     i,j,k,size;
2456  int *rowchoise,*colchoise;
2457  BOOLEAN rowch,colch;
2458  ideal result;
2459  matrix tmp;
2460  poly p,q;
2461
2462  i = binom(a->rows(),ar);
2463  j = binom(a->cols(),ar);
2464
2465  rowchoise=(int *)omAlloc(ar*sizeof(int));
2466  colchoise=(int *)omAlloc(ar*sizeof(int));
2467  if ((i>512) || (j>512) || (i*j >512)) size=512;
2468  else size=i*j;
2469  result=idInit(size,1);
2470  tmp=mpNew(ar,ar);
2471  k = 0; /* the index in result*/
2472  idInitChoise(ar,1,a->rows(),&rowch,rowchoise);
2473  while (!rowch)
2474  {
2475    idInitChoise(ar,1,a->cols(),&colch,colchoise);
2476    while (!colch)
2477    {
2478      for (i=1; i<=ar; i++)
2479      {
2480        for (j=1; j<=ar; j++)
2481        {
2482          MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
2483        }
2484      }
2485      p = mpDetBareiss(tmp);
2486      if (p!=NULL)
2487      {
2488        if (R!=NULL)
2489        {
2490          q = p;
2491          p = kNF(R,currQuotient,q);
2492          pDelete(&q);
2493        }
2494        if (p!=NULL)
2495        {
2496          if (k>=size)
2497          {
2498            pEnlargeSet(&result->m,size,32);
2499            size += 32;
2500          }
2501          result->m[k] = p;
2502          k++;
2503        }
2504      }
2505      idGetNextChoise(ar,a->cols(),&colch,colchoise);
2506    }
2507    idGetNextChoise(ar,a->rows(),&rowch,rowchoise);
2508  }
2509  /*delete the matrix tmp*/
2510  for (i=1; i<=ar; i++)
2511  {
2512    for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
2513  }
2514  idDelete((ideal*)&tmp);
2515  if (k==0)
2516  {
2517    k=1;
2518    result->m[0]=NULL;
2519  }
2520  omFreeSize((ADDRESS)rowchoise,ar*sizeof(int));
2521  omFreeSize((ADDRESS)colchoise,ar*sizeof(int));
2522  pEnlargeSet(&result->m,size,k-size);
2523  IDELEMS(result) = k;
2524  return (result);
2525}
2526#else
2527/*2
2528* compute all ar-minors of the matrix a
2529* the caller of mpRecMin
2530* the elements of the result are not in R (if R!=NULL)
2531*/
2532ideal idMinors(matrix a, int ar, ideal R)
2533{
2534  int elems=0;
2535  int r=a->nrows,c=a->ncols;
2536  int i;
2537  matrix b;
2538  ideal result,h;
2539  ring origR;
2540  sip_sring tmpR;
2541  Exponent_t bound;
2542
2543  if((ar<=0) || (ar>r) || (ar>c))
2544  {
2545    Werror("%d-th minor, matrix is %dx%d",ar,r,c);
2546    return NULL;
2547  }
2548  h = idMatrix2Module(mpCopy(a));
2549  bound = smExpBound(h,c,r,ar);
2550  idDelete(&h);
2551  smRingChange(&origR,tmpR,bound);
2552  b = mpNew(r,c);
2553  for (i=r*c-1;i>=0;i--)
2554  {
2555    if (a->m[i])
2556      b->m[i] = prCopyR(a->m[i],origR);
2557  }
2558  if (R) R = idrCopyR(R,origR);
2559  result=idInit(32,1);
2560  if(ar>1) mpRecMin(ar-1,result,elems,b,r,c,NULL,R);
2561  else mpMinorToResult(result,elems,b,r,c,R);
2562  idDelete((ideal *)&b);
2563  if (R) idDelete(&R);
2564  idSkipZeroes(result);
2565  rChangeCurrRing(origR);
2566  result = idrMoveR(result,&tmpR);
2567  smRingClean(origR,tmpR);
2568  idTest(result);
2569  return result;
2570}
2571#endif
2572
2573/*2
2574*skips all zeroes and double elements, searches also for units
2575*/
2576void idCompactify(ideal id)
2577{
2578  int i,j;
2579  BOOLEAN b=FALSE;
2580
2581  i = IDELEMS(id)-1;
2582  while ((! b) && (i>=0))
2583  {
2584    b=pIsUnit(id->m[i]);
2585    i--;
2586  }
2587  if (b)
2588  {
2589    for(i=IDELEMS(id)-1;i>=0;i--) pDelete(&id->m[i]);
2590    id->m[0]=pOne();
2591  }
2592  else
2593  {
2594    idDelMultiples(id);
2595  }
2596  idSkipZeroes(id);
2597}
2598
2599/*2
2600*returns TRUE if id1 is a submodule of id2
2601*/
2602BOOLEAN idIsSubModule(ideal id1,ideal id2)
2603{
2604  int i;
2605  poly p;
2606
2607  if (idIs0(id1)) return TRUE;
2608  for (i=0;i<IDELEMS(id1);i++)
2609  {
2610    if (id1->m[i] != NULL)
2611    {
2612      p = kNF(id2,currQuotient,id1->m[i]);
2613      if (p != NULL)
2614      {
2615        pDelete(&p);
2616        return FALSE;
2617      }
2618    }
2619  }
2620  return TRUE;
2621}
2622
2623/*2
2624* returns the ideals of initial terms
2625*/
2626ideal idHead(ideal h)
2627{
2628  ideal m = idInit(IDELEMS(h),h->rank);
2629  int i;
2630
2631  for (i=IDELEMS(h)-1;i>=0; i--)
2632  {
2633    if (h->m[i]!=NULL) m->m[i]=pHead(h->m[i]);
2634  }
2635  return m;
2636}
2637
2638ideal idHomogen(ideal h, int varnum)
2639{
2640  ideal m = idInit(IDELEMS(h),h->rank);
2641  int i;
2642
2643  for (i=IDELEMS(h)-1;i>=0; i--)
2644  {
2645    m->m[i]=pHomogen(h->m[i],varnum);
2646  }
2647  return m;
2648}
2649
2650/*------------------type conversions----------------*/
2651ideal idVec2Ideal(poly vec)
2652{
2653   ideal result=idInit(1,1);
2654   omFree((ADDRESS)result->m);
2655   result->m=NULL; // remove later
2656   pVec2Polys(vec, &(result->m), &(IDELEMS(result)));
2657   return result;
2658}
2659
2660#define NEW_STUFF
2661#ifndef NEW_STUFF
2662// converts mat to module, destroys mat
2663ideal idMatrix2Module(matrix mat)
2664{
2665  int mc=MATCOLS(mat);
2666  int mr=MATROWS(mat);
2667  ideal result = idInit(si_max(mc,1),si_max(mr,1));
2668  int i,j;
2669  poly h;
2670
2671  for(j=0;j<mc /*MATCOLS(mat)*/;j++) /* j is also index in result->m */
2672  {
2673    for (i=1;i<=mr /*MATROWS(mat)*/;i++)
2674    {
2675      h = MATELEM(mat,i,j+1);
2676      if (h!=NULL)
2677      {
2678        MATELEM(mat,i,j+1)=NULL;
2679        pSetCompP(h,i);
2680        result->m[j] = pAdd(result->m[j],h);
2681      }
2682    }
2683  }
2684  // obachman: need to clean this up
2685  idDelete((ideal*) &mat);
2686  return result;
2687}
2688#else
2689
2690#include "sbuckets.h"
2691
2692// converts mat to module, destroys mat
2693ideal idMatrix2Module(matrix mat)
2694{
2695  int mc=MATCOLS(mat);
2696  int mr=MATROWS(mat);
2697  ideal result = idInit(si_max(mc,1),si_max(mr,1));
2698  int i,j, l;
2699  poly h;
2700  poly p;
2701  sBucket_pt bucket = sBucketCreate(currRing);
2702
2703  for(j=0;j<mc /*MATCOLS(mat)*/;j++) /* j is also index in result->m */
2704  {
2705    for (i=1;i<=mr /*MATROWS(mat)*/;i++)
2706    {
2707      h = MATELEM(mat,i,j+1);
2708      if (h!=NULL)
2709      {
2710        l=pLength(h);
2711        MATELEM(mat,i,j+1)=NULL;
2712        p_SetCompP(h,i, currRing);
2713        sBucket_Merge_p(bucket, h, l);
2714      }
2715    }
2716    sBucketClearMerge(bucket, &(result->m[j]), &l);
2717  }
2718  sBucketDestroy(&bucket);
2719
2720  // obachman: need to clean this up
2721  idDelete((ideal*) &mat);
2722  return result;
2723}
2724#endif
2725
2726/*2
2727* converts a module into a matrix, destroyes the input
2728*/
2729matrix idModule2Matrix(ideal mod)
2730{
2731  matrix result = mpNew(mod->rank,IDELEMS(mod));
2732  int i,cp;
2733  poly p,h;
2734
2735  for(i=0;i<IDELEMS(mod);i++)
2736  {
2737    p=mod->m[i];
2738    mod->m[i]=NULL;
2739    while (p!=NULL)
2740    {
2741      h=p;
2742      pIter(p);
2743      pNext(h)=NULL;
2744//      cp = si_max(1,pGetComp(h));     // if used for ideals too
2745      cp = pGetComp(h);
2746      pSetComp(h,0);
2747      pSetmComp(h);
2748#ifdef TEST
2749      if (cp>mod->rank)
2750      {
2751        Print("## inv. rank %d -> %d\n",mod->rank,cp);
2752        int k,l,o=mod->rank;
2753        mod->rank=cp;
2754        matrix d=mpNew(mod->rank,IDELEMS(mod));
2755        for (l=1; l<=o; l++)
2756        {
2757          for (k=1; k<=IDELEMS(mod); k++)
2758          {
2759            MATELEM(d,l,k)=MATELEM(result,l,k);
2760            MATELEM(result,l,k)=NULL;
2761          }
2762        }
2763        idDelete((ideal *)&result);
2764        result=d;
2765      }
2766#endif
2767      MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
2768    }
2769  }
2770  // obachman 10/99: added the following line, otherwise memory leack!
2771  idDelete(&mod);
2772  return result;
2773}
2774
2775matrix idModule2formatedMatrix(ideal mod,int rows, int cols)
2776{
2777  matrix result = mpNew(rows,cols);
2778  int i,cp,r=idRankFreeModule(mod),c=IDELEMS(mod);
2779  poly p,h;
2780
2781  if (r>rows) r = rows;
2782  if (c>cols) c = cols;
2783  for(i=0;i<c;i++)
2784  {
2785    p=mod->m[i];
2786    mod->m[i]=NULL;
2787    while (p!=NULL)
2788    {
2789      h=p;
2790      pIter(p);
2791      pNext(h)=NULL;
2792      cp = pGetComp(h);
2793      if (cp<=r)
2794      {
2795        pSetComp(h,0);
2796        pSetmComp(h);
2797        MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
2798      }
2799      else
2800        pDelete(&h);
2801    }
2802  }
2803  idDelete(&mod);
2804  return result;
2805}
2806
2807/*2
2808* substitute the n-th variable by the monomial e in id
2809* destroy id
2810*/
2811ideal  idSubst(ideal id, int n, poly e)
2812{
2813  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
2814  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
2815
2816  res->rank = id->rank;
2817  for(k--;k>=0;k--)
2818  {
2819    res->m[k]=pSubst(id->m[k],n,e);
2820    id->m[k]=NULL;
2821  }
2822  idDelete(&id);
2823  return res;
2824}
2825
2826BOOLEAN idHomModule(ideal m, ideal Q, intvec **w)
2827{
2828  if (w!=NULL) *w=NULL;
2829  if ((Q!=NULL) && (!idHomIdeal(Q,NULL))) return FALSE;
2830  if (idIs0(m))
2831  {
2832    if (w!=NULL) (*w)=new intvec(m->rank);
2833    return TRUE;
2834  }
2835
2836  int i,j,cmax=2,order=0,ord,* diff,* iscom,diffmin=32000;
2837  poly p=NULL;
2838  int length=IDELEMS(m);
2839  polyset P=m->m;
2840  polyset F=(polyset)omAlloc(length*sizeof(poly));
2841  for (i=length-1;i>=0;i--)
2842  {
2843    p=F[i]=P[i];
2844    cmax=si_max(cmax,(int)pMaxComp(p)+1);
2845  }
2846  diff = (int *)omAlloc0(cmax*sizeof(int));
2847  if (w!=NULL) *w=new intvec(cmax-1);
2848  iscom = (int *)omAlloc0(cmax*sizeof(int));
2849  i=0;
2850  while (i<=length)
2851  {
2852    if (i<length)
2853    {
2854      p=F[i];
2855      while ((p!=NULL) && (!iscom[pGetComp(p)])) pIter(p);
2856    }
2857    if ((p==NULL) && (i<length))
2858    {
2859      i++;
2860    }
2861    else
2862    {
2863      if (p==NULL)
2864      {
2865        i=0;
2866        while ((i<length) && (F[i]==NULL)) i++;
2867        if (i>=length) break;
2868        p = F[i];
2869      }
2870      if (pLexOrder && (currRing->order[0]==ringorder_lp))
2871        order=pTotaldegree(p);
2872      else
2873      //  order = p->order;
2874        order = pFDeg(p,currRing);
2875      order += diff[pGetComp(p)];
2876      p = F[i];
2877//Print("Actual p=F[%d]: ",i);pWrite(p);
2878      F[i] = NULL;
2879      i=0;
2880    }
2881    while (p!=NULL)
2882    {
2883      //if (pLexOrder)
2884      //  ord=pTotaldegree(p);
2885      //else
2886      //  ord = p->order;
2887      ord = pFDeg(p,currRing);
2888      if (!iscom[pGetComp(p)])
2889      {
2890        diff[pGetComp(p)] = order-ord;
2891        iscom[pGetComp(p)] = 1;
2892/*
2893*PrintS("new diff: ");
2894*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
2895*PrintLn();
2896*PrintS("new iscom: ");
2897*for (j=0;j<cmax;j++) Print("%d ",iscom[j]);
2898*PrintLn();
2899*Print("new set %d, order %d, ord %d, diff %d\n",pGetComp(p),order,ord,diff[pGetComp(p)]);
2900*/
2901      }
2902      else
2903      {
2904/*
2905*PrintS("new diff: ");
2906*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
2907*PrintLn();
2908*Print("order %d, ord %d, diff %d\n",order,ord,diff[pGetComp(p)]);
2909*/
2910        if (order != ord+diff[pGetComp(p)])
2911        {
2912          omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
2913          omFreeSize((ADDRESS) diff,cmax*sizeof(int));
2914          omFreeSize((ADDRESS) F,length*sizeof(poly));
2915          delete *w;*w=NULL;
2916          return FALSE;
2917        }
2918      }
2919      pIter(p);
2920    }
2921  }
2922  omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
2923  omFreeSize((ADDRESS) F,length*sizeof(poly));
2924  for (i=1;i<cmax;i++) (**w)[i-1]=diff[i];
2925  for (i=1;i<cmax;i++)
2926  {
2927    if (diff[i]<diffmin) diffmin=diff[i];
2928  }
2929  if (w!=NULL)
2930  {
2931    for (i=1;i<cmax;i++)
2932    {
2933      (**w)[i-1]=diff[i]-diffmin;
2934    }
2935  }
2936  omFreeSize((ADDRESS) diff,cmax*sizeof(int));
2937  return TRUE;
2938}
2939
2940BOOLEAN idTestHomModule(ideal m, ideal Q, intvec *w)
2941{
2942  if ((Q!=NULL) && (!idHomIdeal(Q,NULL)))  { PrintS(" Q not hom\n"); return FALSE;}
2943  if (idIs0(m)) return TRUE;
2944
2945  int cmax=-1;
2946  int i;
2947  poly p=NULL;
2948  int length=IDELEMS(m);
2949  polyset P=m->m;
2950  for (i=length-1;i>=0;i--)
2951  {
2952    p=P[i];
2953    if (p!=NULL) cmax=si_max(cmax,(int)pMaxComp(p)+1);
2954  }
2955  if (w != NULL)
2956  if (w->length()+1 < cmax)
2957  {
2958    // Print("length: %d - %d \n", w->length(),cmax);
2959    return FALSE;
2960  }
2961
2962  if(w!=NULL)
2963    pSetModDeg(w);
2964
2965  for (i=length-1;i>=0;i--)
2966  {
2967    p=P[i];
2968    poly q=p;
2969    if (p!=NULL)
2970    {
2971      int d=pFDeg(p,currRing);
2972      loop
2973      {
2974        pIter(p);
2975        if (p==NULL) break;
2976        if (d!=pFDeg(p,currRing))
2977        {
2978          //pWrite(q); wrp(p); Print(" -> %d - %d\n",d,pFDeg(p,currRing));
2979          if(w!=NULL)
2980            pSetModDeg(NULL);
2981          return FALSE;
2982        }
2983      }
2984    }
2985  }
2986
2987  if(w!=NULL)
2988    pSetModDeg(NULL);
2989
2990  return TRUE;
2991}
2992
2993ideal idJet(ideal i,int d)
2994{
2995  ideal r=idInit((i->nrows)*(i->ncols),i->rank);
2996  r->nrows = i-> nrows;
2997  r->ncols = i-> ncols;
2998  //r->rank = i-> rank;
2999  int k;
3000  for(k=(i->nrows)*(i->ncols)-1;k>=0; k--)
3001  {
3002    r->m[k]=ppJet(i->m[k],d);
3003  }
3004  return r;
3005}
3006
3007ideal idJetW(ideal i,int d, intvec * iv)
3008{
3009  ideal r=idInit(IDELEMS(i),i->rank);
3010  if (ecartWeights!=NULL)
3011  {
3012    WerrorS("cannot compute weighted jets now");
3013  }
3014  else
3015  {
3016    short *w=iv2array(iv);
3017    int k;
3018    for(k=0; k<IDELEMS(i); k++)
3019    {
3020      r->m[k]=ppJetW(i->m[k],d,w);
3021    }
3022    omFreeSize((ADDRESS)w,(pVariables+1)*sizeof(short));
3023  }
3024  return r;
3025}
3026
3027int idMinDegW(ideal M,intvec *w)
3028{
3029  int d=-1;
3030  for(int i=0;i<IDELEMS(M);i++)
3031  {
3032    int d0=pMinDeg(M->m[i],w);
3033    if(-1<d0&&(d0<d||d==-1))
3034      d=d0;
3035  }
3036  return d;
3037}
3038
3039ideal idSeries(int n,ideal M,matrix U,intvec *w)
3040{
3041  for(int i=IDELEMS(M)-1;i>=0;i--)
3042  {
3043    if(U==NULL)
3044      M->m[i]=pSeries(n,M->m[i],NULL,w);
3045    else
3046    {
3047      M->m[i]=pSeries(n,M->m[i],MATELEM(U,i+1,i+1),w);
3048      MATELEM(U,i+1,i+1)=NULL;
3049    }
3050  }
3051  if(U!=NULL)
3052    idDelete((ideal*)&U);
3053  return M;
3054}
3055
3056matrix idDiff(matrix i, int k)
3057{
3058  int e=MATCOLS(i)*MATROWS(i);
3059  matrix r=mpNew(MATROWS(i),MATCOLS(i));
3060  int j;
3061  for(j=0; j<e; j++)
3062  {
3063    r->m[j]=pDiff(i->m[j],k);
3064  }
3065  return r;
3066}
3067
3068matrix idDiffOp(ideal I, ideal J,BOOLEAN multiply)
3069{
3070  matrix r=mpNew(IDELEMS(I),IDELEMS(J));
3071  int i,j;
3072  for(i=0; i<IDELEMS(I); i++)
3073  {
3074    for(j=0; j<IDELEMS(J); j++)
3075    {
3076      MATELEM(r,i+1,j+1)=pDiffOp(I->m[i],J->m[j],multiply);
3077    }
3078  }
3079  return r;
3080}
3081
3082/*3
3083*handles for some ideal operations the ring/syzcomp managment
3084*returns all syzygies (componentwise-)shifted by -syzcomp
3085*or -syzcomp-1 (in case of ideals as input)
3086static ideal idHandleIdealOp(ideal arg,int syzcomp,int isIdeal=FALSE)
3087{
3088  ring orig_ring=currRing;
3089  ring syz_ring=rCurrRingAssure_SyzComp();
3090  rSetSyzComp(length);
3091
3092  ideal s_temp;
3093  if (orig_ring!=syz_ring)
3094    s_temp=idrMoveR_NoSort(arg,orig_ring);
3095  else
3096    s_temp=arg;
3097
3098  ideal s_temp1 = kStd(s_temp,currQuotient,testHomog,&w,NULL,length);
3099  if (w!=NULL) delete w;
3100
3101  if (syz_ring!=orig_ring)
3102  {
3103    idDelete(&s_temp);
3104    rChangeCurrRing(orig_ring);
3105  }
3106
3107  idDelete(&temp);
3108  ideal temp1=idRingCopy(s_temp1,syz_ring);
3109
3110  if (syz_ring!=orig_ring)
3111  {
3112    rChangeCurrRing(syz_ring);
3113    idDelete(&s_temp1);
3114    rChangeCurrRing(orig_ring);
3115    rKill(syz_ring);
3116  }
3117
3118  for (i=0;i<IDELEMS(temp1);i++)
3119  {
3120    if ((temp1->m[i]!=NULL)
3121    && (pGetComp(temp1->m[i])<=length))
3122    {
3123      pDelete(&(temp1->m[i]));
3124    }
3125    else
3126    {
3127      pShift(&(temp1->m[i]),-length);
3128    }
3129  }
3130  temp1->rank = rk;
3131  idSkipZeroes(temp1);
3132
3133  return temp1;
3134}
3135*/
3136/*2
3137* represents (h1+h2)/h2=h1/(h1 intersect h2)
3138*/
3139//ideal idModulo (ideal h2,ideal h1)
3140ideal idModulo (ideal h2,ideal h1, tHomog hom, intvec ** w)
3141{
3142  intvec *wtmp=NULL;
3143
3144  int i,j,k,rk,flength=0,slength,length;
3145  poly p,q;
3146
3147  if (idIs0(h2))
3148    return idFreeModule(si_max(1,h2->ncols));
3149  if (!idIs0(h1))
3150    flength = idRankFreeModule(h1);
3151  slength = idRankFreeModule(h2);
3152  length  = si_max(flength,slength);
3153  if (length==0)
3154  {
3155    length = 1;
3156  }
3157  ideal temp = idInit(IDELEMS(h2),length+IDELEMS(h2));
3158  if ((w!=NULL)&&((*w)!=NULL))
3159  {
3160    //Print("input weights:");(*w)->show(1);PrintLn();
3161    int d;
3162    int k;
3163    wtmp=new intvec(length+IDELEMS(h2));
3164    for (i=0;i<length;i++)
3165      ((*wtmp)[i])=(**w)[i];
3166    for (i=0;i<IDELEMS(h2);i++)
3167    {
3168      poly p=h2->m[i];
3169      if (p!=NULL)
3170      {
3171        d = pDeg(p);
3172        k= pGetComp(p);
3173        if (slength>0) k--;
3174        d +=((**w)[k]);
3175        ((*wtmp)[i+length]) = d;
3176      }
3177    }
3178    //Print("weights:");wtmp->show(1);PrintLn();
3179  }
3180  for (i=0;i<IDELEMS(h2);i++)
3181  {
3182    temp->m[i] = pCopy(h2->m[i]);
3183    q = pOne();
3184    pSetComp(q,i+1+length);
3185    pSetmComp(q);
3186    if(temp->m[i]!=NULL)
3187    {
3188      if (slength==0) pShift(&(temp->m[i]),1);
3189      p = temp->m[i];
3190      while (pNext(p)!=NULL) pIter(p);
3191      pNext(p) = q;
3192    }
3193    else
3194      temp->m[i]=q;
3195  }
3196  rk = k = IDELEMS(h2);
3197  if (!idIs0(h1))
3198  {
3199    pEnlargeSet(&(temp->m),IDELEMS(temp),IDELEMS(h1));
3200    IDELEMS(temp) += IDELEMS(h1);
3201    for (i=0;i<IDELEMS(h1);i++)
3202    {
3203      if (h1->m[i]!=NULL)
3204      {
3205        temp->m[k] = pCopy(h1->m[i]);
3206        if (flength==0) pShift(&(temp->m[k]),1);
3207        k++;
3208      }
3209    }
3210  }
3211
3212  ring orig_ring=currRing;
3213  ring syz_ring=rCurrRingAssure_SyzComp();
3214  rSetSyzComp(length);
3215  ideal s_temp;
3216
3217  if (syz_ring != orig_ring)
3218  {
3219    s_temp = idrMoveR_NoSort(temp, orig_ring);
3220  }
3221  else
3222  {
3223    s_temp = temp;
3224  }
3225
3226  idTest(s_temp);
3227  ideal s_temp1 = kStd(s_temp,currQuotient,hom,&wtmp,NULL,length);
3228
3229  //if (wtmp!=NULL)  Print("output weights:");wtmp->show(1);PrintLn();
3230  if ((w!=NULL) && (*w !=NULL) && (wtmp!=NULL))
3231  {
3232    delete *w;
3233    *w=new intvec(IDELEMS(h2));
3234    for (i=0;i<IDELEMS(h2);i++)
3235      ((**w)[i])=(*wtmp)[i+length];
3236  }
3237  if (wtmp!=NULL) delete wtmp;
3238
3239  for (i=0;i<IDELEMS(s_temp1);i++)
3240  {
3241    if ((s_temp1->m[i]!=NULL)
3242    && (pGetComp(s_temp1->m[i])<=length))
3243    {
3244      pDelete(&(s_temp1->m[i]));
3245    }
3246    else
3247    {
3248      pShift(&(s_temp1->m[i]),-length);
3249    }
3250  }
3251  s_temp1->rank = rk;
3252  idSkipZeroes(s_temp1);
3253
3254  if (syz_ring!=orig_ring)
3255  {
3256    rChangeCurrRing(orig_ring);
3257    s_temp1 = idrMoveR_NoSort(s_temp1, syz_ring);
3258    rKill(syz_ring);
3259    // Hmm ... here seems to be a memory leak
3260    // However, simply deleting it causes memory trouble
3261    // idDelete(&s_temp);
3262  }
3263  else
3264  {
3265    idDelete(&temp);
3266  }
3267  idTest(s_temp1);
3268  return s_temp1;
3269}
3270
3271int idElem(const ideal F)
3272{
3273  int i=0,j=IDELEMS(F)-1;
3274
3275  while(j>=0)
3276  {
3277    if ((F->m)[j]!=NULL) i++;
3278    j--;
3279  }
3280  return i;
3281}
3282
3283/*
3284*computes module-weights for liftings of homogeneous modules
3285*/
3286intvec * idMWLift(ideal mod,intvec * weights)
3287{
3288  if (idIs0(mod)) return new intvec(2);
3289  int i=IDELEMS(mod);
3290  while ((i>0) && (mod->m[i-1]==NULL)) i--;
3291  intvec *result = new intvec(i+1);
3292  while (i>0)
3293  {
3294    (*result)[i]=pFDeg(mod->m[i],currRing)+(*weights)[pGetComp(mod->m[i])];
3295  }
3296  return result;
3297}
3298
3299/*2
3300*sorts the kbase for idCoef* in a special way (lexicographically
3301*with x_max,...,x_1)
3302*/
3303ideal idCreateSpecialKbase(ideal kBase,intvec ** convert)
3304{
3305  int i;
3306  ideal result;
3307
3308  if (idIs0(kBase)) return NULL;
3309  result = idInit(IDELEMS(kBase),kBase->rank);
3310  *convert = idSort(kBase,FALSE);
3311  for (i=0;i<(*convert)->length();i++)
3312  {
3313    result->m[i] = pCopy(kBase->m[(**convert)[i]-1]);
3314  }
3315  return result;
3316}
3317
3318/*2
3319*returns the index of a given monom in the list of the special kbase
3320*/
3321int idIndexOfKBase(poly monom, ideal kbase)
3322{
3323  int j=IDELEMS(kbase);
3324
3325  while ((j>0) && (kbase->m[j-1]==NULL)) j--;
3326  if (j==0) return -1;
3327  int i=pVariables;
3328  while (i>0)
3329  {
3330    loop
3331    {
3332      if (pGetExp(monom,i)>pGetExp(kbase->m[j-1],i)) return -1;
3333      if (pGetExp(monom,i)==pGetExp(kbase->m[j-1],i)) break;
3334      j--;
3335      if (j==0) return -1;
3336    }
3337    if (i==1)
3338    {
3339      while(j>0)
3340      {
3341        if (pGetComp(monom)==pGetComp(kbase->m[j-1])) return j-1;
3342        if (pGetComp(monom)>pGetComp(kbase->m[j-1])) return -1;
3343        j--;
3344      }
3345    }
3346    i--;
3347  }
3348  return -1;
3349}
3350
3351/*2
3352*decomposes the monom in a part of coefficients described by the
3353*complement of how and a monom in variables occuring in how, the
3354*index of which in kbase is returned as integer pos (-1 if it don't
3355*exists)
3356*/
3357poly idDecompose(poly monom, poly how, ideal kbase, int * pos)
3358{
3359  int i;
3360  poly coeff=pOne(), base=pOne();
3361
3362  for (i=1;i<=pVariables;i++)
3363  {
3364    if (pGetExp(how,i)>0)
3365    {
3366      pSetExp(base,i,pGetExp(monom,i));
3367    }
3368    else
3369    {
3370      pSetExp(coeff,i,pGetExp(monom,i));
3371    }
3372  }
3373  pSetComp(base,pGetComp(monom));
3374  pSetm(base);
3375  pSetCoeff(coeff,nCopy(pGetCoeff(monom)));
3376  pSetm(coeff);
3377  *pos = idIndexOfKBase(base,kbase);
3378  if (*pos<0)
3379    pDelete(&coeff);
3380  pDelete(&base);
3381  return coeff;
3382}
3383
3384/*2
3385*returns a matrix A of coefficients with kbase*A=arg
3386*if all monomials in variables of how occur in kbase
3387*the other are deleted
3388*/
3389matrix idCoeffOfKBase(ideal arg, ideal kbase, poly how)
3390{
3391  matrix result;
3392  ideal tempKbase;
3393  poly p,q;
3394  intvec * convert;
3395  int i=IDELEMS(kbase),j=IDELEMS(arg),k,pos;
3396#if 0
3397  while ((i>0) && (kbase->m[i-1]==NULL)) i--;
3398  if (idIs0(arg))
3399    return mpNew(i,1);
3400  while ((j>0) && (arg->m[j-1]==NULL)) j--;
3401  result = mpNew(i,j);
3402#else
3403  result = mpNew(i, j);
3404  while ((j>0) && (arg->m[j-1]==NULL)) j--;
3405#endif
3406
3407  tempKbase = idCreateSpecialKbase(kbase,&convert);
3408  for (k=0;k<j;k++)
3409  {
3410    p = arg->m[k];
3411    while (p!=NULL)
3412    {
3413      q = idDecompose(p,how,tempKbase,&pos);
3414      if (pos>=0)
3415      {
3416        MATELEM(result,(*convert)[pos],k+1) =
3417            pAdd(MATELEM(result,(*convert)[pos],k+1),q);
3418      }
3419      else
3420        pDelete(&q);
3421      pIter(p);
3422    }
3423  }
3424  idDelete(&tempKbase);
3425  return result;
3426}
3427
3428/*3
3429* searches for units in the components of the module arg and
3430* returns the first one
3431*/
3432static int idReadOutUnits(ideal arg,int* comp)
3433{
3434  if (idIs0(arg)) return -1;
3435  int i=0,j, generator=-1;
3436  int rk_arg=arg->rank; //idRankFreeModule(arg);
3437  int * componentIsUsed =(int *)omAlloc((rk_arg+1)*sizeof(int));
3438  poly p,q;
3439
3440  while ((generator<0) && (i<IDELEMS(arg)))
3441  {
3442    memset(componentIsUsed,0,(rk_arg+1)*sizeof(int));
3443    p = arg->m[i];
3444    while (p!=NULL)
3445    {
3446      j = pGetComp(p);
3447      if (componentIsUsed[j]==0)
3448      {
3449        if (pLmIsConstantComp(p))
3450        {
3451          generator = i;
3452          componentIsUsed[j] = 1;
3453        }
3454        else
3455        {
3456          componentIsUsed[j] = -1;
3457        }
3458      }
3459      else if (componentIsUsed[j]>0)
3460      {
3461        (componentIsUsed[j])++;
3462      }
3463      pIter(p);
3464    }
3465    i++;
3466  }
3467  i = 0;
3468  *comp = -1;
3469  for (j=0;j<=rk_arg;j++)
3470  {
3471    if (componentIsUsed[j]>0)
3472    {
3473      if ((*comp==-1) || (componentIsUsed[j]<i))
3474      {
3475        *comp = j;
3476        i= componentIsUsed[j];
3477      }
3478    }
3479  }
3480  omFree(componentIsUsed);
3481  return generator;
3482}
3483
3484#if 0
3485static void idDeleteComp(ideal arg,int red_comp)
3486{
3487  int i,j;
3488  poly p;
3489
3490  for (i=IDELEMS(arg)-1;i>=0;i--)
3491  {
3492    p = arg->m[i];
3493    while (p!=NULL)
3494    {
3495      j = pGetComp(p);
3496      if (j>red_comp)
3497      {
3498        pSetComp(p,j-1);
3499        pSetm(p);
3500      }
3501      pIter(p);
3502    }
3503  }
3504  (arg->rank)--;
3505}
3506#endif
3507
3508static void idDeleteComps(ideal arg,int* red_comp,int del)
3509// red_comp is an array [0..args->rank]
3510{
3511  int i,j;
3512  poly p;
3513
3514  for (i=IDELEMS(arg)-1;i>=0;i--)
3515  {
3516    p = arg->m[i];
3517    while (p!=NULL)
3518    {
3519      j = pGetComp(p);
3520      if (red_comp[j]!=j)
3521      {
3522        pSetComp(p,red_comp[j]);
3523        pSetmComp(p);
3524      }
3525      pIter(p);
3526    }
3527  }
3528  (arg->rank) -= del;
3529}
3530
3531/*2
3532* returns the presentation of an isomorphic, minimally
3533* embedded  module (arg represents the quotient!)
3534*/
3535ideal idMinEmbedding(ideal arg,BOOLEAN inPlace, intvec **w)
3536{
3537  if (idIs0(arg)) return idInit(1,arg->rank);
3538  int i,next_gen,next_comp;
3539  ideal res=arg;
3540
3541  if (!inPlace) res = idCopy(arg);
3542  res->rank=si_max(res->rank,idRankFreeModule(res));
3543  int *red_comp=(int*)omAlloc((res->rank+1)*sizeof(int));
3544  for (i=res->rank;i>=0;i--) red_comp[i]=i;
3545
3546  int del=0;
3547  loop
3548  {
3549    next_gen = idReadOutUnits(res,&next_comp);
3550    if (next_gen<0) break;
3551    del++;
3552    syGaussForOne(res,next_gen,next_comp,0,IDELEMS(res));
3553    for(i=next_comp+1;i<=arg->rank;i++) red_comp[i]--;
3554    if ((w !=NULL)&&(*w!=NULL))
3555    {
3556      for(i=next_comp;i<(*w)->length();i++) (**w)[i-1]=(**w)[i];
3557    }
3558  }
3559
3560  idDeleteComps(res,red_comp,del);
3561  idSkipZeroes(res);
3562  omFree(red_comp);
3563
3564  if ((w !=NULL)&&(*w!=NULL) &&(del>0))
3565  {
3566    intvec *wtmp=new intvec((*w)->length()-del);
3567    for(i=0;i<res->rank;i++) (*wtmp)[i]=(**w)[i];
3568    delete *w;
3569    *w=wtmp;
3570  }
3571  return res;
3572}
3573
3574/*2
3575* transpose a module
3576*/
3577ideal idTransp(ideal a)
3578{
3579  int r = a->rank, c = IDELEMS(a);
3580  ideal b =  idInit(r,c);
3581
3582  for (int i=c; i>0; i--)
3583  {
3584    poly p=a->m[i-1];
3585    while(p!=NULL)
3586    {
3587      poly h=pHead(p);
3588      int co=pGetComp(h)-1;
3589      pSetComp(h,i);
3590      pSetmComp(h);
3591      b->m[co]=pAdd(b->m[co],h);
3592      pIter(p);
3593    }
3594  }
3595  return b;
3596}
3597
3598intvec * idQHomWeight(ideal id)
3599{
3600  poly head, tail;
3601  int k;
3602  int in=IDELEMS(id)-1, ready=0, all=0,
3603      coldim=pVariables, rowmax=2*coldim;
3604  if (in<0) return NULL;
3605  intvec *imat=new intvec(rowmax+1,coldim,0);
3606
3607  do
3608  {
3609    head = id->m[in--];
3610    if (head!=NULL)
3611    {
3612      tail = pNext(head);
3613      while (tail!=NULL)
3614      {
3615        all++;
3616        for (k=1;k<=coldim;k++)
3617          IMATELEM(*imat,all,k) = pGetExpDiff(head,tail,k);
3618        if (all==rowmax)
3619        {
3620          ivTriangIntern(imat, ready, all);
3621          if (ready==coldim)
3622          {
3623            delete imat;
3624            return NULL;
3625          }
3626        }
3627        pIter(tail);
3628      }
3629    }
3630  } while (in>=0);
3631  if (all>ready)
3632  {
3633    ivTriangIntern(imat, ready, all);
3634    if (ready==coldim)
3635    {
3636      delete imat;
3637      return NULL;
3638    }
3639  }
3640  intvec *result = ivSolveKern(imat, ready);
3641  delete imat;
3642  return result;
3643}
3644
3645BOOLEAN idIsZeroDim(ideal I)
3646{
3647  BOOLEAN *UsedAxis=(BOOLEAN *)omAlloc0(pVariables*sizeof(BOOLEAN));
3648  int i,n;
3649  poly po;
3650  BOOLEAN res=TRUE;
3651  for(i=IDELEMS(I)-1;i>=0;i--)
3652  {
3653    po=I->m[i];
3654    if ((po!=NULL) &&((n=pIsPurePower(po))!=0)) UsedAxis[n-1]=TRUE;
3655  }
3656  for(i=pVariables-1;i>=0;i--)
3657  {
3658    if(UsedAxis[i]==FALSE) {res=FALSE; break;} // not zero-dim.
3659  }
3660  omFreeSize(UsedAxis,pVariables*sizeof(BOOLEAN));
3661  return res;
3662}
3663
3664void idNormalize(ideal I)
3665{
3666  if (rField_has_simple_inverse()) return; /* Z/p, GF(p,n), R, long R/C */
3667  int i;
3668  poly p;
3669  for(i=IDELEMS(I)-1;i>=0;i--)
3670  {
3671    p=I->m[i] ;
3672    while(p!=NULL)
3673    {
3674      nNormalize(pGetCoeff(p));
3675      pIter(p);
3676    }
3677  }
3678}
3679
3680#include "clapsing.h"
3681
3682poly id_GCD(poly f, poly g, const ring r)
3683{
3684  ring save_r=currRing;
3685  rChangeCurrRing(r);
3686  ideal I=idInit(2,1); I->m[0]=f; I->m[1]=g;
3687  intvec *w = NULL;
3688  ideal S=idSyzygies(I,testHomog,&w);
3689  if (w!=NULL) delete w;
3690  poly gg=pTakeOutComp(&(S->m[0]),2);
3691  idDelete(&S);
3692  poly gcd_p=singclap_pdivide(f,gg);
3693  pDelete(&gg);
3694  rChangeCurrRing(save_r);
3695  return gcd_p;
3696}
3697
3698/*2
3699* xx,q: arrays of length 0..rl-1
3700* xx[i]: SB mod q[i]
3701* assume: char=0
3702* assume: q[i]!=0
3703* destroys xx
3704*/
3705ideal idChineseRemainder(ideal *xx, number *q, int rl)
3706{
3707  ideal result=idInit(IDELEMS(xx[0]),1);
3708  int i,j;
3709  poly r,h,hh,res_p;
3710  number *x=(number *)omAlloc(rl*sizeof(number));
3711  for(i=IDELEMS(result)-1;i>=0;i--)
3712  {
3713    res_p=NULL;
3714    loop
3715    {
3716      r=NULL;
3717      for(j=rl-1;j>=0;j--)
3718      {
3719        h=xx[j]->m[i];
3720        if ((h!=NULL)
3721        &&((r==NULL)||(pLmCmp(r,h)==-1)))
3722          r=h;
3723      }
3724      if (r==NULL) break;
3725      h=pHead(r);
3726      for(j=rl-1;j>=0;j--)
3727      {
3728        hh=xx[j]->m[i];
3729        if ((hh!=NULL) && (pLmCmp(r,hh)==0))
3730        {
3731          x[j]=pGetCoeff(hh);
3732          hh=pLmFreeAndNext(hh);
3733          xx[j]->m[i]=hh;
3734        }
3735        else
3736          x[j]=nlInit(0);
3737      }
3738      number n=nlChineseRemainder(x,q,rl);
3739      for(j=rl-1;j>=0;j--)
3740      {
3741        nlDelete(&(x[j]),currRing);
3742      }
3743      pSetCoeff(h,n);
3744      //Print("new mon:");pWrite(h);
3745      res_p=pAdd(res_p,h);
3746    }
3747    result->m[i]=res_p;
3748  }
3749  omFree(x);
3750  for(i=rl-1;i>=0;i--) idDelete(&(xx[i]));
3751  omFree(xx);
3752  return result;
3753}
3754ideal idChineseRemainder(ideal *xx, intvec *iv)
3755{
3756  int rl=iv->length();
3757  number *q=(number *)omAlloc(rl*sizeof(number));
3758  int i;
3759  for(i=0; i<rl; i++)
3760  {
3761    q[i]=nInit((*iv)[i]);
3762  }
3763  return idChineseRemainder(xx,q,rl);
3764}
Note: See TracBrowser for help on using the repository browser.