source: git/libpolys/polys/simpleideals.cc @ a665eb

spielwiese
Last change on this file since a665eb was a665eb, checked in by Hans Schoenemann <hannes@…>, 13 years ago
id_HomIdeal, id_MaxIdeal
  • Property mode set to 100644
File size: 54.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT - all basic methods to manipulate ideals
7*/
8
9
10/* includes */
11#include "config.h"
12#include <misc/auxiliary.h>
13#include <misc/options.h>
14#include <omalloc/omalloc.h>
15#include <polys/monomials/p_polys.h>
16#include <misc/intvec.h>
17#include <polys/simpleideals.h>
18
19/*2
20* initialise an ideal
21*/
22ideal idInit(int idsize, int rank)
23{
24  /*- initialise an ideal -*/
25  ideal hh = (ideal )omAllocBin(sip_sideal_bin);
26  hh->nrows = 1;
27  hh->rank = rank;
28  IDELEMS(hh) = idsize;
29  if (idsize>0)
30  {
31    hh->m = (poly *)omAlloc0(idsize*sizeof(poly));
32  }
33  else
34    hh->m=NULL;
35  return hh;
36}
37
38#ifdef PDEBUG
39// this is only for outputting an ideal within the debugger
40void idShow(const ideal id, const ring lmRing, const ring tailRing, const int debugPrint)
41{
42  assume( debugPrint >= 0 );
43
44  if( id == NULL )
45    PrintS("(NULL)");
46  else
47  {
48    Print("Module of rank %ld,real rank %ld and %d generators.\n",
49          id->rank,id_RankFreeModule(id, lmRing, tailRing),IDELEMS(id));
50
51    int j = (id->ncols*id->nrows) - 1;
52    while ((j > 0) && (id->m[j]==NULL)) j--;
53    for (int i = 0; i <= j; i++)
54    {
55      Print("generator %d: ",i); p_DebugPrint(id->m[i], lmRing, tailRing, debugPrint);
56    }
57  }
58}
59#endif
60
61/* index of generator with leading term in ground ring (if any);
62   otherwise -1 */
63int id_PosConstant(ideal id, const ring r)
64{
65  int k;
66  for (k = IDELEMS(id)-1; k>=0; k--)
67  {
68    if (p_LmIsConstantComp(id->m[k], r) == TRUE)
69      return k;
70  }
71  return -1;
72}
73
74/*2
75* initialise the maximal ideal (at 0)
76*/
77ideal id_MaxIdeal (const ring r)
78{
79  int l;
80  ideal hh=NULL;
81
82  hh=idInit(rVar(r),1);
83  for (l=0; l<rVar(r); l++)
84  {
85    hh->m[l] = p_One(r);
86    p_SetExp(hh->m[l],l+1,1,r);
87    p_Setm(hh->m[l],r);
88  }
89  return hh;
90}
91
92/*2
93* deletes an ideal/matrix
94*/
95void id_Delete (ideal * h, ring r)
96{
97  int j,elems;
98  if (*h == NULL)
99    return;
100  elems=j=(*h)->nrows*(*h)->ncols;
101  if (j>0)
102  {
103    do
104    {
105      p_Delete(&((*h)->m[--j]), r);
106    }
107    while (j>0);
108    omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*elems);
109  }
110  omFreeBin((ADDRESS)*h, sip_sideal_bin);
111  *h=NULL;
112}
113
114
115/*2
116* Shallowdeletes an ideal/matrix
117*/
118void id_ShallowDelete (ideal *h, ring r)
119{
120  int j,elems;
121  if (*h == NULL)
122    return;
123  elems=j=(*h)->nrows*(*h)->ncols;
124  if (j>0)
125  {
126    do
127    {
128      p_ShallowDelete(&((*h)->m[--j]), r);
129    }
130    while (j>0);
131    omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*elems);
132  }
133  omFreeBin((ADDRESS)*h, sip_sideal_bin);
134  *h=NULL;
135}
136
137/*2
138*gives an ideal the minimal possible size
139*/
140void idSkipZeroes (ideal ide)
141{
142  int k;
143  int j = -1;
144  BOOLEAN change=FALSE;
145  for (k=0; k<IDELEMS(ide); k++)
146  {
147    if (ide->m[k] != NULL)
148    {
149      j++;
150      if (change)
151      {
152        ide->m[j] = ide->m[k];
153      }
154    }
155    else
156    {
157      change=TRUE;
158    }
159  }
160  if (change)
161  {
162    if (j == -1)
163      j = 0;
164    else
165    {
166      for (k=j+1; k<IDELEMS(ide); k++)
167        ide->m[k] = NULL;
168    }
169    pEnlargeSet(&(ide->m),IDELEMS(ide),j+1-IDELEMS(ide));
170    IDELEMS(ide) = j+1;
171  }
172}
173
174/*2
175* copies the first k (>= 1) entries of the given ideal
176* and returns these as a new ideal
177* (Note that the copied polynomials may be zero.)
178*/
179ideal id_CopyFirstK (const ideal ide, const int k,const ring r)
180{
181  ideal newI = idInit(k, 0);
182  for (int i = 0; i < k; i++)
183    newI->m[i] = p_Copy(ide->m[i],r);
184  return newI;
185}
186
187/*2
188* ideal id = (id[i])
189* result is leadcoeff(id[i]) = 1
190*/
191void id_Norm(ideal id, const ring r)
192{
193  for (int i=IDELEMS(id)-1; i>=0; i--)
194  {
195    if (id->m[i] != NULL)
196    {
197      p_Norm(id->m[i],r);
198    }
199  }
200}
201
202/*2
203* ideal id = (id[i]), c any unit
204* if id[i] = c*id[j] then id[j] is deleted for j > i
205*/
206void id_DelMultiples(ideal id, const ring r)
207{
208  int i, j;
209  int k = IDELEMS(id)-1;
210  for (i=k; i>=0; i--)
211  {
212    if (id->m[i]!=NULL)
213    {
214      for (j=k; j>i; j--)
215      {
216        if (id->m[j]!=NULL)
217        {
218#ifdef HAVE_RINGS
219          if (rField_is_Ring(r))
220          {
221            /* if id[j] = c*id[i] then delete id[j].
222               In the below cases of a ground field, we
223               check whether id[i] = c*id[j] and, if so,
224               delete id[j] for historical reasons (so
225               that previous output does not change) */
226            if (p_ComparePolys(id->m[j], id->m[i],r)) p_Delete(&id->m[j],r);
227          }
228          else
229          {
230            if (p_ComparePolys(id->m[i], id->m[j],r)) p_Delete(&id->m[j],r);
231          }
232#else
233          if (p_ComparePolys(id->m[i], id->m[j],r)) p_Delete(&id->m[j],r);
234#endif   
235        }
236      }
237    }
238  }
239}
240
241/*2
242* ideal id = (id[i])
243* if id[i] = id[j] then id[j] is deleted for j > i
244*/
245void id_DelEquals(ideal id, const ring r)
246{
247  int i, j;
248  int k = IDELEMS(id)-1;
249  for (i=k; i>=0; i--)
250  {
251    if (id->m[i]!=NULL)
252    {
253      for (j=k; j>i; j--)
254      {
255        if ((id->m[j]!=NULL)
256        && (p_EqualPolys(id->m[i], id->m[j],r)))
257        {
258          p_Delete(&id->m[j],r);
259        }
260      }
261    }
262  }
263}
264
265//
266// Delete id[j], if Lm(j) == Lm(i) and both LC(j), LC(i) are units and j > i
267//
268void id_DelLmEquals(ideal id, const ring r)
269{
270  int i, j;
271  int k = IDELEMS(id)-1;
272  for (i=k; i>=0; i--)
273  {
274    if (id->m[i] != NULL)
275    {
276      for (j=k; j>i; j--)
277      {
278        if ((id->m[j] != NULL)
279        && p_LmEqual(id->m[i], id->m[j],r)
280#ifdef HAVE_RINGS
281        && n_IsUnit(pGetCoeff(id->m[i]),r->cf) && n_IsUnit(pGetCoeff(id->m[j]),r->cf)
282#endif
283        )
284        {
285          p_Delete(&id->m[j],r);
286        }
287      }
288    }
289  }
290}
291
292//
293// delete id[j], if LT(j) == coeff*mon*LT(i) and vice versa, i.e.,
294// delete id[i], if LT(i) == coeff*mon*LT(j)
295//
296void id_DelDiv(ideal id, const ring r)
297{
298  int i, j;
299  int k = IDELEMS(id)-1;
300  for (i=k; i>=0; i--)
301  {
302    if (id->m[i] != NULL)
303    {
304      for (j=k; j>i; j--)
305      {
306        if (id->m[j]!=NULL)
307        {
308#ifdef HAVE_RINGS
309          if (rField_is_Ring(r))
310          {
311            if (p_DivisibleByRingCase(id->m[i], id->m[j],r))
312            {
313              p_Delete(&id->m[j],r);
314            }
315            else if (p_DivisibleByRingCase(id->m[j], id->m[i],r))
316            {
317              p_Delete(&id->m[i],r);
318              break;
319            }
320          }
321          else
322          {
323#endif
324          /* the case of a ground field: */
325          if (p_DivisibleBy(id->m[i], id->m[j],r))
326          {
327            p_Delete(&id->m[j],r);
328          }
329          else if (p_DivisibleBy(id->m[j], id->m[i],r))
330          {
331            p_Delete(&id->m[i],r);
332            break;
333          }
334#ifdef HAVE_RINGS
335          }
336#endif   
337        }
338      }
339    }
340  }
341}
342
343/*2
344*test if the ideal has only constant polynomials
345*/
346BOOLEAN id_IsConstant(ideal id, const ring r)
347{
348  int k;
349  for (k = IDELEMS(id)-1; k>=0; k--)
350  {
351    if (!p_IsConstantPoly(id->m[k],r))
352      return FALSE;
353  }
354  return TRUE;
355}
356
357/*2
358* copy an ideal
359*/
360#ifdef PDEBUG
361ideal id_DBCopy(ideal h1,const char *f,int l, const ring r)
362{
363  int i;
364  ideal h2;
365
366  id_DBTest(h1,PDEBUG,f,l,r);
367//#ifdef TEST
368  if (h1 == NULL)
369  {
370    h2=idDBInit(1,1,f,l);
371  }
372  else
373//#endif
374  {
375    h2=idDBInit(IDELEMS(h1),h1->rank,f,l);
376    for (i=IDELEMS(h1)-1; i>=0; i--)
377      h2->m[i] = p_Copy(h1->m[i],r);
378  }
379  return h2;
380}
381#else
382ideal id_Copy(ideal h1, const ring r)
383{
384  int i;
385  ideal h2;
386
387//#ifdef TEST
388  if (h1 == NULL)
389  {
390    h2=idInit(1,1);
391  }
392  else
393//#endif
394  {
395    h2=idInit(IDELEMS(h1),h1->rank);
396    for (i=IDELEMS(h1)-1; i>=0; i--)
397      h2->m[i] = p_Copy(h1->m[i],r);
398  }
399  return h2;
400}
401#endif
402
403#ifdef PDEBUG
404void id_DBTest(ideal h1, int level, const char *f,const int l, const ring r)
405{
406  int i;
407
408  if (h1 != NULL)
409  {
410    // assume(IDELEMS(h1) > 0); for ideal/module, does not apply to matrix
411    omCheckAddrSize(h1,sizeof(*h1));
412    omdebugAddrSize(h1->m,h1->ncols*h1->nrows*sizeof(poly));
413    /* to be able to test matrices: */
414    for (i=(h1->ncols*h1->nrows)-1; i>=0; i--)
415      _p_Test(h1->m[i], r, level);
416    int new_rk=id_RankFreeModule(h1,r);
417    if(new_rk > h1->rank)
418    {
419      dReportError("wrong rank %d (should be %d) in %s:%d\n",
420                   h1->rank, new_rk, f,l);
421      omPrintAddrInfo(stderr, h1, " for ideal");
422      h1->rank=new_rk;
423    }
424  }
425}
426#endif
427
428/*3
429* for idSort: compare a and b revlex inclusive module comp.
430*/
431static int p_Comp_RevLex(poly a, poly b,BOOLEAN nolex, const ring R)
432{
433  if (b==NULL) return 1;
434  if (a==NULL) return -1;
435
436  if (nolex)
437  {
438    int r=p_LmCmp(a,b,R);
439    if (r!=0) return r;
440    number h=n_Sub(pGetCoeff(a),pGetCoeff(b),R->cf);
441    r = -1+n_IsZero(h,R->cf)+2*n_GreaterZero(h,R->cf); /* -1: <, 0:==, 1: > */
442    n_Delete(&h, R->cf);
443    return r;
444  }
445  int l=rVar(R);
446  while ((l>0) && (p_GetExp(a,l,R)==p_GetExp(b,l,R))) l--;
447  if (l==0)
448  {
449    if (p_GetComp(a,R)==p_GetComp(b,R))
450    {
451      number h=n_Sub(pGetCoeff(a),pGetCoeff(b),R->cf);
452      int r = -1+n_IsZero(h,R->cf)+2*n_GreaterZero(h,R->cf); /* -1: <, 0:==, 1: > */
453      n_Delete(&h,R->cf);
454      return r;
455    }
456    if (p_GetComp(a,R)>p_GetComp(b,R)) return 1;
457  }
458  else if (p_GetExp(a,l,R)>p_GetExp(b,l,R))
459    return 1;
460  return -1;
461}
462
463/*2
464*sorts the ideal w.r.t. the actual ringordering
465*uses lex-ordering when nolex = FALSE
466*/
467intvec *id_Sort(ideal id,BOOLEAN nolex, const ring r)
468{
469  poly p,q;
470  intvec * result = new intvec(IDELEMS(id));
471  int i, j, actpos=0, newpos, l;
472  int diff, olddiff, lastcomp, newcomp;
473  BOOLEAN notFound;
474
475  for (i=0;i<IDELEMS(id);i++)
476  {
477    if (id->m[i]!=NULL)
478    {
479      notFound = TRUE;
480      newpos = actpos / 2;
481      diff = (actpos+1) / 2;
482      diff = (diff+1) / 2;
483      lastcomp = p_Comp_RevLex(id->m[i],id->m[(*result)[newpos]],nolex,r);
484      if (lastcomp<0)
485      {
486        newpos -= diff;
487      }
488      else if (lastcomp>0)
489      {
490        newpos += diff;
491      }
492      else
493      {
494        notFound = FALSE;
495      }
496      //while ((newpos>=0) && (newpos<actpos) && (notFound))
497      while (notFound && (newpos>=0) && (newpos<actpos))
498      {
499        newcomp = p_Comp_RevLex(id->m[i],id->m[(*result)[newpos]],nolex,r);
500        olddiff = diff;
501        if (diff>1)
502        {
503          diff = (diff+1) / 2;
504          if ((newcomp==1)
505          && (actpos-newpos>1)
506          && (diff>1)
507          && (newpos+diff>=actpos))
508          {
509            diff = actpos-newpos-1;
510          }
511          else if ((newcomp==-1)
512          && (diff>1)
513          && (newpos<diff))
514          {
515            diff = newpos;
516          }
517        }
518        if (newcomp<0)
519        {
520          if ((olddiff==1) && (lastcomp>0))
521            notFound = FALSE;
522          else
523            newpos -= diff;
524        }
525        else if (newcomp>0)
526        {
527          if ((olddiff==1) && (lastcomp<0))
528          {
529            notFound = FALSE;
530            newpos++;
531          }
532          else
533          {
534            newpos += diff;
535          }
536        }
537        else
538        {
539          notFound = FALSE;
540        }
541        lastcomp = newcomp;
542        if (diff==0) notFound=FALSE; /*hs*/
543      }
544      if (newpos<0) newpos = 0;
545      if (newpos>actpos) newpos = actpos;
546      while ((newpos<actpos) && (p_Comp_RevLex(id->m[i],id->m[(*result)[newpos]],nolex,r)==0))
547        newpos++;
548      for (j=actpos;j>newpos;j--)
549      {
550        (*result)[j] = (*result)[j-1];
551      }
552      (*result)[newpos] = i;
553      actpos++;
554    }
555  }
556  for (j=0;j<actpos;j++) (*result)[j]++;
557  return result;
558}
559
560/*2
561* concat the lists h1 and h2 without zeros
562*/
563ideal id_SimpleAdd (ideal h1,ideal h2, const ring R)
564{
565  int i,j,r,l;
566  ideal result;
567
568  if (h1==NULL) return id_Copy(h2,R);
569  if (h2==NULL) return id_Copy(h1,R);
570  j = IDELEMS(h1)-1;
571  while ((j >= 0) && (h1->m[j] == NULL)) j--;
572  i = IDELEMS(h2)-1;
573  while ((i >= 0) && (h2->m[i] == NULL)) i--;
574  r = si_max(h1->rank,h2->rank);
575  if (i+j==(-2))
576    return idInit(1,r);
577  else
578    result=idInit(i+j+2,r);
579  for (l=j; l>=0; l--)
580  {
581    result->m[l] = p_Copy(h1->m[l],R);
582  }
583  r = i+j+1;
584  for (l=i; l>=0; l--, r--)
585  {
586    result->m[r] = p_Copy(h2->m[l],R);
587  }
588  return result;
589}
590
591/*2
592* insert h2 into h1 (if h2 is not the zero polynomial)
593* return TRUE iff h2 was indeed inserted
594*/
595BOOLEAN idInsertPoly (ideal h1, poly h2)
596{
597  if (h2==NULL) return FALSE;
598  int j = IDELEMS(h1)-1;
599  while ((j >= 0) && (h1->m[j] == NULL)) j--;
600  j++;
601  if (j==IDELEMS(h1))
602  {
603    pEnlargeSet(&(h1->m),IDELEMS(h1),16);
604    IDELEMS(h1)+=16;
605  }
606  h1->m[j]=h2;
607  return TRUE;
608}
609
610/*2
611* insert h2 into h1 depending on the two boolean parameters:
612* - if zeroOk is true, then h2 will also be inserted when it is zero
613* - if duplicateOk is true, then h2 will also be inserted when it is
614*   already present in h1
615* return TRUE iff h2 was indeed inserted
616*/
617BOOLEAN id_InsertPolyWithTests (ideal h1, const int validEntries,
618  const poly h2, const bool zeroOk, const bool duplicateOk, const ring r)
619{
620  if ((!zeroOk) && (h2 == NULL)) return FALSE;
621  if (!duplicateOk)
622  {
623    bool h2FoundInH1 = false;
624    int i = 0;
625    while ((i < validEntries) && (!h2FoundInH1))
626    {
627      h2FoundInH1 = p_EqualPolys(h1->m[i], h2,r);
628      i++;
629    }
630    if (h2FoundInH1) return FALSE;
631  }
632  if (validEntries == IDELEMS(h1))
633  {
634    pEnlargeSet(&(h1->m), IDELEMS(h1), 16);
635    IDELEMS(h1) += 16;
636  }
637  h1->m[validEntries] = h2;
638  return TRUE;
639}
640
641/*2
642* h1 + h2
643*/
644ideal id_Add (ideal h1,ideal h2, const ring r)
645{
646  ideal result = id_SimpleAdd(h1,h2,r);
647  id_Compactify(result,r);
648  return result;
649}
650
651/*2
652* h1 * h2
653*/
654ideal  id_Mult (ideal h1,ideal  h2, const ring r)
655{
656  int i,j,k;
657  ideal  hh;
658
659  j = IDELEMS(h1);
660  while ((j > 0) && (h1->m[j-1] == NULL)) j--;
661  i = IDELEMS(h2);
662  while ((i > 0) && (h2->m[i-1] == NULL)) i--;
663  j = j * i;
664  if (j == 0)
665    hh = idInit(1,1);
666  else
667    hh=idInit(j,1);
668  if (h1->rank<h2->rank)
669    hh->rank = h2->rank;
670  else
671    hh->rank = h1->rank;
672  if (j==0) return hh;
673  k = 0;
674  for (i=0; i<IDELEMS(h1); i++)
675  {
676    if (h1->m[i] != NULL)
677    {
678      for (j=0; j<IDELEMS(h2); j++)
679      {
680        if (h2->m[j] != NULL)
681        {
682          hh->m[k] = pp_Mult_qq(h1->m[i],h2->m[j],r);
683          k++;
684        }
685      }
686    }
687  }
688  {
689    id_Compactify(hh,r);
690    return hh;
691  }
692}
693
694/*2
695*returns true if h is the zero ideal
696*/
697BOOLEAN idIs0 (ideal h)
698{
699  int i;
700
701  if (h == NULL) return TRUE;
702  i = IDELEMS(h)-1;
703  while ((i >= 0) && (h->m[i] == NULL))
704  {
705    i--;
706  }
707  if (i < 0)
708    return TRUE;
709  else
710    return FALSE;
711}
712
713/*2
714* return the maximal component number found in any polynomial in s
715*/
716long idRankFreeModule (ideal s, ring lmRing, ring tailRing)
717{
718  if (s!=NULL)
719  {
720    int  j=0;
721
722    if (rRing_has_Comp(tailRing) && rRing_has_Comp(lmRing))
723    {
724      int  l=IDELEMS(s);
725      poly *p=s->m;
726      int  k;
727      for (; l != 0; l--)
728      {
729        if (*p!=NULL)
730        {
731          pp_Test(*p, lmRing, tailRing);
732          k = p_MaxComp(*p, lmRing, tailRing);
733          if (k>j) j = k;
734        }
735        p++;
736      }
737    }
738    return j;
739  }
740  return -1;
741}
742
743BOOLEAN idIsModule(ideal id, ring r)
744{
745  if (id != NULL && rRing_has_Comp(r))
746  {
747    int j, l = IDELEMS(id);
748    for (j=0; j<l; j++)
749    {
750      if (id->m[j] != NULL && p_GetComp(id->m[j], r) > 0) return TRUE;
751    }
752  }
753  return FALSE;
754}
755
756
757/*2
758*returns true if id is homogenous with respect to the aktual weights
759*/
760BOOLEAN id_HomIdeal (ideal id, ideal Q, const ring r)
761{
762  int i;
763  BOOLEAN     b;
764  if ((id == NULL) || (IDELEMS(id) == 0)) return TRUE;
765  i = 0;
766  b = TRUE;
767  while ((i < IDELEMS(id)) && b)
768  {
769    b = p_IsHomogeneous(id->m[i],r);
770    i++;
771  }
772  if ((b) && (Q!=NULL) && (IDELEMS(Q)>0))
773  {
774    i=0;
775    while ((i < IDELEMS(Q)) && b)
776    {
777      b = p_IsHomogeneous(Q->m[i],r);
778      i++;
779    }
780  }
781  return b;
782}
783
784/*2
785*the minimal index of used variables - 1
786*/
787int p_LowVar (poly p, const ring r)
788{
789  int k,l,lex;
790
791  if (p == NULL) return -1;
792
793  k = 32000;/*a very large dummy value*/
794  while (p != NULL)
795  {
796    l = 1;
797    lex = p_GetExp(p,l,r);
798    while ((l < rVar(r)) && (lex == 0))
799    {
800      l++;
801      lex = p_GetExp(p,l,r);
802    }
803    l--;
804    if (l < k) k = l;
805    pIter(p);
806  }
807  return k;
808}
809
810/*3
811*multiplies p with t (!cas) or  (t-1)
812*the index of t is:1, so we have to shift all variables
813*p is NOT in the actual ring, it has no t
814*/
815static poly p_MultWithT (poly p,BOOLEAN cas, const ring r)
816{
817  /*qp is the working pointer in p*/
818  /*result is the result, qresult is the working pointer*/
819  /*pp is p in the actual ring(shifted), qpp the working pointer*/
820  poly result,qp,pp;
821  poly qresult=NULL;
822  poly qpp=NULL;
823  int  i,j,lex;
824  number n;
825
826  pp = NULL;
827  result = NULL;
828  qp = p;
829  while (qp != NULL)
830  {
831    i = 0;
832    if (result == NULL)
833    {/*first monomial*/
834      result = p_Init(r);
835      qresult = result;
836    }
837    else
838    {
839      qresult->next = p_Init(r);
840      pIter(qresult);
841    }
842    for (j=rVar(r)-1; j>0; j--)
843    {
844      lex = p_GetExp(qp,j,r);
845      p_SetExp(qresult,j+1,lex,r);/*copy all variables*/
846    }
847    lex = p_GetComp(qp,r);
848    p_SetComp(qresult,lex,r);
849    n=n_Copy(pGetCoeff(qp),r->cf);
850    pSetCoeff0(qresult,n);
851    qresult->next = NULL;
852    p_Setm(qresult,r);
853    /*qresult is now qp brought into the actual ring*/
854    if (cas)
855    { /*case: mult with t-1*/
856      p_SetExp(qresult,1,0,r);
857      p_Setm(qresult,r);
858      if (pp == NULL)
859      { /*first monomial*/
860        pp = p_Copy(qresult,r);
861        qpp = pp;
862      }
863      else
864      {
865        qpp->next = p_Copy(qresult,r);
866        pIter(qpp);
867      }
868      pGetCoeff(qpp)=n_Neg(pGetCoeff(qpp),r->cf);
869      /*now qpp contains -1*qp*/
870    }
871    p_SetExp(qresult,1,1,r);/*this is mult. by t*/
872    p_Setm(qresult,r);
873    pIter(qp);
874  }
875  /*
876  *now p is processed:
877  *result contains t*p
878  * if cas: pp contains -1*p (in the new ring)
879  */
880  if (cas)  qresult->next = pp;
881  /*  else      qresult->next = NULL;*/
882  return result;
883}
884
885/*2
886* verschiebt die Indizees der Modulerzeugenden um i
887*/
888void pShift (poly * p,int i)
889{
890  poly qp1 = *p,qp2 = *p;/*working pointers*/
891  int     j = pMaxComp(*p),k = pMinComp(*p);
892
893  if (j+i < 0) return ;
894  while (qp1 != NULL)
895  {
896    if ((pGetComp(qp1)+i > 0) || ((j == -i) && (j == k)))
897    {
898      pAddComp(qp1,i);
899      pSetmComp(qp1);
900      qp2 = qp1;
901      pIter(qp1);
902    }
903    else
904    {
905      if (qp2 == *p)
906      {
907        pIter(*p);
908        pLmDelete(&qp2);
909        qp2 = *p;
910        qp1 = *p;
911      }
912      else
913      {
914        qp2->next = qp1->next;
915        if (qp1!=NULL) pLmDelete(&qp1);
916        qp1 = qp2->next;
917      }
918    }
919  }
920}
921
922/*2
923*initialized a field with r numbers between beg and end for the
924*procedure idNextChoise
925*/
926void idInitChoise (int r,int beg,int end,BOOLEAN  *endch,int * choise)
927{
928  /*returns the first choise of r numbers between beg and end*/
929  int i;
930  for (i=0; i<r; i++)
931  {
932    choise[i] = 0;
933  }
934  if (r <= end-beg+1)
935    for (i=0; i<r; i++)
936    {
937      choise[i] = beg+i;
938    }
939  if (r > end-beg+1)
940    *endch = TRUE;
941  else
942    *endch = FALSE;
943}
944
945/*2
946*returns the next choise of r numbers between beg and end
947*/
948void idGetNextChoise (int r,int end,BOOLEAN  *endch,int * choise)
949{
950  int i = r-1,j;
951  while ((i >= 0) && (choise[i] == end))
952  {
953    i--;
954    end--;
955  }
956  if (i == -1)
957    *endch = TRUE;
958  else
959  {
960    choise[i]++;
961    for (j=i+1; j<r; j++)
962    {
963      choise[j] = choise[i]+j-i;
964    }
965    *endch = FALSE;
966  }
967}
968
969/*2
970*takes the field choise of d numbers between beg and end, cancels the t-th
971*entree and searches for the ordinal number of that d-1 dimensional field
972* w.r.t. the algorithm of construction
973*/
974int idGetNumberOfChoise(int t, int d, int begin, int end, int * choise)
975{
976  int * localchoise,i,result=0;
977  BOOLEAN b=FALSE;
978
979  if (d<=1) return 1;
980  localchoise=(int*)omAlloc((d-1)*sizeof(int));
981  idInitChoise(d-1,begin,end,&b,localchoise);
982  while (!b)
983  {
984    result++;
985    i = 0;
986    while ((i<t) && (localchoise[i]==choise[i])) i++;
987    if (i>=t)
988    {
989      i = t+1;
990      while ((i<d) && (localchoise[i-1]==choise[i])) i++;
991      if (i>=d)
992      {
993        omFreeSize((ADDRESS)localchoise,(d-1)*sizeof(int));
994        return result;
995      }
996    }
997    idGetNextChoise(d-1,end,&b,localchoise);
998  }
999  omFreeSize((ADDRESS)localchoise,(d-1)*sizeof(int));
1000  return 0;
1001}
1002
1003/*2
1004*computes the binomial coefficient
1005*/
1006int binom (int n,int r)
1007{
1008  int i,result;
1009
1010  if (r==0) return 1;
1011  if (n-r<r) return binom(n,n-r);
1012  result = n-r+1;
1013  for (i=2;i<=r;i++)
1014  {
1015    result *= n-r+i;
1016    if (result<0)
1017    {
1018      WarnS("overflow in binomials");
1019      return 0;
1020    }
1021    result /= i;
1022  }
1023  return result;
1024}
1025
1026/*2
1027*the free module of rank i
1028*/
1029ideal idFreeModule (int i)
1030{
1031  int j;
1032  ideal h;
1033
1034  h=idInit(i,i);
1035  for (j=0; j<i; j++)
1036  {
1037    h->m[j] = p_One(r);
1038    pSetComp(h->m[j],j+1);
1039    pSetmComp(h->m[j]);
1040  }
1041  return h;
1042}
1043
1044ideal idSectWithElim (ideal h1,ideal h2)
1045// does not destroy h1,h2
1046{
1047  if (TEST_OPT_PROT) PrintS("intersect by elimination method\n");
1048  assume(!idIs0(h1));
1049  assume(!idIs0(h2));
1050  assume(IDELEMS(h1)<=IDELEMS(h2));
1051  assume(idRankFreeModule(h1)==0);
1052  assume(idRankFreeModule(h2)==0);
1053  // add a new variable:
1054  int j;
1055  ring origRing=currRing;
1056  ring r=rCopy0(origRing);
1057  r->N++;
1058  r->block0[0]=1;
1059  r->block1[0]= r->N;
1060  omFree(r->order);
1061  r->order=(int*)omAlloc0(3*sizeof(int*));
1062  r->order[0]=ringorder_dp;
1063  r->order[1]=ringorder_C;
1064  char **names=(char**)omAlloc0(rVar(r) * sizeof(char_ptr));
1065  for (j=0;j<r->N-1;j++) names[j]=r->names[j];
1066  names[r->N-1]=omStrDup("@");
1067  omFree(r->names);
1068  r->names=names;
1069  rComplete(r,TRUE);
1070  // fetch h1, h2
1071  ideal h;
1072  h1=idrCopyR(h1,origRing,r);
1073  h2=idrCopyR(h2,origRing,r);
1074  // switch to temp. ring r
1075  rChangeCurrRing(r);
1076  // create 1-t, t
1077  poly omt=p_One(r);
1078  pSetExp(omt,r->N,1);
1079  poly t=pCopy(omt);
1080  pSetm(omt);
1081  omt=pNeg(omt);
1082  omt=pAdd(omt,p_One(r));
1083  // compute (1-t)*h1
1084  h1=(ideal)mpMultP((matrix)h1,omt);
1085  // compute t*h2
1086  h2=(ideal)mpMultP((matrix)h2,pCopy(t));
1087  // (1-t)h1 + t*h2
1088  h=idInit(IDELEMS(h1)+IDELEMS(h2),1);
1089  int l;
1090  for (l=IDELEMS(h1)-1; l>=0; l--)
1091  {
1092    h->m[l] = h1->m[l];  h1->m[l]=NULL;
1093  }
1094  j=IDELEMS(h1);
1095  for (l=IDELEMS(h2)-1; l>=0; l--)
1096  {
1097    h->m[l+j] = h2->m[l];  h2->m[l]=NULL;
1098  }
1099  idDelete(&h1);
1100  idDelete(&h2);
1101  // eliminate t:
1102
1103  ideal res=idElimination(h,t);
1104  // cleanup
1105  idDelete(&h);
1106  res=idrMoveR(res,r,origRing);
1107  rChangeCurrRing(origRing);
1108  rKill(r);
1109  return res;
1110}
1111
1112/*2
1113*computes recursively all monomials of a certain degree
1114*in every step the actvar-th entry in the exponential
1115*vector is incremented and the other variables are
1116*computed by recursive calls of makemonoms
1117*if the last variable is reached, the difference to the
1118*degree is computed directly
1119*vars is the number variables
1120*actvar is the actual variable to handle
1121*deg is the degree of the monomials to compute
1122*monomdeg is the actual degree of the monomial in consideration
1123*/
1124static void makemonoms(int vars,int actvar,int deg,int monomdeg)
1125{
1126  poly p;
1127  int i=0;
1128
1129  if ((idpowerpoint == 0) && (actvar ==1))
1130  {
1131    idpower[idpowerpoint] = p_One(r);
1132    monomdeg = 0;
1133  }
1134  while (i<=deg)
1135  {
1136    if (deg == monomdeg)
1137    {
1138      pSetm(idpower[idpowerpoint]);
1139      idpowerpoint++;
1140      return;
1141    }
1142    if (actvar == vars)
1143    {
1144      pSetExp(idpower[idpowerpoint],actvar,deg-monomdeg);
1145      pSetm(idpower[idpowerpoint]);
1146      pTest(idpower[idpowerpoint]);
1147      idpowerpoint++;
1148      return;
1149    }
1150    else
1151    {
1152      p = pCopy(idpower[idpowerpoint]);
1153      makemonoms(vars,actvar+1,deg,monomdeg);
1154      idpower[idpowerpoint] = p;
1155    }
1156    monomdeg++;
1157    pSetExp(idpower[idpowerpoint],actvar,pGetExp(idpower[idpowerpoint],actvar)+1);
1158    pSetm(idpower[idpowerpoint]);
1159    pTest(idpower[idpowerpoint]);
1160    i++;
1161  }
1162}
1163
1164/*2
1165*returns the deg-th power of the maximal ideal of 0
1166*/
1167ideal id_MaxIdeal(int deg, const ring r)
1168{
1169  if (deg < 0)
1170  {
1171    WarnS("maxideal: power must be non-negative");
1172  }
1173  if (deg < 1)
1174  {
1175    ideal I=idInit(1,1);
1176    I->m[0]=p_One(r);
1177    return I;
1178  }
1179  if (deg == 1)
1180  {
1181    return idMaxIdeal(r);
1182  }
1183
1184  int vars = rVar(r);
1185  int i = binom(vars+deg-1,deg);
1186  if (i<=0) return idInit(1,1);
1187  ideal id=idInit(i,1);
1188  idpower = id->m;
1189  idpowerpoint = 0;
1190  makemonoms(vars,1,deg,0);
1191  idpower = NULL;
1192  idpowerpoint = 0;
1193  return id;
1194}
1195
1196/*2
1197*computes recursively all generators of a certain degree
1198*of the ideal "givenideal"
1199*elms is the number elements in the given ideal
1200*actelm is the actual element to handle
1201*deg is the degree of the power to compute
1202*gendeg is the actual degree of the generator in consideration
1203*/
1204static void makepotence(int elms,int actelm,int deg,int gendeg)
1205{
1206  poly p;
1207  int i=0;
1208
1209  if ((idpowerpoint == 0) && (actelm ==1))
1210  {
1211    idpower[idpowerpoint] = p_One(r);
1212    gendeg = 0;
1213  }
1214  while (i<=deg)
1215  {
1216    if (deg == gendeg)
1217    {
1218      idpowerpoint++;
1219      return;
1220    }
1221    if (actelm == elms)
1222    {
1223      p=pPower(pCopy(givenideal[actelm-1]),deg-gendeg);
1224      idpower[idpowerpoint]=pMult(idpower[idpowerpoint],p);
1225      idpowerpoint++;
1226      return;
1227    }
1228    else
1229    {
1230      p = pCopy(idpower[idpowerpoint]);
1231      makepotence(elms,actelm+1,deg,gendeg);
1232      idpower[idpowerpoint] = p;
1233    }
1234    gendeg++;
1235    idpower[idpowerpoint]=pMult(idpower[idpowerpoint],pCopy(givenideal[actelm-1]));
1236    i++;
1237  }
1238}
1239
1240/*2
1241*returns the deg-th power of the ideal gid
1242*/
1243//ideal idPower(ideal gid,int deg)
1244//{
1245//  int i;
1246//  ideal id;
1247//
1248//  if (deg < 1) deg = 1;
1249//  i = binom(IDELEMS(gid)+deg-1,deg);
1250//  id=idInit(i,1);
1251//  idpower = id->m;
1252//  givenideal = gid->m;
1253//  idpowerpoint = 0;
1254//  makepotence(IDELEMS(gid),1,deg,0);
1255//  idpower = NULL;
1256//  givenideal = NULL;
1257//  idpowerpoint = 0;
1258//  return id;
1259//}
1260static void idNextPotence(ideal given, ideal result,
1261  int begin, int end, int deg, int restdeg, poly ap)
1262{
1263  poly p;
1264  int i;
1265
1266  p = pPower(pCopy(given->m[begin]),restdeg);
1267  i = result->nrows;
1268  result->m[i] = pMult(pCopy(ap),p);
1269//PrintS(".");
1270  (result->nrows)++;
1271  if (result->nrows >= IDELEMS(result))
1272  {
1273    pEnlargeSet(&(result->m),IDELEMS(result),16);
1274    IDELEMS(result) += 16;
1275  }
1276  if (begin == end) return;
1277  for (i=restdeg-1;i>0;i--)
1278  {
1279    p = pPower(pCopy(given->m[begin]),i);
1280    p = pMult(pCopy(ap),p);
1281    idNextPotence(given, result, begin+1, end, deg, restdeg-i, p);
1282    pDelete(&p);
1283  }
1284  idNextPotence(given, result, begin+1, end, deg, restdeg, ap);
1285}
1286
1287ideal id_Power(ideal given,int exp, const ring r)
1288{
1289  ideal result,temp;
1290  poly p1;
1291  int i;
1292
1293  if (idIs0(given)) return idInit(1,1);
1294  temp = id_Copy(given,r);
1295  idSkipZeroes(temp);
1296  i = binom(IDELEMS(temp)+exp-1,exp);
1297  result = idInit(i,1);
1298  result->nrows = 0;
1299//Print("ideal contains %d elements\n",i);
1300  p1=p_One(r);
1301  idNextPotence(temp,result,0,IDELEMS(temp)-1,exp,exp,p1);
1302  p_Delete(&p1,r);
1303  id_Delete(&temp,r);
1304  result->nrows = 1;
1305  id_DelEquals(result,r);
1306  idSkipZeroes(result);
1307  return result;
1308}
1309
1310/*2
1311* compute the which-th ar-minor of the matrix a
1312*/
1313poly idMinor(matrix a, int ar, unsigned long which, ideal R)
1314{
1315  int     i,j,k,size;
1316  unsigned long curr;
1317  int *rowchoise,*colchoise;
1318  BOOLEAN rowch,colch;
1319  ideal result;
1320  matrix tmp;
1321  poly p,q;
1322
1323  i = binom(a->rows(),ar);
1324  j = binom(a->cols(),ar);
1325
1326  rowchoise=(int *)omAlloc(ar*sizeof(int));
1327  colchoise=(int *)omAlloc(ar*sizeof(int));
1328  if ((i>512) || (j>512) || (i*j >512)) size=512;
1329  else size=i*j;
1330  result=idInit(size,1);
1331  tmp=mpNew(ar,ar);
1332  k = 0; /* the index in result*/
1333  curr = 0; /* index of current minor */
1334  idInitChoise(ar,1,a->rows(),&rowch,rowchoise);
1335  while (!rowch)
1336  {
1337    idInitChoise(ar,1,a->cols(),&colch,colchoise);
1338    while (!colch)
1339    {
1340      if (curr == which)
1341      {
1342        for (i=1; i<=ar; i++)
1343        {
1344          for (j=1; j<=ar; j++)
1345          {
1346            MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
1347          }
1348        }
1349        p = mpDetBareiss(tmp);
1350        if (p!=NULL)
1351        {
1352          if (R!=NULL)
1353          {
1354            q = p;
1355            p = kNF(R,currQuotient,q);
1356            pDelete(&q);
1357          }
1358          /*delete the matrix tmp*/
1359          for (i=1; i<=ar; i++)
1360          {
1361            for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
1362          }
1363          idDelete((ideal*)&tmp);
1364          omFreeSize((ADDRESS)rowchoise,ar*sizeof(int));
1365          omFreeSize((ADDRESS)colchoise,ar*sizeof(int));
1366          return (p);
1367        }
1368      }
1369      curr++;
1370      idGetNextChoise(ar,a->cols(),&colch,colchoise);
1371    }
1372    idGetNextChoise(ar,a->rows(),&rowch,rowchoise);
1373  }
1374  return (poly) 1;
1375}
1376
1377#ifdef WITH_OLD_MINOR
1378/*2
1379* compute all ar-minors of the matrix a
1380*/
1381ideal idMinors(matrix a, int ar, ideal R)
1382{
1383  int     i,j,k,size;
1384  int *rowchoise,*colchoise;
1385  BOOLEAN rowch,colch;
1386  ideal result;
1387  matrix tmp;
1388  poly p,q;
1389
1390  i = binom(a->rows(),ar);
1391  j = binom(a->cols(),ar);
1392
1393  rowchoise=(int *)omAlloc(ar*sizeof(int));
1394  colchoise=(int *)omAlloc(ar*sizeof(int));
1395  if ((i>512) || (j>512) || (i*j >512)) size=512;
1396  else size=i*j;
1397  result=idInit(size,1);
1398  tmp=mpNew(ar,ar);
1399  k = 0; /* the index in result*/
1400  idInitChoise(ar,1,a->rows(),&rowch,rowchoise);
1401  while (!rowch)
1402  {
1403    idInitChoise(ar,1,a->cols(),&colch,colchoise);
1404    while (!colch)
1405    {
1406      for (i=1; i<=ar; i++)
1407      {
1408        for (j=1; j<=ar; j++)
1409        {
1410          MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
1411        }
1412      }
1413      p = mpDetBareiss(tmp);
1414      if (p!=NULL)
1415      {
1416        if (R!=NULL)
1417        {
1418          q = p;
1419          p = kNF(R,currQuotient,q);
1420          pDelete(&q);
1421        }
1422        if (p!=NULL)
1423        {
1424          if (k>=size)
1425          {
1426            pEnlargeSet(&result->m,size,32);
1427            size += 32;
1428          }
1429          result->m[k] = p;
1430          k++;
1431        }
1432      }
1433      idGetNextChoise(ar,a->cols(),&colch,colchoise);
1434    }
1435    idGetNextChoise(ar,a->rows(),&rowch,rowchoise);
1436  }
1437  /*delete the matrix tmp*/
1438  for (i=1; i<=ar; i++)
1439  {
1440    for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
1441  }
1442  idDelete((ideal*)&tmp);
1443  if (k==0)
1444  {
1445    k=1;
1446    result->m[0]=NULL;
1447  }
1448  omFreeSize((ADDRESS)rowchoise,ar*sizeof(int));
1449  omFreeSize((ADDRESS)colchoise,ar*sizeof(int));
1450  pEnlargeSet(&result->m,size,k-size);
1451  IDELEMS(result) = k;
1452  return (result);
1453}
1454#else
1455/*2
1456* compute all ar-minors of the matrix a
1457* the caller of mpRecMin
1458* the elements of the result are not in R (if R!=NULL)
1459*/
1460ideal idMinors(matrix a, int ar, ideal R)
1461{
1462  int elems=0;
1463  int r=a->nrows,c=a->ncols;
1464  int i;
1465  matrix b;
1466  ideal result,h;
1467  ring origR;
1468  ring tmpR;
1469  long bound;
1470
1471  if((ar<=0) || (ar>r) || (ar>c))
1472  {
1473    Werror("%d-th minor, matrix is %dx%d",ar,r,c);
1474    return NULL;
1475  }
1476  h = idMatrix2Module(mpCopy(a));
1477  bound = smExpBound(h,c,r,ar);
1478  idDelete(&h);
1479  tmpR=smRingChange(&origR,bound);
1480  b = mpNew(r,c);
1481  for (i=r*c-1;i>=0;i--)
1482  {
1483    if (a->m[i])
1484      b->m[i] = prCopyR(a->m[i],origR);
1485  }
1486  if (R!=NULL)
1487  {
1488    R = idrCopyR(R,origR);
1489    //if (ar>1) // otherwise done in mpMinorToResult
1490    //{
1491    //  matrix bb=(matrix)kNF(R,currQuotient,(ideal)b);
1492    //  bb->rank=b->rank; bb->nrows=b->nrows; bb->ncols=b->ncols;
1493    //  idDelete((ideal*)&b); b=bb;
1494    //}
1495  }
1496  result=idInit(32,1);
1497  if(ar>1) mpRecMin(ar-1,result,elems,b,r,c,NULL,R);
1498  else mpMinorToResult(result,elems,b,r,c,R);
1499  idDelete((ideal *)&b);
1500  if (R!=NULL) idDelete(&R);
1501  idSkipZeroes(result);
1502  rChangeCurrRing(origR);
1503  result = idrMoveR(result,tmpR);
1504  smKillModifiedRing(tmpR);
1505  idTest(result);
1506  return result;
1507}
1508#endif
1509
1510/*2
1511*skips all zeroes and double elements, searches also for units
1512*/
1513void id_Compactify(ideal id, const ring r)
1514{
1515  int i,j;
1516  BOOLEAN b=FALSE;
1517
1518  i = IDELEMS(id)-1;
1519  while ((! b) && (i>=0))
1520  {
1521    b=p_IsUnit(id->m[i],r);
1522    i--;
1523  }
1524  if (b)
1525  {
1526    for(i=IDELEMS(id)-1;i>=0;i--) p_Delete(&id->m[i],r);
1527    id->m[0]=p_One(r);
1528  }
1529  else
1530  {
1531    id_DelMultiples(id,r);
1532  }
1533  idSkipZeroes(id);
1534}
1535
1536/*2
1537*returns TRUE if id1 is a submodule of id2
1538*/
1539BOOLEAN idIsSubModule(ideal id1,ideal id2)
1540{
1541  int i;
1542  poly p;
1543
1544  if (idIs0(id1)) return TRUE;
1545  for (i=0;i<IDELEMS(id1);i++)
1546  {
1547    if (id1->m[i] != NULL)
1548    {
1549      p = kNF(id2,currQuotient,id1->m[i]);
1550      if (p != NULL)
1551      {
1552        pDelete(&p);
1553        return FALSE;
1554      }
1555    }
1556  }
1557  return TRUE;
1558}
1559
1560/*2
1561* returns the ideals of initial terms
1562*/
1563ideal idHead(ideal h)
1564{
1565  ideal m = idInit(IDELEMS(h),h->rank);
1566  int i;
1567
1568  for (i=IDELEMS(h)-1;i>=0; i--)
1569  {
1570    if (h->m[i]!=NULL) m->m[i]=pHead(h->m[i]);
1571  }
1572  return m;
1573}
1574
1575ideal idHomogen(ideal h, int varnum)
1576{
1577  ideal m = idInit(IDELEMS(h),h->rank);
1578  int i;
1579
1580  for (i=IDELEMS(h)-1;i>=0; i--)
1581  {
1582    m->m[i]=pHomogen(h->m[i],varnum);
1583  }
1584  return m;
1585}
1586
1587/*------------------type conversions----------------*/
1588ideal idVec2Ideal(poly vec)
1589{
1590   ideal result=idInit(1,1);
1591   omFree((ADDRESS)result->m);
1592   result->m=NULL; // remove later
1593   pVec2Polys(vec, &(result->m), &(IDELEMS(result)));
1594   return result;
1595}
1596
1597#define NEW_STUFF
1598#ifndef NEW_STUFF
1599// converts mat to module, destroys mat
1600ideal idMatrix2Module(matrix mat)
1601{
1602  int mc=MATCOLS(mat);
1603  int mr=MATROWS(mat);
1604  ideal result = idInit(si_max(mc,1),si_max(mr,1));
1605  int i,j;
1606  poly h;
1607
1608  for(j=0;j<mc /*MATCOLS(mat)*/;j++) /* j is also index in result->m */
1609  {
1610    for (i=1;i<=mr /*MATROWS(mat)*/;i++)
1611    {
1612      h = MATELEM(mat,i,j+1);
1613      if (h!=NULL)
1614      {
1615        MATELEM(mat,i,j+1)=NULL;
1616        pSetCompP(h,i);
1617        result->m[j] = pAdd(result->m[j],h);
1618      }
1619    }
1620  }
1621  // obachman: need to clean this up
1622  idDelete((ideal*) &mat);
1623  return result;
1624}
1625#else
1626
1627#include "sbuckets.h"
1628
1629// converts mat to module, destroys mat
1630ideal idMatrix2Module(matrix mat)
1631{
1632  int mc=MATCOLS(mat);
1633  int mr=MATROWS(mat);
1634  ideal result = idInit(si_max(mc,1),si_max(mr,1));
1635  int i,j, l;
1636  poly h;
1637  poly p;
1638  sBucket_pt bucket = sBucketCreate(currRing);
1639
1640  for(j=0;j<mc /*MATCOLS(mat)*/;j++) /* j is also index in result->m */
1641  {
1642    for (i=1;i<=mr /*MATROWS(mat)*/;i++)
1643    {
1644      h = MATELEM(mat,i,j+1);
1645      if (h!=NULL)
1646      {
1647        l=pLength(h);
1648        MATELEM(mat,i,j+1)=NULL;
1649        p_SetCompP(h,i, currRing);
1650        sBucket_Merge_p(bucket, h, l);
1651      }
1652    }
1653    sBucketClearMerge(bucket, &(result->m[j]), &l);
1654  }
1655  sBucketDestroy(&bucket);
1656
1657  // obachman: need to clean this up
1658  idDelete((ideal*) &mat);
1659  return result;
1660}
1661#endif
1662
1663/*2
1664* converts a module into a matrix, destroyes the input
1665*/
1666matrix idModule2Matrix(ideal mod)
1667{
1668  matrix result = mpNew(mod->rank,IDELEMS(mod));
1669  int i,cp;
1670  poly p,h;
1671
1672  for(i=0;i<IDELEMS(mod);i++)
1673  {
1674    p=pReverse(mod->m[i]);
1675    mod->m[i]=NULL;
1676    while (p!=NULL)
1677    {
1678      h=p;
1679      pIter(p);
1680      pNext(h)=NULL;
1681//      cp = si_max(1,pGetComp(h));     // if used for ideals too
1682      cp = pGetComp(h);
1683      pSetComp(h,0);
1684      pSetmComp(h);
1685#ifdef TEST
1686      if (cp>mod->rank)
1687      {
1688        Print("## inv. rank %ld -> %d\n",mod->rank,cp);
1689        int k,l,o=mod->rank;
1690        mod->rank=cp;
1691        matrix d=mpNew(mod->rank,IDELEMS(mod));
1692        for (l=1; l<=o; l++)
1693        {
1694          for (k=1; k<=IDELEMS(mod); k++)
1695          {
1696            MATELEM(d,l,k)=MATELEM(result,l,k);
1697            MATELEM(result,l,k)=NULL;
1698          }
1699        }
1700        idDelete((ideal *)&result);
1701        result=d;
1702      }
1703#endif
1704      MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
1705    }
1706  }
1707  // obachman 10/99: added the following line, otherwise memory leack!
1708  idDelete(&mod);
1709  return result;
1710}
1711
1712matrix idModule2formatedMatrix(ideal mod,int rows, int cols)
1713{
1714  matrix result = mpNew(rows,cols);
1715  int i,cp,r=idRankFreeModule(mod),c=IDELEMS(mod);
1716  poly p,h;
1717
1718  if (r>rows) r = rows;
1719  if (c>cols) c = cols;
1720  for(i=0;i<c;i++)
1721  {
1722    p=pReverse(mod->m[i]);
1723    mod->m[i]=NULL;
1724    while (p!=NULL)
1725    {
1726      h=p;
1727      pIter(p);
1728      pNext(h)=NULL;
1729      cp = pGetComp(h);
1730      if (cp<=r)
1731      {
1732        pSetComp(h,0);
1733        pSetmComp(h);
1734        MATELEM(result,cp,i+1) = pAdd(MATELEM(result,cp,i+1),h);
1735      }
1736      else
1737        pDelete(&h);
1738    }
1739  }
1740  idDelete(&mod);
1741  return result;
1742}
1743
1744/*2
1745* substitute the n-th variable by the monomial e in id
1746* destroy id
1747*/
1748ideal  idSubst(ideal id, int n, poly e)
1749{
1750  int k=MATROWS((matrix)id)*MATCOLS((matrix)id);
1751  ideal res=(ideal)mpNew(MATROWS((matrix)id),MATCOLS((matrix)id));
1752
1753  res->rank = id->rank;
1754  for(k--;k>=0;k--)
1755  {
1756    res->m[k]=pSubst(id->m[k],n,e);
1757    id->m[k]=NULL;
1758  }
1759  idDelete(&id);
1760  return res;
1761}
1762
1763BOOLEAN idHomModule(ideal m, ideal Q, intvec **w)
1764{
1765  if (w!=NULL) *w=NULL;
1766  if ((Q!=NULL) && (!idHomIdeal(Q,NULL))) return FALSE;
1767  if (idIs0(m))
1768  {
1769    if (w!=NULL) (*w)=new intvec(m->rank);
1770    return TRUE;
1771  }
1772
1773  long cmax=1,order=0,ord,* diff,diffmin=32000;
1774  int *iscom;
1775  int i,j;
1776  poly p=NULL;
1777  pFDegProc d;
1778  if (pLexOrder && (currRing->order[0]==ringorder_lp))
1779     d=p_Totaldegree;
1780  else
1781     d=pFDeg;
1782  int length=IDELEMS(m);
1783  polyset P=m->m;
1784  polyset F=(polyset)omAlloc(length*sizeof(poly));
1785  for (i=length-1;i>=0;i--)
1786  {
1787    p=F[i]=P[i];
1788    cmax=si_max(cmax,(long)pMaxComp(p));
1789  }
1790  cmax++;
1791  diff = (long *)omAlloc0(cmax*sizeof(long));
1792  if (w!=NULL) *w=new intvec(cmax-1);
1793  iscom = (int *)omAlloc0(cmax*sizeof(int));
1794  i=0;
1795  while (i<=length)
1796  {
1797    if (i<length)
1798    {
1799      p=F[i];
1800      while ((p!=NULL) && (iscom[pGetComp(p)]==0)) pIter(p);
1801    }
1802    if ((p==NULL) && (i<length))
1803    {
1804      i++;
1805    }
1806    else
1807    {
1808      if (p==NULL) /* && (i==length) */
1809      {
1810        i=0;
1811        while ((i<length) && (F[i]==NULL)) i++;
1812        if (i>=length) break;
1813        p = F[i];
1814      }
1815      //if (pLexOrder && (currRing->order[0]==ringorder_lp))
1816      //  order=pTotaldegree(p);
1817      //else
1818      //  order = p->order;
1819      //  order = pFDeg(p,currRing);
1820      order = d(p,currRing) +diff[pGetComp(p)];
1821      //order += diff[pGetComp(p)];
1822      p = F[i];
1823//Print("Actual p=F[%d]: ",i);pWrite(p);
1824      F[i] = NULL;
1825      i=0;
1826    }
1827    while (p!=NULL)
1828    {
1829      if (pLexOrder && (currRing->order[0]==ringorder_lp))
1830        ord=pTotaldegree(p);
1831      else
1832      //  ord = p->order;
1833        ord = pFDeg(p,currRing);
1834      if (iscom[pGetComp(p)]==0)
1835      {
1836        diff[pGetComp(p)] = order-ord;
1837        iscom[pGetComp(p)] = 1;
1838/*
1839*PrintS("new diff: ");
1840*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
1841*PrintLn();
1842*PrintS("new iscom: ");
1843*for (j=0;j<cmax;j++) Print("%d ",iscom[j]);
1844*PrintLn();
1845*Print("new set %d, order %d, ord %d, diff %d\n",pGetComp(p),order,ord,diff[pGetComp(p)]);
1846*/
1847      }
1848      else
1849      {
1850/*
1851*PrintS("new diff: ");
1852*for (j=0;j<cmax;j++) Print("%d ",diff[j]);
1853*PrintLn();
1854*Print("order %d, ord %d, diff %d\n",order,ord,diff[pGetComp(p)]);
1855*/
1856        if (order != (ord+diff[pGetComp(p)]))
1857        {
1858          omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
1859          omFreeSize((ADDRESS) diff,cmax*sizeof(long));
1860          omFreeSize((ADDRESS) F,length*sizeof(poly));
1861          delete *w;*w=NULL;
1862          return FALSE;
1863        }
1864      }
1865      pIter(p);
1866    }
1867  }
1868  omFreeSize((ADDRESS) iscom,cmax*sizeof(int));
1869  omFreeSize((ADDRESS) F,length*sizeof(poly));
1870  for (i=1;i<cmax;i++) (**w)[i-1]=(int)(diff[i]);
1871  for (i=1;i<cmax;i++)
1872  {
1873    if (diff[i]<diffmin) diffmin=diff[i];
1874  }
1875  if (w!=NULL)
1876  {
1877    for (i=1;i<cmax;i++)
1878    {
1879      (**w)[i-1]=(int)(diff[i]-diffmin);
1880    }
1881  }
1882  omFreeSize((ADDRESS) diff,cmax*sizeof(long));
1883  return TRUE;
1884}
1885
1886BOOLEAN idTestHomModule(ideal m, ideal Q, intvec *w)
1887{
1888  if ((Q!=NULL) && (!idHomIdeal(Q,NULL)))  { PrintS(" Q not hom\n"); return FALSE;}
1889  if (idIs0(m)) return TRUE;
1890
1891  int cmax=-1;
1892  int i;
1893  poly p=NULL;
1894  int length=IDELEMS(m);
1895  polyset P=m->m;
1896  for (i=length-1;i>=0;i--)
1897  {
1898    p=P[i];
1899    if (p!=NULL) cmax=si_max(cmax,(int)pMaxComp(p)+1);
1900  }
1901  if (w != NULL)
1902  if (w->length()+1 < cmax)
1903  {
1904    // Print("length: %d - %d \n", w->length(),cmax);
1905    return FALSE;
1906  }
1907
1908  if(w!=NULL)
1909    pSetModDeg(w);
1910
1911  for (i=length-1;i>=0;i--)
1912  {
1913    p=P[i];
1914    poly q=p;
1915    if (p!=NULL)
1916    {
1917      int d=pFDeg(p,currRing);
1918      loop
1919      {
1920        pIter(p);
1921        if (p==NULL) break;
1922        if (d!=pFDeg(p,currRing))
1923        {
1924          //pWrite(q); wrp(p); Print(" -> %d - %d\n",d,pFDeg(p,currRing));
1925          if(w!=NULL)
1926            pSetModDeg(NULL);
1927          return FALSE;
1928        }
1929      }
1930    }
1931  }
1932
1933  if(w!=NULL)
1934    pSetModDeg(NULL);
1935
1936  return TRUE;
1937}
1938
1939ideal idJet(ideal i,int d)
1940{
1941  ideal r=idInit((i->nrows)*(i->ncols),i->rank);
1942  r->nrows = i-> nrows;
1943  r->ncols = i-> ncols;
1944  //r->rank = i-> rank;
1945  int k;
1946  for(k=(i->nrows)*(i->ncols)-1;k>=0; k--)
1947  {
1948    r->m[k]=ppJet(i->m[k],d);
1949  }
1950  return r;
1951}
1952
1953ideal idJetW(ideal i,int d, intvec * iv)
1954{
1955  ideal r=idInit(IDELEMS(i),i->rank);
1956  if (ecartWeights!=NULL)
1957  {
1958    WerrorS("cannot compute weighted jets now");
1959  }
1960  else
1961  {
1962    short *w=iv2array(iv);
1963    int k;
1964    for(k=0; k<IDELEMS(i); k++)
1965    {
1966      r->m[k]=ppJetW(i->m[k],d,w);
1967    }
1968    omFreeSize((ADDRESS)w,(rVar(r)+1)*sizeof(short));
1969  }
1970  return r;
1971}
1972
1973int idMinDegW(ideal M,intvec *w)
1974{
1975  int d=-1;
1976  for(int i=0;i<IDELEMS(M);i++)
1977  {
1978    int d0=pMinDeg(M->m[i],w);
1979    if(-1<d0&&(d0<d||d==-1))
1980      d=d0;
1981  }
1982  return d;
1983}
1984
1985ideal idSeries(int n,ideal M,matrix U,intvec *w)
1986{
1987  for(int i=IDELEMS(M)-1;i>=0;i--)
1988  {
1989    if(U==NULL)
1990      M->m[i]=pSeries(n,M->m[i],NULL,w);
1991    else
1992    {
1993      M->m[i]=pSeries(n,M->m[i],MATELEM(U,i+1,i+1),w);
1994      MATELEM(U,i+1,i+1)=NULL;
1995    }
1996  }
1997  if(U!=NULL)
1998    idDelete((ideal*)&U);
1999  return M;
2000}
2001
2002matrix idDiff(matrix i, int k)
2003{
2004  int e=MATCOLS(i)*MATROWS(i);
2005  matrix r=mpNew(MATROWS(i),MATCOLS(i));
2006  r->rank=i->rank;
2007  int j;
2008  for(j=0; j<e; j++)
2009  {
2010    r->m[j]=pDiff(i->m[j],k);
2011  }
2012  return r;
2013}
2014
2015matrix idDiffOp(ideal I, ideal J,BOOLEAN multiply)
2016{
2017  matrix r=mpNew(IDELEMS(I),IDELEMS(J));
2018  int i,j;
2019  for(i=0; i<IDELEMS(I); i++)
2020  {
2021    for(j=0; j<IDELEMS(J); j++)
2022    {
2023      MATELEM(r,i+1,j+1)=pDiffOp(I->m[i],J->m[j],multiply);
2024    }
2025  }
2026  return r;
2027}
2028
2029int idElem(const ideal F)
2030{
2031  int i=0,j=IDELEMS(F)-1;
2032
2033  while(j>=0)
2034  {
2035    if ((F->m)[j]!=NULL) i++;
2036    j--;
2037  }
2038  return i;
2039}
2040
2041/*
2042*computes module-weights for liftings of homogeneous modules
2043*/
2044intvec * idMWLift(ideal mod,intvec * weights)
2045{
2046  if (idIs0(mod)) return new intvec(2);
2047  int i=IDELEMS(mod);
2048  while ((i>0) && (mod->m[i-1]==NULL)) i--;
2049  intvec *result = new intvec(i+1);
2050  while (i>0)
2051  {
2052    (*result)[i]=pFDeg(mod->m[i],currRing)+(*weights)[pGetComp(mod->m[i])];
2053  }
2054  return result;
2055}
2056
2057/*2
2058*sorts the kbase for idCoef* in a special way (lexicographically
2059*with x_max,...,x_1)
2060*/
2061ideal idCreateSpecialKbase(ideal kBase,intvec ** convert)
2062{
2063  int i;
2064  ideal result;
2065
2066  if (idIs0(kBase)) return NULL;
2067  result = idInit(IDELEMS(kBase),kBase->rank);
2068  *convert = idSort(kBase,FALSE);
2069  for (i=0;i<(*convert)->length();i++)
2070  {
2071    result->m[i] = pCopy(kBase->m[(**convert)[i]-1]);
2072  }
2073  return result;
2074}
2075
2076/*2
2077*returns the index of a given monom in the list of the special kbase
2078*/
2079int idIndexOfKBase(poly monom, ideal kbase)
2080{
2081  int j=IDELEMS(kbase);
2082
2083  while ((j>0) && (kbase->m[j-1]==NULL)) j--;
2084  if (j==0) return -1;
2085  int i=rVar(r);
2086  while (i>0)
2087  {
2088    loop
2089    {
2090      if (pGetExp(monom,i)>pGetExp(kbase->m[j-1],i)) return -1;
2091      if (pGetExp(monom,i)==pGetExp(kbase->m[j-1],i)) break;
2092      j--;
2093      if (j==0) return -1;
2094    }
2095    if (i==1)
2096    {
2097      while(j>0)
2098      {
2099        if (pGetComp(monom)==pGetComp(kbase->m[j-1])) return j-1;
2100        if (pGetComp(monom)>pGetComp(kbase->m[j-1])) return -1;
2101        j--;
2102      }
2103    }
2104    i--;
2105  }
2106  return -1;
2107}
2108
2109/*2
2110*decomposes the monom in a part of coefficients described by the
2111*complement of how and a monom in variables occuring in how, the
2112*index of which in kbase is returned as integer pos (-1 if it don't
2113*exists)
2114*/
2115poly idDecompose(poly monom, poly how, ideal kbase, int * pos)
2116{
2117  int i;
2118  poly coeff=p_One(r), base=p_One(r);
2119
2120  for (i=1;i<=rVar(r);i++)
2121  {
2122    if (pGetExp(how,i)>0)
2123    {
2124      pSetExp(base,i,pGetExp(monom,i));
2125    }
2126    else
2127    {
2128      pSetExp(coeff,i,pGetExp(monom,i));
2129    }
2130  }
2131  pSetComp(base,pGetComp(monom));
2132  pSetm(base);
2133  pSetCoeff(coeff,nCopy(pGetCoeff(monom)));
2134  pSetm(coeff);
2135  *pos = idIndexOfKBase(base,kbase);
2136  if (*pos<0)
2137    pDelete(&coeff);
2138  pDelete(&base);
2139  return coeff;
2140}
2141
2142/*2
2143*returns a matrix A of coefficients with kbase*A=arg
2144*if all monomials in variables of how occur in kbase
2145*the other are deleted
2146*/
2147matrix idCoeffOfKBase(ideal arg, ideal kbase, poly how)
2148{
2149  matrix result;
2150  ideal tempKbase;
2151  poly p,q;
2152  intvec * convert;
2153  int i=IDELEMS(kbase),j=IDELEMS(arg),k,pos;
2154#if 0
2155  while ((i>0) && (kbase->m[i-1]==NULL)) i--;
2156  if (idIs0(arg))
2157    return mpNew(i,1);
2158  while ((j>0) && (arg->m[j-1]==NULL)) j--;
2159  result = mpNew(i,j);
2160#else
2161  result = mpNew(i, j);
2162  while ((j>0) && (arg->m[j-1]==NULL)) j--;
2163#endif
2164
2165  tempKbase = idCreateSpecialKbase(kbase,&convert);
2166  for (k=0;k<j;k++)
2167  {
2168    p = arg->m[k];
2169    while (p!=NULL)
2170    {
2171      q = idDecompose(p,how,tempKbase,&pos);
2172      if (pos>=0)
2173      {
2174        MATELEM(result,(*convert)[pos],k+1) =
2175            pAdd(MATELEM(result,(*convert)[pos],k+1),q);
2176      }
2177      else
2178        pDelete(&q);
2179      pIter(p);
2180    }
2181  }
2182  idDelete(&tempKbase);
2183  return result;
2184}
2185
2186/*3
2187* searches for the next unit in the components of the module arg and
2188* returns the first one;
2189*/
2190static int id_ReadOutPivot(ideal arg,int* comp, const ring r)
2191{
2192  if (idIs0(arg)) return -1;
2193  int i=0,j, generator=-1;
2194  int rk_arg=arg->rank; //idRankFreeModule(arg);
2195  int * componentIsUsed =(int *)omAlloc((rk_arg+1)*sizeof(int));
2196  poly p;
2197
2198  while ((generator<0) && (i<IDELEMS(arg)))
2199  {
2200    memset(componentIsUsed,0,(rk_arg+1)*sizeof(int));
2201    p = arg->m[i];
2202    while (p!=NULL)
2203    {
2204      j = p_GetComp(p,r);
2205      if (componentIsUsed[j]==0)
2206      {
2207#ifdef HAVE_RINGS
2208        if (p_LmIsConstantComp(p,r) &&
2209            (!rField_is_Ring(r) || n_IsUnit(pGetCoeff(p),r->cf)))
2210        {
2211#else
2212        if (p_LmIsConstantComp(p,r))
2213        {
2214#endif
2215          generator = i;
2216          componentIsUsed[j] = 1;
2217        }
2218        else
2219        {
2220          componentIsUsed[j] = -1;
2221        }
2222      }
2223      else if (componentIsUsed[j]>0)
2224      {
2225        (componentIsUsed[j])++;
2226      }
2227      pIter(p);
2228    }
2229    i++;
2230  }
2231  i = 0;
2232  *comp = -1;
2233  for (j=0;j<=rk_arg;j++)
2234  {
2235    if (componentIsUsed[j]>0)
2236    {
2237      if ((*comp==-1) || (componentIsUsed[j]<i))
2238      {
2239        *comp = j;
2240        i= componentIsUsed[j];
2241      }
2242    }
2243  }
2244  omFree(componentIsUsed);
2245  return generator;
2246}
2247
2248#if 0
2249static void idDeleteComp(ideal arg,int red_comp)
2250{
2251  int i,j;
2252  poly p;
2253
2254  for (i=IDELEMS(arg)-1;i>=0;i--)
2255  {
2256    p = arg->m[i];
2257    while (p!=NULL)
2258    {
2259      j = pGetComp(p);
2260      if (j>red_comp)
2261      {
2262        pSetComp(p,j-1);
2263        pSetm(p);
2264      }
2265      pIter(p);
2266    }
2267  }
2268  (arg->rank)--;
2269}
2270#endif
2271
2272static void idDeleteComps(ideal arg,int* red_comp,int del)
2273// red_comp is an array [0..args->rank]
2274{
2275  int i,j;
2276  poly p;
2277
2278  for (i=IDELEMS(arg)-1;i>=0;i--)
2279  {
2280    p = arg->m[i];
2281    while (p!=NULL)
2282    {
2283      j = pGetComp(p);
2284      if (red_comp[j]!=j)
2285      {
2286        pSetComp(p,red_comp[j]);
2287        pSetmComp(p);
2288      }
2289      pIter(p);
2290    }
2291  }
2292  (arg->rank) -= del;
2293}
2294
2295/*2
2296* returns the presentation of an isomorphic, minimally
2297* embedded  module (arg represents the quotient!)
2298*/
2299ideal idMinEmbedding(ideal arg,BOOLEAN inPlace, intvec **w, const ring r)
2300{
2301  if (idIs0(arg)) return idInit(1,arg->rank);
2302  int i,next_gen,next_comp;
2303  ideal res=arg;
2304  if (!inPlace) res = id_Copy(arg,r);
2305  res->rank=si_max(res->rank,id_RankFreeModule(res,r));
2306  int *red_comp=(int*)omAlloc((res->rank+1)*sizeof(int));
2307  for (i=res->rank;i>=0;i--) red_comp[i]=i;
2308
2309  int del=0;
2310  loop
2311  {
2312    next_gen = idReadOutPivot(res,&next_comp);
2313    if (next_gen<0) break;
2314    del++;
2315    syGaussForOne(res,next_gen,next_comp,0,IDELEMS(res));
2316    for(i=next_comp+1;i<=arg->rank;i++) red_comp[i]--;
2317    if ((w !=NULL)&&(*w!=NULL))
2318    {
2319      for(i=next_comp;i<(*w)->length();i++) (**w)[i-1]=(**w)[i];
2320    }
2321  }
2322
2323  idDeleteComps(res,red_comp,del);
2324  idSkipZeroes(res);
2325  omFree(red_comp);
2326
2327  if ((w !=NULL)&&(*w!=NULL) &&(del>0))
2328  {
2329    intvec *wtmp=new intvec((*w)->length()-del);
2330    for(i=0;i<res->rank;i++) (*wtmp)[i]=(**w)[i];
2331    delete *w;
2332    *w=wtmp;
2333  }
2334  return res;
2335}
2336
2337intvec * idQHomWeight(ideal id)
2338{
2339  poly head, tail;
2340  int k;
2341  int in=IDELEMS(id)-1, ready=0, all=0,
2342      coldim=rVar(r), rowmax=2*coldim;
2343  if (in<0) return NULL;
2344  intvec *imat=new intvec(rowmax+1,coldim,0);
2345
2346  do
2347  {
2348    head = id->m[in--];
2349    if (head!=NULL)
2350    {
2351      tail = pNext(head);
2352      while (tail!=NULL)
2353      {
2354        all++;
2355        for (k=1;k<=coldim;k++)
2356          IMATELEM(*imat,all,k) = pGetExpDiff(head,tail,k);
2357        if (all==rowmax)
2358        {
2359          ivTriangIntern(imat, ready, all);
2360          if (ready==coldim)
2361          {
2362            delete imat;
2363            return NULL;
2364          }
2365        }
2366        pIter(tail);
2367      }
2368    }
2369  } while (in>=0);
2370  if (all>ready)
2371  {
2372    ivTriangIntern(imat, ready, all);
2373    if (ready==coldim)
2374    {
2375      delete imat;
2376      return NULL;
2377    }
2378  }
2379  intvec *result = ivSolveKern(imat, ready);
2380  delete imat;
2381  return result;
2382}
2383
2384BOOLEAN idIsZeroDim(ideal I)
2385{
2386  BOOLEAN *UsedAxis=(BOOLEAN *)omAlloc0(rVar(r)*sizeof(BOOLEAN));
2387  int i,n;
2388  poly po;
2389  BOOLEAN res=TRUE;
2390  for(i=IDELEMS(I)-1;i>=0;i--)
2391  {
2392    po=I->m[i];
2393    if ((po!=NULL) &&((n=pIsPurePower(po))!=0)) UsedAxis[n-1]=TRUE;
2394  }
2395  for(i=rVar(r)-1;i>=0;i--)
2396  {
2397    if(UsedAxis[i]==FALSE) {res=FALSE; break;} // not zero-dim.
2398  }
2399  omFreeSize(UsedAxis,rVar(r)*sizeof(BOOLEAN));
2400  return res;
2401}
2402
2403void id_Normalize(ideal I,const ring r)
2404{
2405  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
2406  int i;
2407  for(i=IDELEMS(I)-1;i>=0;i--)
2408  {
2409    p_Normalize(I->m[i],r);
2410  }
2411}
2412
2413// #include <kernel/clapsing.h>
2414
2415#ifdef HAVE_FACTORY
2416poly id_GCD(poly f, poly g, const ring r)
2417{
2418  ring save_r=r;
2419  rChangeCurrRing(r);
2420  ideal I=idInit(2,1); I->m[0]=f; I->m[1]=g;
2421  intvec *w = NULL;
2422  ideal S=idSyzygies(I,testHomog,&w);
2423  if (w!=NULL) delete w;
2424  poly gg=pTakeOutComp(&(S->m[0]),2);
2425  idDelete(&S);
2426  poly gcd_p=singclap_pdivide(f,gg);
2427  pDelete(&gg);
2428  rChangeCurrRing(save_r);
2429  return gcd_p;
2430}
2431#endif
2432
2433/*2
2434* xx,q: arrays of length 0..rl-1
2435* xx[i]: SB mod q[i]
2436* assume: char=0
2437* assume: q[i]!=0
2438* destroys xx
2439*/
2440#ifdef HAVE_FACTORY
2441ideal idChineseRemainder(ideal *xx, number *q, int rl, const ring r)
2442{
2443  int cnt=IDELEMS(xx[0])*xx[0]->nrows;
2444  ideal result=idInit(cnt,xx[0]->rank);
2445  result->nrows=xx[0]->nrows; // for lifting matrices
2446  result->ncols=xx[0]->ncols; // for lifting matrices
2447  int i,j;
2448  poly r,h,hh,res_p;
2449  number *x=(number *)omAlloc(rl*sizeof(number));
2450  for(i=cnt-1;i>=0;i--)
2451  {
2452    res_p=NULL;
2453    loop
2454    {
2455      r=NULL;
2456      for(j=rl-1;j>=0;j--)
2457      {
2458        h=xx[j]->m[i];
2459        if ((h!=NULL)
2460        &&((r==NULL)||(pLmCmp(r,h)==-1)))
2461          r=h;
2462      }
2463      if (r==NULL) break;
2464      h=pHead(r);
2465      for(j=rl-1;j>=0;j--)
2466      {
2467        hh=xx[j]->m[i];
2468        if ((hh!=NULL) && (pLmCmp(r,hh)==0))
2469        {
2470          x[j]=pGetCoeff(hh);
2471          hh=pLmFreeAndNext(hh);
2472          xx[j]->m[i]=hh;
2473        }
2474        else
2475          x[j]=nlInit(0, r);
2476      }
2477      number n=nlChineseRemainder(x,q,rl);
2478      for(j=rl-1;j>=0;j--)
2479      {
2480        x[j]=NULL; // nlInit(0...) takes no memory
2481      }
2482      if (nlIsZero(n)) pDelete(&h);
2483      else
2484      {
2485        pSetCoeff(h,n);
2486        //Print("new mon:");pWrite(h);
2487        res_p=pAdd(res_p,h);
2488      }
2489    }
2490    result->m[i]=res_p;
2491  }
2492  omFree(x);
2493  for(i=rl-1;i>=0;i--) idDelete(&(xx[i]));
2494  omFree(xx);
2495  return result;
2496}
2497#endif
2498/* currently unsed:
2499ideal idChineseRemainder(ideal *xx, intvec *iv)
2500{
2501  int rl=iv->length();
2502  number *q=(number *)omAlloc(rl*sizeof(number));
2503  int i;
2504  for(i=0; i<rl; i++)
2505  {
2506    q[i]=nInit((*iv)[i]);
2507  }
2508  return idChineseRemainder(xx,q,rl);
2509}
2510*/
2511/*
2512 * lift ideal with coeffs over Z (mod N) to Q via Farey
2513 */
2514ideal id_Farey(ideal x, number N, const ring r)
2515{
2516  int cnt=IDELEMS(x)*x->nrows;
2517  ideal result=idInit(cnt,x->rank);
2518  result->nrows=x->nrows; // for lifting matrices
2519  result->ncols=x->ncols; // for lifting matrices
2520
2521  int i;
2522  for(i=cnt-1;i>=0;i--)
2523  {
2524    poly h=p_Copy(x->m[i],r);
2525    result->m[i]=h;
2526    while(h!=NULL)
2527    {
2528      number c=pGetCoeff(h);
2529      pSetCoeff0(h,nlFarey(c,N));
2530      n_Delete(&c,r->cf);
2531      pIter(h);
2532    }
2533    while((result->m[i]!=NULL)&&(n_IsZero(pGetCoeff(result->m[i]),r->cf)))
2534    {
2535      p_LmDelete(&(result->m[i]),r);
2536    }
2537    h=result->m[i];
2538    while((h!=NULL) && (pNext(h)!=NULL))
2539    {
2540      if(n_IsZero(pGetCoeff(pNext(h)),r->cf))
2541      {
2542        p_LmDelete(&pNext(h),r);
2543      }
2544      else pIter(h);
2545    }
2546  }
2547  return result;
2548}
2549
2550/*2
2551* transpose a module
2552*/
2553ideal id_Transp(ideal a, const ring rRing)
2554{
2555  int r = a->rank, c = IDELEMS(a);
2556  ideal b =  idInit(r,c);
2557
2558  for (int i=c; i>0; i--)
2559  {
2560    poly p=a->m[i-1];
2561    while(p!=NULL)
2562    {
2563      poly h=p_Head(p, rRing);
2564      int co=p_GetComp(h, rRing)-1;
2565      p_SetComp(h, i, rRing);
2566      p_Setm(h, rRing);
2567      b->m[co] = p_Add_q(b->m[co], h, rRing);
2568      pIter(p);
2569    }
2570  }
2571  return b;
2572}
2573
2574
2575
2576/*2
2577* The following is needed to compute the image of certain map used in
2578* the computation of cohomologies via BGG
2579* let M = { w_1, ..., w_k }, k = size(M) == ncols(M), n = nvars(currRing).
2580* assuming that nrows(M) <= m*n; the procedure computes:
2581* transpose(M) * transpose( var(1) I_m | ... | var(n) I_m ) :== transpose(module{f_1, ... f_k}),
2582* where f_i = \sum_{j=1}^{m} (w_i, v_j) gen(j),  (w_i, v_j) is a `scalar` multiplication.
2583* that is, if w_i = (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m) then
2584
2585  (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m)
2586*  var_1  ... var_1  |  var_2  ...  var_2  | ... |  var_n  ...  var(n)
2587*  gen_1  ... gen_m  |  gen_1  ...  gen_m  | ... |  gen_1  ...  gen_m
2588+ =>
2589  f_i =
2590
2591   a^1_1 * var(1) * gen(1) + ... + a^1_m * var(1) * gen(m) +
2592   a^2_1 * var(2) * gen(1) + ... + a^2_m * var(2) * gen(m) +
2593                             ...
2594   a^n_1 * var(n) * gen(1) + ... + a^n_m * var(n) * gen(m);
2595
2596   NOTE: for every f_i we run only ONCE along w_i saving partial sums into a temporary array of polys of size m
2597*/
2598ideal id_TensorModuleMult(const int m, const ideal M, const ring rRing)
2599{
2600// #ifdef DEBU
2601//  WarnS("tensorModuleMult!!!!");
2602
2603  assume(m > 0);
2604  assume(M != NULL);
2605
2606  const int n = rRing->N;
2607
2608  assume(M->rank <= m * n);
2609
2610  const int k = IDELEMS(M);
2611
2612  ideal idTemp =  idInit(k,m); // = {f_1, ..., f_k }
2613
2614  for( int i = 0; i < k; i++ ) // for every w \in M
2615  {
2616    poly pTempSum = NULL;
2617
2618    poly w = M->m[i];
2619
2620    while(w != NULL) // for each term of w...
2621    {
2622      poly h = p_Head(w, rRing);
2623
2624      const int  gen = p_GetComp(h, rRing); // 1 ...
2625
2626      assume(gen > 0);
2627      assume(gen <= n*m);
2628
2629      // TODO: write a formula with %, / instead of while!
2630      /*
2631      int c = gen;
2632      int v = 1;
2633      while(c > m)
2634      {
2635        c -= m;
2636        v++;
2637      }
2638      */
2639
2640      int cc = gen % m;
2641      if( cc == 0) cc = m;
2642      int vv = 1 + (gen - cc) / m;
2643
2644//      assume( cc == c );
2645//      assume( vv == v );
2646
2647      //  1<= c <= m
2648      assume( cc > 0 );
2649      assume( cc <= m );
2650
2651      assume( vv > 0 );
2652      assume( vv <= n );
2653
2654      assume( (cc + (vv-1)*m) == gen );
2655
2656      p_IncrExp(h, vv, rRing); // h *= var(j) && //      p_AddExp(h, vv, 1, rRing);
2657      p_SetComp(h, cc, rRing);
2658
2659      p_Setm(h, rRing);         // addjust degree after the previous steps!
2660
2661      pTempSum = p_Add_q(pTempSum, h, rRing); // it is slow since h will be usually put to the back of pTempSum!!!
2662
2663      pIter(w);
2664    }
2665
2666    idTemp->m[i] = pTempSum;
2667  }
2668
2669  // simplify idTemp???
2670
2671  ideal idResult = id_Transp(idTemp, rRing);
2672
2673  id_Delete(&idTemp, rRing);
2674
2675  return(idResult);
2676}
Note: See TracBrowser for help on using the repository browser.