source: git/kernel/syz1.cc @ 0de0509

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