source: git/kernel/syz1.cc @ 22579cf

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