source: git/factory/facBivar.cc @ e3198f

spielwiese
Last change on this file since e3198f was 1a3011e, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: avoid double checking of factors during henselLiftAndEarly chg: lower liftBound in bivariate factorization
  • Property mode set to 100644
File size: 9.9 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facBivar.cc
5 *
6 * bivariate factorization over Q(a)
7 *
8 * @author Martin Lee
9 *
10 * @internal @version \$Id$
11 *
12 **/
13/*****************************************************************************/
14
15#include "config.h"
16
17#include "assert.h"
18#include "debug.h"
19#include "timing.h"
20
21#include "cf_algorithm.h"
22#include "facFqBivar.h"
23#include "facBivar.h"
24
25#ifdef HAVE_NTL
26TIMING_DEFINE_PRINT(fac_uni_factorizer)
27TIMING_DEFINE_PRINT(fac_hensel_lift)
28TIMING_DEFINE_PRINT(fac_factor_recombination)
29
30CFList conv (const CFFList& L)
31{
32  CFList result;
33  for (CFFListIterator i= L; i.hasItem(); i++)
34    result.append (i.getItem().factor());
35  return result;
36}
37
38bool testPoint (const CanonicalForm& F, CanonicalForm& G, int i)
39{
40  G= F (i, 2);
41  if (G.inCoeffDomain() || degree (F, 1) > degree (G, 1))
42    return false;
43
44  if (degree (gcd (deriv (G, G.mvar()), G)) > 0)
45    return false;
46  return true;
47}
48
49CanonicalForm evalPoint (const CanonicalForm& F, int& i)
50{
51  Variable x= Variable (1);
52  Variable y= Variable (2);
53  CanonicalForm result;
54
55  int k;
56
57  if (i == 0)
58  {
59    if (testPoint (F, result, i))
60      return result;
61  }
62  do
63  {
64    if (i > 0)
65      k= 1;
66    else
67      k= 2;
68    while (k < 3)
69    {
70      if (k == 1)
71      {
72        if (testPoint (F, result, i))
73          return result;
74      }
75      else
76      {
77        if (testPoint (F, result, -i))
78          return result;
79      }
80      k++;
81    }
82    i++;
83  } while (1);
84}
85
86CFList biFactorize (const CanonicalForm& F, const Variable& v)
87{
88  if (F.inCoeffDomain())
89    return CFList(F);
90
91  bool extension= (v.level() != 1);
92  CanonicalForm A;
93  CanonicalForm multiplier= 1;
94  if (isOn (SW_RATIONAL))
95    A= F*bCommonDen (F);
96  else
97    A= F;
98
99  if (A.isUnivariate())
100  {
101    CFFList buf;
102    if (extension)
103      buf= factorize (A, v);
104    else
105      buf= factorize (A, true);
106    CFList result= conv (buf);
107    if (result.getFirst().inCoeffDomain())
108      result.removeFirst();
109    return result;
110  }
111
112  CFMap N;
113  A= compress (A, N);
114  Variable y= A.mvar();
115
116  if (y.level() > 2) return CFList (F);
117  Variable x= Variable (1);
118
119  CanonicalForm contentAx= content (A, x);
120  CanonicalForm contentAy= content (A);
121
122  A= A/(contentAx*contentAy);
123  CFFList contentAxFactors, contentAyFactors;
124
125  if (extension)
126  {
127    if (!contentAx.inCoeffDomain())
128    {
129      contentAxFactors= factorize (contentAx, v);
130      if (contentAxFactors.getFirst().factor().inCoeffDomain())
131        contentAxFactors.removeFirst();
132    }
133    if (!contentAy.inCoeffDomain())
134    {
135      contentAyFactors= factorize (contentAy, v);
136      if (contentAyFactors.getFirst().factor().inCoeffDomain())
137        contentAyFactors.removeFirst();
138    }
139  }
140  else
141  {
142    if (!contentAx.inCoeffDomain())
143    {
144      contentAxFactors= factorize (contentAx, true);
145      if (contentAxFactors.getFirst().factor().inCoeffDomain())
146        contentAxFactors.removeFirst();
147    }
148    if (!contentAy.inCoeffDomain())
149    {
150      contentAyFactors= factorize (contentAy, true);
151      if (contentAyFactors.getFirst().factor().inCoeffDomain())
152        contentAyFactors.removeFirst();
153    }
154  }
155
156  //check trivial case
157  if (degree (A) == 1 || degree (A, 1) == 1)
158  {
159    CFList factors;
160    factors.append (A);
161
162    appendSwapDecompress (factors, conv (contentAxFactors),
163                          conv (contentAyFactors), false, false, N);
164
165    normalize (factors);
166    return factors;
167  }
168
169  //trivial case
170  CFList factors;
171  if (A.inCoeffDomain())
172  {
173    append (factors, conv (contentAxFactors));
174    append (factors, conv (contentAyFactors));
175    decompress (factors, N);
176    return factors;
177  }
178  else if (A.isUnivariate())
179  {
180    if (extension)
181      factors= conv (factorize (A, v));
182    else
183      factors= conv (factorize (A, true));
184    append (factors, conv (contentAxFactors));
185    append (factors, conv (contentAyFactors));
186    decompress (factors, N);
187    return factors;
188  }
189
190  bool swap= false;
191  if (degree (A) > degree (A, x))
192  {
193    A= swapvar (A, y, x);
194    swap= true;
195  }
196
197  CanonicalForm Aeval, bufAeval, buf;
198  CFList uniFactors, list, bufUniFactors;
199  DegreePattern degs;
200  DegreePattern bufDegs;
201
202  CanonicalForm Aeval2, bufAeval2;
203  CFList bufUniFactors2, list2, uniFactors2;
204  DegreePattern degs2;
205  DegreePattern bufDegs2;
206  bool swap2= false;
207
208  // several univariate factorizations to obtain more information about the
209  // degree pattern therefore usually less combinations have to be tried during
210  // the recombination process
211  int factorNums= 2;
212  int subCheck1= substituteCheck (A, x);
213  int subCheck2= substituteCheck (A, y);
214  buf= swapvar (A,x,y);
215  int evaluation, evaluation2, bufEvaluation= 0, bufEvaluation2= 0;
216  for (int i= 0; i < factorNums; i++)
217  {
218    bufAeval= A;
219    bufAeval= evalPoint (A, bufEvaluation);
220
221    bufAeval2= buf;
222    bufAeval2= evalPoint (buf, bufEvaluation2);
223
224
225    // univariate factorization
226    TIMING_START (fac_uni_factorizer);
227
228    if (extension)
229      bufUniFactors= conv (factorize (bufAeval, v));
230    else
231      bufUniFactors= conv (factorize (bufAeval, true));
232    TIMING_END_AND_PRINT (fac_uni_factorizer,
233                          "time for univariate factorization: ");
234    DEBOUTLN (cerr, "Lc (bufAeval)*prod (bufUniFactors)== bufAeval " <<
235              prod (bufUniFactors)*Lc (bufAeval) == bufAeval);
236
237    TIMING_START (fac_uni_factorizer);
238    if (extension)
239      bufUniFactors2= conv (factorize (bufAeval2, v));
240    else
241      bufUniFactors2= conv (factorize (bufAeval2, true));
242    TIMING_END_AND_PRINT (fac_uni_factorizer,
243                          "time for univariate factorization in y: ");
244    DEBOUTLN (cerr, "Lc (Aeval2)*prod (uniFactors2)== Aeval2 " <<
245              prod (bufUniFactors2)*Lc (bufAeval2) == bufAeval2);
246
247    if (bufUniFactors.getFirst().inCoeffDomain())
248      bufUniFactors.removeFirst();
249    if (bufUniFactors2.getFirst().inCoeffDomain())
250      bufUniFactors2.removeFirst();
251    if (bufUniFactors.length() == 1 || bufUniFactors2.length() == 1)
252    {
253      factors.append (A);
254
255      appendSwapDecompress (factors, conv (contentAxFactors),
256                            conv (contentAyFactors), swap, swap2, N);
257
258      if (isOn (SW_RATIONAL))
259        normalize (factors);
260      return factors;
261    }
262
263    if (i == 0)
264    {
265      if (subCheck1 > 0)
266      {
267        int subCheck= substituteCheck (bufUniFactors);
268
269        if (subCheck > 1 && (subCheck1%subCheck == 0))
270        {
271          CanonicalForm bufA= A;
272          subst (bufA, bufA, subCheck, x);
273          factors= biFactorize (bufA, v);
274          reverseSubst (factors, subCheck, x);
275          appendSwapDecompress (factors, conv (contentAxFactors),
276                                conv (contentAyFactors), swap, swap2, N);
277          if (isOn (SW_RATIONAL))
278            normalize (factors);
279          return factors;
280        }
281      }
282
283      if (subCheck2 > 0)
284      {
285        int subCheck= substituteCheck (bufUniFactors2);
286
287        if (subCheck > 1 && (subCheck2%subCheck == 0))
288        {
289          CanonicalForm bufA= A;
290          subst (bufA, bufA, subCheck, y);
291          factors= biFactorize (bufA, v);
292          reverseSubst (factors, subCheck, y);
293          appendSwapDecompress (factors, conv (contentAxFactors),
294                                conv (contentAyFactors), swap, swap2, N);
295          if (isOn (SW_RATIONAL))
296            normalize (factors);
297          return factors;
298        }
299      }
300    }
301
302    // degree analysis
303    bufDegs = DegreePattern (bufUniFactors);
304    bufDegs2= DegreePattern (bufUniFactors2);
305
306    if (i == 0)
307    {
308      Aeval= bufAeval;
309      evaluation= bufEvaluation;
310      uniFactors= bufUniFactors;
311      degs= bufDegs;
312      Aeval2= bufAeval2;
313      evaluation2= bufEvaluation2;
314      uniFactors2= bufUniFactors2;
315      degs2= bufDegs2;
316    }
317    else
318    {
319      degs.intersect (bufDegs);
320      degs2.intersect (bufDegs2);
321      if (bufUniFactors2.length() < uniFactors2.length())
322      {
323        uniFactors2= bufUniFactors2;
324        Aeval2= bufAeval2;
325        evaluation2= bufEvaluation2;
326      }
327      if (bufUniFactors.length() < uniFactors.length())
328      {
329        if (evaluation != 0)
330        {
331          uniFactors= bufUniFactors;
332          Aeval= bufAeval;
333          evaluation= bufEvaluation;
334        }
335      }
336    }
337    bufEvaluation++;
338    bufEvaluation2++;
339  }
340
341  if (evaluation != 0 && (uniFactors.length() > uniFactors2.length() ||
342      (uniFactors.length() == uniFactors2.length()
343       && degs.getLength() > degs2.getLength())))
344  {
345    degs= degs2;
346    uniFactors= uniFactors2;
347    evaluation= evaluation2;
348    Aeval= Aeval2;
349    A= buf;
350    swap2= true;
351  }
352
353
354  if (degs.getLength() == 1) // A is irreducible
355  {
356    factors.append (A);
357    appendSwapDecompress (factors, conv (contentAxFactors),
358                          conv (contentAyFactors), swap, swap2, N);
359    if (isOn (SW_RATIONAL))
360      normalize (factors);
361    return factors;
362  }
363
364  A= A (y + evaluation, y);
365
366  int liftBound= degree (A, y) + 1;
367
368  ExtensionInfo dummy= ExtensionInfo (false);
369  bool earlySuccess= false;
370  CFList earlyFactors;
371  TIMING_START (fac_hensel_lift);
372  uniFactors= henselLiftAndEarly
373             (A, earlySuccess, earlyFactors, degs, liftBound,
374              uniFactors, dummy, evaluation);
375  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
376  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
377
378  CanonicalForm MODl= power (y, liftBound);
379
380  factors= factorRecombination (uniFactors, A, MODl, degs, 1,
381                                uniFactors.length()/2);
382
383  if (earlySuccess)
384    factors= Union (earlyFactors, factors);
385  else if (!earlySuccess && degs.getLength() == 1)
386    factors= earlyFactors;
387
388  for (CFListIterator i= factors; i.hasItem(); i++)
389    i.getItem()= i.getItem() (y - evaluation, y);
390
391  appendSwapDecompress (factors, conv (contentAxFactors),
392                        conv (contentAyFactors), swap, swap2, N);
393  if (isOn (SW_RATIONAL))
394    normalize (factors);
395
396  return factors;
397}
398
399#endif
Note: See TracBrowser for help on using the repository browser.