source: git/dyn_modules/callgfanlib/initialReduction.cc @ 528fa78

jengelh-datetimespielwiese
Last change on this file since 528fa78 was 528fa78, checked in by Yue Ren <ren@…>, 9 years ago
new: reduceInitially(ideal I, const number p)
  • Property mode set to 100644
File size: 6.7 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 * assumes that I is generated by elements which are homogeneous in x of the same degree.
199 **/
200bool reduceInitially(ideal I, const number p)
201{
202  sortByLeadmonom(I); int i,j;
203  for (i=1; i<IDELEMS(I); i++)
204    if (pReduce(I->m[i],p)) return true;
205
206  /***
207   * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
208   **/
209  for (i=0; i<IDELEMS(I)-1; i++)
210    for (j=i+1; j<IDELEMS(I); j++)
211      if (reduceInitially(I->m[j], I->m[i]) && pReduce(I->m[j],p)) return true;
212
213  /***
214   * the second pass. removing terms divisible by lt(g_j) out of g_i for i<j
215   **/
216  for (i=0; i<IDELEMS(I)-1; i++)
217    for (j=i+1; j<IDELEMS(I); j++)
218      if (reduceInitially(I->m[i], I->m[j]) && pReduce(I->m[i],p)) return true;
219  return false;
220}
221
222
223#ifndef NDEBUG
224BOOLEAN reduceInitially1(leftv res, leftv args)
225{
226  leftv u = args;
227  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
228  {
229    leftv v = u->next;
230    if ((v != NULL) && (v->Typ() == NUMBER_CMD))
231    {
232      ideal I = (ideal) u->CopyD();
233      number p = (number) v->CopyD();
234      omUpdateInfo();
235      Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
236      (void) reduceInitially(I,p);
237      omUpdateInfo();
238      Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
239      n_Delete(&p,currRing->cf);
240      res->rtyp = IDEAL_CMD;
241      res->data = (char*) I;
242      return FALSE;
243    }
244  }
245  return TRUE;
246}
247#endif //NDEBUG
Note: See TracBrowser for help on using the repository browser.