source: git/Singular/syz1.cc @ 6e56de

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