source: git/Singular/hdegree.cc @ c232af

spielwiese
Last change on this file since c232af was c232af, checked in by Olaf Bachmann <obachman@…>, 24 years ago
* omalloc stuff git-svn-id: file:///usr/local/Singular/svn/trunk@4524 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.20 2000-08-14 12:56:18 obachman 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 (pComp0(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
1050
1051void scComputeHC(ideal S, int ak, poly &hEdge)
1052{
1053  int  i;
1054  Exponent_t  k = ak;
1055  hNvar = pVariables;
1056  hexist = hInit(S, NULL, &hNexist);
1057  if (k!=0)
1058    hComp(hexist, hNexist, k, hexist, &hNstc);
1059  else
1060    hNstc = hNexist;
1061  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
1062  hvar = (varset)omAlloc((hNvar + 1) * sizeof(int));
1063  hpure = (scmon)omAlloc((1 + (hNvar * hNvar)) * sizeof(Exponent_t));
1064  stcmem = hCreate(hNvar - 1);
1065  for (i = hNvar; i>0; i--)
1066    hvar[i] = i;
1067  hStaircase(hexist, &hNstc, hvar, hNvar);
1068  if ((hNvar > 2) && (hNstc > 10))
1069    hOrdSupp(hexist, hNstc, hvar, hNvar);
1070  memset(hpure, 0, (hNvar + 1) * sizeof(Exponent_t));
1071  hPure(hexist, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
1072  hLexS(hexist, hNstc, hvar, hNvar);
1073  if (hEdge!=NULL)
1074    pFree1(hEdge);
1075  hEdge = pInit();
1076  pWork = pInit();
1077  hHedgeStep(hpure, hexist, hNstc, hvar, hNvar,hEdge);
1078  pSetComp(hEdge,ak);
1079  hKill(stcmem, hNvar - 1);
1080  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1081  omFreeSize((ADDRESS)hvar, (hNvar + 1) * sizeof(int));
1082  omFreeSize((ADDRESS)hpure, (1 + (hNvar * hNvar)) * sizeof(Exponent_t));
1083  hDelete(hexist, hNexist);
1084  pFree1(pWork);
1085}
1086
1087
1088//  kbase
1089
1090static poly last;
1091static scmon act;
1092
1093static void scElKbase()
1094{
1095  poly q = pInit();
1096  pSetCoeff0(q,nInit(1));
1097  pSetExpV(q,act);
1098  pNext(q) = NULL;
1099  last = pNext(last) = q;
1100}
1101
1102static Exponent_t scMax( int i, scfmon stc, int Nvar)
1103{
1104  Exponent_t x, y=stc[0][Nvar];
1105  for (; i;)
1106  {
1107    i--;
1108    x = stc[i][Nvar];
1109    if (x > y) y = x;
1110  }
1111  return y;
1112}
1113
1114static Exponent_t scMin( int i, scfmon stc, int Nvar)
1115{
1116  Exponent_t x, y=stc[0][Nvar];
1117  for (; i;)
1118  {
1119    i--;
1120    x = stc[i][Nvar];
1121    if (x < y) y = x;
1122  }
1123  return y;
1124}
1125
1126static Exponent_t scRestrict( int &Nstc, scfmon stc, int Nvar)
1127{
1128  Exponent_t x, y;
1129  int i, j, Istc = Nstc;
1130
1131  y = MAX_EXPONENT;
1132  for (i=Nstc-1; i>=0; i--)
1133  {
1134    j = Nvar-1;
1135    loop
1136    {
1137      if(stc[i][j] != 0) break;
1138      j--;
1139      if (j == 0)
1140      {
1141        Istc--;
1142        x = stc[i][Nvar];
1143        if (x < y) y = x;
1144        stc[i] = NULL;
1145        break;
1146      }
1147    }
1148  }
1149  if (Istc < Nstc)
1150  {
1151    for (i=Nstc-1; i>=0; i--)
1152    {
1153      if (stc[i] && (stc[i][Nvar] >= y))
1154      {
1155        Istc--;
1156        stc[i] = NULL;
1157      }
1158    }
1159    j = 0;
1160    while (stc[j]) j++;
1161    i = j+1;
1162    for(; i<Nstc; i++)
1163    {
1164      if (stc[i])
1165      {
1166        stc[j] = stc[i];
1167        j++;
1168      }
1169    }
1170    Nstc = Istc;
1171    return y;
1172  }
1173  else
1174    return -1;
1175}
1176
1177static void scAll( int Nvar, Exponent_t deg)
1178{
1179  int i;
1180  Exponent_t d = deg;
1181  if (d == 0)
1182  {
1183    for (i=Nvar; i; i--) act[i] = 0;
1184    scElKbase();
1185    return;
1186  }
1187  if (Nvar == 1)
1188  {
1189    act[1] = d;
1190    scElKbase();
1191    return;
1192  }
1193  do
1194  {
1195    act[Nvar] = d;
1196    scAll(Nvar-1, deg-d);
1197    d--;
1198  } while (d >= 0);
1199}
1200
1201static void scAllKbase( int Nvar, Exponent_t ideg, Exponent_t deg)
1202{
1203  do
1204  {
1205    act[Nvar] = ideg;
1206    scAll(Nvar-1, deg-ideg);
1207    ideg--;
1208  } while (ideg >= 0);
1209}
1210
1211static void scDegKbase( scfmon stc, int Nstc, int Nvar, Exponent_t deg)
1212{
1213  int  Ivar, Istc, i, j;
1214  scfmon sn;
1215  Exponent_t x, ideg;
1216
1217  if (deg == 0)
1218  {
1219    for (i=Nvar; i; i--) act[i] = 0;
1220    scElKbase();
1221    return;
1222  }
1223  if (Nvar == 1)
1224  {
1225    for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
1226    act[1] = deg;
1227    scElKbase();
1228    return;
1229  }
1230  Ivar = Nvar-1;
1231  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1232  x = scRestrict(Nstc, sn, Nvar);
1233  if (x <= 0)
1234  {
1235    if (x == 0) return;
1236    ideg = deg;
1237  }
1238  else
1239  {
1240    if (deg < x) ideg = deg;
1241    else ideg = x-1;
1242    if (Nstc == 0)
1243    {
1244      scAllKbase(Nvar, ideg, deg);
1245      return;
1246    }
1247  }
1248  loop
1249  {
1250    x = scMax(Nstc, sn, Nvar);
1251    while (ideg >= x)
1252    {
1253      act[Nvar] = ideg;
1254      scDegKbase(sn, Nstc, Ivar, deg-ideg);
1255      ideg--;
1256    }
1257    if (ideg < 0) return;
1258    Istc = Nstc;
1259    for (i=Nstc-1; i>=0; i--)
1260    {
1261      if (ideg < sn[i][Nvar])
1262      {
1263        Istc--;
1264        sn[i] = NULL;
1265      }
1266    }
1267    if (Istc == 0)
1268    {
1269      scAllKbase(Nvar, ideg, deg);
1270      return;
1271    }
1272    j = 0;
1273    while (sn[j]) j++;
1274    i = j+1;
1275    for (; i<Nstc; i++)
1276    {
1277      if (sn[i])
1278      {
1279        sn[j] = sn[i];
1280        j++;
1281      }
1282    }
1283    Nstc = Istc;
1284  }
1285}
1286
1287static void scInKbase( scfmon stc, int Nstc, int Nvar)
1288{
1289  int  Ivar, Istc, i, j;
1290  scfmon sn;
1291  Exponent_t x, ideg;
1292
1293  if (Nvar == 1)
1294  {
1295    ideg = scMin(Nstc, stc, 1);
1296    while (ideg > 0)
1297    {
1298      ideg--;
1299      act[1] = ideg;
1300      scElKbase();
1301    }
1302    return;
1303  }
1304  Ivar = Nvar-1;
1305  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1306  x = scRestrict(Nstc, sn, Nvar);
1307  if (x == 0) return;
1308  ideg = x-1;
1309  loop
1310  {
1311    x = scMax(Nstc, sn, Nvar);
1312    while (ideg >= x)
1313    {
1314      act[Nvar] = ideg;
1315      scInKbase(sn, Nstc, Ivar);
1316      ideg--;
1317    }
1318    if (ideg < 0) return;
1319    Istc = Nstc;
1320    for (i=Nstc-1; i>=0; i--)
1321    {
1322      if (ideg < sn[i][Nvar])
1323      {
1324        Istc--;
1325        sn[i] = NULL;
1326      }
1327    }
1328    j = 0;
1329    while (sn[j]) j++;
1330    i = j+1;
1331    for (; i<Nstc; i++)
1332    {
1333      if (sn[i])
1334      {
1335        sn[j] = sn[i];
1336        j++;
1337      }
1338    }
1339    Nstc = Istc;
1340  }
1341}
1342
1343static ideal scIdKbase()
1344{
1345  polyset mm;
1346  ideal res;
1347  poly p, q = last;
1348  int i = pLength(q);
1349  res = idInit(i,1);
1350  mm = res->m;
1351  i = 0;
1352  do
1353  {
1354    mm[i] = q;
1355    i++;
1356    p = pNext(q);
1357    pNext(q) = NULL;
1358    q = p;
1359  } while (q!=NULL);
1360  return res;
1361}
1362
1363extern ideal scKBase(int deg, ideal s, ideal Q)
1364{
1365  int  i, di;
1366  poly p;
1367
1368  if (deg < 0)
1369  {
1370    di = scDimInt(s, Q);
1371    if (di != 0)
1372    {
1373      //Werror("KBase not finite");
1374      return idInit(1,s->rank);
1375    }
1376  }
1377  stcmem = hCreate(pVariables - 1);
1378  hexist = hInit(s, Q, &hNexist);
1379  p = last = pInit();
1380  /*pNext(p) = NULL;*/
1381  act = (scmon)omAlloc((pVariables + 1) * sizeof(Exponent_t));
1382  *act = 0;
1383  if (!hNexist)
1384  {
1385    scAll(pVariables, deg);
1386    goto ende;
1387  }
1388  if (!hisModule)
1389  {
1390    if (deg < 0) scInKbase(hexist, hNexist, pVariables);
1391    else scDegKbase(hexist, hNexist, pVariables, deg);
1392  }
1393  else
1394  {
1395    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1396    for (i = 1; i <= hisModule; i++)
1397    {
1398      *act = i;
1399      hComp(hexist, hNexist, i, hstc, &hNstc);
1400      if (hNstc)
1401      {
1402        if (deg < 0) scInKbase(hstc, hNstc, pVariables);
1403        else scDegKbase(hstc, hNstc, pVariables, deg);
1404      }
1405      else
1406        scAll(pVariables, deg);
1407    }
1408    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1409  }
1410ende:
1411  hDelete(hexist, hNexist);
1412  omFreeSize((ADDRESS)act, (pVariables + 1) * sizeof(Exponent_t));
1413  hKill(stcmem, pVariables - 1);
1414  pDelete1(&p);
1415  if (p == NULL)
1416    return idInit(1,s->rank);
1417  else
1418  {
1419    last = p;
1420    ideal res=scIdKbase();
1421    res->rank=s->rank;
1422    return res;
1423  }
1424}
1425
1426
Note: See TracBrowser for help on using the repository browser.