source: git/kernel/hilb.cc @ ebf6a6

spielwiese
Last change on this file since ebf6a6 was ebf6a6, checked in by Hans Schönemann <hannes@…>, 14 years ago
*hannes: malloc overflow git-svn-id: file:///usr/local/Singular/svn/trunk@11947 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 9.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: hilb.cc,v 1.4 2009-07-06 11:54:10 Singular Exp $ */
5/*
6*  ABSTRACT -  Hilbert series
7*/
8
9#include "mod2.h"
10#include "structs.h"
11#include "febase.h"
12#include "omalloc.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      if (z>(INT_MAX)/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=pVariables; 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 = (int *)omAllocBin(int_bin);
240  *p0 = 1;
241  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
242  hvar = (varset)omAlloc((pVariables + 1) * sizeof(int));
243  hpure = (scmon)omAlloc((1 + (pVariables * pVariables)) * sizeof(int));
244  stcmem = hCreate(pVariables - 1);
245  Qpol = (int **)omAlloc((pVariables + 1) * sizeof(int *));
246  Ql = (int *)omAlloc0((pVariables + 1) * sizeof(int));
247  Q0 = (int *)omAlloc((pVariables + 1) * sizeof(int));
248  *Qpol = NULL;
249  hLength = k = j = 0;
250  mc = hisModule;
251  if (mc!=0)
252  {
253    mw = hMinModulweight(modulweight);
254    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
255  }
256  else
257  {
258    mw = 0;
259    hstc = hexist;
260    hNstc = hNexist;
261  }
262  loop
263  {
264    if (mc!=0)
265    {
266      hComp(hexist, hNexist, mc, hstc, &hNstc);
267      if (modulweight != NULL)
268        j = (*modulweight)[mc-1]-mw;
269    }
270    if (hNstc!=0)
271    {
272      hNvar = pVariables;
273      for (i = hNvar; i>=0; i--)
274        hvar[i] = i;
275      //if (notstc) // TODO: no mon divides another
276        hStaircase(hstc, &hNstc, hvar, hNvar);
277      hSupp(hstc, hNstc, hvar, &hNvar);
278      if (hNvar!=0)
279      {
280        if ((hNvar > 2) && (hNstc > 10))
281          hOrdSupp(hstc, hNstc, hvar, hNvar);
282        hHilbEst(hstc, hNstc, hvar, hNvar);
283        memset(hpure, 0, (pVariables + 1) * sizeof(int));
284        hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
285        hLexS(hstc, hNstc, hvar, hNvar);
286        Q0[hNvar] = 0;
287        hHilbStep(hpure, hstc, hNstc, hvar, hNvar, p0, 1);
288      }
289    }
290    else
291    {
292      if(*Qpol!=NULL)
293        (**Qpol)++;
294      else
295      {
296        *Qpol = (int *)omAllocBin(int_bin);
297        hLength = *Ql = **Qpol = 1;
298      }
299    }
300    if (*Qpol!=NULL)
301    {
302      i = hLength;
303      while ((i > 0) && ((*Qpol)[i - 1] == 0))
304        i--;
305      if (i > 0)
306      {
307        l = i + j;
308        if (l > k)
309        {
310          work = new intvec(l);
311          for (ii=0; ii<k; ii++)
312            (*work)[ii] = (*hseries1)[ii];
313          if (hseries1 != NULL)
314            delete hseries1;
315          hseries1 = work;
316          k = l;
317        }
318        while (i > 0)
319        {
320          (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
321          (*Qpol)[i - 1] = 0;
322          i--;
323        }
324      }
325    }
326    mc--;
327    if (mc <= 0)
328      break;
329  }
330  if (k==0)
331  {
332    hseries1=new intvec(2);
333    (*hseries1)[0]=0;
334    (*hseries1)[1]=0;
335  }
336  else
337  {
338    l = k+1;
339    while ((*hseries1)[l-2]==0) l--;
340    if (l!=k)
341    {
342      work = new intvec(l);
343      for (ii=l-2; ii>=0; ii--)
344        (*work)[ii] = (*hseries1)[ii];
345      delete hseries1;
346      hseries1 = work;
347    }
348    (*hseries1)[l-1] = mw;
349  }
350  for (i = 0; i <= pVariables; i++)
351  {
352    if (Ql[i]!=0)
353      omFreeSize((ADDRESS)Qpol[i], Ql[i] * sizeof(int));
354  }
355  omFreeSize((ADDRESS)Q0, (pVariables + 1) * sizeof(int));
356  omFreeSize((ADDRESS)Ql, (pVariables + 1) * sizeof(int));
357  omFreeSize((ADDRESS)Qpol, (pVariables + 1) * sizeof(int *));
358  hKill(stcmem, pVariables - 1);
359  omFreeSize((ADDRESS)hpure, (1 + (pVariables * pVariables)) * sizeof(int));
360  omFreeSize((ADDRESS)hvar, (pVariables + 1) * sizeof(int));
361  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
362  hDelete(hexist, hNexist);
363  omFreeBin((ADDRESS)p0,  int_bin);
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(pVariables+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.