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

spielwiese
Last change on this file since 1f637e was 1f637e, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
FIX: replaced pVariables by currRing->N... ADD: added 'void rChangeCurrRing(ring r);' to polys/polys.h CHG: started changing gr_kstd2.cc and syz*.cc
  • Property mode set to 100644
File size: 28.6 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6*  ABSTRACT -  dimension, multiplicity, HC, kbase
7*/
8
9#include <kernel/mod2.h>
10
11
12#include <kernel/structs.h>
13#include <kernel/febase.h>
14#include <omalloc/omalloc.h>
15#include <kernel/ideals.h>
16#include <polys/polys.h>
17#include <misc/intvec.h>
18#include <coeffs/numbers.h>
19#include <kernel/hutil.h>
20#include <kernel/stairc.h>
21
22int  hCo, hMu, hMu2;
23omBin indlist_bin = omGetSpecBin(sizeof(indlist));
24
25/*0 implementation*/
26
27// dimension
28
29void hDimSolve(scmon pure, int Npure, scfmon rad, int Nrad,
30 varset var, int Nvar)
31{
32  int  dn, iv, rad0, b, c, x;
33  scmon pn;
34  scfmon rn;
35  if (Nrad < 2)
36  {
37    dn = Npure + Nrad;
38    if (dn < hCo)
39      hCo = dn;
40    return;
41  }
42  if (Npure+1 >= hCo)
43    return;
44  iv = Nvar;
45  while(pure[var[iv]]) iv--;
46  hStepR(rad, Nrad, var, iv, &rad0);
47  if (rad0!=0)
48  {
49    iv--;
50    if (rad0 < Nrad)
51    {
52      pn = hGetpure(pure);
53      rn = hGetmem(Nrad, rad, radmem[iv]);
54      hDimSolve(pn, Npure + 1, rn, rad0, var, iv);
55      b = rad0;
56      c = Nrad;
57      hElimR(rn, &rad0, b, c, var, iv);
58      hPure(rn, b, &c, var, iv, pn, &x);
59      hLex2R(rn, rad0, b, c, var, iv, hwork);
60      rad0 += (c - b);
61      hDimSolve(pn, Npure + x, rn, rad0, var, iv);
62    }
63    else
64    {
65      hDimSolve(pure, Npure, rad, Nrad, var, iv);
66    }
67  }
68  else
69    hCo = Npure + 1;
70}
71
72int  scDimInt(ideal S, ideal Q)
73{
74  int  mc;
75  hexist = hInit(S, Q, &hNexist);
76  if (!hNexist)
77    return (currRing->N);
78  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
79  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
80  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
81  mc = hisModule;
82  if (!mc)
83  {
84    hrad = hexist;
85    hNrad = hNexist;
86  }
87  else
88    hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
89  radmem = hCreate((currRing->N) - 1);
90  hCo = (currRing->N) + 1;
91  loop
92  {
93    if (mc)
94      hComp(hexist, hNexist, mc, hrad, &hNrad);
95    if (hNrad)
96    {
97      hNvar = (currRing->N);
98      hRadical(hrad, &hNrad, hNvar);
99      hSupp(hrad, hNrad, hvar, &hNvar);
100      if (hNvar)
101      {
102        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
103        hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
104        hLexR(hrad, hNrad, hvar, hNvar);
105        hDimSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
106      }
107    }
108    else
109    {
110      hCo = 0;
111      break;
112    }
113    mc--;
114    if (mc <= 0)
115      break;
116  }
117  hKill(radmem, (currRing->N) - 1);
118  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
119  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
120  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
121  hDelete(hexist, hNexist);
122  if (hisModule)
123    omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
124  return (currRing->N) - hCo;
125}
126
127// independent set
128static scmon hInd;
129
130static void hIndSolve(scmon pure, int Npure, scfmon rad, int Nrad,
131 varset var, int Nvar)
132{
133  int  dn, iv, rad0, b, c, x;
134  scmon pn;
135  scfmon rn;
136  if (Nrad < 2)
137  {
138    dn = Npure + Nrad;
139    if (dn < hCo)
140    {
141      hCo = dn;
142      for (iv=(currRing->N); iv; iv--)
143      {
144        if (pure[iv])
145          hInd[iv] = 0;
146        else
147          hInd[iv] = 1;
148      }
149      if (Nrad)
150      {
151        pn = *rad;
152        iv = Nvar;
153        loop
154        {
155          x = var[iv];
156          if (pn[x])
157          {
158            hInd[x] = 0;
159            break;
160          }
161          iv--;
162        }
163      }
164    }
165    return;
166  }
167  if (Npure+1 >= hCo)
168    return;
169  iv = Nvar;
170  while(pure[var[iv]]) iv--;
171  hStepR(rad, Nrad, var, iv, &rad0);
172  if (rad0)
173  {
174    iv--;
175    if (rad0 < Nrad)
176    {
177      pn = hGetpure(pure);
178      rn = hGetmem(Nrad, rad, radmem[iv]);
179      pn[var[iv + 1]] = 1;
180      hIndSolve(pn, Npure + 1, rn, rad0, var, iv);
181      pn[var[iv + 1]] = 0;
182      b = rad0;
183      c = Nrad;
184      hElimR(rn, &rad0, b, c, var, iv);
185      hPure(rn, b, &c, var, iv, pn, &x);
186      hLex2R(rn, rad0, b, c, var, iv, hwork);
187      rad0 += (c - b);
188      hIndSolve(pn, Npure + x, rn, rad0, var, iv);
189    }
190    else
191    {
192      hIndSolve(pure, Npure, rad, Nrad, var, iv);
193    }
194  }
195  else
196  {
197    hCo = Npure + 1;
198    for (x=(currRing->N); x; x--)
199    {
200      if (pure[x])
201        hInd[x] = 0;
202      else
203        hInd[x] = 1;
204    }
205    hInd[var[iv]] = 0;
206  }
207}
208
209intvec * scIndIntvec(ideal S, ideal Q)
210{
211  intvec *Set=new intvec((currRing->N));
212  int  mc,i;
213  hexist = hInit(S, Q, &hNexist);
214  if (hNexist==0)
215  {
216    for(i=0; i<(currRing->N); i++)
217      (*Set)[i]=1;
218    return Set;
219  }
220  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
221  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
222  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
223  hInd = (scmon)omAlloc((1 + (currRing->N)) * sizeof(int));
224  mc = hisModule;
225  if (mc==0)
226  {
227    hrad = hexist;
228    hNrad = hNexist;
229  }
230  else
231    hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
232  radmem = hCreate((currRing->N) - 1);
233  hCo = (currRing->N) + 1;
234  loop
235  {
236    if (mc!=0)
237      hComp(hexist, hNexist, mc, hrad, &hNrad);
238    if (hNrad!=0)
239    {
240      hNvar = (currRing->N);
241      hRadical(hrad, &hNrad, hNvar);
242      hSupp(hrad, hNrad, hvar, &hNvar);
243      if (hNvar!=0)
244      {
245        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
246        hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
247        hLexR(hrad, hNrad, hvar, hNvar);
248        hIndSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
249      }
250    }
251    else
252    {
253      hCo = 0;
254      break;
255    }
256    mc--;
257    if (mc <= 0)
258      break;
259  }
260  for(i=0; i<(currRing->N); i++)
261    (*Set)[i] = hInd[i+1];
262  hKill(radmem, (currRing->N) - 1);
263  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
264  omFreeSize((ADDRESS)hInd, (1 + (currRing->N)) * sizeof(int));
265  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
266  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
267  hDelete(hexist, hNexist);
268  if (hisModule)
269    omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
270  return Set;
271}
272
273indset ISet, JSet;
274
275static BOOLEAN hNotZero(scfmon rad, int Nrad, varset var, int Nvar)
276{
277  int  k1, i;
278  k1 = var[Nvar];
279  i = 0;
280  loop
281  {
282    if (rad[i][k1]==0)
283      return FALSE;
284    i++;
285    if (i == Nrad)
286      return TRUE;
287  }
288}
289
290static void hIndep(scmon pure)
291{
292  int iv;
293  intvec *Set;
294
295  Set = ISet->set = new intvec((currRing->N));
296  for (iv=(currRing->N); iv!=0 ; iv--)
297  {
298    if (pure[iv])
299      (*Set)[iv-1] = 0;
300    else
301      (*Set)[iv-1] = 1;
302  }
303  ISet = ISet->nx = (indset)omAlloc0Bin(indlist_bin);
304  hMu++;
305}
306
307void hIndMult(scmon pure, int Npure, scfmon rad, int Nrad,
308 varset var, int Nvar)
309{
310  int  dn, iv, rad0, b, c, x;
311  scmon pn;
312  scfmon rn;
313  if (Nrad < 2)
314  {
315    dn = Npure + Nrad;
316    if (dn == hCo)
317    {
318      if (Nrad==0)
319        hIndep(pure);
320      else
321      {
322        pn = *rad;
323        for (iv = Nvar; iv!=0; iv--)
324        {
325          x = var[iv];
326          if (pn[x])
327          {
328            pure[x] = 1;
329            hIndep(pure);
330            pure[x] = 0;
331          }
332        }
333      }
334    }
335    return;
336  }
337  iv = Nvar;
338  dn = Npure+1;
339  if (dn >= hCo)
340  {
341    if (dn > hCo)
342      return;
343    loop
344    {
345      if(!pure[var[iv]])
346      {
347        if(hNotZero(rad, Nrad, var, iv))
348        {
349          pure[var[iv]] = 1;
350          hIndep(pure);
351          pure[var[iv]] = 0;
352        }
353      }
354      iv--;
355      if (!iv)
356        return;
357    }
358  }
359  while(pure[var[iv]]) iv--;
360  hStepR(rad, Nrad, var, iv, &rad0);
361  iv--;
362  if (rad0 < Nrad)
363  {
364    pn = hGetpure(pure);
365    rn = hGetmem(Nrad, rad, radmem[iv]);
366    pn[var[iv + 1]] = 1;
367    hIndMult(pn, Npure + 1, rn, rad0, var, iv);
368    pn[var[iv + 1]] = 0;
369    b = rad0;
370    c = Nrad;
371    hElimR(rn, &rad0, b, c, var, iv);
372    hPure(rn, b, &c, var, iv, pn, &x);
373    hLex2R(rn, rad0, b, c, var, iv, hwork);
374    rad0 += (c - b);
375    hIndMult(pn, Npure + x, rn, rad0, var, iv);
376  }
377  else
378  {
379    hIndMult(pure, Npure, rad, Nrad, var, iv);
380  }
381}
382
383/*3
384* consider indset x := !pure
385* (for all i) (if(sm(i) > x) return FALSE)
386* else return TRUE
387*/
388static BOOLEAN hCheck1(indset sm, scmon pure)
389{
390  int iv;
391  intvec *Set;
392  while (sm->nx != NULL)
393  {
394    Set = sm->set;
395    iv=(currRing->N);
396    loop
397    {
398      if (((*Set)[iv-1] == 0) && (pure[iv] == 0))
399        break;
400      iv--;
401      if (iv == 0)
402        return FALSE;
403    }
404    sm = sm->nx;
405  }
406  return TRUE;
407}
408
409/*3
410* consider indset x := !pure
411* (for all i) if(x > sm(i)) delete sm(i))
412* return (place for x)
413*/
414static indset hCheck2(indset sm, scmon pure)
415{
416  int iv;
417  intvec *Set;
418  indset be, a1 = NULL;
419  while (sm->nx != NULL)
420  {
421    Set = sm->set;
422    iv=(currRing->N);
423    loop
424    {
425      if ((pure[iv] == 1) && ((*Set)[iv-1] == 1))
426        break;
427      iv--;
428      if (iv == 0)
429      {
430        if (a1 == NULL)
431        {
432          a1 = sm;
433        }
434        else
435        {
436          hMu2--;
437          be->nx = sm->nx;
438          delete Set;
439          omFreeBin((ADDRESS)sm, indlist_bin);
440          sm = be;
441        }
442        break;
443      }
444    }
445    be = sm;
446    sm = sm->nx;
447  }
448  if (a1 != NULL)
449  {
450    return a1;
451  }
452  else
453  {
454    hMu2++;
455    sm->set = new intvec((currRing->N));
456    sm->nx = (indset)omAlloc0Bin(indlist_bin);
457    return sm;
458  }
459}
460
461/*2
462*  definition x >= y
463*      x(i) == 0 => y(i) == 0
464*      > ex. j with x(j) == 1 and y(j) == 0
465*/
466static void hCheckIndep(scmon pure)
467{
468  intvec *Set;
469  indset res;
470  int iv;
471  if (hCheck1(ISet, pure))
472  {
473    if (hCheck1(JSet, pure))
474    {
475      res = hCheck2(JSet,pure);
476      if (res == NULL)
477        return;
478      Set = res->set;
479      for (iv=(currRing->N); iv; iv--)
480      {
481        if (pure[iv])
482          (*Set)[iv-1] = 0;
483        else
484          (*Set)[iv-1] = 1;
485      }
486    }
487  }
488}
489
490void hIndAllMult(scmon pure, int Npure, scfmon rad, int Nrad,
491 varset var, int Nvar)
492{
493  int  dn, iv, rad0, b, c, x;
494  scmon pn;
495  scfmon rn;
496  if (Nrad < 2)
497  {
498    dn = Npure + Nrad;
499    if (dn > hCo)
500    {
501      if (!Nrad)
502        hCheckIndep(pure);
503      else
504      {
505        pn = *rad;
506        for (iv = Nvar; iv; iv--)
507        {
508          x = var[iv];
509          if (pn[x])
510          {
511            pure[x] = 1;
512            hCheckIndep(pure);
513            pure[x] = 0;
514          }
515        }
516      }
517    }
518    return;
519  }
520  iv = Nvar;
521  while(pure[var[iv]]) iv--;
522  hStepR(rad, Nrad, var, iv, &rad0);
523  iv--;
524  if (rad0 < Nrad)
525  {
526    pn = hGetpure(pure);
527    rn = hGetmem(Nrad, rad, radmem[iv]);
528    pn[var[iv + 1]] = 1;
529    hIndAllMult(pn, Npure + 1, rn, rad0, var, iv);
530    pn[var[iv + 1]] = 0;
531    b = rad0;
532    c = Nrad;
533    hElimR(rn, &rad0, b, c, var, iv);
534    hPure(rn, b, &c, var, iv, pn, &x);
535    hLex2R(rn, rad0, b, c, var, iv, hwork);
536    rad0 += (c - b);
537    hIndAllMult(pn, Npure + x, rn, rad0, var, iv);
538  }
539  else
540  {
541    hIndAllMult(pure, Npure, rad, Nrad, var, iv);
542  }
543}
544
545// multiplicity
546
547static int hZeroMult(scmon pure, scfmon stc, int Nstc, varset var, int Nvar)
548{
549  int  iv = Nvar -1, sum, a, a0, a1, b, i;
550  int  x, x0;
551  scmon pn;
552  scfmon sn;
553  if (!iv)
554    return pure[var[1]];
555  else if (!Nstc)
556  {
557    sum = 1;
558    for (i = Nvar; i; i--)
559      sum *= pure[var[i]];
560    return sum;
561  }
562  x = a = 0;
563  pn = hGetpure(pure);
564  sn = hGetmem(Nstc, stc, stcmem[iv]);
565  hStepS(sn, Nstc, var, Nvar, &a, &x);
566  if (a == Nstc)
567    return pure[var[Nvar]] * hZeroMult(pn, sn, a, var, iv);
568  else
569    sum = x * hZeroMult(pn, sn, a, var, iv);
570  b = a;
571  loop
572  {
573    a0 = a;
574    x0 = x;
575    hStepS(sn, Nstc, var, Nvar, &a, &x);
576    hElimS(sn, &b, a0, a, var, iv);
577    a1 = a;
578    hPure(sn, a0, &a1, var, iv, pn, &i);
579    hLex2S(sn, b, a0, a1, var, iv, hwork);
580    b += (a1 - a0);
581    if (a < Nstc)
582    {
583      sum += (x - x0) * hZeroMult(pn, sn, b, var, iv);
584    }
585    else
586    {
587      sum += (pure[var[Nvar]] - x0) * hZeroMult(pn, sn, b, var, iv);
588      return sum;
589    }
590  }
591}
592
593static void hProject(scmon pure, varset sel)
594{
595  int  i, i0, k;
596  i0 = 0;
597  for (i = 1; i <= (currRing->N); i++)
598  {
599    if (pure[i])
600    {
601      i0++;
602      sel[i0] = i;
603    }
604  }
605  i = hNstc;
606  memcpy(hwork, hstc, i * sizeof(scmon));
607  hStaircase(hwork, &i, sel, i0);
608  if ((i0 > 2) && (i > 10))
609    hOrdSupp(hwork, i, sel, i0);
610  memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
611  hPure(hwork, 0, &i, sel, i0, hpur0, &k);
612  hLexS(hwork, i, sel, i0);
613  hMu += hZeroMult(hpur0, hwork, i, sel, i0);
614}
615
616static void hDimMult(scmon pure, int Npure, scfmon rad, int Nrad,
617 varset var, int Nvar)
618{
619  int  dn, iv, rad0, b, c, x;
620  scmon pn;
621  scfmon rn;
622  if (Nrad < 2)
623  {
624    dn = Npure + Nrad;
625    if (dn == hCo)
626    {
627      if (!Nrad)
628        hProject(pure, hsel);
629      else
630      {
631        pn = *rad;
632        for (iv = Nvar; iv; iv--)
633        {
634          x = var[iv];
635          if (pn[x])
636          {
637            pure[x] = 1;
638            hProject(pure, hsel);
639            pure[x] = 0;
640          }
641        }
642      }
643    }
644    return;
645  }
646  iv = Nvar;
647  dn = Npure+1;
648  if (dn >= hCo)
649  {
650    if (dn > hCo)
651      return;
652    loop
653    {
654      if(!pure[var[iv]])
655      {
656        if(hNotZero(rad, Nrad, var, iv))
657        {
658          pure[var[iv]] = 1;
659          hProject(pure, hsel);
660          pure[var[iv]] = 0;
661        }
662      }
663      iv--;
664      if (!iv)
665        return;
666    }
667  }
668  while(pure[var[iv]]) iv--;
669  hStepR(rad, Nrad, var, iv, &rad0);
670  iv--;
671  if (rad0 < Nrad)
672  {
673    pn = hGetpure(pure);
674    rn = hGetmem(Nrad, rad, radmem[iv]);
675    pn[var[iv + 1]] = 1;
676    hDimMult(pn, Npure + 1, rn, rad0, var, iv);
677    pn[var[iv + 1]] = 0;
678    b = rad0;
679    c = Nrad;
680    hElimR(rn, &rad0, b, c, var, iv);
681    hPure(rn, b, &c, var, iv, pn, &x);
682    hLex2R(rn, rad0, b, c, var, iv, hwork);
683    rad0 += (c - b);
684    hDimMult(pn, Npure + x, rn, rad0, var, iv);
685  }
686  else
687  {
688    hDimMult(pure, Npure, rad, Nrad, var, iv);
689  }
690}
691
692static void hDegree(ideal S, ideal Q)
693{
694  int  di;
695  int  mc;
696  hexist = hInit(S, Q, &hNexist);
697  if (!hNexist)
698  {
699    hCo = 0;
700    hMu = 1;
701    return;
702  }
703  //hWeight();
704  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
705  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
706  hsel = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
707  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
708  hpur0 = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
709  mc = hisModule;
710  hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
711  if (!mc)
712  {
713    memcpy(hrad, hexist, hNexist * sizeof(scmon));
714    hstc = hexist;
715    hNrad = hNstc = hNexist;
716  }
717  else
718    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
719  radmem = hCreate((currRing->N) - 1);
720  stcmem = hCreate((currRing->N) - 1);
721  hCo = (currRing->N) + 1;
722  di = hCo + 1;
723  loop
724  {
725    if (mc)
726    {
727      hComp(hexist, hNexist, mc, hrad, &hNrad);
728      hNstc = hNrad;
729      memcpy(hstc, hrad, hNrad * sizeof(scmon));
730    }
731    if (hNrad)
732    {
733      hNvar = (currRing->N);
734      hRadical(hrad, &hNrad, hNvar);
735      hSupp(hrad, hNrad, hvar, &hNvar);
736      if (hNvar)
737      {
738        hCo = hNvar;
739        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
740        hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
741        hLexR(hrad, hNrad, hvar, hNvar);
742        hDimSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
743      }
744    }
745    else
746    {
747      hNvar = 1;
748      hCo = 0;
749    }
750    if (hCo < di)
751    {
752      di = hCo;
753      hMu = 0;
754    }
755    if (hNvar && (hCo == di))
756    {
757      if (di && (di < (currRing->N)))
758        hDimMult(hpure, hNpure, hrad, hNrad, hvar, hNvar);
759      else if (!di)
760        hMu++;
761      else
762      {
763        hStaircase(hstc, &hNstc, hvar, hNvar);
764        if ((hNvar > 2) && (hNstc > 10))
765          hOrdSupp(hstc, hNstc, hvar, hNvar);
766        memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
767        hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
768        hLexS(hstc, hNstc, hvar, hNvar);
769        hMu += hZeroMult(hpur0, hstc, hNstc, hvar, hNvar);
770      }
771    }
772    mc--;
773    if (mc <= 0)
774      break;
775  }
776  hCo = di;
777  hKill(stcmem, (currRing->N) - 1);
778  hKill(radmem, (currRing->N) - 1);
779  omFreeSize((ADDRESS)hpur0, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
780  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
781  omFreeSize((ADDRESS)hsel, ((currRing->N) + 1) * sizeof(int));
782  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
783  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
784  omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
785  hDelete(hexist, hNexist);
786  if (hisModule)
787    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
788}
789
790int  scMultInt(ideal S, ideal Q)
791{
792  hDegree(S, Q);
793  return hMu;
794}
795
796void scPrintDegree(int co, int mu)
797{
798  int di = (currRing->N)-co;
799  if (pOrdSgn == 1)
800  {
801    if (di>0)
802      Print("// dimension (proj.)  = %d\n// degree (proj.)   = %d\n", di-1, mu);
803    else
804      Print("// dimension (affine) = 0\n// degree (affine)  = %d\n",       mu);
805  }
806  else
807    Print("// dimension (local)   = %d\n// multiplicity = %d\n", di, mu);
808}
809
810void scDegree(ideal S, intvec *modulweight, ideal Q)
811{
812  int co, mu, l;
813  intvec *hseries2;
814  intvec *hseries1 = hFirstSeries(S, modulweight, Q);
815  l = hseries1->length()-1;
816  if (l > 1)
817    hseries2 = hSecondSeries(hseries1);
818  else
819    hseries2 = hseries1;
820  hDegreeSeries(hseries1, hseries2, &co, &mu);
821  if ((l == 1) &&(mu == 0))
822    scPrintDegree((currRing->N)+1, 0);
823  else
824    scPrintDegree(co, mu);
825  if (l>1)
826    delete hseries1;
827  delete hseries2;
828}
829
830static void hDegree0(ideal S, ideal Q)
831{
832  int  mc;
833  hexist = hInit(S, Q, &hNexist);
834  if (!hNexist)
835  {
836    hMu = -1;
837    return;
838  }
839  else
840    hMu = 0;
841  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
842  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
843  hpur0 = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
844  mc = hisModule;
845  if (!mc)
846  {
847    hstc = hexist;
848    hNstc = hNexist;
849  }
850  else
851    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
852  stcmem = hCreate((currRing->N) - 1);
853  loop
854  {
855    if (mc)
856    {
857      hComp(hexist, hNexist, mc, hstc, &hNstc);
858      if (!hNstc)
859      {
860        hMu = -1;
861        break;
862      }
863    }
864    hNvar = (currRing->N);
865    for (int i = hNvar; i; i--)
866      hvar[i] = i;
867    hStaircase(hstc, &hNstc, hvar, hNvar);
868    hSupp(hstc, hNstc, hvar, &hNvar);
869    if ((hNvar == (currRing->N)) && (hNstc >= (currRing->N)))
870    {
871      if ((hNvar > 2) && (hNstc > 10))
872        hOrdSupp(hstc, hNstc, hvar, hNvar);
873      memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
874      hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
875      if (hNpure == hNvar)
876      {
877        hLexS(hstc, hNstc, hvar, hNvar);
878        hMu += hZeroMult(hpur0, hstc, hNstc, hvar, hNvar);
879      }
880      else
881        hMu = -1;
882    }
883    else if (hNvar)
884      hMu = -1;
885    mc--;
886    if (mc <= 0 || hMu < 0)
887      break;
888  }
889  hKill(stcmem, (currRing->N) - 1);
890  omFreeSize((ADDRESS)hpur0, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
891  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
892  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
893  hDelete(hexist, hNexist);
894  if (hisModule)
895    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
896}
897
898int  scMult0Int(ideal S, ideal Q)
899{
900  hDegree0(S, Q);
901  return hMu;
902}
903
904
905// HC
906
907static poly pWork;
908
909static void hHedge(poly hEdge)
910{
911  pSetm(pWork);
912  if (pLmCmp(pWork, hEdge) == pOrdSgn)
913  {
914    for (int i = hNvar; i>0; i--)
915      pSetExp(hEdge,i, pGetExp(pWork,i));
916    pSetm(hEdge);
917  }
918}
919
920
921static void hHedgeStep(scmon pure, scfmon stc,
922                       int Nstc, varset var, int Nvar,poly hEdge)
923{
924  int  iv = Nvar -1, k = var[Nvar], a, a0, a1, b, i;
925  int  x, x0;
926  scmon pn;
927  scfmon sn;
928  if (iv==0)
929  {
930    pSetExp(pWork, k, pure[k]);
931    hHedge(hEdge);
932    return;
933  }
934  else if (Nstc==0)
935  {
936    for (i = Nvar; i>0; i--)
937      pSetExp(pWork, var[i], pure[var[i]]);
938    hHedge(hEdge);
939    return;
940  }
941  x = a = 0;
942  pn = hGetpure(pure);
943  sn = hGetmem(Nstc, stc, stcmem[iv]);
944  hStepS(sn, Nstc, var, Nvar, &a, &x);
945  if (a == Nstc)
946  {
947    pSetExp(pWork, k, pure[k]);
948    hHedgeStep(pn, sn, a, var, iv,hEdge);
949    return;
950  }
951  else
952  {
953    pSetExp(pWork, k, x);
954    hHedgeStep(pn, sn, a, var, iv,hEdge);
955  }
956  b = a;
957  loop
958  {
959    a0 = a;
960    x0 = x;
961    hStepS(sn, Nstc, var, Nvar, &a, &x);
962    hElimS(sn, &b, a0, a, var, iv);
963    a1 = a;
964    hPure(sn, a0, &a1, var, iv, pn, &i);
965    hLex2S(sn, b, a0, a1, var, iv, hwork);
966    b += (a1 - a0);
967    if (a < Nstc)
968    {
969      pSetExp(pWork, k, x);
970      hHedgeStep(pn, sn, b, var, iv,hEdge);
971    }
972    else
973    {
974      pSetExp(pWork, k, pure[k]);
975      hHedgeStep(pn, sn, b, var, iv,hEdge);
976      return;
977    }
978  }
979}
980
981void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge, ring tailRing)
982{
983  int  i;
984  int  k = ak;
985  hNvar = (currRing->N);
986  hexist = hInit(S, Q, &hNexist, tailRing);
987  if (k!=0)
988    hComp(hexist, hNexist, k, hexist, &hNstc);
989  else
990    hNstc = hNexist;
991  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
992  hvar = (varset)omAlloc((hNvar + 1) * sizeof(int));
993  hpure = (scmon)omAlloc((1 + (hNvar * hNvar)) * sizeof(int));
994  stcmem = hCreate(hNvar - 1);
995  for (i = hNvar; i>0; i--)
996    hvar[i] = i;
997  hStaircase(hexist, &hNstc, hvar, hNvar);
998  if ((hNvar > 2) && (hNstc > 10))
999    hOrdSupp(hexist, hNstc, hvar, hNvar);
1000  memset(hpure, 0, (hNvar + 1) * sizeof(int));
1001  hPure(hexist, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
1002  hLexS(hexist, hNstc, hvar, hNvar);
1003  if (hEdge!=NULL)
1004    pLmFree(hEdge);
1005  hEdge = pInit();
1006  pWork = pInit();
1007  hHedgeStep(hpure, hexist, hNstc, hvar, hNvar,hEdge);
1008  pSetComp(hEdge,ak);
1009  hKill(stcmem, hNvar - 1);
1010  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1011  omFreeSize((ADDRESS)hvar, (hNvar + 1) * sizeof(int));
1012  omFreeSize((ADDRESS)hpure, (1 + (hNvar * hNvar)) * sizeof(int));
1013  hDelete(hexist, hNexist);
1014  pLmFree(pWork);
1015}
1016
1017
1018
1019//  kbase
1020
1021static poly last;
1022static scmon act;
1023
1024static void scElKbase()
1025{
1026  poly q = pInit();
1027  pSetCoeff0(q,nInit(1));
1028  pSetExpV(q,act);
1029  pNext(q) = NULL;
1030  last = pNext(last) = q;
1031}
1032
1033static int scMax( int i, scfmon stc, int Nvar)
1034{
1035  int x, y=stc[0][Nvar];
1036  for (; i;)
1037  {
1038    i--;
1039    x = stc[i][Nvar];
1040    if (x > y) y = x;
1041  }
1042  return y;
1043}
1044
1045static int scMin( int i, scfmon stc, int Nvar)
1046{
1047  int x, y=stc[0][Nvar];
1048  for (; i;)
1049  {
1050    i--;
1051    x = stc[i][Nvar];
1052    if (x < y) y = x;
1053  }
1054  return y;
1055}
1056
1057static int scRestrict( int &Nstc, scfmon stc, int Nvar)
1058{
1059  int x, y;
1060  int i, j, Istc = Nstc;
1061
1062  y = INT_MAX;
1063  for (i=Nstc-1; i>=0; i--)
1064  {
1065    j = Nvar-1;
1066    loop
1067    {
1068      if(stc[i][j] != 0) break;
1069      j--;
1070      if (j == 0)
1071      {
1072        Istc--;
1073        x = stc[i][Nvar];
1074        if (x < y) y = x;
1075        stc[i] = NULL;
1076        break;
1077      }
1078    }
1079  }
1080  if (Istc < Nstc)
1081  {
1082    for (i=Nstc-1; i>=0; i--)
1083    {
1084      if (stc[i] && (stc[i][Nvar] >= y))
1085      {
1086        Istc--;
1087        stc[i] = NULL;
1088      }
1089    }
1090    j = 0;
1091    while (stc[j]) j++;
1092    i = j+1;
1093    for(; i<Nstc; i++)
1094    {
1095      if (stc[i])
1096      {
1097        stc[j] = stc[i];
1098        j++;
1099      }
1100    }
1101    Nstc = Istc;
1102    return y;
1103  }
1104  else
1105    return -1;
1106}
1107
1108static void scAll( int Nvar, int deg)
1109{
1110  int i;
1111  int d = deg;
1112  if (d == 0)
1113  {
1114    for (i=Nvar; i; i--) act[i] = 0;
1115    scElKbase();
1116    return;
1117  }
1118  if (Nvar == 1)
1119  {
1120    act[1] = d;
1121    scElKbase();
1122    return;
1123  }
1124  do
1125  {
1126    act[Nvar] = d;
1127    scAll(Nvar-1, deg-d);
1128    d--;
1129  } while (d >= 0);
1130}
1131
1132static void scAllKbase( int Nvar, int ideg, int deg)
1133{
1134  do
1135  {
1136    act[Nvar] = ideg;
1137    scAll(Nvar-1, deg-ideg);
1138    ideg--;
1139  } while (ideg >= 0);
1140}
1141
1142static void scDegKbase( scfmon stc, int Nstc, int Nvar, int deg)
1143{
1144  int  Ivar, Istc, i, j;
1145  scfmon sn;
1146  int x, ideg;
1147
1148  if (deg == 0)
1149  {
1150    for (i=Nstc-1; i>=0; i--)
1151    {
1152      for (j=Nvar;j;j--){ if(stc[i][j]) break; }
1153      if (j==0){return;}
1154    }
1155    for (i=Nvar; i; i--) act[i] = 0;
1156    scElKbase();
1157    return;
1158  }
1159  if (Nvar == 1)
1160  {
1161    for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
1162    act[1] = deg;
1163    scElKbase();
1164    return;
1165  }
1166  Ivar = Nvar-1;
1167  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1168  x = scRestrict(Nstc, sn, Nvar);
1169  if (x <= 0)
1170  {
1171    if (x == 0) return;
1172    ideg = deg;
1173  }
1174  else
1175  {
1176    if (deg < x) ideg = deg;
1177    else ideg = x-1;
1178    if (Nstc == 0)
1179    {
1180      scAllKbase(Nvar, ideg, deg);
1181      return;
1182    }
1183  }
1184  loop
1185  {
1186    x = scMax(Nstc, sn, Nvar);
1187    while (ideg >= x)
1188    {
1189      act[Nvar] = ideg;
1190      scDegKbase(sn, Nstc, Ivar, deg-ideg);
1191      ideg--;
1192    }
1193    if (ideg < 0) return;
1194    Istc = Nstc;
1195    for (i=Nstc-1; i>=0; i--)
1196    {
1197      if (ideg < sn[i][Nvar])
1198      {
1199        Istc--;
1200        sn[i] = NULL;
1201      }
1202    }
1203    if (Istc == 0)
1204    {
1205      scAllKbase(Nvar, ideg, deg);
1206      return;
1207    }
1208    j = 0;
1209    while (sn[j]) j++;
1210    i = j+1;
1211    for (; i<Nstc; i++)
1212    {
1213      if (sn[i])
1214      {
1215        sn[j] = sn[i];
1216        j++;
1217      }
1218    }
1219    Nstc = Istc;
1220  }
1221}
1222
1223static void scInKbase( scfmon stc, int Nstc, int Nvar)
1224{
1225  int  Ivar, Istc, i, j;
1226  scfmon sn;
1227  int x, ideg;
1228
1229  if (Nvar == 1)
1230  {
1231    ideg = scMin(Nstc, stc, 1);
1232    while (ideg > 0)
1233    {
1234      ideg--;
1235      act[1] = ideg;
1236      scElKbase();
1237    }
1238    return;
1239  }
1240  Ivar = Nvar-1;
1241  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1242  x = scRestrict(Nstc, sn, Nvar);
1243  if (x == 0) return;
1244  ideg = x-1;
1245  loop
1246  {
1247    x = scMax(Nstc, sn, Nvar);
1248    while (ideg >= x)
1249    {
1250      act[Nvar] = ideg;
1251      scInKbase(sn, Nstc, Ivar);
1252      ideg--;
1253    }
1254    if (ideg < 0) return;
1255    Istc = Nstc;
1256    for (i=Nstc-1; i>=0; i--)
1257    {
1258      if (ideg < sn[i][Nvar])
1259      {
1260        Istc--;
1261        sn[i] = NULL;
1262      }
1263    }
1264    j = 0;
1265    while (sn[j]) j++;
1266    i = j+1;
1267    for (; i<Nstc; i++)
1268    {
1269      if (sn[i])
1270      {
1271        sn[j] = sn[i];
1272        j++;
1273      }
1274    }
1275    Nstc = Istc;
1276  }
1277}
1278
1279static ideal scIdKbase()
1280{
1281  polyset mm;
1282  ideal res;
1283  poly p, q = last;
1284  int i = pLength(q);
1285  res = idInit(i,1);
1286  mm = res->m;
1287  i = 0;
1288  do
1289  {
1290    mm[i] = q;
1291    i++;
1292    p = pNext(q);
1293    pNext(q) = NULL;
1294    q = p;
1295  } while (q!=NULL);
1296  return res;
1297}
1298
1299ideal scKBase(int deg, ideal s, ideal Q, intvec * mv)
1300{
1301  int  i, di;
1302  poly p;
1303
1304  if (deg < 0)
1305  {
1306    di = scDimInt(s, Q);
1307    if (di != 0)
1308    {
1309      //Werror("KBase not finite");
1310      return idInit(1,s->rank);
1311    }
1312  }
1313  stcmem = hCreate((currRing->N) - 1);
1314  hexist = hInit(s, Q, &hNexist);
1315  p = last = pInit();
1316  /*pNext(p) = NULL;*/
1317  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1318  *act = 0;
1319  if (!hNexist)
1320  {
1321    scAll((currRing->N), deg);
1322    goto ende;
1323  }
1324  if (!hisModule)
1325  {
1326    if (deg < 0) scInKbase(hexist, hNexist, (currRing->N));
1327    else scDegKbase(hexist, hNexist, (currRing->N), deg);
1328  }
1329  else
1330  {
1331    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1332    for (i = 1; i <= hisModule; i++)
1333    {
1334      *act = i;
1335      hComp(hexist, hNexist, i, hstc, &hNstc);
1336      int deg_ei=deg;
1337      if (mv!=NULL) deg_ei -= (*mv)[i-1];
1338      if ((deg < 0) || (deg_ei>=0))
1339      {
1340        if (hNstc)
1341        {
1342          if (deg < 0) scInKbase(hstc, hNstc, (currRing->N));
1343          else scDegKbase(hstc, hNstc, (currRing->N), deg_ei);
1344        }
1345        else
1346          scAll((currRing->N), deg_ei);
1347      }
1348    }
1349    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1350  }
1351ende:
1352  hDelete(hexist, hNexist);
1353  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1354  hKill(stcmem, (currRing->N) - 1);
1355  pLmDelete(&p);
1356  if (p == NULL)
1357    return idInit(1,s->rank);
1358  else
1359  {
1360    last = p;
1361    ideal res=scIdKbase();
1362    res->rank=s->rank;
1363    return res;
1364  }
1365}
1366
1367#if 0 //-- alternative implementation of scComputeHC
1368void scComputeHCw(ideal ss, ideal Q, int ak, poly &hEdge, ring tailRing)
1369{
1370  int  i, di;
1371  poly p;
1372
1373  if (hEdge!=NULL)
1374    pLmFree(hEdge);
1375
1376  ideal s=idInit(IDELEMS(ss),ak);
1377  for(i=IDELEMS(ss)-1;i>=0;i--)
1378  {
1379    if (ss->m[i]!=NULL) s->m[i]=pHead(ss->m[i]);
1380  }
1381  di = scDimInt(s, Q);
1382  stcmem = hCreate((currRing->N) - 1);
1383  hexist = hInit(s, Q, &hNexist);
1384  p = last = pInit();
1385  /*pNext(p) = NULL;*/
1386  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1387  *act = 0;
1388  if (!hNexist)
1389  {
1390    scAll((currRing->N), -1);
1391    goto ende;
1392  }
1393  if (!hisModule)
1394  {
1395    scInKbase(hexist, hNexist, (currRing->N));
1396  }
1397  else
1398  {
1399    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1400    for (i = 1; i <= hisModule; i++)
1401    {
1402      *act = i;
1403      hComp(hexist, hNexist, i, hstc, &hNstc);
1404      if (hNstc)
1405      {
1406        scInKbase(hstc, hNstc, (currRing->N));
1407      }
1408      else
1409        scAll((currRing->N), -1);
1410    }
1411    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1412  }
1413ende:
1414  hDelete(hexist, hNexist);
1415  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1416  hKill(stcmem, (currRing->N) - 1);
1417  pDeleteLm(&p);
1418  idDelete(&s);
1419  if (p == NULL)
1420  {
1421    return; // no HEdge
1422  }
1423  else
1424  {
1425    last = p;
1426    ideal res=scIdKbase();
1427    res->rank=ss->rank;
1428    poly p_ind=res->m[0]; int ind=0;
1429    for(i=IDELEMS(res)-1;i>0;i--)
1430    {
1431      if (pCmp(res->m[i],p_ind)==-1) { p_ind=res->m[i]; ind=i; }
1432    }
1433    assume(p_ind!=NULL);
1434    assume(res->m[ind]==p_ind);
1435    hEdge=p_ind;
1436    res->m[ind]=NULL;
1437    nDelete(&pGetCoeff(hEdge));
1438    pGetCoeff(hEdge)=NULL;
1439    for(i=(currRing->N);i>0;i--)
1440      pIncrExp(hEdge,i);
1441    pSetm(hEdge);
1442   
1443    idDelete(&res);
1444    return;
1445  }
1446}
1447#endif
Note: See TracBrowser for help on using the repository browser.