source: git/Singular/hdegree.cc @ d2b2a7

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