source: git/Singular/syz1.cc @ 28e0ac8

fieker-DuValspielwiese
Last change on this file since 28e0ac8 was 16acb0, checked in by Hans Schönemann <hannes@…>, 22 years ago
* hannes: more ring-independence in changing rings, deleting data etc. (see lists::Clean(r), sleftv::CleanUp(r), syKillComputation(r), nap_Delete) git-svn-id: file:///usr/local/Singular/svn/trunk@5708 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 68.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: syz1.cc,v 1.74 2002-01-10 12:33:24 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      coefgcd =
901        nGcd(pGetCoeff(tso.p1),pGetCoeff(tso.p2),currRing);
902      tso.syz = pHead(tso.lcm);
903      p = tso.syz;
904      pSetCoeff(p,nDiv(pGetCoeff(tso.p1),coefgcd));
905      pGetCoeff(p) = nNeg(pGetCoeff(p));
906      pSetComp(p,tso.ind2+1);
907      p_Setm_Syz(p, currRing, Components, ShiftedComponents); // actueller index
908      pNext(p) = pHead(tso.lcm);
909      pIter(p);
910      pSetComp(p,tso.ind1+1);
911      p_Setm_Syz(p, currRing, Components, ShiftedComponents); // actueller index
912      pSetCoeff(p,nDiv(pGetCoeff(tso.p2),coefgcd));
913      nDelete(&coefgcd);
914      if (tso.p != NULL)
915      {
916        kBucketInit(syzstr->bucket,tso.p,-1);
917        q = kBucketGetLm(syzstr->bucket);
918        j = Fin[pGetComp(q)]-1;
919        int pos = j+Hin[pGetComp(q)];
920        loop
921        {
922          if (j<0) break;
923          if (pLmDivisibleByNoComp(redset[j],q))
924          {
925            pNext(p) = pHead(q);
926            pIter(p);
927            pSetComp(p,bin[j]+1);
928            p_Setm_Syz(p, currRing, Components, ShiftedComponents); // actueller index
929//if (pLength(redset[j])!=syzstr->elemLength[index][bin[j]])
930//Print("Halt");
931//if (pLength(redset[j])!=syzstr->elemLength[index][bin[j]])
932//Print("Halt");
933            pGetCoeff(p) = nNeg(pGetCoeff(p));
934            number up = kBucketPolyRed(syzstr->bucket,redset[j],elL[bin[j]],
935                                       NULL);
936            // Thomas: Check whether you need number here
937            nDelete(&up);
938            q = kBucketGetLm(syzstr->bucket);
939            if (q==NULL) break;
940            j = Fin[pGetComp(q)]-1;
941            pos = j+Hin[pGetComp(q)];
942          }
943          else
944          {
945            j++;
946            if (j==pos) break;
947          }
948        }
949        int lb;
950        kBucketClear(syzstr->bucket,&tso.p,&lb);
951      }
952      if (tso.p != NULL)
953      {
954        if (TEST_OPT_PROT) Print("g");
955        if (k==IDELEMS((syzstr->res)[index]))
956        {
957          syEnlargeFields(syzstr,index);
958          bin=syzstr->backcomponents[index];
959          elL=syzstr->elemLength[index];
960          redset=syzstr->orderedRes[index]->m;
961          Components = syzstr->truecomponents[index];
962          ShiftedComponents = syzstr->ShiftedComponents[index];
963        }
964        pNext(p) = pHead(tso.p);
965        pIter(p);
966
967        assume(p!= NULL);
968        k++;
969        syzstr->res[index]->m[k-1] = tso.p;
970        syzstr->elemLength[index][k-1] = pLength(tso.p);
971        pNorm(syzstr->res[index]->m[k-1]);
972        need_reset = syOrder(syzstr->res[index]->m[k-1],syzstr,index,k);
973        pSetComp(p,k); // actueller index
974        p_Setm_Syz(p, currRing, Components, ShiftedComponents);
975        pGetCoeff(p) = nNeg(pGetCoeff(p));
976
977        tso.isNotMinimal = p;
978        tso.p = NULL;
979      }
980      else
981      {
982        if (TEST_OPT_PROT) Print(".");
983        //if (index % 2==0)
984          //euler++;
985        //else
986          //euler--;
987      }
988      if (ks==IDELEMS(syzstr->res[index+1]))
989      {
990        syEnlargeFields(syzstr,index+1);
991      }
992      syzstr->res[index+1]->m[ks] = tso.syz;
993      syzstr->elemLength[index+1][ks] = pLength(tso.syz);
994      pNorm(syzstr->res[index+1]->m[ks]);
995      tso.syz =NULL;
996      tso.syzind = ks;
997      if (need_reset)
998        syResetShiftedComponents(syzstr, index+1);
999      if (syOrder(syzstr->res[index+1]->m[ks],syzstr,index+1,ks+1))
1000         syResetShiftedComponents(syzstr, index+2);
1001      ks++;
1002      p = NULL;
1003      nextPairs[(*spl1)[i]-1] = tso;
1004    }
1005    i++;
1006  }
1007  delete spl1;
1008}
1009
1010/*3
1011* reduces the generators of the module index in degree deg
1012* (which are actual syzygies of the module index-1)
1013* wrt. the ideal generated by elements of lower degrees
1014*/
1015static void syRedGenerOfCurrDeg(syStrategy syzstr, int deg, int index)
1016{
1017  ideal res=syzstr->res[index];
1018  int i=0,j,k=IDELEMS(res),kk;
1019  SSet sPairs=syzstr->resPairs[index-1];
1020
1021  while ((k>0) && (res->m[k-1]==NULL)) k--;
1022  while ((i<(*syzstr->Tl)[index-1]) && (((sPairs)[i].syz==NULL) ||
1023          ((sPairs)[i].order<deg)))
1024    i++;
1025  if ((i>=(*syzstr->Tl)[index-1]) || ((sPairs)[i].order>deg)) return;
1026  while ((i<(*syzstr->Tl)[index-1]) && (((sPairs)[i].syz==NULL) ||
1027         ((sPairs)[i].order==deg)))
1028  {
1029    if ((sPairs)[i].syz!=NULL)
1030    {
1031      j = k-1;
1032      while ((j>=0) && (res->m[j]!=NULL) &&
1033             ((sPairs)[i].syz!=NULL))
1034      {
1035        if (pLmDivisibleBy(res->m[j],(sPairs)[i].syz))
1036        {
1037          (sPairs)[i].syz =
1038            ksOldSpolyRed(res->m[j],(sPairs)[i].syz);
1039            //sySPolyRed((sPairs)[i].syz,res->m[j]);
1040          j = k-1;
1041        }
1042        else
1043        {
1044          j--;
1045        }
1046      }
1047      if ((sPairs)[i].syz != NULL)
1048      {
1049        if (k==IDELEMS(res))
1050        {
1051          syEnlargeFields(syzstr,index);
1052          res=syzstr->res[index];
1053        }
1054        if (TEST_OPT_DEBUG)
1055        {
1056          if ((sPairs)[i].isNotMinimal==NULL)
1057          {
1058            PrintLn();
1059            Print("minimal generator: ");pWrite((syzstr->resPairs[index-1])[i].syz);
1060            Print("comes from: ");pWrite((syzstr->resPairs[index-1])[i].p1);
1061            Print("and: ");pWrite((syzstr->resPairs[index-1])[i].p2);
1062          }
1063        }
1064        //res->m[k] = (sPairs)[i].syz;
1065        res->m[k] = syRedtail((sPairs)[i].syz,syzstr,index);
1066        (sPairs)[i].syzind = k;
1067        syzstr->elemLength[index][k] = pLength((sPairs)[i].syz);
1068        pNorm(res->m[k]);
1069  //      (sPairs)[i].syz = NULL;
1070        k++;
1071        if (syOrder(res->m[k-1],syzstr,index,k))
1072          syResetShiftedComponents(syzstr, index);
1073        //euler++;
1074      }
1075      else
1076        (sPairs)[i].syzind = -1;
1077    }
1078    i++;
1079  }
1080}
1081
1082/*3
1083* puts a pair into the right place in resPairs
1084*/
1085void syEnterPair(SSet sPairs, SObject * so, int * sPlength,int index)
1086{
1087  int ll,k,no=(*so).order,sP=*sPlength,i;
1088  poly p=(*so).lcm;
1089
1090  if ((sP==0) || (sPairs[sP-1].order<=no))
1091    ll = sP;
1092  else if (sP==1)
1093    ll = 0;
1094  else
1095  {
1096    int an=0,en=sP-1;
1097    loop
1098    {
1099      if (an>=en-1)
1100      {
1101        if ((sPairs[an].order<=no) && (sPairs[an+1].order>no))
1102        {
1103          ll = an+1;
1104          break;
1105        }
1106        else if ((sPairs[en].order<=no) && (sPairs[en+1].order>no))
1107        {
1108          ll = en+1;
1109          break;
1110        }
1111        else if (sPairs[an].order>no)
1112        {
1113          ll = an;
1114          break;
1115        }
1116        else
1117        {
1118          Print("Hier ist was faul!\n");
1119          break;
1120        }
1121      }
1122      i=(an+en) / 2;
1123      if (sPairs[i].order <= no)
1124        an=i;
1125      else
1126        en=i;
1127    }
1128  }
1129  for (k=(*sPlength);k>ll;k--)
1130  {
1131    syCopyPair(&sPairs[k-1],&sPairs[k]);
1132  }
1133  syCopyPair(so,&sPairs[ll]);
1134  (*sPlength)++;
1135}
1136void syEnterPair(syStrategy syzstr, SObject * so, int * sPlength,int index)
1137{
1138  int ll;
1139
1140  if (*sPlength>=(*syzstr->Tl)[index])
1141  {
1142    SSet temp = (SSet)omAlloc0(((*syzstr->Tl)[index]+16)*sizeof(SObject));
1143    for (ll=0;ll<(*syzstr->Tl)[index];ll++)
1144    {
1145      temp[ll].p = (syzstr->resPairs[index])[ll].p;
1146      temp[ll].p1 = (syzstr->resPairs[index])[ll].p1;
1147      temp[ll].p2 = (syzstr->resPairs[index])[ll].p2;
1148      temp[ll].syz = (syzstr->resPairs[index])[ll].syz;
1149      temp[ll].lcm = (syzstr->resPairs[index])[ll].lcm;
1150      temp[ll].ind1 = (syzstr->resPairs[index])[ll].ind1;
1151      temp[ll].ind2 = (syzstr->resPairs[index])[ll].ind2;
1152      temp[ll].syzind = (syzstr->resPairs[index])[ll].syzind;
1153      temp[ll].order = (syzstr->resPairs[index])[ll].order;
1154      temp[ll].isNotMinimal = (syzstr->resPairs[index])[ll].isNotMinimal;
1155      temp[ll].length = (syzstr->resPairs[index])[ll].length;
1156      temp[ll].reference = (syzstr->resPairs[index])[ll].reference;
1157    }
1158    if (syzstr->resPairs[index] != NULL) // OB: ?????
1159      omFreeSize((ADDRESS)syzstr->resPairs[index],(*syzstr->Tl)[index]*sizeof(SObject));
1160    (*syzstr->Tl)[index] += 16;
1161    syzstr->resPairs[index] = temp;
1162  }
1163  syEnterPair(syzstr->resPairs[index],so,sPlength,index);
1164}
1165
1166/*3
1167* computes pairs from the new elements (beginning with the element newEl)
1168* in the module index
1169*/
1170static void syCreateNewPairs(syStrategy syzstr, int index, int newEl)
1171{
1172  SSet temp;
1173  SObject tso;
1174  int i,ii,j,k=IDELEMS(syzstr->res[index]),l=(*syzstr->Tl)[index],ll;
1175  int qc,first,pos,jj,j1;
1176  int * bci=syzstr->backcomponents[index];
1177  poly p,q;
1178  polyset rs=syzstr->res[index]->m,nPm;
1179
1180
1181  while ((k>0) && (rs[k-1]==NULL)) k--;
1182  if (newEl>=k) return;
1183
1184  long * ShiftedComponents = syzstr->ShiftedComponents[index];
1185  int* Components = syzstr->truecomponents[index];
1186
1187  ideal nP=idInit(k,syzstr->res[index]->rank);
1188  nPm=nP->m;
1189  while ((l>0) && ((syzstr->resPairs[index])[l-1].p1==NULL)) l--;
1190  for (j=newEl;j<k;j++)
1191  {
1192    q = rs[j];
1193    qc = pGetComp(q);
1194    first = syzstr->Firstelem[index-1][pGetComp(q)]-1;
1195    pos = first+syzstr->Howmuch[index-1][pGetComp(q)];
1196    for (i=first;i<pos;i++)
1197    {
1198      jj = bci[i];
1199      if (jj>=j) break;
1200      p = pOne();
1201      pLcm(rs[jj],q,p);
1202      pSetComp(p,j+1);
1203      p_Setm_Syz(p, currRing, Components, ShiftedComponents);
1204      ii = first;
1205      loop
1206      {
1207        j1 = bci[ii];
1208        if (nPm[j1]!=NULL)
1209        {
1210          if (pLmDivisibleByNoComp(nPm[j1],p))
1211          {
1212            pDelete(&p);
1213            break;
1214          }
1215          else if (pLmDivisibleByNoComp(p,nPm[j1]))
1216          {
1217            pDelete(&(nPm[j1]));
1218            //break;
1219          }
1220        }
1221        ii++;
1222        if (ii>=pos) break;
1223      }
1224      if (p!=NULL)
1225      {
1226        nPm[jj] = p;
1227      }
1228    }
1229    for (i=first;i<pos;i++)
1230    {
1231      ii = bci[i];
1232      if (nPm[ii]!=NULL)
1233      {
1234        if (l>=(*syzstr->Tl)[index])
1235        {
1236          temp = (SSet)omAlloc0(((*syzstr->Tl)[index]+16)*sizeof(SObject));
1237          for (ll=0;ll<(*syzstr->Tl)[index];ll++)
1238          {
1239            temp[ll].p = (syzstr->resPairs[index])[ll].p;
1240            temp[ll].p1 = (syzstr->resPairs[index])[ll].p1;
1241            temp[ll].p2 = (syzstr->resPairs[index])[ll].p2;
1242            temp[ll].syz = (syzstr->resPairs[index])[ll].syz;
1243            temp[ll].lcm = (syzstr->resPairs[index])[ll].lcm;
1244            temp[ll].ind1 = (syzstr->resPairs[index])[ll].ind1;
1245            temp[ll].ind2 = (syzstr->resPairs[index])[ll].ind2;
1246            temp[ll].syzind = (syzstr->resPairs[index])[ll].syzind;
1247            temp[ll].order = (syzstr->resPairs[index])[ll].order;
1248            temp[ll].isNotMinimal = (syzstr->resPairs[index])[ll].isNotMinimal;
1249          }
1250          if (syzstr->resPairs[index] != NULL) // OB: ????
1251            omFreeSize((ADDRESS)syzstr->resPairs[index],(*syzstr->Tl)[index]*sizeof(SObject));
1252          (*syzstr->Tl)[index] += 16;
1253          syzstr->resPairs[index] = temp;
1254        }
1255        tso.lcm = p = nPm[ii];
1256        nPm[ii] = NULL;
1257        //#ifdef HAVE_SHIFTED_EXPONENTS
1258        //tso.order = pTotaldegree(p);
1259        //p->exp[currRing->pOrdIndex]=tso.order+0x40000000;
1260        //#else
1261        tso.order = pTotaldegree(p);
1262        pSetOrder(p, tso.order);
1263        //#endif
1264        if ((syzstr->cw!=NULL) && (index>0) && (pGetComp(q)>0))
1265        {
1266          int ii=index-1,jj=pGetComp(q);
1267          while (ii>0)
1268          {
1269            jj = pGetComp(syzstr->res[ii]->m[jj-1]);
1270            ii--;
1271          }
1272          tso.order += (*syzstr->cw)[jj-1];
1273        }
1274        tso.p1 = rs[ii];
1275        tso.p2 = q;
1276        tso.ind1 = ii;
1277        tso.ind2 = j;
1278        tso.syzind = -1;
1279        tso.isNotMinimal = NULL;
1280        tso.p = NULL;
1281        tso.syz = NULL;
1282        syEnterPair(syzstr->resPairs[index],&tso,&l,index);
1283      }
1284    }
1285  }
1286  idDelete(&nP);
1287}
1288
1289static SSet syChosePairsPutIn(syStrategy syzstr, int *index,
1290               int *howmuch, int * actdeg, int an, int en)
1291{
1292  int newdeg=*actdeg,newindex=-1,i,t,sldeg;
1293  poly p;
1294  SSet result;
1295  SRes resPairs=syzstr->resPairs;
1296
1297  if (an>syzstr->length) return NULL;
1298  if (en>syzstr->length) en=syzstr->length;
1299  while (*index<en)
1300  {
1301    if (resPairs[*index]!=NULL)
1302    {
1303      sldeg = (*actdeg)+*index;
1304      i = 0;
1305      if (*index!=0)
1306      {
1307        while ((i<(*syzstr->Tl)[*index]))
1308        {
1309          if ((resPairs[*index])[i].lcm!=NULL)
1310          {
1311            if ((resPairs[*index])[i].order == sldeg)
1312            {
1313              result = &(resPairs[*index])[i];
1314              *howmuch =1;
1315              i++;
1316              while ((i<(*syzstr->Tl)[*index]) && ((resPairs[*index])[i].lcm!=NULL)
1317                      && ((resPairs[*index])[i].order == sldeg))
1318              {
1319                i++;
1320                (*howmuch)++;
1321              }
1322              return result;
1323            }
1324          }
1325          i++;
1326        }
1327      }
1328      else
1329      {
1330        while ((i<(*syzstr->Tl)[*index]))
1331        {
1332          if ((resPairs[*index])[i].syz!=NULL)
1333          {
1334            if ((resPairs[*index])[i].order == sldeg)
1335            {
1336              result = &(resPairs[*index])[i];
1337              (*howmuch) =1;
1338              i++;
1339              while ((i<(*syzstr->Tl)[*index]) && ((resPairs[*index])[i].syz!=NULL)
1340                      && ((resPairs[*index])[i].order == *actdeg))
1341              {
1342                i++;
1343                (*howmuch)++;
1344              }
1345              return result;
1346            }
1347          }
1348          i++;
1349        }
1350      }
1351    }
1352    (*index)++;
1353  }
1354  *index = an;
1355  //if (TEST_OPT_PROT) Print("(Euler:%d)",euler);
1356  while (*index<en)
1357  {
1358    if (resPairs[*index]!=NULL)
1359    {
1360      i = 0;
1361      while ((i<(*syzstr->Tl)[*index]))
1362      {
1363        t = *actdeg+*index;
1364        if (((resPairs[*index])[i].lcm!=NULL) ||
1365              ((resPairs[*index])[i].syz!=NULL))
1366        {
1367          if ((resPairs[*index])[i].order > t)
1368            t = (resPairs[*index])[i].order;
1369        }
1370        if ((t>*actdeg+*index) && ((newdeg==*actdeg) || (t<newdeg+*index)))
1371        {
1372          newdeg = t-*index;
1373          newindex = *index;
1374          break;
1375        }
1376        i++;
1377      }
1378    }
1379    (*index)++;
1380  }
1381  if (newdeg>*actdeg)
1382  {
1383    *actdeg = newdeg;
1384    *index = newindex;
1385    return syChosePairsPutIn(syzstr,index,howmuch,actdeg,an,en);
1386  }
1387  else return NULL;
1388}
1389
1390/*3
1391* FOR THE HOMOGENEOUS CASE ONLY!
1392* looks through the pair set and the given module for
1393* remaining pairs or generators to consider
1394* returns a pointer to the first pair and the number of them in the given module
1395* works with slanted degree (i.e. deg=realdeg-index)
1396*/
1397SSet syChosePairs(syStrategy syzstr, int *index, int *howmuch, int * actdeg)
1398{
1399  return syChosePairsPutIn(syzstr,index,howmuch,actdeg,0,syzstr->length);
1400}
1401
1402/*3
1403* FOR THE INHOMOGENEOUS CASE ONLY!
1404* looks through the pair set and the given module for
1405* remaining pairs or generators to consider
1406* returns a pointer to the first pair and the number of them in the given module
1407* works with slanted degree (i.e. deg=realdeg-index)
1408* looks first through the 0 and 1 module then through the other
1409*/
1410static SSet syChosePairsIH(syStrategy syzstr, int *index,
1411               int *howmuch, int * actdeg, int mindeg)
1412{
1413  SSet result=NULL;
1414
1415  result = syChosePairsPutIn(syzstr,index,howmuch,actdeg,0,2);
1416  if (result == NULL)
1417  {
1418    *actdeg = mindeg;
1419    result = syChosePairsPutIn(syzstr,index,howmuch,actdeg,2,syzstr->length);
1420  }
1421  return result;
1422}
1423
1424/*3
1425* looks through the pair set and the given module for
1426* remaining pairs or generators to consider
1427* returns a pointer to the first pair and the number of them in the given module
1428* works deg by deg
1429*/
1430/*
1431*static SSet syChosePairs1(SRes resPairs,intvec * Tl, int *index, int *howmuch,
1432*                   int length,int * actdeg)
1433*{
1434*  int newdeg=*actdeg,newindex=-1,i,t;
1435*  SSet result;
1436*
1437*  while (*index>=0)
1438*  {
1439*    if (resPairs[*index]!=NULL)
1440*    {
1441*      i = 0;
1442*      if (*index!=0)
1443*      {
1444*        while ((i<(*Tl)[*index]))
1445*        {
1446*          if ((resPairs[*index])[i].lcm!=NULL)
1447*          {
1448*            if (pGetOrder((resPairs[*index])[i].lcm) == *actdeg)
1449*            {
1450*              result = &(resPairs[*index])[i];
1451*              *howmuch =1;
1452*              i++;
1453*              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].lcm!=NULL)
1454*                      && (pGetOrder((resPairs[*index])[i].lcm) == *actdeg))
1455*              {
1456*                i++;
1457*                (*howmuch)++;
1458*              }
1459*              return result;
1460*            }
1461*          }
1462*          i++;
1463*        }
1464*      }
1465*      else
1466*      {
1467*        while ((i<(*Tl)[*index]))
1468*        {
1469*          if ((resPairs[*index])[i].syz!=NULL)
1470*          {
1471*            if ((resPairs[*index])[i].order == *actdeg)
1472*            {
1473*              result = &(resPairs[*index])[i];
1474*              (*howmuch) =1;
1475*              i++;
1476*              while ((i<(*Tl)[*index]) && ((resPairs[*index])[i].syz!=NULL)
1477*                      && ((resPairs[*index])[i].order == *actdeg))
1478*              {
1479*                i++;
1480*                (*howmuch)++;
1481*              }
1482*              return result;
1483*            }
1484*          }
1485*          i++;
1486*        }
1487*      }
1488*    }
1489*    (*index)--;
1490*  }
1491*  *index = length-1;
1492*  while (*index>=0)
1493*  {
1494*    if (resPairs[*index]!=NULL)
1495*    {
1496*      i = 0;
1497*      while ((i<(*Tl)[*index]))
1498*      {
1499*        t = *actdeg;
1500*        if ((resPairs[*index])[i].lcm!=NULL)
1501*        {
1502*          if (pGetOrder((resPairs[*index])[i].lcm) > *actdeg)
1503*            t = pGetOrder((resPairs[*index])[i].lcm);
1504*        }
1505*        else if ((resPairs[*index])[i].syz!=NULL)
1506*        {
1507*          if ((resPairs[*index])[i].order > *actdeg)
1508*            t = (resPairs[*index])[i].order;
1509*        }
1510*        if ((t>*actdeg) && ((newdeg==*actdeg) || (t<newdeg)))
1511*        {
1512*          newdeg = t;
1513*          newindex = *index;
1514*          break;
1515*        }
1516*        i++;
1517*      }
1518*    }
1519*    (*index)--;
1520*  }
1521*  if (newdeg>*actdeg)
1522*  {
1523*    *actdeg = newdeg;
1524*    *index = newindex;
1525*    return syChosePairs1(resPairs,Tl,index,howmuch,length,actdeg);
1526*  }
1527*  else return NULL;
1528*}
1529*/
1530/*3
1531* statistics of the resolution
1532*/
1533static void syStatistics(resolvente res,int length)
1534{
1535  int i,j=1,k;
1536  Order_t deg = 0;
1537
1538  PrintLn();
1539  while ((j<length) && (res[j]!=NULL))
1540  {
1541    Print("In module %d: \n",j);
1542    k = 0;
1543    while ((k<IDELEMS(res[j])) && (res[j]->m[k]!=NULL))
1544    {
1545      i = 1;
1546      deg = pGetOrder(res[j]->m[k]);
1547      k++;
1548      while ((k<IDELEMS(res[j])) && (res[j]->m[k]!=NULL) &&
1549              (pGetOrder(res[j]->m[k])==deg))
1550      {
1551        i++;
1552        k++;
1553      }
1554      Print("%d elements of degree %d\n",i,deg);
1555    }
1556    j++;
1557  }
1558}
1559
1560/*3
1561* initialize a module
1562*/
1563int syInitSyzMod(syStrategy syzstr, int index, int init)
1564{
1565  int result;
1566
1567  if (syzstr->res[index]==NULL)
1568  {
1569    syzstr->res[index] = idInit(init-1,1);
1570    syzstr->truecomponents[index] = (int*)omAlloc0(init*sizeof(int));
1571    syzstr->ShiftedComponents[index] = (long*)omAlloc0(init*sizeof(long));
1572    if (index==0)
1573    {
1574      for (int i=0;i<init;i++)
1575      {
1576        syzstr->truecomponents[0][i] = i;
1577        syzstr->ShiftedComponents[0][i] = (i)*SYZ_SHIFT_BASE;
1578      }
1579    }
1580    syzstr->backcomponents[index] = (int*)omAlloc0(init*sizeof(int));
1581    syzstr->Howmuch[index] = (int*)omAlloc0(init*sizeof(int));
1582    syzstr->Firstelem[index] = (int*)omAlloc0(init*sizeof(int));
1583    syzstr->elemLength[index] = (int*)omAlloc0(init*sizeof(int));
1584    syzstr->orderedRes[index] = idInit(init-1,1);
1585    syzstr->sev[index] = (unsigned long*) omAlloc0(init*sizeof(unsigned long));
1586    result = 0;
1587  }
1588  else
1589  {
1590    result = IDELEMS(syzstr->res[index]);
1591    while ((result>0) && (syzstr->res[index]->m[result-1]==NULL)) result--;
1592  }
1593  return result;
1594}
1595
1596/*3
1597* deletes a resolution
1598*/
1599void syKillComputation(syStrategy syzstr, ring r)
1600{
1601  if (syzstr->references>0)
1602  {
1603    (syzstr->references)--;
1604  }
1605  else
1606  {
1607    int i,j;
1608    if (syzstr->minres!=NULL)
1609    {
1610      for (i=0;i<syzstr->length;i++)
1611      {
1612        if (syzstr->minres[i]!=NULL)
1613        {
1614          for (j=0;j<IDELEMS(syzstr->minres[i]);j++)
1615          {
1616            if (syzstr->minres[i]->m[j]!=NULL)
1617              p_Delete(&(syzstr->minres[i]->m[j]),r);
1618          }
1619        }
1620        id_Delete(&(syzstr->minres[i]),r);
1621      }
1622      omFreeSize((ADDRESS)syzstr->minres,(syzstr->length+1)*sizeof(ideal));
1623    }
1624    if (syzstr->fullres!=NULL)
1625    {
1626      for (i=0;i<syzstr->length;i++)
1627      {
1628        if (syzstr->fullres[i]!=NULL)
1629        {
1630          for (j=0;j<IDELEMS(syzstr->fullres[i]);j++)
1631          {
1632            if (syzstr->fullres[i]->m[j]!=NULL)
1633              p_Delete(&(syzstr->fullres[i]->m[j]),r);
1634          }
1635        }
1636        id_Delete(&(syzstr->fullres[i]),r);
1637      }
1638      omFreeSize((ADDRESS)syzstr->fullres,(syzstr->length+1)*sizeof(ideal));
1639    }
1640    if (syzstr->weights!=0)
1641    {
1642      for (i=0;i<syzstr->length;i++)
1643      {
1644        if (syzstr->weights[i]!=NULL)
1645        {
1646          delete syzstr->weights[i];
1647        }
1648      }
1649      omFreeSize((ADDRESS)syzstr->weights,syzstr->length*sizeof(intvec*));
1650    }
1651
1652    ring sr=syzstr->syRing;
1653    if (sr==NULL) sr=r;
1654
1655    if (syzstr->resPairs!=NULL)
1656    {
1657      for (i=0;i<syzstr->length;i++)
1658      {
1659        for (j=0;j<(*syzstr->Tl)[i];j++)
1660        {
1661          if ((syzstr->resPairs[i])[j].lcm!=NULL)
1662            p_Delete(&((syzstr->resPairs[i])[j].lcm),sr);
1663          if ((i>0) && ((syzstr->resPairs[i])[j].syz!=NULL))
1664            p_Delete(&((syzstr->resPairs[i])[j].syz),sr);
1665        }
1666        if (syzstr->orderedRes[i]!=NULL)
1667        {
1668          for (j=0;j<IDELEMS(syzstr->orderedRes[i]);j++)
1669          {
1670            syzstr->orderedRes[i]->m[j] = NULL;
1671          }
1672        }
1673        id_Delete(&(syzstr->orderedRes[i]),sr);
1674        if (syzstr->truecomponents[i]!=NULL)
1675        {
1676          omFreeSize((ADDRESS)syzstr->truecomponents[i],(IDELEMS(syzstr->res[i])+1)*sizeof(int));
1677          syzstr->truecomponents[i]=NULL;
1678          omFreeSize((ADDRESS)syzstr->ShiftedComponents[i],(IDELEMS(syzstr->res[i])+1)*sizeof(long));
1679          syzstr->ShiftedComponents[i]=NULL;
1680        }
1681        if (syzstr->backcomponents[i]!=NULL)
1682        {
1683          omFreeSize((ADDRESS)syzstr->backcomponents[i],(IDELEMS(syzstr->res[i])+1)*sizeof(int));
1684          syzstr->backcomponents[i]=NULL;
1685        }
1686        if (syzstr->Howmuch[i]!=NULL)
1687        {
1688          omFreeSize((ADDRESS)syzstr->Howmuch[i],(IDELEMS(syzstr->res[i])+1)*sizeof(int));
1689          syzstr->Howmuch[i]=NULL;
1690        }
1691        if (syzstr->Firstelem[i]!=NULL)
1692        {
1693          omFreeSize((ADDRESS)syzstr->Firstelem[i],(IDELEMS(syzstr->res[i])+1)*sizeof(int));
1694          syzstr->Firstelem[i]=NULL;
1695        }
1696        if (syzstr->elemLength[i]!=NULL)
1697        {
1698          omFreeSize((ADDRESS)syzstr->elemLength[i],(IDELEMS(syzstr->res[i])+1)*sizeof(int));
1699          syzstr->elemLength[i]=NULL;
1700        }
1701        if (syzstr->res[i]!=NULL)
1702        {
1703          for (j=0;j<IDELEMS(syzstr->res[i]);j++)
1704          {
1705            if (syzstr->res[i]->m[j]!=NULL)
1706              p_Delete(&(syzstr->res[i]->m[j]),sr);
1707          }
1708        }
1709        if ((syzstr->hilb_coeffs!=NULL)
1710        && (syzstr->hilb_coeffs[i]!=NULL))
1711          delete syzstr->hilb_coeffs[i];
1712        if (syzstr->sev[i] != NULL)
1713          omFreeSize((ADDRESS)syzstr->sev[i], (IDELEMS(syzstr->res[i])+1)*sizeof(unsigned long));
1714        id_Delete(&(syzstr->res[i]),sr);
1715        if (syzstr->resPairs[i] != NULL) // OB: ????
1716          omFreeSize((ADDRESS)syzstr->resPairs[i],(*syzstr->Tl)[i]*sizeof(SObject));
1717      }
1718      omFreeSize((ADDRESS)syzstr->resPairs,syzstr->length*sizeof(SObject*));
1719      omFreeSize((ADDRESS)syzstr->res,(syzstr->length+1)*sizeof(ideal));
1720      omFreeSize((ADDRESS)syzstr->orderedRes,(syzstr->length+1)*sizeof(ideal));
1721      omFreeSize((ADDRESS)syzstr->elemLength,(syzstr->length+1)*sizeof(int*));
1722      omFreeSize((ADDRESS)syzstr->truecomponents,(syzstr->length+1)*sizeof(int*));
1723      omFreeSize((ADDRESS)syzstr->ShiftedComponents,(syzstr->length+1)*sizeof(long*));
1724      if (syzstr->sev != NULL)
1725        omFreeSize(((ADDRESS)syzstr->sev), (syzstr->length+1)*sizeof(unsigned long*));
1726      omFreeSize((ADDRESS)syzstr->backcomponents,(syzstr->length+1)*sizeof(int*));
1727      omFreeSize((ADDRESS)syzstr->Howmuch,(syzstr->length+1)*sizeof(int*));
1728      omFreeSize((ADDRESS)syzstr->Firstelem,(syzstr->length+1)*sizeof(int*));
1729      if (syzstr->hilb_coeffs!=NULL)
1730        omFreeSize((ADDRESS)syzstr->hilb_coeffs,(syzstr->length+1)*sizeof(intvec*));
1731    }
1732    if (syzstr->cw!=NULL)
1733      delete syzstr->cw;
1734    if (syzstr->betti!=NULL)
1735      delete syzstr->betti;
1736    if (syzstr->resolution!=NULL)
1737      delete syzstr->resolution;
1738    if (syzstr->Tl!=NULL)
1739      delete syzstr->Tl;
1740    if ((syzstr->syRing != NULL) && (syzstr->syRing != r))
1741    {
1742      rKill(syzstr->syRing);
1743    }
1744    omFreeSize((ADDRESS)syzstr, sizeof(ssyStrategy));
1745  }
1746}
1747
1748/*3
1749* read out the Betti numbers from resolution
1750* (if not LaScala calls the traditional Betti procedure)
1751*/
1752intvec * syBettiOfComputation(syStrategy syzstr, BOOLEAN minim,int * row_shift)
1753{
1754  int dummy;
1755  if (syzstr->betti!=NULL)
1756  {
1757    if (minim || (syzstr->resPairs!=NULL))
1758      return ivCopy(syzstr->betti);
1759  }
1760  intvec *result;
1761  if (syzstr->resPairs!=NULL)
1762  {
1763    int i,j=-1,jj=-1,l,sh=0;
1764    SRes rP=syzstr->resPairs;
1765
1766    if (syzstr->hilb_coeffs!=NULL) sh = 1;
1767    l = syzstr->length;
1768    while ((l>0) && (rP[l-1]==NULL)) l--;
1769    if (l==0) return NULL;
1770    l--;
1771    while (l>=sh)
1772    {
1773      i = 0;
1774      while ((i<(*syzstr->Tl)[l]) &&
1775        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)))
1776      {
1777        if (rP[l][i].isNotMinimal==NULL)
1778        {
1779          if (j<rP[l][i].order-l)
1780            j = rP[l][i].order-l;
1781          if (jj<l)
1782            jj = l;
1783        }
1784        i++;
1785      }
1786      l--;
1787    }
1788    j = j+sh;
1789    jj = jj+2;
1790    result=new intvec(j,jj-sh,0);
1791    IMATELEM(*result,1,1)
1792      = max(1,idRankFreeModule(syzstr->res[1],
1793                               (syzstr->syRing!=NULL?syzstr->syRing:currRing)));
1794    for (i=sh;i<jj;i++)
1795    {
1796      j = 0;
1797      while (j<(*syzstr->Tl)[i])
1798      {
1799        if ((rP[i][j].isNotMinimal==NULL)&&
1800           ((rP[i][j].lcm!=NULL) || (rP[i][j].syz!=NULL)))
1801          IMATELEM(*result,rP[i][j].order-i+sh,i+2-sh)++;
1802        j++;
1803      }
1804    }
1805  }
1806  else if (syzstr->fullres!=NULL)
1807    result = syBetti(syzstr->fullres,syzstr->length,&dummy,NULL,minim,row_shift);
1808  else
1809    result = syBetti(syzstr->minres,syzstr->length,&dummy,NULL,minim,row_shift);
1810  if ((minim) || (syzstr->resPairs!=NULL))
1811    syzstr->betti = ivCopy(result);
1812  return result;
1813}
1814
1815/*2
1816* read out the Betti numbers from resolution
1817* (interpreter interface)
1818*/
1819BOOLEAN syBetti2(leftv res, leftv u, leftv w)
1820{
1821  syStrategy syzstr=(syStrategy)u->Data();
1822  BOOLEAN minim=(int)w->Data();
1823  int row_shift=0;
1824
1825  res->data=(void *)syBettiOfComputation(syzstr,minim,&row_shift);
1826  atSet(res,omStrDup("rowShift"),(void*)row_shift,INT_CMD);
1827  return FALSE;
1828}
1829BOOLEAN syBetti1(leftv res, leftv u)
1830{
1831  syStrategy syzstr=(syStrategy)u->Data();
1832
1833  res->data=(void *)syBettiOfComputation(syzstr);
1834  return FALSE;
1835}
1836
1837/*3
1838* computes the allocated length of the resolution
1839*/
1840int syLength(syStrategy syzstr)
1841{
1842  return syzstr->length;
1843}
1844
1845/*3
1846* computes the real length of the resolution
1847*/
1848int sySize(syStrategy syzstr)
1849{
1850  resolvente r=syzstr->res;
1851  if (r==NULL)
1852    r = syzstr->fullres;
1853  if (r==NULL)
1854    r = syzstr->minres;
1855  if (r==NULL)
1856  {
1857    WerrorS("No resolution found");
1858    return 0;
1859  }
1860  int i=syzstr->length;
1861  while ((i>0) && (r[i-1]==NULL)) i--;
1862  return i;
1863}
1864
1865/*3
1866* computes the cohomological dimension of res[1]
1867*/
1868int syDim(syStrategy syzstr)
1869{
1870  int i,j=-1,l;
1871  if (syzstr->resPairs!=NULL)
1872  {
1873    SRes rP=syzstr->resPairs;
1874
1875    l = syzstr->length;
1876    while ((l>0) && (rP[l-1]==NULL)) l--;
1877    if (l==0) return -1;
1878    l--;
1879    while (l>=0)
1880    {
1881      i = 0;
1882      while ((i<(*syzstr->Tl)[l]) &&
1883        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)) &&
1884        (rP[l][i].isNotMinimal!=NULL))
1885      {
1886        i++;
1887      }
1888      if ((i<(*syzstr->Tl)[l]) &&
1889        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)) &&
1890        (rP[l][i].isNotMinimal==NULL))
1891        return l;
1892      l--;
1893    }
1894    return l;
1895  }
1896  else
1897    return sySize(syzstr);
1898}
1899
1900/*3
1901* copies the resolution (by increment the reference counter)
1902*/
1903syStrategy syCopy(syStrategy syzstr)
1904{
1905  syStrategy result=syzstr;
1906  (result->references)++;
1907  return result;
1908}
1909
1910/*2
1911* local print procedure used in syPrint
1912*/
1913static void syPrintEmptySpaces(int i)
1914{
1915  if (i!=0)
1916  {
1917    Print(" ");
1918    syPrintEmptySpaces(i/10);
1919  }
1920}
1921
1922/*2
1923* local print procedure used in syPrint
1924*/
1925static void syPrintEmptySpaces1(int i)
1926{
1927  if (i!=0)
1928  {
1929    Print(" ");
1930    syPrintEmptySpaces1(i-1);
1931  }
1932}
1933
1934/*2
1935* local print procedure used in syPrint
1936*/
1937static int syLengthInt(int i)
1938{
1939  int j=0;
1940
1941  if (i==0) return 1;
1942  while (i!=0)
1943  {
1944    j++;
1945    i = i/10;
1946  }
1947  return j;
1948}
1949
1950/*3
1951* prints the resolution as sequence of free modules
1952*/
1953void syPrint(syStrategy syzstr)
1954{
1955  if ((syzstr->resPairs==NULL) && (syzstr->fullres==NULL)
1956     && (syzstr->minres==NULL))
1957  {
1958    Print("No resolution defined\n");
1959    return;
1960  }
1961  int l=0;
1962  if (syzstr->resolution==NULL)
1963  {
1964    int j;
1965    if (syzstr->resPairs!=NULL)
1966    {
1967      syzstr->resolution = new intvec(syzstr->length+1);
1968      SRes rP=syzstr->resPairs;
1969      (*syzstr->resolution)[0]
1970        = max(1,idRankFreeModule(syzstr->res[1],
1971                                 (syzstr->syRing != NULL ? syzstr->syRing : currRing)));
1972      while ((l<syzstr->length) && (rP[l]!=NULL))
1973      {
1974        j=0;
1975        while ((j<(*syzstr->Tl)[l]) &&
1976          ((rP[l][j].lcm!=NULL) || (rP[l][j].syz!=NULL)))
1977        {
1978          if (rP[l][j].isNotMinimal==NULL)
1979            ((*syzstr->resolution)[l+1])++;
1980          j++;
1981        }
1982        l++;
1983      }
1984    }
1985    else
1986    {
1987      resolvente rr;
1988      syzstr->resolution = new intvec(syzstr->length+2);
1989      if (syzstr->minres!=NULL)
1990        rr = syzstr->minres;
1991      else
1992        rr = syzstr->fullres;
1993      (*syzstr->resolution)[0]
1994        = max(1,idRankFreeModule(rr[0],
1995                                 (syzstr->syRing != NULL ? syzstr->syRing : currRing)));
1996      while ((l<syzstr->length) && (rr[l]!=NULL))
1997      {
1998        j = IDELEMS(rr[l]);
1999        while ((j>0) && (rr[l]->m[j-1]==NULL)) j--;
2000        ((*syzstr->resolution)[l+1]) = j;
2001        l++;
2002      }
2003    }
2004  }
2005  char *sn=currRingHdl->id;
2006  int sl=strlen(sn);
2007  syPrintEmptySpaces1(sl);
2008  l = 0;
2009  loop
2010  {
2011    Print("%d",(*syzstr->resolution)[l]);
2012    if ((l>=syzstr->resolution->length()) || ((*syzstr->resolution)[l]==0))
2013      break;
2014    syPrintEmptySpaces1(sl+5);
2015    l++;
2016  }
2017  PrintLn();
2018  l = 0;
2019  loop
2020  {
2021    Print(sn);
2022    if ((l>=syzstr->resolution->length()) || ((*syzstr->resolution)[l]==0))
2023      break;
2024    syPrintEmptySpaces((*syzstr->resolution)[l]);
2025    Print(" <-- ");
2026    l++;
2027  }
2028  PrintLn();
2029  PrintLn();
2030  l = 0;
2031  loop
2032  {
2033    Print("%d",l);
2034    if ((l>=syzstr->resolution->length()) || ((*syzstr->resolution)[l]==0))
2035      break;
2036    syPrintEmptySpaces1(sl+5+syLengthInt((*syzstr->resolution)[l])-
2037                         syLengthInt(l));
2038    l++;
2039  }
2040  PrintLn();
2041  if (syzstr->minres==NULL)
2042  {
2043    Print("resolution not minimized yet");
2044    PrintLn();
2045  }
2046}
2047
2048/*2
2049* divides out the weight monomials (given by the Schreyer-ordering)
2050* from the LaScala-resolution
2051*/
2052static resolvente syReorder(resolvente res,int length,
2053        syStrategy syzstr,BOOLEAN toCopy=TRUE,resolvente totake=NULL)
2054{
2055  int i,j,l;
2056  poly p,q,tq;
2057  polyset ri1;
2058  resolvente fullres;
2059  ring origR=syzstr->syRing;
2060  fullres = (resolvente)omAlloc0((length+1)*sizeof(ideal));
2061  if (totake==NULL)
2062    totake = res;
2063  for (i=length-1;i>0;i--)
2064  {
2065    if (res[i]!=NULL)
2066    {
2067      if (i>1)
2068      {
2069        j = IDELEMS(res[i-1]);
2070        while ((j>0) && (res[i-1]->m[j-1]==NULL)) j--;
2071        fullres[i-1] = idInit(IDELEMS(res[i]),j);
2072        ri1 = totake[i-1]->m;
2073        for (j=IDELEMS(res[i])-1;j>=0;j--)
2074        {
2075          p = res[i]->m[j];
2076          q = NULL;
2077          while (p!=NULL)
2078          {
2079            if (toCopy)
2080            {
2081              if (origR!=NULL)
2082                tq = prHeadR(p,origR);
2083              else
2084                tq = pHead(p);
2085              pIter(p);
2086            }
2087            else
2088            {
2089              res[i]->m[j] = NULL;
2090              if (origR!=NULL)
2091              {
2092                poly pp=p;
2093                pIter(p);
2094                pNext(pp)=NULL;
2095                tq = prMoveR(pp, origR);
2096              }
2097              else
2098              {
2099                tq = p;
2100                pIter(p);
2101                pNext(tq) = NULL;
2102              }
2103            }
2104//            pWrite(tq);
2105            pTest(tq);
2106            for (l=pVariables;l>0;l--)
2107            {
2108              if (origR!=NULL)
2109                pSubExp(tq,l, p_GetExp(ri1[pGetComp(tq)-1],l,origR));
2110              else
2111                pSubExp(tq,l, pGetExp(ri1[pGetComp(tq)-1],l));
2112            }
2113            pSetm(tq);
2114            pTest(tq);
2115            q = pAdd(q,tq);
2116            pTest(q);
2117          }
2118          fullres[i-1]->m[j] = q;
2119        }
2120      }
2121      else
2122      {
2123        if (origR!=NULL)
2124        {
2125          fullres[i-1] = idInit(IDELEMS(res[i]),res[i]->rank);
2126          for (j=IDELEMS(res[i])-1;j>=0;j--)
2127          {
2128            if (toCopy)
2129              fullres[i-1]->m[j] = prCopyR(res[i]->m[j], origR);
2130            else
2131            {
2132              fullres[i-1]->m[j] = prMoveR(res[i]->m[j], origR);
2133              res[i]->m[j] = NULL;
2134            }
2135          }
2136        }
2137        else
2138        {
2139          if (toCopy)
2140            fullres[i-1] = idCopy(res[i]);
2141          else
2142          {
2143            fullres[i-1] = res[i];
2144            res[i] = NULL;
2145          }
2146        }
2147        for (j=IDELEMS(fullres[i-1])-1;j>=0;j--)
2148          fullres[i-1]->m[j] = pSortCompCorrect(fullres[i-1]->m[j]);
2149      }
2150      if (!toCopy)
2151      {
2152        if (res[i]!=NULL) idDelete(&res[i]);
2153      }
2154    }
2155  }
2156  if (!toCopy)
2157    omFreeSize((ADDRESS)res,(length+1)*sizeof(ideal));
2158  //syzstr->length = length;
2159  return fullres;
2160}
2161
2162/*3
2163* converts a resolution into a list of modules
2164*/
2165lists syConvRes(syStrategy syzstr,BOOLEAN toDel)
2166{
2167  if ((syzstr->fullres==NULL) && (syzstr->minres==NULL))
2168  {
2169    if (syzstr->hilb_coeffs==NULL)
2170    {
2171      syzstr->fullres = syReorder(syzstr->res,syzstr->length,syzstr);
2172    }
2173    else
2174    {
2175      syzstr->minres = syReorder(syzstr->orderedRes,syzstr->length,syzstr);
2176      syKillEmptyEntres(syzstr->minres,syzstr->length);
2177    }
2178  }
2179  resolvente tr;
2180  int typ0=IDEAL_CMD;
2181  if (syzstr->minres!=NULL)
2182    tr = syzstr->minres;
2183  else
2184    tr = syzstr->fullres;
2185  resolvente trueres=NULL;
2186  intvec ** w=NULL;
2187  if (syzstr->length>0)
2188  {
2189    trueres=(resolvente)omAlloc0((syzstr->length)*sizeof(ideal));
2190    for (int i=(syzstr->length)-1;i>=0;i--)
2191    {
2192      if (tr[i]!=NULL)
2193      {
2194        trueres[i] = idCopy(tr[i]);
2195      }
2196    }
2197    if (idRankFreeModule(trueres[0]) > 0)
2198      typ0 = MODUL_CMD;
2199    if (syzstr->weights!=NULL)
2200    {
2201      w = (intvec**)omAlloc0((syzstr->length)*sizeof(intvec*));
2202      for (int i=(syzstr->length)-1;i>=0;i--)
2203      {
2204        if (syzstr->weights[i]!=NULL) w[i] = ivCopy(syzstr->weights[i]);
2205      }
2206    }
2207  }
2208  lists li = liMakeResolv(trueres,syzstr->length,syzstr->list_length,typ0,w);
2209  if (w != NULL) omFreeSize(w, (syzstr->length)*sizeof(intvec*));
2210  if (toDel) syKillComputation(syzstr);
2211  return li;
2212}
2213
2214/*3
2215* converts a list of modules into a resolution
2216*/
2217syStrategy syConvList(lists li,BOOLEAN toDel)
2218{
2219  int typ0;
2220  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2221
2222  resolvente fr = liFindRes(li,&(result->length),&typ0,&(result->weights));
2223  if (fr != NULL)
2224  {
2225
2226    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
2227    for (int i=result->length-1;i>=0;i--)
2228    {
2229      if (fr[i]!=NULL)
2230        result->fullres[i] = idCopy(fr[i]);
2231    }
2232    result->list_length=result->length;
2233    omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
2234  }
2235  else
2236  {
2237    omFreeSize(result, sizeof(ssyStrategy));
2238    result = NULL;
2239  }
2240  if (toDel) li->Clean();
2241  return result;
2242}
2243
2244/*3
2245* converts a list of modules into a minimal resolution
2246*/
2247syStrategy syForceMin(lists li)
2248{
2249  int typ0;
2250  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2251
2252  resolvente fr = liFindRes(li,&(result->length),&typ0);
2253  result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
2254  for (int i=result->length-1;i>=0;i--)
2255  {
2256    if (fr[i]!=NULL)
2257      result->minres[i] = idCopy(fr[i]);
2258  }
2259  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
2260  return result;
2261}
2262
2263/*2
2264* deleting all monomials the component of which correspond
2265* to non-minimal generators
2266*/
2267static poly syStripOut(poly p,intvec * toStrip)
2268{
2269  if (toStrip==NULL) return p;
2270  poly pp=p;
2271
2272  while ((pp!=NULL) && ((*toStrip)[pGetComp(pp)]!=0))
2273    pDeleteLm(&pp);
2274  p = pp;
2275  if (pp!=NULL)
2276  {
2277    while (pNext(pp)!=NULL)
2278    {
2279      if ((*toStrip)[pGetComp(pNext(pp))]!=0)
2280        pDeleteLm(&pNext(pp));
2281      else
2282        pIter(pp);
2283    }
2284  }
2285  return p;
2286}
2287
2288/*2
2289* copies only those monomials the component of which correspond
2290* to minimal generators
2291*/
2292static poly syStripOutCopy(poly p,intvec * toStrip)
2293{
2294  if (toStrip==NULL) return pCopy(p);
2295  poly result=NULL,pp;
2296
2297  while (p!=NULL)
2298  {
2299    if ((*toStrip)[pGetComp(p)]==0)
2300    {
2301      if (result==NULL)
2302      {
2303        result = pp = pHead(p);
2304      }
2305      else
2306      {
2307        pNext(pp) = pHead(p);
2308        pIter(pp);
2309      }
2310    }
2311    pIter(p);
2312  }
2313  return result;
2314}
2315
2316/*2
2317* minimizes toMin
2318*/
2319static poly syMinimizeP(int toMin,syStrategy syzstr,intvec * ordn,int index,
2320                        intvec * toStrip)
2321{
2322  int ii=0,i,j,tc;
2323  poly p,pp,q=NULL,tq,pisN;
2324  SSet sPairs=syzstr->resPairs[index];
2325  poly tempStripped=NULL;
2326
2327  //pp=pCopy(syzstr->res[index+1]->m[toMin]);
2328  pp = syStripOutCopy(syzstr->res[index+1]->m[toMin],toStrip);
2329  while ((ii<ordn->length()) && ((*ordn)[ii]!=-1) &&
2330             (sPairs[(*ordn)[ii]].syzind!=toMin))
2331  {
2332    ii++;
2333  }
2334  while (ii>=0)
2335  {
2336    i = (*ordn)[ii];
2337    if (sPairs[i].isNotMinimal!=NULL)
2338    {
2339      tempStripped =
2340        syStripOutCopy(syzstr->res[index+1]->m[sPairs[i].syzind],toStrip);
2341      pisN = sPairs[i].isNotMinimal;
2342      tc = pGetComp(pisN);
2343      p = pp;
2344      while (p!=NULL)
2345      {
2346        if (pGetComp(p)==tc)
2347        {
2348          tq = pInit();
2349          for(j=pVariables; j>0; j--)
2350            pSetExp(tq,j, pGetExp(p,j)-pGetExp(pisN,j));
2351          pSetComp(tq, 0);
2352          pSetCoeff0(tq,nDiv(pGetCoeff(p),pGetCoeff(pisN)));
2353          pGetCoeff(tq) = nNeg(pGetCoeff(tq));
2354          pSetm(tq);
2355          q = pAdd(q,pMult_mm(pCopy(tempStripped),tq));
2356          pDelete(&tq);
2357        }
2358        pIter(p);
2359      }
2360      if (q!=NULL)
2361      {
2362        pp = pAdd(pp,q);
2363        q = NULL;
2364      }
2365      pDelete(&tempStripped);
2366    }
2367    ii--;
2368  }
2369  return pp;
2370}
2371
2372/*2
2373* minimizes toMin
2374*/
2375static poly syMinimizeP1(int toMin,syStrategy syzstr,intvec * ordn,int index,
2376                        intvec * toStrip)
2377{
2378  int ii=0,i,j,tc,lp,ltS=-1;
2379  poly p,mp=NULL,pp,q=NULL,tq,pisN;
2380  SSet sPairs=syzstr->resPairs[index];
2381  poly tempStripped=NULL;
2382
2383  pp = syStripOutCopy(syzstr->res[index+1]->m[toMin],toStrip);
2384  kBucketInit(syzstr->bucket,pp,-1);
2385  while ((ii<ordn->length()) && ((*ordn)[ii]!=-1) &&
2386             (sPairs[(*ordn)[ii]].syzind!=toMin))
2387  {
2388    ii++;
2389  }
2390  while (ii>=0)
2391  {
2392    i = (*ordn)[ii];
2393    if (sPairs[i].isNotMinimal!=NULL)
2394    {
2395      tempStripped =
2396        syStripOutCopy(syzstr->res[index+1]->m[sPairs[i].syzind],toStrip);
2397      tc = pGetComp(sPairs[i].isNotMinimal);
2398      //p = pTakeOutComp1(&tempStripped,tc);
2399      int lu;
2400      pTakeOutComp(&tempStripped,tc,&p,&lu);
2401      kBucketTakeOutComp(syzstr->bucket,tc,&mp,&lp);
2402      mp = pDivideM(mp,p);
2403      while (mp!=NULL)
2404      {
2405        p = pNext(mp);
2406        pNext(mp) = NULL;
2407        ltS = -1;
2408        kBucket_Minus_m_Mult_p(syzstr->bucket,mp,tempStripped,&ltS);
2409        mp = p;
2410      }
2411      pDelete(&mp);
2412      pDelete(&tempStripped);
2413    }
2414    ii--;
2415  }
2416  kBucketClear(syzstr->bucket,&pp,&lp);
2417  return pp;
2418}
2419
2420/*2
2421* deletes empty components after minimization
2422*/
2423void syKillEmptyEntres(resolvente res,int length)
2424{
2425  int i,j,jj,k,rj;
2426  intvec * changes;
2427  poly p;
2428  ideal ri;
2429
2430  for (i=0;i<length;i++)
2431  {
2432    ri = res[i];
2433    if (ri!=NULL)
2434    {
2435      rj = IDELEMS(ri);
2436      changes = new intvec(rj+1,1,-1);
2437      while ((rj>0) && (ri->m[rj-1]==NULL)) rj--;
2438      j = k = 0;
2439      while (j+k<rj)
2440      {
2441        if (ri->m[j+k]!=NULL)
2442        {
2443          ri->m[j] = ri->m[j+k];
2444          (*changes)[j+k+1] = j+1;
2445          j++;
2446        }
2447        else
2448        {
2449          k++;
2450        }
2451      }
2452      for (jj=j;jj<rj;jj++)
2453        ri->m[jj] = NULL;
2454      if (res[i+1]!=NULL)
2455      {
2456        ri = res[i+1];
2457        for (j=IDELEMS(ri)-1;j>=0;j--)
2458        {
2459          p = ri->m[j];
2460          while (p!=NULL)
2461          {
2462            pSetComp(p,(*changes)[pGetComp(p)]);
2463            pSetm(p);
2464            pIter(p);
2465          }
2466        }
2467      }
2468      delete changes;
2469    }
2470  }
2471}
2472
2473/*2
2474* determines the components for minimization
2475*/
2476static intvec * syToStrip(syStrategy syzstr, int index)
2477{
2478  intvec * result=NULL;
2479
2480  if ((syzstr->resPairs[index-1]!=NULL) && (!idIs0(syzstr->res[index])))
2481  {
2482    result=new intvec(IDELEMS(syzstr->res[index])+1);
2483    for (int i=(*syzstr->Tl)[index-1]-1;i>=0;i--)
2484    {
2485      if (syzstr->resPairs[index-1][i].isNotMinimal!=NULL)
2486      {
2487        (*result)[syzstr->resPairs[index-1][i].syzind+1] = 1;
2488      }
2489    }
2490  }
2491  return result;
2492}
2493
2494/*2
2495* re-computes the order of pairs during the algorithm
2496* this ensures to procede with a triangular matrix
2497*/
2498static intvec * syOrdPairs(SSet sPairs, int length)
2499{
2500  intvec * result=new intvec(length,1,-1);
2501  int i,j=0,k=-1,l,ii;
2502
2503  loop
2504  {
2505    l = -1;
2506    for(i=0;i<length;i++)
2507    {
2508      if (sPairs[i].syzind>k)
2509      {
2510        if (l==-1)
2511        {
2512          l = sPairs[i].syzind;
2513          ii = i;
2514        }
2515        else
2516        {
2517          if (sPairs[i].syzind<l)
2518          {
2519            l = sPairs[i].syzind;
2520            ii = i;
2521          }
2522        }
2523      }
2524    }
2525    if (l==-1) break;
2526    (*result)[j] = ii;
2527    j++;
2528    k = l;
2529  }
2530  return result;
2531}
2532
2533/*2
2534* minimizes the output of LaScala
2535*/
2536static resolvente syReadOutMinimalRes(syStrategy syzstr,
2537           BOOLEAN computeStd=FALSE)
2538{
2539  intvec * Strip, * ordn;
2540  resolvente tres=(resolvente)omAlloc0((syzstr->length+1)*sizeof(ideal));
2541  ring tmpR = NULL;
2542  ring origR = currRing;
2543
2544//Print("Hier ");
2545  if (computeStd)
2546  {
2547    tres[0] = syzstr->res[1];
2548    syzstr->res[1] = idInit(IDELEMS(tres[0]),tres[0]->rank);
2549    return tres;
2550  }
2551  int i,j,l,index,ii,i1;
2552  poly p;
2553  ideal rs;
2554  SSet sPairs;
2555  int * ord,*b0,*b1;
2556
2557  assume(syzstr->syRing != NULL);
2558  rChangeCurrRing(syzstr->syRing);
2559//Print("laeufts ");
2560  syzstr->bucket = kBucketCreate();
2561  for (index=syzstr->length-1;index>0;index--)
2562  {
2563    if (syzstr->resPairs[index]!=NULL)
2564    {
2565//Print("ideal %d: \n",index);
2566      currcomponents = syzstr->truecomponents[index];
2567      currShiftedComponents = syzstr->ShiftedComponents[index];
2568      rChangeSComps(currcomponents, currShiftedComponents,
2569                    IDELEMS(syzstr->res[index]));
2570      sPairs = syzstr->resPairs[index];
2571      Strip = syToStrip(syzstr,index);
2572      tres[index+1] = idInit(IDELEMS(syzstr->res[index+1]),syzstr->res[index+1]->rank);
2573      i1 = (*syzstr->Tl)[index];
2574//Print("i1= %d\n",i1);
2575      ordn = syOrdPairs(sPairs,i1);
2576      for (i=0;i<i1;i++)
2577      {
2578        if ((sPairs[i].isNotMinimal==NULL) && (sPairs[i].lcm!=NULL))
2579        {
2580          l = sPairs[i].syzind;
2581//Print("Minimiere Poly %d: ",l);pWrite(syzstr->res[index+1]->m[l]);
2582          tres[index+1]->m[l] =
2583            syMinimizeP1(l,syzstr,ordn,index,Strip);
2584        }
2585      }
2586      delete Strip;
2587      Strip = NULL;
2588    }
2589  }
2590  currcomponents = syzstr->truecomponents[0];
2591  currShiftedComponents = syzstr->ShiftedComponents[0];
2592  rChangeSComps( currcomponents, currShiftedComponents,
2593                 IDELEMS(syzstr->res[0]));
2594  tres[1] = idInit(IDELEMS(syzstr->res[1]),syzstr->res[1]->rank);
2595  sPairs = syzstr->resPairs[0];
2596  for (i=(*syzstr->Tl)[0]-1;i>=0;i--)
2597  {
2598    if (sPairs[i].syzind>=0)
2599    {
2600      tres[1]->m[sPairs[i].syzind] = pCopy(syzstr->res[1]->m[sPairs[i].syzind]);
2601    }
2602  }
2603/*--- changes to the original ring------------------*/
2604  kBucketDestroy(&syzstr->bucket);
2605  if (syzstr->syRing != NULL)
2606  {
2607    rChangeCurrRing(origR);
2608    // Thomas: now make sure that all data which you need is pFetchCopied
2609    // maybe incoporate it into syReorder ??
2610  }
2611  tres = syReorder(tres,syzstr->length,syzstr,FALSE,syzstr->res);
2612  syKillEmptyEntres(tres,syzstr->length);
2613  idSkipZeroes(tres[0]);
2614  return tres;
2615}
2616
2617/*3
2618* minimizes any kind of resolution
2619*/
2620syStrategy syMinimize(syStrategy syzstr)
2621{
2622  if (syzstr->minres==NULL)
2623  {
2624    if (syzstr->resPairs!=NULL)
2625    {
2626      if (syzstr->hilb_coeffs==NULL)
2627      {
2628        syzstr->minres = syReadOutMinimalRes(syzstr);
2629      }
2630      else
2631      {
2632        syzstr->minres = syReorder(syzstr->orderedRes,syzstr->length,syzstr);
2633      }
2634    }
2635    else if (syzstr->fullres!=NULL)
2636    {
2637      syMinimizeResolvente(syzstr->fullres,syzstr->length,1);
2638      syzstr->minres = syzstr->fullres;
2639      syzstr->fullres = NULL;
2640    }
2641  }
2642  (syzstr->references)++;
2643  return syzstr;
2644}
2645
2646/*2
2647* implementation of LaScala's algorithm
2648* assumes that the given module is homogeneous
2649* works with slanted degree, uses syChosePairs
2650*/
2651syStrategy syLaScala3(ideal arg,int * length)
2652{
2653  BOOLEAN noPair=FALSE;
2654  int i,j,actdeg=32000,index=0,reg=-1;
2655  int startdeg,howmuch;
2656  poly p;
2657  ideal temp;
2658  SSet nextPairs;
2659  syStrategy syzstr=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2660  ring origR = currRing;
2661
2662  if ((idIs0(arg)) ||
2663      ((idRankFreeModule(arg)>0) && (!idHomModule(arg,NULL,&(syzstr->cw)))))
2664  {
2665    syzstr->minres = (resolvente)omAlloc0Bin(ideal_bin);
2666    syzstr->length = 1;
2667    syzstr->minres[0] = idInit(1,arg->rank);
2668    return syzstr;
2669  }
2670
2671  //crit = 0;
2672  //euler = -1;
2673  redpol = pInit();
2674  syzstr->length = *length = pVariables+2;
2675
2676  // Creare dp,S ring and change to it
2677  syzstr->syRing = rCurrRingAssure_dp_S();
2678  assume(syzstr->syRing != origR);
2679
2680  // set initial ShiftedComps
2681  currcomponents = (int*)omAlloc0((arg->rank+1)*sizeof(int));
2682  currShiftedComponents = (long*)omAlloc0((arg->rank+1)*sizeof(long));
2683  for (i=0;i<=arg->rank;i++)
2684  {
2685    currShiftedComponents[i] = (i)*SYZ_SHIFT_BASE;
2686    currcomponents[i] = i;
2687  }
2688  rChangeSComps(currcomponents, currShiftedComponents, arg->rank);
2689/*--- initializes the data structures---------------*/
2690  syzstr->Tl = new intvec(*length);
2691  temp = idInit(IDELEMS(arg),arg->rank);
2692  for (i=0;i<IDELEMS(arg);i++)
2693  {
2694    temp->m[i] = prCopyR( arg->m[i], origR);
2695    if (temp->m[i]!=NULL)
2696    {
2697      j = pTotaldegree(temp->m[i]);
2698      if (j<actdeg) actdeg = j;
2699    }
2700  }
2701  idTest(temp);
2702  idSkipZeroes(temp);
2703  idTest(temp);
2704  syzstr->resPairs = syInitRes(temp,length,syzstr->Tl,syzstr->cw);
2705  omFreeSize((ADDRESS)currcomponents,(arg->rank+1)*sizeof(int));
2706  omFreeSize((ADDRESS)currShiftedComponents,(arg->rank+1)*sizeof(long));
2707  syzstr->res = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
2708  syzstr->orderedRes = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
2709  syzstr->elemLength = (int**)omAlloc0((*length+1)*sizeof(int*));
2710  syzstr->truecomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
2711  syzstr->ShiftedComponents = (long**)omAlloc0((*length+1)*sizeof(long*));
2712  syzstr->backcomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
2713  syzstr->Howmuch = (int**)omAlloc0((*length+1)*sizeof(int*));
2714  syzstr->Firstelem = (int**)omAlloc0((*length+1)*sizeof(int*));
2715  syzstr->sev = (unsigned long **) omAlloc0((*length+1)*sizeof(unsigned long *));
2716  syzstr->bucket = kBucketCreate();
2717  int len0=idRankFreeModule(temp)+1;
2718
2719  startdeg = actdeg;
2720  nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2721  //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2722/*--- computes the resolution ----------------------*/
2723  while (nextPairs!=NULL)
2724  {
2725    if (TEST_OPT_PROT) Print("%d",actdeg);
2726    if (TEST_OPT_PROT) Print("(m%d)",index);
2727    if (index==0)
2728      i = syInitSyzMod(syzstr,index,len0);
2729    else
2730      i = syInitSyzMod(syzstr,index);
2731    currcomponents = syzstr->truecomponents[max(index-1,0)];
2732    currShiftedComponents = syzstr->ShiftedComponents[max(index-1,0)];
2733    rChangeSComps(currcomponents, currShiftedComponents,
2734                  IDELEMS(syzstr->res[max(index-1,0)]));
2735    j = syInitSyzMod(syzstr,index+1);
2736    if (index>0)
2737    {
2738      syRedNextPairs(nextPairs,syzstr,howmuch,index);
2739      syCompactifyPairSet(syzstr->resPairs[index],(*syzstr->Tl)[index],0);
2740    }
2741    else
2742      syRedGenerOfCurrDeg(syzstr,actdeg,index+1);
2743/*--- creates new pairs -----------------------------*/
2744    syCreateNewPairs(syzstr,index,i);
2745    if (index<(*length)-1)
2746    {
2747      syCreateNewPairs(syzstr,index+1,j);
2748    }
2749    index++;
2750    nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2751    //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2752  }
2753  if (temp!=NULL) idDelete(&temp);
2754  kBucketDestroy(&(syzstr->bucket));
2755  if (origR != syzstr->syRing)
2756    rChangeCurrRing(origR);
2757  pDeleteLm(&redpol);
2758  if (TEST_OPT_PROT) PrintLn();
2759  return syzstr;
2760}
2761
Note: See TracBrowser for help on using the repository browser.