source: git/kernel/ideals.cc @ b83062

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