source: git/kernel/ideals.cc @ 3060a7

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