source: git/kernel/syz1.cc @ e585b6

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