source: git/kernel/ideals.cc @ fda6e6e

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