source: git/Singular/hilb.cc @ 82716e

spielwiese
Last change on this file since 82716e was e78cce, checked in by Olaf Bachmann <obachman@…>, 26 years ago
Mon Dec 15 23:31:56 1997 Olaf Bachmann <obachman@mathematik.uni-kl.de> * polys-impl.h: Various small bug - fixes and changes in preparation of COMP_FAST merge: hilbert series, maps are ok; sres and syz are still buggy; small things still need to be checked; git-svn-id: file:///usr/local/Singular/svn/trunk@978 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 9.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: hilb.cc,v 1.7 1997-12-15 22:46:24 obachman Exp $ */
5/*
6*  ABSTRACT -  Hilbert series
7*/
8
9#include "mod2.h"
10#include "tok.h"
11#include "febase.h"
12#include "mmemory.h"
13#include "polys.h"
14#include "intvec.h"
15#include "hutil.h"
16#include "stairc.h"
17
18static int  **Qpol;
19static int  *Q0, *Ql;
20static int  hLength;
21
22static int hMinModulweight(intvec *modulweight)
23{
24  int i,j,k;
25
26  if(modulweight==NULL) return 0;
27  j=(*modulweight)[0];
28  for(i=modulweight->rows()-1;i!=0;i--)
29  {
30    k=(*modulweight)[i];
31    if(k<j) j=k;
32  }
33  return j;
34}
35
36static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
37{
38  int  i, j;
39  int  x, y, z = 1;
40  int  *p;
41  for (i = Nvar; i>0; i--)
42  {
43    x = 0;
44    for (j = 0; j < Nstc; j++)
45    {
46      y = stc[j][var[i]];
47      if (y > x)
48        x = y;
49    }
50    z += x;
51    j = i - 1;
52    if (z > Ql[j])
53    {
54      p = (int *)Alloc(z * sizeof(int));
55      if (Ql[j]!=0)
56      {
57        if (j==0)
58          memcpy(p, Qpol[j], Ql[j] * sizeof(int));
59        Free((ADDRESS)Qpol[j], Ql[j] * sizeof(int));
60      }
61      if (j==0)
62      {
63        for (x = Ql[j]; x < z; x++)
64          p[x] = 0;
65      }
66      Ql[j] = z;
67      Qpol[j] = p;
68    }
69  }
70}
71
72static int  *hAddHilb(int Nv, int x, int *pol, int *lp)
73{
74  int  l = *lp, ln, i;
75  int  *pon;
76  *lp = ln = l + x;
77  pon = Qpol[Nv];
78  memcpy(pon, pol, l * sizeof(int));
79  if (l > x)
80  {
81    for (i = x; i < l; i++)
82      pon[i] -= pol[i - x];
83    for (i = l; i < ln; i++)
84      pon[i] = -pol[i - x];
85  }
86  else
87  {
88    for (i = l; i < x; i++)
89      pon[i] = 0;
90    for (i = x; i < ln; i++)
91      pon[i] = -pol[i - x];
92  }
93  return pon;
94}
95
96static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
97{
98  int  l = lp, x, i, j;
99  int  *p, *pl;
100  p = pol;
101  for (i = Nv; i>0; i--)
102  {
103    x = pure[var[i + 1]];
104    if (x!=0)
105      p = hAddHilb(i, x, p, &l);
106  }
107  pl = *Qpol;
108  j = Q0[Nv + 1];
109  for (i = 0; i < l; i++)
110    pl[i + j] += p[i];
111  x = pure[var[1]];
112  if (x!=0)
113  {
114    j += x;
115    for (i = 0; i < l; i++)
116      pl[i + j] -= p[i];
117  }
118  j += l;
119  if (j > hLength)
120    hLength = j;
121}
122
123static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var,
124 int Nvar, int *pol, int Lpol)
125{
126  int  iv = Nvar -1, ln, a, a0, a1, b, i;
127  Exponent_t  x, x0;
128  scmon pn;
129  scfmon sn;
130  int  *pon;
131  if (Nstc==0)
132  {
133    hLastHilb(pure, iv, var, pol, Lpol);
134    return;
135  }
136  x = a = 0;
137  pn = hGetpure(pure);
138  sn = hGetmem(Nstc, stc, stcmem[iv]);
139  hStepS(sn, Nstc, var, Nvar, &a, &x);
140  Q0[iv] = Q0[Nvar];
141  ln = Lpol;
142  pon = pol;
143  if (a == Nstc)
144  {
145    x = pure[var[Nvar]];
146    if (x!=0)
147      pon = hAddHilb(iv, x, pon, &ln);
148    hHilbStep(pn, sn, a, var, iv, pon, ln);
149    return;
150  }
151  else
152  {
153    pon = hAddHilb(iv, x, pon, &ln);
154    hHilbStep(pn, sn, a, var, iv, pon, ln);
155  }
156  b = a;
157  x0 = 0;
158  loop
159  {
160    Q0[iv] += (x - x0);
161    a0 = a;
162    x0 = x;
163    hStepS(sn, Nstc, var, Nvar, &a, &x);
164    hElimS(sn, &b, a0, a, var, iv);
165    a1 = a;
166    hPure(sn, a0, &a1, var, iv, pn, &i);
167    hLex2S(sn, b, a0, a1, var, iv, hwork);
168    b += (a1 - a0);
169    ln = Lpol;
170    if (a < Nstc)
171    {
172      pon = hAddHilb(iv, x - x0, pol, &ln);
173      hHilbStep(pn, sn, b, var, iv, pon, ln);
174    }
175    else
176    {
177      x = pure[var[Nvar]];
178      if (x!=0)
179        pon = hAddHilb(iv, x - x0, pol, &ln);
180      else
181        pon = pol;
182      hHilbStep(pn, sn, b, var, iv, pon, ln);
183      return;
184    }
185  }
186}
187
188/*
189*basic routines
190*/
191
192static intvec * hSeries(ideal S, intvec *modulweight, ideal Q, int notstc)
193{
194  intvec *work, *hseries1=NULL;
195  Exponent_t  mc;
196  int  *p0;
197  int  i, j, k, l, ii, mw;
198  hexist = hInit(S, Q, &hNexist);
199  if (hNexist==0)
200  {
201    hseries1=new intvec(2);
202    (*hseries1)[0]=1;
203    (*hseries1)[1]=0;
204    return hseries1;
205  }
206  p0 = (int *)Alloc(sizeof(int));
207  *p0 = 1;
208  hwork = (scfmon)Alloc(hNexist * sizeof(scmon));
209  hvar = (varset)Alloc((pVariables + 1) * sizeof(int));
210  hpure = (scmon)Alloc((1 + (pVariables * pVariables)) * sizeof(Exponent_t));
211  stcmem = hCreate(pVariables - 1);
212  Qpol = (int **)Alloc((pVariables + 1) * sizeof(int *));
213  Ql = (int *)Alloc0((pVariables + 1) * sizeof(int));
214  Q0 = (int *)Alloc((pVariables + 1) * sizeof(int));
215  *Qpol = NULL;
216  hLength = k = j = 0;
217  mc = hisModule;
218  if (mc!=0)
219  {
220    mw = hMinModulweight(modulweight);
221    hstc = (scfmon)Alloc(hNexist * sizeof(scmon));
222  }
223  else
224  {
225    mw = 0;
226    hstc = hexist;
227    hNstc = hNexist;
228  }
229  loop
230  {
231    if (mc!=0)
232    {
233      hComp(hexist, hNexist, mc, hstc, &hNstc);
234      if (modulweight != NULL)
235        j = (*modulweight)[mc-1]-mw;
236    }
237    if (hNstc!=0)
238    {
239      hNvar = pVariables;
240      for (i = hNvar; i; i--)
241        hvar[i] = i;
242      if (notstc)
243        hStaircase(hstc, &hNstc, hvar, hNvar);
244      hSupp(hstc, hNstc, hvar, &hNvar);
245      if (hNvar!=0)
246      {
247        if ((hNvar > 2) && (hNstc > 10))
248          hOrdSupp(hstc, hNstc, hvar, hNvar);
249        hHilbEst(hstc, hNstc, hvar, hNvar);
250        memset(hpure, 0, (pVariables + 1) * sizeof(Exponent_t));
251        hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
252        hLexS(hstc, hNstc, hvar, hNvar);
253        Q0[hNvar] = 0;
254        hHilbStep(hpure, hstc, hNstc, hvar, hNvar, p0, 1);
255      }
256    }
257    else
258    {
259      if(*Qpol!=NULL)
260        (**Qpol)++;
261      else
262      {
263        *Qpol = (int *)Alloc(sizeof(int));
264        hLength = *Ql = **Qpol = 1;
265      }
266    }
267    if (*Qpol!=NULL)
268    {
269      i = hLength;
270      while ((i > 0) && ((*Qpol)[i - 1] == 0))
271        i--;
272      if (i > 0)
273      {
274        l = i + j;
275        if (l > k)
276        {
277          work = new intvec(l);
278          for (ii=0; ii<k; ii++)
279            (*work)[ii] = (*hseries1)[ii];
280          if (hseries1 != NULL)
281            delete hseries1;
282          hseries1 = work;
283          k = l;
284        }
285        while (i > 0)
286        {
287          (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
288          (*Qpol)[i - 1] = 0;
289          i--;
290        }
291      }
292    }
293    mc--;
294    if (mc <= 0)
295      break;
296  }
297  if (k==0)
298  {
299    hseries1=new intvec(2);
300    (*hseries1)[0]=0;
301    (*hseries1)[1]=0;
302  }
303  else
304  {
305    l = k+1;
306    while ((*hseries1)[l-2]==0) l--;
307    if (l!=k)
308    {
309      work = new intvec(l);
310      for (ii=l-2; ii>=0; ii--)
311        (*work)[ii] = (*hseries1)[ii];
312      delete hseries1;
313      hseries1 = work;
314    }
315    (*hseries1)[l-1] = mw;
316  }
317  for (i = 0; i <= pVariables; i++)
318  {
319    if (Ql[i]!=0)
320      Free((ADDRESS)Qpol[i], Ql[i] * sizeof(int));
321  }
322  Free((ADDRESS)Q0, (pVariables + 1) * sizeof(int));
323  Free((ADDRESS)Ql, (pVariables + 1) * sizeof(int));
324  Free((ADDRESS)Qpol, (pVariables + 1) * sizeof(int *));
325  hKill(stcmem, pVariables - 1);
326  Free((ADDRESS)hpure, (1 + (pVariables * pVariables)) * sizeof(Exponent_t));
327  Free((ADDRESS)hvar, (pVariables + 1) * sizeof(int));
328  Free((ADDRESS)hwork, hNexist * sizeof(scmon));
329  hDelete(hexist, hNexist);
330  Free((ADDRESS)p0, sizeof(int));
331  if (hisModule!=0)
332    Free((ADDRESS)hstc, hNexist * sizeof(scmon));
333  return hseries1;
334}
335
336
337intvec * hHstdSeries(ideal S, intvec *modulweight, ideal Q)
338{
339  return hSeries(S, modulweight, Q, 0);
340}
341
342intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q)
343{
344  return hSeries(S, modulweight, Q, 1);
345}
346
347intvec * hSecondSeries(intvec *hseries1)
348{
349  intvec *work, *hseries2;
350  int i, j, k, s, t, l;
351  if (hseries1 == NULL)
352    return NULL;
353  work = new intvec(hseries1);
354  k = l = work->length()-1;
355  s = 0;
356  for (i = k-1; i >= 0; i--)
357    s += (*work)[i];
358  loop
359  {
360    if ((s != 0) || (k == 1))
361      break;
362    s = 0;
363    t = (*work)[k-1];
364    k--;
365    for (i = k-1; i >= 0; i--)
366    {
367      j = (*work)[i];
368      (*work)[i] = -t;
369      s += t;
370      t += j;
371    }
372  }
373  hseries2 = new intvec(k+1);
374  for (i = k-1; i >= 0; i--)
375    (*hseries2)[i] = (*work)[i];
376  (*hseries2)[k] = (*work)[l];
377  delete work;
378  return hseries2;
379}
380
381static void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
382{
383  int m, i, j, k;
384  *co = *mu = 0;
385  if ((s1 == NULL) || (s2 == NULL))
386    return;
387  i = s1->length();
388  j = s2->length();
389  if (j > i)
390    return;
391  m = 0;
392  for(k=j-2; k>=0; k--)
393    m += (*s2)[k];
394  *mu = m;
395  *co = i - j;
396}
397
398static void hPrintHilb(intvec *hseries)
399{
400  int  i, j, l, k;
401  if (hseries == NULL)
402    return;
403  l = hseries->length()-1;
404  k = (*hseries)[l];
405  for (i = 0; i < l; i++)
406  {
407    j = (*hseries)[i];
408    if (j != 0)
409    {
410      Print("//  %8d t^%d\n", j, i+k);
411    }
412  }
413}
414
415static void hPrintDegree(int co, int mu)
416{
417  if (pOrdSgn == 1)
418    Print("// codimension = %d\n// dimension   = %d\n// degree      = %d\n",
419      co, pVariables - co, mu);
420  else
421    Print("// codimension  = %d\n// dimension    = %d\n// multiplicity = %d\n",
422      co, pVariables - co, mu);
423}
424
425/*
426*caller
427*/
428void hLookSeries(ideal S, intvec *modulweight, ideal Q)
429{
430  int co, mu, l;
431  intvec *hseries2;
432  intvec *hseries1 = hFirstSeries(S, modulweight, Q);
433  hPrintHilb(hseries1);
434  l = hseries1->length()-1;
435  if (l > 1)
436    hseries2 = hSecondSeries(hseries1);
437  else
438    hseries2 = hseries1;
439  hDegreeSeries(hseries1, hseries2, &co, &mu);
440  PrintLn();
441  hPrintHilb(hseries2);
442  if ((l == 1) &&(mu == 0))
443    hPrintDegree(pVariables+1, 0);
444  else
445    hPrintDegree(co, mu);
446  if (l>1)
447    delete hseries1;
448  delete hseries2;
449}
450
Note: See TracBrowser for help on using the repository browser.