source: git/kernel/ideals.cc @ 230d82

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