source: git/Singular/hdegree.cc @ 907274

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