source: git/kernel/hdegree.cc @ 87beab7

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