source: git/kernel/syz1.cc @ 83192d

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