source: git/Singular/hdegree.cc @ 126cfa

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