source: git/kernel/hilb.cc @ fbc7cb

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