source: git/factory/facBivar.cc @ e4fe2b

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