source: git/Singular/ideals.cc @ c4bbf1f

spielwiese
Last change on this file since c4bbf1f was 4b5c87, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: several changes concerning stdfac git-svn-id: file:///usr/local/Singular/svn/trunk@1416 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 61.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: ideals.cc,v 1.22 1998-04-22 07:48:54 Singular Exp $ */
5/*
6* ABSTRACT - all basic methods to manipulate ideals
7*/
8
9/* includes */
10#include "mod2.h"
11#include "tok.h"
12#include "mmemory.h"
13#include "febase.h"
14#include "numbers.h"
15#include "polys.h"
16#include "ipid.h"
17#include "ring.h"
18#include "kstd1.h"
19#include "matpol.h"
20#include "weight.h"
21#include "intvec.h"
22#include "syz.h"
23#include "ideals.h"
24#include "lists.h"
25
26
27static poly * idpower;
28/*collects the monomials in makemonoms, must be allocated befor*/
29static int idpowerpoint;
30/*index of the actual monomial in idpower*/
31static poly * givenideal;
32/*the ideal from which a power is computed*/
33
34/*0 implementation*/
35
36/*2
37* initialise an ideal
38*/
39#ifdef MDEBUG
40ideal idDBInit(int idsize, int rank, char *f, int l)
41#else
42ideal idInit(int idsize, int rank)
43#endif
44{
45  /*- initialise an ideal -*/
46#ifdef MDEBUG
47  ideal hh = (ideal )mmDBAllocBlock(sizeof(*hh),f,l);
48#else
49  ideal hh = (ideal )Alloc(sizeof(*hh));
50#endif
51  hh->nrows = 1;
52  hh->rank = rank;
53  IDELEMS(hh) = idsize;
54  if (idsize>0)
55  {
56#ifdef MDEBUG
57    hh->m = (poly *)mmDBAllocBlock0(idsize*sizeof(poly),f,l);
58#else
59    hh->m = (poly *)Alloc0(idsize*sizeof(poly));
60#endif
61  }
62  else
63    hh->m=NULL;
64  return hh;
65}
66
67/*2
68* initialise the maximal ideal (at 0)
69*/
70ideal idMaxIdeal (void)
71{
72  int l;
73  ideal hh=NULL;
74
75  hh=idInit(pVariables,1);
76  for (l=0; l<pVariables; l++)
77  {
78    hh->m[l] = pOne();
79    pSetExp(hh->m[l],l+1,1);
80    pSetm(hh->m[l]);
81  }
82  return hh;
83}
84
85/*2
86* deletes an ideal/matrix
87*/
88#ifdef MDEBUG
89void idDBDelete (ideal* h, char *f, int l)
90#else
91void idDelete (ideal * h)
92#endif
93{
94  int j,elems;
95  if (*h == NULL)
96    return;
97  elems=j=(*h)->nrows*(*h)->ncols;
98  if (j>0)
99  {
100    do
101    {
102      #ifdef MDEBUG
103      pDBDelete(&((*h)->m[--j]),f,l);
104      #else
105      pDelete(&((*h)->m[--j]));
106      #endif
107    }
108    while (j>0);
109    #ifdef MDEBUG
110    mmDBFreeBlock((ADDRESS)((*h)->m),sizeof(poly)*elems,f,l);
111    #else
112    Free((ADDRESS)((*h)->m),sizeof(poly)*elems);
113    #endif
114  }
115  #ifdef MDEBUG
116  mmDBFreeBlock((ADDRESS)(*h),sizeof(**h),f,l);
117  #else
118  Free((ADDRESS)*h,sizeof(**h));
119  #endif
120  *h=NULL;
121}
122
123/*2
124*gives an ideal the minimal possible size
125*/
126void idSkipZeroes (ideal ide)
127{
128  int k;
129  int j = -1;
130  BOOLEAN change=FALSE;
131  for (k=0; k<IDELEMS(ide); k++)
132  {
133    if (ide->m[k] != NULL)
134    {
135      j++;
136      if (change)
137      {
138        ide->m[j] = ide->m[k];
139      }
140    }
141    else
142    {
143      change=TRUE;
144    }
145  }
146  if (change)
147  {
148    if (j == -1)
149      j = 0;
150    else
151    {
152      for (k=j+1; k<IDELEMS(ide); k++)
153        ide->m[k] = NULL;
154    }
155    pEnlargeSet(&(ide->m),IDELEMS(ide),j+1-IDELEMS(ide));
156    IDELEMS(ide) = j+1;
157  }
158}
159
160/*2
161* ideal id = (id[i])
162* result is leadcoeff(id[i]) = 1
163*/
164void idNorm(ideal id)
165{
166  for (int i=0; i<IDELEMS(id); i++)
167  {
168    if (id->m[i] != NULL)
169    {
170      pNorm(id->m[i]);
171    }
172  }
173}
174
175/*2
176* ideal id = (id[i]), c any number
177* if id[i] = c*id[j] then id[j] is deleted for j > i
178*/
179void idDelMultiples(ideal id)
180{
181  int i, j, t;
182  int k = IDELEMS(id), l = k;
183  for (i=k-2; i>=0; i--)
184  {
185    if (id->m[i]!=NULL)
186    {
187      for (j=l-1; j>i; j--)
188      {
189        if ((id->m[j]!=NULL)
190        && (pComparePolys(id->m[i], id->m[j])))
191        {
192          pDelete(&id->m[j]);
193          l--;
194          for(t=j; t<l; t++)
195          {
196            id->m[t] = id->m[t+1];
197          }
198        }
199      }
200    }
201  }
202  if (l != k)
203  {
204    pEnlargeSet(&id->m, k, l-k);
205    IDELEMS(id) = l;
206  }
207}
208
209/*2
210* ideal id = (id[i])
211* if id[i] = id[j] then id[j] is deleted for j > i
212*/
213void idDelEquals(ideal id)
214{
215  int i, j, t;
216  int k = IDELEMS(id), l = k;
217  for (i=k-2; i>=0; i--)
218  {
219    for (j=l-1; j>i; j--)
220    {
221      if (pEqualPolys(id->m[i], id->m[j]))
222      {
223        pDelete(&id->m[j]);
224        l--;
225        for(t=j; t<l; t++)
226        {
227          id->m[t] = id->m[t+1];
228        }
229      }
230    }
231  }
232  if (l != k)
233  {
234    pEnlargeSet(&id->m, k, l-k);
235    IDELEMS(id) = l;
236  }
237}
238
239/*2
240* copy an ideal
241*/
242#ifdef MDEBUG
243ideal idDBCopy(ideal h1,char *f,int l)
244#else
245ideal idCopy (ideal h1)
246#endif
247{
248  int i;
249  ideal h2;
250
251//#ifdef TEST
252  if (h1 == NULL)
253  {
254#ifdef MDEBUG
255    h2=idDBInit(1,1,f,l);
256#else
257    h2=idInit(1,1);
258#endif
259  }
260  else
261//#endif
262  {
263#ifdef MDEBUG
264    h2=idDBInit(IDELEMS(h1),h1->rank,f,l);
265#else
266    h2=idInit(IDELEMS(h1),h1->rank);
267#endif
268    for (i=IDELEMS(h1)-1; i>=0; i--)
269#ifdef MDEBUG
270      h2->m[i] = pDBCopy(h1->m[i],f,l);
271#else
272      h2->m[i] = pCopy(h1->m[i]);
273#endif
274  }
275  return h2;
276}
277
278#ifdef PDEBUG
279void idDBTest(ideal h1,char *f,int l)
280{
281  int i;
282
283  if (h1 != NULL)
284  {
285    #ifdef MDEBUG
286    mmDBTestBlock(h1,sizeof(*h1),f,l);
287    #endif
288    for (i=IDELEMS(h1)-1; i>=0; i--)
289      pDBTest(h1->m[i],f,l);
290  }
291}
292#endif
293
294/*3
295* for idSort: compare a and b revlex inclusive module comp.
296*/
297static int pComp_RevLex(poly a, poly b)
298{
299 int l=pVariables;
300 while ((l>=0) && (pGetExp(a,l)==pGetExp(b,l))) l--;
301 if (l<0)
302   return 0;
303 else if (pGetExp(a,l)>pGetExp(b,l))
304   return 1;
305 else
306   return -1;
307}
308
309/*2
310*sorts the ideal w.r.t. the actual ringordering
311*uses lex-ordering when nolex = FALSE
312*/
313intvec *idSort(ideal id,BOOLEAN nolex)
314{
315  poly p,q;
316  intvec * result = new intvec(IDELEMS(id));
317  int i, j, actpos=0, newpos, l;
318  int diff, olddiff, lastcomp, newcomp;
319  BOOLEAN notFound;
320
321  pCompProc oldComp=pComp0;
322
323  if (!nolex) pComp0=pComp_RevLex;
324
325  for (i=0;i<IDELEMS(id);i++)
326  {
327    if (id->m[i]!=NULL)
328    {
329      notFound = TRUE;
330      newpos = actpos / 2;
331      diff = (actpos+1) / 2;
332      diff = (diff+1) / 2;
333      lastcomp = pComp0(id->m[i],id->m[(*result)[newpos]]);
334      if (lastcomp<0)
335      {
336        newpos -= diff;
337      }
338      else if (lastcomp>0)
339      {
340        newpos += diff;
341      }
342      else
343      {
344        notFound = FALSE;
345      }
346      while ((newpos>=0) && (newpos<actpos) && (notFound))
347      {
348        newcomp = pComp0(id->m[i],id->m[(*result)[newpos]]);
349        olddiff = diff;
350        if (diff>1)
351        {
352          diff = (diff+1) / 2;
353          if ((newcomp==1)
354          && (actpos-newpos>1)
355          && (diff>1)
356          && (newpos+diff>=actpos))
357          {
358            diff = actpos-newpos-1;
359          }
360          else if ((newcomp==-1)
361          && (diff>1)
362          && (newpos<diff))
363          {
364            diff = newpos;
365          }
366        }
367        if (newcomp<0)
368        {
369          if ((olddiff==1) && (lastcomp>0))
370            notFound = FALSE;
371          else
372            newpos -= diff;
373        }
374        else if (newcomp>0)
375        {
376          if ((olddiff==1) && (lastcomp<0))
377          {
378            notFound = FALSE;
379            newpos++;
380          }
381          else
382          {
383            newpos += diff;
384          }
385        }
386        else
387        {
388          notFound = FALSE;
389        }
390        lastcomp = newcomp;
391        if (diff==0) notFound=FALSE; /*hs*/
392      }
393      if (newpos<0) newpos = 0;
394      if (newpos>=actpos)
395      {
396        (*result)[actpos] = i;
397      }
398      else
399      {
400        for (j=actpos;j>newpos;j--)
401        {
402          (*result)[j] = (*result)[j-1];
403        }
404        (*result)[newpos] = i;
405      }
406      actpos++;
407    }
408  }
409  for (j=0;j<actpos;j++) (*result)[j]++;
410  pComp0=oldComp;
411  return result;
412}
413
414/*2
415* concat the lists h1 and h2 without zeros
416*/
417ideal idSimpleAdd (ideal h1,ideal h2)
418{
419  int i,j,r,l;
420  ideal result;
421
422  if (h1==NULL) return idCopy(h2);
423  if (h2==NULL) return idCopy(h1);
424  j = IDELEMS(h1)-1;
425  while ((j >= 0) && (h1->m[j] == NULL)) j--;
426  i = IDELEMS(h2)-1;
427  while ((i >= 0) && (h2->m[i] == NULL)) i--;
428  r = max(h1->rank,h2->rank);
429  if (i+j==(-2))
430    return idInit(1,r);
431  else
432    result=idInit(i+j+2,r);
433  for (l=j; l>=0; l--)
434  {
435    result->m[l] = pCopy(h1->m[l]);
436  }
437  r = i+j+1;
438  for (l=i; l>=0; l--, r--)
439  {
440    result->m[r] = pCopy(h2->m[l]);
441  }
442  return result;
443}
444
445/*2
446* h1 + h2
447*/
448ideal idAdd (ideal h1,ideal h2)
449{
450  ideal result = idSimpleAdd(h1,h2);
451  ideal tmp = idCompactify(result);
452
453  idDelete(&result);
454  return tmp;
455}
456
457/*2
458* h1 * h2
459*/
460ideal  idMult (ideal h1,ideal  h2)
461{
462  int i,j,k;
463  ideal  hh;
464
465  j = IDELEMS(h1);
466  while ((j > 0) && (h1->m[j-1] == NULL)) j--;
467  i = IDELEMS(h2);
468  while ((i > 0) && (h2->m[i-1] == NULL)) i--;
469  j = j * i;
470  if (j == 0)
471    hh = idInit(1,1);
472  else
473    hh=idInit(j,1);
474  if (h1->rank<h2->rank)
475    hh->rank = h2->rank;
476  else
477    hh->rank = h1->rank;
478  if (j==0) return hh;
479  k = 0;
480  for (i=0; i<IDELEMS(h1); i++)
481  {
482    if (h1->m[i] != NULL)
483    {
484      for (j=0; j<IDELEMS(h2); j++)
485      {
486        if (h2->m[j] != NULL)
487        {
488          hh->m[k] = pMult(pCopy(h1->m[i]),pCopy(h2->m[j]));
489          k++;
490        }
491      }
492    }
493  }
494  {
495    ideal tmp = idCompactify(hh);
496    idDelete(&hh);
497    return tmp;
498    return hh;
499  }
500}
501
502/*2
503*returns true if h is the zero ideal
504*/
505BOOLEAN idIs0 (ideal h)
506{
507  int i;
508
509  if (h == NULL) return TRUE;
510  i = IDELEMS(h);
511  while ((i > 0) && (h->m[i-1] == NULL))
512  {
513    i--;
514  }
515  if (i == 0)
516    return TRUE;
517  else
518    return FALSE;
519}
520
521/*2
522* return the maximal component number found in any polynomial in s
523*/
524int idRankFreeModule (ideal s)
525{
526  if (s!=NULL)
527  {
528    int  j=0;
529    int  l=IDELEMS(s);
530    poly *p=s->m;
531    int  k;
532
533    for (; l != 0; l--)
534    {
535      if (*p!=NULL)
536      {
537        k = pMaxComp(*p);
538        if (k>j) j = k;
539      }
540      p++;
541    }
542    return j;
543  }
544  return -1;
545}
546
547/*2
548*returns true if id is homogenous with respect to the aktual weights
549*/
550BOOLEAN idHomIdeal (ideal id, ideal Q)
551{
552  int i;
553  BOOLEAN     b;
554  if ((id == NULL) || (IDELEMS(id) == 0)) return TRUE;
555  i = 0;
556  b = TRUE;
557  while ((i < IDELEMS(id)) && b)
558  {
559    b = pIsHomogeneous(id->m[i]);
560    i++;
561  }
562  if ((b) && (Q!=NULL) && (IDELEMS(Q)>0))
563  {
564    i=0;
565    while ((i < IDELEMS(Q)) && b)
566    {
567      b = pIsHomogeneous(Q->m[i]);
568      i++;
569    }
570  }
571  return b;
572}
573
574/*2
575*returns a minimized set of generators of h1
576*/
577ideal idMinBase (ideal h1)
578{
579  ideal h2, h3,h4,e;
580  int j,k;
581  int i,l,ll;
582  intvec * wth;
583  BOOLEAN homog;
584
585  homog = idHomModule(h1,currQuotient,&wth);
586  if ((currRing->OrdSgn == 1) && (!homog))
587  {
588    Warn("minbase applies only to the local or homogeneous case");
589    e=idCopy(h1);
590    return e;
591  }
592  if ((currRing->OrdSgn == 1) && (homog))
593  {
594    lists re=min_std(h1,currQuotient,(tHomog)homog,&wth,NULL,0,2);
595    h2 = (ideal)re->m[1].data;
596    re->m[1].data = NULL;
597    re->m[1].rtyp = NONE;
598    re->Clean();
599    return h2;
600  }
601  e=idInit(1,h1->rank);
602  if (idIs0(h1))
603  {
604    return e;
605  }
606  pEnlargeSet(&(e->m),IDELEMS(e),15);
607  IDELEMS(e) = 16;
608  h2 = std(h1,currQuotient,isNotHomog,NULL);
609  h3 = idMaxIdeal();
610  h4=idMult(h2,h3);
611  idDelete(&h3);
612  h3=std(h4,currQuotient,isNotHomog,NULL);
613  k = IDELEMS(h3);
614  while ((k > 0) && (h3->m[k-1] == NULL)) k--;
615  j = -1;
616  l = IDELEMS(h2);
617  while ((l > 0) && (h2->m[l-1] == NULL)) l--;
618  for (i=l-1; i>=0; i--)
619  {
620    if (h2->m[i] != NULL)
621    {
622      ll = 0;
623      while ((ll < k) && ((h3->m[ll] == NULL)
624      || !pDivisibleBy(h3->m[ll],h2->m[i])))
625        ll++;
626      if (ll >= k)
627      {
628        j++;
629        if (j > IDELEMS(e)-1)
630        {
631          pEnlargeSet(&(e->m),IDELEMS(e),16);
632          IDELEMS(e) += 16;
633        }
634        e->m[j] = pCopy(h2->m[i]);
635      }
636    }
637  }
638  idDelete(&h2);
639  idDelete(&h3);
640  idDelete(&h4);
641  if (currQuotient!=NULL)
642  {
643    h3=idInit(1,e->rank);
644    h2=kNF(h3,currQuotient,e);
645    idDelete(&h3);
646    idDelete(&e);
647    e=h2;
648  }
649  idSkipZeroes(e);
650  return e;
651}
652
653/*2
654*the minimal index of used variables - 1
655*/
656int pLowVar (poly p)
657{
658  int k,l,lex;
659
660  if (p == NULL) return -1;
661
662  k = 32000;/*a very large dummy value*/
663  while (p != NULL)
664  {
665    l = 1;
666    lex = pGetExp(p,l);
667    while ((l <= pVariables) && (lex == 0))
668    {
669      l++;
670      lex = pGetExp(p,l);
671    }
672    l--;
673    if (l < k) k = l;
674    pIter(p);
675  }
676  return k;
677}
678
679/*3
680*multiplies p with t (!cas) or  (t-1)
681*the index of t is:1, so we have to shift all variables
682*p is NOT in the actual ring, it has no t
683*/
684static poly pMultWithT (poly p,BOOLEAN cas)
685{
686  /*qp is the working pointer in p*/
687  /*result is the result, qresult is the working pointer*/
688  /*pp is p in the actual ring(shifted), qpp the working pointer*/
689  poly result,qp,pp;
690  poly qresult=NULL;
691  poly qpp=NULL;
692  int  i,j,lex;
693  number n;
694
695  pp = NULL;
696  result = NULL;
697  qp = p;
698  while (qp != NULL)
699  {
700    i = 0;
701    if (result == NULL)
702    {/*first monomial*/
703      result = pInit();
704      qresult = result;
705    }
706    else
707    {
708      qresult->next = pInit();
709      pIter(qresult);
710    }
711    for (j=pVariables-1; j>0; j--)
712    {
713      lex = pGetExp(qp,j);
714      pSetExp(qresult,j+1,lex);/*copy all variables*/
715    }
716    lex = pGetComp(qp);
717    pSetComp(qresult,lex);
718    n=nCopy(pGetCoeff(qp));
719    pSetCoeff0(qresult,n);
720    qresult->next = NULL;
721    pSetm(qresult);
722    /*qresult is now qp brought into the actual ring*/
723    if (cas)
724    { /*case: mult with t-1*/
725      pSetExp(qresult,1,0);
726      pSetm(qresult);
727      if (pp == NULL)
728      { /*first monomial*/
729        pp = pCopy(qresult);
730        qpp = pp;
731      }
732      else
733      {
734        qpp->next = pCopy(qresult);
735        pIter(qpp);
736      }
737      pGetCoeff(qpp)=nNeg(pGetCoeff(qpp));
738      /*now qpp contains -1*qp*/
739    }
740    pSetExp(qresult,1,1);/*this is mult. by t*/
741    pSetm(qresult);
742    pIter(qp);
743  }
744  /*
745  *now p is processed:
746  *result contains t*p
747  * if cas: pp contains -1*p (in the new ring)
748  */
749  if (cas)  qresult->next = pp;
750  /*  else      qresult->next = NULL;*/
751  return result;
752}
753
754/*3
755*deletes the place of t in p (t: variable with index 1)
756*p is NOT in the actual ring: it has pVariables+1 variables
757*/
758static poly pDivByT (poly * p,int size)
759{
760
761  poly result=NULL,
762       resultp=NULL , /** working pointer in result **/
763       pp;
764  int  i,j;
765
766  while (*p != NULL)
767  {
768    i = 0;
769    if (result == NULL)
770    {/*the first monomial*/
771      result = pInit();
772      resultp = result;
773      resultp->next = NULL;
774    }
775    else
776    {
777      resultp->next = pInit();
778      pIter(resultp);
779      resultp->next = NULL;
780    }
781    for (j=1; j<=pVariables; j++)
782    {
783      pSetExp(resultp,j,pGetExp(*p,j+1));
784    }
785    pSetComp(resultp,pGetComp(*p));
786    pSetCoeff0(resultp,pGetCoeff(*p));
787    pSetm(resultp);
788    pp = (*p)->next;
789    Free((ADDRESS)*p,size);
790    *p = pp;
791  }
792  return result;
793}
794
795/*2
796*dehomogenized the generators of the ideal id1 with the leading
797*monomial of p replaced by n
798*/
799ideal idDehomogen (ideal id1,poly p,number n)
800{
801  int i;
802  ideal result;
803
804  if (idIs0(id1))
805  {
806    return idInit(1,id1->rank);
807  }
808  result=idInit(IDELEMS(id1),id1->rank);
809  for (i=0; i<IDELEMS(id1); i++)
810  {
811    result->m[i] = pDehomogen(id1->m[i],p,n);
812  }
813  return result;
814}
815
816/*2
817* verschiebt die Indizees der Modulerzeugenden um i
818*/
819void pShift (poly * p,int i)
820{
821  poly qp1 = *p,qp2 = *p;/*working pointers*/
822  int     j = pMaxComp(*p),k = pMinComp(*p);
823
824  if (j+i < 0) return ;
825  while (qp1 != NULL)
826  {
827    if ((pGetComp(qp1)+i > 0) || ((j == -i) && (j == k)))
828    {
829      pSetComp(qp1,pGetComp(qp1)+i);
830      qp2 = qp1;
831      pIter(qp1);
832    }
833    else
834    {
835      if (qp2 == *p)
836      {
837        pIter(*p);
838        qp2->next = NULL;
839        pDelete(&qp2);
840        qp2 = *p;
841        qp1 = *p;
842      }
843      else
844      {
845        qp2->next = qp1->next;
846        qp1->next = NULL;
847        pDelete(&qp1);
848        qp1 = qp2->next;
849      }
850    }
851  }
852}
853
854/*2
855*initialized a field with r numbers between beg and end for the
856*procedure idNextChoise
857*/
858void idInitChoise (int r,int beg,int end,BOOLEAN  *endch,int * choise)
859{
860  /*returns the first choise of r numbers between beg and end*/
861  int i;
862  for (i=0; i<r; i++)
863  {
864    choise[i] = 0;
865  }
866  if (r <= end-beg+1)
867    for (i=0; i<r; i++)
868    {
869      choise[i] = beg+i;
870    }
871  if (r > end-beg+1)
872    *endch = TRUE;
873  else
874    *endch = FALSE;
875}
876
877/*2
878*returns the next choise of r numbers between beg and end
879*/
880void idGetNextChoise (int r,int end,BOOLEAN  *endch,int * choise)
881{
882  int i = r-1,j;
883  while ((i >= 0) && (choise[i] == end))
884  {
885    i--;
886    end--;
887  }
888  if (i == -1)
889    *endch = TRUE;
890  else
891  {
892    choise[i]++;
893    for (j=i+1; j<r; j++)
894    {
895      choise[j] = choise[i]+j-i;
896    }
897    *endch = FALSE;
898  }
899}
900
901/*2
902*takes the field choise of d numbers between beg and end, cancels the t-th
903*entree and searches for the ordinal number of that d-1 dimensional field
904* w.r.t. the algorithm of construction
905*/
906int idGetNumberOfChoise(int t, int d, int begin, int end, int * choise)
907{
908  int * localchoise,i,result=0;
909  BOOLEAN b=FALSE;
910
911  if (d<=1) return 1;
912  localchoise=(int*)Alloc((d-1)*sizeof(int));
913  idInitChoise(d-1,begin,end,&b,localchoise);
914  while (!b)
915  {
916    result++;
917    i = 0;
918    while ((i<t) && (localchoise[i]==choise[i])) i++;
919    if (i>=t)
920    {
921      i = t+1;
922      while ((i<d) && (localchoise[i-1]==choise[i])) i++;
923      if (i>=d)
924      {
925        Free((ADDRESS)localchoise,(d-1)*sizeof(int));
926        return result;
927      }
928    }
929    idGetNextChoise(d-1,end,&b,localchoise);
930  }
931  Free((ADDRESS)localchoise,(d-1)*sizeof(int));
932  return 0;
933}
934
935/*2
936*computes the binomial coefficient
937*/
938int binom (int n,int r)
939{
940  int i,result;
941
942  if (r==0) return 1;
943  if (n-r<r) return binom(n,n-r);
944  result = n-r+1;
945  for (i=2;i<=r;i++)
946  {
947    result *= n-r+i;
948    result /= i;
949  }
950  return result;
951}
952
953/*2
954*the free module of rank i
955*/
956ideal idFreeModule (int i)
957{
958  int j;
959  ideal h;
960
961  h=idInit(i,i);
962  for (j=0; j<i; j++)
963  {
964    h->m[j] = pOne();
965    pSetComp(h->m[j],j+1);
966  }
967  return h;
968}
969
970/*2
971* h3 := h1 intersect h2
972*/
973ideal idSect (ideal h1,ideal h2)
974{
975  ideal first=h2,second=h1,temp,temp1,result;
976  int i,j,k,flength,slength,length,rank=min(h1->rank,h2->rank);
977  intvec *w;
978  poly p,q;
979
980  if ((idIs0(h1)) && (idIs0(h2)))  return idInit(1,rank);
981  if (IDELEMS(h1)<IDELEMS(h2))
982  {
983    first = h1;
984    second = h2;
985  }
986  flength = idRankFreeModule(first);
987  slength = idRankFreeModule(second);
988  length  = max(flength,slength);
989  if (length==0)
990  {
991    length = 1;
992  }
993  temp = idInit(IDELEMS(first),1);
994  j = IDELEMS(temp);
995  while ((j>0) && (first->m[j-1]==NULL)) j--;
996  k = 0;
997  for (i=0;i<j;i++)
998  {
999    if (first->m[i]!=NULL)
1000    {
1001      temp->m[k] = pCopy(first->m[i]);
1002      q = pOne();
1003      pSetComp(q,i+1+length);
1004      if (flength==0) pShift(&(temp->m[k]),1);
1005      p = temp->m[k];
1006      while (pNext(p)) pIter(p);
1007      pNext(p) = q;
1008      k++;
1009    }
1010  }
1011  pEnlargeSet(&(temp->m),IDELEMS(temp),j+IDELEMS(second)-IDELEMS(temp));
1012  IDELEMS(temp) = j+IDELEMS(second);
1013  for (i=0;i<IDELEMS(second);i++)
1014  {
1015    if (second->m[i]!=NULL)
1016    {
1017      temp->m[k] = pCopy(second->m[i]);
1018      if (slength==0) pShift(&(temp->m[k]),1);
1019      k++;
1020    }
1021  }
1022  pSetSyzComp(length);
1023  temp1 = std(temp,currQuotient,testHomog,&w,NULL,length);
1024  if (w!=NULL) delete w;
1025  pSetSyzComp(0);
1026  idDelete(&temp);
1027  result = idInit(IDELEMS(temp1),rank);
1028  j = 0;
1029  for (i=0;i<IDELEMS(temp1);i++)
1030  {
1031    if ((temp1->m[i]!=NULL)
1032    && (pGetComp(temp1->m[i])>length))
1033    {
1034      p = temp1->m[i];
1035//PrintS("die Syzygie ist: ");pWrite(p);
1036      temp1->m[i] = NULL;
1037      while (p!=NULL)
1038      {
1039        q = pNext(p);
1040        pNext(p) = NULL;
1041        k = pGetComp(p)-1-length;
1042        pSetComp(p,0);
1043//PrintS("das %d-te Element: ",k);pWrite(first->m[k]);
1044        result->m[j] = pAdd(result->m[j],pMult(pCopy(first->m[k]),p));
1045        p = q;
1046      }
1047//PrintS("Generator ist: ");pWrite(result->m[j]);
1048      j++;
1049    }
1050  }
1051  idSkipZeroes(result);
1052  return result;
1053}
1054
1055/*2
1056* ideal/module intersection for a list of objects
1057* given as 'resolvente'
1058*/
1059ideal idMultSect(resolvente arg, int length)
1060{
1061  int i,j=0,k=0,syzComp,l,maxrk=-1,realrki;
1062  ideal bigmat,tempstd,result;
1063  poly p;
1064  int isIdeal=0;
1065  intvec * w=NULL;
1066
1067  /* find 0-ideals and max rank -----------------------------------*/
1068  for (i=0;i<length;i++)
1069  {
1070    if (!idIs0(arg[i]))
1071    {
1072      realrki=idRankFreeModule(arg[i]);
1073      k++;
1074      j += IDELEMS(arg[i]);
1075      if (realrki>maxrk) maxrk = realrki;
1076    }
1077    else
1078    {
1079      if (arg[i]!=NULL)
1080      {
1081        return idInit(1,arg[i]->rank);
1082      }
1083    }
1084  }
1085  if (maxrk == 0)
1086  {
1087    isIdeal = 1;
1088    maxrk = 1;
1089  }
1090  /* init -----------------------------------------------------------*/
1091  j += maxrk;
1092  bigmat = idInit(j,(k+1)*maxrk);
1093  syzComp = k*maxrk;
1094  pSetSyzComp(syzComp);
1095  /* create unit matrices ------------------------------------------*/
1096  for (i=0;i<maxrk;i++)
1097  {
1098    for (j=0;j<=k;j++)
1099    {
1100      p = pOne();
1101      pSetComp(p,i+1+j*maxrk);
1102      pSetm(p);
1103      bigmat->m[i] = pAdd(bigmat->m[i],p);
1104    }
1105  }
1106  /* enter given ideals ------------------------------------------*/
1107  i = maxrk;
1108  k = 0;
1109  for (j=0;j<length;j++)
1110  {
1111    if (arg[j]!=NULL)
1112    {
1113      for (l=0;l<IDELEMS(arg[j]);l++)
1114      {
1115        if (arg[j]->m[l]!=NULL)
1116        {
1117          bigmat->m[i] = pCopy(arg[j]->m[l]);
1118          pShift(&(bigmat->m[i]),k*maxrk+isIdeal);
1119          i++;
1120        }
1121      }
1122      k++;
1123    }
1124  }
1125  /* std computation --------------------------------------------*/
1126  tempstd = std(bigmat,currQuotient,testHomog,&w,NULL,syzComp);
1127  if (w!=NULL) delete w;
1128  idDelete(&bigmat);
1129  pSetSyzComp(0);
1130  /* interprete result ----------------------------------------*/
1131  result = idInit(8,maxrk);
1132  k = 0;
1133  for (j=0;j<IDELEMS(tempstd);j++)
1134  {
1135    if ((tempstd->m[j]!=NULL) && (pGetComp(tempstd->m[j])>syzComp))
1136    {
1137      if (k>=IDELEMS(result))
1138      {
1139        pEnlargeSet(&(result->m),IDELEMS(result),8);
1140        IDELEMS(result) += 8;
1141      }
1142      p = tempstd->m[j];
1143      tempstd->m[j] = NULL;
1144      pShift(&p,-syzComp-isIdeal);
1145      result->m[k] = p;
1146      k++;
1147    }
1148  }
1149  /* clean up ----------------------------------------------------*/
1150  idDelete(&tempstd);
1151  idSkipZeroes(result);
1152  return result;
1153}
1154
1155/*2
1156*computes the rank of the free module in which h1 embeds
1157*/
1158int lengthFreeModule (ideal  h1)
1159{
1160  int     i,j,k;
1161
1162  if (idIs0(h1)) return 0;
1163  k = -1;
1164  for (i=0; i<IDELEMS(h1); i++)
1165  {
1166    j = pMaxComp(h1->m[i]);
1167    if (j>k) k = j;
1168  }
1169  return k;
1170}
1171
1172/*2
1173*computes syzygies of h1,
1174*if quot != NULL it computes in the quotient ring modulo "quot"
1175*/
1176ideal  idPrepare (ideal  h1,ideal  quot, tHomog h,
1177  int* syzcomp, int *quotgen, int *quotdim, intvec **w)
1178{
1179  ideal   h2, h3;
1180  int     i;
1181  int     j,k;
1182  poly    p,q;
1183  BOOLEAN orderChanged=FALSE;
1184
1185  if (idIs0(h1)) return NULL;
1186  k = lengthFreeModule(h1);
1187  if (*syzcomp<k) *syzcomp = k;
1188  h2 = NULL;
1189  h2=idCopy(h1);
1190  i = IDELEMS(h2)-1;
1191  //while ((i >= 0) && (h2->m[i] == NULL)) i--;
1192  if (k == 0)
1193  {
1194    for (j=0; j<=i; j++) pShift(&(h2->m[j]),1);
1195    *syzcomp = 1;
1196  }
1197  h2->rank = *syzcomp+i+1;
1198  for (j=0; j<=i; j++)
1199  {
1200    p = h2->m[j];
1201    q = pOne();
1202#ifdef DRING
1203    if (pDRING) { pdDFlag(q)=1; pSetm(q); }
1204#endif
1205    pSetComp(q,*syzcomp+1+j);
1206    if (p!=NULL)
1207    {
1208      while (pNext(p)) pIter(p);
1209      p->next = q;
1210    }
1211    else
1212      h2->m[j]=q;
1213  }
1214  if ((pOrdSgn!=1) || (h!=TRUE) || (*syzcomp>1)|| (pLexOrder))
1215  {
1216    pSetSyzComp(*syzcomp);
1217    orderChanged=TRUE;
1218  }
1219  else
1220  {
1221    for(j=0;(j<IDELEMS(h2) && (!orderChanged));j++)
1222    {
1223      if (h2->m[j] != NULL)
1224      {
1225        p=h2->m[j];
1226        while ((p!=NULL)
1227        && (pGetComp(p)<=*syzcomp))
1228        {
1229          if (pIsConstantComp(p))
1230          {
1231            pSetSyzComp(*syzcomp);
1232            orderChanged=TRUE;
1233            break;
1234          }
1235          pIter(p);
1236        }
1237      }
1238    }
1239  }
1240//  if (orderChanged) PrintS("order changed\n");
1241//  else PrintS("order not changed\n");
1242#ifdef PDEBUG
1243  for(j=0;j<IDELEMS(h2);j++) pTest(h2->m[j]);
1244#endif
1245  h3=std(h2,quot,h,w,NULL,*syzcomp);
1246  //h3->rank = h2->rank; done by std -> initBuchMora -> initS
1247  h3->rank-=*syzcomp;
1248  idDelete(&h2);
1249  if (orderChanged) pSetSyzComp(0);
1250  return h3;
1251}
1252
1253ideal idSyzygies (ideal  h1,ideal  quot, tHomog h,intvec **w)
1254{
1255  int d;
1256  return idSyzygies(h1,quot,h,w,FALSE,d);
1257}
1258
1259ideal idSyzygies (ideal  h1,ideal  quot, tHomog h,intvec **w,
1260                  BOOLEAN setRegularity, int &deg)
1261{
1262  ideal e, h3;
1263  poly  p;
1264  int   i, j, k=0, quotdim, quotgen,length=0,reg;
1265  BOOLEAN isMonomial=TRUE;
1266
1267#ifdef PDEBUG
1268  int ii;
1269  for(ii=0;ii<IDELEMS(h1);ii++) pTest(h1->m[ii]);
1270  if (quot!=NULL)
1271  {
1272    for(ii=0;ii<IDELEMS(quot);ii++) pTest(quot->m[ii]);
1273  }
1274#endif
1275  if (idIs0(h1))
1276    return idFreeModule(IDELEMS(h1));
1277  h3=idPrepare(h1,quot,h,&k,&quotgen,&quotdim,w);
1278  if (h3==NULL)
1279    return idFreeModule(IDELEMS(h1));
1280  i = -1;
1281  e=idInit(16,h3->rank);
1282  for (j=0; j<IDELEMS(h3); j++)
1283  {
1284    if (h3->m[j] != NULL)
1285    {
1286      if (pMinComp(h3->m[j]) > k)
1287      {
1288        p = h3->m[j];
1289        h3->m[j]=NULL;
1290        pShift(&p,-k);
1291        if (p!=NULL)
1292        {
1293          i++;
1294          if (i+1 >= IDELEMS(e))
1295          {
1296            pEnlargeSet(&(e->m),IDELEMS(e),16);
1297            IDELEMS(e) += 16;
1298          }
1299          e->m[i] = p;
1300        }
1301      }
1302      else
1303      {
1304        isMonomial=isMonomial && (pNext(h3->m[j])==NULL);
1305        pDelete(&pNext(h3->m[j]));
1306      }
1307    }
1308  }
1309  if ((!isMonomial) && (!TEST_OPT_NOTREGULARITY) && (setRegularity) && (h==isHomog))
1310  {
1311    idSkipZeroes(h3);
1312    resolvente res = sySchreyerResolvente(h3,-1,&length,TRUE);
1313    intvec * dummy = syBetti(res,length,&reg, *w);
1314    deg = reg+2;
1315    delete dummy;
1316    for (j=0;j<length;j++)
1317    {
1318      if (res[j]!=NULL) idDelete(&(res[j]));
1319    }
1320    Free((ADDRESS)res,length*sizeof(ideal));
1321  }
1322  idDelete(&h3);
1323  idSkipZeroes(e);
1324  return e;
1325}
1326
1327/*2
1328* computes syzygies and minimizes the input (h1)
1329* ONLY used in syMinRes
1330*/
1331ideal idSyzMin (ideal h1,ideal  quot, tHomog h,intvec **w,
1332                  BOOLEAN setRegularity, int &deg)
1333{
1334  ideal e, h3;
1335  poly  p;
1336  intvec * reord;
1337  int   i, l=0, j, k=0, quotdim, quotgen,length=0,reg;
1338  BOOLEAN isMonomial=TRUE;
1339
1340  if (idIs0(h1))
1341    return idFreeModule(1);
1342//PrintS("h1 vorher\n");
1343//for (i=0;i<IDELEMS(h1);i++)
1344//{
1345//Print("Element %d: ",i);pWrite(h1->m[i]);
1346//}
1347  idSkipZeroes(h1);
1348  h3=idPrepare(h1,quot,h,&k,&quotgen,&quotdim,w);
1349//PrintS("h1 nachher\n");
1350//for (i=0;i<IDELEMS(h3);i++)
1351//{
1352//Print("Element %d: ",i);pWrite(h3->m[i]);
1353//}
1354  if (h3==NULL)
1355    return idFreeModule(1);
1356  for (i=IDELEMS(h1);i!=0;i--)
1357    pDelete(&(h1->m[i-1]));
1358  reord = new intvec(IDELEMS(h1)+1);
1359  i = -1;
1360  e=idInit(16,h3->rank);
1361  for (j=0; j<IDELEMS(h3); j++)
1362  {
1363    if (h3->m[j] != NULL)
1364    {
1365      p = h3->m[j];
1366      if (pGetComp(p) > k)
1367      {
1368        h3->m[j]=NULL;
1369        pShift(&p,-k);
1370        if (p!=NULL)
1371        {
1372          i++;
1373          if (i+1 >= IDELEMS(e))
1374          {
1375            pEnlargeSet(&(e->m),IDELEMS(e),16);
1376            IDELEMS(e) += 16;
1377          }
1378          e->m[i] = p;
1379        }
1380      }
1381      else
1382      {
1383        while ((pNext(p)!=NULL) && (pGetComp(pNext(p))<=k)) pIter(p);
1384        if (pIsConstantComp(pNext(p)))
1385        {
1386          (*reord)[pGetComp(pNext(p))-k] = l+1;
1387//Print("Element %d mit Comp %d: ",j,pGetComp(pNext(p))-k);
1388//pWrite(h3->m[j]);
1389          h1->m[l] = h3->m[j];
1390          h3->m[j] = pCopy(h1->m[l]);
1391          pDelete(&pNext(p));
1392          l++;
1393        }
1394        isMonomial=isMonomial && (pGetComp(pNext(h3->m[j]))>k);
1395        pDelete(&pNext(h3->m[j]));
1396      }
1397    }
1398  }
1399  if ((!isMonomial) && (!TEST_OPT_NOTREGULARITY) && (setRegularity) && (h==isHomog))
1400  {
1401    idSkipZeroes(h3);
1402    resolvente res = sySchreyerResolvente(h3,0,&length);
1403    intvec * dummy = syBetti(res,length,&reg, *w);
1404    deg = reg+2;
1405    delete dummy;
1406    for (j=0;j<length;j++)
1407    {
1408      if (res[j]!=NULL) idDelete(&(res[j]));
1409    }
1410    Free((ADDRESS)res,length*sizeof(ideal));
1411  }
1412  idDelete(&h3);
1413//PrintS("Komponententransformation: ");
1414//reord->show();
1415//PrintLn();
1416  for (i=IDELEMS(e);i!=0;i--)
1417  {
1418    if (e->m[i-1]!=NULL)
1419    {
1420      p = e->m[i-1];
1421      while (p!=NULL)
1422      {
1423        j = pGetComp(p);
1424        pSetComp(p,(*reord)[j]);
1425        pIter(p);
1426      }
1427      e->m[i-1] = pOrdPolyMerge(e->m[i-1]);
1428    }
1429  }
1430  idSkipZeroes(e);
1431  delete reord;
1432  return e;
1433}
1434
1435/*
1436*computes a standard basis for h1 and stores the transformation matrix
1437* in ma
1438*/
1439ideal idLiftStd (ideal  h1,ideal  quot, matrix* ma, tHomog h)
1440{
1441  int   i, j, k=0, t, quotgen, inputIsIdeal=lengthFreeModule(h1);
1442  ideal e, h3;
1443  poly  p=NULL, q, qq;
1444  intvec *w=NULL;
1445
1446  idDelete((ideal*)ma);
1447  *ma=mpNew(1,0);
1448  if (idIs0(h1))
1449    return idInit(1,h1->rank);
1450  h3=idPrepare(h1,quot,h,&k,&quotgen,&i,&w);
1451  if (w!=NULL) delete w;
1452  i = 0;
1453  for (j=0;j<IDELEMS(h3);j++)
1454  {
1455    if ((h3->m[j] != NULL) && (pMinComp(h3->m[j]) <= k))
1456      i++;
1457  }
1458  j = IDELEMS(h1);
1459  idDelete((ideal*)ma);
1460  *ma = mpNew(j,i);
1461  i = -1;
1462  e=idInit(16,h1->rank);
1463  for (j=0; j<IDELEMS(h3); j++)
1464  {
1465    if ((h3->m[j] != NULL) && (pMinComp(h3->m[j]) <= k))
1466    {
1467      q = h3->m[j];
1468      qq=q;
1469      h3->m[j]=NULL;
1470      while (pNext(q) != NULL)
1471      {
1472        if (pGetComp(pNext(q)) > k)
1473        {
1474          p = pNext(q);
1475          pNext(q) = pNext(pNext(q));
1476          pNext(p) = NULL;
1477          t=pGetComp(p);
1478          pSetComp(p,0);
1479          MATELEM(*ma,t-k,i+2) = pAdd(MATELEM(*ma,t-k,i+2),p);
1480        }
1481        else
1482        {
1483          pIter(q);
1484        }
1485      }
1486      if (!inputIsIdeal) pShift(&qq,-1);
1487      //if (quotgen != 0)
1488      //{
1489      //  q = kNF(quot,qq);
1490      //  pDelete(&qq);
1491      //}
1492      //else
1493        q = qq;
1494      if (q !=NULL)
1495      {
1496        i++;
1497        if (i+1 >= IDELEMS(e))
1498        {
1499          pEnlargeSet(&(e->m),IDELEMS(e),16);
1500          IDELEMS(e) += 16;
1501        }
1502        e->m[i] = q;
1503      }
1504    }
1505  }
1506  idDelete(&h3);
1507  idSkipZeroes(e);
1508  return e;
1509}
1510
1511/*2
1512*computes a representation of the generators of submod with respect to those
1513* of mod
1514*/
1515ideal idLiftNonStB (ideal  mod, ideal submod)
1516{
1517  int   lsmod =idRankFreeModule(submod), i, j, k, quotgen;
1518  ideal result, h3, temp;
1519
1520  if (idIs0(mod))
1521    return idInit(1,mod->rank);
1522
1523  if  ((idRankFreeModule(mod)!=0) && (lsmod==0)) lsmod=1;
1524  k=lsmod;
1525  //idSkipZeroes(mod);
1526  h3=idPrepare(mod,currQuotient,(tHomog)FALSE,&k,&quotgen,&i,NULL);
1527  for (j=0;j<IDELEMS(h3);j++)
1528  {
1529    if ((h3->m[j] != NULL) && (pMinComp(h3->m[j]) > k))
1530      pDelete(&(h3->m[j]));
1531  }
1532  idSkipZeroes(h3);
1533  if (lsmod==0)
1534  {
1535    temp = idCopy(submod);
1536    for (j=IDELEMS(temp);j>0;j--)
1537    {
1538      if (temp->m[j-1]!=NULL)
1539        pShift(&(temp->m[j-1]),1);
1540    }
1541  }
1542  else
1543  {
1544    temp = submod;
1545  }
1546  pSetSyzComp(k);
1547  result = kNF(h3,currQuotient,temp,k);
1548  result->rank = h3->rank;
1549  idDelete(&h3);
1550  if (lsmod==0)
1551    idDelete(&temp);
1552  pSetSyzComp(0);
1553  for (j=0;j<IDELEMS(result);j++)
1554  {
1555    if (result->m[j]!=NULL)
1556    {
1557      if (pGetComp(result->m[j])<=k)
1558      {
1559        WerrorS("2nd module lies not in the first");
1560        idDelete(&result);
1561        return idInit(1,1);
1562      }
1563      else
1564      {
1565        pShift(&(result->m[j]),-k);
1566        pNeg(result->m[j]);
1567      }
1568    }
1569  }
1570  return result;
1571}
1572
1573/*2
1574*computes a representation of the generators of submod with respect to those
1575* of mod which is given as standardbasis,
1576* uses currQuotient as the quotient ideal (if not NULL)
1577*/
1578ideal  idLift (ideal  mod,ideal submod)
1579{
1580  ideal temp, result;
1581  int   j,k;
1582  poly  p,q;
1583  BOOLEAN reported=FALSE;
1584
1585  if (idIs0(mod)) return idInit(1,mod->rank);
1586  k = lengthFreeModule(mod);
1587  temp=idCopy(mod);
1588  if (k == 0)
1589  {
1590    for (j=0; j<IDELEMS(temp); j++)
1591    {
1592      if (temp->m[j]!=NULL) pSetCompP(temp->m[j],1);
1593    }
1594    k = 1;
1595  }
1596  for (j=0; j<IDELEMS(temp); j++)
1597  {
1598    if (temp->m[j]!=NULL)
1599    {
1600      p = temp->m[j];
1601      q = pOne();
1602      pGetCoeff(q)=nNeg(pGetCoeff(q));   //set q to -1
1603      pSetComp(q,k+1+j);
1604      while (pNext(p)) pIter(p);
1605      pNext(p) = q;
1606    }
1607  }
1608  result=idInit(IDELEMS(submod),submod->rank);
1609  pSetSyzComp(k);
1610  for (j=0; j<IDELEMS(submod); j++)
1611  {
1612    if (submod->m[j]!=NULL)
1613    {
1614      p = pCopy(submod->m[j]);
1615      if (pGetComp(p)==0) pSetCompP(p,1);
1616      q = kNF(temp,currQuotient,p);
1617      pDelete(&p);
1618      if (q!=NULL)
1619      {
1620        if (pMinComp(q)<=k)
1621        {
1622          if (!reported)
1623          {
1624            Warn("first module not a standardbasis");
1625            Warn("or second not a proper submodule");
1626            reported=TRUE;
1627          }
1628          pDelete(&q);
1629        }
1630        else
1631        {
1632          pShift(&q,-k);
1633          result->m[j] = q;
1634        }
1635      }
1636    }
1637  }
1638  pSetSyzComp(0);
1639  //idSkipZeroes(result);
1640  idDelete(&temp);
1641  return result;
1642}
1643
1644/*2
1645*computes the quotient of h1,h2
1646*/
1647#ifdef OLD_QUOT
1648ideal idQuot (ideal  h1, ideal h2, BOOLEAN h1IsSB)
1649{
1650  int i,j = 0,l,ll,k,kkk,k1,k2,kmax;
1651  ideal h3,h4;
1652  poly     p,q = NULL;
1653  BOOLEAN     b = FALSE;
1654
1655  k1 = lengthFreeModule(h1);
1656  k2 = lengthFreeModule(h2);
1657  k=max(k1,k2);
1658  if (k==0) { k = 1; b=TRUE; }
1659  h4 = idInit(1,1);
1660  for (i=0; i<IDELEMS(h2); i++)
1661  {
1662    if (h2->m[i] != NULL)
1663    {
1664      p = pCopy(h2->m[i]);
1665      if (k2 == 0)
1666        pShift(&p,j*k+1);
1667      else
1668        pShift(&p,j*k);
1669      q = pAdd(q,p);
1670      j++;
1671    }
1672  }
1673  kmax = j*k+1;
1674  p = pOne();
1675  pSetComp(p,kmax);
1676  pSetSyzComp(kmax-1);
1677  q = pAdd(q,p);
1678  h4->m[0] = q;
1679  if (k2 == 0)
1680  {
1681    if (k > IDELEMS(h4))
1682    {
1683      pEnlargeSet(&(h4->m),IDELEMS(h4),k-IDELEMS(h4));
1684      IDELEMS(h4) = k;
1685    }
1686    for (i=1; i<k; i++)
1687    {
1688      p = pCopy(h4->m[i-1]);
1689      pShift(&p,1);
1690      h4->m[i] = p;
1691    }
1692  }
1693  kkk = IDELEMS(h4);
1694  i = IDELEMS(h1);
1695  while (h1->m[i-1] == NULL) i--;
1696  for (l=0; l<i; l++)
1697  {
1698    if(h1->m[l]!=NULL)
1699    {
1700      for (ll=0; ll<j; ll++)
1701      {
1702        p = pCopy(h1->m[l]);
1703        if (k1 == 0)
1704          pShift(&p,ll*k+1);
1705        else
1706          pShift(&p,ll*k);
1707        if (kkk >= IDELEMS(h4))
1708        {
1709          pEnlargeSet(&(h4->m),IDELEMS(h4),16);
1710          IDELEMS(h4) += 16;
1711        }
1712        h4->m[kkk] = p;
1713        kkk++;
1714      }
1715    }
1716  }
1717  h3 = std(h4,currQuotient,(tHomog)FALSE,NULL,NULL,kmax-1);
1718  pSetSyzComp(0);
1719  idDelete(&h4);
1720  for (i=0;i<IDELEMS(h3);i++)
1721  {
1722    if ((h3->m[i]!=NULL) && (pGetComp(h3->m[i])>=kmax))
1723    {
1724      if (b)
1725        pShift(&h3->m[i],-kmax);
1726      else
1727        pShift(&h3->m[i],-kmax+1);
1728    }
1729    else
1730      pDelete(&h3->m[i]);
1731  }
1732  if (b)
1733    h3->rank = 1;
1734  else
1735    h3->rank = h1->rank;
1736  h4=idCompactify(h3);
1737  idDelete(&h3);
1738  return h4;
1739}
1740#else
1741ideal idQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb)
1742{
1743  int i,j = 0,l,ll,k,kkk,k1,k2,kmax;
1744  intvec * weights,* weights1;
1745  ideal h3,h4,temph1;
1746  poly     p,q = NULL;
1747  BOOLEAN  resultIsIdeal = FALSE,addOnlyOne=TRUE;
1748  tHomog   hom=isNotHomog;
1749  BITSET old_test=test;
1750
1751  k1 = lengthFreeModule(h1);
1752  k2 = lengthFreeModule(h2);
1753  if (idIs0(h2))
1754  {
1755    if (k1==0)
1756    {
1757      h3 = idInit(1,1);
1758      h3->m[0] = pOne();
1759    }
1760    else
1761      h3 = idFreeModule(h1->rank);
1762    return h3;
1763  }
1764  k=max(k1,k2);
1765  if (k==0)
1766  {
1767    k = 1;
1768    resultIsIdeal=TRUE;
1769  }
1770  else if ((k1>0)&&(k2>0))
1771  {
1772    resultIsIdeal=TRUE;
1773  }
1774  hom = (tHomog)idHomModule(h1,currQuotient,&weights) ;
1775  h4 = idInit(1,1);
1776  for (i=0; i<IDELEMS(h2); i++)
1777  {
1778    if (h2->m[i] != NULL)
1779    {
1780      p = pCopy(h2->m[i]);
1781      if (k2 == 0)
1782        pShift(&p,j*k+1);
1783      else
1784        pShift(&p,j*k);
1785      q = pAdd(q,p);
1786      j++;
1787    }
1788  }
1789  kmax = j*k+1;
1790  p = pOne();
1791  pSetComp(p,kmax);
1792  pSetSyzComp(kmax-1);
1793  q = pAdd(q,p);
1794  h4->m[0] = q;
1795  if (k2 == 0)
1796  {
1797    if (k > IDELEMS(h4))
1798    {
1799      pEnlargeSet(&(h4->m),IDELEMS(h4),k-IDELEMS(h4));
1800      IDELEMS(h4) = k;
1801    }
1802    if (k>1) addOnlyOne = FALSE;
1803    for (i=1; i<k; i++)
1804    {
1805      p = pCopy(h4->m[i-1]);
1806      pShift(&p,1);
1807      h4->m[i] = p;
1808    }
1809  }
1810  kkk = IDELEMS(h4);
1811  if (addOnlyOne && (!h1IsStb))
1812    temph1 = std(h1,currQuotient,hom,&weights,NULL);
1813  else
1814    temph1 = h1;
1815  idTest(temph1);
1816  i = IDELEMS(temph1);
1817  while ((i>0) && (temph1->m[i-1]==NULL)) i--;
1818  for (l=0; l<i; l++)
1819  {
1820    if(temph1->m[l]!=NULL)
1821    {
1822      for (ll=0; ll<j; ll++)
1823      {
1824        p = pCopy(temph1->m[l]);
1825        if (k1 == 0)
1826          pShift(&p,ll*k+1);
1827        else
1828          pShift(&p,ll*k);
1829        if (kkk >= IDELEMS(h4))
1830        {
1831          pEnlargeSet(&(h4->m),IDELEMS(h4),16);
1832          IDELEMS(h4) += 16;
1833        }
1834        h4->m[kkk] = p;
1835        kkk++;
1836      }
1837    }
1838  }
1839  idTest(h4);
1840  if (addOnlyOne)
1841  {
1842    p = h4->m[0];
1843    for (i=0;i<IDELEMS(h4)-1;i++)
1844    {
1845      h4->m[i] = h4->m[i+1];
1846    }
1847    h4->m[IDELEMS(h4)-1] = p;
1848    idSkipZeroes(h4);
1849    test |= Sy_bit(OPT_SB_1);
1850  }
1851  idTest(h4);
1852  hom = (tHomog)idHomModule(h4,currQuotient,&weights1);
1853  if (addOnlyOne)
1854  {
1855    h3 = std(h4,currQuotient,hom,&weights1,NULL,kmax-1,IDELEMS(h4)-1);
1856  }
1857  else
1858  {
1859    h3 = std(h4,currQuotient,hom,&weights1,NULL,kmax-1);
1860  }
1861  idTest(h3);
1862  idDelete(&h4);
1863  pSetSyzComp(0);
1864  for (i=0;i<IDELEMS(h3);i++)
1865  {
1866    if ((h3->m[i]!=NULL) && (pGetComp(h3->m[i])>=kmax))
1867    {
1868      if (resultIsIdeal)
1869        pShift(&h3->m[i],-kmax);
1870      else
1871        pShift(&h3->m[i],-kmax+1);
1872    }
1873    else
1874      pDelete(&h3->m[i]);
1875  }
1876  if (resultIsIdeal)
1877    h3->rank = 1;
1878  else
1879    h3->rank = h1->rank;
1880  idTest(h3);
1881  idSkipZeroes(h3);
1882  //h4=idCompactify(h3);
1883  //idDelete(&h3);
1884  if ((addOnlyOne) && (!h1IsStb))
1885    idDelete(&temph1);
1886  if (weights!=NULL) delete weights;
1887  if (weights1!=NULL) delete weights1;
1888  test = old_test;
1889  return h3;
1890}
1891#endif
1892
1893/*2
1894*computes recursively all monomials of a certain degree
1895*in every step the actvar-th entry in the exponential
1896*vector is incremented and the other variables are
1897*computed by recursive calls of makemonoms
1898*if the last variable is reached, the difference to the
1899*degree is computed directly
1900*vars is the number variables
1901*actvar is the actual variable to handle
1902*deg is the degree of the monomials to compute
1903*monomdeg is the actual degree of the monomial in consideration
1904*/
1905static void makemonoms(int vars,int actvar,int deg,int monomdeg)
1906{
1907  poly p;
1908  int i=0;
1909
1910  if ((idpowerpoint == 0) && (actvar ==1))
1911  {
1912    idpower[idpowerpoint] = pOne();
1913    monomdeg = 0;
1914  }
1915  while (i<=deg)
1916  {
1917    if (deg == monomdeg)
1918    {
1919      pSetm(idpower[idpowerpoint]);
1920      idpowerpoint++;
1921      return;
1922    }
1923    if (actvar == vars)
1924    {
1925      pSetExp(idpower[idpowerpoint],actvar,deg-monomdeg);
1926      pSetm(idpower[idpowerpoint]);
1927      idpowerpoint++;
1928      return;
1929    }
1930    else
1931    {
1932      p = pCopy(idpower[idpowerpoint]);
1933      makemonoms(vars,actvar+1,deg,monomdeg);
1934      idpower[idpowerpoint] = p;
1935    }
1936    monomdeg++;
1937    pSetExp(idpower[idpowerpoint],actvar,pGetExp(idpower[idpowerpoint],actvar)+1);
1938    pSetm(idpower[idpowerpoint]);
1939    i++;
1940  }
1941}
1942
1943/*2
1944*returns the deg-th power of the maximal ideal of 0
1945*/
1946ideal idMaxIdeal(int deg)
1947{
1948  if (deg < 0)
1949  {
1950    WarnS("maxideal: power must be non-negative");
1951  }
1952  if (deg < 1)
1953  {
1954    ideal I=idInit(1,1);
1955    I->m[0]=pOne();
1956    return I;
1957  }
1958  if (deg == 1)
1959  {
1960    return idMaxIdeal();
1961  }
1962
1963  int vars = currRing->N;
1964  int i = binom(vars+deg-1,deg);
1965  ideal id=idInit(i,1);
1966  idpower = id->m;
1967  idpowerpoint = 0;
1968  makemonoms(vars,1,deg,0);
1969  idpower = NULL;
1970  idpowerpoint = 0;
1971  return id;
1972}
1973
1974/*2
1975*computes recursively all generators of a certain degree
1976*of the ideal "givenideal"
1977*elms is the number elements in the given ideal
1978*actelm is the actual element to handle
1979*deg is the degree of the power to compute
1980*gendeg is the actual degree of the generator in consideration
1981*/
1982static void makepotence(int elms,int actelm,int deg,int gendeg)
1983{
1984  poly p;
1985  int i=0;
1986
1987  if ((idpowerpoint == 0) && (actelm ==1))
1988  {
1989    idpower[idpowerpoint] = pOne();
1990    gendeg = 0;
1991  }
1992  while (i<=deg)
1993  {
1994    if (deg == gendeg)
1995    {
1996      idpowerpoint++;
1997      return;
1998    }
1999    if (actelm == elms)
2000    {
2001      p=pPower(pCopy(givenideal[actelm-1]),deg-gendeg);
2002      idpower[idpowerpoint]=pMult(idpower[idpowerpoint],p);
2003      idpowerpoint++;
2004      return;
2005    }
2006    else
2007    {
2008      p = pCopy(idpower[idpowerpoint]);
2009      makepotence(elms,actelm+1,deg,gendeg);
2010      idpower[idpowerpoint] = p;
2011    }
2012    gendeg++;
2013    idpower[idpowerpoint]=pMult(idpower[idpowerpoint],pCopy(givenideal[actelm-1]));
2014    i++;
2015  }
2016}
2017
2018/*2
2019*returns the deg-th power of the ideal gid
2020*/
2021//ideal idPower(ideal gid,int deg)
2022//{
2023//  int i;
2024//  ideal id;
2025//
2026//  if (deg < 1) deg = 1;
2027//  i = binom(IDELEMS(gid)+deg-1,deg);
2028//  id=idInit(i,1);
2029//  idpower = id->m;
2030//  givenideal = gid->m;
2031//  idpowerpoint = 0;
2032//  makepotence(IDELEMS(gid),1,deg,0);
2033//  idpower = NULL;
2034//  givenideal = NULL;
2035//  idpowerpoint = 0;
2036//  return id;
2037//}
2038static void idNextPotence(ideal given, ideal result,
2039  int begin, int end, int deg, int restdeg, poly ap)
2040{
2041  poly p;
2042  int i;
2043
2044  p = pPower(pCopy(given->m[begin]),restdeg);
2045  i = result->nrows;
2046  result->m[i] = pMult(pCopy(ap),p);
2047//PrintS(".");
2048  (result->nrows)++;
2049  if (result->nrows >= IDELEMS(result))
2050  {
2051    pEnlargeSet(&(result->m),IDELEMS(result),16);
2052    IDELEMS(result) += 16;
2053  }
2054  if (begin == end) return;
2055  for (i=restdeg-1;i>=0;i--)
2056  {
2057    p = pPower(pCopy(given->m[begin]),i);
2058    p = pMult(pCopy(ap),p);
2059    idNextPotence(given, result, begin+1, end, deg, restdeg-i, p);
2060    pDelete(&p);
2061  }
2062}
2063
2064ideal idPower(ideal given,int exp)
2065{
2066  ideal result,temp;
2067  poly p1;
2068  int i;
2069
2070  if (idIs0(given)) return idInit(1,1);
2071  temp = idCopy(given);
2072  idSkipZeroes(temp);
2073  i = binom(IDELEMS(temp)+exp-1,exp);
2074  result = idInit(i,1);
2075  result->nrows = 0;
2076//Print("ideal contains %d elements\n",i);
2077  p1=pOne();
2078  idNextPotence(temp,result,0,IDELEMS(temp)-1,exp,exp,p1);
2079  pDelete(&p1);
2080  idDelete(&temp);
2081  idSkipZeroes(result);
2082  result->nrows = 1;
2083  return result;
2084}
2085
2086/*2
2087* eliminate delVar (product of vars) in h1
2088*/
2089ideal idElimination (ideal h1,poly delVar,intvec *hilb)
2090{
2091  int    i,j=0,k,l;
2092  ideal  h,hh, h3;
2093  int    *ord,*block0,*block1;
2094  int    ordersize=2;
2095  short  **wv;
2096  tHomog hom;
2097  intvec * w;
2098  sip_sring tmpR;
2099  ring origR = currRing;
2100
2101  if (delVar==NULL)
2102  {
2103    return idCopy(h1);
2104  }
2105  if (currQuotient!=NULL)
2106  {
2107    WerrorS("cannot eliminate in a qring");
2108    return idCopy(h1);
2109  }
2110  if (idIs0(h1)) return idInit(1,h1->rank);
2111  hom=(tHomog)idHomModule(h1,NULL,&w); //sets w to weight vector or NULL
2112  h3=idInit(16,h1->rank);
2113  for (k=0;; k++)
2114  {
2115    if (currRing->order[k]!=0) ordersize++;
2116    else break;
2117  }
2118  ord=(int*)Alloc0(ordersize*sizeof(int));
2119  block0=(int*)Alloc(ordersize*sizeof(int));
2120  block1=(int*)Alloc(ordersize*sizeof(int));
2121  for (k=0;; k++)
2122  {
2123    if (currRing->order[k]!=0)
2124    {
2125      block0[k+1] = currRing->block0[k];
2126      block1[k+1] = currRing->block1[k];
2127      ord[k+1] = currRing->order[k];
2128    }
2129    else
2130      break;
2131  }
2132  block0[0] = 1;
2133  block1[0] = pVariables;
2134  wv=(short**) Alloc0(ordersize*sizeof(short**));
2135  memcpy4(wv+1,currRing->wvhdl,(ordersize-1)*sizeof(short**));
2136  wv[0]=(short*)Alloc0((pVariables+1)*sizeof(short));
2137  for (j=0;j<pVariables;j++)
2138    if (pGetExp(delVar,j+1)!=0) wv[0][j]=1;
2139  ord[0] = ringorder_a;
2140
2141  // fill in tmp ring to get back the data later on
2142  tmpR  = *origR;
2143  tmpR.order = ord;
2144  tmpR.block0 = block0;
2145  tmpR.block1 = block1;
2146  tmpR.wvhdl = wv;
2147  rComplete(&tmpR);
2148
2149  // change into the new ring
2150  pChangeRing(pVariables,currRing->OrdSgn,ord,block0,block1,wv);
2151  //  rChangeCurrRing(&tmpR, FALSE);
2152  currRing = &tmpR;
2153  h = idInit(IDELEMS(h1),1);
2154  // fetch data from the old ring
2155  for (k=0;k<IDELEMS(h1);k++) h->m[k] = pFetchCopy(origR, h1->m[k]);
2156  // compute std
2157  hh = std(h,NULL,hom,&w,hilb);
2158  idDelete(&h);
2159
2160  // go back to the original ring
2161  rChangeCurrRing(origR,FALSE);
2162  i = IDELEMS(hh)-1;
2163  while ((i >= 0) && (hh->m[i] == NULL)) i--;
2164  j = -1;
2165  // fetch data from temp ring
2166  for (k=0; k<=i; k++)
2167  {
2168    l=pVariables;
2169    while ((l>0) && (pRingGetExp(&tmpR, hh->m[k],l)*pGetExp(delVar,l)==0)) l--;
2170    if (l==0)
2171    {
2172      j++;
2173      if (j >= IDELEMS(h3))
2174      {
2175        pEnlargeSet(&(h3->m),IDELEMS(h3),16);
2176        IDELEMS(h3) += 16;
2177      }
2178      h3->m[j] = pFetchCopy(&tmpR, hh->m[k]);
2179    }
2180  }
2181  idDelete(&hh);
2182  idSkipZeroes(h3);
2183  Free((ADDRESS)wv[0],(pVariables+1)*sizeof(short));
2184  Free((ADDRESS)wv,ordersize*sizeof(short**));
2185  Free((ADDRESS)ord,ordersize*sizeof(int));
2186  Free((ADDRESS)block0,ordersize*sizeof(int));
2187  Free((ADDRESS)block1,ordersize*sizeof(int));
2188  if (w!=NULL)
2189    delete w;
2190  return h3;
2191}
2192
2193/*3
2194* produces recursively the ideal of all arxar-minors of a
2195*/
2196static void idRecMin(matrix a,int ar,poly *barDiv,ideal result,
2197              int * nextPlace, int rowToChose=0)
2198{
2199//Print("Level: %d\n",ar);
2200/*--- there is no non-zero minor coming from a------------*/
2201  if((ar<0) || (ar>min(a->ncols,a->nrows)-1))
2202  {
2203    return;
2204  }
2205
2206/*--- initializing temporary structures-------------------*/
2207  int i,j,r=rowToChose,c,newi,newp,k;
2208  poly p=NULL,pp;
2209
2210  if (ar==0)
2211  {
2212/*--- ar is 0, the matrix-entres are minors---------------*/
2213    for (i=a->nrows;i>0;i--)
2214    {
2215      for (j=a->ncols;j>0;j--)
2216      {
2217        p = MATELEM(a,i,j);
2218        if (p!=NULL)
2219        {
2220          if (*nextPlace>=IDELEMS(result))
2221          {
2222            pEnlargeSet(&(result->m),IDELEMS(result),16);
2223            IDELEMS(result) += 16;
2224          }
2225          MATELEM(a,i,j) = NULL;
2226          result->m[*nextPlace] = p;
2227          (*nextPlace)++;
2228        }
2229      }
2230    }
2231    idTest(result);
2232    idDelete((ideal*)&a);
2233    return;
2234  }
2235/*--- ar>0, we perform one step of the Bareiss algorithm--*/
2236  p = pCopy(*barDiv);   //we had to store barDiv for the remaining loops
2237  pp = pCopy(p);   //we had to store barDiv for the remaining loops
2238  matrix nextStep = mpOneStepBareiss(a,barDiv,&r,&c);
2239//Print("next row is: %d, next col: %d\n",r,c);
2240/*--- there is no pivot - the matrix is zero -------------*/
2241  if (r*c==0)
2242  {
2243    idDelete((ideal*)&a);
2244    return;
2245  }
2246/*--- we read out the r-1 x c-1 matrix for the next step--*/
2247  if ((a->nrows-1)*(a->ncols-1)>0)
2248  {
2249    matrix next = mpNew(a->nrows-1,a->ncols-1);
2250    for (i=a->nrows-1;i>0;i--)
2251    {
2252      for (j=a->ncols-1;j>0;j--)
2253      {
2254        MATELEM(next,i,j) = MATELEM(nextStep,i,j);
2255        MATELEM(nextStep,i,j) =NULL;
2256      }
2257    }
2258    idDelete((ideal*)&nextStep);
2259/*--- we call the next Step------------------------------*/
2260    idRecMin(next,ar-1,barDiv,result,nextPlace);
2261    next = NULL;
2262  }
2263/*--- now we have to take out the r-th row...------------*/
2264  if (((a->nrows)>1) && (rowToChose==0))
2265  {
2266    nextStep = mpNew(a->nrows-1,a->ncols);
2267    for (i=r-1;i>0;i--)
2268    {
2269      for (j=a->ncols;j>0;j--)
2270      {
2271        MATELEM(nextStep,i,j) = pCopy(MATELEM(a,i,j));
2272      }
2273    }
2274    for (i=a->nrows;i>r;i--)
2275    {
2276      for (j=a->ncols;j>0;j--)
2277      {
2278        MATELEM(nextStep,i-1,j) = pCopy(MATELEM(a,i,j));
2279      }
2280    }
2281/*--- and to perform the algorithm with the rest---------*/
2282    idRecMin(nextStep,ar,&p,result,nextPlace);
2283    nextStep = NULL;
2284  }
2285/*--- now we have to take out the c-th col...------------*/
2286  if ((a->nrows)>1)
2287  {
2288    nextStep = mpNew(a->nrows,a->ncols-1);
2289    for (i=a->nrows;i>0;i--)
2290    {
2291      for (j=c-1;j>0;j--)
2292      {
2293        MATELEM(nextStep,i,j) = MATELEM(a,i,j);
2294        MATELEM(a,i,j) = NULL;
2295      }
2296    }
2297    for (i=a->nrows;i>0;i--)
2298    {
2299      for (j=a->ncols;j>c;j--)
2300      {
2301        MATELEM(nextStep,i,j-1) = MATELEM(a,i,j);
2302        MATELEM(a,i,j) = NULL;
2303      }
2304    }
2305/*--- and to perform the algorithm with the rest---------*/
2306    idDelete((ideal*)&a);
2307    idRecMin(nextStep,ar,&p,result,nextPlace,r);
2308    nextStep = NULL;
2309  }
2310/*--- deleting temporary structures and returns----------*/
2311  pDelete(barDiv);
2312  pDelete(&p);
2313  pDelete(&pp);
2314  return;
2315}
2316
2317/*2
2318* compute all ar-minors of the matrix a
2319* the caller of idRecMin
2320*/
2321ideal idMinors(matrix a, int ar)
2322{
2323  if((ar<=0) || (ar>min(a->ncols,a->nrows)))
2324  {
2325    Werror("%d-th minor ",ar);
2326    return NULL;
2327  }
2328  int i=0;
2329  poly barDiv=NULL;
2330  ideal result=idInit(16,1);
2331  idTest(result);
2332
2333  idRecMin(mpCopy(a),ar-1,&barDiv,result,&i);
2334  idSkipZeroes(result);
2335  return result;
2336}
2337
2338/*2
2339*returns TRUE if p is a unit element in the current ring
2340*/
2341BOOLEAN pIsUnit(poly p)
2342{
2343  int i;
2344
2345  if (p == NULL) return FALSE;
2346  i = 1;
2347  while (i<=pVariables && pGetExp(p,i) == 0) i++;
2348  if (i > pVariables && (pGetComp(p) == 0))
2349  {
2350    if (currRing->OrdSgn == 1 && pNext(p) !=NULL) return FALSE;
2351    return TRUE;
2352  }
2353  return FALSE;
2354}
2355
2356/*2
2357*skips all zeroes and double elements, searches also for units
2358*/
2359ideal idCompactify(ideal id)
2360{
2361  ideal result = NULL;
2362  int i,j;
2363  BOOLEAN b=FALSE;
2364
2365  result=idCopy(id);
2366  i = IDELEMS(result)-1;
2367  while ((! b) && (i>=0))
2368  {
2369    b=pIsUnit(result->m[i]);
2370    i--;
2371  }
2372  if (b)
2373  {
2374    for (i=IDELEMS(result)-1;i>=0;i--)
2375      pDelete(&result->m[i]);
2376    result->m[0]=pOne();
2377  }
2378  else
2379  {
2380    for (i=1;i<IDELEMS(result);i++)
2381    {
2382      if (result->m[i]!=NULL)
2383      {
2384        for (j=0;j<i;j++)
2385        {
2386          if ((result->m[j]!=NULL)
2387          && (pComparePolys(result->m[i],result->m[j])))
2388          {
2389            pDelete(&(result->m[j]));
2390          }
2391        }
2392      }
2393    }
2394  }
2395  idSkipZeroes(result);
2396  return result;
2397}
2398
2399/*2
2400*returns TRUE if id1 is a submodule of id2
2401*/
2402BOOLEAN idIsSubModule(ideal id1,ideal id2)
2403{
2404  int i;
2405  poly p;
2406
2407  if (idIs0(id1)) return TRUE;
2408  for (i=0;i<IDELEMS(id1);i++)
2409  {
2410    if (id1->m[i] != NULL)
2411    {
2412      p = kNF(id2,currQuotient,id1->m[i]);
2413      if (p != NULL)
2414      {
2415        pDelete(&p);
2416        return FALSE;
2417      }
2418    }
2419  }
2420  return TRUE;
2421}
2422
2423/*2
2424* returns the ideals of initial terms
2425*/
2426ideal idHead(ideal h)
2427{
2428  ideal m = idInit(IDELEMS(h),h->rank);
2429  int i;
2430
2431  for (i=IDELEMS(h)-1;i>=0; i--)
2432  {
2433    if (h->m[i]!=NULL) m->m[i]=pHead(h->m[i]);
2434  }
2435  return m;
2436}
2437
2438ideal idHomogen(ideal h, int varnum)
2439{
2440  ideal m = idInit(IDELEMS(h),h->rank);
2441  int i;
2442
2443  for (i=IDELEMS(h)-1;i>=0; i--)
2444  {
2445    m->m[i]=pHomogen(h->m[i],varnum);
2446  }
2447  return m;
2448}
2449
2450/*------------------type conversions----------------*/
2451ideal idVec2Ideal(poly vec)
2452{
2453   ideal result=idInit(1,1);
2454   Free((ADDRESS)result->m,sizeof(poly));
2455   result->m=NULL; // remove later
2456   pVec2Polys(vec, &(result->m), &(IDELEMS(result)));
2457   return result;
2458}
2459
2460ideal idMatrix2Module(matrix mat)
2461{
2462  ideal result = idInit(MATCOLS(mat),MATROWS(mat));
2463  int i,j;
2464  poly h;
2465#ifdef DRING
2466  poly p;
2467#endif
2468
2469  for(j=0;j<MATCOLS(mat);j++) /* j is also index in result->m */
2470  {
2471    for (i=1;i<=MATROWS(mat);i++)
2472    {
2473      h = MATELEM(mat,i,j+1);
2474      if (h!=NULL)
2475      {
2476        MATELEM(mat,i,j+1)=NULL;
2477        pSetCompP(h,i);
2478#ifdef DRING
2479        pdSetDFlag(h,0);
2480#endif
2481        result->m[j] = pAdd(result->m[j],h);
2482      }
2483    }
2484  }
2485  return result;
2486}
2487
2488/*2
2489* converts a module into a matrix, destroyes the input
2490*/
2491matrix idModule2Matrix(ideal mod)
2492{
2493  matrix result = mpNew(mod->rank,IDELEMS(mod));
2494  int i,cp;
2495  poly p,h;
2496
2497  for(i=0;i<IDELEMS(mod);i++)
2498  {
2499    p=mod->m[i];
2500    mod->m[i]=NULL;
2501    while (p!=NULL)
2502    {
2503      h=p;
2504      pIter(p);
2505      pNext(h)=NULL;
2506//      cp = max(1,pGetComp(h));     // if used for ideals too
2507      cp = pGetComp(h);
2508      pSetComp(h,0);
2509#ifdef TEST
2510      if (cp>mod->rank)
2511      {
2512        Print("## inv. rank %d -> %d\n",mod->rank,cp);
2513        int k,l,o=mod->rank;
2514        mod->rank=cp;
2515        matrix d=mpNew(mod->rank,IDELEMS(mod));
2516        for (l=1; l<=o; l++)
2517        {
2518          for (k=1; k<=IDELEMS(mod); k++)
2519          {
2520            MATELEM(d,l,k)=MATELEM(result,l,k);
2521            MATELEM(result,l,k)=NULL;
2522          }
2523        }
2524        idDelete((ideal *)&result);
2525        result=d;
2526      }
2527#endif
2528      MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
2529    }
2530  }
2531  return result;
2532}
2533
2534matrix idModule2formatedMatrix(ideal mod,int rows, int cols)
2535{
2536  matrix result = mpNew(rows,cols);
2537  int i,cp,r=idRankFreeModule(mod),c=IDELEMS(mod);
2538  poly p,h;
2539
2540  if (r>rows) r = rows;
2541  if (c>cols) c = cols;
2542  for(i=0;i<c;i++)
2543  {
2544    p=mod->m[i];
2545    mod->m[i]=NULL;
2546    while (p!=NULL)
2547    {
2548      h=p;
2549      pIter(p);
2550      pNext(h)=NULL;
2551      cp = pGetComp(h);
2552      if (cp<=r)
2553      {
2554        pSetComp(h,0);
2555        MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
2556      }
2557      else
2558        pDelete(&h);
2559    }
2560  }
2561  idDelete(&mod);
2562  return result;
2563}
2564
2565/*2
2566* substitute the n-th variable by the monomial e in id
2567* destroy id
2568*/
2569ideal  idSubst(ideal id, int n, poly e)
2570{
2571  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
2572  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
2573
2574  res->rank = id->rank;
2575  for(k--;k>=0;k--)
2576  {
2577    res->m[k]=pSubst(id->m[k],n,e);
2578    id->m[k]=NULL;
2579  }
2580  idDelete(&id);
2581  return res;
2582}
2583
2584BOOLEAN idHomModule(ideal m, ideal Q, intvec **w)
2585{
2586  if (w!=NULL) *w=NULL;
2587  if ((Q!=NULL) && (!idHomIdeal(Q,NULL))) return FALSE;
2588  if (idIs0(m)) return TRUE;
2589
2590  int i,j,cmax=2,order=0,ord,* diff,* iscom,diffmin=32000;
2591  poly p=NULL;
2592  int length=IDELEMS(m);
2593  polyset P=m->m;
2594  polyset F=(polyset)Alloc(length*sizeof(poly));
2595  for (i=length-1;i>=0;i--)
2596  {
2597    p=F[i]=P[i];
2598    cmax=max(cmax,pMaxComp(p)+1);
2599  }
2600  diff = (int *)Alloc0(cmax*sizeof(int));
2601  if (w!=NULL) *w=new intvec(cmax-1);
2602  iscom = (int *)Alloc0(cmax*sizeof(int));
2603  i=0;
2604  while (i<=length)
2605  {
2606    if (i<length)
2607    {
2608      p=F[i];
2609      while ((p!=NULL) && (!iscom[pGetComp(p)])) pIter(p);
2610    }
2611    if ((p==NULL) && (i<length))
2612    {
2613      i++;
2614    }
2615    else
2616    {
2617      if (p==NULL)
2618      {
2619        i=0;
2620        while ((i<length) && (F[i]==NULL)) i++;
2621        if (i>=length) break;
2622        p = F[i];
2623      }
2624      if (pLexOrder)
2625        order=pTotaldegree(p);
2626      else
2627      //  order = p->order;
2628        order = pFDeg(p);
2629      order += diff[pGetComp(p)];
2630      p = F[i];
2631//Print("Actual p=F[%d]: ",i);pWrite(p);
2632      F[i] = NULL;
2633      i=0;
2634    }
2635    while (p!=NULL)
2636    {
2637      //if (pLexOrder)
2638      //  ord=pTotaldegree(p);
2639      //else
2640      //  ord = p->order;
2641      ord = pFDeg(p);
2642      if (!iscom[pGetComp(p)])
2643      {
2644        diff[pGetComp(p)] = order-ord;
2645        iscom[pGetComp(p)] = 1;
2646/*
2647*PrintS("new diff: ");
2648*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
2649*PrintLn();
2650*PrintS("new iscom: ");
2651*for (j=0;j<cmax;j++) Print("%d ",iscom[j]);
2652*PrintLn();
2653*Print("new set %d, order %d, ord %d, diff %d\n",pGetComp(p),order,ord,diff[pGetComp(p)]);
2654*/
2655      }
2656      else
2657      {
2658/*
2659*PrintS("new diff: ");
2660*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
2661*PrintLn();
2662*Print("order %d, ord %d, diff %d\n",order,ord,diff[pGetComp(p)]);
2663*/
2664        if (order != ord+diff[pGetComp(p)])
2665        {
2666          Free((ADDRESS) iscom,cmax*sizeof(int));
2667          Free((ADDRESS) diff,cmax*sizeof(int));
2668          Free((ADDRESS) F,length*sizeof(poly));
2669          delete *w;*w=NULL;
2670          return FALSE;
2671        }
2672      }
2673      pIter(p);
2674    }
2675  }
2676  Free((ADDRESS) iscom,cmax*sizeof(int));
2677  Free((ADDRESS) F,length*sizeof(poly));
2678  for (i=1;i<cmax;i++) (**w)[i-1]=diff[i];
2679  for (i=1;i<cmax;i++)
2680  {
2681    if (diff[i]<diffmin) diffmin=diff[i];
2682  }
2683  for (i=1;i<cmax;i++)
2684  {
2685    (**w)[i-1]=diff[i]-diffmin;
2686  }
2687  Free((ADDRESS) diff,cmax*sizeof(int));
2688  return TRUE;
2689}
2690
2691ideal idJet(ideal i,int d)
2692{
2693  ideal r=idInit(IDELEMS(i),i->rank);
2694  int k;
2695  for(k=0; k<IDELEMS(i); k++)
2696  {
2697    r->m[k]=pJet(i->m[k],d);
2698  }
2699  return r;
2700}
2701
2702ideal idJetW(ideal i,int d, intvec * iv)
2703{
2704  ideal r=idInit(IDELEMS(i),i->rank);
2705  if (ecartWeights!=NULL)
2706  {
2707    WerrorS("cannot compute weighted jets now");
2708  }
2709  else
2710  {
2711    short *w=iv2array(iv);
2712    int k;
2713    for(k=0; k<IDELEMS(i); k++)
2714    {
2715      r->m[k]=pJetW(i->m[k],d,w);
2716    }
2717    Free((ADDRESS)w,(pVariables+1)*sizeof(short));
2718  }
2719  return r;
2720}
2721
2722matrix idDiff(matrix i, int k)
2723{
2724  int e=MATCOLS(i)*MATROWS(i);
2725  matrix r=mpNew(MATROWS(i),MATCOLS(i));
2726  int j;
2727  for(j=0; j<e; j++)
2728  {
2729    r->m[j]=pDiff(i->m[j],k);
2730  }
2731  return r;
2732}
2733
2734matrix idDiffOp(ideal I, ideal J,BOOLEAN multiply)
2735{
2736  matrix r=mpNew(IDELEMS(I),IDELEMS(J));
2737  int i,j;
2738  for(i=0; i<IDELEMS(I); i++)
2739  {
2740    for(j=0; j<IDELEMS(J); j++)
2741    {
2742      MATELEM(r,i+1,j+1)=pDiffOp(I->m[i],J->m[j],multiply);
2743    }
2744  }
2745  return r;
2746}
2747
2748/*2
2749* represents (h1+h2)/h2=h1/(h1 intersect h2)
2750*/
2751ideal idModulo (ideal h2,ideal h1)
2752{
2753  ideal temp,temp1;
2754  int i,j,k,rk,flength=0,slength,length;
2755  intvec * w;
2756  poly p,q;
2757
2758  if (idIs0(h2))
2759    return idFreeModule(max(1,h2->ncols));
2760  if (!idIs0(h1))
2761    flength = idRankFreeModule(h1);
2762  slength = idRankFreeModule(h2);
2763  length  = max(flength,slength);
2764  if (length==0)
2765  {
2766    length = 1;
2767  }
2768  temp = idInit(IDELEMS(h2),1);
2769  for (i=0;i<IDELEMS(h2);i++)
2770  {
2771    temp->m[i] = pCopy(h2->m[i]);
2772    q = pOne();
2773    pSetComp(q,i+1+length);
2774    if(temp->m[i]!=NULL)
2775    {
2776      if (slength==0) pShift(&(temp->m[i]),1);
2777      p = temp->m[i];
2778      while (pNext(p)!=NULL) pIter(p);
2779      pNext(p) = q;
2780    }
2781    else
2782      temp->m[i]=q;
2783  }
2784  rk = k = IDELEMS(h2);
2785  if (!idIs0(h1))
2786  {
2787    pEnlargeSet(&(temp->m),IDELEMS(temp),IDELEMS(h1));
2788    IDELEMS(temp) += IDELEMS(h1);
2789    for (i=0;i<IDELEMS(h1);i++)
2790    {
2791      if (h1->m[i]!=NULL)
2792      {
2793        temp->m[k] = pCopy(h1->m[i]);
2794        if (flength==0) pShift(&(temp->m[k]),1);
2795        k++;
2796      }
2797    }
2798  }
2799  pSetSyzComp(length);
2800  temp1 = std(temp,currQuotient,testHomog,&w,NULL,length);
2801  pSetSyzComp(0);
2802  idDelete(&temp);
2803  if (w!=NULL) delete w;
2804  for (i=0;i<IDELEMS(temp1);i++)
2805  {
2806    if ((temp1->m[i]!=NULL)
2807    && (pGetComp(temp1->m[i])<=length))
2808    {
2809      pDelete(&(temp1->m[i]));
2810    }
2811    else
2812    {
2813      pShift(&(temp1->m[i]),-length);
2814    }
2815  }
2816  idSkipZeroes(temp1);
2817  temp1->rank = rk;
2818  return temp1;
2819}
2820
2821int idElem(ideal F)
2822{
2823  int i=0,j=0;
2824
2825  while(j<IDELEMS(F))
2826  {
2827   if ((F->m)[j]!=NULL) i++;
2828   j++;
2829  }
2830  return i;
2831}
2832
2833/*
2834*computes module-weights for liftings of homogeneous modules
2835*/
2836intvec * idMWLift(ideal mod,intvec * weights)
2837{
2838  if (idIs0(mod)) return new intvec(2);
2839  int i=IDELEMS(mod);
2840  while ((i>0) && (mod->m[i-1]==NULL)) i--;
2841  intvec *result = new intvec(i+1);
2842  while (i>0)
2843  {
2844    (*result)[i]=pFDeg(mod->m[i])+(*weights)[pGetComp(mod->m[i])];
2845  }
2846  return result;
2847}
2848
2849/*2
2850*sorts the kbase for idCoef* in a special way (lexicographically
2851*with x_max,...,x_1)
2852*/
2853ideal idCreateSpecialKbase(ideal kBase,intvec ** convert)
2854{
2855  int i;
2856  ideal result;
2857
2858  if (idIs0(kBase)) return NULL;
2859  result = idInit(IDELEMS(kBase),kBase->rank);
2860  *convert = idSort(kBase,FALSE);
2861  for (i=0;i<(*convert)->length();i++)
2862  {
2863    result->m[i] = pCopy(kBase->m[(**convert)[i]-1]);
2864  }
2865  return result;
2866}
2867
2868/*2
2869*returns the index of a given monom in the list of the special kbase
2870*/
2871int idIndexOfKBase(poly monom, ideal kbase)
2872{
2873  int j=IDELEMS(kbase);
2874
2875  while ((j>0) && (kbase->m[j-1]==NULL)) j--;
2876  if (j==0) return -1;
2877  int i=pVariables;
2878  while (i>=0)
2879  {
2880    while ((j>0) && (pGetExp(monom,i)<pGetExp(kbase->m[j-1],i))) j--;
2881    if ((j==0) || (pGetExp(monom,i)>pGetExp(kbase->m[j-1],i))) return -1;
2882    if (i==0)
2883      return j-1;
2884    else
2885      i--;
2886  }
2887  return -1;
2888}
2889
2890/*2
2891*decomposes the monom in a part of coefficients described by the
2892*complement of how and a monom in variables occuring in how, the
2893*index of which in kbase is returned as integer pos (-1 if it don't
2894*exists)
2895*/
2896poly idDecompose(poly monom, poly how, ideal kbase, int * pos)
2897{
2898  int i;
2899  poly coeff=pOne(), base=pOne();
2900
2901  for (i=1;i<=pVariables;i++)
2902  {
2903    if (pGetExp(how,i)>0)
2904    {
2905      pSetExp(base,i,pGetExp(monom,i));
2906    }
2907    else
2908    {
2909      pSetExp(coeff,i,pGetExp(monom,i));
2910    }
2911  }
2912  pSetComp(base,pGetComp(monom));
2913  pSetm(base);
2914  pSetCoeff(coeff,nCopy(pGetCoeff(monom)));
2915  pSetm(coeff);
2916  *pos = idIndexOfKBase(base,kbase);
2917  if (*pos<0)
2918    pDelete(&coeff);
2919  pDelete(&base);
2920  return coeff;
2921}
2922
2923/*2
2924*returns a matrix A of coefficients with kbase*A=arg
2925*if all monomials in variables of how occur in kbase
2926*the other are deleted
2927*/
2928matrix idCoeffOfKBase(ideal arg, ideal kbase, poly how)
2929{
2930  matrix result;
2931  ideal tempKbase;
2932  poly p,q;
2933  intvec * convert;
2934  int i=IDELEMS(kbase),j=IDELEMS(arg),k,pos;
2935
2936  while ((i>0) && (kbase->m[i-1]==NULL)) i--;
2937  if (idIs0(arg))
2938    return mpNew(i,1);
2939  while ((j>0) && (arg->m[j-1]==NULL)) j--;
2940  result = mpNew(i,j);
2941  tempKbase = idCreateSpecialKbase(kbase,&convert);
2942  for (k=0;k<j;k++)
2943  {
2944    p = arg->m[k];
2945    while (p!=NULL)
2946    {
2947      q = idDecompose(p,how,tempKbase,&pos);
2948      if (pos>=0)
2949      {
2950        MATELEM(result,(*convert)[pos],k+1) =
2951            pAdd(MATELEM(result,(*convert)[pos],k+1),q);
2952      }
2953      else
2954        pDelete(&q);
2955      p = pIter(p);
2956    }
2957  }
2958  idDelete(&tempKbase);
2959  return result;
2960}
2961
2962intvec * idQHomWeights(ideal id)
2963{
2964  intvec * imat=new intvec(2*pVariables,pVariables,0);
2965  poly actHead=NULL,wPoint=NULL;
2966  int actIndex,i=-1,j=1,k;
2967  BOOLEAN notReady=TRUE;
2968
2969  while (notReady)
2970  {
2971    if (wPoint==NULL)
2972    {
2973      i++;
2974      while ((i<IDELEMS(id))
2975      && ((id->m[i]==NULL) || (pNext(id->m[i])==NULL)))
2976        i++;
2977      if (i<IDELEMS(id))
2978      {
2979        actHead = id->m[i];
2980        wPoint = pNext(actHead);
2981      }
2982    }
2983    while ((wPoint!=NULL) && (j<=2*pVariables))
2984    {
2985      for (k=1;k<=pVariables;k++)
2986        IMATELEM(*imat,j,k) += pGetExp(actHead,k)-pGetExp(wPoint,k);
2987      pIter(wPoint);
2988      j++;
2989    }
2990    if ((i>=IDELEMS(id)) || (j>2*pVariables))
2991    {
2992      ivTriangMat(imat,1,1);
2993      j = ivFirstEmptyRow(imat);
2994      if ((i>=IDELEMS(id)) || (j>pVariables)) notReady=FALSE;
2995    }
2996  }
2997  intvec *result=NULL;
2998  if (j<=pVariables)
2999  {
3000    result=ivSolveIntMat(imat);
3001  }
3002  //else
3003  //{
3004  //  WerrorS("not homogeneous");
3005  //}
3006  delete imat;
3007  return result;
3008}
3009
3010/*2
3011* returns the presentation of an isomorphic, minimally
3012* embedded  module
3013*/
3014ideal idMinEmbedding(ideal arg)
3015{
3016  if (idIs0(arg)) return idInit(1,arg->rank);
3017
3018  int i,j,k,pC;
3019  poly p,q;
3020  int rg=arg->rank;
3021  ideal res = idCopy(arg);
3022  intvec *indexMap=new intvec(rg+1);
3023  intvec *toKill=new intvec(rg+1);
3024
3025  loop
3026  {
3027    k = 0;
3028    for (i=indexMap->length()-1;i>0;i--)
3029    {
3030      (*indexMap)[i] = i;
3031      (*toKill)[i] = 0;
3032    }
3033    for (j=IDELEMS(res)-1;j>=0;j--)
3034    {
3035      if ((res->m[j]!=NULL) && (pIsConstantComp(res->m[j])) &&
3036           (pNext(res->m[j])==NULL))
3037      {
3038        pC = pGetComp(res->m[j]);
3039        if ((*toKill)[pC]==0)
3040        {
3041          rg--;
3042          (*toKill)[pC] = 1;
3043          for (i=indexMap->length()-1;i>=pC;i--)
3044            (*indexMap)[i]--;
3045        }
3046        pDelete(&(res->m[j]));
3047        k++;
3048      }
3049    }
3050    idSkipZeroes(res);
3051    if (k==0) break;
3052    if (rg>0)
3053    {
3054      res->rank=rg;
3055      for (j=IDELEMS(res)-1;j>=0;j--)
3056      {
3057        while ((res->m[j]!=NULL) && ((*toKill)[pGetComp(res->m[j])]==1))
3058          pDelete1(&res->m[j]);
3059        p = res->m[j];
3060        while ((p!=NULL) && (pNext(p)!=NULL))
3061        {
3062          pSetComp(p,(*indexMap)[pGetComp(p)]);
3063          while ((pNext(p)!=NULL) && ((*toKill)[pGetComp(pNext(p))]==1))
3064            pDelete1(&pNext(p));
3065          pIter(p);
3066        }
3067        if (p!=NULL) pSetComp(p,(*indexMap)[pGetComp(p)]);
3068      }
3069      idSkipZeroes(res);
3070    }
3071    else
3072    {
3073      idDelete(&res);
3074      res=idFreeModule(1);
3075      break;
3076    }
3077  }
3078  delete toKill;
3079  delete indexMap;
3080  return res;
3081}
3082
Note: See TracBrowser for help on using the repository browser.