source: git/kernel/hdegree.cc @ 98938c

spielwiese
Last change on this file since 98938c was 30d73e, checked in by Wilfred Pohl <pohl@…>, 20 years ago
fix error in kbase git-svn-id: file:///usr/local/Singular/svn/trunk@7513 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 26.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: hdegree.cc,v 1.3 2004-10-05 11:25:39 pohl 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
809int scDegree(ideal S, ideal Q)
810{
811  hDegree(S, Q);
812  int di = pVariables-hCo;
813  if (di!=0) return hMu;
814  else /*variety is empty or not projective*/ return 0;
815}
816
817static void hDegree0(ideal S, ideal Q)
818{
819  int  mc;
820  hexist = hInit(S, Q, &hNexist);
821  if (!hNexist)
822  {
823    hMu = -1;
824    return;
825  }
826  else
827    hMu = 0;
828  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
829  hvar = (varset)omAlloc((pVariables + 1) * sizeof(int));
830  hpur0 = (scmon)omAlloc((1 + (pVariables * pVariables)) * sizeof(int));
831  mc = hisModule;
832  if (!mc)
833  {
834    hstc = hexist;
835    hNstc = hNexist;
836  }
837  else
838    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
839  stcmem = hCreate(pVariables - 1);
840  loop
841  {
842    if (mc)
843    {
844      hComp(hexist, hNexist, mc, hstc, &hNstc);
845      if (!hNstc)
846      {
847        hMu = -1;
848        break;
849      }
850    }
851    hNvar = pVariables;
852    for (int i = hNvar; i; i--)
853      hvar[i] = i;
854    hStaircase(hstc, &hNstc, hvar, hNvar);
855    hSupp(hstc, hNstc, hvar, &hNvar);
856    if ((hNvar == pVariables) && (hNstc >= pVariables))
857    {
858      if ((hNvar > 2) && (hNstc > 10))
859        hOrdSupp(hstc, hNstc, hvar, hNvar);
860      memset(hpur0, 0, (pVariables + 1) * sizeof(int));
861      hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
862      if (hNpure == hNvar)
863      {
864        hLexS(hstc, hNstc, hvar, hNvar);
865        hMu += hZeroMult(hpur0, hstc, hNstc, hvar, hNvar);
866      }
867      else
868        hMu = -1;
869    }
870    else if (hNvar)
871      hMu = -1;
872    mc--;
873    if (mc <= 0 || hMu < 0)
874      break;
875  }
876  hKill(stcmem, pVariables - 1);
877  omFreeSize((ADDRESS)hpur0, (1 + (pVariables * pVariables)) * sizeof(int));
878  omFreeSize((ADDRESS)hvar, (pVariables + 1) * sizeof(int));
879  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
880  hDelete(hexist, hNexist);
881  if (hisModule)
882    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
883}
884
885int  scMult0Int(ideal S, ideal Q)
886{
887  hDegree0(S, Q);
888  return hMu;
889}
890
891
892// HC
893
894static poly pWork;
895
896static void hHedge(poly hEdge)
897{
898  pSetm(pWork);
899  if (pLmCmp(pWork, hEdge) == pOrdSgn)
900  {
901    for (int i = hNvar; i>0; i--)
902      pSetExp(hEdge,i, pGetExp(pWork,i));
903    pSetm(hEdge);
904  }
905}
906
907
908static void hHedgeStep(scmon pure, scfmon stc,
909                       int Nstc, varset var, int Nvar,poly hEdge)
910{
911  int  iv = Nvar -1, k = var[Nvar], a, a0, a1, b, i;
912  int  x, x0;
913  scmon pn;
914  scfmon sn;
915  if (iv==0)
916  {
917    pSetExp(pWork, k, pure[k]);
918    hHedge(hEdge);
919    return;
920  }
921  else if (Nstc==0)
922  {
923    for (i = Nvar; i>0; i--)
924      pSetExp(pWork, var[i], pure[var[i]]);
925    hHedge(hEdge);
926    return;
927  }
928  x = a = 0;
929  pn = hGetpure(pure);
930  sn = hGetmem(Nstc, stc, stcmem[iv]);
931  hStepS(sn, Nstc, var, Nvar, &a, &x);
932  if (a == Nstc)
933  {
934    pSetExp(pWork, k, pure[k]);
935    hHedgeStep(pn, sn, a, var, iv,hEdge);
936    return;
937  }
938  else
939  {
940    pSetExp(pWork, k, x);
941    hHedgeStep(pn, sn, a, var, iv,hEdge);
942  }
943  b = a;
944  loop
945  {
946    a0 = a;
947    x0 = x;
948    hStepS(sn, Nstc, var, Nvar, &a, &x);
949    hElimS(sn, &b, a0, a, var, iv);
950    a1 = a;
951    hPure(sn, a0, &a1, var, iv, pn, &i);
952    hLex2S(sn, b, a0, a1, var, iv, hwork);
953    b += (a1 - a0);
954    if (a < Nstc)
955    {
956      pSetExp(pWork, k, x);
957      hHedgeStep(pn, sn, b, var, iv,hEdge);
958    }
959    else
960    {
961      pSetExp(pWork, k, pure[k]);
962      hHedgeStep(pn, sn, b, var, iv,hEdge);
963      return;
964    }
965  }
966}
967
968void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge, ring tailRing)
969{
970  int  i;
971  int  k = ak;
972  hNvar = pVariables;
973  hexist = hInit(S, Q, &hNexist, tailRing);
974  if (k!=0)
975    hComp(hexist, hNexist, k, hexist, &hNstc);
976  else
977    hNstc = hNexist;
978  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
979  hvar = (varset)omAlloc((hNvar + 1) * sizeof(int));
980  hpure = (scmon)omAlloc((1 + (hNvar * hNvar)) * sizeof(int));
981  stcmem = hCreate(hNvar - 1);
982  for (i = hNvar; i>0; i--)
983    hvar[i] = i;
984  hStaircase(hexist, &hNstc, hvar, hNvar);
985  if ((hNvar > 2) && (hNstc > 10))
986    hOrdSupp(hexist, hNstc, hvar, hNvar);
987  memset(hpure, 0, (hNvar + 1) * sizeof(int));
988  hPure(hexist, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
989  hLexS(hexist, hNstc, hvar, hNvar);
990  if (hEdge!=NULL)
991    pLmFree(hEdge);
992  hEdge = pInit();
993  pWork = pInit();
994  hHedgeStep(hpure, hexist, hNstc, hvar, hNvar,hEdge);
995  pSetComp(hEdge,ak);
996  hKill(stcmem, hNvar - 1);
997  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
998  omFreeSize((ADDRESS)hvar, (hNvar + 1) * sizeof(int));
999  omFreeSize((ADDRESS)hpure, (1 + (hNvar * hNvar)) * sizeof(int));
1000  hDelete(hexist, hNexist);
1001  pLmFree(pWork);
1002}
1003
1004
1005//  kbase
1006
1007static poly last;
1008static scmon act;
1009
1010static void scElKbase()
1011{
1012  poly q = pInit();
1013  pSetCoeff0(q,nInit(1));
1014  pSetExpV(q,act);
1015  pNext(q) = NULL;
1016  last = pNext(last) = q;
1017}
1018
1019static int scMax( int i, scfmon stc, int Nvar)
1020{
1021  int x, y=stc[0][Nvar];
1022  for (; i;)
1023  {
1024    i--;
1025    x = stc[i][Nvar];
1026    if (x > y) y = x;
1027  }
1028  return y;
1029}
1030
1031static int scMin( 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 scRestrict( int &Nstc, scfmon stc, int Nvar)
1044{
1045  int x, y;
1046  int i, j, Istc = Nstc;
1047
1048  y = MAX_EXPONENT;
1049  for (i=Nstc-1; i>=0; i--)
1050  {
1051    j = Nvar-1;
1052    loop
1053    {
1054      if(stc[i][j] != 0) break;
1055      j--;
1056      if (j == 0)
1057      {
1058        Istc--;
1059        x = stc[i][Nvar];
1060        if (x < y) y = x;
1061        stc[i] = NULL;
1062        break;
1063      }
1064    }
1065  }
1066  if (Istc < Nstc)
1067  {
1068    for (i=Nstc-1; i>=0; i--)
1069    {
1070      if (stc[i] && (stc[i][Nvar] >= y))
1071      {
1072        Istc--;
1073        stc[i] = NULL;
1074      }
1075    }
1076    j = 0;
1077    while (stc[j]) j++;
1078    i = j+1;
1079    for(; i<Nstc; i++)
1080    {
1081      if (stc[i])
1082      {
1083        stc[j] = stc[i];
1084        j++;
1085      }
1086    }
1087    Nstc = Istc;
1088    return y;
1089  }
1090  else
1091    return -1;
1092}
1093
1094static void scAll( int Nvar, int deg)
1095{
1096  int i;
1097  int d = deg;
1098  if (d == 0)
1099  {
1100    for (i=Nvar; i; i--) act[i] = 0;
1101    scElKbase();
1102    return;
1103  }
1104  if (Nvar == 1)
1105  {
1106    act[1] = d;
1107    scElKbase();
1108    return;
1109  }
1110  do
1111  {
1112    act[Nvar] = d;
1113    scAll(Nvar-1, deg-d);
1114    d--;
1115  } while (d >= 0);
1116}
1117
1118static void scAllKbase( int Nvar, int ideg, int deg)
1119{
1120  do
1121  {
1122    act[Nvar] = ideg;
1123    scAll(Nvar-1, deg-ideg);
1124    ideg--;
1125  } while (ideg >= 0);
1126}
1127
1128static void scDegKbase( scfmon stc, int Nstc, int Nvar, int deg)
1129{
1130  int  Ivar, Istc, i, j;
1131  scfmon sn;
1132  int x, ideg;
1133
1134  if (deg == 0)
1135  {
1136    for (i=Nstc-1; i>=0; i--)
1137    {
1138      for (j=Nvar;j;j--){ if(stc[i][j]) break; }
1139      if (j==0){return;}
1140    }
1141    for (i=Nvar; i; i--) act[i] = 0;
1142    scElKbase();
1143    return;
1144  }
1145  if (Nvar == 1)
1146  {
1147    for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
1148    act[1] = deg;
1149    scElKbase();
1150    return;
1151  }
1152  Ivar = Nvar-1;
1153  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1154  x = scRestrict(Nstc, sn, Nvar);
1155  if (x <= 0)
1156  {
1157    if (x == 0) return;
1158    ideg = deg;
1159  }
1160  else
1161  {
1162    if (deg < x) ideg = deg;
1163    else ideg = x-1;
1164    if (Nstc == 0)
1165    {
1166      scAllKbase(Nvar, ideg, deg);
1167      return;
1168    }
1169  }
1170  loop
1171  {
1172    x = scMax(Nstc, sn, Nvar);
1173    while (ideg >= x)
1174    {
1175      act[Nvar] = ideg;
1176      scDegKbase(sn, Nstc, Ivar, deg-ideg);
1177      ideg--;
1178    }
1179    if (ideg < 0) return;
1180    Istc = Nstc;
1181    for (i=Nstc-1; i>=0; i--)
1182    {
1183      if (ideg < sn[i][Nvar])
1184      {
1185        Istc--;
1186        sn[i] = NULL;
1187      }
1188    }
1189    if (Istc == 0)
1190    {
1191      scAllKbase(Nvar, ideg, deg);
1192      return;
1193    }
1194    j = 0;
1195    while (sn[j]) j++;
1196    i = j+1;
1197    for (; i<Nstc; i++)
1198    {
1199      if (sn[i])
1200      {
1201        sn[j] = sn[i];
1202        j++;
1203      }
1204    }
1205    Nstc = Istc;
1206  }
1207}
1208
1209static void scInKbase( scfmon stc, int Nstc, int Nvar)
1210{
1211  int  Ivar, Istc, i, j;
1212  scfmon sn;
1213  int x, ideg;
1214
1215  if (Nvar == 1)
1216  {
1217    ideg = scMin(Nstc, stc, 1);
1218    while (ideg > 0)
1219    {
1220      ideg--;
1221      act[1] = ideg;
1222      scElKbase();
1223    }
1224    return;
1225  }
1226  Ivar = Nvar-1;
1227  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1228  x = scRestrict(Nstc, sn, Nvar);
1229  if (x == 0) return;
1230  ideg = x-1;
1231  loop
1232  {
1233    x = scMax(Nstc, sn, Nvar);
1234    while (ideg >= x)
1235    {
1236      act[Nvar] = ideg;
1237      scInKbase(sn, Nstc, Ivar);
1238      ideg--;
1239    }
1240    if (ideg < 0) return;
1241    Istc = Nstc;
1242    for (i=Nstc-1; i>=0; i--)
1243    {
1244      if (ideg < sn[i][Nvar])
1245      {
1246        Istc--;
1247        sn[i] = NULL;
1248      }
1249    }
1250    j = 0;
1251    while (sn[j]) j++;
1252    i = j+1;
1253    for (; i<Nstc; i++)
1254    {
1255      if (sn[i])
1256      {
1257        sn[j] = sn[i];
1258        j++;
1259      }
1260    }
1261    Nstc = Istc;
1262  }
1263}
1264
1265static ideal scIdKbase()
1266{
1267  polyset mm;
1268  ideal res;
1269  poly p, q = last;
1270  int i = pLength(q);
1271  res = idInit(i,1);
1272  mm = res->m;
1273  i = 0;
1274  do
1275  {
1276    mm[i] = q;
1277    i++;
1278    p = pNext(q);
1279    pNext(q) = NULL;
1280    q = p;
1281  } while (q!=NULL);
1282  return res;
1283}
1284
1285extern ideal scKBase(int deg, ideal s, ideal Q)
1286{
1287  int  i, di;
1288  poly p;
1289
1290  if (deg < 0)
1291  {
1292    di = scDimInt(s, Q);
1293    if (di != 0)
1294    {
1295      //Werror("KBase not finite");
1296      return idInit(1,s->rank);
1297    }
1298  }
1299  stcmem = hCreate(pVariables - 1);
1300  hexist = hInit(s, Q, &hNexist);
1301  p = last = pInit();
1302  /*pNext(p) = NULL;*/
1303  act = (scmon)omAlloc((pVariables + 1) * sizeof(int));
1304  *act = 0;
1305  if (!hNexist)
1306  {
1307    scAll(pVariables, deg);
1308    goto ende;
1309  }
1310  if (!hisModule)
1311  {
1312    if (deg < 0) scInKbase(hexist, hNexist, pVariables);
1313    else scDegKbase(hexist, hNexist, pVariables, deg);
1314  }
1315  else
1316  {
1317    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1318    for (i = 1; i <= hisModule; i++)
1319    {
1320      *act = i;
1321      hComp(hexist, hNexist, i, hstc, &hNstc);
1322      if (hNstc)
1323      {
1324        if (deg < 0) scInKbase(hstc, hNstc, pVariables);
1325        else scDegKbase(hstc, hNstc, pVariables, deg);
1326      }
1327      else
1328        scAll(pVariables, deg);
1329    }
1330    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1331  }
1332ende:
1333  hDelete(hexist, hNexist);
1334  omFreeSize((ADDRESS)act, (pVariables + 1) * sizeof(int));
1335  hKill(stcmem, pVariables - 1);
1336  pDeleteLm(&p);
1337  if (p == NULL)
1338    return idInit(1,s->rank);
1339  else
1340  {
1341    last = p;
1342    ideal res=scIdKbase();
1343    res->rank=s->rank;
1344    return res;
1345  }
1346}
1347
1348
Note: See TracBrowser for help on using the repository browser.