source: git/kernel/hdegree.cc @ 0d0ba9

spielwiese
Last change on this file since 0d0ba9 was 0d0ba9, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: EXPONENT is not used any more git-svn-id: file:///usr/local/Singular/svn/trunk@8769 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 26.4 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: hdegree.cc,v 1.6 2005-10-31 17:35:36 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//  kbase
1018
1019static poly last;
1020static scmon act;
1021
1022static void scElKbase()
1023{
1024  poly q = pInit();
1025  pSetCoeff0(q,nInit(1));
1026  pSetExpV(q,act);
1027  pNext(q) = NULL;
1028  last = pNext(last) = q;
1029}
1030
1031static int scMax( int i, scfmon stc, int Nvar)
1032{
1033  int x, y=stc[0][Nvar];
1034  for (; i;)
1035  {
1036    i--;
1037    x = stc[i][Nvar];
1038    if (x > y) y = x;
1039  }
1040  return y;
1041}
1042
1043static int scMin( int i, scfmon stc, int Nvar)
1044{
1045  int x, y=stc[0][Nvar];
1046  for (; i;)
1047  {
1048    i--;
1049    x = stc[i][Nvar];
1050    if (x < y) y = x;
1051  }
1052  return y;
1053}
1054
1055static int scRestrict( int &Nstc, scfmon stc, int Nvar)
1056{
1057  int x, y;
1058  int i, j, Istc = Nstc;
1059
1060  y = INT_MAX;
1061  for (i=Nstc-1; i>=0; i--)
1062  {
1063    j = Nvar-1;
1064    loop
1065    {
1066      if(stc[i][j] != 0) break;
1067      j--;
1068      if (j == 0)
1069      {
1070        Istc--;
1071        x = stc[i][Nvar];
1072        if (x < y) y = x;
1073        stc[i] = NULL;
1074        break;
1075      }
1076    }
1077  }
1078  if (Istc < Nstc)
1079  {
1080    for (i=Nstc-1; i>=0; i--)
1081    {
1082      if (stc[i] && (stc[i][Nvar] >= y))
1083      {
1084        Istc--;
1085        stc[i] = NULL;
1086      }
1087    }
1088    j = 0;
1089    while (stc[j]) j++;
1090    i = j+1;
1091    for(; i<Nstc; i++)
1092    {
1093      if (stc[i])
1094      {
1095        stc[j] = stc[i];
1096        j++;
1097      }
1098    }
1099    Nstc = Istc;
1100    return y;
1101  }
1102  else
1103    return -1;
1104}
1105
1106static void scAll( int Nvar, int deg)
1107{
1108  int i;
1109  int d = deg;
1110  if (d == 0)
1111  {
1112    for (i=Nvar; i; i--) act[i] = 0;
1113    scElKbase();
1114    return;
1115  }
1116  if (Nvar == 1)
1117  {
1118    act[1] = d;
1119    scElKbase();
1120    return;
1121  }
1122  do
1123  {
1124    act[Nvar] = d;
1125    scAll(Nvar-1, deg-d);
1126    d--;
1127  } while (d >= 0);
1128}
1129
1130static void scAllKbase( int Nvar, int ideg, int deg)
1131{
1132  do
1133  {
1134    act[Nvar] = ideg;
1135    scAll(Nvar-1, deg-ideg);
1136    ideg--;
1137  } while (ideg >= 0);
1138}
1139
1140static void scDegKbase( scfmon stc, int Nstc, int Nvar, int deg)
1141{
1142  int  Ivar, Istc, i, j;
1143  scfmon sn;
1144  int x, ideg;
1145
1146  if (deg == 0)
1147  {
1148    for (i=Nstc-1; i>=0; i--)
1149    {
1150      for (j=Nvar;j;j--){ if(stc[i][j]) break; }
1151      if (j==0){return;}
1152    }
1153    for (i=Nvar; i; i--) act[i] = 0;
1154    scElKbase();
1155    return;
1156  }
1157  if (Nvar == 1)
1158  {
1159    for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
1160    act[1] = deg;
1161    scElKbase();
1162    return;
1163  }
1164  Ivar = Nvar-1;
1165  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1166  x = scRestrict(Nstc, sn, Nvar);
1167  if (x <= 0)
1168  {
1169    if (x == 0) return;
1170    ideg = deg;
1171  }
1172  else
1173  {
1174    if (deg < x) ideg = deg;
1175    else ideg = x-1;
1176    if (Nstc == 0)
1177    {
1178      scAllKbase(Nvar, ideg, deg);
1179      return;
1180    }
1181  }
1182  loop
1183  {
1184    x = scMax(Nstc, sn, Nvar);
1185    while (ideg >= x)
1186    {
1187      act[Nvar] = ideg;
1188      scDegKbase(sn, Nstc, Ivar, deg-ideg);
1189      ideg--;
1190    }
1191    if (ideg < 0) return;
1192    Istc = Nstc;
1193    for (i=Nstc-1; i>=0; i--)
1194    {
1195      if (ideg < sn[i][Nvar])
1196      {
1197        Istc--;
1198        sn[i] = NULL;
1199      }
1200    }
1201    if (Istc == 0)
1202    {
1203      scAllKbase(Nvar, ideg, deg);
1204      return;
1205    }
1206    j = 0;
1207    while (sn[j]) j++;
1208    i = j+1;
1209    for (; i<Nstc; i++)
1210    {
1211      if (sn[i])
1212      {
1213        sn[j] = sn[i];
1214        j++;
1215      }
1216    }
1217    Nstc = Istc;
1218  }
1219}
1220
1221static void scInKbase( scfmon stc, int Nstc, int Nvar)
1222{
1223  int  Ivar, Istc, i, j;
1224  scfmon sn;
1225  int x, ideg;
1226
1227  if (Nvar == 1)
1228  {
1229    ideg = scMin(Nstc, stc, 1);
1230    while (ideg > 0)
1231    {
1232      ideg--;
1233      act[1] = ideg;
1234      scElKbase();
1235    }
1236    return;
1237  }
1238  Ivar = Nvar-1;
1239  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1240  x = scRestrict(Nstc, sn, Nvar);
1241  if (x == 0) return;
1242  ideg = x-1;
1243  loop
1244  {
1245    x = scMax(Nstc, sn, Nvar);
1246    while (ideg >= x)
1247    {
1248      act[Nvar] = ideg;
1249      scInKbase(sn, Nstc, Ivar);
1250      ideg--;
1251    }
1252    if (ideg < 0) return;
1253    Istc = Nstc;
1254    for (i=Nstc-1; i>=0; i--)
1255    {
1256      if (ideg < sn[i][Nvar])
1257      {
1258        Istc--;
1259        sn[i] = NULL;
1260      }
1261    }
1262    j = 0;
1263    while (sn[j]) j++;
1264    i = j+1;
1265    for (; i<Nstc; i++)
1266    {
1267      if (sn[i])
1268      {
1269        sn[j] = sn[i];
1270        j++;
1271      }
1272    }
1273    Nstc = Istc;
1274  }
1275}
1276
1277static ideal scIdKbase()
1278{
1279  polyset mm;
1280  ideal res;
1281  poly p, q = last;
1282  int i = pLength(q);
1283  res = idInit(i,1);
1284  mm = res->m;
1285  i = 0;
1286  do
1287  {
1288    mm[i] = q;
1289    i++;
1290    p = pNext(q);
1291    pNext(q) = NULL;
1292    q = p;
1293  } while (q!=NULL);
1294  return res;
1295}
1296
1297extern ideal scKBase(int deg, ideal s, ideal Q)
1298{
1299  int  i, di;
1300  poly p;
1301
1302  if (deg < 0)
1303  {
1304    di = scDimInt(s, Q);
1305    if (di != 0)
1306    {
1307      //Werror("KBase not finite");
1308      return idInit(1,s->rank);
1309    }
1310  }
1311  stcmem = hCreate(pVariables - 1);
1312  hexist = hInit(s, Q, &hNexist);
1313  p = last = pInit();
1314  /*pNext(p) = NULL;*/
1315  act = (scmon)omAlloc((pVariables + 1) * sizeof(int));
1316  *act = 0;
1317  if (!hNexist)
1318  {
1319    scAll(pVariables, deg);
1320    goto ende;
1321  }
1322  if (!hisModule)
1323  {
1324    if (deg < 0) scInKbase(hexist, hNexist, pVariables);
1325    else scDegKbase(hexist, hNexist, pVariables, deg);
1326  }
1327  else
1328  {
1329    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1330    for (i = 1; i <= hisModule; i++)
1331    {
1332      *act = i;
1333      hComp(hexist, hNexist, i, hstc, &hNstc);
1334      if (hNstc)
1335      {
1336        if (deg < 0) scInKbase(hstc, hNstc, pVariables);
1337        else scDegKbase(hstc, hNstc, pVariables, deg);
1338      }
1339      else
1340        scAll(pVariables, deg);
1341    }
1342    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1343  }
1344ende:
1345  hDelete(hexist, hNexist);
1346  omFreeSize((ADDRESS)act, (pVariables + 1) * sizeof(int));
1347  hKill(stcmem, pVariables - 1);
1348  pDeleteLm(&p);
1349  if (p == NULL)
1350    return idInit(1,s->rank);
1351  else
1352  {
1353    last = p;
1354    ideal res=scIdKbase();
1355    res->rank=s->rank;
1356    return res;
1357  }
1358}
1359
1360
Note: See TracBrowser for help on using the repository browser.