source: git/dyn_modules/callgfanlib/initialReduction.cc @ edb60f

jengelh-datetimespielwiese
Last change on this file since edb60f was edb60f, checked in by Yue Ren <ren@…>, 9 years ago
new: reduceInitially(ideal I, const number p, poly g)
  • Property mode set to 100644
File size: 9.4 KB
Line 
1#include <kernel/polys.h>
2#include <libpolys/polys/monomials/p_polys.h>
3#include <singularWishlist.h>
4
5#include <Singular/ipid.h>
6
7/***
8 * changes a polynomial g with the help p-t such that
9 * 1) each term of g has a distinct monomial in x
10 * 2) no term of g has a coefficient divisible by p
11 * in particular, this means that all g_\alpha can be obtained
12 * by reading the coefficients and that g is initially reduced
13 * with respect to p-t
14 **/
15static bool pReduce(poly g, const number p)
16{
17  poly toBeChecked = pNext(g);
18  pNext(g) = NULL; poly gEnd = g;
19  poly gCache;
20
21  number coeff, pPower; int power; poly subst;
22  while(toBeChecked)
23  {
24    for (gCache = g; gCache; pIter(gCache))
25      if (p_LeadmonomDivisibleBy(gCache,toBeChecked,currRing)) break;
26    if (gCache)
27    {
28      n_Power(p,p_GetExp(toBeChecked,1,currRing)-p_GetExp(gCache,1,currRing),&pPower,currRing->cf);
29      coeff = n_Mult(p_GetCoeff(toBeChecked,currRing),pPower,currRing->cf);
30      p_SetCoeff(gCache,n_Add(p_GetCoeff(gCache,currRing),coeff,currRing->cf),currRing);
31      n_Delete(&pPower,currRing->cf); n_Delete(&coeff,currRing->cf);
32      toBeChecked=p_LmDeleteAndNext(toBeChecked,currRing);
33    }
34    else
35    {
36      if (n_DivBy(p_GetCoeff(toBeChecked,currRing),p,currRing->cf))
37      {
38        coeff=n_Div(p_GetCoeff(toBeChecked,currRing),p,currRing->cf);
39        power=1;
40        while (n_DivBy(coeff,p,currRing->cf))
41        {
42          coeff=n_Div(p_GetCoeff(pNext(g),currRing),p,currRing->cf);
43          power++;
44          if (power<1)
45          {
46            WerrorS("pReduce: overflow in exponent");
47            return true;
48          }
49        }
50        subst=p_LmInit(toBeChecked,currRing);
51        p_AddExp(subst,1,power,currRing);
52        p_SetCoeff(subst,coeff,currRing);
53        p_Setm(subst,currRing); pTest(subst);
54        toBeChecked=p_LmDeleteAndNext(toBeChecked,currRing);
55        toBeChecked=p_Add_q(toBeChecked,subst,currRing);
56        pTest(toBeChecked);
57      }
58      else
59      {
60        pNext(gEnd)=toBeChecked;
61        pIter(gEnd); pIter(toBeChecked);
62        pNext(gEnd)=NULL;
63      }
64    }
65  }
66  return false;
67}
68
69
70#ifndef NDEBUG
71BOOLEAN pReduce(leftv res, leftv args)
72{
73  leftv u = args;
74  if ((u != NULL) && (u->Typ() == POLY_CMD))
75  {
76    poly g; number p = n_Init(3,currRing->cf);
77    omUpdateInfo();
78    Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
79    g = (poly) u->CopyD();
80    (void) pReduce(g,p);
81    p_Delete(&g,currRing);
82    omUpdateInfo();
83    Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
84    g = (poly) u->CopyD();
85    (void) pReduce(g,p);
86    n_Delete(&p,currRing->cf);
87    res->rtyp = POLY_CMD;
88    res->data = (char*) g;
89    return FALSE;
90  }
91  return TRUE;
92}
93#endif //NDEBUG
94
95
96/***
97 * returns, if it exists, a pointer to the first term in g
98 * that is divisible by (the leading monomial of) m or, if it does not exist, the NULL pointer.
99 * if g is homogeneous in x with the same degree as m,
100 * then it returns the first term with the same monomial in x as m,
101 * if the t-degree of the term is higher than the t-degree of m, or NULL otherwise.
102 **/
103static inline poly firstTermDivisibleBy(const poly g, const poly m)
104{
105  poly gCache = NULL;
106  for (gCache=g; gCache; pIter(gCache))
107    if (p_LeadmonomDivisibleBy(m,gCache,currRing)) break;
108  return gCache;
109}
110
111
112/***
113 * reduces h initially with respect to g,
114 * returns false if h was initially reduced in the first place,
115 * returns true if reductions have taken place.
116 * assumes that h and g are in pReduced form and homogeneous in x of the same degree
117 **/
118bool reduceInitially(poly &h, const poly g)
119{
120  poly hCache;
121  for (hCache=h; hCache; pIter(hCache))
122    if (p_LeadmonomDivisibleBy(g,hCache,currRing)) break;
123  if (hCache)
124  {
125    number gAlpha = p_GetCoeff(g,currRing);
126    poly hAlphaT = p_Init(currRing);
127    p_SetCoeff(hAlphaT,n_Copy(p_GetCoeff(hCache,currRing),currRing->cf),currRing);
128    p_SetExp(hAlphaT,1,p_GetExp(hCache,1,currRing)-p_GetExp(g,1,currRing),currRing);
129    for (int i=2; i<=currRing->N; i++)
130      p_SetExp(hAlphaT,i,0,currRing);
131    p_Setm(hAlphaT,currRing); pTest(hAlphaT);
132    h = p_Add_q(p_Mult_nn(h,gAlpha,currRing),
133                p_Neg(p_Mult_q(p_Copy(g,currRing),hAlphaT,currRing),currRing),
134                currRing);
135    pTest(h);
136    return true;
137  }
138  return false;
139}
140
141
142#ifndef NDEBUG
143BOOLEAN reduceInitially0(leftv res, leftv args)
144{
145  leftv u = args;
146  if ((u != NULL) && (u->Typ() == POLY_CMD))
147  {
148    leftv v = u->next;
149    if ((v != NULL) && (v->Typ() == POLY_CMD))
150    {
151      poly g,h;
152      omUpdateInfo();
153      Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
154      h = (poly) u->CopyD();
155      g = (poly) v->CopyD();
156      (void)reduceInitially(h,g);
157      p_Delete(&h,currRing);
158      p_Delete(&g,currRing);
159      omUpdateInfo();
160      Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
161      h = (poly) u->CopyD();
162      g = (poly) v->CopyD();
163      (void)reduceInitially(h,g);
164      p_Delete(&g,currRing);
165      res->rtyp = POLY_CMD;
166      res->data = (char*) h;
167      return FALSE;
168    }
169  }
170  return TRUE;
171}
172#endif //NDEBUG
173
174
175static inline void sortByLeadmonom(ideal I)
176{
177  poly cache; int i, n=IDELEMS(I), newn;
178  do
179  {
180    newn=0;
181    for (i=1; i<n; i++)
182    {
183      if (pLmCmp(I->m[i-1],I->m[i])<0)
184      {
185        cache=I->m[i-1];
186        I->m[i-1]=I->m[i];
187        I->m[i]=cache;
188        newn = i;
189      }
190    }
191    n=newn;
192  } while(n);
193}
194
195
196/***
197 * reduces I initially with respect to itself and with respect to p-t.
198 * also sorts the generators of I with respect to the leading monomials in descending order.
199 * assumes that I is generated by elements which are homogeneous in x of the same degree.
200 **/
201bool reduceInitially(ideal I, const number p)
202{
203  poly cache; int i,j,n=IDELEMS(I);
204  do
205  {
206    j=0;
207    for (i=1; i<n; i++)
208    {
209      if (pLmCmp(I->m[i-1],I->m[i])<0)
210      {
211        cache=I->m[i-1];
212        I->m[i-1]=I->m[i];
213        I->m[i]=cache;
214        j = i;
215      }
216    }
217    n=j;
218  } while(n);
219  for (i=1; i<IDELEMS(I); i++)
220    if (pReduce(I->m[i],p)) return true;
221
222  /***
223   * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
224   **/
225  for (i=0; i<IDELEMS(I)-1; i++)
226    for (j=i+1; j<IDELEMS(I); j++)
227      if (reduceInitially(I->m[j], I->m[i]) && pReduce(I->m[j],p)) return true;
228
229  /***
230   * the second pass. removing terms divisible by lt(g_j) out of g_i for i<j
231   **/
232  for (i=0; i<IDELEMS(I)-1; i++)
233    for (j=i+1; j<IDELEMS(I); j++)
234      if (reduceInitially(I->m[i], I->m[j]) && pReduce(I->m[i],p)) return true;
235  return false;
236}
237
238
239#ifndef NDEBUG
240BOOLEAN reduceInitially1(leftv res, leftv args)
241{
242  leftv u = args;
243  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
244  {
245    leftv v = u->next;
246    if ((v != NULL) && (v->Typ() == NUMBER_CMD))
247    {
248      ideal I; number p;
249      omUpdateInfo();
250      Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
251      I = (ideal) u->CopyD();
252      p = (number) v->CopyD();
253      (void) reduceInitially(I,p);
254      id_Delete(&I,currRing);
255      n_Delete(&p,currRing->cf);
256      omUpdateInfo();
257      Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
258      I = (ideal) u->CopyD();
259      p = (number) v->CopyD();
260      (void) reduceInitially(I,p);
261      n_Delete(&p,currRing->cf);
262      res->rtyp = IDEAL_CMD;
263      res->data = (char*) I;
264      return FALSE;
265    }
266  }
267  return TRUE;
268}
269#endif //NDEBUG
270
271
272/***
273 * inserts g into I and reduces I with respect to itself and p-t
274 * assumes that I was already sorted and initially reduced in the first place
275 **/
276bool reduceInitially(ideal I, const number p, poly g)
277{
278  pEnlargeSet(&(I->m),IDELEMS(I),1);
279  IDELEMS(I)++; int i,j,k;
280  for (j=IDELEMS(I)-1; j>0; j--)
281  {
282    if (pLmCmp(I->m[j-1],g)<0)
283      I->m[j] = I->m[j-1];
284    else
285    {
286      I->m[j] = g;
287      break;
288    }
289  }
290  if (j<1) I->m[0] = g;
291
292  /***
293   * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
294   * removing terms with the same monomials in x as lt(g_j) out of g_k for j<k
295   **/
296  for (i=0; i<j; i++)
297    if (reduceInitially(I->m[j], I->m[i]) && pReduce(I->m[j],p)) return true;
298  for (k=j+1; k<IDELEMS(I); k++)
299    if (reduceInitially(I->m[k], I->m[j]) && pReduce(I->m[k],p)) return true;
300
301  /***
302   * the second pass. removing terms divisible by lt(g_j) and lt(g_k) out of g_i for i<j<k
303   * removing terms divisible by lt(g_k) out of g_j for j<k
304   **/
305  for (i=0; i<j; i++)
306    for (k=j; k<IDELEMS(I); k++)
307      if (reduceInitially(I->m[i], I->m[k]) && pReduce(I->m[i],p)) return true;
308  for (k=j+1; k<IDELEMS(I); k++)
309    if (reduceInitially(I->m[j],I->m[k]) && pReduce(I->m[j],p)) return true;
310
311  return false;
312}
313
314
315#ifndef NDEBUG
316BOOLEAN reduceInitially2(leftv res, leftv args)
317{
318  leftv u = args;
319  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
320  {
321    leftv v = u->next;
322    if ((v != NULL) && (v->Typ() == NUMBER_CMD))
323    {
324      leftv w = v->next;
325      if ((w != NULL) && (w->Typ() == POLY_CMD))
326      {
327        ideal I; number p; poly g;
328        omUpdateInfo();
329        Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
330        I = (ideal) u->CopyD();
331        p = (number) v->CopyD();
332        g = (poly) w->CopyD();
333        (void) reduceInitially(I,p,g);
334        id_Delete(&I,currRing);
335        n_Delete(&p,currRing->cf);
336        omUpdateInfo();
337        Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
338        I = (ideal) u->CopyD();
339        p = (number) v->CopyD();
340        g = (poly) w->CopyD();
341        (void) reduceInitially(I,p,g);
342        n_Delete(&p,currRing->cf);
343        res->rtyp = IDEAL_CMD;
344        res->data = (char*) I;
345        return FALSE;
346      }
347    }
348  }
349  return TRUE;
350}
351#endif //NDEBUG
Note: See TracBrowser for help on using the repository browser.