source: git/kernel/syz1.cc @ cdec33

spielwiese
Last change on this file since cdec33 was fd1b0f2, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: 'hres' uses syRing and returns the result in re/orderedRes... From: Oleksandr Motsak <motsak@mathematik.uni-kl.de> git-svn-id: file:///usr/local/Singular/svn/trunk@14036 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 67.8 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      if(syzstr->syRing->typ[1].ord_typ == ro_syzcomp)
1643        rNChangeSComps(NULL, NULL, syzstr->syRing);
1644     
1645      rKill(syzstr->syRing);
1646    }
1647    omFreeSize((ADDRESS)syzstr, sizeof(ssyStrategy));
1648  }
1649}
1650
1651/*2
1652* divides out the weight monomials (given by the Schreyer-ordering)
1653* from the LaScala-resolution
1654*/
1655resolvente syReorder(resolvente res,int length,
1656        syStrategy syzstr,BOOLEAN toCopy,resolvente totake)
1657{
1658  int i,j,l;
1659  poly p,q,tq;
1660  polyset ri1;
1661  resolvente fullres;
1662  ring origR=syzstr->syRing;
1663  fullres = (resolvente)omAlloc0((length+1)*sizeof(ideal));
1664  if (totake==NULL)
1665    totake = res;
1666  for (i=length-1;i>0;i--)
1667  {
1668    if (res[i]!=NULL)
1669    {
1670      if (i>1)
1671      {
1672        j = IDELEMS(res[i-1]);
1673        while ((j>0) && (res[i-1]->m[j-1]==NULL)) j--;
1674        fullres[i-1] = idInit(IDELEMS(res[i]),j);
1675        ri1 = totake[i-1]->m;
1676        for (j=IDELEMS(res[i])-1;j>=0;j--)
1677        {
1678          p = res[i]->m[j];
1679          q = NULL;
1680          while (p!=NULL)
1681          {
1682            if (toCopy)
1683            {
1684              if (origR!=NULL)
1685                tq = prHeadR(p,origR);
1686              else
1687                tq = pHead(p);
1688              pIter(p);
1689            }
1690            else
1691            {
1692              res[i]->m[j] = NULL;
1693              if (origR!=NULL)
1694              {
1695                poly pp=p;
1696                pIter(p);
1697                pNext(pp)=NULL;
1698                tq = prMoveR(pp, origR);
1699              }
1700              else
1701              {
1702                tq = p;
1703                pIter(p);
1704                pNext(tq) = NULL;
1705              }
1706            }
1707//            pWrite(tq);
1708            pTest(tq);
1709            for (l=pVariables;l>0;l--)
1710            {
1711              if (origR!=NULL)
1712                pSubExp(tq,l, p_GetExp(ri1[pGetComp(tq)-1],l,origR));
1713              else
1714                pSubExp(tq,l, pGetExp(ri1[pGetComp(tq)-1],l));
1715            }
1716            pSetm(tq);
1717            pTest(tq);
1718            q = pAdd(q,tq);
1719            pTest(q);
1720          }
1721          fullres[i-1]->m[j] = q;
1722        }
1723      }
1724      else
1725      {
1726        if (origR!=NULL)
1727        {
1728          fullres[i-1] = idInit(IDELEMS(res[i]),res[i]->rank);
1729          for (j=IDELEMS(res[i])-1;j>=0;j--)
1730          {
1731            if (toCopy)
1732              fullres[i-1]->m[j] = prCopyR(res[i]->m[j], origR);
1733            else
1734            {
1735              fullres[i-1]->m[j] = prMoveR(res[i]->m[j], origR);
1736              res[i]->m[j] = NULL;
1737            }
1738          }
1739        }
1740        else
1741        {
1742          if (toCopy)
1743            fullres[i-1] = idCopy(res[i]);
1744          else
1745          {
1746            fullres[i-1] = res[i];
1747            res[i] = NULL;
1748          }
1749        }
1750        for (j=IDELEMS(fullres[i-1])-1;j>=0;j--)
1751          fullres[i-1]->m[j] = pSortCompCorrect(fullres[i-1]->m[j]);
1752      }
1753      if (!toCopy)
1754      {
1755        if (res[i]!=NULL) idDelete(&res[i]);
1756      }
1757    }
1758  }
1759  if (!toCopy)
1760    omFreeSize((ADDRESS)res,(length+1)*sizeof(ideal));
1761  //syzstr->length = length;
1762  return fullres;
1763}
1764
1765/*3
1766* read out the Betti numbers from resolution
1767* (if not LaScala calls the traditional Betti procedure)
1768*/
1769intvec * syBettiOfComputation(syStrategy syzstr, BOOLEAN minim,int * row_shift,
1770                              intvec* weights)
1771{
1772  int dummy;
1773  BOOLEAN std_weights=TRUE;
1774  if ((weights!=NULL)
1775  && (syzstr->betti!=NULL)
1776  && (syzstr->weights!=NULL) && (syzstr->weights[0]!=NULL))
1777  {
1778    int i;
1779    for(i=weights->length()-1; i>=0; i--)
1780    {
1781      //Print("test %d: %d - %d\n",i,(*weights)[i], (*(syzstr->weights[0]))[i]);
1782      if ((*weights)[i]!=(*(syzstr->weights[0]))[i])
1783      {
1784        std_weights=FALSE;
1785        break;
1786      }
1787    }
1788  }
1789  if ((syzstr->betti!=NULL)
1790  && (std_weights))
1791  {
1792    if (minim || (syzstr->resPairs!=NULL))
1793      return ivCopy(syzstr->betti);
1794  }
1795
1796  resolvente fullres = syzstr->fullres;
1797  resolvente minres = syzstr->minres;
1798  const int length = syzstr->length;
1799 
1800  if ((fullres==NULL) && (minres==NULL))
1801  {
1802     if (syzstr->hilb_coeffs==NULL)
1803     {
1804        fullres = syReorder(syzstr->res, length, syzstr);
1805     }
1806     else
1807     {
1808        minres = syReorder(syzstr->orderedRes, length, syzstr);
1809        syKillEmptyEntres(minres, length);
1810     }
1811  }
1812 
1813  intvec *result=NULL;
1814 
1815  if (fullres!=NULL)
1816    result = syBetti(fullres,length,&dummy,weights,minim,row_shift);
1817  else
1818    result = syBetti(minres,length,&dummy,weights,minim,row_shift);
1819
1820
1821  return result; /// Don't change the syzstr???
1822
1823  // TODO: cleanup thses!
1824  if( fullres != NULL && syzstr->fullres == NULL )
1825    syzstr->fullres = fullres;
1826  if( minres != NULL && syzstr->minres == NULL )
1827    syzstr->minres = minres;
1828
1829  if ((result!=NULL)
1830  && ((minim) || (syzstr->resPairs!=NULL))
1831  && std_weights)
1832  {
1833    syzstr->betti = ivCopy(result); // cache the result...
1834  }
1835 
1836  return result;
1837}
1838
1839/*3
1840* computes the real length of the resolution
1841*/
1842int sySize(syStrategy syzstr)
1843{
1844  resolvente r=syzstr->res;
1845  if (r==NULL)
1846    r = syzstr->fullres;
1847  if (r==NULL)
1848    r = syzstr->minres;
1849  if (r==NULL)
1850  {
1851    WerrorS("No resolution found");
1852    return 0;
1853  }
1854  int i=syzstr->length;
1855  while ((i>0) && (r[i-1]==NULL)) i--;
1856  return i;
1857}
1858
1859/*3
1860* computes the cohomological dimension of res[1]
1861*/
1862int syDim(syStrategy syzstr)
1863{
1864  int i,j=-1,l;
1865  if (syzstr->resPairs!=NULL)
1866  {
1867    SRes rP=syzstr->resPairs;
1868
1869    l = syzstr->length;
1870    while ((l>0) && (rP[l-1]==NULL)) l--;
1871    if (l==0) return -1;
1872    l--;
1873    while (l>=0)
1874    {
1875      i = 0;
1876      while ((i<(*syzstr->Tl)[l]) &&
1877        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)) &&
1878        (rP[l][i].isNotMinimal!=NULL))
1879      {
1880        i++;
1881      }
1882      if ((i<(*syzstr->Tl)[l]) &&
1883        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)) &&
1884        (rP[l][i].isNotMinimal==NULL))
1885        return l;
1886      l--;
1887    }
1888    return l;
1889  }
1890  else
1891    return sySize(syzstr);
1892}
1893
1894/*3
1895* copies the resolution (by increment the reference counter)
1896*/
1897syStrategy syCopy(syStrategy syzstr)
1898{
1899  syStrategy result=syzstr;
1900  (result->references)++;
1901  return result;
1902}
1903
1904/*2
1905* local print procedure used in syPrint
1906*/
1907static void syPrintEmptySpaces(int i)
1908{
1909  if (i!=0)
1910  {
1911    PrintS(" ");
1912    syPrintEmptySpaces(i/10);
1913  }
1914}
1915
1916/*2
1917* local print procedure used in syPrint
1918*/
1919static void syPrintEmptySpaces1(int i)
1920{
1921  if (i!=0)
1922  {
1923    PrintS(" ");
1924    syPrintEmptySpaces1(i-1);
1925  }
1926}
1927
1928/*2
1929* local print procedure used in syPrint
1930*/
1931static int syLengthInt(int i)
1932{
1933  int j=0;
1934
1935  if (i==0) return 1;
1936  while (i!=0)
1937  {
1938    j++;
1939    i = i/10;
1940  }
1941  return j;
1942}
1943
1944/*3
1945* prints the resolution as sequence of free modules
1946*/
1947void syPrint(syStrategy syzstr)
1948{
1949  if ( (syzstr->resPairs==NULL) &&
1950       (syzstr->fullres==NULL)  &&
1951       (syzstr->minres==NULL) &&
1952       (syzstr->resolution == NULL) )
1953  {
1954    PrintS("No resolution defined\n");
1955    return;
1956  }
1957
1958  intvec* resolution = syzstr->resolution;
1959 
1960  if (resolution==NULL)
1961  {
1962    if (syzstr->resPairs!=NULL)
1963    {
1964      resolution = new intvec(syzstr->length+1);
1965      SRes rP = syzstr->resPairs;
1966//      assume(idRankFreeModule(syzstr->res[1], (syzstr->syRing != NULL ? syzstr->syRing : currRing))==syzstr->res[1]->rank);
1967      (*resolution)[0] = syzstr->res[1]->rank;
1968      int k=0;
1969      while ((k<syzstr->length) && (rP[k]!=NULL))
1970      {
1971        int j = 0;
1972        while ((j<(*syzstr->Tl)[k]) &&
1973          ((rP[k][j].lcm!=NULL) || (rP[k][j].syz!=NULL)))
1974        {
1975          if (rP[k][j].isNotMinimal==NULL)
1976            ((*resolution)[k+1])++;
1977          j++;
1978        }
1979        k++;
1980      }
1981    }
1982    else
1983    {
1984      resolution = new intvec(syzstr->length+2);
1985      resolvente rr;
1986      if (syzstr->minres!=NULL)
1987        rr = syzstr->minres;
1988      else
1989        rr = syzstr->fullres;
1990      (*resolution)[0]
1991        = si_max(1,(int)idRankFreeModule(rr[0],
1992                                 (syzstr->syRing != NULL ? syzstr->syRing : currRing)));
1993      int k=0;
1994      while ((k<syzstr->length) && (rr[k]!=NULL))
1995      {
1996        (*resolution)[k+1] = idSize(rr[k]);
1997        k++;
1998      }
1999    }
2000  }
2001  const char *sn=currRingHdl->id;
2002  int sl=strlen(sn);
2003  syPrintEmptySpaces1(sl);
2004  int k = 0;
2005  loop
2006  {
2007    if ((k>=resolution->length()) || ((*resolution)[k]==0))
2008      break;
2009    Print("%d",(*resolution)[k]);
2010    syPrintEmptySpaces1(sl+5);
2011    k++;
2012  }
2013  PrintLn();
2014  k = 0;
2015  loop
2016  {
2017    if ((k>=resolution->length()) || ((*resolution)[k]==0))
2018      break;
2019    PrintS(sn);
2020    if (((k+1)>=resolution->length()) || ((*resolution)[(k+1)]==0))
2021      break;
2022    PrintS(" <-- ");
2023    syPrintEmptySpaces((*resolution)[k]);
2024    k++;
2025  }
2026  PrintLn();
2027  PrintLn();
2028  k = 0;
2029  loop
2030  {
2031    if ((k>=resolution->length()) || ((*resolution)[k]==0))
2032      break;
2033    Print("%d",k);
2034    syPrintEmptySpaces1(sl+5+syLengthInt((*resolution)[k])-
2035                         syLengthInt(k));
2036    k++;
2037  }
2038  PrintLn();
2039  if (syzstr->minres==NULL)
2040  {
2041    PrintS("resolution not minimized yet");
2042    PrintLn();
2043  }
2044
2045  if (syzstr->resolution == NULL) syzstr->resolution = resolution; 
2046}
2047
2048/*2
2049* deleting all monomials the component of which correspond
2050* to non-minimal generators
2051*/
2052static poly syStripOut(poly p,intvec * toStrip)
2053{
2054  if (toStrip==NULL) return p;
2055  poly pp=p;
2056
2057  while ((pp!=NULL) && ((*toStrip)[pGetComp(pp)]!=0))
2058    pLmDelete(&pp);
2059  p = pp;
2060  if (pp!=NULL)
2061  {
2062    while (pNext(pp)!=NULL)
2063    {
2064      if ((*toStrip)[pGetComp(pNext(pp))]!=0)
2065        pLmDelete(&pNext(pp));
2066      else
2067        pIter(pp);
2068    }
2069  }
2070  return p;
2071}
2072
2073/*2
2074* copies only those monomials the component of which correspond
2075* to minimal generators
2076*/
2077static poly syStripOutCopy(poly p,intvec * toStrip)
2078{
2079  if (toStrip==NULL) return pCopy(p);
2080  poly result=NULL,pp;
2081
2082  while (p!=NULL)
2083  {
2084    if ((*toStrip)[pGetComp(p)]==0)
2085    {
2086      if (result==NULL)
2087      {
2088        result = pp = pHead(p);
2089      }
2090      else
2091      {
2092        pNext(pp) = pHead(p);
2093        pIter(pp);
2094      }
2095    }
2096    pIter(p);
2097  }
2098  return result;
2099}
2100
2101/*2
2102* minimizes toMin
2103*/
2104static poly syMinimizeP(int toMin,syStrategy syzstr,intvec * ordn,int index,
2105                        intvec * toStrip)
2106{
2107  int ii=0,i,j,tc;
2108  poly p,pp,q=NULL,tq,pisN;
2109  SSet sPairs=syzstr->resPairs[index];
2110  poly tempStripped=NULL;
2111
2112  //pp=pCopy(syzstr->res[index+1]->m[toMin]);
2113  pp = syStripOutCopy(syzstr->res[index+1]->m[toMin],toStrip);
2114  while ((ii<ordn->length()) && ((*ordn)[ii]!=-1) &&
2115             (sPairs[(*ordn)[ii]].syzind!=toMin))
2116  {
2117    ii++;
2118  }
2119  while (ii>=0)
2120  {
2121    i = (*ordn)[ii];
2122    if (sPairs[i].isNotMinimal!=NULL)
2123    {
2124      tempStripped =
2125        syStripOutCopy(syzstr->res[index+1]->m[sPairs[i].syzind],toStrip);
2126      pisN = sPairs[i].isNotMinimal;
2127      tc = pGetComp(pisN);
2128      p = pp;
2129      while (p!=NULL)
2130      {
2131        if (pGetComp(p)==tc)
2132        {
2133          tq = pInit();
2134          for(j=pVariables; j>0; j--)
2135            pSetExp(tq,j, pGetExp(p,j)-pGetExp(pisN,j));
2136          pSetComp(tq, 0);
2137          pSetCoeff0(tq,nDiv(pGetCoeff(p),pGetCoeff(pisN)));
2138          pGetCoeff(tq) = nNeg(pGetCoeff(tq));
2139          pSetm(tq);
2140          q = pAdd(q,pMult_mm(pCopy(tempStripped),tq));
2141          pDelete(&tq);
2142        }
2143        pIter(p);
2144      }
2145      if (q!=NULL)
2146      {
2147        pp = pAdd(pp,q);
2148        q = NULL;
2149      }
2150      pDelete(&tempStripped);
2151    }
2152    ii--;
2153  }
2154  return pp;
2155}
2156
2157/*2
2158* minimizes toMin
2159*/
2160static poly syMinimizeP1(int toMin,syStrategy syzstr,intvec * ordn,int index,
2161                        intvec * toStrip)
2162{
2163  int ii=0,i,j,tc,lp,ltS=-1;
2164  poly p,mp=NULL,pp,q=NULL,tq,pisN;
2165  SSet sPairs=syzstr->resPairs[index];
2166  poly tempStripped=NULL;
2167
2168  pp = syStripOutCopy(syzstr->res[index+1]->m[toMin],toStrip);
2169  kBucketInit(syzstr->bucket,pp,-1);
2170  while ((ii<ordn->length()) && ((*ordn)[ii]!=-1) &&
2171             (sPairs[(*ordn)[ii]].syzind!=toMin))
2172  {
2173    ii++;
2174  }
2175  while (ii>=0)
2176  {
2177    i = (*ordn)[ii];
2178    if (sPairs[i].isNotMinimal!=NULL)
2179    {
2180      tempStripped =
2181        syStripOutCopy(syzstr->res[index+1]->m[sPairs[i].syzind],toStrip);
2182      tc = pGetComp(sPairs[i].isNotMinimal);
2183      //p = pTakeOutComp1(&tempStripped,tc);
2184      int lu;
2185      pTakeOutComp(&tempStripped,tc,&p,&lu);
2186      kBucketTakeOutComp(syzstr->bucket,tc,&mp,&lp);
2187      mp = pDivideM(mp,p);
2188      while (mp!=NULL)
2189      {
2190        p = pNext(mp);
2191        pNext(mp) = NULL;
2192        ltS = -1;
2193        kBucket_Minus_m_Mult_p(syzstr->bucket,mp,tempStripped,&ltS);
2194        mp = p;
2195      }
2196      pDelete(&mp);
2197      pDelete(&tempStripped);
2198    }
2199    ii--;
2200  }
2201  kBucketClear(syzstr->bucket,&pp,&lp);
2202  return pp;
2203}
2204
2205/*2
2206* deletes empty components after minimization
2207*/
2208void syKillEmptyEntres(resolvente res,int length)
2209{
2210  int i,j,jj,k,rj;
2211  intvec * changes;
2212  poly p;
2213  ideal ri;
2214
2215  for (i=0;i<length;i++)
2216  {
2217    ri = res[i];
2218    if (ri!=NULL)
2219    {
2220      rj = IDELEMS(ri);
2221      changes = new intvec(rj+1,1,-1);
2222      while ((rj>0) && (ri->m[rj-1]==NULL)) rj--;
2223      j = k = 0;
2224      while (j+k<rj)
2225      {
2226        if (ri->m[j+k]!=NULL)
2227        {
2228          ri->m[j] = ri->m[j+k];
2229          (*changes)[j+k+1] = j+1;
2230          j++;
2231        }
2232        else
2233        {
2234          k++;
2235        }
2236      }
2237      for (jj=j;jj<rj;jj++)
2238        ri->m[jj] = NULL;
2239      if (res[i+1]!=NULL)
2240      {
2241        ri = res[i+1];
2242        for (j=IDELEMS(ri)-1;j>=0;j--)
2243        {
2244          p = ri->m[j];
2245          while (p!=NULL)
2246          {
2247            pSetComp(p,(*changes)[pGetComp(p)]);
2248            pSetm(p);
2249            pIter(p);
2250          }
2251        }
2252      }
2253      delete changes;
2254    }
2255  }
2256}
2257
2258/*2
2259* determines the components for minimization
2260*/
2261static intvec * syToStrip(syStrategy syzstr, int index)
2262{
2263  intvec * result=NULL;
2264
2265  if ((syzstr->resPairs[index-1]!=NULL) && (!idIs0(syzstr->res[index])))
2266  {
2267    result=new intvec(IDELEMS(syzstr->res[index])+1);
2268    for (int i=(*syzstr->Tl)[index-1]-1;i>=0;i--)
2269    {
2270      if (syzstr->resPairs[index-1][i].isNotMinimal!=NULL)
2271      {
2272        (*result)[syzstr->resPairs[index-1][i].syzind+1] = 1;
2273      }
2274    }
2275  }
2276  return result;
2277}
2278
2279/*2
2280* re-computes the order of pairs during the algorithm
2281* this ensures to procede with a triangular matrix
2282*/
2283static intvec * syOrdPairs(SSet sPairs, int length)
2284{
2285  intvec * result=new intvec(length,1,-1);
2286  int i,j=0,k=-1,l,ii;
2287
2288  loop
2289  {
2290    l = -1;
2291    for(i=0;i<length;i++)
2292    {
2293      if (sPairs[i].syzind>k)
2294      {
2295        if (l==-1)
2296        {
2297          l = sPairs[i].syzind;
2298          ii = i;
2299        }
2300        else
2301        {
2302          if (sPairs[i].syzind<l)
2303          {
2304            l = sPairs[i].syzind;
2305            ii = i;
2306          }
2307        }
2308      }
2309    }
2310    if (l==-1) break;
2311    (*result)[j] = ii;
2312    j++;
2313    k = l;
2314  }
2315  return result;
2316}
2317
2318/*2
2319* minimizes the output of LaScala
2320*/
2321static resolvente syReadOutMinimalRes(syStrategy syzstr,
2322           BOOLEAN computeStd=FALSE)
2323{
2324  intvec * Strip, * ordn;
2325  resolvente tres=(resolvente)omAlloc0((syzstr->length+1)*sizeof(ideal));
2326  ring tmpR = NULL;
2327  ring origR = currRing;
2328
2329//Print("Hier ");
2330  if (computeStd)
2331  {
2332    tres[0] = syzstr->res[1];
2333    syzstr->res[1] = idInit(IDELEMS(tres[0]),tres[0]->rank);
2334    return tres;
2335  }
2336  int i,j,l,index,ii,i1;
2337  poly p;
2338  ideal rs;
2339  SSet sPairs;
2340  int * ord,*b0,*b1;
2341
2342  assume(syzstr->syRing != NULL);
2343  rChangeCurrRing(syzstr->syRing);
2344//Print("laeufts ");
2345  syzstr->bucket = kBucketCreate();
2346  for (index=syzstr->length-1;index>0;index--)
2347  {
2348    if (syzstr->resPairs[index]!=NULL)
2349    {
2350//Print("ideal %d: \n",index);
2351      currcomponents = syzstr->truecomponents[index];
2352      currShiftedComponents = syzstr->ShiftedComponents[index];
2353      rChangeSComps(currcomponents, currShiftedComponents,
2354                    IDELEMS(syzstr->res[index]));
2355      sPairs = syzstr->resPairs[index];
2356      Strip = syToStrip(syzstr,index);
2357      tres[index+1] = idInit(IDELEMS(syzstr->res[index+1]),syzstr->res[index+1]->rank);
2358      i1 = (*syzstr->Tl)[index];
2359//Print("i1= %d\n",i1);
2360      ordn = syOrdPairs(sPairs,i1);
2361      for (i=0;i<i1;i++)
2362      {
2363        if ((sPairs[i].isNotMinimal==NULL) && (sPairs[i].lcm!=NULL))
2364        {
2365          l = sPairs[i].syzind;
2366//Print("Minimiere Poly %d: ",l);pWrite(syzstr->res[index+1]->m[l]);
2367          tres[index+1]->m[l] =
2368            syMinimizeP1(l,syzstr,ordn,index,Strip);
2369        }
2370      }
2371      delete Strip;
2372      Strip = NULL;
2373    }
2374  }
2375  currcomponents = syzstr->truecomponents[0];
2376  currShiftedComponents = syzstr->ShiftedComponents[0];
2377  rChangeSComps( currcomponents, currShiftedComponents,
2378                 IDELEMS(syzstr->res[0]));
2379  tres[1] = idInit(IDELEMS(syzstr->res[1]),syzstr->res[1]->rank);
2380  sPairs = syzstr->resPairs[0];
2381  for (i=(*syzstr->Tl)[0]-1;i>=0;i--)
2382  {
2383    if (sPairs[i].syzind>=0)
2384    {
2385      tres[1]->m[sPairs[i].syzind] = pCopy(syzstr->res[1]->m[sPairs[i].syzind]);
2386    }
2387  }
2388/*--- changes to the original ring------------------*/
2389  kBucketDestroy(&syzstr->bucket);
2390  if (syzstr->syRing != NULL)
2391  {
2392    rChangeCurrRing(origR);
2393    // Thomas: now make sure that all data which you need is pFetchCopied
2394    // maybe incoporate it into syReorder ??
2395  }
2396  tres = syReorder(tres,syzstr->length,syzstr,FALSE,syzstr->res);
2397  syKillEmptyEntres(tres,syzstr->length);
2398  idSkipZeroes(tres[0]);
2399  return tres;
2400}
2401
2402/*3
2403* minimizes any kind of resolution
2404*/
2405syStrategy syMinimize(syStrategy syzstr)
2406{
2407  if (syzstr->minres==NULL)
2408  {
2409    if (syzstr->resPairs!=NULL)
2410    {
2411      if (syzstr->hilb_coeffs==NULL)
2412      {
2413        syzstr->minres = syReadOutMinimalRes(syzstr);
2414      }
2415      else
2416      {
2417        syzstr->minres = syReorder(syzstr->orderedRes,syzstr->length,syzstr);
2418      }
2419    }
2420    else if (syzstr->fullres!=NULL)
2421    {
2422      syMinimizeResolvente(syzstr->fullres,syzstr->length,1);
2423      syzstr->minres = syzstr->fullres;
2424      syzstr->fullres = NULL;
2425    }
2426  }
2427  (syzstr->references)++;
2428  return syzstr;
2429}
2430
2431/*2
2432* implementation of LaScala's algorithm
2433* assumes that the given module is homogeneous
2434* works with slanted degree, uses syChosePairs
2435*/
2436syStrategy syLaScala3(ideal arg,int * length)
2437{
2438  BOOLEAN noPair=FALSE;
2439  int i,j,actdeg=32000,index=0,reg=-1;
2440  int startdeg,howmuch;
2441  poly p;
2442  ideal temp;
2443  SSet nextPairs;
2444  syStrategy syzstr=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2445  ring origR = currRing;
2446
2447  if ((idIs0(arg)) ||
2448      ((idRankFreeModule(arg)>0) && (!idHomModule(arg,NULL,&(syzstr->cw)))))
2449  {
2450    syzstr->minres = (resolvente)omAlloc0Bin(char_ptr_bin);
2451    syzstr->length = 1;
2452    syzstr->minres[0] = idInit(1,arg->rank);
2453    return syzstr;
2454  }
2455
2456  //crit = 0;
2457  //euler = -1;
2458  redpol = pInit();
2459  syzstr->length = *length = pVariables+2;
2460
2461  // Creare dp,S ring and change to it
2462  syzstr->syRing = rCurrRingAssure_dp_S();
2463  assume(syzstr->syRing != origR);
2464
2465  // set initial ShiftedComps
2466  currcomponents = (int*)omAlloc0((arg->rank+1)*sizeof(int));
2467  currShiftedComponents = (long*)omAlloc0((arg->rank+1)*sizeof(long));
2468  for (i=0;i<=arg->rank;i++)
2469  {
2470    currShiftedComponents[i] = (i)*SYZ_SHIFT_BASE;
2471    currcomponents[i] = i;
2472  }
2473  rChangeSComps(currcomponents, currShiftedComponents, arg->rank);
2474/*--- initializes the data structures---------------*/
2475  syzstr->Tl = new intvec(*length);
2476  temp = idInit(IDELEMS(arg),arg->rank);
2477  for (i=0;i<IDELEMS(arg);i++)
2478  {
2479    temp->m[i] = prCopyR( arg->m[i], origR);
2480    if (temp->m[i]!=NULL)
2481    {
2482      j = pTotaldegree(temp->m[i]);
2483      if (j<actdeg) actdeg = j;
2484    }
2485  }
2486  idTest(temp);
2487  idSkipZeroes(temp);
2488  idTest(temp);
2489  syzstr->resPairs = syInitRes(temp,length,syzstr->Tl,syzstr->cw);
2490  omFreeSize((ADDRESS)currcomponents,(arg->rank+1)*sizeof(int));
2491  omFreeSize((ADDRESS)currShiftedComponents,(arg->rank+1)*sizeof(long));
2492  syzstr->res = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
2493  syzstr->orderedRes = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
2494  syzstr->elemLength = (int**)omAlloc0((*length+1)*sizeof(int*));
2495  syzstr->truecomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
2496  syzstr->ShiftedComponents = (long**)omAlloc0((*length+1)*sizeof(long*));
2497  syzstr->backcomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
2498  syzstr->Howmuch = (int**)omAlloc0((*length+1)*sizeof(int*));
2499  syzstr->Firstelem = (int**)omAlloc0((*length+1)*sizeof(int*));
2500  syzstr->sev = (unsigned long **) omAlloc0((*length+1)*sizeof(unsigned long *));
2501  syzstr->bucket = kBucketCreate();
2502  int len0=idRankFreeModule(temp)+1;
2503
2504  startdeg = actdeg;
2505  nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2506  //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2507/*--- computes the resolution ----------------------*/
2508  while (nextPairs!=NULL)
2509  {
2510    if (TEST_OPT_PROT) Print("%d",actdeg);
2511    if (TEST_OPT_PROT) Print("(m%d)",index);
2512    if (index==0)
2513      i = syInitSyzMod(syzstr,index,len0);
2514    else
2515      i = syInitSyzMod(syzstr,index);
2516    currcomponents = syzstr->truecomponents[si_max(index-1,0)];
2517    currShiftedComponents = syzstr->ShiftedComponents[si_max(index-1,0)];
2518    rChangeSComps(currcomponents, currShiftedComponents,
2519                  IDELEMS(syzstr->res[si_max(index-1,0)]));
2520    j = syInitSyzMod(syzstr,index+1);
2521    if (index>0)
2522    {
2523      syRedNextPairs(nextPairs,syzstr,howmuch,index);
2524      syCompactifyPairSet(syzstr->resPairs[index],(*syzstr->Tl)[index],0);
2525    }
2526    else
2527      syRedGenerOfCurrDeg(syzstr,actdeg,index+1);
2528/*--- creates new pairs -----------------------------*/
2529    syCreateNewPairs(syzstr,index,i);
2530    if (index<(*length)-1)
2531    {
2532      syCreateNewPairs(syzstr,index+1,j);
2533    }
2534    index++;
2535    nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2536    //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2537  }
2538  if (temp!=NULL) idDelete(&temp);
2539  kBucketDestroy(&(syzstr->bucket));
2540  if (origR != syzstr->syRing)
2541    rChangeCurrRing(origR);
2542  pLmDelete(&redpol);
2543
2544  if (TEST_OPT_PROT) PrintLn();
2545
2546  if (syzstr->minres==NULL)
2547      if (syzstr->resPairs!=NULL)
2548        if (syzstr->hilb_coeffs==NULL)
2549          syzstr->minres = syReadOutMinimalRes(syzstr);
2550 
2551  return syzstr;
2552}
2553
2554
2555
2556/*2
2557* more general implementation of LaScala's algorithm
2558* assumes that the given module is (quasi-)homogeneous
2559* works with slanted degree, uses syChosePairs
2560*/
2561syStrategy syLaScala(ideal arg, int& maxlength, intvec* weights)
2562{
2563  BOOLEAN noPair=FALSE;
2564  int i,j,actdeg=32000,index=0,reg=-1;
2565  int startdeg,howmuch;
2566  poly p;
2567  ideal temp;
2568  SSet nextPairs;
2569  syStrategy syzstr=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2570  ring origR = currRing;
2571
2572  if(weights!= NULL) 
2573    syzstr->cw = new intvec(weights);
2574  else
2575    syzstr->cw = NULL;
2576
2577  if ((idIs0(arg)) ||
2578      ((idRankFreeModule(arg)>0) && (!idTestHomModule(arg, NULL, syzstr->cw))))
2579  {
2580    syzstr->minres = (resolvente)omAlloc0Bin(char_ptr_bin);
2581    syzstr->length = 1;
2582    syzstr->minres[0] = idInit(1,arg->rank);
2583    return syzstr;
2584  }
2585
2586
2587  //crit = 0;
2588  //euler = -1;
2589  redpol = pInit();
2590 
2591  if( maxlength > 0 )
2592    syzstr->length = maxlength; //  = pVariables+2;
2593  else
2594    syzstr->length = maxlength = pVariables+2; 
2595
2596  // Creare dp,S ring and change to it
2597  syzstr->syRing = rCurrRingAssure_dp_S();
2598  assume(syzstr->syRing != origR);
2599  assume(syzstr->syRing->typ[1].ord_typ == ro_syzcomp); 
2600
2601
2602  // set initial ShiftedComps
2603  currcomponents = (int*)omAlloc0((arg->rank+1)*sizeof(int));
2604  currShiftedComponents = (long*)omAlloc0((arg->rank+1)*sizeof(long));
2605  for (i=0;i<=arg->rank;i++)
2606  {
2607    currShiftedComponents[i] = (i)*SYZ_SHIFT_BASE;
2608    currcomponents[i] = i;
2609  }
2610  rChangeSComps(currcomponents, currShiftedComponents, arg->rank);
2611/*--- initializes the data structures---------------*/
2612  syzstr->Tl = new intvec(maxlength);
2613  temp = idInit(IDELEMS(arg),arg->rank);
2614  for (i=0;i<IDELEMS(arg);i++)
2615  {
2616    temp->m[i] = prCopyR( arg->m[i], origR);
2617    if (temp->m[i]!=NULL)
2618    {
2619      j = pTotaldegree(temp->m[i]);
2620      if (j<actdeg) actdeg = j;
2621    }
2622  }
2623  idTest(temp);
2624  idSkipZeroes(temp);
2625  idTest(temp);
2626  syzstr->resPairs = syInitRes(temp,&maxlength,syzstr->Tl,syzstr->cw);
2627  omFreeSize((ADDRESS)currcomponents,(arg->rank+1)*sizeof(int));
2628  omFreeSize((ADDRESS)currShiftedComponents,(arg->rank+1)*sizeof(long));
2629 
2630  syzstr->res = (resolvente)omAlloc0((maxlength+1)*sizeof(ideal));
2631  syzstr->orderedRes = (resolvente)omAlloc0((maxlength+1)*sizeof(ideal));
2632  syzstr->elemLength = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2633
2634  syzstr->truecomponents = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2635  syzstr->ShiftedComponents = (long**)omAlloc0((maxlength+1)*sizeof(long*));
2636
2637  syzstr->backcomponents = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2638  syzstr->Howmuch = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2639  syzstr->Firstelem = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2640  syzstr->sev = (unsigned long **) omAlloc0((maxlength+1)*sizeof(unsigned long *));
2641
2642  assume( syzstr->length == maxlength );
2643 
2644  syzstr->bucket = kBucketCreate();
2645  int len0=idRankFreeModule(temp)+1;
2646
2647  startdeg = actdeg;
2648  nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2649  //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2650/*--- computes the resolution ----------------------*/
2651  while (nextPairs!=NULL)
2652  {
2653    if (TEST_OPT_PROT) Print("%d",actdeg);
2654    if (TEST_OPT_PROT) Print("(m%d)",index);
2655    if (index==0)
2656      i = syInitSyzMod(syzstr,index,len0);
2657    else
2658      i = syInitSyzMod(syzstr,index);
2659    currcomponents = syzstr->truecomponents[si_max(index-1,0)];
2660    currShiftedComponents = syzstr->ShiftedComponents[si_max(index-1,0)];
2661    rChangeSComps(currcomponents, currShiftedComponents,
2662                  IDELEMS(syzstr->res[si_max(index-1,0)]));
2663    j = syInitSyzMod(syzstr,index+1);
2664    if (index>0)
2665    {
2666      syRedNextPairs(nextPairs,syzstr,howmuch,index);
2667      syCompactifyPairSet(syzstr->resPairs[index],(*syzstr->Tl)[index],0);
2668    }
2669    else
2670      syRedGenerOfCurrDeg(syzstr,actdeg,index+1);
2671/*--- creates new pairs -----------------------------*/
2672    syCreateNewPairs(syzstr,index,i);
2673    if (index<(maxlength-1))
2674    {
2675      syCreateNewPairs(syzstr,index+1,j);
2676    }
2677    index++;
2678    nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2679    //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2680  }
2681  if (temp!=NULL) idDelete(&temp);
2682  kBucketDestroy(&(syzstr->bucket));
2683  if (origR != syzstr->syRing)
2684    rChangeCurrRing(origR);
2685  pLmDelete(&redpol);
2686  if (TEST_OPT_PROT) PrintLn();
2687  return syzstr;
2688}
2689
Note: See TracBrowser for help on using the repository browser.