source: git/Singular/syz1.cc @ 0224c8

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