source: git/libpolys/polys/ideals.cc @ 68e548

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