source: git/kernel/ideals.cc @ 86016d

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