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

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