source: git/kernel/GBEngine/syz1.cc @ 823db5

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