source: git/kernel/syz1.cc @ 762407

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