source: git/kernel/syz1.cc @ 1f637e

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