source: git/Singular/syz1.cc @ c1489f2

fieker-DuValspielwiese
Last change on this file since c1489f2 was c1489f2, checked in by Olaf Bachmann <obachman@…>, 25 years ago
* stuff from Thomas * various bug fixes git-svn-id: file:///usr/local/Singular/svn/trunk@3730 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 67.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: syz1.cc,v 1.43 1999-10-19 12:42:49 obachman Exp $ */
5/*
6* ABSTRACT: resolutions
7*/
8
9#include "mod2.h"
10#include "tok.h"
11#include "mmemory.h"
12#include "polys.h"
13#include "febase.h"
14#include "kstd1.h"
15#include "kutil.h"
16#include "stairc.h"
17#include "ipid.h"
18#include "cntrlc.h"
19#include "ipid.h"
20#include "intvec.h"
21#include "ipshell.h"
22#include "limits.h"
23#include "numbers.h"
24#include "modulop.h"
25#include "ideals.h"
26#include "intvec.h"
27#include "ring.h"
28#include "lists.h"
29#include "syz.h"
30#include "kbuckets.h"
31#include "polys-comp.h"
32#include <limits.h>
33
34extern void rSetmS(poly p, int* Components, long* ShiftedComponents);
35
36/*--------------static variables------------------------*/
37/*---points to the real components, shifted of the actual module-*/
38int *  currcomponents=NULL;
39long *  currShiftedComponents=NULL;
40
41
42/*---head-term-polynomials for the reduction------------*/
43static poly redpol=NULL;
44/*---counts number of applications of GM-criteria-------*/
45//static int crit;
46//static int euler;
47
48/*3
49* deletes all entres of a pair
50*/
51void syDeletePair(SObject * so)
52{
53  pDelete(&(*so).p);
54  pDelete(&(*so).lcm);
55  pDelete(&(*so).syz);
56  (*so).p1 = NULL;
57  (*so).p2 = NULL;
58  (*so).ind1 = 0;
59  (*so).ind2 = 0;
60  (*so).syzind = -1;
61  (*so).order = 0;
62  (*so).isNotMinimal = NULL;
63  (*so).length = -1;
64  (*so).reference = -1;
65}
66
67/*3
68* initializes all entres of a pair
69*/
70void syInitializePair(SObject * so)
71{
72  (*so).p = NULL;
73  (*so).lcm = NULL;
74  (*so).syz = NULL;
75  (*so).p1 = NULL;
76  (*so).p2 = NULL;
77  (*so).ind1 = 0;
78  (*so).ind2 = 0;
79  (*so).syzind = -1;
80  (*so).order = 0;
81  (*so).isNotMinimal = NULL;
82  (*so).length = -1;
83  (*so).reference = -1;
84}
85
86/*3
87* puts all entres of a pair to another
88*/
89void syCopyPair(SObject * argso, SObject * imso)
90{
91  *imso=*argso;
92  (*argso).p = NULL;
93  (*argso).p1 = NULL;
94  (*argso).p2 = NULL;
95  (*argso).lcm = NULL;
96  (*argso).syz = NULL;
97  (*argso).ind1 = 0;
98  (*argso).ind2 = 0;
99  (*argso).syzind = -1;
100  (*argso).order = 0;
101  (*argso).isNotMinimal = NULL;
102  (*argso).length = -1;
103  (*argso).reference = -1;
104}
105
106/*3
107* deletes empty objects from a pair set beginning with
108* pair first
109* assumes a pair to be empty if .lcm does so
110*/
111void syCompactifyPairSet(SSet sPairs, int sPlength, int first)
112{
113  int k=first,kk=0;
114
115  while (k+kk<sPlength)
116  {
117    if (sPairs[k+kk].lcm!=NULL)
118    {
119      if (kk>0) syCopyPair(&sPairs[k+kk],&sPairs[k]);
120      k++;
121    }
122    else
123    {
124      kk++;
125    }
126  }
127  while (k<sPlength)
128  {
129    syInitializePair(&sPairs[k]);
130    k++;
131  }
132}
133
134/*3
135* deletes empty objects from a pair set beginning with
136* pair first
137* assumes a pair to be empty if .lcm does so
138*/
139void syCompactify1(SSet sPairs, int* sPlength, int first)
140{
141  int k=first,kk=0;
142
143  while (k+kk<*sPlength)
144  {
145    if (sPairs[k+kk].lcm!=NULL)
146    {
147      if (kk>0) syCopyPair(&sPairs[k+kk],&sPairs[k]);
148      k++;
149    }
150    else
151    {
152      kk++;
153    }
154  }
155  while (k<*sPlength)
156  {
157    syInitializePair(&sPairs[k]);
158    k++;
159  }
160  *sPlength -= kk;
161}
162
163static pCompProc syOldComp1;
164/*3
165* replaces comp1dpc during homogeneous syzygy-computations
166* compares with components of currcomponents instead of the
167* exp[0]
168*/
169static int syzcomp0dpc(poly p1, poly p2)
170{
171  /*4 compare monomials by order then revlex*/
172    int i = pVariables;
173    if ((pGetExp(p1,i) == pGetExp(p2,i)))
174    {
175      do
176      {
177        i--;
178        if (i <= 1)
179        {
180           /*4 handle module case:*/
181           if (pGetComp(p1)==pGetComp(p2)) return 0;
182           else if
183              (currcomponents[pGetComp(p1)]>currcomponents[pGetComp(p2)])
184                return 1;
185           else return -1;
186        }
187      } while ((pGetExp(p1,i) == pGetExp(p2,i)));
188    }
189    if (pGetExp(p1,i) < pGetExp(p2,i)) return 1;
190    return -1;
191}
192
193/*3
194* replaces comp1dpc during homogeneous syzygy-computations
195* compares with components of currcomponents instead of the
196* exp[0]
197*/
198int syzcomp1dpc(poly p1, poly p2)
199{
200  int i = pVariables;
201  while ((i>1) && (pGetExp(p1,i)==pGetExp(p2,i)))
202    i--;
203  if (i>1)
204  {
205    if (pGetExp(p1,i) < pGetExp(p2,i)) return 1;
206    return -1;
207  }
208  int o1=pGetComp(p1);
209  int o2=pGetComp(p2);
210  if (o1==o2/*pGetComp(p1)==pGetComp(p2)*/) return 0;
211  if (currcomponents[o1]>currcomponents[o2]) return 1;
212  return -1;
213
214}
215
216/*3
217* replaces comp1dpc during homogeneous syzygy-computations
218* compares with components of currcomponents instead of the
219* exp[0]
220*/
221
222#ifdef PDEBUG
223static int syzcomp2dpc_test(poly p1, poly p2)
224{
225  long c1, c2, cc1, cc2, ccc1, ccc2, ec1, ec2;
226  c1 = pGetComp(p1);
227  c2 = pGetComp(p2);
228  cc1 = currcomponents[c1];
229  cc2 = currcomponents[c2];
230  ccc1 = currShiftedComponents[cc1];
231  ccc2 = currShiftedComponents[cc2];
232  ec1 = p1->exp.l[currRing->typ[1].data.syzcomp.place];
233  ec2 = p2->exp.l[currRing->typ[1].data.syzcomp.place];
234
235  if (ec1 != ccc1)
236  {
237    Warn("Shifted comp of p1 out of sync. should %d, is %d", ccc1, ec1);
238    //mmDBInfoBlock(p1);
239  }
240  if (ec2 != ccc2)
241  {
242    Warn("Shifted comp of p2 out of sync. should %d, is %d", ccc2, ec2);
243    //mmDBInfoBlock(p2);
244  }
245
246  if (c1 == c2)
247  {
248    assume(ccc1 == ccc2);
249  }
250  else if (cc1 > cc2)
251  {
252    assume(ccc1 > ccc2);
253  }
254  else
255  {
256    assume (cc1 < cc2);
257    assume (ccc1 < ccc2);
258  }
259  int o1=pGetOrder(p1), o2=pGetOrder(p2);
260  if (o1 > o2) return 1;
261  if (o1 < o2) return -1;
262
263  //if (o1>0)
264  {
265    int i = pVariables;
266    while ((i>1) && (pGetExp(p1,i)==pGetExp(p2,i)))
267      i--;
268    //(*orderingdepth)[pVariables-i]++;
269    if (i>1)
270    {
271      if (pGetExp(p1,i) < pGetExp(p2,i)) return 1;
272      return -1;
273    }
274  }
275  o1=pGetComp(p1);
276  o2=pGetComp(p2);
277  if (o1==o2/*pGetComp(p1)==pGetComp(p2)*/) return 0;
278  if (currcomponents[o1]>currcomponents[o2]) return 1;
279  return -1;
280}
281int syzcomp2dpc(poly p1, poly p2)
282{
283  int test = syzcomp2dpc_test(p1, p2);
284  int rc = rComp0(p1, p2);
285  if (test != rc)
286  {
287    assume(0);
288    syzcomp2dpc_test(p1, p2);
289    rComp0(p1, p2);
290  }
291  return test;
292}
293#endif // PDEBUG
294
295/*3
296* compares only the monomial without component
297*/
298static int syzcompmonomdp(poly p1, poly p2)
299{
300  int i;
301
302  /*4 compare monomials by order then revlex*/
303  if (pGetOrder(p1) == pGetOrder(p2))
304  {
305    i = pVariables;
306    if ((pGetExp(p1,i) == pGetExp(p2,i)))
307    {
308      do
309      {
310        i--;
311        if (i <= 1)
312          return 0;
313      } while ((pGetExp(p1,i) == pGetExp(p2,i)));
314    }
315    if (pGetExp(p1,i) < pGetExp(p2,i)) return 1;
316    return -1;
317  }
318  else if (pGetOrder(p1) > pGetOrder(p2)) return 1;
319  return -1;
320}
321
322poly syOrdPolySchreyer(poly p)
323{
324  return pOrdPolyMerge(p);
325}
326
327static int syLength(poly p)
328{
329  int result=0;
330
331  while (p!=NULL)
332  {
333    result++;
334    pIter(p);
335  }
336  return result;
337}
338
339poly syRedtail (poly p, syStrategy syzstr, int index)
340{
341  poly h, hn;
342  int j,pos;
343  ideal redWith=syzstr->orderedRes[index];
344
345  h = p;
346  hn = pNext(h);
347  while(hn != NULL)
348  {
349    j = syzstr->Firstelem[index-1][pGetComp(hn)]-1;
350    if (j>=0)
351    {
352      pos = j+syzstr->Howmuch[index-1][pGetComp(hn)];
353      while (j < pos)
354      {
355        if (pDivisibleBy2(redWith->m[j], hn))
356        {
357          //hn = sySPolyRed(hn,redWith->m[j]);
358          hn = ksOldSpolyRed(redWith->m[j],hn);
359          if (hn == NULL)
360          {
361            pNext(h) = NULL;
362            return p;
363          }
364          j = syzstr->Firstelem[index-1][pGetComp(hn)]-1;
365          pos = j+syzstr->Howmuch[index-1][pGetComp(hn)];
366        }
367        else
368        {
369          j++;
370        }
371      }
372    }
373    h = pNext(h) = hn;
374    hn = pNext(h);
375  }
376  return p;
377}
378
379
380/*3
381* local procedure for of syInitRes for the module case
382*/
383static int syChMin(intvec * iv)
384{
385  int i,j=-1,r=-1;
386
387  for (i=iv->length()-1;i>=0;i--)
388  {
389    if ((*iv)[i]>=0)
390    {
391      if ((j<0) || ((*iv)[i]<j))
392      {
393        j = (*iv)[i];
394        r = i;
395      }
396    }
397  }
398  return r;
399}
400
401/*3
402* initialize the resolution and puts in the argument as
403* zeroth entre, length must be > 0
404* assumes that the basering is degree-compatible
405*/
406SRes syInitRes(ideal arg,int * length, intvec * Tl, intvec * cw)
407{
408  if (idIs0(arg)) return NULL;
409  SRes resPairs = (SRes)Alloc0(*length*sizeof(SSet));
410  resPairs[0] = (SSet)Alloc0(IDELEMS(arg)*sizeof(SObject));
411  intvec * iv=NULL;
412  int i,j;
413
414  if (idRankFreeModule(arg)==0)
415  {
416    iv = idSort(arg);
417    for (i=0;i<IDELEMS(arg);i++)
418    {
419      (resPairs[0])[i].syz = /*pCopy*/(arg->m[(*iv)[i]-1]);
420      arg->m[(*iv)[i]-1] = NULL;
421      (resPairs[0])[i].order = pTotaldegree((resPairs[0])[i].syz);
422    }
423  }
424  else
425  {
426    iv = new intvec(IDELEMS(arg),1,-1);
427    for (i=0;i<IDELEMS(arg);i++)
428    {
429      (*iv)[i] = pTotaldegree(arg->m[i])+(*cw)[pGetComp(arg->m[i])-1];
430    }
431    for (i=0;i<IDELEMS(arg);i++)
432    {
433      j = syChMin(iv);
434      if (j<0) break;
435      (resPairs[0])[i].syz = arg->m[j];
436      arg->m[j] = NULL;
437      (resPairs[0])[i].order = (*iv)[j];
438      (*iv)[j] = -1;
439    }
440  }
441  if (iv!=NULL)  delete iv;
442  (*Tl)[0] = IDELEMS(arg);
443  return resPairs;
444}
445
446// rearrange shifted components
447long syReorderShiftedComponents(long * sc, int n)
448{
449  long holes = 0;
450  int i;
451  long new_comps = 0, new_space, max;
452
453  // count number of holes
454  for (i=1; i<n; i++)
455  {
456    if (sc[i-1] + 1 < sc[i]) holes++;
457  }
458
459  if (LONG_MAX - SYZ_SHIFT_BASE <= sc[n-1])
460  {
461    // need new components
462    new_comps = (1 << SYZ_SHIFT_MAX_NEW_COMP_ESTIMATE) - 1;
463    max = LONG_MAX;
464  }
465  else
466  {
467    max = sc[n-1] + SYZ_SHIFT_BASE;
468  }
469
470  // no we arrange things such that
471  // (n - holes) + holes*new_space + new_comps*SYZ_SHIFT_BASE= LONG_MAX
472  new_space = (max - n + holes - new_comps*SYZ_SHIFT_BASE) / holes;
473
474  assume(new_space < SYZ_SHIFT_BASE && new_space >= 4);
475
476  long* tc = ( long*) Alloc(n*sizeof(long));
477  tc[0] = sc[0];
478  // rearrange things
479  for (i=1; i<n; i++)
480  {
481    if (sc[i-1] + 1 < sc[i])
482    {
483      tc[i] = tc[i-1] + new_space;
484    }
485    else
486    {
487      tc[i] = tc[i-1] + 1;
488    }
489    assume(tc[i] > tc[i-1]);
490  }
491
492  assume(LONG_MAX -  SYZ_SHIFT_BASE > tc[n-1]);
493#ifdef HAVE_ASSUME
494  for (i=1; i<n; i++)
495  {
496    assume(tc[i] >= 0);
497    assume(tc[i-1] + 1 <= tc[i]);
498  }
499#endif
500
501  memcpyW(sc, tc, n);
502  Free(tc, n*sizeof(long));
503  return new_space;
504}
505
506// this make a Setm on p
507static void pResetSetm(poly p)
508{
509#ifdef PDEBUG
510  poly q = p;
511#endif
512  while (p!= NULL)
513  {
514    pSetm(p);
515    pIter(p);
516  }
517#ifdef PDEBUG
518  if (! pTest(q))
519  {
520    pSetm(p);
521  }
522#endif
523}
524
525void syResetShiftedComponents(syStrategy syzstr, int index,int hilb)
526{
527  assume(index > 0);
528  int i;
529  if (syzstr->res[index] != NULL)
530  {
531    long * prev_s;
532    int* prev_c;
533    int p_length;
534    rGetSComps(&prev_c, &prev_s, &p_length);
535    currcomponents = syzstr->truecomponents[index-1];
536    currShiftedComponents = syzstr->ShiftedComponents[index-1];
537    rChangeSComps(currcomponents,
538                  currShiftedComponents,
539                  IDELEMS(syzstr->res[index-1]));
540    if (hilb==0)
541    {
542      ideal id = syzstr->res[index];
543      for (i=0; i<IDELEMS(id); i++)
544      {
545        pResetSetm(id->m[i]);
546      }
547    }
548    else if (hilb==1)
549    {
550      assume (index>1);
551      assume (syzstr->resPairs[index-1]!=NULL);
552      SSet Pairs=syzstr->resPairs[index-1];
553      SSet Pairs1=syzstr->resPairs[index]; 
554      int till=(*syzstr->Tl)[index-1];
555      for (i=0;i<till;i++)
556      {
557        if (Pairs[i].syz!=NULL)
558          pResetSetm(Pairs[i].syz);
559      }
560      till=(*syzstr->Tl)[index];
561      for (i=0;i<till;i++)
562      {
563        if (Pairs1[i].p!=NULL)
564          pResetSetm(Pairs1[i].p);
565      }
566    }
567    currcomponents  = prev_c;
568    currShiftedComponents = prev_s;
569    rChangeSComps(prev_c, prev_s, p_length);
570  }
571}
572
573/*3
574* determines the place of a polynomial in the right ordered resolution
575* set the vectors of truecomponents
576*/
577static BOOLEAN syOrder(poly p,syStrategy syzstr,int index,
578                    int realcomp)
579{
580  int i=IDELEMS(syzstr->res[index-1])+1,j=0,k,tc,orc,ie=realcomp-1;
581  int *trind1=syzstr->truecomponents[index-1];
582  int *trind=syzstr->truecomponents[index];
583  long *shind=syzstr->ShiftedComponents[index];
584  int *bc=syzstr->backcomponents[index];
585  int *F1=syzstr->Firstelem[index-1];
586  int *H1=syzstr->Howmuch[index-1];
587  poly pp;
588  polyset o_r=syzstr->orderedRes[index]->m;
589  polyset or1=syzstr->orderedRes[index-1]->m;
590  BOOLEAN ret = FALSE;
591
592  // if != 0, then new element can go into same component
593  // i.e., we do not need to leave space in shifted components
594  long same_comp = 0;
595
596  if (p==NULL) return FALSE;
597  if (realcomp==0) realcomp=1;
598
599  if (index>1)
600    tc = trind1[pGetComp(p)]-1;
601  else
602    tc = pGetComp(p)-1;
603  loop         //while ((j<ie) && (trind1[orc]<=tc+1))
604  {
605    if (j>=ie)
606      break;
607    else
608    {
609      orc = pGetComp(o_r[j]);
610      if (trind1[orc]>tc+1) break;
611      else if (trind1[orc] == tc+1)
612      {
613        same_comp = 1;
614      }
615      else
616      {
617        assume(same_comp == 0);
618      }
619      j += H1[orc];
620    }
621  }
622  if (j>ie)
623  {
624    WerrorS("orderedRes to small");
625    return FALSE;
626  }
627  ie++;
628  if (j == (ie -1))
629  {
630    // new element is the last in ordered module
631    if (same_comp == 0)
632      same_comp = SYZ_SHIFT_BASE;
633
634    // test wheter we have enough space for new shifted component
635    if ((LONG_MAX - same_comp) <= shind[ie-1])
636    {
637      long new_space = syReorderShiftedComponents(shind, ie);
638      assume((LONG_MAX - same_comp) > shind[ie-1]);
639      ret = TRUE;
640      if (TEST_OPT_PROT) Print("(T%u)", new_space);
641    }
642
643    // yes, then set new shifted component
644    assume(ie == 1 || shind[ie-1] > 0);
645    shind[ie] = shind[ie-1] + same_comp;
646  }
647  else
648  {
649    // new element must come in between
650    // i.e. at place j+1
651    long prev, next;
652
653    // test whether new component can get shifted value
654    prev = shind[j];
655    next = shind[j+1];
656    assume(next > prev);
657    if ((same_comp && prev + 2 >= next) || (!same_comp && next - prev < 4))
658    {
659       long new_space = syReorderShiftedComponents(shind, ie);
660      prev = shind[j];
661      next = shind[j+1];
662      assume((same_comp && prev + 2 < next) || (!same_comp && next - prev >= 4));
663      ret = TRUE;
664     if (TEST_OPT_PROT) Print("(B%u)", new_space);
665    }
666
667    // make room for insertion of j+1 shifted component
668    for (k=ie; k > j+1; k--) shind[k] = shind[k-1];
669
670    if (same_comp)
671    {
672      // can simply add one
673      shind[j+1] = prev + 1;
674      assume(shind[j+1] + 1 < shind[j+2]);
675    }
676    else
677    {
678      // need to leave more breathing room - i.e. value goes in
679      // between
680      shind[j+1]  = prev + ((next - prev) >> 1);
681      assume (shind[j] + 1 < shind[j+1] && shind[j+1] + 1 < shind[j+2]);
682    }
683  }
684
685  if (o_r[j]!=NULL)
686  {
687    for (k=ie-1;k>j;k--)
688    {
689      o_r[k] = o_r[k-1];
690      bc[k] = bc[k-1];
691    }
692  }
693  o_r[j] = p;
694  bc[j] = realcomp-1;
695  (H1[pGetComp(p)])++;
696  for (k=0;k<i;k++)
697  {
698    if (F1[k]>j)
699      (F1[k])++;
700  }
701  if (F1[pGetComp(p)]==0)
702    F1[pGetComp(p)]=j+1;
703  for (k=0;k<IDELEMS((syzstr->res)[index]);k++)
704  {
705    if (trind[k]>j)
706      trind[k] += 1;
707  }
708  for (k=IDELEMS((syzstr->res)[index])-1;k>realcomp;k--)
709    trind[k] = trind[k-1];
710  trind[realcomp] = j+1;
711  return ret;
712}
713
714//#define OLD_PAIR_ORDER
715#ifdef OLD_PAIR_ORDER
716static intvec* syLinStrat(SSet nextPairs, syStrategy syzstr,
717                          int howmuch, int index)
718{
719  int i=howmuch-1,i1=0,l,ll;
720  int ** Fin=syzstr->Firstelem;
721  int ** Hin=syzstr->Howmuch;
722  int ** bin=syzstr->backcomponents;
723  ideal o_r=syzstr->orderedRes[index+1];
724  intvec *result=new intvec(howmuch+1);
725  BOOLEAN isDivisible;
726  SObject tso;
727
728  while (i>=0)
729  {
730    tso = nextPairs[i];
731    isDivisible = FALSE;
732    if (syzstr->res[index+1]!=NULL)
733    {
734      l = Fin[index][pGetComp(tso.lcm)]-1;
735      if (l>=0)
736      {
737        ll = l+Hin[index][pGetComp(tso.lcm)];
738        while ((l<ll) && (!isDivisible))
739        {
740          if (o_r->m[l]!=NULL)
741          {
742            isDivisible = isDivisible ||
743              pDivisibleBy2(o_r->m[l],tso.lcm);
744          }
745          l++;
746        }
747      }
748    }
749    if (isDivisible)
750    {
751      syDeletePair(&nextPairs[i]);
752      //crit++;
753    }
754    else
755    {
756      nextPairs[i].p =
757        //sySPoly(tso.p1, tso.p2,tso.lcm);
758        spSpolyCreate(tso.p2, tso.p1,NULL,spSpolyLoop_General);
759      (*result)[i1] = i+1;
760      i1++;
761    }
762    i--;
763  }
764  return result;
765}
766#else
767static intvec* syLinStrat(SSet nextPairs, syStrategy syzstr,
768                          int howmuch, int index)
769{
770  int i=howmuch-1,i1=0,i2,i3,l,ll;
771  int ** Fin=syzstr->Firstelem;
772  int ** Hin=syzstr->Howmuch;
773  int ** bin=syzstr->backcomponents;
774  ideal o_r=syzstr->orderedRes[index+1];
775  intvec *result=new intvec(howmuch+1);
776  intvec *spl=new intvec(howmuch,1,-1);
777  BOOLEAN isDivisible;
778  SObject tso;
779
780  while (i>=0)
781  {
782    tso = nextPairs[i];
783    isDivisible = FALSE;
784    if (syzstr->res[index+1]!=NULL)
785    {
786      l = Fin[index][pGetComp(tso.lcm)]-1;
787      if (l>=0)
788      {
789        ll = l+Hin[index][pGetComp(tso.lcm)];
790        while ((l<ll) && (!isDivisible))
791        {
792          if (o_r->m[l]!=NULL)
793          {
794            isDivisible = isDivisible ||
795              pDivisibleBy2(o_r->m[l],tso.lcm);
796          }
797          l++;
798        }
799      }
800    }
801    if (isDivisible)
802    {
803      syDeletePair(&nextPairs[i]);
804      //crit++;
805    }
806    else
807    {
808      pTest(tso.p2);
809      pTest(tso.p1);
810      nextPairs[i].p =
811        ksOldCreateSpoly(tso.p2, tso.p1,NULL);
812      (*spl)[i] = pLength(nextPairs[i].p);
813    }
814    i--;
815  }
816  i3 = 0;
817  loop
818  {
819    i2 = -1;
820    for (i1=0;i1<howmuch;i1++)
821    {
822      if (i2==-1)
823      {
824        if ((*spl)[i1]!=-1)
825        {
826          i2 = i1;
827        }
828      }
829      else
830      {
831        if (((*spl)[i1]>=0) && ((*spl)[i1]<(*spl)[i2]))
832        {
833          i2 = i1;
834        }
835      }
836    }
837    if (i2>=0)
838    {
839      (*result)[i3] = i2+1;
840      (*spl)[i2] = -1;
841      i3++;
842    }
843    else
844    {
845      break;
846    }
847  }
848  return result;
849}
850#endif
851
852void syEnlargeFields(syStrategy syzstr,int index)
853{
854  pEnlargeSet(&(syzstr->res[index]->m),IDELEMS(syzstr->res[index]),16);
855  syzstr->truecomponents[index]=(int*)ReAlloc0((ADDRESS)syzstr->truecomponents[index],
856                               (IDELEMS(syzstr->res[index])+1)*sizeof(int),
857                               (IDELEMS(syzstr->res[index])+17)*sizeof(int));
858  syzstr->ShiftedComponents[index]
859    =(long*)ReAlloc0((ADDRESS)syzstr->ShiftedComponents[index],
860                              (IDELEMS(syzstr->res[index])+1)*sizeof(long),
861                              (IDELEMS(syzstr->res[index])+17)*sizeof(long));
862  syzstr->backcomponents[index]=(int*)ReAlloc0((ADDRESS)syzstr->backcomponents[index],
863                               (IDELEMS(syzstr->res[index])+1)*sizeof(int),
864                               (IDELEMS(syzstr->res[index])+17)*sizeof(int));
865  syzstr->Howmuch[index]=(int*)ReAlloc0((ADDRESS)syzstr->Howmuch[index],
866                               (IDELEMS(syzstr->res[index])+1)*sizeof(int),
867                               (IDELEMS(syzstr->res[index])+17)*sizeof(int));
868  syzstr->Firstelem[index]=(int*)ReAlloc0((ADDRESS)syzstr->Firstelem[index],
869                               (IDELEMS(syzstr->res[index])+1)*sizeof(int),
870                               (IDELEMS(syzstr->res[index])+17)*sizeof(int));
871  syzstr->elemLength[index]=(int*)ReAlloc0((ADDRESS)syzstr->elemLength[index],
872                               (IDELEMS(syzstr->res[index])+1)*sizeof(int),
873                               (IDELEMS(syzstr->res[index])+17)*sizeof(int));
874//Print("sort %d has now size %d\n",index,IDELEMS(res[index])+17);
875  IDELEMS(syzstr->res[index]) += 16;
876  pEnlargeSet(&(syzstr->orderedRes[index]->m),IDELEMS(syzstr->orderedRes[index]),16);
877  IDELEMS(syzstr->orderedRes[index]) += 16;
878}
879/*3
880* reduces all pairs of degree deg in the module index
881* put the reduced generators to the resolvente which contains
882* the truncated kStd
883*/
884static void syRedNextPairs(SSet nextPairs, syStrategy syzstr,
885               int howmuch, int index)
886{
887  int i,j,k=IDELEMS(syzstr->res[index]);
888  int ks=IDELEMS(syzstr->res[index+1]),kk,l,ll;
889  int * Fin=syzstr->Firstelem[index-1];
890  int * Hin=syzstr->Howmuch[index-1];
891  int * bin=syzstr->backcomponents[index];
892  int * elL=syzstr->elemLength[index];
893  number coefgcd,n;
894  polyset redset=syzstr->orderedRes[index]->m;
895  poly p=NULL,q;
896  intvec *spl1=new intvec(howmuch+1);
897  SObject tso;
898  long * ShiftedComponents = syzstr->ShiftedComponents[index];
899  int* Components = syzstr->truecomponents[index];
900  assume(Components != NULL && ShiftedComponents != NULL);
901  BOOLEAN need_reset;
902
903  if ((nextPairs==NULL) || (howmuch==0)) return;
904  while ((k>0) && (syzstr->res[index]->m[k-1]==NULL)) k--;
905  while ((ks>0) && (syzstr->res[index+1]->m[ks-1]==NULL)) ks--;
906  spl1 = syLinStrat(nextPairs,syzstr,howmuch,index);
907  i=0;
908  while ((*spl1)[i]>0)
909  {
910    need_reset = FALSE;
911    tso = nextPairs[(*spl1)[i]-1];
912    if ((tso.p1!=NULL) && (tso.p2!=NULL))
913    {
914      coefgcd =
915        nGcd(pGetCoeff(tso.p1),pGetCoeff(tso.p2));
916      tso.syz = pHead(tso.lcm);
917      p = tso.syz;
918      pSetCoeff(p,nDiv(pGetCoeff(tso.p1),coefgcd));
919      pGetCoeff(p) = nNeg(pGetCoeff(p));
920      pSetComp(p,tso.ind2+1);
921      rSetmS(p, Components, ShiftedComponents); // actueller index
922      pNext(p) = pHead(tso.lcm);
923      pIter(p);
924      pSetComp(p,tso.ind1+1);
925      rSetmS(p, Components, ShiftedComponents); // actueller index
926      pSetCoeff(p,nDiv(pGetCoeff(tso.p2),coefgcd));
927      nDelete(&coefgcd);
928      if (tso.p != NULL)
929      {
930        kBucketInit(syzstr->bucket,tso.p,-1);
931        q = kBucketGetLm(syzstr->bucket);
932        j = Fin[pGetComp(q)]-1;
933        int pos = j+Hin[pGetComp(q)];
934        loop
935        {
936          if (j<0) break;
937          if (pDivisibleBy2(redset[j],q))
938          {
939            pNext(p) = pHead(q);
940            pIter(p);
941            pSetComp(p,bin[j]+1);
942            rSetmS(p, Components, ShiftedComponents); // actueller index
943//if (pLength(redset[j])!=syzstr->elemLength[index][bin[j]])
944//Print("Halt");
945//if (pLength(redset[j])!=syzstr->elemLength[index][bin[j]])
946//Print("Halt");
947            pGetCoeff(p) = nNeg(pGetCoeff(p));
948            number up = kBucketPolyRed(syzstr->bucket,redset[j],elL[bin[j]],
949                                       NULL);
950            // Thomas: Check whether you need number here
951            nDelete(&up);
952            q = kBucketGetLm(syzstr->bucket);
953            if (q==NULL) break;
954            j = Fin[pGetComp(q)]-1;
955            pos = j+Hin[pGetComp(q)];
956          }
957          else
958          {
959            j++;
960            if (j==pos) break;
961          }
962        }
963        int lb;
964        kBucketClear(syzstr->bucket,&tso.p,&lb);
965      }
966      if (tso.p != NULL)
967      {
968        if (TEST_OPT_PROT) Print("g");
969        if (k==IDELEMS((syzstr->res)[index]))
970        {
971          syEnlargeFields(syzstr,index);
972          bin=syzstr->backcomponents[index];
973          elL=syzstr->elemLength[index];
974          redset=syzstr->orderedRes[index]->m;
975          Components = syzstr->truecomponents[index];
976          ShiftedComponents = syzstr->ShiftedComponents[index];
977        }
978        pNext(p) = pHead(tso.p);
979        pIter(p);
980
981        assume(p!= NULL);
982        k++;
983        syzstr->res[index]->m[k-1] = tso.p;
984        syzstr->elemLength[index][k-1] = pLength(tso.p);
985        pNorm(syzstr->res[index]->m[k-1]);
986        need_reset = syOrder(syzstr->res[index]->m[k-1],syzstr,index,k);
987        pSetComp(p,k); // actueller index
988        rSetmS(p, Components, ShiftedComponents);
989        pGetCoeff(p) = nNeg(pGetCoeff(p));
990
991        tso.isNotMinimal = p;
992        tso.p = NULL;
993      }
994      else
995      {
996        if (TEST_OPT_PROT) Print(".");
997        //if (index % 2==0)
998          //euler++;
999        //else
1000          //euler--;
1001      }
1002      if (ks==IDELEMS(syzstr->res[index+1]))
1003      {
1004        syEnlargeFields(syzstr,index+1);
1005      }
1006      syzstr->res[index+1]->m[ks] = tso.syz;
1007      syzstr->elemLength[index+1][ks] = pLength(tso.syz);
1008      pNorm(syzstr->res[index+1]->m[ks]);
1009      tso.syz =NULL;
1010      tso.syzind = ks;
1011      if (need_reset)
1012        syResetShiftedComponents(syzstr, index+1);
1013      if (syOrder(syzstr->res[index+1]->m[ks],syzstr,index+1,ks+1))
1014         syResetShiftedComponents(syzstr, index+2);
1015      ks++;
1016      p = NULL;
1017      nextPairs[(*spl1)[i]-1] = tso;
1018    }
1019    i++;
1020  }
1021  delete spl1;
1022}
1023
1024/*3
1025* reduces the generators of the module index in degree deg
1026* (which are actual syzygies of the module index-1)
1027* wrt. the ideal generated by elements of lower degrees
1028*/
1029static void syRedGenerOfCurrDeg(syStrategy syzstr, int deg, int index)
1030{
1031  ideal res=syzstr->res[index];
1032  int i=0,j,k=IDELEMS(res),kk;
1033  SSet sPairs=syzstr->resPairs[index-1];
1034
1035  while ((k>0) && (res->m[k-1]==NULL)) k--;
1036  while ((i<(*syzstr->Tl)[index-1]) && (((sPairs)[i].syz==NULL) ||
1037          ((sPairs)[i].order<deg)))
1038    i++;
1039  if ((i>=(*syzstr->Tl)[index-1]) || ((sPairs)[i].order>deg)) return;
1040  while ((i<(*syzstr->Tl)[index-1]) && (((sPairs)[i].syz==NULL) ||
1041         ((sPairs)[i].order==deg)))
1042  {
1043    if ((sPairs)[i].syz!=NULL)
1044    {
1045      j = k-1;
1046      while ((j>=0) && (res->m[j]!=NULL) &&
1047             ((sPairs)[i].syz!=NULL))
1048      {
1049        if (pDivisibleBy1(res->m[j],(sPairs)[i].syz))
1050        {
1051          (sPairs)[i].syz =
1052            ksOldSpolyRed(res->m[j],(sPairs)[i].syz);
1053            //sySPolyRed((sPairs)[i].syz,res->m[j]);
1054          j = k-1;
1055        }
1056        else
1057        {
1058          j--;
1059        }
1060      }
1061      if ((sPairs)[i].syz != NULL)
1062      {
1063        if (k==IDELEMS(res))
1064        {
1065          syEnlargeFields(syzstr,index);
1066          res=syzstr->res[index];
1067        }
1068        if (BTEST1(6))
1069        {
1070          if ((sPairs)[i].isNotMinimal==NULL)
1071          {
1072            PrintLn();
1073            Print("minimal generator: ");pWrite((syzstr->resPairs[index-1])[i].syz);
1074            Print("comes from: ");pWrite((syzstr->resPairs[index-1])[i].p1);
1075            Print("and: ");pWrite((syzstr->resPairs[index-1])[i].p2);
1076          }
1077        }
1078        //res->m[k] = (sPairs)[i].syz;
1079        res->m[k] = syRedtail((sPairs)[i].syz,syzstr,index);
1080        (sPairs)[i].syzind = k;
1081        syzstr->elemLength[index][k] = pLength((sPairs)[i].syz);
1082        pNorm(res->m[k]);
1083  //      (sPairs)[i].syz = NULL;
1084        k++;
1085        if (syOrder(res->m[k-1],syzstr,index,k))
1086          syResetShiftedComponents(syzstr, index);
1087        //euler++;
1088      }
1089      else
1090        (sPairs)[i].syzind = -1;
1091    }
1092    i++;
1093  }
1094}
1095
1096/*3
1097* puts a pair into the right place in resPairs
1098*/
1099static void syEnterPair(SSet sPairs, SObject * so, int * sPlength,int index)
1100{
1101  int ll,k,no=(*so).order,sP=*sPlength,i;
1102  poly p=(*so).lcm;
1103
1104  if ((sP==0) || (sPairs[sP-1].order<=no))
1105    ll = sP;
1106  else if (sP==1)
1107    ll = 0;
1108  else
1109  {
1110    int an=0,en=sP-1;
1111    loop
1112    {
1113      if (an>=en-1)
1114      {
1115        if ((sPairs[an].order<=no) && (sPairs[an+1].order>no))
1116        {
1117          ll = an+1;
1118          break;
1119        }
1120        else if ((sPairs[en].order<=no) && (sPairs[en+1].order>no))
1121        {
1122          ll = en+1;
1123          break;
1124        }
1125        else if (sPairs[an].order>no)
1126        {
1127          ll = an;
1128          break;
1129        }
1130        else
1131        {
1132          Print("Hier ist was faul!\n");
1133          break;
1134        }
1135      }
1136      i=(an+en) / 2;
1137      if (sPairs[i].order <= no)
1138        an=i;
1139      else
1140        en=i;
1141    }
1142  }
1143  for (k=(*sPlength);k>ll;k--)
1144  {
1145    syCopyPair(&sPairs[k-1],&sPairs[k]);
1146  }
1147  syCopyPair(so,&sPairs[ll]);
1148  (*sPlength)++;
1149}
1150void syEnterPair(syStrategy syzstr, SObject * so, int * sPlength,int index)
1151{
1152  int ll;
1153
1154  if (*sPlength>=(*syzstr->Tl)[index])
1155  {
1156    SSet temp = (SSet)Alloc0(((*syzstr->Tl)[index]+16)*sizeof(SObject));
1157    for (ll=0;ll<(*syzstr->Tl)[index];ll++)
1158    {
1159      temp[ll].p = (syzstr->resPairs[index])[ll].p;
1160      temp[ll].p1 = (syzstr->resPairs[index])[ll].p1;
1161      temp[ll].p2 = (syzstr->resPairs[index])[ll].p2;
1162      temp[ll].syz = (syzstr->resPairs[index])[ll].syz;
1163      temp[ll].lcm = (syzstr->resPairs[index])[ll].lcm;
1164      temp[ll].ind1 = (syzstr->resPairs[index])[ll].ind1;
1165      temp[ll].ind2 = (syzstr->resPairs[index])[ll].ind2;
1166      temp[ll].syzind = (syzstr->resPairs[index])[ll].syzind;
1167      temp[ll].order = (syzstr->resPairs[index])[ll].order;
1168      temp[ll].isNotMinimal = (syzstr->resPairs[index])[ll].isNotMinimal;
1169      temp[ll].length = (syzstr->resPairs[index])[ll].length;
1170      temp[ll].reference = (syzstr->resPairs[index])[ll].reference;
1171    }
1172    Free((ADDRESS)syzstr->resPairs[index],(*syzstr->Tl)[index]*sizeof(SObject));
1173    (*syzstr->Tl)[index] += 16;
1174    syzstr->resPairs[index] = temp;
1175  }
1176  syEnterPair(syzstr->resPairs[index],so,sPlength,index);
1177}
1178
1179/*3
1180* computes pairs from the new elements (beginning with the element newEl)
1181* in the module index
1182*/
1183static void syCreateNewPairs(syStrategy syzstr, int index, int newEl)
1184{
1185  SSet temp;
1186  SObject tso;
1187  int i,ii,j,k=IDELEMS(syzstr->res[index]),l=(*syzstr->Tl)[index],ll;
1188  int qc,first,pos,jj,j1;
1189  int * bci=syzstr->backcomponents[index];
1190  poly p,q;
1191  polyset rs=syzstr->res[index]->m,nPm;
1192
1193
1194  while ((k>0) && (rs[k-1]==NULL)) k--;
1195  if (newEl>=k) return;
1196
1197  long * ShiftedComponents = syzstr->ShiftedComponents[index];
1198  int* Components = syzstr->truecomponents[index];
1199
1200  ideal nP=idInit(k,syzstr->res[index]->rank);
1201  nPm=nP->m;
1202  while ((l>0) && ((syzstr->resPairs[index])[l-1].p1==NULL)) l--;
1203  for (j=newEl;j<k;j++)
1204  {
1205    q = rs[j];
1206    qc = pGetComp(q);
1207    first = syzstr->Firstelem[index-1][pGetComp(q)]-1;
1208    pos = first+syzstr->Howmuch[index-1][pGetComp(q)];
1209    for (i=first;i<pos;i++)
1210    {
1211      jj = bci[i];
1212      if (jj>=j) break;
1213      p = pOne();
1214      pLcm(rs[jj],q,p);
1215      pSetComp(p,j+1);
1216      rSetmS(p, Components, ShiftedComponents);
1217      ii = first;
1218      loop
1219      {
1220        j1 = bci[ii];
1221        if (nPm[j1]!=NULL)
1222        {
1223          if (pDivisibleBy2(nPm[j1],p))
1224          {
1225            pDelete(&p);
1226            break;
1227          }
1228          else if (pDivisibleBy2(p,nPm[j1]))
1229          {
1230            pDelete(&(nPm[j1]));
1231            //break;
1232          }
1233        }
1234        ii++;
1235        if (ii>=pos) break;
1236      }
1237      if (p!=NULL)
1238      {
1239        nPm[jj] = p;
1240      }
1241    }
1242    for (i=first;i<pos;i++)
1243    {
1244      ii = bci[i];
1245      if (nPm[ii]!=NULL)
1246      {
1247        if (l>=(*syzstr->Tl)[index])
1248        {
1249          temp = (SSet)Alloc0(((*syzstr->Tl)[index]+16)*sizeof(SObject));
1250          for (ll=0;ll<(*syzstr->Tl)[index];ll++)
1251          {
1252            temp[ll].p = (syzstr->resPairs[index])[ll].p;
1253            temp[ll].p1 = (syzstr->resPairs[index])[ll].p1;
1254            temp[ll].p2 = (syzstr->resPairs[index])[ll].p2;
1255            temp[ll].syz = (syzstr->resPairs[index])[ll].syz;
1256            temp[ll].lcm = (syzstr->resPairs[index])[ll].lcm;
1257            temp[ll].ind1 = (syzstr->resPairs[index])[ll].ind1;
1258            temp[ll].ind2 = (syzstr->resPairs[index])[ll].ind2;
1259            temp[ll].syzind = (syzstr->resPairs[index])[ll].syzind;
1260            temp[ll].order = (syzstr->resPairs[index])[ll].order;
1261            temp[ll].isNotMinimal = (syzstr->resPairs[index])[ll].isNotMinimal;
1262          }
1263          Free((ADDRESS)syzstr->resPairs[index],(*syzstr->Tl)[index]*sizeof(SObject));
1264          (*syzstr->Tl)[index] += 16;
1265          syzstr->resPairs[index] = temp;
1266        }
1267        tso.lcm = p = nPm[ii];
1268        nPm[ii] = NULL;
1269        tso.order = pGetOrder(p) = pTotaldegree(p);
1270        if ((syzstr->cw!=NULL) && (index>0) && (pGetComp(q)>0))
1271        {
1272          int ii=index-1,jj=pGetComp(q);
1273          while (ii>0)
1274          {
1275            jj = pGetComp(syzstr->res[ii]->m[jj-1]);
1276            ii--;
1277          }
1278          tso.order += (*syzstr->cw)[jj-1];
1279        }
1280        tso.p1 = rs[ii];
1281        tso.p2 = q;
1282        tso.ind1 = ii;
1283        tso.ind2 = j;
1284        tso.syzind = -1;
1285        tso.isNotMinimal = NULL;
1286        tso.p = NULL;
1287        tso.syz = NULL;
1288        syEnterPair(syzstr->resPairs[index],&tso,&l,index);
1289      }
1290    }
1291  }
1292  idDelete(&nP);
1293}
1294
1295static SSet syChosePairsPutIn(syStrategy syzstr, int *index,
1296               int *howmuch, int * actdeg, int an, int en)
1297{
1298  int newdeg=*actdeg,newindex=-1,i,t,sldeg;
1299  poly p;
1300  SSet result;
1301  SRes resPairs=syzstr->resPairs;
1302
1303  if (an>syzstr->length) return NULL;
1304  if (en>syzstr->length) en=syzstr->length;
1305  while (*index<en)
1306  {
1307    if (resPairs[*index]!=NULL)
1308    {
1309      sldeg = (*actdeg)+*index;
1310      i = 0;
1311      if (*index!=0)
1312      {
1313        while ((i<(*syzstr->Tl)[*index]))
1314        {
1315          if ((resPairs[*index])[i].lcm!=NULL)
1316          {
1317            if ((resPairs[*index])[i].order == sldeg)
1318            {
1319              result = &(resPairs[*index])[i];
1320              *howmuch =1;
1321              i++;
1322              while ((i<(*syzstr->Tl)[*index]) && ((resPairs[*index])[i].lcm!=NULL)
1323                      && ((resPairs[*index])[i].order == sldeg))
1324              {
1325                i++;
1326                (*howmuch)++;
1327              }
1328              return result;
1329            }
1330          }
1331          i++;
1332        }
1333      }
1334      else
1335      {
1336        while ((i<(*syzstr->Tl)[*index]))
1337        {
1338          if ((resPairs[*index])[i].syz!=NULL)
1339          {
1340            if ((resPairs[*index])[i].order == sldeg)
1341            {
1342              result = &(resPairs[*index])[i];
1343              (*howmuch) =1;
1344              i++;
1345              while ((i<(*syzstr->Tl)[*index]) && ((resPairs[*index])[i].syz!=NULL)
1346                      && ((resPairs[*index])[i].order == *actdeg))
1347              {
1348                i++;
1349                (*howmuch)++;
1350              }
1351              return result;
1352            }
1353          }
1354          i++;
1355        }
1356      }
1357    }
1358    (*index)++;
1359  }
1360  *index = an;
1361  //if (TEST_OPT_PROT) Print("(Euler:%d)",euler);
1362  while (*index<en)
1363  {
1364    if (resPairs[*index]!=NULL)
1365    {
1366      i = 0;
1367      while ((i<(*syzstr->Tl)[*index]))
1368      {
1369        t = *actdeg+*index;
1370        if (((resPairs[*index])[i].lcm!=NULL) ||
1371              ((resPairs[*index])[i].syz!=NULL))
1372        {
1373          if ((resPairs[*index])[i].order > t)
1374            t = (resPairs[*index])[i].order;
1375        }
1376        if ((t>*actdeg+*index) && ((newdeg==*actdeg) || (t<newdeg+*index)))
1377        {
1378          newdeg = t-*index;
1379          newindex = *index;
1380          break;
1381        }
1382        i++;
1383      }
1384    }
1385    (*index)++;
1386  }
1387  if (newdeg>*actdeg)
1388  {
1389    *actdeg = newdeg;
1390    *index = newindex;
1391    return syChosePairsPutIn(syzstr,index,howmuch,actdeg,an,en);
1392  }
1393  else return NULL;
1394}
1395
1396/*3
1397* FOR THE HOMOGENEOUS CASE ONLY!
1398* looks through the pair set and the given module for
1399* remaining pairs or generators to consider
1400* returns a pointer to the first pair and the number of them in the given module
1401* works with slanted degree (i.e. deg=realdeg-index)
1402*/
1403SSet syChosePairs(syStrategy syzstr, int *index, int *howmuch, int * actdeg)
1404{
1405  return syChosePairsPutIn(syzsts,index,howmuch,actdeg,0,syzstr->length);
1406}
1407
1408/*3
1409* FOR THE INHOMOGENEOUS CASE ONLY!
1410* looks through the pair set and the given module for
1411* remaining pairs or generators to consider
1412* returns a pointer to the first pair and the number of them in the given module
1413* works with slanted degree (i.e. deg=realdeg-index)
1414* looks first through the 0 and 1 module then through the other
1415*/
1416static SSet syChosePairsIH(syStrategy syzstr, int *index,
1417               int *howmuch, int * actdeg, int mindeg)
1418{
1419  SSet result=NULL;
1420
1421  result = syChosePairsPutIn(syzstr,index,howmuch,actdeg,0,2);
1422  if (result == NULL)
1423  {
1424    *actdeg = mindeg;
1425    result = syChosePairsPutIn(syzstr,index,howmuch,actdeg,2,syzstr->length);
1426  }
1427  return result;
1428}
1429
1430/*3
1431* looks through the pair set and the given module for
1432* remaining pairs or generators to consider
1433* returns a pointer to the first pair and the number of them in the given module
1434* works deg by deg
1435*/
1436/*
1437*static SSet syChosePairs1(SRes resPairs,intvec * Tl, int *index, int *howmuch,
1438*                   int length,int * actdeg)
1439*{
1440*  int newdeg=*actdeg,newindex=-1,i,t;
1441*  SSet result;
1442*
1443*  while (*index>=0)
1444*  {
1445*    if (resPairs[*index]!=NULL)
1446*    {
1447*      i = 0;
1448*      if (*index!=0)
1449*      {
1450*        while ((i<(*Tl)[*index]))
1451*        {
1452*          if ((resPairs[*index])[i].lcm!=NULL)
1453*          {
1454*            if (pGetOrder((resPairs[*index])[i].lcm) == *actdeg)
1455*            {
1456*              result = &(resPairs[*index])[i];
1457*              *howmuch =1;
1458*              i++;
1459*              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].lcm!=NULL)
1460*                      && (pGetOrder((resPairs[*index])[i].lcm) == *actdeg))
1461*              {
1462*                i++;
1463*                (*howmuch)++;
1464*              }
1465*              return result;
1466*            }
1467*          }
1468*          i++;
1469*        }
1470*      }
1471*      else
1472*      {
1473*        while ((i<(*Tl)[*index]))
1474*        {
1475*          if ((resPairs[*index])[i].syz!=NULL)
1476*          {
1477*            if ((resPairs[*index])[i].order == *actdeg)
1478*            {
1479*              result = &(resPairs[*index])[i];
1480*              (*howmuch) =1;
1481*              i++;
1482*              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].syz!=NULL)
1483*                      && ((resPairs[*index])[i].order == *actdeg))
1484*              {
1485*                i++;
1486*                (*howmuch)++;
1487*              }
1488*              return result;
1489*            }
1490*          }
1491*          i++;
1492*        }
1493*      }
1494*    }
1495*    (*index)--;
1496*  }
1497*  *index = length-1;
1498*  while (*index>=0)
1499*  {
1500*    if (resPairs[*index]!=NULL)
1501*    {
1502*      i = 0;
1503*      while ((i<(*Tl)[*index]))
1504*      {
1505*        t = *actdeg;
1506*        if ((resPairs[*index])[i].lcm!=NULL)
1507*        {
1508*          if (pGetOrder((resPairs[*index])[i].lcm) > *actdeg)
1509*            t = pGetOrder((resPairs[*index])[i].lcm);
1510*        }
1511*        else if ((resPairs[*index])[i].syz!=NULL)
1512*        {
1513*          if ((resPairs[*index])[i].order > *actdeg)
1514*            t = (resPairs[*index])[i].order;
1515*        }
1516*        if ((t>*actdeg) && ((newdeg==*actdeg) || (t<newdeg)))
1517*        {
1518*          newdeg = t;
1519*          newindex = *index;
1520*          break;
1521*        }
1522*        i++;
1523*      }
1524*    }
1525*    (*index)--;
1526*  }
1527*  if (newdeg>*actdeg)
1528*  {
1529*    *actdeg = newdeg;
1530*    *index = newindex;
1531*    return syChosePairs1(resPairs,Tl,index,howmuch,length,actdeg);
1532*  }
1533*  else return NULL;
1534*}
1535*/
1536/*3
1537* statistics of the resolution
1538*/
1539static void syStatistics(resolvente res,int length)
1540{
1541  int i,j=1,k;
1542  Order_t deg = 0;
1543
1544  PrintLn();
1545  while ((j<length) && (res[j]!=NULL))
1546  {
1547    Print("In module %d: \n",j);
1548    k = 0;
1549    while ((k<IDELEMS(res[j])) && (res[j]->m[k]!=NULL))
1550    {
1551      i = 1;
1552      deg = pGetOrder(res[j]->m[k]);
1553      k++;
1554      while ((k<IDELEMS(res[j])) && (res[j]->m[k]!=NULL) &&
1555              (pGetOrder(res[j]->m[k])==deg))
1556      {
1557        i++;
1558        k++;
1559      }
1560      Print("%d elements of degree %d\n",i,deg);
1561    }
1562    j++;
1563  }
1564}
1565
1566/*3
1567* initialize a module
1568*/
1569int syInitSyzMod(syStrategy syzstr, int index, int init)
1570{
1571  int result;
1572
1573  if (syzstr->res[index]==NULL)
1574  {
1575    syzstr->res[index] = idInit(init-1,1);
1576    syzstr->truecomponents[index] = (int*)Alloc0(init*sizeof(int));
1577    syzstr->ShiftedComponents[index] = (long*)Alloc0(init*sizeof(long));
1578    if (index==0)
1579    {
1580      for (int i=0;i<init;i++)
1581      {
1582        syzstr->truecomponents[0][i] = i;
1583        syzstr->ShiftedComponents[0][i] = (i)*SYZ_SHIFT_BASE;
1584      }
1585    }
1586    syzstr->backcomponents[index] = (int*)Alloc0(init*sizeof(int));
1587    syzstr->Howmuch[index] = (int*)Alloc0(init*sizeof(int));
1588    syzstr->Firstelem[index] = (int*)Alloc0(init*sizeof(int));
1589    syzstr->elemLength[index] = (int*)Alloc0(init*sizeof(int));
1590    syzstr->orderedRes[index] = idInit(init-1,1);
1591    result = 0;
1592  }
1593  else
1594  {
1595    result = IDELEMS(syzstr->res[index]);
1596    while ((result>0) && (syzstr->res[index]->m[result-1]==NULL)) result--;
1597  }
1598  return result;
1599}
1600
1601/*3
1602* deletes a resolution
1603*/
1604void syKillComputation(syStrategy syzstr)
1605{
1606  if (syzstr->references>0)
1607  {
1608    (syzstr->references)--;
1609  }
1610  else
1611  {
1612    int i,j;
1613    if (syzstr->minres!=NULL)
1614    {
1615      for (i=0;i<syzstr->length;i++)
1616      {
1617        if (syzstr->minres[i]!=NULL)
1618        {
1619          for (j=0;j<IDELEMS(syzstr->minres[i]);j++)
1620          {
1621            if (syzstr->minres[i]->m[j]!=NULL)
1622              pDelete(&(syzstr->minres[i]->m[j]));
1623          }
1624        }
1625        idDelete(&(syzstr->minres[i]));
1626      }
1627      Free((ADDRESS)syzstr->minres,(syzstr->length+1)*sizeof(ideal));
1628    }
1629    if (syzstr->fullres!=NULL)
1630    {
1631      for (i=0;i<syzstr->length;i++)
1632      {
1633        if (syzstr->fullres[i]!=NULL)
1634        {
1635          for (j=0;j<IDELEMS(syzstr->fullres[i]);j++)
1636          {
1637            if (syzstr->fullres[i]->m[j]!=NULL)
1638              pDelete(&(syzstr->fullres[i]->m[j]));
1639          }
1640        }
1641        idDelete(&(syzstr->fullres[i]));
1642      }
1643      Free((ADDRESS)syzstr->fullres,(syzstr->length+1)*sizeof(ideal));
1644    }
1645    ring origR = currRing;
1646
1647    if (syzstr->syRing != NULL)
1648      rChangeCurrRing(syzstr->syRing, FALSE);
1649
1650    if (syzstr->resPairs!=NULL)
1651    {
1652      for (i=0;i<syzstr->length;i++)
1653      {
1654        for (j=0;j<(*syzstr->Tl)[i];j++)
1655        {
1656          if ((syzstr->resPairs[i])[j].lcm!=NULL)
1657            pDelete(&((syzstr->resPairs[i])[j].lcm));
1658        }
1659        if (syzstr->orderedRes[i]!=NULL)
1660        {
1661          for (j=0;j<IDELEMS(syzstr->orderedRes[i]);j++)
1662          {
1663            syzstr->orderedRes[i]->m[j] = NULL;
1664          }
1665        }
1666        idDelete(&(syzstr->orderedRes[i]));
1667        if (syzstr->truecomponents[i]!=NULL)
1668        {
1669          Free((ADDRESS)syzstr->truecomponents[i],(IDELEMS(syzstr->res[i])+1)*sizeof(int));
1670          syzstr->truecomponents[i]=NULL;
1671          Free((ADDRESS)syzstr->ShiftedComponents[i],(IDELEMS(syzstr->res[i])+1)*sizeof(long));
1672          syzstr->ShiftedComponents[i]=NULL;
1673        }
1674        if (syzstr->backcomponents[i]!=NULL)
1675        {
1676          Free((ADDRESS)syzstr->backcomponents[i],(IDELEMS(syzstr->res[i])+1)*sizeof(int));
1677          syzstr->backcomponents[i]=NULL;
1678        }
1679        if (syzstr->Howmuch[i]!=NULL)
1680        {
1681          Free((ADDRESS)syzstr->Howmuch[i],(IDELEMS(syzstr->res[i])+1)*sizeof(int));
1682          syzstr->Howmuch[i]=NULL;
1683        }
1684        if (syzstr->Firstelem[i]!=NULL)
1685        {
1686          Free((ADDRESS)syzstr->Firstelem[i],(IDELEMS(syzstr->res[i])+1)*sizeof(int));
1687          syzstr->Firstelem[i]=NULL;
1688        }
1689        if (syzstr->res[i]!=NULL)
1690        {
1691          for (j=0;j<IDELEMS(syzstr->res[i]);j++)
1692          {
1693            if (syzstr->res[i]->m[j]!=NULL)
1694              pDelete(&(syzstr->res[i]->m[j]));
1695          }
1696        }
1697        if (syzstr->weights!=0)
1698        {
1699          if (syzstr->weights[i]!=NULL)
1700          {
1701            delete syzstr->weights[i];
1702          }
1703          Free((ADDRESS)syzstr->weights,syzstr->length*sizeof(intvec*));
1704        }
1705        idDelete(&(syzstr->res[i]));
1706        Free((ADDRESS)syzstr->resPairs[i],(*syzstr->Tl)[i]*sizeof(SObject));
1707      }
1708      Free((ADDRESS)syzstr->resPairs,syzstr->length*sizeof(SObject*));
1709      Free((ADDRESS)syzstr->res,(syzstr->length+1)*sizeof(ideal));
1710      Free((ADDRESS)syzstr->orderedRes,(syzstr->length+1)*sizeof(ideal));
1711      Free((ADDRESS)syzstr->elemLength,(syzstr->length+1)*sizeof(int*));
1712      Free((ADDRESS)syzstr->truecomponents,(syzstr->length+1)*sizeof(int*));
1713      Free((ADDRESS)syzstr->ShiftedComponents,(syzstr->length+1)*sizeof(long*));
1714      Free((ADDRESS)syzstr->backcomponents,(syzstr->length+1)*sizeof(int*));
1715      Free((ADDRESS)syzstr->Howmuch,(syzstr->length+1)*sizeof(int*));
1716      Free((ADDRESS)syzstr->Firstelem,(syzstr->length+1)*sizeof(int*));
1717    }
1718    if (syzstr->cw!=NULL)
1719      delete syzstr->cw;
1720    if (syzstr->betti!=NULL)
1721      delete syzstr->betti;
1722    if (syzstr->resolution!=NULL)
1723      delete syzstr->resolution;
1724    if (syzstr->syRing != NULL)
1725    {
1726      rChangeCurrRing(origR, FALSE);
1727      rKill(syzstr->syRing);
1728    }
1729    FreeSizeOf((ADDRESS)syzstr,ssyStrategy);
1730  }
1731}
1732
1733/*3
1734* read out the Betti numbers from resolution
1735* (if not LaScala calls the traditional Betti procedure)
1736*/
1737intvec * syBettiOfComputation(syStrategy syzstr, BOOLEAN minim)
1738{
1739  int dummy;
1740  if (syzstr->betti!=NULL)
1741  {
1742    return ivCopy(syzstr->betti);
1743  }
1744  intvec *result;
1745  if (syzstr->resPairs!=NULL)
1746  {
1747    int i,j=-1,jj=-1,l,sh=0;
1748    SRes rP=syzstr->resPairs;
1749
1750    if (syzstr->hilb_coeffs!=NULL) sh = 1;
1751    l = syzstr->length;
1752    while ((l>0) && (rP[l-1]==NULL)) l--;
1753    if (l==0) return NULL;
1754    l--;
1755    while (l>=sh)
1756    {
1757      i = 0;
1758      while ((i<(*syzstr->Tl)[l]) &&
1759        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)))
1760      {
1761        if (rP[l][i].isNotMinimal==NULL)
1762        {
1763          if (j<rP[l][i].order-l)
1764            j = rP[l][i].order-l;
1765          if (jj<l)
1766            jj = l;
1767        }
1768        i++;
1769      }
1770      l--;
1771    }
1772    j = j+sh;
1773    jj = jj+2;
1774    result=new intvec(j,jj-sh,0);
1775    IMATELEM(*result,1,1) = max(1,idRankFreeModule(syzstr->res[1]));
1776    for (i=sh;i<jj;i++)
1777    {
1778      j = 0;
1779      while (j<(*syzstr->Tl)[i])
1780      {
1781        if ((rP[i][j].isNotMinimal==NULL)&&
1782           ((rP[i][j].lcm!=NULL) || (rP[i][j].syz!=NULL)))
1783          IMATELEM(*result,rP[i][j].order-i+sh,i+2-sh)++;
1784        j++;
1785      }
1786    }
1787    syzstr->betti = result;
1788  }
1789  else if (syzstr->fullres!=NULL)
1790    syzstr->betti = syBetti(syzstr->fullres,syzstr->length,&dummy,NULL,minim);
1791  else
1792    syzstr->betti = syBetti(syzstr->minres,syzstr->length,&dummy,NULL,minim);
1793  return ivCopy(syzstr->betti);
1794}
1795
1796/*2
1797* read out the Betti numbers from resolution
1798* (interpreter interface)
1799*/
1800BOOLEAN syBetti2(leftv res, leftv u, leftv w)
1801{
1802  syStrategy syzstr=(syStrategy)u->Data();
1803  BOOLEAN minim=(int)w->Data();
1804
1805  res->data=(void *)syBettiOfComputation(syzstr,minim);
1806  return FALSE;
1807}
1808BOOLEAN syBetti1(leftv res, leftv u)
1809{
1810  syStrategy syzstr=(syStrategy)u->Data();
1811
1812  res->data=(void *)syBettiOfComputation(syzstr);
1813  return FALSE;
1814}
1815
1816/*3
1817* computes the allocated length of the resolution
1818*/
1819int syLength(syStrategy syzstr)
1820{
1821  return syzstr->length;
1822}
1823
1824/*3
1825* computes the real length of the resolution
1826*/
1827int sySize(syStrategy syzstr)
1828{
1829  resolvente r=syzstr->res;
1830  if (r==NULL)
1831    r = syzstr->fullres;
1832  if (r==NULL)
1833    r = syzstr->minres;
1834  if (r==NULL)
1835  {
1836    WerrorS("No resolution found");
1837    return 0;
1838  }
1839  int i=syzstr->length;
1840  while ((i>0) && (r[i-1]==NULL)) i--;
1841  return i;
1842}
1843
1844/*3
1845* computes the cohomological dimension of res[1]
1846*/
1847int syDim(syStrategy syzstr)
1848{
1849  int i,j=-1,l;
1850  if (syzstr->resPairs!=NULL)
1851  {
1852    SRes rP=syzstr->resPairs;
1853
1854    l = syzstr->length;
1855    while ((l>0) && (rP[l-1]==NULL)) l--;
1856    if (l==0) return -1;
1857    l--;
1858    while (l>=0)
1859    {
1860      i = 0;
1861      while ((i<(*syzstr->Tl)[l]) &&
1862        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)) &&
1863        (rP[l][i].isNotMinimal!=NULL))
1864      {
1865        i++;
1866      }
1867      if ((i<(*syzstr->Tl)[l]) &&
1868        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)) &&
1869        (rP[l][i].isNotMinimal==NULL))
1870        return l;
1871      l--;
1872    }
1873    return l;
1874  }
1875  else
1876    return sySize(syzstr);
1877}
1878
1879/*3
1880* copies the resolution (by increment the reference counter)
1881*/
1882syStrategy syCopy(syStrategy syzstr)
1883{
1884  syStrategy result=syzstr;
1885  (result->references)++;
1886  return result;
1887}
1888
1889/*2
1890* local print procedure used in syPrint
1891*/
1892static void syPrintEmptySpaces(int i)
1893{
1894  if (i!=0)
1895  {
1896    Print(" ");
1897    syPrintEmptySpaces(i/10);
1898  }
1899}
1900
1901/*2
1902* local print procedure used in syPrint
1903*/
1904static void syPrintEmptySpaces1(int i)
1905{
1906  if (i!=0)
1907  {
1908    Print(" ");
1909    syPrintEmptySpaces1(i-1);
1910  }
1911}
1912
1913/*2
1914* local print procedure used in syPrint
1915*/
1916static int syLengthInt(int i)
1917{
1918  int j=0;
1919
1920  if (i==0) return 1;
1921  while (i!=0)
1922  {
1923    j++;
1924    i = i/10;
1925  }
1926  return j;
1927}
1928
1929/*3
1930* prints the resolution as sequence of free modules
1931*/
1932void syPrint(syStrategy syzstr)
1933{
1934  if ((syzstr->resPairs==NULL) && (syzstr->fullres==NULL)
1935     && (syzstr->minres==NULL))
1936  {
1937    Print("No resolution defined\n");
1938    return;
1939  }
1940  int l=0;
1941  if (syzstr->resolution==NULL)
1942  {
1943    int j;
1944    if (syzstr->resPairs!=NULL)
1945    {
1946      syzstr->resolution = new intvec(syzstr->length+1);
1947      SRes rP=syzstr->resPairs;
1948      (*syzstr->resolution)[0] = max(1,idRankFreeModule(syzstr->res[1]));
1949      while ((l<syzstr->length) && (rP[l]!=NULL))
1950      {
1951        j=0;
1952        while ((j<(*syzstr->Tl)[l]) &&
1953          ((rP[l][j].lcm!=NULL) || (rP[l][j].syz!=NULL)))
1954        {
1955          if (rP[l][j].isNotMinimal==NULL)
1956            ((*syzstr->resolution)[l+1])++;
1957          j++;
1958        }
1959        l++;
1960      }
1961    }
1962    else
1963    {
1964      resolvente rr;
1965      syzstr->resolution = new intvec(syzstr->length+2);
1966      if (syzstr->minres!=NULL)
1967        rr = syzstr->minres;
1968      else
1969        rr = syzstr->fullres;
1970      (*syzstr->resolution)[0] = max(1,idRankFreeModule(rr[0]));
1971      while ((l<syzstr->length) && (rr[l]!=NULL))
1972      {
1973        j = IDELEMS(rr[l]);
1974        while ((l>0) && (rr[l]->m[j-1]==NULL)) j--;
1975        ((*syzstr->resolution)[l+1]) = j;
1976        l++;
1977      }
1978    }
1979  }
1980  char *sn=currRingHdl->id;
1981  int sl=strlen(sn);
1982  syPrintEmptySpaces1(sl);
1983  l = 0;
1984  loop
1985  {
1986    Print("%d",(*syzstr->resolution)[l]);
1987    if ((l>=syzstr->resolution->length()) || ((*syzstr->resolution)[l]==0))
1988      break;
1989    syPrintEmptySpaces1(sl+5);
1990    l++;
1991  }
1992  PrintLn();
1993  l = 0;
1994  loop
1995  {
1996    Print(sn);
1997    if ((l>=syzstr->resolution->length()) || ((*syzstr->resolution)[l]==0))
1998      break;
1999    syPrintEmptySpaces((*syzstr->resolution)[l]);
2000    Print(" <-- ");
2001    l++;
2002  }
2003  PrintLn();
2004  PrintLn();
2005  l = 0;
2006  loop
2007  {
2008    Print("%d",l);
2009    if ((l>=syzstr->resolution->length()) || ((*syzstr->resolution)[l]==0))
2010      break;
2011    syPrintEmptySpaces1(sl+5+syLengthInt((*syzstr->resolution)[l])-
2012                         syLengthInt(l));
2013    l++;
2014  }
2015  PrintLn();
2016  if (syzstr->minres==NULL)
2017  {
2018    Print("resolution not minimized yet");
2019    PrintLn();
2020  }
2021}
2022
2023/*2
2024* divides out the weight monomials (given by the Schreyer-ordering)
2025* from the LaScala-resolution
2026*/
2027static resolvente syReorder(resolvente res,int length,
2028        syStrategy syzstr,BOOLEAN toCopy=TRUE,resolvente totake=NULL)
2029{
2030  int i,j,l;
2031  poly p,q,tq;
2032  polyset ri1;
2033  resolvente fullres;
2034  ring origR=syzstr->syRing;
2035  fullres = (resolvente)Alloc0((length+1)*sizeof(ideal));
2036  if (totake==NULL)
2037    totake = res;
2038  for (i=length-1;i>0;i--)
2039  {
2040    if (res[i]!=NULL)
2041    {
2042      if (i>1)
2043      {
2044        j = IDELEMS(res[i-1]);
2045        while ((j>0) && (res[i-1]->m[j-1]==NULL)) j--;
2046        fullres[i-1] = idInit(IDELEMS(res[i]),j);
2047        ri1 = totake[i-1]->m;
2048        for (j=IDELEMS(res[i])-1;j>=0;j--)
2049        {
2050          p = res[i]->m[j];
2051          q = NULL;
2052          while (p!=NULL)
2053          {
2054            if (toCopy)
2055            {
2056              if (origR!=NULL)
2057                tq = pFetchHead(origR,p);
2058              else
2059                tq = pHead(p);
2060              pIter(p);
2061            }
2062            else
2063            {
2064              res[i]->m[j] = NULL;
2065              if (origR!=NULL)
2066              {
2067                poly pp=p;
2068                pIter(p);
2069                pNext(pp)=NULL;
2070                tq = pFetchCopyDelete(origR,pp);
2071              }
2072              else
2073              {
2074                tq = p;
2075                pIter(p);
2076                pNext(tq) = NULL;
2077              }
2078            }
2079            for (l=pVariables;l>0;l--)
2080            {
2081              if (origR!=NULL)
2082                pSubExp(tq,l, pRingGetExp(origR,ri1[pGetComp(tq)-1],l));
2083              else
2084                pSubExp(tq,l, pGetExp(ri1[pGetComp(tq)-1],l));
2085            }
2086            pSetm(tq);
2087            q = pAdd(q,tq);
2088          }
2089          fullres[i-1]->m[j] = q;
2090        }
2091      }
2092      else
2093      {
2094        if (origR!=NULL)
2095        {
2096          fullres[i-1] = idInit(IDELEMS(res[i]),res[i]->rank);
2097          for (j=IDELEMS(res[i])-1;j>=0;j--)
2098          {
2099            if (toCopy)
2100              fullres[i-1]->m[j] = pFetchCopy(origR,res[i]->m[j]);
2101            else
2102            {
2103              fullres[i-1]->m[j] = pFetchCopyDelete(origR,res[i]->m[j]);
2104              res[i]->m[j] = NULL;
2105            }
2106          }
2107        }
2108        else
2109        {
2110          if (toCopy)
2111            fullres[i-1] = idCopy(res[i]);
2112          else
2113          {
2114            fullres[i-1] = res[i];
2115            res[i] = NULL;
2116          }
2117        }
2118        for (j=IDELEMS(fullres[i-1])-1;j>=0;j--)
2119          fullres[i-1]->m[j] = pOrdPolyInsertSetm(fullres[i-1]->m[j]);
2120      }
2121      if (!toCopy)
2122      {
2123        if (res[i]!=NULL) idDelete(&res[i]);
2124      }
2125    }
2126  }
2127  if (!toCopy)
2128    Free((ADDRESS)res,(length+1)*sizeof(ideal));
2129  //syzstr->length = length;
2130  return fullres;
2131}
2132
2133/*3
2134* converts a resolution into a list of modules
2135*/
2136lists syConvRes(syStrategy syzstr,BOOLEAN toDel)
2137{
2138  if ((syzstr->fullres==NULL) && (syzstr->minres==NULL))
2139  {
2140    if (syzstr->hilb_coeffs==NULL)
2141    {
2142      syzstr->fullres = syReorder(syzstr->res,syzstr->length,syzstr);
2143    }
2144    else
2145    {
2146      syzstr->minres = syReorder(syzstr->orderedRes,syzstr->length,syzstr);
2147      syKillEmptyEntres(syzstr->minres,syzstr->length);
2148    }
2149  }
2150  resolvente tr;
2151  int typ0=IDEAL_CMD;
2152  if (syzstr->minres!=NULL)
2153    tr = syzstr->minres;
2154  else
2155    tr = syzstr->fullres;
2156  resolvente trueres=NULL;
2157  intvec ** w=NULL;
2158  if (syzstr->length>0)
2159  {
2160    trueres=(resolvente)Alloc0((syzstr->length)*sizeof(ideal));
2161    for (int i=(syzstr->length)-1;i>=0;i--)
2162    {
2163      if (tr[i]!=NULL)
2164      {
2165        trueres[i] = idCopy(tr[i]);
2166      }
2167    }
2168    if (idRankFreeModule(trueres[0]) > 0)
2169      typ0 = MODUL_CMD;
2170    if (syzstr->weights!=NULL)
2171    {
2172      w = (intvec**)Alloc0((syzstr->length)*sizeof(intvec*));
2173      for (int i=(syzstr->length)-1;i>=0;i--)
2174      {
2175        if (syzstr->weights[i]!=NULL) w[i] = ivCopy(syzstr->weights[i]);
2176      }
2177    }
2178  }
2179  lists li = liMakeResolv(trueres,syzstr->length,syzstr->list_length,typ0,w);
2180  if (toDel) syKillComputation(syzstr);
2181  return li;
2182}
2183
2184/*3
2185* converts a list of modules into a resolution
2186*/
2187syStrategy syConvList(lists li,BOOLEAN toDel)
2188{
2189  int typ0;
2190  syStrategy result=(syStrategy)Alloc0SizeOf(ssyStrategy);
2191
2192  resolvente fr = liFindRes(li,&(result->length),&typ0,&(result->weights));
2193  if (fr != NULL)
2194  {
2195
2196    result->fullres = (resolvente)Alloc0((result->length+1)*sizeof(ideal));
2197    for (int i=result->length-1;i>=0;i--)
2198    {
2199      if (fr[i]!=NULL)
2200        result->fullres[i] = idCopy(fr[i]);
2201    }
2202    result->list_length=result->length;
2203    Free((ADDRESS)fr,(result->length)*sizeof(ideal));
2204  }
2205  else
2206  {
2207    FreeSizeOf (result, ssyStrategy);
2208    result = NULL;
2209  }
2210  if (toDel) li->Clean();
2211  return result;
2212}
2213
2214/*3
2215* converts a list of modules into a minimal resolution
2216*/
2217syStrategy syForceMin(lists li)
2218{
2219  int typ0;
2220  syStrategy result=(syStrategy)Alloc0SizeOf(ssyStrategy);
2221
2222  resolvente fr = liFindRes(li,&(result->length),&typ0);
2223  result->minres = (resolvente)Alloc0((result->length+1)*sizeof(ideal));
2224  for (int i=result->length-1;i>=0;i--)
2225  {
2226    if (fr[i]!=NULL)
2227      result->minres[i] = idCopy(fr[i]);
2228  }
2229  Free((ADDRESS)fr,(result->length)*sizeof(ideal));
2230  return result;
2231}
2232
2233/*2
2234* deleting all monomials the component of which correspond
2235* to non-minimal generators
2236*/
2237static poly syStripOut(poly p,intvec * toStrip)
2238{
2239  if (toStrip==NULL) return p;
2240  poly pp=p;
2241
2242  while ((pp!=NULL) && ((*toStrip)[pGetComp(pp)]!=0))
2243    pDelete1(&pp);
2244  p = pp;
2245  if (pp!=NULL)
2246  {
2247    while (pNext(pp)!=NULL)
2248    {
2249      if ((*toStrip)[pGetComp(pNext(pp))]!=0)
2250        pDelete1(&pNext(pp));
2251      else
2252        pIter(pp);
2253    }
2254  }
2255  return p;
2256}
2257
2258/*2
2259* copies only those monomials the component of which correspond
2260* to minimal generators
2261*/
2262static poly syStripOutCopy(poly p,intvec * toStrip)
2263{
2264  if (toStrip==NULL) return pCopy(p);
2265  poly result=NULL,pp;
2266
2267  while (p!=NULL)
2268  {
2269    if ((*toStrip)[pGetComp(p)]==0)
2270    {
2271      if (result==NULL)
2272      {
2273        result = pp = pHead(p);
2274      }
2275      else
2276      {
2277        pNext(pp) = pHead(p);
2278        pIter(pp);
2279      }
2280    }
2281    pIter(p);
2282  }
2283  return result;
2284}
2285
2286/*2
2287* minimizes toMin
2288*/
2289static poly syMinimizeP(int toMin,syStrategy syzstr,intvec * ordn,int index,
2290                        intvec * toStrip)
2291{
2292  int ii=0,i,j,tc;
2293  poly p,pp,q=NULL,tq,pisN;
2294  SSet sPairs=syzstr->resPairs[index];
2295  poly tempStripped=NULL;
2296
2297  //pp=pCopy(syzstr->res[index+1]->m[toMin]);
2298  pp = syStripOutCopy(syzstr->res[index+1]->m[toMin],toStrip);
2299  while ((ii<ordn->length()) && ((*ordn)[ii]!=-1) &&
2300             (sPairs[(*ordn)[ii]].syzind!=toMin))
2301  {
2302    ii++;
2303  }
2304  while (ii>=0)
2305  {
2306    i = (*ordn)[ii];
2307    if (sPairs[i].isNotMinimal!=NULL)
2308    {
2309      tempStripped =
2310        syStripOutCopy(syzstr->res[index+1]->m[sPairs[i].syzind],toStrip);
2311      pisN = sPairs[i].isNotMinimal;
2312      tc = pGetComp(pisN);
2313      p = pp;
2314      while (p!=NULL)
2315      {
2316        if (pGetComp(p)==tc)
2317        {
2318          tq = pInit();
2319          for(j=pVariables; j>0; j--)
2320            pSetExp(tq,j, pGetExp(p,j)-pGetExp(pisN,j));
2321          pSetComp(tq, 0);
2322          pSetCoeff0(tq,nDiv(pGetCoeff(p),pGetCoeff(pisN)));
2323          pGetCoeff(tq) = nNeg(pGetCoeff(tq));
2324          pSetm(tq);
2325          q = pAdd(q,pMultT(pCopy(tempStripped),tq));
2326          pDelete(&tq);
2327        }
2328        pIter(p);
2329      }
2330      if (q!=NULL)
2331      {
2332        pp = pAdd(pp,q);
2333        q = NULL;
2334      }
2335      pDelete(&tempStripped);
2336    }
2337    ii--;
2338  }
2339  return pp;
2340}
2341
2342/*2
2343* minimizes toMin
2344*/
2345static poly syMinimizeP1(int toMin,syStrategy syzstr,intvec * ordn,int index,
2346                        intvec * toStrip)
2347{
2348  int ii=0,i,j,tc,lp,ltS=-1;
2349  poly p,mp=NULL,pp,q=NULL,tq,pisN;
2350  SSet sPairs=syzstr->resPairs[index];
2351  poly tempStripped=NULL;
2352
2353  pp = syStripOutCopy(syzstr->res[index+1]->m[toMin],toStrip);
2354  kBucketInit(syzstr->bucket,pp,-1);
2355  while ((ii<ordn->length()) && ((*ordn)[ii]!=-1) &&
2356             (sPairs[(*ordn)[ii]].syzind!=toMin))
2357  {
2358    ii++;
2359  }
2360  while (ii>=0)
2361  {
2362    i = (*ordn)[ii];
2363    if (sPairs[i].isNotMinimal!=NULL)
2364    {
2365      tempStripped =
2366        syStripOutCopy(syzstr->res[index+1]->m[sPairs[i].syzind],toStrip);
2367      tc = pGetComp(sPairs[i].isNotMinimal);
2368      //p = pTakeOutComp1(&tempStripped,tc);
2369      int lu;
2370      pTakeOutComp(&tempStripped,tc,&p,&lu);
2371      kBucketTakeOutComp(syzstr->bucket,tc,&mp,&lp);
2372      mp = pDivideM(mp,p);
2373      while (mp!=NULL)
2374      {
2375        p = pNext(mp);
2376        pNext(mp) = NULL;
2377        ltS = -1;
2378        kBucket_Minus_m_Mult_p(syzstr->bucket,mp,tempStripped,&ltS);
2379        mp = p;
2380      }
2381      pDelete(&mp);
2382      pDelete(&tempStripped);
2383    }
2384    ii--;
2385  }
2386  kBucketClear(syzstr->bucket,&pp,&lp);
2387  return pp;
2388}
2389
2390/*2
2391* deletes empty components after minimization
2392*/
2393void syKillEmptyEntres(resolvente res,int length)
2394{
2395  int i,j,jj,k,rj;
2396  intvec * changes;
2397  poly p;
2398  ideal ri;
2399
2400  for (i=0;i<length;i++)
2401  {
2402    ri = res[i];
2403    if (ri!=NULL)
2404    {
2405      rj = IDELEMS(ri);
2406      changes = new intvec(rj+1,1,-1);
2407      while ((rj>0) && (ri->m[rj-1]==NULL)) rj--;
2408      j = k = 0;
2409      while (j+k<rj)
2410      {
2411        if (ri->m[j+k]!=NULL)
2412        {
2413          ri->m[j] = ri->m[j+k];
2414          (*changes)[j+k+1] = j+1;
2415          j++;
2416        }
2417        else
2418        {
2419          k++;
2420        }
2421      }
2422      for (jj=j;jj<rj;jj++)
2423        ri->m[jj] = NULL;
2424      if (res[i+1]!=NULL)
2425      {
2426        ri = res[i+1];
2427        for (j=IDELEMS(ri)-1;j>=0;j--)
2428        {
2429          p = ri->m[j];
2430          while (p!=NULL)
2431          {
2432            pSetComp(p,(*changes)[pGetComp(p)]);
2433            pSetm(p);
2434            pIter(p);
2435          }
2436        }
2437      }
2438    }
2439  }
2440}
2441
2442/*2
2443* determines the components for minimization
2444*/
2445static intvec * syToStrip(syStrategy syzstr, int index)
2446{
2447  intvec * result=NULL;
2448
2449  if ((syzstr->resPairs[index-1]!=NULL) && (!idIs0(syzstr->res[index])))
2450  {
2451    result=new intvec(IDELEMS(syzstr->res[index])+1);
2452    for (int i=(*syzstr->Tl)[index-1]-1;i>=0;i--)
2453    {
2454      if (syzstr->resPairs[index-1][i].isNotMinimal!=NULL)
2455      {
2456        (*result)[syzstr->resPairs[index-1][i].syzind+1] = 1;
2457      }
2458    }
2459  }
2460  return result;
2461}
2462
2463/*2
2464* re-computes the order of pairs during the algorithm
2465* this ensures to procede with a triangular matrix
2466*/
2467static intvec * syOrdPairs(SSet sPairs, int length)
2468{
2469  intvec * result=new intvec(length,1,-1);
2470  int i,j=0,k=-1,l,ii;
2471
2472  loop
2473  {
2474    l = -1;
2475    for(i=0;i<length;i++)
2476    {
2477      if (sPairs[i].syzind>k)
2478      {
2479        if (l==-1)
2480        {
2481          l = sPairs[i].syzind;
2482          ii = i;
2483        }
2484        else
2485        {
2486          if (sPairs[i].syzind<l)
2487          {
2488            l = sPairs[i].syzind;
2489            ii = i;
2490          }
2491        }
2492      }
2493    }
2494    if (l==-1) break;
2495    (*result)[j] = ii;
2496    j++;
2497    k = l;
2498  }
2499  return result;
2500}
2501
2502/*2
2503* minimizes the output of LaScala
2504*/
2505static resolvente syReadOutMinimalRes(syStrategy syzstr,
2506           BOOLEAN computeStd=FALSE)
2507{
2508  intvec * Strip, * ordn;
2509  resolvente tres=(resolvente)Alloc0((syzstr->length+1)*sizeof(ideal));
2510  ring tmpR = NULL;
2511  ring origR = currRing;
2512
2513//Print("Hier ");
2514  if (computeStd)
2515  {
2516    tres[0] = syzstr->res[1];
2517    syzstr->res[1] = idInit(IDELEMS(tres[0]),tres[0]->rank);
2518    return tres;
2519  }
2520  int i,j,l,index,ii,i1;
2521  poly p;
2522  ideal rs;
2523  SSet sPairs;
2524  int * ord,*b0,*b1;
2525  pSetmProc oldSetm=pSetm;
2526  pCompProc oldComp0=pComp0;
2527
2528  assume(syzstr->syRing != NULL);
2529  rChangeCurrRing(syzstr->syRing, TRUE);
2530//Print("laeufts ");
2531  pComp0 = syzcomp2dpc;
2532  syzstr->bucket = kBucketCreate();
2533  for (index=syzstr->length-1;index>0;index--)
2534  {
2535    if (syzstr->resPairs[index]!=NULL)
2536    {
2537//Print("ideal %d: \n",index);
2538      currcomponents = syzstr->truecomponents[index];
2539      currShiftedComponents = syzstr->ShiftedComponents[index];
2540      rChangeSComps(currcomponents, currShiftedComponents,
2541                    IDELEMS(syzstr->res[index]));
2542      sPairs = syzstr->resPairs[index];
2543      Strip = syToStrip(syzstr,index);
2544      tres[index+1] = idInit(IDELEMS(syzstr->res[index+1]),syzstr->res[index+1]->rank);
2545      i1 = (*syzstr->Tl)[index];
2546//Print("i1= %d\n",i1);
2547      ordn = syOrdPairs(sPairs,i1);
2548      for (i=0;i<i1;i++)
2549      {
2550        if ((sPairs[i].isNotMinimal==NULL) && (sPairs[i].lcm!=NULL))
2551        {
2552          l = sPairs[i].syzind;
2553//Print("Minimiere Poly %d: ",l);pWrite(syzstr->res[index+1]->m[l]);
2554          tres[index+1]->m[l] =
2555            syMinimizeP1(l,syzstr,ordn,index,Strip);
2556        }
2557      }
2558      delete Strip;
2559      Strip = NULL;
2560    }
2561  }
2562  currcomponents = syzstr->truecomponents[0];
2563  currShiftedComponents = syzstr->ShiftedComponents[0];
2564  rChangeSComps( currcomponents, currShiftedComponents,
2565                 IDELEMS(syzstr->res[0]));
2566  tres[1] = idInit(IDELEMS(syzstr->res[1]),syzstr->res[1]->rank);
2567  sPairs = syzstr->resPairs[0];
2568  for (i=(*syzstr->Tl)[0]-1;i>=0;i--)
2569  {
2570    if (sPairs[i].syzind>=0)
2571    {
2572      tres[1]->m[sPairs[i].syzind] = pCopy(syzstr->res[1]->m[sPairs[i].syzind]);
2573    }
2574  }
2575/*--- changes to the original ring------------------*/
2576  kBucketDestroy(&syzstr->bucket);
2577  if (syzstr->syRing != NULL)
2578  {
2579    rChangeCurrRing(origR,TRUE);
2580    // Thomas: now make sure that all data which you need is pFetchCopied
2581    // maybe incoporate it into syReorder ??
2582  }
2583  else
2584  {
2585    pSetm=oldSetm;
2586    pComp0=oldComp0;
2587  }
2588  tres = syReorder(tres,syzstr->length,syzstr,FALSE,syzstr->res);
2589  syKillEmptyEntres(tres,syzstr->length);
2590  idSkipZeroes(tres[0]);
2591  return tres;
2592}
2593
2594/*3
2595* minimizes any kind of resolution
2596*/
2597syStrategy syMinimize(syStrategy syzstr)
2598{
2599  if (syzstr->minres==NULL)
2600  {
2601    if (syzstr->resPairs!=NULL)
2602    {
2603      if (syzstr->hilb_coeffs==NULL)
2604      {
2605        syzstr->minres = syReadOutMinimalRes(syzstr);
2606      }
2607      else
2608      {
2609        syzstr->minres = syReorder(syzstr->orderedRes,syzstr->length,syzstr);
2610      }
2611    }
2612    else if (syzstr->fullres!=NULL)
2613    {
2614      syMinimizeResolvente(syzstr->fullres,syzstr->length,1);
2615      syzstr->minres = syzstr->fullres;
2616      syzstr->fullres = NULL;
2617    }
2618  }
2619  (syzstr->references)++;
2620  return syzstr;
2621}
2622
2623/*2
2624* implementation of LaScala's algorithm
2625* assumes that the given module is homogeneous
2626* works with slanted degree, uses syChosePairs
2627*/
2628syStrategy syLaScala3(ideal arg,int * length)
2629{
2630  BOOLEAN noPair=FALSE;
2631  int i,j,* ord,*b0,*b1,actdeg=32000,index=0,reg=-1;
2632  int startdeg,howmuch;
2633  poly p;
2634  ideal temp;
2635  SSet nextPairs;
2636  syStrategy syzstr=(syStrategy)Alloc0SizeOf(ssyStrategy);
2637  ring tmpR = NULL;
2638  ring origR = currRing;
2639
2640  if ((idIs0(arg)) ||
2641      ((idRankFreeModule(arg)>0) && (!idHomModule(arg,NULL,&(syzstr->cw)))))
2642  {
2643    syzstr->minres = (resolvente)Alloc0SizeOf(ideal);
2644    syzstr->length = 1;
2645    syzstr->minres[0] = idInit(1,arg->rank);
2646    return syzstr;
2647  }
2648  pSetmProc oldSetm=pSetm;
2649  pCompProc oldComp0=pComp0;
2650
2651  //crit = 0;
2652  //euler = -1;
2653  redpol = pInit();
2654  syzstr->length = *length = pVariables+2;
2655
2656  // construct new ring
2657  tmpR = (ring) Alloc0SizeOf(sip_sring);
2658  tmpR->wvhdl = (int **)Alloc0(3 * sizeof(int *));
2659  ord = (int*)Alloc0(3*sizeof(int));
2660  b0 = (int*)Alloc0(3*sizeof(int));
2661  b1 = (int*)Alloc0(3*sizeof(int));
2662  ord[0] = ringorder_dp;
2663  ord[1] = ringorder_S;
2664  b0[0] = 1;
2665  b1[0] = currRing->N;
2666  tmpR->OrdSgn = 1;
2667  tmpR->N = currRing->N;
2668  tmpR->ch = currRing->ch;
2669  tmpR->order = ord;
2670  tmpR->block0 = b0;
2671  tmpR->block1 = b1;
2672  tmpR->P = currRing->P;
2673
2674  if (currRing->parameter!=NULL)
2675  {
2676    tmpR->minpoly=nCopy(currRing->minpoly);
2677    tmpR->parameter=(char **)Alloc(rPar(currRing)*sizeof(char *));
2678    for(i=0;i<tmpR->P;i++)
2679    {
2680      tmpR->parameter[i]=mstrdup(currRing->parameter[i]);
2681    }
2682  }
2683  tmpR->names   = (char **)Alloc(tmpR->N * sizeof(char *));
2684  for (i=0; i<tmpR->N; i++)
2685  {
2686    tmpR->names[i] = mstrdup(currRing->names[i]);
2687  }
2688  currcomponents = (int*)Alloc0((arg->rank+1)*sizeof(int));
2689  currShiftedComponents = (long*)Alloc0((arg->rank+1)*sizeof(long));
2690  for (i=0;i<=arg->rank;i++)
2691  {
2692    currShiftedComponents[i] = (i)*SYZ_SHIFT_BASE;
2693    currcomponents[i] = i;
2694  }
2695  rComplete(tmpR, 1);
2696  rChangeCurrRing(tmpR, TRUE);
2697  rChangeSComps(currcomponents, currShiftedComponents, arg->rank);
2698  syzstr->syRing = tmpR;
2699  pComp0 = syzcomp2dpc;
2700/*--- initializes the data structures---------------*/
2701  syzstr->Tl = new intvec(*length);
2702  temp = idInit(IDELEMS(arg),arg->rank);
2703  for (i=0;i<IDELEMS(arg);i++)
2704  {
2705    temp->m[i] = pFetchCopy(origR, arg->m[i]);
2706    if (temp->m[i]!=NULL)
2707    {
2708      j = pTotaldegree(temp->m[i]);
2709      if (j<actdeg) actdeg = j;
2710    }
2711  }
2712  idTest(temp);
2713  idSkipZeroes(temp);
2714  idTest(temp);
2715  syzstr->resPairs = syInitRes(temp,length,syzstr->Tl,syzstr->cw);
2716  Free((ADDRESS)currcomponents,(arg->rank+1)*sizeof(int));
2717  Free((ADDRESS)currShiftedComponents,(arg->rank+1)*sizeof(long));
2718  syzstr->res = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2719  syzstr->orderedRes = (resolvente)Alloc0((*length+1)*sizeof(ideal));
2720  syzstr->elemLength = (int**)Alloc0((*length+1)*sizeof(int*));
2721  syzstr->truecomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2722  syzstr->ShiftedComponents = (long**)Alloc0((*length+1)*sizeof(long*));
2723  syzstr->backcomponents = (int**)Alloc0((*length+1)*sizeof(int*));
2724  syzstr->Howmuch = (int**)Alloc0((*length+1)*sizeof(int*));
2725  syzstr->Firstelem = (int**)Alloc0((*length+1)*sizeof(int*));
2726  syzstr->bucket = kBucketCreate();
2727  int len0=idRankFreeModule(arg)+1;
2728  startdeg = actdeg;
2729  nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2730  //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2731/*--- computes the resolution ----------------------*/
2732  while (nextPairs!=NULL)
2733  {
2734    if (TEST_OPT_PROT) Print("%d",actdeg);
2735    if (TEST_OPT_PROT) Print("(m%d)",index);
2736    if (index==0)
2737      i = syInitSyzMod(syzstr,index,len0);
2738    else
2739      i = syInitSyzMod(syzstr,index);
2740    currcomponents = syzstr->truecomponents[max(index-1,0)];
2741    currShiftedComponents = syzstr->ShiftedComponents[max(index-1,0)];
2742    rChangeSComps(currcomponents, currShiftedComponents,
2743                  IDELEMS(syzstr->res[max(index-1,0)]));
2744    j = syInitSyzMod(syzstr,index+1);
2745    if (index>0)
2746    {
2747      syRedNextPairs(nextPairs,syzstr,howmuch,index);
2748      syCompactifyPairSet(syzstr->resPairs[index],(*syzstr->Tl)[index],0);
2749    }
2750    else
2751      syRedGenerOfCurrDeg(syzstr,actdeg,index+1);
2752/*--- creates new pairs -----------------------------*/
2753    syCreateNewPairs(syzstr,index,i);
2754    if (index<(*length)-1)
2755    {
2756      syCreateNewPairs(syzstr,index+1,j);
2757    }
2758    index++;
2759    nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2760    //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2761  }
2762  if (temp!=NULL) idDelete(&temp);
2763  kBucketDestroy(&(syzstr->bucket));
2764  if (ord!=NULL)
2765  {
2766    rChangeCurrRing(origR,TRUE);
2767  }
2768  else
2769  {
2770    pSetm=oldSetm;
2771    pComp0=oldComp0;
2772  }
2773  pDelete1(&redpol);
2774  if (TEST_OPT_PROT) PrintLn();
2775  return syzstr;
2776}
2777
Note: See TracBrowser for help on using the repository browser.