source: git/Singular/syz1.cc @ 35aab3

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