source: git/factory/facBivar.cc @ d990001

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