source: git/kernel/hdegree.cc @ fbc7cb

jengelh-datetimespielwiese
Last change on this file since fbc7cb was a9c298, checked in by Hans Schoenemann <hannes@…>, 9 years ago
format stuff
  • Property mode set to 100644
File size: 29.0 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 "singularconfig.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
988  #if HAVE_RINGS
989  if (rField_is_Ring(currRing) && (currRing->OrdSgn == -1))
990  {
991  //consider just monic generators (over rings with zero-divisors)
992  ideal SS=id_Copy(S,tailRing);
993  for(i=0;i<=idElem(SS);i++)
994        {
995        if(pIsPurePower(SS->m[i])==0)
996                p_Delete(&SS->m[i],tailRing);
997        }
998        S=id_Copy(SS,tailRing);
999  }
1000  #endif
1001
1002  hNvar = (currRing->N);
1003  hexist = hInit(S, Q, &hNexist, tailRing);
1004  if (k!=0)
1005    hComp(hexist, hNexist, k, hexist, &hNstc);
1006  else
1007    hNstc = hNexist;
1008  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
1009  hvar = (varset)omAlloc((hNvar + 1) * sizeof(int));
1010  hpure = (scmon)omAlloc((1 + (hNvar * hNvar)) * sizeof(int));
1011  stcmem = hCreate(hNvar - 1);
1012  for (i = hNvar; i>0; i--)
1013    hvar[i] = i;
1014  hStaircase(hexist, &hNstc, hvar, hNvar);
1015  if ((hNvar > 2) && (hNstc > 10))
1016    hOrdSupp(hexist, hNstc, hvar, hNvar);
1017  memset(hpure, 0, (hNvar + 1) * sizeof(int));
1018  hPure(hexist, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
1019  hLexS(hexist, hNstc, hvar, hNvar);
1020  if (hEdge!=NULL)
1021    pLmFree(hEdge);
1022  hEdge = pInit();
1023  pWork = pInit();
1024  hHedgeStep(hpure, hexist, hNstc, hvar, hNvar,hEdge);
1025  pSetComp(hEdge,ak);
1026  hKill(stcmem, hNvar - 1);
1027  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1028  omFreeSize((ADDRESS)hvar, (hNvar + 1) * sizeof(int));
1029  omFreeSize((ADDRESS)hpure, (1 + (hNvar * hNvar)) * sizeof(int));
1030  hDelete(hexist, hNexist);
1031  pLmFree(pWork);
1032}
1033
1034
1035
1036//  kbase
1037
1038static poly last;
1039static scmon act;
1040
1041static void scElKbase()
1042{
1043  poly q = pInit();
1044  pSetCoeff0(q,nInit(1));
1045  pSetExpV(q,act);
1046  pNext(q) = NULL;
1047  last = pNext(last) = q;
1048}
1049
1050static int scMax( int i, scfmon stc, int Nvar)
1051{
1052  int x, y=stc[0][Nvar];
1053  for (; i;)
1054  {
1055    i--;
1056    x = stc[i][Nvar];
1057    if (x > y) y = x;
1058  }
1059  return y;
1060}
1061
1062static int scMin( int i, scfmon stc, int Nvar)
1063{
1064  int x, y=stc[0][Nvar];
1065  for (; i;)
1066  {
1067    i--;
1068    x = stc[i][Nvar];
1069    if (x < y) y = x;
1070  }
1071  return y;
1072}
1073
1074static int scRestrict( int &Nstc, scfmon stc, int Nvar)
1075{
1076  int x, y;
1077  int i, j, Istc = Nstc;
1078
1079  y = MAX_INT_VAL;
1080  for (i=Nstc-1; i>=0; i--)
1081  {
1082    j = Nvar-1;
1083    loop
1084    {
1085      if(stc[i][j] != 0) break;
1086      j--;
1087      if (j == 0)
1088      {
1089        Istc--;
1090        x = stc[i][Nvar];
1091        if (x < y) y = x;
1092        stc[i] = NULL;
1093        break;
1094      }
1095    }
1096  }
1097  if (Istc < Nstc)
1098  {
1099    for (i=Nstc-1; i>=0; i--)
1100    {
1101      if (stc[i] && (stc[i][Nvar] >= y))
1102      {
1103        Istc--;
1104        stc[i] = NULL;
1105      }
1106    }
1107    j = 0;
1108    while (stc[j]) j++;
1109    i = j+1;
1110    for(; i<Nstc; i++)
1111    {
1112      if (stc[i])
1113      {
1114        stc[j] = stc[i];
1115        j++;
1116      }
1117    }
1118    Nstc = Istc;
1119    return y;
1120  }
1121  else
1122    return -1;
1123}
1124
1125static void scAll( int Nvar, int deg)
1126{
1127  int i;
1128  int d = deg;
1129  if (d == 0)
1130  {
1131    for (i=Nvar; i; i--) act[i] = 0;
1132    scElKbase();
1133    return;
1134  }
1135  if (Nvar == 1)
1136  {
1137    act[1] = d;
1138    scElKbase();
1139    return;
1140  }
1141  do
1142  {
1143    act[Nvar] = d;
1144    scAll(Nvar-1, deg-d);
1145    d--;
1146  } while (d >= 0);
1147}
1148
1149static void scAllKbase( int Nvar, int ideg, int deg)
1150{
1151  do
1152  {
1153    act[Nvar] = ideg;
1154    scAll(Nvar-1, deg-ideg);
1155    ideg--;
1156  } while (ideg >= 0);
1157}
1158
1159static void scDegKbase( scfmon stc, int Nstc, int Nvar, int deg)
1160{
1161  int  Ivar, Istc, i, j;
1162  scfmon sn;
1163  int x, ideg;
1164
1165  if (deg == 0)
1166  {
1167    for (i=Nstc-1; i>=0; i--)
1168    {
1169      for (j=Nvar;j;j--){ if(stc[i][j]) break; }
1170      if (j==0){return;}
1171    }
1172    for (i=Nvar; i; i--) act[i] = 0;
1173    scElKbase();
1174    return;
1175  }
1176  if (Nvar == 1)
1177  {
1178    for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
1179    act[1] = deg;
1180    scElKbase();
1181    return;
1182  }
1183  Ivar = Nvar-1;
1184  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1185  x = scRestrict(Nstc, sn, Nvar);
1186  if (x <= 0)
1187  {
1188    if (x == 0) return;
1189    ideg = deg;
1190  }
1191  else
1192  {
1193    if (deg < x) ideg = deg;
1194    else ideg = x-1;
1195    if (Nstc == 0)
1196    {
1197      scAllKbase(Nvar, ideg, deg);
1198      return;
1199    }
1200  }
1201  loop
1202  {
1203    x = scMax(Nstc, sn, Nvar);
1204    while (ideg >= x)
1205    {
1206      act[Nvar] = ideg;
1207      scDegKbase(sn, Nstc, Ivar, deg-ideg);
1208      ideg--;
1209    }
1210    if (ideg < 0) return;
1211    Istc = Nstc;
1212    for (i=Nstc-1; i>=0; i--)
1213    {
1214      if (ideg < sn[i][Nvar])
1215      {
1216        Istc--;
1217        sn[i] = NULL;
1218      }
1219    }
1220    if (Istc == 0)
1221    {
1222      scAllKbase(Nvar, ideg, deg);
1223      return;
1224    }
1225    j = 0;
1226    while (sn[j]) j++;
1227    i = j+1;
1228    for (; i<Nstc; i++)
1229    {
1230      if (sn[i])
1231      {
1232        sn[j] = sn[i];
1233        j++;
1234      }
1235    }
1236    Nstc = Istc;
1237  }
1238}
1239
1240static void scInKbase( scfmon stc, int Nstc, int Nvar)
1241{
1242  int  Ivar, Istc, i, j;
1243  scfmon sn;
1244  int x, ideg;
1245
1246  if (Nvar == 1)
1247  {
1248    ideg = scMin(Nstc, stc, 1);
1249    while (ideg > 0)
1250    {
1251      ideg--;
1252      act[1] = ideg;
1253      scElKbase();
1254    }
1255    return;
1256  }
1257  Ivar = Nvar-1;
1258  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1259  x = scRestrict(Nstc, sn, Nvar);
1260  if (x == 0) return;
1261  ideg = x-1;
1262  loop
1263  {
1264    x = scMax(Nstc, sn, Nvar);
1265    while (ideg >= x)
1266    {
1267      act[Nvar] = ideg;
1268      scInKbase(sn, Nstc, Ivar);
1269      ideg--;
1270    }
1271    if (ideg < 0) return;
1272    Istc = Nstc;
1273    for (i=Nstc-1; i>=0; i--)
1274    {
1275      if (ideg < sn[i][Nvar])
1276      {
1277        Istc--;
1278        sn[i] = NULL;
1279      }
1280    }
1281    j = 0;
1282    while (sn[j]) j++;
1283    i = j+1;
1284    for (; i<Nstc; i++)
1285    {
1286      if (sn[i])
1287      {
1288        sn[j] = sn[i];
1289        j++;
1290      }
1291    }
1292    Nstc = Istc;
1293  }
1294}
1295
1296static ideal scIdKbase()
1297{
1298  polyset mm;
1299  ideal res;
1300  poly p, q = last;
1301  int i = pLength(q);
1302  res = idInit(i,1);
1303  mm = res->m;
1304  i = 0;
1305  do
1306  {
1307    mm[i] = q;
1308    i++;
1309    p = pNext(q);
1310    pNext(q) = NULL;
1311    q = p;
1312  } while (q!=NULL);
1313  return res;
1314}
1315
1316ideal scKBase(int deg, ideal s, ideal Q, intvec * mv)
1317{
1318  int  i, di;
1319  poly p;
1320
1321  if (deg < 0)
1322  {
1323    di = scDimInt(s, Q);
1324    if (di != 0)
1325    {
1326      //Werror("KBase not finite");
1327      return idInit(1,s->rank);
1328    }
1329  }
1330  stcmem = hCreate((currRing->N) - 1);
1331  hexist = hInit(s, Q, &hNexist);
1332  p = last = pInit();
1333  /*pNext(p) = NULL;*/
1334  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1335  *act = 0;
1336  if (!hNexist)
1337  {
1338    scAll((currRing->N), deg);
1339    goto ende;
1340  }
1341  if (!hisModule)
1342  {
1343    if (deg < 0) scInKbase(hexist, hNexist, (currRing->N));
1344    else scDegKbase(hexist, hNexist, (currRing->N), deg);
1345  }
1346  else
1347  {
1348    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1349    for (i = 1; i <= hisModule; i++)
1350    {
1351      *act = i;
1352      hComp(hexist, hNexist, i, hstc, &hNstc);
1353      int deg_ei=deg;
1354      if (mv!=NULL) deg_ei -= (*mv)[i-1];
1355      if ((deg < 0) || (deg_ei>=0))
1356      {
1357        if (hNstc)
1358        {
1359          if (deg < 0) scInKbase(hstc, hNstc, (currRing->N));
1360          else scDegKbase(hstc, hNstc, (currRing->N), deg_ei);
1361        }
1362        else
1363          scAll((currRing->N), deg_ei);
1364      }
1365    }
1366    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1367  }
1368ende:
1369  hDelete(hexist, hNexist);
1370  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1371  hKill(stcmem, (currRing->N) - 1);
1372  pLmDelete(&p);
1373  if (p == NULL)
1374    return idInit(1,s->rank);
1375  else
1376  {
1377    last = p;
1378    ideal res=scIdKbase();
1379    res->rank=s->rank;
1380    return res;
1381  }
1382}
1383
1384#if 0 //-- alternative implementation of scComputeHC
1385void scComputeHCw(ideal ss, ideal Q, int ak, poly &hEdge, ring tailRing)
1386{
1387  int  i, di;
1388  poly p;
1389
1390  if (hEdge!=NULL)
1391    pLmFree(hEdge);
1392
1393  ideal s=idInit(IDELEMS(ss),ak);
1394  for(i=IDELEMS(ss)-1;i>=0;i--)
1395  {
1396    if (ss->m[i]!=NULL) s->m[i]=pHead(ss->m[i]);
1397  }
1398  di = scDimInt(s, Q);
1399  stcmem = hCreate((currRing->N) - 1);
1400  hexist = hInit(s, Q, &hNexist);
1401  p = last = pInit();
1402  /*pNext(p) = NULL;*/
1403  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1404  *act = 0;
1405  if (!hNexist)
1406  {
1407    scAll((currRing->N), -1);
1408    goto ende;
1409  }
1410  if (!hisModule)
1411  {
1412    scInKbase(hexist, hNexist, (currRing->N));
1413  }
1414  else
1415  {
1416    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1417    for (i = 1; i <= hisModule; i++)
1418    {
1419      *act = i;
1420      hComp(hexist, hNexist, i, hstc, &hNstc);
1421      if (hNstc)
1422      {
1423        scInKbase(hstc, hNstc, (currRing->N));
1424      }
1425      else
1426        scAll((currRing->N), -1);
1427    }
1428    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1429  }
1430ende:
1431  hDelete(hexist, hNexist);
1432  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1433  hKill(stcmem, (currRing->N) - 1);
1434  pDeleteLm(&p);
1435  idDelete(&s);
1436  if (p == NULL)
1437  {
1438    return; // no HEdge
1439  }
1440  else
1441  {
1442    last = p;
1443    ideal res=scIdKbase();
1444    res->rank=ss->rank;
1445    poly p_ind=res->m[0]; int ind=0;
1446    for(i=IDELEMS(res)-1;i>0;i--)
1447    {
1448      if (pCmp(res->m[i],p_ind)==-1) { p_ind=res->m[i]; ind=i; }
1449    }
1450    assume(p_ind!=NULL);
1451    assume(res->m[ind]==p_ind);
1452    hEdge=p_ind;
1453    res->m[ind]=NULL;
1454    nDelete(&pGetCoeff(hEdge));
1455    pGetCoeff(hEdge)=NULL;
1456    for(i=(currRing->N);i>0;i--)
1457      pIncrExp(hEdge,i);
1458    pSetm(hEdge);
1459
1460    idDelete(&res);
1461    return;
1462  }
1463}
1464#endif
Note: See TracBrowser for help on using the repository browser.