source: git/kernel/hilb.cc @ abe5c8

fieker-DuValspielwiese
Last change on this file since abe5c8 was 6ce030f, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
removal of the $Id$ svn tag from everywhere NOTE: the git SHA1 may be used instead (only on special places) NOTE: the libraries Singular/LIB/*.lib still contain the marker due to our current use of svn
  • Property mode set to 100644
File size: 9.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5*  ABSTRACT -  Hilbert series
6*/
7
8#include "config.h"
9#include <kernel/mod2.h>
10#include <kernel/structs.h>
11#include <kernel/febase.h>
12#include <omalloc/omalloc.h>
13#include <kernel/polys.h>
14#include <misc/intvec.h>
15#include <kernel/hutil.h>
16#include <kernel/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      if (z>(MAX_INT_VAL)/2)
55      {
56       Werror("interal arrays too big");
57       return;
58      }
59      p = (int *)omAlloc((unsigned long)z * sizeof(int));
60      if (Ql[j]!=0)
61      {
62        if (j==0)
63          memcpy(p, Qpol[j], Ql[j] * sizeof(int));
64        omFreeSize((ADDRESS)Qpol[j], Ql[j] * sizeof(int));
65      }
66      if (j==0)
67      {
68        for (x = Ql[j]; x < z; x++)
69          p[x] = 0;
70      }
71      Ql[j] = z;
72      Qpol[j] = p;
73    }
74  }
75}
76
77static int  *hAddHilb(int Nv, int x, int *pol, int *lp)
78{
79  int  l = *lp, ln, i;
80  int  *pon;
81  *lp = ln = l + x;
82  pon = Qpol[Nv];
83  memcpy(pon, pol, l * sizeof(int));
84  if (l > x)
85  {
86    for (i = x; i < l; i++)
87      pon[i] -= pol[i - x];
88    for (i = l; i < ln; i++)
89      pon[i] = -pol[i - x];
90  }
91  else
92  {
93    for (i = l; i < x; i++)
94      pon[i] = 0;
95    for (i = x; i < ln; i++)
96      pon[i] = -pol[i - x];
97  }
98  return pon;
99}
100
101static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
102{
103  int  l = lp, x, i, j;
104  int  *p, *pl;
105  p = pol;
106  for (i = Nv; i>0; i--)
107  {
108    x = pure[var[i + 1]];
109    if (x!=0)
110      p = hAddHilb(i, x, p, &l);
111  }
112  pl = *Qpol;
113  j = Q0[Nv + 1];
114  for (i = 0; i < l; i++)
115    pl[i + j] += p[i];
116  x = pure[var[1]];
117  if (x!=0)
118  {
119    j += x;
120    for (i = 0; i < l; i++)
121      pl[i + j] -= p[i];
122  }
123  j += l;
124  if (j > hLength)
125    hLength = j;
126}
127
128static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var,
129 int Nvar, int *pol, int Lpol)
130{
131  int  iv = Nvar -1, ln, a, a0, a1, b, i;
132  int  x, x0;
133  scmon pn;
134  scfmon sn;
135  int  *pon;
136  if (Nstc==0)
137  {
138    hLastHilb(pure, iv, var, pol, Lpol);
139    return;
140  }
141  x = a = 0;
142  pn = hGetpure(pure);
143  sn = hGetmem(Nstc, stc, stcmem[iv]);
144  hStepS(sn, Nstc, var, Nvar, &a, &x);
145  Q0[iv] = Q0[Nvar];
146  ln = Lpol;
147  pon = pol;
148  if (a == Nstc)
149  {
150    x = pure[var[Nvar]];
151    if (x!=0)
152      pon = hAddHilb(iv, x, pon, &ln);
153    hHilbStep(pn, sn, a, var, iv, pon, ln);
154    return;
155  }
156  else
157  {
158    pon = hAddHilb(iv, x, pon, &ln);
159    hHilbStep(pn, sn, a, var, iv, pon, ln);
160  }
161  b = a;
162  x0 = 0;
163  loop
164  {
165    Q0[iv] += (x - x0);
166    a0 = a;
167    x0 = x;
168    hStepS(sn, Nstc, var, Nvar, &a, &x);
169    hElimS(sn, &b, a0, a, var, iv);
170    a1 = a;
171    hPure(sn, a0, &a1, var, iv, pn, &i);
172    hLex2S(sn, b, a0, a1, var, iv, hwork);
173    b += (a1 - a0);
174    ln = Lpol;
175    if (a < Nstc)
176    {
177      pon = hAddHilb(iv, x - x0, pol, &ln);
178      hHilbStep(pn, sn, b, var, iv, pon, ln);
179    }
180    else
181    {
182      x = pure[var[Nvar]];
183      if (x!=0)
184        pon = hAddHilb(iv, x - x0, pol, &ln);
185      else
186        pon = pol;
187      hHilbStep(pn, sn, b, var, iv, pon, ln);
188      return;
189    }
190  }
191}
192
193/*
194*basic routines
195*/
196static void hWDegree(intvec *wdegree)
197{
198  int i, k;
199  int x;
200
201  for (i=(currRing->N); i; i--)
202  {
203    x = (*wdegree)[i-1];
204    if (x != 1)
205    {
206      for (k=hNexist-1; k>=0; k--)
207      {
208        hexist[k][i] *= x;
209      }
210    }
211  }
212}
213
214static intvec * hSeries(ideal S, intvec *modulweight,
215                int notstc, intvec *wdegree, ideal Q, ring tailRing)
216{
217  intvec *work, *hseries1=NULL;
218  int  mc;
219  int  p0;
220  int  i, j, k, l, ii, mw;
221  hexist = hInit(S, Q, &hNexist, tailRing);
222  if (hNexist==0)
223  {
224    hseries1=new intvec(2);
225    (*hseries1)[0]=1;
226    (*hseries1)[1]=0;
227    return hseries1;
228  }
229
230  #if 0
231  if (wdegree == NULL)
232    hWeight();
233  else
234    hWDegree(wdegree);
235  #else
236  if (wdegree != NULL) hWDegree(wdegree);
237  #endif
238
239  p0 = 1;
240  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
241  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
242  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
243  stcmem = hCreate((currRing->N) - 1);
244  Qpol = (int **)omAlloc(((currRing->N) + 1) * sizeof(int *));
245  Ql = (int *)omAlloc0(((currRing->N) + 1) * sizeof(int));
246  Q0 = (int *)omAlloc(((currRing->N) + 1) * sizeof(int));
247  *Qpol = NULL;
248  hLength = k = j = 0;
249  mc = hisModule;
250  if (mc!=0)
251  {
252    mw = hMinModulweight(modulweight);
253    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
254  }
255  else
256  {
257    mw = 0;
258    hstc = hexist;
259    hNstc = hNexist;
260  }
261  loop
262  {
263    if (mc!=0)
264    {
265      hComp(hexist, hNexist, mc, hstc, &hNstc);
266      if (modulweight != NULL)
267        j = (*modulweight)[mc-1]-mw;
268    }
269    if (hNstc!=0)
270    {
271      hNvar = (currRing->N);
272      for (i = hNvar; i>=0; i--)
273        hvar[i] = i;
274      //if (notstc) // TODO: no mon divides another
275        hStaircase(hstc, &hNstc, hvar, hNvar);
276      hSupp(hstc, hNstc, hvar, &hNvar);
277      if (hNvar!=0)
278      {
279        if ((hNvar > 2) && (hNstc > 10))
280          hOrdSupp(hstc, hNstc, hvar, hNvar);
281        hHilbEst(hstc, hNstc, hvar, hNvar);
282        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
283        hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
284        hLexS(hstc, hNstc, hvar, hNvar);
285        Q0[hNvar] = 0;
286        hHilbStep(hpure, hstc, hNstc, hvar, hNvar, &p0, 1);
287      }
288    }
289    else
290    {
291      if(*Qpol!=NULL)
292        (**Qpol)++;
293      else
294      {
295        *Qpol = (int *)omAlloc(sizeof(int));
296        hLength = *Ql = **Qpol = 1;
297      }
298    }
299    if (*Qpol!=NULL)
300    {
301      i = hLength;
302      while ((i > 0) && ((*Qpol)[i - 1] == 0))
303        i--;
304      if (i > 0)
305      {
306        l = i + j;
307        if (l > k)
308        {
309          work = new intvec(l);
310          for (ii=0; ii<k; ii++)
311            (*work)[ii] = (*hseries1)[ii];
312          if (hseries1 != NULL)
313            delete hseries1;
314          hseries1 = work;
315          k = l;
316        }
317        while (i > 0)
318        {
319          (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
320          (*Qpol)[i - 1] = 0;
321          i--;
322        }
323      }
324    }
325    mc--;
326    if (mc <= 0)
327      break;
328  }
329  if (k==0)
330  {
331    hseries1=new intvec(2);
332    (*hseries1)[0]=0;
333    (*hseries1)[1]=0;
334  }
335  else
336  {
337    l = k+1;
338    while ((*hseries1)[l-2]==0) l--;
339    if (l!=k)
340    {
341      work = new intvec(l);
342      for (ii=l-2; ii>=0; ii--)
343        (*work)[ii] = (*hseries1)[ii];
344      delete hseries1;
345      hseries1 = work;
346    }
347    (*hseries1)[l-1] = mw;
348  }
349  for (i = 0; i <= (currRing->N); i++)
350  {
351    if (Ql[i]!=0)
352      omFreeSize((ADDRESS)Qpol[i], Ql[i] * sizeof(int));
353  }
354  omFreeSize((ADDRESS)Q0, ((currRing->N) + 1) * sizeof(int));
355  omFreeSize((ADDRESS)Ql, ((currRing->N) + 1) * sizeof(int));
356  omFreeSize((ADDRESS)Qpol, ((currRing->N) + 1) * sizeof(int *));
357  hKill(stcmem, (currRing->N) - 1);
358  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
359  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
360  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
361  hDelete(hexist, hNexist);
362  if (hisModule!=0)
363    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
364  return hseries1;
365}
366
367
368intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
369{
370  return hSeries(S, modulweight, 0, wdegree, Q, tailRing);
371}
372
373intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
374{
375  return hSeries(S, modulweight, 1, wdegree, Q, tailRing);
376}
377
378intvec * hSecondSeries(intvec *hseries1)
379{
380  intvec *work, *hseries2;
381  int i, j, k, s, t, l;
382  if (hseries1 == NULL)
383    return NULL;
384  work = new intvec(hseries1);
385  k = l = work->length()-1;
386  s = 0;
387  for (i = k-1; i >= 0; i--)
388    s += (*work)[i];
389  loop
390  {
391    if ((s != 0) || (k == 1))
392      break;
393    s = 0;
394    t = (*work)[k-1];
395    k--;
396    for (i = k-1; i >= 0; i--)
397    {
398      j = (*work)[i];
399      (*work)[i] = -t;
400      s += t;
401      t += j;
402    }
403  }
404  hseries2 = new intvec(k+1);
405  for (i = k-1; i >= 0; i--)
406    (*hseries2)[i] = (*work)[i];
407  (*hseries2)[k] = (*work)[l];
408  delete work;
409  return hseries2;
410}
411
412void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
413{
414  int m, i, j, k;
415  *co = *mu = 0;
416  if ((s1 == NULL) || (s2 == NULL))
417    return;
418  i = s1->length();
419  j = s2->length();
420  if (j > i)
421    return;
422  m = 0;
423  for(k=j-2; k>=0; k--)
424    m += (*s2)[k];
425  *mu = m;
426  *co = i - j;
427}
428
429static void hPrintHilb(intvec *hseries)
430{
431  int  i, j, l, k;
432  if (hseries == NULL)
433    return;
434  l = hseries->length()-1;
435  k = (*hseries)[l];
436  for (i = 0; i < l; i++)
437  {
438    j = (*hseries)[i];
439    if (j != 0)
440    {
441      Print("//  %8d t^%d\n", j, i+k);
442    }
443  }
444}
445
446/*
447*caller
448*/
449void hLookSeries(ideal S, intvec *modulweight, ideal Q)
450{
451  int co, mu, l;
452  intvec *hseries2;
453  intvec *hseries1 = hFirstSeries(S, modulweight, Q);
454  hPrintHilb(hseries1);
455  l = hseries1->length()-1;
456  if (l > 1)
457    hseries2 = hSecondSeries(hseries1);
458  else
459    hseries2 = hseries1;
460  hDegreeSeries(hseries1, hseries2, &co, &mu);
461  PrintLn();
462  hPrintHilb(hseries2);
463  if ((l == 1) &&(mu == 0))
464    scPrintDegree((currRing->N)+1, 0);
465  else
466    scPrintDegree(co, mu);
467  if (l>1)
468    delete hseries1;
469  delete hseries2;
470}
471
Note: See TracBrowser for help on using the repository browser.