source: git/Singular/syz1.cc @ 2f12a6f

spielwiese
Last change on this file since 2f12a6f was 2f12a6f, checked in by Hans Schönemann <hannes@…>, 27 years ago
* hannes: ipshell.cc: fixed a bug in iiExport (multiple exports in rings) subexpr.cc: implemented sleftv::Copy for ring/qring syz1.cc: added missing idSkipZeros to LaScala* git-svn-id: file:///usr/local/Singular/svn/trunk@298 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 59.6 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: syz1.cc,v 1.7 1997-05-18 11:13:14 Singular Exp $ */
5/*
6* ABSTRACT: resolutions
7*/
8
9
10#include "mod2.h"
11#include "mmemory.h"
12#include "polys.h"
13#include "febase.h"
14#include "kstd1.h"
15#include "kutil.h"
16#include "spolys.h"
17#include "spolys0.h"
18#include "stairc.h"
19#include "ipid.h"
20#include "cntrlc.h"
21#include "ipid.h"
22#include "intvec.h"
23#include "ipshell.h"
24#include "limits.h"
25#include "tok.h"
26#include "numbers.h"
27#include "ideals.h"
28#include "intvec.h"
29#include "ring.h"
30#include "syz.h"
31
32/*--------------static variables------------------------*/
33/*---contains the real components wrt. cdp--------------*/
34static int ** truecomponents=NULL;
35static int ** backcomponents=NULL;
36static int ** Howmuch=NULL;
37static int ** Firstelem=NULL;
38/*---points to the real components of the actual module-*/
39static int *  currcomponents=NULL;
40/*---head-term-polynomials for the reduction------------*/
41static poly redpol=NULL;
42/*---counts number of applications of GM-criteria-------*/
43//static int crit;
44//static int zeroRed;
45//static int dsim;
46//static int simple;
47static int euler;
48/*---controls the ordering------------------------------*/
49static intvec * orderedcomp;
50static int *binomials;
51static int highdeg;
52static int highdeg_1;
53
54/*3
55* deletes all entres of a pair
56*/
57static void syDeletePair(SObject * so)
58{
59  pDelete(&(*so).p);
60  pDelete(&(*so).lcm);
61  pDelete(&(*so).syz);
62  (*so).p1 = NULL;
63  (*so).p2 = NULL;
64  (*so).ind1 = 0;
65  (*so).ind2 = 0;
66  (*so).order = 0;
67  (*so).isNotMinimal = 0;
68}
69
70/*3
71* puts all entres of a pair to another
72*/
73static void syCopyPair(SObject * argso, SObject * imso)
74{
75  (*imso).p = (*argso).p;
76  (*imso).p1 = (*argso).p1;
77  (*imso).p2 = (*argso).p2;
78  (*imso).lcm = (*argso).lcm;
79  (*imso).syz = (*argso).syz;
80  (*imso).ind1 = (*argso).ind1;
81  (*imso).ind2 = (*argso).ind2;
82  (*imso).order = (*argso).order;
83  (*imso).isNotMinimal = (*argso).isNotMinimal;
84  (*argso).p = NULL;
85  (*argso).p1 = NULL;
86  (*argso).p2 = NULL;
87  (*argso).lcm = NULL;
88  (*argso).syz = NULL;
89  (*argso).ind1 = 0;
90  (*argso).ind2 = 0;
91  (*argso).isNotMinimal = 0;
92}
93
94/*3
95* deletes empty objects from a pair set beginning with
96* pair first
97* assumes a pair to be empty if .lcm does so
98*/
99static void syCompactifyPairSet(SSet sPairs, int sPlength, int first)
100{
101  int k=first,kk=0;
102
103  while (k+kk<sPlength)
104  {
105    if (sPairs[k+kk].lcm!=NULL)
106    {
107      if (kk>0) syCopyPair(&sPairs[k+kk],&sPairs[k]);
108      k++;
109    }
110    else
111    {
112      kk++;
113    }
114  }
115}
116
117/*3
118* deletes empty objects from a pair set beginning with
119* pair first
120* assumes a pair to be empty if .lcm does so
121*/
122static void syCompactify1(SSet sPairs, int* sPlength, int first)
123{
124  int k=first,kk=0;
125
126  while (k+kk<=*sPlength)
127  {
128    if (sPairs[k+kk].lcm!=NULL)
129    {
130      if (kk>0) syCopyPair(&sPairs[k+kk],&sPairs[k]);
131      k++;
132    }
133    else
134    {
135      kk++;
136    }
137  }
138  *sPlength -= kk;
139}
140
141/*3
142* replaces comp1dpc during homogeneous syzygy-computations
143* compares with components of currcomponents instead of the
144* exp[0]
145*/
146static int syzcomp1dpc(poly p1, poly p2)
147{
148  /*4 compare monomials by order then revlex*/
149    int i = pVariables;
150    if ((p1->exp[i] == p2->exp[i]))
151    {
152      do
153      {
154        i--;
155        if (i <= 1)
156        {
157          //(*orderingdepth)[pVariables-i]++;
158           /*4 handle module case:*/
159           if (p1->exp[0]==p2->exp[0]) return 0;
160           else if 
161              (currcomponents[p1->exp[0]]>currcomponents[p2->exp[0]]) 
162                return 1;
163           else return -1;
164        }
165      } while ((p1->exp[i] == p2->exp[i]));
166    }
167    //(*orderingdepth)[pVariables-i]++;
168    if (p1->exp[i] < p2->exp[i]) return 1;
169    return -1;
170}
171
172/*3
173* replaces comp1dpc during homogeneous syzygy-computations
174* compares with components of currcomponents instead of the
175* exp[0]
176*/
177static int syzcomp2dpc(poly p1, poly p2)
178{
179  int o1=pGetOrder(p1), o2=pGetOrder(p2);
180  if (o1 > o2) return 1;
181  if (o1 < o2) return -1;
182
183  if (o1>0)
184  {
185    int i = pVariables;
186    while ((i>1) && (p1->exp[i]==p2->exp[i]))
187      i--;
188    //(*orderingdepth)[pVariables-i]++;
189    if (i>1)
190    {
191      if (p1->exp[i] < p2->exp[i]) return 1;
192      return -1;
193    }
194  }
195  o1=p1->exp[0];
196  o2=p2->exp[0];
197  if (o1==o2/*p1->exp[0]==p2->exp[0]*/) return 0;
198  if (currcomponents[o1]>currcomponents[o2]) return 1;
199  return -1;
200}
201
202/*3
203* compares only the monomial without component
204*/
205static int syzcompmonomdp(poly p1, poly p2)
206{
207  int i;
208
209  /*4 compare monomials by order then revlex*/
210  if (pGetOrder(p1) == pGetOrder(p2))
211  {
212    i = pVariables;
213    if ((p1->exp[i] == p2->exp[i]))
214    {
215      do
216      {
217        i--;
218        if (i <= 1)
219          return 0;
220      } while ((p1->exp[i] == p2->exp[i]));
221    }
222    if (p1->exp[i] < p2->exp[i]) return 1;
223    return -1;
224  }
225  else if (pGetOrder(p1) > pGetOrder(p2)) return 1;
226  return -1;
227}
228
229/*
230* (i,j) in binomials[j-1,i-j+1]
231* table size is: pVariables * (highdeg+1)
232* highdeg_1==highdeg+1
233*/
234static void syBinomSet()
235{
236  highdeg_1=highdeg+1;
237  for(int j=1;j<=highdeg;j++)
238  {
239    binomials[j/*0,j*/] = j;
240    for (int i=1;i<pVariables;i++)
241    {
242      binomials[i*(highdeg_1)+j/*i,j*/]
243      = binomials[(i-1)*(highdeg_1)+j/*i-1,j*/]*(j+i)/(i+1);
244    }
245  }
246  for (int i=0;i<pVariables;i++)
247  { 
248    binomials[i*(highdeg_1)/*i,0*/]=0;
249  }
250}
251
252static inline int syBinom(int i,int j)
253{
254  return binomials[(j-1)*(highdeg_1)+i-j+1/*j-1,i-j*/];
255}
256
257static void syzSetm(poly p)
258{
259  int i = 0,j;
260  for (j=pVariables;j>0;j--)
261    i += p->exp[j];
262
263  if (i<=highdeg)
264  {
265    i=1;
266    j = -INT_MAX;
267    int *ip=binomials+pGetExp(p,1);
268    loop
269    {
270      j += (*ip);
271      if (i==pVariables) break;
272      i++;
273      ip+=highdeg_1+pGetExp(p,i);
274    }
275    pGetOrder(p) = j;
276  }
277  else
278    pGetOrder(p)=i;
279}
280
281static inline poly syMultT(poly p,poly m)
282{
283  poly q,result=q=pNew();
284  int j;
285 
286  if (pGetOrder(p)>0)
287  {
288    loop
289    {
290      spMemcpy(q,p);
291      for (j=pVariables;j>0;j--)
292        pGetExp(q,j) += pGetExp(m,j);
293      pSetCoeff0(q,nCopy(pGetCoeff(p)));
294      syzSetm(q);
295      pIter(p);
296      if (p!=NULL)
297      {
298        pNext(q) = pNew();
299        pIter(q);
300      }
301      else break;
302    }
303  }
304  else
305  {
306    poly lastmon=NULL;
307    int i=0;
308    loop
309    {
310      if (pGetOrder(p)!=i)
311      {
312        spMemcpy(q,p);
313        for (j=pVariables;j>0;j--)
314          pGetExp(q,j) += pGetExp(m,j);
315        syzSetm(q);
316        lastmon = q;
317        i = pGetOrder(p);
318      }
319      else
320      {
321        spMemcpy(q,lastmon);
322        pSetComp(q,pGetComp(p));
323      }
324      pSetCoeff0(q,nCopy(pGetCoeff(p)));
325      pIter(p);
326      if (p!=NULL)
327      {
328        pNext(q) = pNew();
329        pIter(q);
330      }
331      else break;
332    }
333  }
334  pNext(q) = NULL;
335  return result;
336}
337
338static inline poly syMultTNeg(poly p,poly m)
339{
340  poly q,result=q=pNew();
341  int j;
342 
343  if (pGetOrder(p)>0)
344  {
345    loop
346    {
347      spMemcpy(q,p);
348      for (j=pVariables;j>0;j--)
349        pGetExp(q,j) += pGetExp(m,j);
350      pSetCoeff0(q,nCopy(pGetCoeff(p)));
351      pGetCoeff(q) = nNeg(pGetCoeff(q));
352      syzSetm(q);
353      pIter(p);
354      if (p!=NULL)
355      {
356        pNext(q) = pNew();
357        pIter(q);
358      }
359      else break;
360    }
361  }
362  else
363  {
364    poly lastmon=NULL;
365    int i=0;
366    loop
367    {
368      if (pGetOrder(p)!=i)
369      {
370        spMemcpy(q,p);
371        for (j=pVariables;j>0;j--)
372          pGetExp(q,j) += pGetExp(m,j);
373        syzSetm(q);
374        lastmon = q;
375        i = pGetOrder(p);
376      }
377      else
378      {
379        spMemcpy(q,lastmon);
380        pSetComp(q,pGetComp(p));
381      }
382      pSetCoeff0(q,nCopy(pGetCoeff(p)));
383      pGetCoeff(q) = nNeg(pGetCoeff(q));
384      pIter(p);
385      if (p!=NULL)
386      {
387        pNext(q) = pNew();
388        pIter(q);
389      }
390      else break;
391    }
392  }
393  pNext(q) = NULL;
394  return result;
395}
396
397static poly syMultT1(poly p,poly m)
398{
399  poly q,result=q=pNew();
400  int j;
401 
402  if (pGetOrder(p)>0)
403  {
404    loop
405    {
406      spMemcpy(q,p);
407      for (j=pVariables;j>0;j--)
408        pGetExp(q,j) += pGetExp(m,j);
409      pSetCoeff0(q,nMult(pGetCoeff(p),pGetCoeff(m)));
410      syzSetm(q);
411      pIter(p);
412      if (p!=NULL)
413      {
414        pNext(q) = pNew();
415        pIter(q);
416      }
417      else break;
418    }
419  }
420  else
421  {
422    poly lastmon=NULL;
423    int i=0;
424    loop
425    {
426      if (pGetOrder(p)!=i)
427      {
428        spMemcpy(q,p);
429        for (j=pVariables;j>0;j--)
430          pGetExp(q,j) += pGetExp(m,j);
431        syzSetm(q);
432        lastmon = q;
433        i = pGetOrder(p);
434      }
435      else
436      {
437        spMemcpy(q,lastmon);
438        pSetComp(q,pGetComp(p));
439      }
440      pSetCoeff0(q,nMult(pGetCoeff(p),pGetCoeff(m)));
441      pIter(p);
442      if (p!=NULL)
443      {
444        pNext(q) = pNew();
445        pIter(q);
446      }
447      else break;
448    }
449  }
450  pNext(q) = NULL;
451  return result;
452}
453
454static inline poly syAdd(poly m1,poly m2)
455{
456  if (m1==NULL)
457    return m2;
458  if (m2==NULL)
459    return m1;
460  if ((pGetOrder(m1)>0) && (pGetOrder(m2)>0))
461  {
462    return pAdd(m1,m2);
463  }
464  else
465  {
466    poly p,result=p=pNew();
467    number n;
468    loop
469    {
470      if (pGetOrder(m1)>pGetOrder(m2))
471      {
472        pNext(p) = m1;
473        pIter(p);
474        pIter(m1);
475      }
476      else if (pGetOrder(m1)<pGetOrder(m2))
477      {
478        pNext(p) = m2;
479        pIter(p);
480        pIter(m2);
481      }
482      else if (currcomponents[m1->exp[0]]>currcomponents[m2->exp[0]])
483      {
484        pNext(p) = m1;
485        pIter(p);
486        pIter(m1);
487      }
488      else if (currcomponents[m1->exp[0]]<currcomponents[m2->exp[0]])
489      {
490        pNext(p) = m2;
491        pIter(p);
492        pIter(m2);
493      }
494      else
495      {
496        n = nAdd(pGetCoeff(m1),pGetCoeff(m2));
497        if (nIsZero(n))
498        {
499          pDelete1(&m1);
500        }
501        else
502        {
503          pNext(p) = m1;
504          pIter(p);
505          pIter(m1);
506          pSetCoeff(p,n);
507        }
508        pDelete1(&m2);
509      }
510      if (m1==NULL)
511      {
512        pNext(p) = m2;
513        break;
514      }
515      if (m2==NULL)
516      {
517        pNext(p) = m1;
518        break;
519      }
520    }
521    pDelete1(&result);
522    return result;
523  }
524}
525
526poly syOrdPolySchreyer(poly p)
527{
528  return pOrdPolySchreyer(p);
529}
530
531static poly sySPAdd(poly m1,poly m2,poly m)
532{
533  int i2=-INT_MAX;
534  int i1;
535  int j;
536  //register poly m1=m11;
537  poly p=redpol;
538  poly tm2=pNew();
539  number multwith=pGetCoeff(m);
540
541  pDelete1(&m1);
542  i1 = pGetOrder(m1);
543  pIter(m2);
544  spMemcpy(tm2,m2);
545  j = 1;
546  pGetExp(tm2,1)+=pGetExp(m,1);
547  int *ip=binomials+pGetExp(tm2,1);
548  loop
549  {
550    i2 += (*ip);
551    if (j==pVariables) break;
552    j++;
553    pGetExp(tm2,j)+=pGetExp(m,j);
554    ip+=highdeg_1+pGetExp(tm2,j);
555  }
556  pGetOrder(tm2) = i2;
557  pSetCoeff0(tm2,nMult(pGetCoeff(m2),multwith));
558  loop
559  {
560    j = i1-i2;
561    if (j>0/*m1->order>tm2->order*/)
562    {
563      p = pNext(p) = m1;
564      pIter(m1);
565      if (m1!=NULL)
566        i1 = pGetOrder(m1);
567      else
568        break;
569    }
570    else if (j<0/*m1->order<tm2->order*/)
571    {
572      p = pNext(p) = tm2;
573      j = pGetOrder(m2);
574      pIter(m2);
575      if (m2!=NULL)
576      {
577        j -=pGetOrder(m2);
578        tm2 = pNew();
579        if (j!=0)
580        {
581          spMemcpy(tm2,m2);
582          j = 1;
583          pGetExp(tm2,1)+=pGetExp(m,1);
584          int *ip=binomials+pGetExp(tm2,1);
585          i2=-INT_MAX;
586          loop
587          {
588            i2 += (*ip);
589            if (j==pVariables) break;
590            j++;
591            pGetExp(tm2,j)+=pGetExp(m,j);
592            ip+=highdeg_1+pGetExp(tm2,j);
593          }
594          pGetOrder(tm2) = i2;
595          //simple++;
596        }
597        else
598        {
599          spMemcpy(tm2,p);
600          pSetComp(tm2,pGetComp(m2));
601          //dsim++;
602        }
603        //pNext(tm2) = NULL;
604        pSetCoeff0(tm2,nMult(pGetCoeff(m2),multwith));
605      }
606      else
607        break;
608    }
609    else 
610    {
611      j = currcomponents[m1->exp[0]]-currcomponents[m2->exp[0]];
612      if (j>0/*currcomponents[m1->exp[0]]>currcomponents[m2->exp[0]]*/)
613      {
614        p = pNext(p) = m1;
615        pIter(m1);
616        if (m1!=NULL)
617          i1 = pGetOrder(m1);
618        else
619          break;
620      }
621      else 
622      {
623        if (j<0/*currcomponents[m1->exp[0]]<currcomponents[m2->exp[0]]*/)
624        {
625          p = pNext(p) = tm2;
626          j = pGetOrder(m2);
627        }
628        else
629        {
630          number n = nAdd(pGetCoeff(m1),pGetCoeff(tm2));
631          if (nIsZero(n))
632          {
633            pDelete1(&tm2);
634            nDelete(&n);
635            j = 0;
636          }
637          else
638          {
639            p = pNext(p) = tm2;
640            pSetCoeff(p,n);
641            j = pGetOrder(m2);
642          }
643          pDelete1(&m1);
644          if (m1==NULL)
645          {
646            pIter(m2);
647            break;
648          }
649          else
650            i1 = pGetOrder(m1);
651        }
652        pIter(m2);
653        if (m2!=NULL)
654        {
655          j -=pGetOrder(m2);
656          tm2 = pNew();
657          if (j!=0)
658          {
659            spMemcpy(tm2,m2);
660            j = 1;
661            pGetExp(tm2,1)+=pGetExp(m,1);
662            int *ip=binomials+pGetExp(tm2,1);
663            i2=-INT_MAX;
664            loop
665            {
666              i2 += (*ip);
667              if (j==pVariables) break;
668              j++;
669              pGetExp(tm2,j)+=pGetExp(m,j);
670              ip+=highdeg_1+pGetExp(tm2,j);
671            }
672            pGetOrder(tm2) = i2;
673            //simple++;
674          }
675          else
676          {
677            spMemcpy(tm2,p);
678            pSetComp(tm2,pGetComp(m2));
679            //dsim++;
680          }
681          //pNext(tm2) = NULL;
682          pSetCoeff0(tm2,nMult(pGetCoeff(m2),multwith));
683        }
684        else
685          break;
686      }
687    }
688  }
689  if (m2==NULL)
690  {
691    pNext(p) = m1;
692  }
693  else
694  {
695    pNext(p) = syMultT1(m2,m);
696  }
697  return pNext(redpol);
698}
699
700static inline poly sySPoly(poly m1,poly m2,poly lcm)
701{
702  poly a=pDivide(lcm,m1);
703  poly b1=NULL,b2=NULL;
704  if (pNext(m1)!=NULL)
705    b1 = syMultT(pNext(m1),a);
706  pDelete(&a);
707  a=pDivide(lcm,m2);
708  if (pNext(m2)!=NULL)
709    b2 = syMultTNeg(pNext(m2),a);
710  pDelete(&a);
711  return syAdd(b1,b2);
712}
713
714static poly sySPolyRed(poly m1,poly m2)
715{
716  if (pNext(m2)==NULL)
717  {
718    pDelete1(&m1);
719    return m1;
720  }
721  poly a=pDivide(m1,m2);
722  pSetCoeff0(a,nCopy(pGetCoeff(m1)));
723  pGetCoeff(a) = nNeg(pGetCoeff(a));
724  //return spSpolyRed(m2,m1,NULL);
725  if (pNext(m2)==NULL)
726  {
727    pDelete1(&m1);
728    return m1;
729  }
730  if (pNext(m1)==NULL)
731    return syMultT1(pNext(m2),a);
732  poly res;
733  if (pGetOrder(m1)>0)
734  {
735    res = spSpolyRed(m2,m1,NULL);
736  }
737  else
738  {
739    res = sySPAdd(m1,m2,a);
740  }
741  pDelete(&a);
742  return res;
743}
744
745BOOLEAN syDivisibleBy(poly a, poly b)
746{
747  if (a->exp[0]==b->exp[0])
748  {
749    int i=pVariables-1;
750    short *e1=&(a->exp[1]);
751    short *e2=&(b->exp[1]);
752    if ((*e1) > (*e2)) return FALSE;
753    do
754    {
755      i--;
756      e1++;
757      e2++;
758      if ((*e1) > (*e2)) return FALSE;
759    } while (i>0);
760    return TRUE;
761  }
762  else 
763    return FALSE;
764}
765
766BOOLEAN syDivisibleBy1(poly a, poly b)
767{
768  int i=pVariables-1;
769  short *e1=&(a->exp[1]);
770  short *e2=&(b->exp[1]);
771  if ((*e1) > (*e2)) return FALSE;
772  do
773  {
774    i--;
775    e1++;
776    e2++;
777    if ((*e1) > (*e2)) return FALSE;
778  } while (i>0);
779  return TRUE;
780}
781/*
782*int syDivisibleBy2(poly a, poly b)
783*{
784*  int i=pVariables-1;
785*  short *e1=&(a->exp[1]);
786*  short *e2=&(b->exp[1]);
787*  int result=(*e2)-(*e1);
788*  do
789*  {
790*    i--;
791*    e1++;
792*    e2++;
793*  } while (i>0);
794*  return 0;
795*}
796*/
797
798static int syLength(poly p)
799{
800  int result=0;
801
802  while (p!=NULL)
803  {
804    result++;
805    pIter(p);
806  }
807  return result;
808}
809
810poly syRedtail (poly p, ideal redWith, int index)
811{
812  poly h, hn;
813  int j,pos;
814
815  h = p;
816  hn = pNext(h);
817  while(hn != NULL)
818  {
819    j = Firstelem[index-1][pGetComp(hn)]-1;
820    if (j>=0)
821    {
822      pos = j+Howmuch[index-1][pGetComp(hn)];
823      while (j < pos)
824      {
825        if (syDivisibleBy1(redWith->m[j], hn))
826        {
827           //int syL=syLength(redWith->m[j]);
828            //Print("r");
829//for(int jj=j+1;jj<pos;jj++)
830//{
831  //if (syDivisibleBy1(redWith->m[jj],hn))
832  //{
833    //int syL1=syLength(redWith->m[jj]);
834    //if (syL1<syL)
835    //{
836      //j = jj;
837      //syL = syL1;
838      //Print("take poly %d with %d monomials\n",j,syL);
839      //Print("have poly %d with %d monomials\n",jj,syL1);
840    //}
841  //}
842//}
843          //if (pGetComp(redWith->m[j])!=pGetComp(hn))
844          //{
845            //Print("Hilfe!!!!!!!\n");
846            //Print("Fehler in Modul %d bei Elem %d\n",index,j);
847            //Print("Poly p: ");pWrite(hn);
848            //Print("Poly redWith: ");pWrite(redWith->m[j]);
849          //}
850          hn = sySPolyRed(hn,redWith->m[j]);
851          if (hn == NULL)
852          {
853            pNext(h) = NULL;
854            return p;
855          }
856          j = Firstelem[index-1][pGetComp(hn)]-1;
857          pos = j+Howmuch[index-1][pGetComp(hn)];
858        }
859        else
860        {
861          j++;
862        }
863      }
864    }
865    h = pNext(h) = hn;
866    hn = pNext(h);
867  }
868  return p;
869}
870
871/*3
872* initialize the resolution and puts in the argument as
873* zeroth entre, length must be > 0
874* assumes that the basering is degree-compatible
875*/
876static SRes syInitRes(ideal arg,int * length, intvec * Tl)
877{
878  if (idIs0(arg)) return NULL;
879  SRes resPairs = (SRes)Alloc0(*length*sizeof(SSet));
880  resPairs[0] = (SSet)Alloc0(IDELEMS(arg)*sizeof(SObject));
881  intvec * iv = idSort(arg);
882  int i;
883
884  for (i=0;i<IDELEMS(arg);i++)
885  {
886    (resPairs[0])[i].syz = pCopy(arg->m[(*iv)[i]-1]);
887  } 
888  delete iv;
889  (*Tl)[0] = IDELEMS(arg);
890  return resPairs;
891}
892
893/*3
894* determines the place of a polynomial in the right ordered resolution
895* set the vectors of truecomponents
896*/
897static void syOrder(poly p,resolvente res,resolvente orderedRes,int index,
898                    int realcomp)
899{
900  int i=IDELEMS(res[index-1])+1,j=0,k,tc,orc,ie=realcomp-1;
901  int *trind1=truecomponents[index-1],*trind=truecomponents[index];
902  int *bc=backcomponents[index],*F1=Firstelem[index-1];
903  int *H1=Howmuch[index-1];
904  poly pp;
905  polyset or=orderedRes[index]->m,or1=orderedRes[index-1]->m;
906
907  if (p==NULL) return;
908  if (realcomp==0) realcomp=1;
909 
910  if (index>1)
911    tc = trind1[pGetComp(p)]-1;
912  else
913    tc = pGetComp(p);
914  loop         //while ((j<ie) && (trind1[orc]<=tc+1))
915  {
916    if (j>=ie)
917      break;
918    else
919    {
920      orc = pGetComp(or[j]);
921      if (trind1[orc]>tc+1) break;
922      j += H1[orc];
923    }
924  }
925  if (j>ie)
926  {
927    WerrorS("orderedRes to small");
928    return;
929  }
930  ie++;
931  if (or[j]!=NULL)
932  {
933    for (k=ie-1;k>j;k--)
934    {
935      or[k] = or[k-1];
936      bc[k] = bc[k-1];
937    }
938  }
939  or[j] = p;
940  bc[j] = realcomp-1;
941  (H1[pGetComp(p)])++;
942  for (k=0;k<i;k++)
943  {
944    if (F1[k]>j)
945      (F1[k])++;
946  }
947  if (F1[pGetComp(p)]==0)
948    F1[pGetComp(p)]=j+1;
949//Print("write in sort %d till %d\n",index-1,i-1);
950//Print("poly: ");pWrite(p);
951//Print("in module %d as %d -th element\n",index,j);
952  for (k=0;k<IDELEMS(res[index]);k++)
953  {
954    if (trind[k]>j)
955      trind[k] += 1;
956  }
957  for (k=IDELEMS(res[index])-1;k>realcomp;k--)
958    trind[k] = trind[k-1];
959  trind[realcomp] = j+1;
960}
961
962/*3
963* reduces all pairs of degree deg in the module index
964* put the reduced generators to the resolvente which contains
965* the truncated std
966*/
967static void syRedPairsOfCurrDeg(SRes resPairs, resolvente res, resolvente orderedRes,
968               int deg, int index, intvec * Tl, intvec ** minGen)
969{
970  int i=0,j,k=IDELEMS(res[index]),kk;
971  number coefgcd,n;
972  poly p=NULL;
973
974  if (resPairs[index]==NULL) return;
975  while ((k>0) && (res[index]->m[k-1]==NULL)) k--;
976  while ((i<(*Tl)[index]) && ((resPairs[index])[i].lcm!=NULL)
977          && (pGetOrder((resPairs[index])[i].lcm)<deg)) i++;
978  if ((i>=(*Tl)[index]) || ((resPairs[index])[i].lcm==NULL)
979          || (pGetOrder((resPairs[index])[i].lcm)>deg)) return;
980  while ((i<(*Tl)[index]) && ((resPairs[index])[i].lcm!=NULL)
981          && (pGetOrder((resPairs[index])[i].lcm)==deg)) 
982  {
983    (resPairs[index])[i].p = 
984      spSpolyCreate((resPairs[index])[i].p2, (resPairs[index])[i].p1,NULL);
985    coefgcd = 
986      nGcd(pGetCoeff((resPairs[index])[i].p1),pGetCoeff((resPairs[index])[i].p1));
987    if ((*minGen[index])[(resPairs[index])[i].ind2]==0)
988    {
989      (resPairs[index])[i].syz = pHead((resPairs[index])[i].lcm);
990      p = (resPairs[index])[i].syz;
991      pSetCoeff(p,nDiv(pGetCoeff((resPairs[index])[i].p1),coefgcd));
992      pGetCoeff(p) = nNeg(pGetCoeff(p));
993      pSetComp(p,(resPairs[index])[i].ind2+1);
994    }
995    if ((*minGen[index])[(resPairs[index])[i].ind1]==0)
996    {
997      if ((resPairs[index])[i].syz==NULL)
998      {
999        (resPairs[index])[i].syz = pHead((resPairs[index])[i].lcm);
1000        p = (resPairs[index])[i].syz;
1001      }
1002      else
1003      {
1004        pNext(p) = pHead((resPairs[index])[i].lcm);
1005        pIter(p);
1006      }
1007      pSetComp(p,(resPairs[index])[i].ind1+1);
1008      pSetCoeff(p,nDiv(pGetCoeff((resPairs[index])[i].p2),coefgcd));
1009    }
1010    nDelete(&coefgcd);
1011    j = k-1;
1012    while ((j>=0) && (res[index]->m[j]!=NULL) 
1013           && ((resPairs[index])[i].p!=NULL))
1014    {
1015      if (syDivisibleBy(res[index]->m[j],(resPairs[index])[i].p))
1016      {
1017        if ((*minGen[index])[j]==0)
1018        {
1019          if ((resPairs[index])[i].syz==NULL)
1020          {
1021            (resPairs[index])[i].syz = pHead((resPairs[index])[i].p);
1022            p = (resPairs[index])[i].syz;
1023          }
1024          else
1025          {
1026            pNext(p) = pHead((resPairs[index])[i].p);
1027            pIter(p);
1028          }
1029          pSetComp(p,j+1);
1030          pGetCoeff(p) = nNeg(pGetCoeff(p));
1031        }
1032        (resPairs[index])[i].p = 
1033          spSpolyRed(res[index]->m[j],(resPairs[index])[i].p,NULL);
1034        j = k-1;
1035      }
1036      else
1037      {
1038        j--;
1039      }
1040    }
1041    if ((resPairs[index])[i].p != NULL)
1042    {
1043      if (TEST_OPT_PROT) Print("g");
1044      if (k==IDELEMS(res[index]))
1045      {
1046        pEnlargeSet(&(res[index]->m),IDELEMS(res[index]),16);
1047        truecomponents[index]=(int*)ReAlloc((ADDRESS)truecomponents[index],
1048                                       (IDELEMS(res[index])+1)*sizeof(int),
1049                                       (IDELEMS(res[index])+17)*sizeof(int));
1050        backcomponents[index]=(int*)ReAlloc((ADDRESS)backcomponents[index],
1051                                       (IDELEMS(res[index])+1)*sizeof(int),
1052                                       (IDELEMS(res[index])+17)*sizeof(int));
1053        Howmuch[index]=(int*)ReAlloc((ADDRESS)Howmuch[index],
1054                                       (IDELEMS(res[index])+1)*sizeof(int),
1055                                       (IDELEMS(res[index])+17)*sizeof(int));
1056        Firstelem[index]=(int*)ReAlloc((ADDRESS)Firstelem[index],
1057                                       (IDELEMS(res[index])+1)*sizeof(int),
1058                                       (IDELEMS(res[index])+17)*sizeof(int));
1059        int mk;
1060        for (mk=IDELEMS(res[index])+1;mk<IDELEMS(res[index])+17;mk++)
1061        {
1062          Howmuch[index][mk] = 0;
1063          Firstelem[index][mk] = 0;
1064        }
1065//Print("sort %d has now size %d\n",index,IDELEMS(res[index])+17);
1066        IDELEMS(res[index]) += 16;
1067        pEnlargeSet(&(orderedRes[index]->m),IDELEMS(orderedRes[index]),16);
1068        IDELEMS(orderedRes[index]) += 16;
1069        intvec * temp = new intvec(minGen[index]->length()+16);
1070        for (kk=minGen[index]->length()-1;kk>=0;kk--)
1071          (*temp)[kk] = (*minGen[index])[kk];
1072        delete minGen[index];
1073        minGen[index] = temp;
1074      }
1075      if ((resPairs[index])[i].syz==NULL)
1076      {
1077        (resPairs[index])[i].syz = pHead((resPairs[index])[i].p);
1078        p = (resPairs[index])[i].syz;
1079      }
1080      else
1081      {
1082        pNext(p) = pHead(resPairs[index][i].p);
1083        pIter(p);
1084      }
1085      if (p!=NULL)
1086      {
1087        k++;
1088        pSetComp(p,k);
1089        pGetCoeff(p) = nNeg(pGetCoeff(p));
1090      }
1091      if (resPairs[index][i].p!=NULL)
1092      {
1093        res[index]->m[k] = resPairs[index][i].p;
1094        pNorm(res[index]->m[k]);
1095        syOrder(res[index]->m[k-1],res,orderedRes,index,k);
1096      }
1097      resPairs[index][i].isNotMinimal = 1;
1098      resPairs[index][i].p =NULL;
1099    }
1100    else
1101      if (TEST_OPT_PROT) Print(".");
1102    p = NULL;
1103    i++;
1104  }
1105} 
1106
1107/*3
1108* reduces all pairs of degree deg in the module index
1109* put the reduced generators to the resolvente which contains
1110* the truncated std
1111*/
1112static void syRedNextPairs(SSet nextPairs, resolvente res, resolvente orderedRes,
1113               int howmuch, int index)
1114{
1115  int i=howmuch-1,j,k=IDELEMS(res[index]),ks=IDELEMS(res[index+1]),kk,l,ll;
1116  number coefgcd,n;
1117  polyset redset=orderedRes[index]->m;
1118  poly p=NULL,q;
1119  BOOLEAN isDivisible;
1120
1121  if ((nextPairs==NULL) || (howmuch==0)) return;
1122  while ((k>0) && (res[index]->m[k-1]==NULL)) k--;
1123  while ((ks>0) && (res[index+1]->m[ks-1]==NULL)) ks--;
1124  while (i>=0)
1125  {
1126    isDivisible = FALSE;
1127    if (res[index+1]!=NULL)
1128    {
1129      l = Firstelem[index][pGetComp(nextPairs[i].lcm)]-1;
1130      if (l>=0)
1131      {
1132        ll = l+Howmuch[index][pGetComp(nextPairs[i].lcm)];
1133        while ((l<ll) && (!isDivisible))
1134        {
1135          if (res[index+1]->m[l]!=NULL)
1136          {
1137            isDivisible = isDivisible || 
1138              syDivisibleBy1(orderedRes[index+1]->m[l],nextPairs[i].lcm);
1139          }
1140          l++;
1141        }
1142      }
1143    }
1144    if (!isDivisible)
1145    {
1146      nextPairs[i].p = 
1147        sySPoly(nextPairs[i].p1, nextPairs[i].p2,nextPairs[i].lcm);
1148      coefgcd = 
1149        nGcd(pGetCoeff(nextPairs[i].p1),pGetCoeff(nextPairs[i].p1));
1150      nextPairs[i].syz = pHead(nextPairs[i].lcm);
1151      pSetm(nextPairs[i].syz);
1152      p = nextPairs[i].syz;
1153      pSetCoeff(p,nDiv(pGetCoeff(nextPairs[i].p1),coefgcd));
1154      pGetCoeff(p) = nNeg(pGetCoeff(p));
1155      pSetComp(p,nextPairs[i].ind2+1);
1156      pNext(p) = pHead(nextPairs[i].lcm);
1157      pIter(p);
1158      pSetm(p);
1159      pSetComp(p,nextPairs[i].ind1+1);
1160      pSetCoeff(p,nDiv(pGetCoeff(nextPairs[i].p2),coefgcd));
1161      nDelete(&coefgcd);
1162      if (nextPairs[i].p != NULL)
1163      {
1164        q = nextPairs[i].p;
1165        j = Firstelem[index-1][pGetComp(q)]-1;
1166        int pos = j+Howmuch[index-1][pGetComp(q)];
1167        loop
1168        {
1169          if (j<0) break;
1170          if (syDivisibleBy1(redset[j],q))
1171          {
1172            //int syL=syLength(redset[j]);
1173            //Print("r");
1174//for(int jj=j+1;jj<k;jj++)
1175//{
1176  //if (syDivisibleBy(redset[jj],q))
1177  //{
1178    //int syL1=syLength(redset[jj]);
1179    //if (syL1<syL)
1180    //{
1181      //j = jj;
1182      //syL = syL1;
1183      //Print("take poly %d with %d monomials\n",j,syL);
1184      //Print("have poly %d with %d monomials\n",jj,syL1);
1185    //}
1186  //}
1187//}
1188            pNext(p) = pHead(q);
1189            pIter(p);
1190            pSetComp(p,backcomponents[index][j]+1);
1191            pGetCoeff(p) = nNeg(pGetCoeff(p));
1192            q = sySPolyRed(q,redset[j]);
1193            if (q==NULL) break;
1194            j = Firstelem[index-1][pGetComp(q)]-1;
1195            pos = j+Howmuch[index-1][pGetComp(q)];
1196          }
1197          else
1198          {
1199            j++;
1200            if (j==pos) break;
1201          }
1202        }
1203        nextPairs[i].p = q;
1204      }
1205      if (nextPairs[i].p != NULL)
1206      {
1207        if (TEST_OPT_PROT) Print("g");
1208        if (k==IDELEMS(res[index]))
1209        {
1210          pEnlargeSet(&(res[index]->m),IDELEMS(res[index]),16);
1211          truecomponents[index]=(int*)ReAlloc((ADDRESS)truecomponents[index],
1212                                       (IDELEMS(res[index])+1)*sizeof(int),
1213                                       (IDELEMS(res[index])+17)*sizeof(int));
1214          backcomponents[index]=(int*)ReAlloc((ADDRESS)backcomponents[index],
1215                                       (IDELEMS(res[index])+1)*sizeof(int),
1216                                       (IDELEMS(res[index])+17)*sizeof(int));
1217          Howmuch[index]=(int*)ReAlloc((ADDRESS)Howmuch[index],
1218                                       (IDELEMS(res[index])+1)*sizeof(int),
1219                                       (IDELEMS(res[index])+17)*sizeof(int));
1220          Firstelem[index]=(int*)ReAlloc((ADDRESS)Firstelem[index],
1221                                       (IDELEMS(res[index])+1)*sizeof(int),
1222                                       (IDELEMS(res[index])+17)*sizeof(int));
1223        int mk;
1224        for (mk=IDELEMS(res[index])+1;mk<IDELEMS(res[index])+17;mk++)
1225        {
1226          Howmuch[index][mk] = 0;
1227          Firstelem[index][mk] = 0;
1228        }
1229//Print("sort %d has now size %d\n",index,IDELEMS(res[index])+17);
1230          IDELEMS(res[index]) += 16;
1231          pEnlargeSet(&(orderedRes[index]->m),IDELEMS(orderedRes[index]),16);
1232          IDELEMS(orderedRes[index]) += 16;
1233          redset=orderedRes[index]->m;
1234        }
1235        pNext(p) = pHead(nextPairs[i].p);
1236        pIter(p);
1237        if (p!=NULL)
1238        {
1239          k++;
1240          pSetComp(p,k);
1241          pGetCoeff(p) = nNeg(pGetCoeff(p));
1242        }
1243        if (nextPairs[i].p!=NULL)
1244        {
1245          res[index]->m[k-1] = nextPairs[i].p;
1246          pNorm(res[index]->m[k-1]);
1247          syOrder(res[index]->m[k-1],res,orderedRes,index,k);
1248        }
1249        nextPairs[i].isNotMinimal = 1;
1250        nextPairs[i].p =NULL;
1251      }
1252      else
1253      {
1254        if (TEST_OPT_PROT) Print(".");
1255        if (index % 2==0)
1256          euler++;
1257        else 
1258          euler--;
1259      }
1260      if (ks==IDELEMS(res[index+1]))
1261      {
1262        pEnlargeSet(&(res[index+1]->m),IDELEMS(res[index+1]),16);
1263        truecomponents[index+1]=(int*)ReAlloc((ADDRESS)truecomponents[index+1],
1264                                       (IDELEMS(res[index+1])+1)*sizeof(int),
1265                                       (IDELEMS(res[index+1])+17)*sizeof(int));
1266        backcomponents[index+1]=(int*)ReAlloc((ADDRESS)backcomponents[index+1],
1267                                       (IDELEMS(res[index+1])+1)*sizeof(int),
1268                                       (IDELEMS(res[index+1])+17)*sizeof(int));
1269        Howmuch[index+1]=(int*)ReAlloc((ADDRESS)Howmuch[index+1],
1270                                       (IDELEMS(res[index+1])+1)*sizeof(int),
1271                                       (IDELEMS(res[index+1])+17)*sizeof(int));
1272        Firstelem[index+1]=(int*)ReAlloc((ADDRESS)Firstelem[index+1],
1273                                       (IDELEMS(res[index+1])+1)*sizeof(int),
1274                                       (IDELEMS(res[index+1])+17)*sizeof(int));
1275        int mk;
1276        for (mk=IDELEMS(res[index+1])+1;mk<IDELEMS(res[index+1])+17;mk++)
1277        {
1278          Howmuch[index+1][mk] = 0;
1279          Firstelem[index+1][mk] = 0;
1280        }
1281//Print("sort %d has now size %d\n",index+1,IDELEMS(res[index+1])+17);
1282        IDELEMS(res[index+1]) += 16;
1283        pEnlargeSet(&(orderedRes[index+1]->m),IDELEMS(orderedRes[index+1]),16);
1284        IDELEMS(orderedRes[index+1]) += 16;
1285      }
1286      currcomponents = truecomponents[index];
1287      res[index+1]->m[ks] = syRedtail(nextPairs[i].syz,orderedRes[index+1],index+1);
1288      currcomponents = truecomponents[index-1];
1289      //res[index+1]->m[ks] = nextPairs[i].syz;
1290      pNorm(res[index+1]->m[ks]);
1291      nextPairs[i].syz =NULL;
1292      syOrder(res[index+1]->m[ks],res,orderedRes,index+1,ks+1);
1293      ks++;
1294      p = NULL;
1295    }
1296    else
1297    {
1298      syDeletePair(&nextPairs[i]);
1299      //crit++;
1300    }
1301    i--;
1302  }
1303} 
1304
1305/*3
1306* reduces the generators of the module index in degree deg
1307* (which are actual syzygies of the module index-1)
1308* wrt. the ideal generated by elements of lower degrees
1309*/
1310static void syRedGenerOfCurrDeg(SRes resPairs,resolvente res,resolvente orderedRes,
1311         int deg, int index,intvec * Tl)
1312{
1313  int i=0,j,k=IDELEMS(res[index]),kk;
1314
1315  while ((k>0) && (res[index]->m[k-1]==NULL)) k--;
1316  while ((i<(*Tl)[index-1]) && (((resPairs[index-1])[i].syz==NULL) || 
1317          (pTotaldegree((resPairs[index-1])[i].syz)<deg)))
1318    i++;
1319  if ((i>=(*Tl)[index-1]) || (pTotaldegree((resPairs[index-1])[i].syz)>deg)) return;
1320  while ((i<(*Tl)[index-1]) && (((resPairs[index-1])[i].syz==NULL) ||
1321         (pTotaldegree((resPairs[index-1])[i].syz)==deg))) 
1322  {
1323    if ((resPairs[index-1])[i].syz!=NULL)
1324    {
1325      j = k-1;
1326      while ((j>=0) && (res[index]->m[j]!=NULL) && 
1327             ((resPairs[index-1])[i].syz!=NULL))
1328      {
1329        if (syDivisibleBy(res[index]->m[j],(resPairs[index-1])[i].syz))
1330        {
1331          //Print("r");
1332          (resPairs[index-1])[i].syz = 
1333            //spSpolyRed(res[index]->m[j],(resPairs[index-1])[i].syz,NULL);
1334            sySPolyRed((resPairs[index-1])[i].syz,res[index]->m[j]);
1335          j = k-1;
1336        }
1337        else
1338        {
1339          j--;
1340        }
1341      }
1342      if ((resPairs[index-1])[i].syz != NULL)
1343      {
1344        if (k==IDELEMS(res[index]))
1345        {
1346          pEnlargeSet(&(res[index]->m),IDELEMS(res[index]),16);
1347          truecomponents[index]=(int*)ReAlloc((ADDRESS)truecomponents[index],
1348                                       (IDELEMS(res[index])+1)*sizeof(int),
1349                                       (IDELEMS(res[index])+17)*sizeof(int));
1350          backcomponents[index]=(int*)ReAlloc((ADDRESS)backcomponents[index],
1351                                       (IDELEMS(res[index])+1)*sizeof(int),
1352                                       (IDELEMS(res[index])+17)*sizeof(int));
1353          Howmuch[index]=(int*)ReAlloc((ADDRESS)Howmuch[index],
1354                                       (IDELEMS(res[index])+1)*sizeof(int),
1355                                       (IDELEMS(res[index])+17)*sizeof(int));
1356          Firstelem[index]=(int*)ReAlloc((ADDRESS)Firstelem[index],
1357                                       (IDELEMS(res[index])+1)*sizeof(int),
1358                                       (IDELEMS(res[index])+17)*sizeof(int));
1359        int mk;
1360        for (mk=IDELEMS(res[index])+1;mk<IDELEMS(res[index])+17;mk++)
1361        {
1362          Howmuch[index][mk] = 0;
1363          Firstelem[index][mk] = 0;
1364        }
1365//Print("sort %d has now size %d\n",index,IDELEMS(res[index])+17);
1366          IDELEMS(res[index]) += 16;
1367          pEnlargeSet(&(orderedRes[index]->m),IDELEMS(orderedRes[index]),16);
1368          IDELEMS(orderedRes[index]) += 16;
1369        }
1370        if (BTEST1(6))
1371        {
1372          if ((resPairs[index-1])[i].isNotMinimal==0)
1373          {
1374            PrintLn();
1375            Print("minimal generator: ");pWrite((resPairs[index-1])[i].syz);
1376            Print("comes from: ");pWrite((resPairs[index-1])[i].p1);
1377            Print("and: ");pWrite((resPairs[index-1])[i].p2);
1378          }
1379        }
1380        //res[index]->m[k] = (resPairs[index-1])[i].syz;
1381        res[index]->m[k] = syRedtail((resPairs[index-1])[i].syz,orderedRes[index],index);
1382        pNorm(res[index]->m[k]);
1383  //      (resPairs[index-1])[i].syz = NULL;
1384        k++;
1385        syOrder(res[index]->m[k-1],res,orderedRes,index,k);
1386        euler++;
1387      }
1388      //else
1389      //{
1390        //zeroRed++;
1391      //}
1392    }
1393    i++;
1394  }
1395}
1396
1397/*3
1398* puts a pair into the right place in resPairs
1399*/
1400static void syEnterPair(SSet sPairs, SObject * so, int * sPlength,int index)
1401{
1402  int ll,k,no=pGetOrder((*so).lcm),sP=*sPlength,i;
1403  poly p=(*so).lcm;
1404
1405  if ((sP==0) || (sPairs[sP-1].order<=no))
1406    ll = sP;
1407  else if (sP==1)
1408    ll = 0;
1409  else
1410  {
1411    int an=0,en=sP-1;
1412    loop
1413    {
1414      if (an>=en-1)
1415      {
1416        if ((sPairs[an].order<=no) && (sPairs[an+1].order>no))
1417        {
1418          ll = an+1;
1419          break;
1420        }
1421        else if ((sPairs[en].order<=no) && (sPairs[en+1].order>no))
1422        {
1423          ll = en+1;
1424          break;
1425        }
1426        else if (sPairs[an].order>no)
1427        {
1428          ll = an;
1429          break;
1430        }
1431        else
1432        {
1433          Print("Hier ist was faul!\n");
1434          break;
1435        }
1436      }
1437      i=(an+en) / 2;
1438      if (sPairs[i].order <= no)
1439        an=i;
1440      else
1441        en=i;
1442    }
1443  }
1444  for (k=(*sPlength);k>ll;k--)
1445  {
1446    syCopyPair(&sPairs[k-1],&sPairs[k]);
1447  }
1448  syCopyPair(so,&sPairs[ll]);
1449  (*sPlength)++;
1450}
1451
1452/*3
1453* applies GM-criteria to triples (i,j,k) where j is fixed
1454*/
1455static void syCrit(SRes resPairs, intvec * Tl, int index, int j)
1456{
1457  int i,k,l,ll,ini=max(1,(*Tl)[index]);
1458  ideal temp=idInit(ini,1);
1459   
1460  for (l=(*Tl)[index]-1;l>=0;l--)
1461  {
1462    if ((resPairs[index][l].lcm!=NULL) && 
1463        (pGetComp(resPairs[index][l].lcm)==j))
1464    {
1465      temp->m[l] = pDivide(resPairs[index][l].p1,resPairs[index][l].lcm);
1466      pSetComp(temp->m[l],resPairs[index][l].ind1);
1467    }
1468  }
1469  idSkipZeroes(temp);
1470  k = IDELEMS(temp)-1;
1471  while ((k>=0) && (temp->m[k]==NULL)) k--;
1472  while (k>0)
1473  {
1474    for (i=0;i<k;i++)
1475    {
1476      l = pVariables;
1477      while ((l>0) && (pGetExp(temp->m[i],l)*pGetExp(temp->m[k],l)==0))
1478        l--;
1479      if (l==0)
1480      {
1481        ll = 0;
1482        while ((ll<(*Tl)[index]) && 
1483          (((pGetComp(temp->m[i])!=resPairs[index][ll].ind1) ||
1484          (pGetComp(temp->m[k])!=resPairs[index][ll].ind2)) &&
1485          ((pGetComp(temp->m[k])!=resPairs[index][ll].ind1) ||
1486          (pGetComp(temp->m[i])!=resPairs[index][ll].ind2))))
1487          ll++;
1488        if ((ll<(*Tl)[index]) && (resPairs[index][ll].lcm!=NULL))
1489        {
1490          syDeletePair(&resPairs[index][ll]);
1491          //crit++;
1492        }
1493      }
1494    }
1495    k--;
1496  }
1497  syCompactifyPairSet(resPairs[index],(*Tl)[index],0);
1498  idDelete(&temp);
1499}
1500
1501/*3
1502* computes pairs from the new elements (beginning with the element newEl)
1503* in the module index
1504*/
1505static void syCreateNewPairs(SRes resPairs, intvec * Tl, resolvente res,
1506                int index, int newEl)
1507{
1508  SSet temp;
1509  SObject tso;
1510  int i,ii,j,k=IDELEMS(res[index]),l=(*Tl)[index],ll;
1511  int qc,first,pos,jj,j1;
1512  poly p,q;
1513  polyset rs=res[index]->m,nPm;
1514
1515  while ((k>0) && (res[index]->m[k-1]==NULL)) k--;
1516  if (newEl>=k) return;
1517  ideal nP=idInit(k,res[index]->rank);
1518  nPm=nP->m;
1519  while ((l>0) && ((resPairs[index])[l-1].p1==NULL)) l--;
1520//Print("new pairs in module %d\n",index);
1521  for (j=newEl;j<k;j++)
1522  {
1523    q = rs[j];
1524    qc = pGetComp(q);
1525    first = Firstelem[index-1][pGetComp(q)]-1;
1526    pos = first+Howmuch[index-1][pGetComp(q)];
1527    for (i=first;i<pos;i++)
1528    {
1529      jj = backcomponents[index][i];
1530      if (jj>=j) break;
1531      p = pOne();
1532      pLcm(rs[jj],q,p);
1533      pSetComp(p,j+1);
1534      ii = first;
1535      loop
1536      {
1537        j1 = backcomponents[index][ii];
1538        if (nPm[j1]!=NULL)
1539        {
1540          if (syDivisibleBy1(nPm[j1],p))
1541          {
1542            pDelete(&p);
1543            break;
1544          }
1545          else if (syDivisibleBy1(p,nPm[j1]))
1546          {
1547            pDelete(&(nPm[j1]));
1548            break;
1549          }
1550        }
1551        ii++;
1552        if (ii>=pos) break;
1553      }
1554      if (p!=NULL)
1555      {
1556        nPm[jj] = p;
1557      }
1558    }
1559    for (i=0;i<j;i++)
1560    {
1561      if (nPm[i]!=NULL)
1562      {
1563        if (l>=(*Tl)[index])
1564        {
1565          temp = (SSet)Alloc0(((*Tl)[index]+16)*sizeof(SObject));
1566          for (ll=0;ll<(*Tl)[index];ll++)
1567          {
1568            temp[ll].p = (resPairs[index])[ll].p;
1569            temp[ll].p1 = (resPairs[index])[ll].p1;
1570            temp[ll].p2 = (resPairs[index])[ll].p2;
1571            temp[ll].syz = (resPairs[index])[ll].syz;
1572            temp[ll].lcm = (resPairs[index])[ll].lcm;
1573            temp[ll].ind1 = (resPairs[index])[ll].ind1;
1574            temp[ll].ind2 = (resPairs[index])[ll].ind2;
1575            temp[ll].order = (resPairs[index])[ll].order;
1576            temp[ll].isNotMinimal = (resPairs[index])[ll].isNotMinimal;
1577          }
1578          Free((ADDRESS)resPairs[index],(*Tl)[index]*sizeof(SObject));
1579          (*Tl)[index] += 16;
1580          resPairs[index] = temp;
1581        }
1582        tso.lcm = p = nPm[i];
1583        nPm[i] = NULL;
1584        tso.order = pGetOrder(p) = pTotaldegree(p);
1585        tso.p1 = rs[i];
1586        tso.p2 = q;
1587        tso.ind1 = i;
1588        tso.ind2 = j;
1589        tso.isNotMinimal = 0;
1590        tso.p = NULL;
1591        tso.syz = NULL;
1592        syEnterPair(resPairs[index],&tso,&l,index);
1593      }
1594    }
1595  }
1596  idDelete(&nP);
1597}
1598
1599static void sySyzTail(SRes resPairs, intvec * Tl, resolvente orderedRes,
1600                resolvente res, int index, int newEl)
1601{
1602  int j,ll,k=IDELEMS(res[index]);
1603  while ((k>0) && (res[index]->m[k-1]==NULL)) k--;
1604  for (j=newEl;j<k;j++)
1605  {
1606    ll = 0;
1607    while ((ll<(*Tl)[index]) && (resPairs[index][ll].p1!=NULL) &&
1608           (resPairs[index][ll].ind1!=j) && (resPairs[index][ll].ind2!=j))
1609      ll++;
1610    if ((ll<(*Tl)[index]) && (resPairs[index][ll].p1!=NULL))
1611      res[index]->m[j] = syRedtail(res[index]->m[j],orderedRes[index],index);
1612  }
1613}
1614
1615/*3
1616* looks through the pair set and the given module for
1617* remaining pairs or generators to consider
1618*/
1619static BOOLEAN syCheckPairs(SRes resPairs,intvec * Tl,
1620                   int length,int * actdeg)
1621{
1622  int newdeg=*actdeg,i,index=0,t;
1623 
1624  while ((index<length) && (resPairs[index]!=NULL))
1625  {
1626    i = 0;
1627    while ((i<(*Tl)[index]))
1628    {
1629      t = *actdeg;
1630      if ((resPairs[index])[i].lcm!=NULL)
1631      {
1632        if (pGetOrder((resPairs[index])[i].lcm) > *actdeg)
1633          t = pGetOrder((resPairs[index])[i].lcm);
1634      }
1635      else if ((resPairs[index])[i].syz!=NULL)
1636      {
1637        if (pGetOrder((resPairs[index])[i].syz) > *actdeg)
1638          t = pGetOrder((resPairs[index])[i].syz);
1639      }
1640      if ((t>*actdeg) && ((newdeg==*actdeg) || (t<newdeg)))
1641      {
1642        newdeg = t;
1643        break;
1644      }
1645      i++;
1646    } 
1647    index++;
1648  }
1649  if (newdeg>*actdeg)
1650  {
1651    *actdeg = newdeg;
1652    return TRUE;
1653  }
1654  else return FALSE;
1655}
1656
1657/*3
1658* looks through the pair set and the given module for
1659* remaining pairs or generators to consider
1660* returns a pointer to the first pair and the number of them in the given module
1661* works with slanted degree (i.e. deg=realdeg-index)
1662*/
1663static SSet syChosePairs(SRes resPairs,intvec * Tl, int *index, int *howmuch,
1664                   int length,int * actdeg)
1665{
1666  int newdeg=*actdeg,newindex=-1,i,t,sldeg;
1667  SSet result;
1668 
1669  while (*index<length)
1670  {
1671    if (resPairs[*index]!=NULL)
1672    {
1673      sldeg = (*actdeg)+*index;
1674      i = 0;
1675      if (*index!=0)
1676      {
1677        while ((i<(*Tl)[*index]))
1678        {
1679          if ((resPairs[*index])[i].lcm!=NULL)
1680          {
1681            if (pGetOrder((resPairs[*index])[i].lcm) == sldeg)
1682            {
1683              result = &(resPairs[*index])[i];
1684              *howmuch =1;
1685              i++;
1686              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].lcm!=NULL)
1687                      && (pGetOrder((resPairs[*index])[i].lcm) == sldeg))
1688              {
1689                i++;
1690                (*howmuch)++;
1691              }
1692              return result;
1693            }
1694          }
1695          i++;
1696        }
1697      }
1698      else
1699      {
1700        while ((i<(*Tl)[*index]))
1701        {
1702          if ((resPairs[*index])[i].syz!=NULL)
1703          {
1704            if (pTotaldegree((resPairs[*index])[i].syz) == sldeg)
1705            {
1706              result = &(resPairs[*index])[i];
1707              (*howmuch) =1;
1708              i++;
1709              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].syz!=NULL)
1710                      && (pTotaldegree((resPairs[*index])[i].syz) == *actdeg))
1711              {
1712                i++;
1713                (*howmuch)++;
1714              }
1715              return result;
1716            }
1717          }
1718          i++;
1719        }
1720      }
1721    }
1722    (*index)++;
1723  }
1724  *index = 0;
1725  //if (TEST_OPT_PROT) Print("(Euler:%d)",euler);
1726  while (*index<length)
1727  {
1728    if (resPairs[*index]!=NULL)
1729    {
1730      i = 0;
1731      while ((i<(*Tl)[*index]))
1732      {
1733        t = *actdeg+*index;
1734        if ((resPairs[*index])[i].lcm!=NULL)
1735        {
1736          if (pGetOrder((resPairs[*index])[i].lcm) > *actdeg+*index)
1737            t = pGetOrder((resPairs[*index])[i].lcm);
1738        }
1739        else if ((resPairs[*index])[i].syz!=NULL)
1740        {
1741          if (pTotaldegree((resPairs[*index])[i].syz) > *actdeg+*index)
1742            t = pTotaldegree((resPairs[*index])[i].syz);
1743        }
1744        if ((t>*actdeg+*index) && ((newdeg==*actdeg) || (t<newdeg+*index)))
1745        {
1746          newdeg = t-*index;
1747          newindex = *index;
1748          break;
1749        }
1750        i++;
1751      } 
1752    }
1753    (*index)++;
1754  }
1755  if (newdeg>*actdeg)
1756  {
1757    *actdeg = newdeg;
1758    *index = newindex;
1759    return syChosePairs(resPairs,Tl,index,howmuch,length,actdeg);
1760  }
1761  else return NULL;
1762}
1763
1764/*3
1765* looks through the pair set and the given module for
1766* remaining pairs or generators to consider
1767* returns a pointer to the first pair and the number of them in the given module
1768* works deg by deg
1769*/
1770static SSet syChosePairs1(SRes resPairs,intvec * Tl, int *index, int *howmuch,
1771                   int length,int * actdeg)
1772{
1773  int newdeg=*actdeg,newindex=-1,i,t;
1774  SSet result;
1775 
1776  while (*index>=0)
1777  {
1778    if (resPairs[*index]!=NULL)
1779    {
1780      i = 0;
1781      if (*index!=0)
1782      {
1783        while ((i<(*Tl)[*index]))
1784        {
1785          if ((resPairs[*index])[i].lcm!=NULL)
1786          {
1787            if (pGetOrder((resPairs[*index])[i].lcm) == *actdeg)
1788            {
1789              result = &(resPairs[*index])[i];
1790              *howmuch =1;
1791              i++;
1792              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].lcm!=NULL)
1793                      && (pGetOrder((resPairs[*index])[i].lcm) == *actdeg))
1794              {
1795                i++;
1796                (*howmuch)++;
1797              }
1798              return result;
1799            }
1800          }
1801          i++;
1802        }
1803      }
1804      else
1805      {
1806        while ((i<(*Tl)[*index]))
1807        {
1808          if ((resPairs[*index])[i].syz!=NULL)
1809          {
1810            if (pTotaldegree((resPairs[*index])[i].syz) == *actdeg)
1811            {
1812              result = &(resPairs[*index])[i];
1813              (*howmuch) =1;
1814              i++;
1815              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].syz!=NULL)
1816                      && (pTotaldegree((resPairs[*index])[i].syz) == *actdeg))
1817              {
1818                i++;
1819                (*howmuch)++;
1820              }
1821              return result;
1822            }
1823          }
1824          i++;
1825        }
1826      }
1827    }
1828    (*index)--;
1829  }
1830  *index = length-1;
1831  while (*index>=0)
1832  {
1833    if (resPairs[*index]!=NULL)
1834    {
1835      i = 0;
1836      while ((i<(*Tl)[*index]))
1837      {
1838        t = *actdeg;
1839        if ((resPairs[*index])[i].lcm!=NULL)
1840        {
1841          if (pGetOrder((resPairs[*index])[i].lcm) > *actdeg)
1842            t = pGetOrder((resPairs[*index])[i].lcm);
1843        }
1844        else if ((resPairs[*index])[i].syz!=NULL)
1845        {
1846          if (pTotaldegree((resPairs[*index])[i].syz) > *actdeg)
1847            t = pTotaldegree((resPairs[*index])[i].syz);
1848        }
1849        if ((t>*actdeg) && ((newdeg==*actdeg) || (t<newdeg)))
1850        {
1851          newdeg = t;
1852          newindex = *index;
1853          break;
1854        }
1855        i++;
1856      } 
1857    }
1858    (*index)--;
1859  }
1860  if (newdeg>*actdeg)
1861  {
1862    *actdeg = newdeg;
1863    *index = newindex;
1864    return syChosePairs1(resPairs,Tl,index,howmuch,length,actdeg);
1865  }
1866  else return NULL;
1867}
1868
1869/*3
1870* statistics of the resolution
1871*/
1872static void syStatistics(resolvente res,int length)
1873{
1874  int i,j=1,k,deg=0;
1875
1876  PrintLn();
1877  while ((j<length) && (res[j]!=NULL))
1878  {
1879    Print("In module %d: \n",j);
1880    k = 0;
1881    while ((k<IDELEMS(res[j])) && (res[j]->m[k]!=NULL))
1882    {
1883      i = 1;
1884      deg = pGetOrder(res[j]->m[k]);
1885      k++;
1886      while ((k<IDELEMS(res[j])) && (res[j]->m[k]!=NULL) &&
1887              (pGetOrder(res[j]->m[k])==deg))
1888      {
1889        i++;
1890        k++;
1891      }
1892      Print("%d elements of degree %d\n",i,deg);
1893    }
1894    j++;
1895  }
1896}
1897
1898/*3
1899* sets regularity computed from the leading terms of arg
1900*/
1901static int sySetRegularity(ideal arg)
1902{
1903  if (idIs0(arg)) return -1;
1904  ideal temp=idInit(IDELEMS(arg),arg->rank);
1905  int i,len,reg=-1;
1906//  intvec * w=new intvec(2);
1907
1908  for (i=IDELEMS(arg)-1;i>=0;i--)
1909  {
1910    if (arg->m[i]!=NULL)
1911      temp->m[i] = pHead(arg->m[i]);
1912  }
1913  idSkipZeroes(temp);
1914  resolvente res = sySchreyerResolvente(temp,-1,&len,TRUE,TRUE);
1915  intvec * dummy = syBetti(res,len,&reg, NULL);
1916  delete dummy;
1917  for (i=0;i<len;i++)
1918  {
1919    if (res[i]!=NULL) idDelete(&(res[i]));
1920  }
1921  Free((ADDRESS)res,len*sizeof(ideal));
1922  idDelete(&temp);
1923  return reg;
1924}
1925
1926/*3
1927* initialize a module
1928*/
1929static int syInitSyzMod(resolvente res,resolvente orderedRes,
1930                         int index)
1931{
1932  int result;
1933
1934  if (res[index]==NULL) 
1935  {
1936    res[index] = idInit(16,1);
1937    truecomponents[index] = (int*)Alloc0(17*sizeof(int));
1938    backcomponents[index] = (int*)Alloc0(17*sizeof(int));
1939    Howmuch[index] = (int*)Alloc0(17*sizeof(int));
1940    Firstelem[index] = (int*)Alloc0(17*sizeof(int));
1941//Print("sort %d has now size %d\n",index,17);
1942    orderedRes[index] = idInit(16,1);
1943    result = 0;
1944  }
1945  else
1946  {
1947    result = IDELEMS(res[index]);
1948    while ((result>0) && (res[index]->m[result-1]==NULL)) result--;
1949  }
1950  return result;
1951}
1952
1953/*2
1954* implementation of LaScala's algorithm
1955* assumes that the given module is homogeneous
1956* works with slanted degree, uses syChosePairs
1957*/
1958resolvente syLaScala1(ideal arg,int * length)
1959{
1960  BOOLEAN noPair=FALSE;
1961  int i,j,* ord,*b0,*b1,actdeg=32000,index=0,reg=-1;
1962  int startdeg,howmuch;
1963  intvec * Tl;
1964  poly p;
1965  ideal temp;
1966  resolvente res,orderedRes;
1967  SSet nextPairs;
1968  SRes resPairs;
1969  pSetmProc oldSetm=pSetm;
1970
1971  //crit = 0;
1972  //zeroRed = 0;
1973  //simple = 0;
1974  //dsim = 0;
1975  euler = -1;
1976  redpol = pNew();
1977  //orderingdepth = new intvec(pVariables+1);
1978  if (*length<=0) *length = pVariables+2;
1979  if (idIs0(arg)) return NULL;
1980  if ((currRing->order[0]==ringorder_dp)
1981  &&  (currRing->order[1]==ringorder_C)
1982  &&  (currRing->order[2]==0))
1983  {
1984    ord=NULL;
1985  } 
1986/*--- changes to a dpC-ring with special comp0------*/
1987  else
1988  {
1989    ord = (int*)Alloc0(3*sizeof(int));
1990    b0 = (int*)Alloc0(3*sizeof(int));
1991    b1 = (int*)Alloc0(3*sizeof(int));
1992    ord[0] = ringorder_dp;
1993    ord[1] = ringorder_C;
1994    b0[1] = 1;
1995    b1[0] = pVariables;
1996    pChangeRing(pVariables,1,ord,b0,b1,currRing->wvhdl);
1997  } 
1998/*--- initializes the data structures---------------*/
1999  Tl = new intvec(*length);
2000  temp = idInit(IDELEMS(arg),arg->rank);
2001  for (i=0;i<IDELEMS(arg);i++)
2002  {
2003    temp->m[i] = pOrdPoly(pCopy(arg->m[i]));
2004    if (temp->m[i]!=NULL) 
2005    {
2006      j = pTotaldegree(temp->m[i]);
2007      if (j<actdeg) actdeg = j;
2008    }
2009  }
2010  {
2011    highdeg=1;
2012    long long t=1;
2013    long long h_n=1+pVariables;
2014    while ((t=(((long long)t*(long long)h_n)/(long long)highdeg))<INT_MAX)
2015    {
2016      highdeg++;
2017      h_n++;
2018    }
2019    highdeg--;
2020    //Print("max deg=%d\n",highdeg);
2021  } 
2022  binomials = (int*)Alloc(pVariables*(highdeg+1)*sizeof(int));
2023  syBinomSet();
2024  pSetm =syzSetm;
2025  for (i=0;i<IDELEMS(arg);i++)
2026  {
2027    p = temp->m[i];
2028    while (p!=NULL)
2029    {
2030      pSetm(p);
2031      pIter(p);
2032    }
2033  }
2034  pComp0 = syzcomp2dpc;
2035  resPairs = syInitRes(temp,length,Tl);
2036  res = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2037  orderedRes = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2038  truecomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2039  backcomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2040  Howmuch = (int**)Alloc0((*length+1)*sizeof(int*));
2041  Firstelem = (int**)Alloc0((*length+1)*sizeof(int*));
2042  int len0=idRankFreeModule(arg)+1;
2043  truecomponents[0] = (int*)Alloc(len0*sizeof(int));
2044  backcomponents[0] = (int*)Alloc(len0*sizeof(int));
2045  Howmuch[0] = (int*)Alloc(len0*sizeof(int));
2046  Firstelem[0] = (int*)Alloc(len0*sizeof(int));
2047//Print("sort %d has now size %d\n",0,len0);
2048  for (i=0;i<len0;i++)
2049    truecomponents[0][i] = i;
2050  startdeg = actdeg;
2051  nextPairs = syChosePairs(resPairs,Tl,&index,&howmuch,*length,&actdeg);
2052  //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2053/*--- computes the resolution ----------------------*/ 
2054  while (nextPairs!=NULL)
2055  {
2056    if (TEST_OPT_PROT) Print("%d",actdeg);
2057    if (TEST_OPT_PROT) Print("(m%d)",index);
2058    currcomponents = truecomponents[max(index-1,0)];
2059    i = syInitSyzMod(res,orderedRes,index);
2060    j = syInitSyzMod(res,orderedRes,index+1);
2061    if (index>0)
2062    {
2063      syRedNextPairs(nextPairs,res,orderedRes,howmuch,index);
2064      syCompactifyPairSet(resPairs[index],(*Tl)[index],0);
2065    }
2066    else
2067      syRedGenerOfCurrDeg(resPairs,res,orderedRes,actdeg,index+1,Tl);
2068/*--- creates new pairs -----------------------------*/     
2069    syCreateNewPairs(resPairs,Tl,res,index,i);
2070    if (index<(*length)-1) 
2071    {
2072      syCreateNewPairs(resPairs,Tl,res,index+1,j);
2073      //currcomponents = truecomponents[index];
2074      //sySyzTail(resPairs,Tl,orderedRes,res,index+1,j);
2075      //currcomponents = truecomponents[index-1];
2076    }
2077    index++;
2078    nextPairs = syChosePairs(resPairs,Tl,&index,&howmuch,*length,&actdeg);
2079    //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2080  }
2081/*--- deletes temporary data structures-------------*/
2082  idDelete(&temp);
2083  for (i=0;i<*length;i++)
2084  {
2085    for (j=0;j<(*Tl)[i];j++)
2086    {
2087      if ((resPairs[i])[j].lcm!=NULL)
2088        pDelete(&(resPairs[i])[j].lcm);
2089    }
2090    if (orderedRes[i]!=NULL)
2091    {
2092      for (j=0;j<IDELEMS(orderedRes[i]);j++)
2093        orderedRes[i]->m[j] = NULL;
2094    }
2095    idDelete(&orderedRes[i]);
2096    if (truecomponents[i]!=NULL)
2097    {
2098      Free((ADDRESS)truecomponents[i],(IDELEMS(res[i])+1)*sizeof(int));
2099      truecomponents[i]=NULL;
2100    }
2101    if (backcomponents[i]!=NULL)
2102    {
2103      Free((ADDRESS)backcomponents[i],(IDELEMS(res[i])+1)*sizeof(int));
2104      backcomponents[i]=NULL;
2105    }
2106    if (Howmuch[i]!=NULL)
2107    {
2108      Free((ADDRESS)Howmuch[i],(IDELEMS(res[i])+1)*sizeof(int));
2109      Howmuch[i]=NULL;
2110    }
2111    if (Firstelem[i]!=NULL)
2112    {
2113      Free((ADDRESS)Firstelem[i],(IDELEMS(res[i])+1)*sizeof(int));
2114      Firstelem[i]=NULL;
2115    }
2116    Free((ADDRESS)resPairs[i],(*Tl)[i]*sizeof(SObject));
2117  }
2118  Free((ADDRESS)resPairs,*length*sizeof(SObject*));
2119  Free((ADDRESS)orderedRes,(*length+1)*sizeof(ideal));
2120  Free((ADDRESS)truecomponents,(*length+1)*sizeof(int*));
2121  Free((ADDRESS)backcomponents,(*length+1)*sizeof(int*));
2122  Free((ADDRESS)Howmuch,(*length+1)*sizeof(int*));
2123  Free((ADDRESS)Firstelem,(*length+1)*sizeof(int*));
2124  truecomponents = NULL;
2125  backcomponents = NULL;
2126  Howmuch = NULL;
2127  Firstelem = NULL;
2128  if (BTEST1(6)) syStatistics(res,(*length+1));
2129  (*length)++;
2130/*--- changes to the original ring------------------*/
2131  if (ord!=NULL)
2132  {
2133    pChangeRing(pVariables,currRing->OrdSgn,currRing->order,
2134    currRing->block0,currRing->block1,currRing->wvhdl);
2135  }
2136  else
2137  {
2138    pSetm=oldSetm;
2139  }
2140  syReOrderResolventFB(res,*length,2);
2141  for (i=0;i<*length-1;i++)
2142  {
2143    res[i] = res[i+1];
2144    if (res[i]!=NULL)
2145    {
2146      if (i>0)
2147      {
2148        for (j=0;j<IDELEMS(res[i]);j++)
2149        {
2150          res[i]->m[j] = syOrdPolySchreyer(res[i]->m[j]);
2151        }
2152        idSkipZeroes(res[i]);
2153      }
2154      else
2155      {
2156        for (j=0;j<IDELEMS(res[i]);j++)
2157        {
2158          res[i]->m[j] = pOrdPoly(res[i]->m[j]);
2159        }
2160        idSkipZeroes(res[i]);
2161      }
2162    }
2163  }
2164  res[*length-1] = NULL;
2165  if (ord!=NULL)
2166  {
2167    Free((ADDRESS)ord,3*sizeof(int));
2168    Free((ADDRESS)b0,3*sizeof(int));
2169    Free((ADDRESS)b1,3*sizeof(int));
2170  } 
2171  Free((ADDRESS)binomials,pVariables*(highdeg_1)*sizeof(int));
2172  pDelete1(&redpol);
2173  if (TEST_OPT_PROT) 
2174  {
2175    //Print("simple: %d\n",simple);
2176    //Print("dsim: %d\n",dsim);
2177    //Print("crit %d-times used \n",crit);
2178    //Print("%d reductions to zero \n",zeroRed);
2179  }
2180  //delete orderingdepth;
2181  return res;
2182}
2183
2184/*2
2185* second implementation of LaScala's algorithm
2186* assumes that the given module is homogeneous
2187* works deg by deg, uses syChosePairs1
2188*/
2189resolvente syLaScala2(ideal arg,int * length)
2190{
2191  BOOLEAN noPair=FALSE;
2192  int i,j,* ord,*b0,*b1,actdeg=32000,index=1,reg=-1;
2193  int startdeg,howmuch;
2194  intvec * Tl;
2195  poly p;
2196  ideal temp;
2197  resolvente res,orderedRes;
2198  SSet nextPairs;
2199  SRes resPairs;
2200
2201  //crit = 0;
2202  //zeroRed = 0;
2203  //simple = 0;
2204  //dsim = 0;
2205  redpol = pNew();
2206  //orderingdepth = new intvec(pVariables+1);
2207  if (*length<=0) *length = pVariables+2;
2208  if (idIs0(arg)) return NULL;
2209/*--- changes to a dpc-ring with special comp0------*/
2210  ord = (int*)Alloc0(3*sizeof(int));
2211  b0 = (int*)Alloc0(3*sizeof(int));
2212  b1 = (int*)Alloc0(3*sizeof(int));
2213  ord[0] = ringorder_dp;
2214  ord[1] = ringorder_C;
2215  b0[1] = 1;
2216  b1[0] = pVariables;
2217  pChangeRing(pVariables,1,ord,b0,b1,currRing->wvhdl);
2218/*--- initializes the data structures---------------*/
2219  Tl = new intvec(*length);
2220  temp = idInit(IDELEMS(arg),arg->rank);
2221  for (i=0;i<IDELEMS(arg);i++)
2222  {
2223    temp->m[i] = pOrdPoly(pCopy(arg->m[i]));
2224    if (temp->m[i]!=NULL) 
2225    {
2226      j = pTotaldegree(temp->m[i]);
2227      if (j<actdeg) actdeg = j;
2228    }
2229  }
2230  {
2231    highdeg=1;
2232    long long t=1;
2233    long long h_n=1+pVariables;
2234    while ((t=(((long long)t*(long long)h_n)/(long long)highdeg))<INT_MAX)
2235    {
2236      highdeg++;
2237      h_n++;
2238    }
2239    highdeg--;
2240    //Print("max deg=%d\n",highdeg);
2241  } 
2242  binomials = (int*)Alloc(pVariables*(highdeg+1)*sizeof(int));
2243  syBinomSet();
2244  pSetm =syzSetm;
2245  for (i=0;i<IDELEMS(arg);i++)
2246  {
2247    p = temp->m[i];
2248    while (p!=NULL)
2249    {
2250      pSetm(p);
2251      pIter(p);
2252    }
2253  }
2254  pComp0 = syzcomp2dpc;
2255  resPairs = syInitRes(temp,length,Tl);
2256  res = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2257  orderedRes = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2258  truecomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2259  backcomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2260  Howmuch = (int**)Alloc0((*length+1)*sizeof(int*));
2261  Firstelem = (int**)Alloc0((*length+1)*sizeof(int*));
2262  int len0=idRankFreeModule(arg)+1;
2263  truecomponents[0] = (int*)Alloc(len0*sizeof(int));
2264  backcomponents[0] = (int*)Alloc(len0*sizeof(int));
2265  Howmuch[0] = (int*)Alloc(len0*sizeof(int));
2266  Firstelem[0] = (int*)Alloc(len0*sizeof(int));
2267//Print("sort %d has now size %d\n",0,len0);
2268  for (i=0;i<len0;i++)
2269    truecomponents[0][i] = i;
2270  startdeg = actdeg;
2271  nextPairs = syChosePairs1(resPairs,Tl,&index,&howmuch,*length,&actdeg);
2272  if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2273/*--- computes the resolution ----------------------*/ 
2274  while (nextPairs!=NULL)
2275  {
2276    if (TEST_OPT_PROT) Print("%d",actdeg);
2277    if (TEST_OPT_PROT) Print("(m%d)",index);
2278    currcomponents = truecomponents[max(index-1,0)];
2279    i = syInitSyzMod(res,orderedRes,index);
2280    j = syInitSyzMod(res,orderedRes,index+1);
2281    if (index>0)
2282    {
2283      syRedNextPairs(nextPairs,res,orderedRes,howmuch,index);
2284      syCompactifyPairSet(resPairs[index],(*Tl)[index],0);
2285    }
2286    else
2287      syRedGenerOfCurrDeg(resPairs,res,orderedRes,actdeg,index+1,Tl);
2288/*--- creates new pairs -----------------------------*/     
2289    syCreateNewPairs(resPairs,Tl,res,index,i);
2290    if (index<(*length)-1) 
2291    {
2292      syCreateNewPairs(resPairs,Tl,res,index+1,j);
2293      currcomponents = truecomponents[index];
2294      sySyzTail(resPairs,Tl,orderedRes,res,index+1,j);
2295      currcomponents = truecomponents[index-1];
2296    }
2297    index--;
2298    nextPairs = syChosePairs1(resPairs,Tl,&index,&howmuch,*length,&actdeg);
2299    if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2300  }
2301/*--- deletes temporary data structures-------------*/
2302  idDelete(&temp);
2303  for (i=0;i<*length;i++)
2304  {
2305    for (j=0;j<(*Tl)[i];j++)
2306    {
2307      if ((resPairs[i])[j].lcm!=NULL)
2308        pDelete(&(resPairs[i])[j].lcm);
2309    }
2310    if (orderedRes[i]!=NULL)
2311    {
2312      for (j=0;j<IDELEMS(orderedRes[i]);j++)
2313        orderedRes[i]->m[j] = NULL;
2314    }
2315    idDelete(&orderedRes[i]);
2316    if (truecomponents[i]!=NULL)
2317    {
2318      Free((ADDRESS)truecomponents[i],(IDELEMS(res[i])+1)*sizeof(int));
2319      truecomponents[i]=NULL;
2320    }
2321    if (backcomponents[i]!=NULL)
2322    {
2323      Free((ADDRESS)backcomponents[i],(IDELEMS(res[i])+1)*sizeof(int));
2324      backcomponents[i]=NULL;
2325    }
2326    if (Howmuch[i]!=NULL)
2327    {
2328      Free((ADDRESS)Howmuch[i],(IDELEMS(res[i])+1)*sizeof(int));
2329      Howmuch[i]=NULL;
2330    }
2331    if (Firstelem[i]!=NULL)
2332    {
2333      Free((ADDRESS)Firstelem[i],(IDELEMS(res[i])+1)*sizeof(int));
2334      Firstelem[i]=NULL;
2335    }
2336    Free((ADDRESS)resPairs[i],(*Tl)[i]*sizeof(SObject));
2337  }
2338  Free((ADDRESS)resPairs,*length*sizeof(SObject*));
2339  Free((ADDRESS)orderedRes,(*length+1)*sizeof(ideal));
2340  Free((ADDRESS)truecomponents,(*length+1)*sizeof(int*));
2341  Free((ADDRESS)backcomponents,(*length+1)*sizeof(int*));
2342  Free((ADDRESS)Howmuch,(*length+1)*sizeof(int*));
2343  Free((ADDRESS)Firstelem,(*length+1)*sizeof(int*));
2344  truecomponents = NULL;
2345  backcomponents = NULL;
2346  Howmuch = NULL;
2347  Firstelem = NULL;
2348  if (BTEST1(6)) syStatistics(res,(*length+1));
2349  (*length)++;
2350/*--- changes to the original ring------------------*/
2351  pChangeRing(pVariables,currRing->OrdSgn,currRing->order,
2352    currRing->block0,currRing->block1,currRing->wvhdl);
2353  syReOrderResolventFB(res,*length,2);
2354  for (i=0;i<*length-1;i++)
2355  {
2356    res[i] = res[i+1];
2357    if (res[i]!=NULL)
2358    {
2359      if (i>0)
2360      {
2361        for (j=0;j<IDELEMS(res[i]);j++)
2362        {
2363          res[i]->m[j] = syOrdPolySchreyer(res[i]->m[j]);
2364        }
2365        idSkipZeroes(res[i]);
2366      }
2367      else
2368      {
2369        for (j=0;j<IDELEMS(res[i]);j++)
2370        {
2371          p = res[i]->m[j];
2372          while (p!=NULL)
2373          {
2374            pSetm(p);
2375            pIter(p);
2376          }
2377        }
2378        idSkipZeroes(res[i]);
2379      }
2380    }
2381  }
2382  res[*length-1] = NULL;
2383  Free((ADDRESS)ord,3*sizeof(int));
2384  Free((ADDRESS)b0,3*sizeof(int));
2385  Free((ADDRESS)b1,3*sizeof(int));
2386  Free((ADDRESS)binomials,pVariables*(highdeg_1)*sizeof(int));
2387  pDelete1(&redpol);
2388  if (TEST_OPT_PROT) 
2389  {
2390    //Print("simple: %d\n",simple);
2391    //Print("dsim: %d\n",dsim);
2392    //Print("crit %d-times used \n",crit);
2393    //Print("%d reductions to zero \n",zeroRed);
2394  }
2395  //delete orderingdepth;
2396  return res;
2397}
Note: See TracBrowser for help on using the repository browser.