source: git/kernel/hdegree.cc @ 16f511

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