source: git/kernel/syz1.cc @ f2b6f0b

spielwiese
Last change on this file since f2b6f0b was 305203, checked in by Hans Schoenemann <hannes@…>, 13 years ago
syStatistics: only for debugging git-svn-id: file:///usr/local/Singular/svn/trunk@13647 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 67.2 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 <kernel/options.h>
12#include <omalloc/omalloc.h>
13#include <kernel/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 <kernel/intvec.h>
20#include <kernel/numbers.h>
21#include <kernel/modulop.h>
22#include <kernel/ideals.h>
23#include <kernel/intvec.h>
24#include <kernel/ring.h>
25#include <kernel/syz.h>
26#include <kernel/kbuckets.h>
27#include <kernel/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 = pVariables;
210    while ((i>1) && (pGetExp(p1,i)==pGetExp(p2,i)))
211      i--;
212    //(*orderingdepth)[pVariables-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      rKill(syzstr->syRing);
1643    }
1644    omFreeSize((ADDRESS)syzstr, sizeof(ssyStrategy));
1645  }
1646}
1647
1648/*2
1649* divides out the weight monomials (given by the Schreyer-ordering)
1650* from the LaScala-resolution
1651*/
1652resolvente syReorder(resolvente res,int length,
1653        syStrategy syzstr,BOOLEAN toCopy,resolvente totake)
1654{
1655  int i,j,l;
1656  poly p,q,tq;
1657  polyset ri1;
1658  resolvente fullres;
1659  ring origR=syzstr->syRing;
1660  fullres = (resolvente)omAlloc0((length+1)*sizeof(ideal));
1661  if (totake==NULL)
1662    totake = res;
1663  for (i=length-1;i>0;i--)
1664  {
1665    if (res[i]!=NULL)
1666    {
1667      if (i>1)
1668      {
1669        j = IDELEMS(res[i-1]);
1670        while ((j>0) && (res[i-1]->m[j-1]==NULL)) j--;
1671        fullres[i-1] = idInit(IDELEMS(res[i]),j);
1672        ri1 = totake[i-1]->m;
1673        for (j=IDELEMS(res[i])-1;j>=0;j--)
1674        {
1675          p = res[i]->m[j];
1676          q = NULL;
1677          while (p!=NULL)
1678          {
1679            if (toCopy)
1680            {
1681              if (origR!=NULL)
1682                tq = prHeadR(p,origR);
1683              else
1684                tq = pHead(p);
1685              pIter(p);
1686            }
1687            else
1688            {
1689              res[i]->m[j] = NULL;
1690              if (origR!=NULL)
1691              {
1692                poly pp=p;
1693                pIter(p);
1694                pNext(pp)=NULL;
1695                tq = prMoveR(pp, origR);
1696              }
1697              else
1698              {
1699                tq = p;
1700                pIter(p);
1701                pNext(tq) = NULL;
1702              }
1703            }
1704//            pWrite(tq);
1705            pTest(tq);
1706            for (l=pVariables;l>0;l--)
1707            {
1708              if (origR!=NULL)
1709                pSubExp(tq,l, p_GetExp(ri1[pGetComp(tq)-1],l,origR));
1710              else
1711                pSubExp(tq,l, pGetExp(ri1[pGetComp(tq)-1],l));
1712            }
1713            pSetm(tq);
1714            pTest(tq);
1715            q = pAdd(q,tq);
1716            pTest(q);
1717          }
1718          fullres[i-1]->m[j] = q;
1719        }
1720      }
1721      else
1722      {
1723        if (origR!=NULL)
1724        {
1725          fullres[i-1] = idInit(IDELEMS(res[i]),res[i]->rank);
1726          for (j=IDELEMS(res[i])-1;j>=0;j--)
1727          {
1728            if (toCopy)
1729              fullres[i-1]->m[j] = prCopyR(res[i]->m[j], origR);
1730            else
1731            {
1732              fullres[i-1]->m[j] = prMoveR(res[i]->m[j], origR);
1733              res[i]->m[j] = NULL;
1734            }
1735          }
1736        }
1737        else
1738        {
1739          if (toCopy)
1740            fullres[i-1] = idCopy(res[i]);
1741          else
1742          {
1743            fullres[i-1] = res[i];
1744            res[i] = NULL;
1745          }
1746        }
1747        for (j=IDELEMS(fullres[i-1])-1;j>=0;j--)
1748          fullres[i-1]->m[j] = pSortCompCorrect(fullres[i-1]->m[j]);
1749      }
1750      if (!toCopy)
1751      {
1752        if (res[i]!=NULL) idDelete(&res[i]);
1753      }
1754    }
1755  }
1756  if (!toCopy)
1757    omFreeSize((ADDRESS)res,(length+1)*sizeof(ideal));
1758  //syzstr->length = length;
1759  return fullres;
1760}
1761
1762/*3
1763* read out the Betti numbers from resolution
1764* (if not LaScala calls the traditional Betti procedure)
1765*/
1766intvec * syBettiOfComputation(syStrategy syzstr, BOOLEAN minim,int * row_shift,
1767                              intvec* weights)
1768{
1769  int dummy;
1770  BOOLEAN std_weights=TRUE;
1771  if ((weights!=NULL)
1772  && (syzstr->betti!=NULL)
1773  && (syzstr->weights!=NULL) && (syzstr->weights[0]!=NULL))
1774  {
1775    int i;
1776    for(i=weights->length()-1; i>=0; i--)
1777    {
1778      //Print("test %d: %d - %d\n",i,(*weights)[i], (*(syzstr->weights[0]))[i]);
1779      if ((*weights)[i]!=(*(syzstr->weights[0]))[i])
1780      {
1781        std_weights=FALSE;
1782        break;
1783      }
1784    }
1785  }
1786  if ((syzstr->betti!=NULL)
1787  && (std_weights))
1788  {
1789    if (minim || (syzstr->resPairs!=NULL))
1790      return ivCopy(syzstr->betti);
1791  }
1792  if ((syzstr->fullres==NULL) && (syzstr->minres==NULL))
1793  {
1794     if (syzstr->hilb_coeffs==NULL)
1795     {
1796        syzstr->fullres = syReorder(syzstr->res,syzstr->length,syzstr);
1797     }
1798     else
1799     {
1800        syzstr->minres = syReorder(syzstr->orderedRes,syzstr->length,syzstr);
1801        syKillEmptyEntres(syzstr->minres,syzstr->length);
1802     }
1803   }
1804  intvec *result=NULL;
1805  if (syzstr->fullres!=NULL)
1806    result = syBetti(syzstr->fullres,syzstr->length,&dummy,weights,minim,row_shift);
1807  else
1808    result = syBetti(syzstr->minres,syzstr->length,&dummy,weights,minim,row_shift);
1809  if ((result!=NULL)
1810  && ((minim) || (syzstr->resPairs!=NULL))
1811  && std_weights)
1812  {
1813    syzstr->betti = ivCopy(result);
1814  }
1815  return result;
1816}
1817
1818/*3
1819* computes the real length of the resolution
1820*/
1821int sySize(syStrategy syzstr)
1822{
1823  resolvente r=syzstr->res;
1824  if (r==NULL)
1825    r = syzstr->fullres;
1826  if (r==NULL)
1827    r = syzstr->minres;
1828  if (r==NULL)
1829  {
1830    WerrorS("No resolution found");
1831    return 0;
1832  }
1833  int i=syzstr->length;
1834  while ((i>0) && (r[i-1]==NULL)) i--;
1835  return i;
1836}
1837
1838/*3
1839* computes the cohomological dimension of res[1]
1840*/
1841int syDim(syStrategy syzstr)
1842{
1843  int i,j=-1,l;
1844  if (syzstr->resPairs!=NULL)
1845  {
1846    SRes rP=syzstr->resPairs;
1847
1848    l = syzstr->length;
1849    while ((l>0) && (rP[l-1]==NULL)) l--;
1850    if (l==0) return -1;
1851    l--;
1852    while (l>=0)
1853    {
1854      i = 0;
1855      while ((i<(*syzstr->Tl)[l]) &&
1856        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)) &&
1857        (rP[l][i].isNotMinimal!=NULL))
1858      {
1859        i++;
1860      }
1861      if ((i<(*syzstr->Tl)[l]) &&
1862        ((rP[l][i].lcm!=NULL) || (rP[l][i].syz!=NULL)) &&
1863        (rP[l][i].isNotMinimal==NULL))
1864        return l;
1865      l--;
1866    }
1867    return l;
1868  }
1869  else
1870    return sySize(syzstr);
1871}
1872
1873/*3
1874* copies the resolution (by increment the reference counter)
1875*/
1876syStrategy syCopy(syStrategy syzstr)
1877{
1878  syStrategy result=syzstr;
1879  (result->references)++;
1880  return result;
1881}
1882
1883/*2
1884* local print procedure used in syPrint
1885*/
1886static void syPrintEmptySpaces(int i)
1887{
1888  if (i!=0)
1889  {
1890    PrintS(" ");
1891    syPrintEmptySpaces(i/10);
1892  }
1893}
1894
1895/*2
1896* local print procedure used in syPrint
1897*/
1898static void syPrintEmptySpaces1(int i)
1899{
1900  if (i!=0)
1901  {
1902    PrintS(" ");
1903    syPrintEmptySpaces1(i-1);
1904  }
1905}
1906
1907/*2
1908* local print procedure used in syPrint
1909*/
1910static int syLengthInt(int i)
1911{
1912  int j=0;
1913
1914  if (i==0) return 1;
1915  while (i!=0)
1916  {
1917    j++;
1918    i = i/10;
1919  }
1920  return j;
1921}
1922
1923/*3
1924* prints the resolution as sequence of free modules
1925*/
1926void syPrint(syStrategy syzstr)
1927{
1928  if ((syzstr->resPairs==NULL) && (syzstr->fullres==NULL)
1929     && (syzstr->minres==NULL))
1930  {
1931    PrintS("No resolution defined\n");
1932    return;
1933  }
1934  int l=0;
1935  if (syzstr->resolution==NULL)
1936  {
1937    int j;
1938    if (syzstr->resPairs!=NULL)
1939    {
1940      syzstr->resolution = new intvec(syzstr->length+1);
1941      SRes rP=syzstr->resPairs;
1942      assume(idRankFreeModule(syzstr->res[1],
1943                                 (syzstr->syRing != NULL ? syzstr->syRing : currRing))==syzstr->res[1]->rank);
1944      (*syzstr->resolution)[0] = syzstr->res[1]->rank;
1945      while ((l<syzstr->length) && (rP[l]!=NULL))
1946      {
1947        j=0;
1948        while ((j<(*syzstr->Tl)[l]) &&
1949          ((rP[l][j].lcm!=NULL) || (rP[l][j].syz!=NULL)))
1950        {
1951          if (rP[l][j].isNotMinimal==NULL)
1952            ((*syzstr->resolution)[l+1])++;
1953          j++;
1954        }
1955        l++;
1956      }
1957    }
1958    else
1959    {
1960      resolvente rr;
1961      syzstr->resolution = new intvec(syzstr->length+2);
1962      if (syzstr->minres!=NULL)
1963        rr = syzstr->minres;
1964      else
1965        rr = syzstr->fullres;
1966      (*syzstr->resolution)[0]
1967        = si_max(1,(int)idRankFreeModule(rr[0],
1968                                 (syzstr->syRing != NULL ? syzstr->syRing : currRing)));
1969      while ((l<syzstr->length) && (rr[l]!=NULL))
1970      {
1971        j = IDELEMS(rr[l]);
1972        while ((j>0) && (rr[l]->m[j-1]==NULL)) j--;
1973        ((*syzstr->resolution)[l+1]) = j;
1974        l++;
1975      }
1976    }
1977  }
1978  const char *sn=currRingHdl->id;
1979  int sl=strlen(sn);
1980  syPrintEmptySpaces1(sl);
1981  l = 0;
1982  loop
1983  {
1984    if ((l>=syzstr->resolution->length()) || ((*syzstr->resolution)[l]==0))
1985      break;
1986    Print("%d",(*syzstr->resolution)[l]);
1987    syPrintEmptySpaces1(sl+5);
1988    l++;
1989  }
1990  PrintLn();
1991  l = 0;
1992  loop
1993  {
1994    if ((l>=syzstr->resolution->length()) || ((*syzstr->resolution)[l]==0))
1995      break;
1996    PrintS(sn);
1997    if (((l+1)>=syzstr->resolution->length()) || ((*syzstr->resolution)[(l+1)]==0))
1998      break;
1999    PrintS(" <-- ");
2000    syPrintEmptySpaces((*syzstr->resolution)[l]);
2001    l++;
2002  }
2003  PrintLn();
2004  PrintLn();
2005  l = 0;
2006  loop
2007  {
2008    if ((l>=syzstr->resolution->length()) || ((*syzstr->resolution)[l]==0))
2009      break;
2010    Print("%d",l);
2011    syPrintEmptySpaces1(sl+5+syLengthInt((*syzstr->resolution)[l])-
2012                         syLengthInt(l));
2013    l++;
2014  }
2015  PrintLn();
2016  if (syzstr->minres==NULL)
2017  {
2018    PrintS("resolution not minimized yet");
2019    PrintLn();
2020  }
2021}
2022
2023/*2
2024* deleting all monomials the component of which correspond
2025* to non-minimal generators
2026*/
2027static poly syStripOut(poly p,intvec * toStrip)
2028{
2029  if (toStrip==NULL) return p;
2030  poly pp=p;
2031
2032  while ((pp!=NULL) && ((*toStrip)[pGetComp(pp)]!=0))
2033    pLmDelete(&pp);
2034  p = pp;
2035  if (pp!=NULL)
2036  {
2037    while (pNext(pp)!=NULL)
2038    {
2039      if ((*toStrip)[pGetComp(pNext(pp))]!=0)
2040        pLmDelete(&pNext(pp));
2041      else
2042        pIter(pp);
2043    }
2044  }
2045  return p;
2046}
2047
2048/*2
2049* copies only those monomials the component of which correspond
2050* to minimal generators
2051*/
2052static poly syStripOutCopy(poly p,intvec * toStrip)
2053{
2054  if (toStrip==NULL) return pCopy(p);
2055  poly result=NULL,pp;
2056
2057  while (p!=NULL)
2058  {
2059    if ((*toStrip)[pGetComp(p)]==0)
2060    {
2061      if (result==NULL)
2062      {
2063        result = pp = pHead(p);
2064      }
2065      else
2066      {
2067        pNext(pp) = pHead(p);
2068        pIter(pp);
2069      }
2070    }
2071    pIter(p);
2072  }
2073  return result;
2074}
2075
2076/*2
2077* minimizes toMin
2078*/
2079static poly syMinimizeP(int toMin,syStrategy syzstr,intvec * ordn,int index,
2080                        intvec * toStrip)
2081{
2082  int ii=0,i,j,tc;
2083  poly p,pp,q=NULL,tq,pisN;
2084  SSet sPairs=syzstr->resPairs[index];
2085  poly tempStripped=NULL;
2086
2087  //pp=pCopy(syzstr->res[index+1]->m[toMin]);
2088  pp = syStripOutCopy(syzstr->res[index+1]->m[toMin],toStrip);
2089  while ((ii<ordn->length()) && ((*ordn)[ii]!=-1) &&
2090             (sPairs[(*ordn)[ii]].syzind!=toMin))
2091  {
2092    ii++;
2093  }
2094  while (ii>=0)
2095  {
2096    i = (*ordn)[ii];
2097    if (sPairs[i].isNotMinimal!=NULL)
2098    {
2099      tempStripped =
2100        syStripOutCopy(syzstr->res[index+1]->m[sPairs[i].syzind],toStrip);
2101      pisN = sPairs[i].isNotMinimal;
2102      tc = pGetComp(pisN);
2103      p = pp;
2104      while (p!=NULL)
2105      {
2106        if (pGetComp(p)==tc)
2107        {
2108          tq = pInit();
2109          for(j=pVariables; j>0; j--)
2110            pSetExp(tq,j, pGetExp(p,j)-pGetExp(pisN,j));
2111          pSetComp(tq, 0);
2112          pSetCoeff0(tq,nDiv(pGetCoeff(p),pGetCoeff(pisN)));
2113          pGetCoeff(tq) = nNeg(pGetCoeff(tq));
2114          pSetm(tq);
2115          q = pAdd(q,pMult_mm(pCopy(tempStripped),tq));
2116          pDelete(&tq);
2117        }
2118        pIter(p);
2119      }
2120      if (q!=NULL)
2121      {
2122        pp = pAdd(pp,q);
2123        q = NULL;
2124      }
2125      pDelete(&tempStripped);
2126    }
2127    ii--;
2128  }
2129  return pp;
2130}
2131
2132/*2
2133* minimizes toMin
2134*/
2135static poly syMinimizeP1(int toMin,syStrategy syzstr,intvec * ordn,int index,
2136                        intvec * toStrip)
2137{
2138  int ii=0,i,j,tc,lp,ltS=-1;
2139  poly p,mp=NULL,pp,q=NULL,tq,pisN;
2140  SSet sPairs=syzstr->resPairs[index];
2141  poly tempStripped=NULL;
2142
2143  pp = syStripOutCopy(syzstr->res[index+1]->m[toMin],toStrip);
2144  kBucketInit(syzstr->bucket,pp,-1);
2145  while ((ii<ordn->length()) && ((*ordn)[ii]!=-1) &&
2146             (sPairs[(*ordn)[ii]].syzind!=toMin))
2147  {
2148    ii++;
2149  }
2150  while (ii>=0)
2151  {
2152    i = (*ordn)[ii];
2153    if (sPairs[i].isNotMinimal!=NULL)
2154    {
2155      tempStripped =
2156        syStripOutCopy(syzstr->res[index+1]->m[sPairs[i].syzind],toStrip);
2157      tc = pGetComp(sPairs[i].isNotMinimal);
2158      //p = pTakeOutComp1(&tempStripped,tc);
2159      int lu;
2160      pTakeOutComp(&tempStripped,tc,&p,&lu);
2161      kBucketTakeOutComp(syzstr->bucket,tc,&mp,&lp);
2162      mp = pDivideM(mp,p);
2163      while (mp!=NULL)
2164      {
2165        p = pNext(mp);
2166        pNext(mp) = NULL;
2167        ltS = -1;
2168        kBucket_Minus_m_Mult_p(syzstr->bucket,mp,tempStripped,&ltS);
2169        mp = p;
2170      }
2171      pDelete(&mp);
2172      pDelete(&tempStripped);
2173    }
2174    ii--;
2175  }
2176  kBucketClear(syzstr->bucket,&pp,&lp);
2177  return pp;
2178}
2179
2180/*2
2181* deletes empty components after minimization
2182*/
2183void syKillEmptyEntres(resolvente res,int length)
2184{
2185  int i,j,jj,k,rj;
2186  intvec * changes;
2187  poly p;
2188  ideal ri;
2189
2190  for (i=0;i<length;i++)
2191  {
2192    ri = res[i];
2193    if (ri!=NULL)
2194    {
2195      rj = IDELEMS(ri);
2196      changes = new intvec(rj+1,1,-1);
2197      while ((rj>0) && (ri->m[rj-1]==NULL)) rj--;
2198      j = k = 0;
2199      while (j+k<rj)
2200      {
2201        if (ri->m[j+k]!=NULL)
2202        {
2203          ri->m[j] = ri->m[j+k];
2204          (*changes)[j+k+1] = j+1;
2205          j++;
2206        }
2207        else
2208        {
2209          k++;
2210        }
2211      }
2212      for (jj=j;jj<rj;jj++)
2213        ri->m[jj] = NULL;
2214      if (res[i+1]!=NULL)
2215      {
2216        ri = res[i+1];
2217        for (j=IDELEMS(ri)-1;j>=0;j--)
2218        {
2219          p = ri->m[j];
2220          while (p!=NULL)
2221          {
2222            pSetComp(p,(*changes)[pGetComp(p)]);
2223            pSetm(p);
2224            pIter(p);
2225          }
2226        }
2227      }
2228      delete changes;
2229    }
2230  }
2231}
2232
2233/*2
2234* determines the components for minimization
2235*/
2236static intvec * syToStrip(syStrategy syzstr, int index)
2237{
2238  intvec * result=NULL;
2239
2240  if ((syzstr->resPairs[index-1]!=NULL) && (!idIs0(syzstr->res[index])))
2241  {
2242    result=new intvec(IDELEMS(syzstr->res[index])+1);
2243    for (int i=(*syzstr->Tl)[index-1]-1;i>=0;i--)
2244    {
2245      if (syzstr->resPairs[index-1][i].isNotMinimal!=NULL)
2246      {
2247        (*result)[syzstr->resPairs[index-1][i].syzind+1] = 1;
2248      }
2249    }
2250  }
2251  return result;
2252}
2253
2254/*2
2255* re-computes the order of pairs during the algorithm
2256* this ensures to procede with a triangular matrix
2257*/
2258static intvec * syOrdPairs(SSet sPairs, int length)
2259{
2260  intvec * result=new intvec(length,1,-1);
2261  int i,j=0,k=-1,l,ii;
2262
2263  loop
2264  {
2265    l = -1;
2266    for(i=0;i<length;i++)
2267    {
2268      if (sPairs[i].syzind>k)
2269      {
2270        if (l==-1)
2271        {
2272          l = sPairs[i].syzind;
2273          ii = i;
2274        }
2275        else
2276        {
2277          if (sPairs[i].syzind<l)
2278          {
2279            l = sPairs[i].syzind;
2280            ii = i;
2281          }
2282        }
2283      }
2284    }
2285    if (l==-1) break;
2286    (*result)[j] = ii;
2287    j++;
2288    k = l;
2289  }
2290  return result;
2291}
2292
2293/*2
2294* minimizes the output of LaScala
2295*/
2296static resolvente syReadOutMinimalRes(syStrategy syzstr,
2297           BOOLEAN computeStd=FALSE)
2298{
2299  intvec * Strip, * ordn;
2300  resolvente tres=(resolvente)omAlloc0((syzstr->length+1)*sizeof(ideal));
2301  ring tmpR = NULL;
2302  ring origR = currRing;
2303
2304//Print("Hier ");
2305  if (computeStd)
2306  {
2307    tres[0] = syzstr->res[1];
2308    syzstr->res[1] = idInit(IDELEMS(tres[0]),tres[0]->rank);
2309    return tres;
2310  }
2311  int i,j,l,index,ii,i1;
2312  poly p;
2313  ideal rs;
2314  SSet sPairs;
2315  int * ord,*b0,*b1;
2316
2317  assume(syzstr->syRing != NULL);
2318  rChangeCurrRing(syzstr->syRing);
2319//Print("laeufts ");
2320  syzstr->bucket = kBucketCreate();
2321  for (index=syzstr->length-1;index>0;index--)
2322  {
2323    if (syzstr->resPairs[index]!=NULL)
2324    {
2325//Print("ideal %d: \n",index);
2326      currcomponents = syzstr->truecomponents[index];
2327      currShiftedComponents = syzstr->ShiftedComponents[index];
2328      rChangeSComps(currcomponents, currShiftedComponents,
2329                    IDELEMS(syzstr->res[index]));
2330      sPairs = syzstr->resPairs[index];
2331      Strip = syToStrip(syzstr,index);
2332      tres[index+1] = idInit(IDELEMS(syzstr->res[index+1]),syzstr->res[index+1]->rank);
2333      i1 = (*syzstr->Tl)[index];
2334//Print("i1= %d\n",i1);
2335      ordn = syOrdPairs(sPairs,i1);
2336      for (i=0;i<i1;i++)
2337      {
2338        if ((sPairs[i].isNotMinimal==NULL) && (sPairs[i].lcm!=NULL))
2339        {
2340          l = sPairs[i].syzind;
2341//Print("Minimiere Poly %d: ",l);pWrite(syzstr->res[index+1]->m[l]);
2342          tres[index+1]->m[l] =
2343            syMinimizeP1(l,syzstr,ordn,index,Strip);
2344        }
2345      }
2346      delete Strip;
2347      Strip = NULL;
2348    }
2349  }
2350  currcomponents = syzstr->truecomponents[0];
2351  currShiftedComponents = syzstr->ShiftedComponents[0];
2352  rChangeSComps( currcomponents, currShiftedComponents,
2353                 IDELEMS(syzstr->res[0]));
2354  tres[1] = idInit(IDELEMS(syzstr->res[1]),syzstr->res[1]->rank);
2355  sPairs = syzstr->resPairs[0];
2356  for (i=(*syzstr->Tl)[0]-1;i>=0;i--)
2357  {
2358    if (sPairs[i].syzind>=0)
2359    {
2360      tres[1]->m[sPairs[i].syzind] = pCopy(syzstr->res[1]->m[sPairs[i].syzind]);
2361    }
2362  }
2363/*--- changes to the original ring------------------*/
2364  kBucketDestroy(&syzstr->bucket);
2365  if (syzstr->syRing != NULL)
2366  {
2367    rChangeCurrRing(origR);
2368    // Thomas: now make sure that all data which you need is pFetchCopied
2369    // maybe incoporate it into syReorder ??
2370  }
2371  tres = syReorder(tres,syzstr->length,syzstr,FALSE,syzstr->res);
2372  syKillEmptyEntres(tres,syzstr->length);
2373  idSkipZeroes(tres[0]);
2374  return tres;
2375}
2376
2377/*3
2378* minimizes any kind of resolution
2379*/
2380syStrategy syMinimize(syStrategy syzstr)
2381{
2382  if (syzstr->minres==NULL)
2383  {
2384    if (syzstr->resPairs!=NULL)
2385    {
2386      if (syzstr->hilb_coeffs==NULL)
2387      {
2388        syzstr->minres = syReadOutMinimalRes(syzstr);
2389      }
2390      else
2391      {
2392        syzstr->minres = syReorder(syzstr->orderedRes,syzstr->length,syzstr);
2393      }
2394    }
2395    else if (syzstr->fullres!=NULL)
2396    {
2397      syMinimizeResolvente(syzstr->fullres,syzstr->length,1);
2398      syzstr->minres = syzstr->fullres;
2399      syzstr->fullres = NULL;
2400    }
2401  }
2402  (syzstr->references)++;
2403  return syzstr;
2404}
2405
2406/*2
2407* implementation of LaScala's algorithm
2408* assumes that the given module is homogeneous
2409* works with slanted degree, uses syChosePairs
2410*/
2411syStrategy syLaScala3(ideal arg,int * length)
2412{
2413  BOOLEAN noPair=FALSE;
2414  int i,j,actdeg=32000,index=0,reg=-1;
2415  int startdeg,howmuch;
2416  poly p;
2417  ideal temp;
2418  SSet nextPairs;
2419  syStrategy syzstr=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2420  ring origR = currRing;
2421
2422  if ((idIs0(arg)) ||
2423      ((idRankFreeModule(arg)>0) && (!idHomModule(arg,NULL,&(syzstr->cw)))))
2424  {
2425    syzstr->minres = (resolvente)omAlloc0Bin(char_ptr_bin);
2426    syzstr->length = 1;
2427    syzstr->minres[0] = idInit(1,arg->rank);
2428    return syzstr;
2429  }
2430
2431  //crit = 0;
2432  //euler = -1;
2433  redpol = pInit();
2434  syzstr->length = *length = pVariables+2;
2435
2436  // Creare dp,S ring and change to it
2437  syzstr->syRing = rCurrRingAssure_dp_S();
2438  assume(syzstr->syRing != origR);
2439
2440  // set initial ShiftedComps
2441  currcomponents = (int*)omAlloc0((arg->rank+1)*sizeof(int));
2442  currShiftedComponents = (long*)omAlloc0((arg->rank+1)*sizeof(long));
2443  for (i=0;i<=arg->rank;i++)
2444  {
2445    currShiftedComponents[i] = (i)*SYZ_SHIFT_BASE;
2446    currcomponents[i] = i;
2447  }
2448  rChangeSComps(currcomponents, currShiftedComponents, arg->rank);
2449/*--- initializes the data structures---------------*/
2450  syzstr->Tl = new intvec(*length);
2451  temp = idInit(IDELEMS(arg),arg->rank);
2452  for (i=0;i<IDELEMS(arg);i++)
2453  {
2454    temp->m[i] = prCopyR( arg->m[i], origR);
2455    if (temp->m[i]!=NULL)
2456    {
2457      j = pTotaldegree(temp->m[i]);
2458      if (j<actdeg) actdeg = j;
2459    }
2460  }
2461  idTest(temp);
2462  idSkipZeroes(temp);
2463  idTest(temp);
2464  syzstr->resPairs = syInitRes(temp,length,syzstr->Tl,syzstr->cw);
2465  omFreeSize((ADDRESS)currcomponents,(arg->rank+1)*sizeof(int));
2466  omFreeSize((ADDRESS)currShiftedComponents,(arg->rank+1)*sizeof(long));
2467  syzstr->res = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
2468  syzstr->orderedRes = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
2469  syzstr->elemLength = (int**)omAlloc0((*length+1)*sizeof(int*));
2470  syzstr->truecomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
2471  syzstr->ShiftedComponents = (long**)omAlloc0((*length+1)*sizeof(long*));
2472  syzstr->backcomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
2473  syzstr->Howmuch = (int**)omAlloc0((*length+1)*sizeof(int*));
2474  syzstr->Firstelem = (int**)omAlloc0((*length+1)*sizeof(int*));
2475  syzstr->sev = (unsigned long **) omAlloc0((*length+1)*sizeof(unsigned long *));
2476  syzstr->bucket = kBucketCreate();
2477  int len0=idRankFreeModule(temp)+1;
2478
2479  startdeg = actdeg;
2480  nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2481  //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2482/*--- computes the resolution ----------------------*/
2483  while (nextPairs!=NULL)
2484  {
2485    if (TEST_OPT_PROT) Print("%d",actdeg);
2486    if (TEST_OPT_PROT) Print("(m%d)",index);
2487    if (index==0)
2488      i = syInitSyzMod(syzstr,index,len0);
2489    else
2490      i = syInitSyzMod(syzstr,index);
2491    currcomponents = syzstr->truecomponents[si_max(index-1,0)];
2492    currShiftedComponents = syzstr->ShiftedComponents[si_max(index-1,0)];
2493    rChangeSComps(currcomponents, currShiftedComponents,
2494                  IDELEMS(syzstr->res[si_max(index-1,0)]));
2495    j = syInitSyzMod(syzstr,index+1);
2496    if (index>0)
2497    {
2498      syRedNextPairs(nextPairs,syzstr,howmuch,index);
2499      syCompactifyPairSet(syzstr->resPairs[index],(*syzstr->Tl)[index],0);
2500    }
2501    else
2502      syRedGenerOfCurrDeg(syzstr,actdeg,index+1);
2503/*--- creates new pairs -----------------------------*/
2504    syCreateNewPairs(syzstr,index,i);
2505    if (index<(*length)-1)
2506    {
2507      syCreateNewPairs(syzstr,index+1,j);
2508    }
2509    index++;
2510    nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2511    //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2512  }
2513  if (temp!=NULL) idDelete(&temp);
2514  kBucketDestroy(&(syzstr->bucket));
2515  if (origR != syzstr->syRing)
2516    rChangeCurrRing(origR);
2517  pLmDelete(&redpol);
2518  if (TEST_OPT_PROT) PrintLn();
2519  return syzstr;
2520}
2521
2522
2523
2524/*2
2525* more general implementation of LaScala's algorithm
2526* assumes that the given module is (quasi-)homogeneous
2527* works with slanted degree, uses syChosePairs
2528*/
2529syStrategy syLaScala(ideal arg, int& maxlength, intvec* weights)
2530{
2531  BOOLEAN noPair=FALSE;
2532  int i,j,actdeg=32000,index=0,reg=-1;
2533  int startdeg,howmuch;
2534  poly p;
2535  ideal temp;
2536  SSet nextPairs;
2537  syStrategy syzstr=(syStrategy)omAlloc0(sizeof(ssyStrategy));
2538  ring origR = currRing;
2539
2540  if(weights!= NULL) 
2541    syzstr->cw = new intvec(weights);
2542  else
2543    syzstr->cw = NULL;
2544
2545  if ((idIs0(arg)) ||
2546      ((idRankFreeModule(arg)>0) && (!idTestHomModule(arg, NULL, syzstr->cw))))
2547  {
2548    syzstr->minres = (resolvente)omAlloc0Bin(char_ptr_bin);
2549    syzstr->length = 1;
2550    syzstr->minres[0] = idInit(1,arg->rank);
2551    return syzstr;
2552  }
2553
2554
2555  //crit = 0;
2556  //euler = -1;
2557  redpol = pInit();
2558 
2559  if( maxlength > 0 )
2560    syzstr->length = maxlength; //  = pVariables+2;
2561  else
2562    syzstr->length = maxlength = pVariables+2; 
2563
2564  // Creare dp,S ring and change to it
2565  syzstr->syRing = rCurrRingAssure_dp_S();
2566  assume(syzstr->syRing != origR);
2567
2568  // set initial ShiftedComps
2569  currcomponents = (int*)omAlloc0((arg->rank+1)*sizeof(int));
2570  currShiftedComponents = (long*)omAlloc0((arg->rank+1)*sizeof(long));
2571  for (i=0;i<=arg->rank;i++)
2572  {
2573    currShiftedComponents[i] = (i)*SYZ_SHIFT_BASE;
2574    currcomponents[i] = i;
2575  }
2576  rChangeSComps(currcomponents, currShiftedComponents, arg->rank);
2577/*--- initializes the data structures---------------*/
2578  syzstr->Tl = new intvec(maxlength);
2579  temp = idInit(IDELEMS(arg),arg->rank);
2580  for (i=0;i<IDELEMS(arg);i++)
2581  {
2582    temp->m[i] = prCopyR( arg->m[i], origR);
2583    if (temp->m[i]!=NULL)
2584    {
2585      j = pTotaldegree(temp->m[i]);
2586      if (j<actdeg) actdeg = j;
2587    }
2588  }
2589  idTest(temp);
2590  idSkipZeroes(temp);
2591  idTest(temp);
2592  syzstr->resPairs = syInitRes(temp,&maxlength,syzstr->Tl,syzstr->cw);
2593  omFreeSize((ADDRESS)currcomponents,(arg->rank+1)*sizeof(int));
2594  omFreeSize((ADDRESS)currShiftedComponents,(arg->rank+1)*sizeof(long));
2595  syzstr->res = (resolvente)omAlloc0((maxlength+1)*sizeof(ideal));
2596  syzstr->orderedRes = (resolvente)omAlloc0((maxlength+1)*sizeof(ideal));
2597  syzstr->elemLength = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2598  syzstr->truecomponents = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2599  syzstr->ShiftedComponents = (long**)omAlloc0((maxlength+1)*sizeof(long*));
2600  syzstr->backcomponents = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2601  syzstr->Howmuch = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2602  syzstr->Firstelem = (int**)omAlloc0((maxlength+1)*sizeof(int*));
2603  syzstr->sev = (unsigned long **) omAlloc0((maxlength+1)*sizeof(unsigned long *));
2604  syzstr->bucket = kBucketCreate();
2605  int len0=idRankFreeModule(temp)+1;
2606
2607  startdeg = actdeg;
2608  nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2609  //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2610/*--- computes the resolution ----------------------*/
2611  while (nextPairs!=NULL)
2612  {
2613    if (TEST_OPT_PROT) Print("%d",actdeg);
2614    if (TEST_OPT_PROT) Print("(m%d)",index);
2615    if (index==0)
2616      i = syInitSyzMod(syzstr,index,len0);
2617    else
2618      i = syInitSyzMod(syzstr,index);
2619    currcomponents = syzstr->truecomponents[si_max(index-1,0)];
2620    currShiftedComponents = syzstr->ShiftedComponents[si_max(index-1,0)];
2621    rChangeSComps(currcomponents, currShiftedComponents,
2622                  IDELEMS(syzstr->res[si_max(index-1,0)]));
2623    j = syInitSyzMod(syzstr,index+1);
2624    if (index>0)
2625    {
2626      syRedNextPairs(nextPairs,syzstr,howmuch,index);
2627      syCompactifyPairSet(syzstr->resPairs[index],(*syzstr->Tl)[index],0);
2628    }
2629    else
2630      syRedGenerOfCurrDeg(syzstr,actdeg,index+1);
2631/*--- creates new pairs -----------------------------*/
2632    syCreateNewPairs(syzstr,index,i);
2633    if (index<(maxlength-1))
2634    {
2635      syCreateNewPairs(syzstr,index+1,j);
2636    }
2637    index++;
2638    nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
2639    //if (TEST_OPT_PROT) Print("(%d,%d)",howmuch,index);
2640  }
2641  if (temp!=NULL) idDelete(&temp);
2642  kBucketDestroy(&(syzstr->bucket));
2643  if (origR != syzstr->syRing)
2644    rChangeCurrRing(origR);
2645  pLmDelete(&redpol);
2646  if (TEST_OPT_PROT) PrintLn();
2647  return syzstr;
2648}
2649
Note: See TracBrowser for help on using the repository browser.