source: git/kernel/syz1.cc @ c68d10

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