source: git/kernel/GBEngine/syz1.cc @ f8aa89

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