# source:git/factory/facFqSquarefree.cc@803956

jengelh-datetimespielwiese
Last change on this file since 803956 was 803956, checked in by Martin Lee <martinlee84@…>, 10 years ago
fix: wrong computation of pth root due to int overflow
• Property mode set to 100644
File size: 6.5 KB
Line
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facFqSquarefree.cc
5 *
6 * This file provides functions for squarefrees factorizing over
7 * \f$F_{p} \f$ , \f$F_{p}(\alpha ) \f$ or GF, as decribed in "Factoring
8 * multivariate polynomials over a finite field" by L. Bernardin
9 *
10 * @author Martin Lee
11 *
12 **/
13/*****************************************************************************/
14
15#include "config.h"
16
17#include "canonicalform.h"
18
19#include "cf_gcd_smallp.h"
20#include "cf_iter.h"
21#include "cf_map.h"
22#include "cf_util.h"
23#include "cf_factory.h"
24#include "facFqSquarefree.h"
25
26#ifdef HAVE_NTL
27#include "NTLconvert.h"
28#endif
29
30static inline
31CanonicalForm
32pthRoot (const CanonicalForm & F, const int & q)
33{
34  CanonicalForm A= F;
35  int p= getCharacteristic ();
36  if (A.inCoeffDomain())
37  {
38    A= power (A, q/p);
39    return A;
40  }
41  else
42  {
43    CanonicalForm buf= 0;
44    for (CFIterator i= A; i.hasTerms(); i++)
45      buf= buf + power(A.mvar(), i.exp()/p)*pthRoot (i.coeff(), q);
46    return buf;
47  }
48}
49
50#ifdef HAVE_NTL
51CanonicalForm
52pthRoot (const CanonicalForm & F, const ZZ& q, const Variable& alpha)
53{
54  CanonicalForm A= F;
55  int p= getCharacteristic ();
56  if (A.inCoeffDomain())
57  {
58    zz_p::init (p);
59    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
60    zz_pE::init (NTLMipo);
61    zz_pX NTLA= convertFacCF2NTLzzpX (A);
62    zz_pE NTLA2= to_zz_pE (NTLA);
63    power (NTLA2, NTLA2, q/p);
64    A= convertNTLzzpE2CF (NTLA2, alpha);
65    return A;
66  }
67  else
68  {
69    CanonicalForm buf= 0;
70    for (CFIterator i= A; i.hasTerms(); i++)
71      buf= buf + power(A.mvar(), i.exp()/p)*pthRoot (i.coeff(), q, alpha);
72    return buf;
73  }
74}
75#endif
76
77CanonicalForm
78maxpthRoot (const CanonicalForm & F, const int & q, int& l)
79{
80  CanonicalForm result= F;
81  bool derivZero= true;
82  l= 0;
83  while (derivZero)
84  {
85    for (int i= 1; i <= result.level(); i++)
86    {
87      if (!deriv (result, Variable (i)).isZero())
88      {
89        derivZero= false;
90        break;
91      }
92    }
93    if (!derivZero)
94      break;
95    result= pthRoot (result, q);
96    l++;
97  }
98  return result;
99}
100
101static inline
102CFFList
103sqrfPosDer (const CanonicalForm & F, const Variable & x,
104            CanonicalForm & c)
105{
106  CanonicalForm b= deriv (F, x);
107  c= gcd (F, b);
108  CanonicalForm w= F/c;
109  CanonicalForm v= b/c;
110  CanonicalForm u= v - deriv (w, x);
111  int j= 1;
112  int p= getCharacteristic();
113  CanonicalForm g;
114  CFFList result;
115  while (j < p - 1 && degree(u) >= 0)
116  {
117    g= gcd (w, u);
118    if (!g.inCoeffDomain())
119      result.append (CFFactor (g, j));
120    w= w/g;
121    c= c/w;
122    v= u/g;
123    u= v - deriv (w, x);
124    j++;
125  }
126  if (!w.inCoeffDomain())
127    result.append (CFFactor (w, j));
128  return result;
129}
130
131CFFList
132squarefreeFactorization (const CanonicalForm & F, const Variable & alpha)
133{
134  int p= getCharacteristic();
135  CanonicalForm A= F;
136  CFMap M;
137  A= compress (A, M);
138  Variable x= A.mvar();
139  int l= x.level();
140  int k;
141  if (CFFactory::gettype() == GaloisFieldDomain)
142    k= getGFDegree();
143  else if (alpha.level() != 1)
144    k= degree (getMipo (alpha));
145  else
146    k= 1;
147  Variable buf, buf2;
148  CanonicalForm tmp;
149
150  CFFList tmp1, tmp2;
151  bool found;
152  for (int i= l; i > 0; i--)
153  {
154    buf= Variable (i);
155    if (degree (deriv (A, buf)) >= 0)
156    {
157      tmp1= sqrfPosDer (A, buf, tmp);
158      A= tmp;
159      for (CFFListIterator j= tmp1; j.hasItem(); j++)
160      {
161        found= false;
162        CFFListIterator k= tmp2;
163        if (!k.hasItem() && !j.getItem().factor().inCoeffDomain()) tmp2.append (j.getItem());
164        else
165        {
166          for (; k.hasItem(); k++)
167          {
168            if (k.getItem().exp() == j.getItem().exp())
169            {
170              k.getItem()= CFFactor (k.getItem().factor()*j.getItem().factor(),
171                                     j.getItem().exp());
172              found= true;
173            }
174          }
175          if (found == false && !j.getItem().factor().inCoeffDomain())
176            tmp2.append(j.getItem());
177        }
178      }
179    }
180  }
181
182  bool degcheck= false;;
183  for (int i= l; i > 0; i--)
184  if (degree (A, Variable(i)) >= p)
185    degcheck= true;
186
187  if (degcheck == false && tmp1.isEmpty() && tmp2.isEmpty())
188    return CFFList (CFFactor (F/Lc(F), 1));
189
190  CanonicalForm buffer;
191#ifdef HAVE_NTL
192  if (alpha.level() == 1)
193#endif
194    buffer= pthRoot (A, ipower (p, k));
195#ifdef HAVE_NTL
196  else
197  {
198    ZZ q;
199    power (q, p, k);
200    buffer= pthRoot (A, q, alpha);
201  }
202#endif
203
204  tmp1= squarefreeFactorization (buffer, alpha);
205
206  CFFList result;
207  buf= alpha;
208  for (CFFListIterator i= tmp2; i.hasItem(); i++)
209  {
210    for (CFFListIterator j= tmp1; j.hasItem(); j++)
211    {
212      tmp= gcd (i.getItem().factor(), j.getItem().factor());
213      i.getItem()= CFFactor (i.getItem().factor()/tmp, i.getItem().exp());
214      j.getItem()= CFFactor (j.getItem().factor()/tmp, j.getItem().exp());
215      if (!tmp.inCoeffDomain())
216      {
217        tmp= M (tmp);
218        result.append (CFFactor (tmp/Lc(tmp),
219                       j.getItem().exp()*p + i.getItem().exp()));
220      }
221    }
222  }
223  for (CFFListIterator i= tmp2; i.hasItem(); i++)
224  {
225    if (!i.getItem().factor().inCoeffDomain())
226    {
227      tmp= M (i.getItem().factor());
228      result.append (CFFactor (tmp/Lc(tmp), i.getItem().exp()));
229    }
230  }
231  for (CFFListIterator j= tmp1; j.hasItem(); j++)
232  {
233    if (!j.getItem().factor().inCoeffDomain())
234    {
235      tmp= M (j.getItem().factor());
236      result.append (CFFactor (tmp/Lc(tmp), j.getItem().exp()*p));
237    }
238  }
239  return result;
240}
241
242CanonicalForm
243sqrfPart (const CanonicalForm& F, CanonicalForm& pthPower,
244          const Variable& alpha)
245{
246  if (F.inCoeffDomain())
247  {
248    pthPower= 1;
249    return F;
250  }
251  CFMap M;
252  CanonicalForm A= compress (F, M);
253  Variable vBuf= alpha;
254  CanonicalForm w, v, b;
255  pthPower= 1;
256  CanonicalForm result;
257  int i= 1;
258  bool allZero= true;
259  for (; i <= A.level(); i++)
260  {
261    if (!deriv (A, Variable (i)).isZero())
262    {
263      allZero= false;
264      break;
265    }
266  }
267  if (allZero)
268  {
269    pthPower= F;
270    return 1;
271  }
272  w= gcd (A, deriv (A, Variable (i)));
273
274  b= A/w;
275  result= b;
276  if (degree (w) < 1)
277    return M (result);
278  i++;
279  for (; i <= A.level(); i++)
280  {
281    if (!deriv (w, Variable (i)).isZero())
282    {
283      b= w;
284      w= gcd (w, deriv (w, Variable (i)));
285      b /= w;
286      if (degree (b) < 1)
287        break;
288      CanonicalForm g;
289      g= gcd (b, result);
290      if (degree (g) > 0)
291        result *= b/g;
292      if (degree (g) <= 0)
293        result *= b;
294    }
295  }
296  result= M (result);
297  return result;
298}
299
Note: See TracBrowser for help on using the repository browser.