source: git/Singular/syz1.cc @ c3c413

spielwiese
Last change on this file since c3c413 was c3c413, checked in by Hans Schönemann <hannes@…>, 27 years ago
* hannes/siebert: minor bug fixes in k*.cc, syz1.cc added degBound to sres git-svn-id: file:///usr/local/Singular/svn/trunk@374 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.8 1997-06-09 12:21:29 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 (TEST_OPT_DEGBOUND && (newdeg > Kstd1_deg))
1756    return NULL;
1757  if (newdeg>*actdeg)
1758  {
1759    *actdeg = newdeg;
1760    *index = newindex;
1761    return syChosePairs(resPairs,Tl,index,howmuch,length,actdeg);
1762  }
1763  else return NULL;
1764}
1765
1766/*3
1767* looks through the pair set and the given module for
1768* remaining pairs or generators to consider
1769* returns a pointer to the first pair and the number of them in the given module
1770* works deg by deg
1771*/
1772static SSet syChosePairs1(SRes resPairs,intvec * Tl, int *index, int *howmuch,
1773                   int length,int * actdeg)
1774{
1775  int newdeg=*actdeg,newindex=-1,i,t;
1776  SSet result;
1777 
1778  while (*index>=0)
1779  {
1780    if (resPairs[*index]!=NULL)
1781    {
1782      i = 0;
1783      if (*index!=0)
1784      {
1785        while ((i<(*Tl)[*index]))
1786        {
1787          if ((resPairs[*index])[i].lcm!=NULL)
1788          {
1789            if (pGetOrder((resPairs[*index])[i].lcm) == *actdeg)
1790            {
1791              result = &(resPairs[*index])[i];
1792              *howmuch =1;
1793              i++;
1794              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].lcm!=NULL)
1795                      && (pGetOrder((resPairs[*index])[i].lcm) == *actdeg))
1796              {
1797                i++;
1798                (*howmuch)++;
1799              }
1800              return result;
1801            }
1802          }
1803          i++;
1804        }
1805      }
1806      else
1807      {
1808        while ((i<(*Tl)[*index]))
1809        {
1810          if ((resPairs[*index])[i].syz!=NULL)
1811          {
1812            if (pTotaldegree((resPairs[*index])[i].syz) == *actdeg)
1813            {
1814              result = &(resPairs[*index])[i];
1815              (*howmuch) =1;
1816              i++;
1817              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].syz!=NULL)
1818                      && (pTotaldegree((resPairs[*index])[i].syz) == *actdeg))
1819              {
1820                i++;
1821                (*howmuch)++;
1822              }
1823              return result;
1824            }
1825          }
1826          i++;
1827        }
1828      }
1829    }
1830    (*index)--;
1831  }
1832  *index = length-1;
1833  while (*index>=0)
1834  {
1835    if (resPairs[*index]!=NULL)
1836    {
1837      i = 0;
1838      while ((i<(*Tl)[*index]))
1839      {
1840        t = *actdeg;
1841        if ((resPairs[*index])[i].lcm!=NULL)
1842        {
1843          if (pGetOrder((resPairs[*index])[i].lcm) > *actdeg)
1844            t = pGetOrder((resPairs[*index])[i].lcm);
1845        }
1846        else if ((resPairs[*index])[i].syz!=NULL)
1847        {
1848          if (pTotaldegree((resPairs[*index])[i].syz) > *actdeg)
1849            t = pTotaldegree((resPairs[*index])[i].syz);
1850        }
1851        if ((t>*actdeg) && ((newdeg==*actdeg) || (t<newdeg)))
1852        {
1853          newdeg = t;
1854          newindex = *index;
1855          break;
1856        }
1857        i++;
1858      } 
1859    }
1860    (*index)--;
1861  }
1862  if (newdeg>*actdeg)
1863  {
1864    *actdeg = newdeg;
1865    *index = newindex;
1866    return syChosePairs1(resPairs,Tl,index,howmuch,length,actdeg);
1867  }
1868  else return NULL;
1869}
1870
1871/*3
1872* statistics of the resolution
1873*/
1874static void syStatistics(resolvente res,int length)
1875{
1876  int i,j=1,k,deg=0;
1877
1878  PrintLn();
1879  while ((j<length) && (res[j]!=NULL))
1880  {
1881    Print("In module %d: \n",j);
1882    k = 0;
1883    while ((k<IDELEMS(res[j])) && (res[j]->m[k]!=NULL))
1884    {
1885      i = 1;
1886      deg = pGetOrder(res[j]->m[k]);
1887      k++;
1888      while ((k<IDELEMS(res[j])) && (res[j]->m[k]!=NULL) &&
1889              (pGetOrder(res[j]->m[k])==deg))
1890      {
1891        i++;
1892        k++;
1893      }
1894      Print("%d elements of degree %d\n",i,deg);
1895    }
1896    j++;
1897  }
1898}
1899
1900/*3
1901* sets regularity computed from the leading terms of arg
1902*/
1903static int sySetRegularity(ideal arg)
1904{
1905  if (idIs0(arg)) return -1;
1906  ideal temp=idInit(IDELEMS(arg),arg->rank);
1907  int i,len,reg=-1;
1908//  intvec * w=new intvec(2);
1909
1910  for (i=IDELEMS(arg)-1;i>=0;i--)
1911  {
1912    if (arg->m[i]!=NULL)
1913      temp->m[i] = pHead(arg->m[i]);
1914  }
1915  idSkipZeroes(temp);
1916  resolvente res = sySchreyerResolvente(temp,-1,&len,TRUE,TRUE);
1917  intvec * dummy = syBetti(res,len,&reg, NULL);
1918  delete dummy;
1919  for (i=0;i<len;i++)
1920  {
1921    if (res[i]!=NULL) idDelete(&(res[i]));
1922  }
1923  Free((ADDRESS)res,len*sizeof(ideal));
1924  idDelete(&temp);
1925  return reg;
1926}
1927
1928/*3
1929* initialize a module
1930*/
1931static int syInitSyzMod(resolvente res,resolvente orderedRes,
1932                         int index)
1933{
1934  int result;
1935
1936  if (res[index]==NULL) 
1937  {
1938    res[index] = idInit(16,1);
1939    truecomponents[index] = (int*)Alloc0(17*sizeof(int));
1940    backcomponents[index] = (int*)Alloc0(17*sizeof(int));
1941    Howmuch[index] = (int*)Alloc0(17*sizeof(int));
1942    Firstelem[index] = (int*)Alloc0(17*sizeof(int));
1943//Print("sort %d has now size %d\n",index,17);
1944    orderedRes[index] = idInit(16,1);
1945    result = 0;
1946  }
1947  else
1948  {
1949    result = IDELEMS(res[index]);
1950    while ((result>0) && (res[index]->m[result-1]==NULL)) result--;
1951  }
1952  return result;
1953}
1954
1955/*2
1956* implementation of LaScala's algorithm
1957* assumes that the given module is homogeneous
1958* works with slanted degree, uses syChosePairs
1959*/
1960resolvente syLaScala1(ideal arg,int * length)
1961{
1962  BOOLEAN noPair=FALSE;
1963  int i,j,* ord,*b0,*b1,actdeg=32000,index=0,reg=-1;
1964  int startdeg,howmuch;
1965  intvec * Tl;
1966  poly p;
1967  ideal temp;
1968  resolvente res,orderedRes;
1969  SSet nextPairs;
1970  SRes resPairs;
1971  pSetmProc oldSetm=pSetm;
1972
1973  //crit = 0;
1974  //zeroRed = 0;
1975  //simple = 0;
1976  //dsim = 0;
1977  euler = -1;
1978  redpol = pNew();
1979  //orderingdepth = new intvec(pVariables+1);
1980  if (*length<=0) *length = pVariables+2;
1981  if (idIs0(arg)) return NULL;
1982  if ((currRing->order[0]==ringorder_dp)
1983  &&  (currRing->order[1]==ringorder_C)
1984  &&  (currRing->order[2]==0))
1985  {
1986    ord=NULL;
1987  } 
1988/*--- changes to a dpC-ring with special comp0------*/
1989  else
1990  {
1991    ord = (int*)Alloc0(3*sizeof(int));
1992    b0 = (int*)Alloc0(3*sizeof(int));
1993    b1 = (int*)Alloc0(3*sizeof(int));
1994    ord[0] = ringorder_dp;
1995    ord[1] = ringorder_C;
1996    b0[1] = 1;
1997    b1[0] = pVariables;
1998    pChangeRing(pVariables,1,ord,b0,b1,currRing->wvhdl);
1999  } 
2000/*--- initializes the data structures---------------*/
2001  Tl = new intvec(*length);
2002  temp = idInit(IDELEMS(arg),arg->rank);
2003  for (i=0;i<IDELEMS(arg);i++)
2004  {
2005    temp->m[i] = pOrdPoly(pCopy(arg->m[i]));
2006    if (temp->m[i]!=NULL) 
2007    {
2008      j = pTotaldegree(temp->m[i]);
2009      if (j<actdeg) actdeg = j;
2010    }
2011  }
2012  {
2013    highdeg=1;
2014    long long t=1;
2015    long long h_n=1+pVariables;
2016    while ((t=(((long long)t*(long long)h_n)/(long long)highdeg))<INT_MAX)
2017    {
2018      highdeg++;
2019      h_n++;
2020    }
2021    highdeg--;
2022    //Print("max deg=%d\n",highdeg);
2023  } 
2024  binomials = (int*)Alloc(pVariables*(highdeg+1)*sizeof(int));
2025  syBinomSet();
2026  pSetm =syzSetm;
2027  for (i=0;i<IDELEMS(arg);i++)
2028  {
2029    p = temp->m[i];
2030    while (p!=NULL)
2031    {
2032      pSetm(p);
2033      pIter(p);
2034    }
2035  }
2036  pComp0 = syzcomp2dpc;
2037  resPairs = syInitRes(temp,length,Tl);
2038  res = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2039  orderedRes = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2040  truecomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2041  backcomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2042  Howmuch = (int**)Alloc0((*length+1)*sizeof(int*));
2043  Firstelem = (int**)Alloc0((*length+1)*sizeof(int*));
2044  int len0=idRankFreeModule(arg)+1;
2045  truecomponents[0] = (int*)Alloc(len0*sizeof(int));
2046  backcomponents[0] = (int*)Alloc(len0*sizeof(int));
2047  Howmuch[0] = (int*)Alloc(len0*sizeof(int));
2048  Firstelem[0] = (int*)Alloc(len0*sizeof(int));
2049//Print("sort %d has now size %d\n",0,len0);
2050  for (i=0;i<len0;i++)
2051    truecomponents[0][i] = i;
2052  startdeg = actdeg;
2053  nextPairs = syChosePairs(resPairs,Tl,&index,&howmuch,*length,&actdeg);
2054  //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2055/*--- computes the resolution ----------------------*/ 
2056  while (nextPairs!=NULL)
2057  {
2058    if (TEST_OPT_PROT) Print("%d",actdeg);
2059    if (TEST_OPT_PROT) Print("(m%d)",index);
2060    currcomponents = truecomponents[max(index-1,0)];
2061    i = syInitSyzMod(res,orderedRes,index);
2062    j = syInitSyzMod(res,orderedRes,index+1);
2063    if (index>0)
2064    {
2065      syRedNextPairs(nextPairs,res,orderedRes,howmuch,index);
2066      syCompactifyPairSet(resPairs[index],(*Tl)[index],0);
2067    }
2068    else
2069      syRedGenerOfCurrDeg(resPairs,res,orderedRes,actdeg,index+1,Tl);
2070/*--- creates new pairs -----------------------------*/     
2071    syCreateNewPairs(resPairs,Tl,res,index,i);
2072    if (index<(*length)-1) 
2073    {
2074      syCreateNewPairs(resPairs,Tl,res,index+1,j);
2075      //currcomponents = truecomponents[index];
2076      //sySyzTail(resPairs,Tl,orderedRes,res,index+1,j);
2077      //currcomponents = truecomponents[index-1];
2078    }
2079    index++;
2080    nextPairs = syChosePairs(resPairs,Tl,&index,&howmuch,*length,&actdeg);
2081    //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2082  }
2083/*--- deletes temporary data structures-------------*/
2084  idDelete(&temp);
2085  for (i=0;i<*length;i++)
2086  {
2087    for (j=0;j<(*Tl)[i];j++)
2088    {
2089      if ((resPairs[i])[j].lcm!=NULL)
2090        pDelete(&(resPairs[i])[j].lcm);
2091    }
2092    if (orderedRes[i]!=NULL)
2093    {
2094      for (j=0;j<IDELEMS(orderedRes[i]);j++)
2095        orderedRes[i]->m[j] = NULL;
2096    }
2097    idDelete(&orderedRes[i]);
2098    if (truecomponents[i]!=NULL)
2099    {
2100      Free((ADDRESS)truecomponents[i],(IDELEMS(res[i])+1)*sizeof(int));
2101      truecomponents[i]=NULL;
2102    }
2103    if (backcomponents[i]!=NULL)
2104    {
2105      Free((ADDRESS)backcomponents[i],(IDELEMS(res[i])+1)*sizeof(int));
2106      backcomponents[i]=NULL;
2107    }
2108    if (Howmuch[i]!=NULL)
2109    {
2110      Free((ADDRESS)Howmuch[i],(IDELEMS(res[i])+1)*sizeof(int));
2111      Howmuch[i]=NULL;
2112    }
2113    if (Firstelem[i]!=NULL)
2114    {
2115      Free((ADDRESS)Firstelem[i],(IDELEMS(res[i])+1)*sizeof(int));
2116      Firstelem[i]=NULL;
2117    }
2118    Free((ADDRESS)resPairs[i],(*Tl)[i]*sizeof(SObject));
2119  }
2120  Free((ADDRESS)resPairs,*length*sizeof(SObject*));
2121  Free((ADDRESS)orderedRes,(*length+1)*sizeof(ideal));
2122  Free((ADDRESS)truecomponents,(*length+1)*sizeof(int*));
2123  Free((ADDRESS)backcomponents,(*length+1)*sizeof(int*));
2124  Free((ADDRESS)Howmuch,(*length+1)*sizeof(int*));
2125  Free((ADDRESS)Firstelem,(*length+1)*sizeof(int*));
2126  truecomponents = NULL;
2127  backcomponents = NULL;
2128  Howmuch = NULL;
2129  Firstelem = NULL;
2130  if (BTEST1(6)) syStatistics(res,(*length+1));
2131  (*length)++;
2132/*--- changes to the original ring------------------*/
2133  if (ord!=NULL)
2134  {
2135    pChangeRing(pVariables,currRing->OrdSgn,currRing->order,
2136    currRing->block0,currRing->block1,currRing->wvhdl);
2137  }
2138  else
2139  {
2140    pSetm=oldSetm;
2141  }
2142  syReOrderResolventFB(res,*length,2);
2143  for (i=0;i<*length-1;i++)
2144  {
2145    res[i] = res[i+1];
2146    if (res[i]!=NULL)
2147    {
2148      if (i>0)
2149      {
2150        for (j=0;j<IDELEMS(res[i]);j++)
2151        {
2152          res[i]->m[j] = syOrdPolySchreyer(res[i]->m[j]);
2153        }
2154        idSkipZeroes(res[i]);
2155      }
2156      else
2157      {
2158        for (j=0;j<IDELEMS(res[i]);j++)
2159        {
2160          res[i]->m[j] = pOrdPoly(res[i]->m[j]);
2161        }
2162        idSkipZeroes(res[i]);
2163      }
2164    }
2165  }
2166  res[*length-1] = NULL;
2167  if (ord!=NULL)
2168  {
2169    Free((ADDRESS)ord,3*sizeof(int));
2170    Free((ADDRESS)b0,3*sizeof(int));
2171    Free((ADDRESS)b1,3*sizeof(int));
2172  } 
2173  Free((ADDRESS)binomials,pVariables*(highdeg_1)*sizeof(int));
2174  pDelete1(&redpol);
2175  if (TEST_OPT_PROT) 
2176  {
2177    //Print("simple: %d\n",simple);
2178    //Print("dsim: %d\n",dsim);
2179    //Print("crit %d-times used \n",crit);
2180    //Print("%d reductions to zero \n",zeroRed);
2181  }
2182  //delete orderingdepth;
2183  return res;
2184}
2185
2186/*2
2187* second implementation of LaScala's algorithm
2188* assumes that the given module is homogeneous
2189* works deg by deg, uses syChosePairs1
2190*/
2191resolvente syLaScala2(ideal arg,int * length)
2192{
2193  BOOLEAN noPair=FALSE;
2194  int i,j,* ord,*b0,*b1,actdeg=32000,index=1,reg=-1;
2195  int startdeg,howmuch;
2196  intvec * Tl;
2197  poly p;
2198  ideal temp;
2199  resolvente res,orderedRes;
2200  SSet nextPairs;
2201  SRes resPairs;
2202
2203  //crit = 0;
2204  //zeroRed = 0;
2205  //simple = 0;
2206  //dsim = 0;
2207  redpol = pNew();
2208  //orderingdepth = new intvec(pVariables+1);
2209  if (*length<=0) *length = pVariables+2;
2210  if (idIs0(arg)) return NULL;
2211/*--- changes to a dpc-ring with special comp0------*/
2212  ord = (int*)Alloc0(3*sizeof(int));
2213  b0 = (int*)Alloc0(3*sizeof(int));
2214  b1 = (int*)Alloc0(3*sizeof(int));
2215  ord[0] = ringorder_dp;
2216  ord[1] = ringorder_C;
2217  b0[1] = 1;
2218  b1[0] = pVariables;
2219  pChangeRing(pVariables,1,ord,b0,b1,currRing->wvhdl);
2220/*--- initializes the data structures---------------*/
2221  Tl = new intvec(*length);
2222  temp = idInit(IDELEMS(arg),arg->rank);
2223  for (i=0;i<IDELEMS(arg);i++)
2224  {
2225    temp->m[i] = pOrdPoly(pCopy(arg->m[i]));
2226    if (temp->m[i]!=NULL) 
2227    {
2228      j = pTotaldegree(temp->m[i]);
2229      if (j<actdeg) actdeg = j;
2230    }
2231  }
2232  {
2233    highdeg=1;
2234    long long t=1;
2235    long long h_n=1+pVariables;
2236    while ((t=(((long long)t*(long long)h_n)/(long long)highdeg))<INT_MAX)
2237    {
2238      highdeg++;
2239      h_n++;
2240    }
2241    highdeg--;
2242    //Print("max deg=%d\n",highdeg);
2243  } 
2244  binomials = (int*)Alloc(pVariables*(highdeg+1)*sizeof(int));
2245  syBinomSet();
2246  pSetm =syzSetm;
2247  for (i=0;i<IDELEMS(arg);i++)
2248  {
2249    p = temp->m[i];
2250    while (p!=NULL)
2251    {
2252      pSetm(p);
2253      pIter(p);
2254    }
2255  }
2256  pComp0 = syzcomp2dpc;
2257  resPairs = syInitRes(temp,length,Tl);
2258  res = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2259  orderedRes = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2260  truecomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2261  backcomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2262  Howmuch = (int**)Alloc0((*length+1)*sizeof(int*));
2263  Firstelem = (int**)Alloc0((*length+1)*sizeof(int*));
2264  int len0=idRankFreeModule(arg)+1;
2265  truecomponents[0] = (int*)Alloc(len0*sizeof(int));
2266  backcomponents[0] = (int*)Alloc(len0*sizeof(int));
2267  Howmuch[0] = (int*)Alloc(len0*sizeof(int));
2268  Firstelem[0] = (int*)Alloc(len0*sizeof(int));
2269//Print("sort %d has now size %d\n",0,len0);
2270  for (i=0;i<len0;i++)
2271    truecomponents[0][i] = i;
2272  startdeg = actdeg;
2273  nextPairs = syChosePairs1(resPairs,Tl,&index,&howmuch,*length,&actdeg);
2274  if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2275/*--- computes the resolution ----------------------*/ 
2276  while (nextPairs!=NULL)
2277  {
2278    if (TEST_OPT_PROT) Print("%d",actdeg);
2279    if (TEST_OPT_PROT) Print("(m%d)",index);
2280    currcomponents = truecomponents[max(index-1,0)];
2281    i = syInitSyzMod(res,orderedRes,index);
2282    j = syInitSyzMod(res,orderedRes,index+1);
2283    if (index>0)
2284    {
2285      syRedNextPairs(nextPairs,res,orderedRes,howmuch,index);
2286      syCompactifyPairSet(resPairs[index],(*Tl)[index],0);
2287    }
2288    else
2289      syRedGenerOfCurrDeg(resPairs,res,orderedRes,actdeg,index+1,Tl);
2290/*--- creates new pairs -----------------------------*/     
2291    syCreateNewPairs(resPairs,Tl,res,index,i);
2292    if (index<(*length)-1) 
2293    {
2294      syCreateNewPairs(resPairs,Tl,res,index+1,j);
2295      currcomponents = truecomponents[index];
2296      sySyzTail(resPairs,Tl,orderedRes,res,index+1,j);
2297      currcomponents = truecomponents[index-1];
2298    }
2299    index--;
2300    nextPairs = syChosePairs1(resPairs,Tl,&index,&howmuch,*length,&actdeg);
2301    if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2302  }
2303/*--- deletes temporary data structures-------------*/
2304  idDelete(&temp);
2305  for (i=0;i<*length;i++)
2306  {
2307    for (j=0;j<(*Tl)[i];j++)
2308    {
2309      if ((resPairs[i])[j].lcm!=NULL)
2310        pDelete(&(resPairs[i])[j].lcm);
2311    }
2312    if (orderedRes[i]!=NULL)
2313    {
2314      for (j=0;j<IDELEMS(orderedRes[i]);j++)
2315        orderedRes[i]->m[j] = NULL;
2316    }
2317    idDelete(&orderedRes[i]);
2318    if (truecomponents[i]!=NULL)
2319    {
2320      Free((ADDRESS)truecomponents[i],(IDELEMS(res[i])+1)*sizeof(int));
2321      truecomponents[i]=NULL;
2322    }
2323    if (backcomponents[i]!=NULL)
2324    {
2325      Free((ADDRESS)backcomponents[i],(IDELEMS(res[i])+1)*sizeof(int));
2326      backcomponents[i]=NULL;
2327    }
2328    if (Howmuch[i]!=NULL)
2329    {
2330      Free((ADDRESS)Howmuch[i],(IDELEMS(res[i])+1)*sizeof(int));
2331      Howmuch[i]=NULL;
2332    }
2333    if (Firstelem[i]!=NULL)
2334    {
2335      Free((ADDRESS)Firstelem[i],(IDELEMS(res[i])+1)*sizeof(int));
2336      Firstelem[i]=NULL;
2337    }
2338    Free((ADDRESS)resPairs[i],(*Tl)[i]*sizeof(SObject));
2339  }
2340  Free((ADDRESS)resPairs,*length*sizeof(SObject*));
2341  Free((ADDRESS)orderedRes,(*length+1)*sizeof(ideal));
2342  Free((ADDRESS)truecomponents,(*length+1)*sizeof(int*));
2343  Free((ADDRESS)backcomponents,(*length+1)*sizeof(int*));
2344  Free((ADDRESS)Howmuch,(*length+1)*sizeof(int*));
2345  Free((ADDRESS)Firstelem,(*length+1)*sizeof(int*));
2346  truecomponents = NULL;
2347  backcomponents = NULL;
2348  Howmuch = NULL;
2349  Firstelem = NULL;
2350  if (BTEST1(6)) syStatistics(res,(*length+1));
2351  (*length)++;
2352/*--- changes to the original ring------------------*/
2353  pChangeRing(pVariables,currRing->OrdSgn,currRing->order,
2354    currRing->block0,currRing->block1,currRing->wvhdl);
2355  syReOrderResolventFB(res,*length,2);
2356  for (i=0;i<*length-1;i++)
2357  {
2358    res[i] = res[i+1];
2359    if (res[i]!=NULL)
2360    {
2361      if (i>0)
2362      {
2363        for (j=0;j<IDELEMS(res[i]);j++)
2364        {
2365          res[i]->m[j] = syOrdPolySchreyer(res[i]->m[j]);
2366        }
2367        idSkipZeroes(res[i]);
2368      }
2369      else
2370      {
2371        for (j=0;j<IDELEMS(res[i]);j++)
2372        {
2373          p = res[i]->m[j];
2374          while (p!=NULL)
2375          {
2376            pSetm(p);
2377            pIter(p);
2378          }
2379        }
2380        idSkipZeroes(res[i]);
2381      }
2382    }
2383  }
2384  res[*length-1] = NULL;
2385  Free((ADDRESS)ord,3*sizeof(int));
2386  Free((ADDRESS)b0,3*sizeof(int));
2387  Free((ADDRESS)b1,3*sizeof(int));
2388  Free((ADDRESS)binomials,pVariables*(highdeg_1)*sizeof(int));
2389  pDelete1(&redpol);
2390  if (TEST_OPT_PROT) 
2391  {
2392    //Print("simple: %d\n",simple);
2393    //Print("dsim: %d\n",dsim);
2394    //Print("crit %d-times used \n",crit);
2395    //Print("%d reductions to zero \n",zeroRed);
2396  }
2397  //delete orderingdepth;
2398  return res;
2399}
Note: See TracBrowser for help on using the repository browser.