source: git/kernel/GBEngine/syz1.cc @ 922890

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