source: git/Singular/hdegree.cc @ 5ca9807

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