source: git/kernel/ideals.cc @ 920804

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