source: git/factory/facAlgExt.cc @ 9e74d36

spielwiese
Last change on this file since 9e74d36 was 9e74d36, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: break factor recovery earlier chg: make sure factor's not a constant from the alg. func. field
  • Property mode set to 100644
File size: 10.1 KB
Line 
1// -*- c++ -*-
2//*****************************************************************************
3/** @file facAlgExt.cc
4 *
5 * @author Martin Lee
6 * @date
7 *
8 * Univariate factorization over algebraic extension of Q using Trager's
9 * algorithm
10 *
11 * @par Copyright:
12 *   (c) by The SINGULAR Team, see LICENSE file
13 *
14**/
15//*****************************************************************************
16
17
18#include "config.h"
19
20
21#include "cf_assert.h"
22#include "debug.h"
23#include "timing.h"
24
25#include "canonicalform.h"
26#include "cf_random.h"
27#include "cf_algorithm.h"
28#include "facFqBivarUtil.h"
29#include "facAlgExt.h"
30#include "cfModResultant.h"
31#include "fac_sqrfree.h"
32
33TIMING_DEFINE_PRINT(fac_alg_resultant)
34TIMING_DEFINE_PRINT(fac_alg_norm)
35TIMING_DEFINE_PRINT(fac_alg_factor_norm)
36TIMING_DEFINE_PRINT(fac_alg_gcd)
37TIMING_DEFINE_PRINT(fac_alg_sqrf)
38TIMING_DEFINE_PRINT(fac_alg_factor_sqrf)
39TIMING_DEFINE_PRINT(fac_alg_time_shift)
40
41// squarefree part of F
42CanonicalForm
43uniSqrfPart (const CanonicalForm& F)
44{
45  ASSERT (F.isUnivariate(), "univariate input expected");
46  ASSERT (getCharacteristic() == 0, "characteristic 0 expected");
47  CanonicalForm G= deriv (F, F.mvar());
48  G= gcd (F, G);
49  return F/G;
50}
51
52CanonicalForm Norm (const CanonicalForm& F, const Variable& alpha)
53{
54  Variable x= Variable (F.level() + 1);
55  Variable y= F.mvar();
56  CanonicalForm g= F (x, alpha);
57  CanonicalForm mipo= getMipo (alpha);
58  mipo= mipo (x, alpha);
59  mipo *= bCommonDen (mipo);
60
61  int degg= degree (g);
62  int degmipo= degree (mipo);
63  CanonicalForm norm;
64  TIMING_START (fac_alg_resultant);
65  if (degg >= 8 || degmipo >= 8)
66    norm= resultantZ (g, mipo, x);
67  else
68    norm= resultant (g, mipo, x);
69  TIMING_END_AND_PRINT (fac_alg_resultant, "time to compute resultant0: ");
70  return norm;
71}
72
73// i is an integer such that Norm (F (x-i*alpha)) is squarefree
74CanonicalForm sqrfNorm (const CanonicalForm& F, const Variable& alpha, int& i)
75{
76  Variable x= Variable (F.level() + 1);
77  Variable y= F.mvar();
78  CanonicalForm g= F (x, alpha);
79  CanonicalForm mipo= getMipo (alpha);
80  mipo= mipo (x, alpha);
81  mipo *= bCommonDen (mipo);
82
83  int degg= degree (g);
84  int degmipo= degree (mipo);
85  CanonicalForm norm;
86  TIMING_START (fac_alg_resultant);
87  if (degg >= 8 || degmipo >= 8)
88    norm= resultantZ (g, mipo, x);
89  else
90    norm= resultant (g, mipo, x);
91  TIMING_END_AND_PRINT (fac_alg_resultant, "time to compute resultant0: ");
92
93  i= 0;
94  int k;
95  if (degree (gcd (deriv (norm, y), norm)) <= 0)
96    return norm;
97  i= 1;
98  do
99  {
100    k= 1;
101    while (k < 3)
102    {
103      if (k == 1)
104      {
105        g= F (y - i*alpha, y);
106        g *= bCommonDen (g);
107        TIMING_START (fac_alg_resultant);
108        if (degg >= 8 || degmipo >= 8)
109          norm= resultantZ (g (x, alpha), mipo, x);
110        else
111          norm= resultant (g (x, alpha), mipo, x);
112        TIMING_END_AND_PRINT (fac_alg_resultant,"time to compute resultant1: ");
113      }
114      else
115      {
116        g= F (y + i*alpha, y);
117        g *= bCommonDen (g);
118        TIMING_START (fac_alg_resultant);
119        if (degg >= 8 || degmipo >= 8)
120          norm= resultantZ (g (x, alpha), mipo, x);
121        else
122          norm= resultant (g (x, alpha), mipo, x);
123        TIMING_END_AND_PRINT (fac_alg_resultant,"time to compute resultant2: ");
124      }
125      if (degree (gcd (deriv (norm, y), norm)) <= 0)
126      {
127        if (k == 2)
128          i= -i;
129        return norm;
130      }
131      k++;
132    }
133    i++;
134  } while (1);
135}
136
137CFList
138AlgExtSqrfFactorize (const CanonicalForm& F, const Variable& alpha)
139{
140  ASSERT (F.isUnivariate(), "univariate input expected");
141  ASSERT (getCharacteristic() == 0, "characteristic 0 expected");
142
143  bool save_rat=!isOn (SW_RATIONAL);
144  On (SW_RATIONAL);
145  CanonicalForm f= F*bCommonDen (F);
146  Variable y= f.mvar();
147  int shift= 0, k= 0, count= 0;
148  CanonicalForm norm, buf, factor, oldF;
149  CFFList normFactors;
150  bool save_sort= !isOn (SW_USE_NTL_SORT);
151  CFList factors, tmp, tmp2;
152  CFFListIterator i;
153  CFListIterator iter;
154  bool shiftBuf= false;
155
156  tmp.append (f);
157  do
158  {
159    tmp2= CFList();
160    for (iter= tmp; iter.hasItem(); iter++)
161    {
162      oldF= iter.getItem()*bCommonDen (iter.getItem());
163      if (shift == 0)
164        f= oldF;
165      else
166        f= oldF (y - shift*alpha, y);
167      TIMING_START (fac_alg_norm);
168      norm= Norm (f, alpha);
169      TIMING_END_AND_PRINT (fac_alg_norm, "time to compute sqrf norm: ");
170      ASSERT (degree (norm, alpha) <= 0, "wrong norm computed");
171
172      TIMING_START (fac_alg_factor_norm);
173      On (SW_USE_NTL_SORT);
174      normFactors= factorize (norm);
175      if (save_sort)
176        Off (SW_USE_NTL_SORT);
177      TIMING_END_AND_PRINT (fac_alg_factor_norm, "time to factor norm: ");
178
179      if (normFactors.getFirst().factor().inCoeffDomain())
180        normFactors.removeFirst();
181      if (normFactors.length() < 2 && normFactors.getLast().exp() == 1)
182      {
183        factors.append (oldF);
184        continue;
185      }
186
187      i= normFactors;
188      shiftBuf= false;
189      if (!(normFactors.length() == 2 &&
190            degree (i.getItem().factor()) <= degree (f)))
191      {
192        TIMING_START (fac_alg_time_shift);
193        if (shift != 0)
194          buf= f;
195        else
196          buf= oldF;
197        shiftBuf= true;
198        TIMING_END_AND_PRINT (fac_alg_time_shift, "time to shift: ");
199      }
200      else
201        buf= oldF;
202
203      count= 0;
204      for (; i.hasItem(); i++)
205      {
206        TIMING_START (fac_alg_gcd);
207        if (shiftBuf)
208          factor= gcd (buf, i.getItem().factor());
209        else
210        {
211          if (shift == 0)
212            factor= gcd (buf, i.getItem().factor());
213          else
214            factor= gcd (buf, i.getItem().factor() (y + shift*alpha, y));
215        }
216        buf /= factor;
217        if (shiftBuf)
218        {
219          if (shift != 0)
220            factor= factor (y + shift*alpha, y);
221        }
222        TIMING_END_AND_PRINT (fac_alg_gcd, "time to recover factors: ");
223        if (i.getItem().exp() == 1 || degree (factor) == 1)
224          factors.append (factor);
225        else
226          tmp2.append (factor);
227        if (buf.inCoeffDomain())
228          break;
229        count++;
230        if (normFactors.length() - 1 == count)
231        {
232          if (shiftBuf)
233          {
234            if (normFactors.getLast().exp() == 1)
235              factors.append (buf (y + shift*alpha, y));
236            else
237              tmp2.append (buf (y + shift*alpha, y));
238          }
239          else
240          {
241            if (normFactors.getLast().exp() == 1)
242              factors.append (buf);
243            else
244              tmp2.append (buf);
245          }
246          buf= 1;
247          break;
248        }
249      }
250    }
251    k++;
252    if (shift == 0)
253    {
254      shift++;
255      k= 1;
256    }
257    if (k == 2)
258      shift= -shift;
259    if (k == 3)
260    {
261      shift= -shift;
262      shift++;
263      k= 1;
264    }
265    tmp= tmp2;
266  }
267  while (!tmp.isEmpty());
268
269  if (save_rat) Off(SW_RATIONAL);
270  return factors;
271}
272
273
274/*CFList
275AlgExtSqrfFactorize (const CanonicalForm& F, const Variable& alpha)
276{
277  ASSERT (F.isUnivariate(), "univariate input expected");
278  ASSERT (getCharacteristic() == 0, "characteristic 0 expected");
279
280  bool save_rat=!isOn (SW_RATIONAL);
281  On (SW_RATIONAL);
282  CanonicalForm f= F*bCommonDen (F);
283  int shift;
284  TIMING_START (fac_alg_norm);
285  CanonicalForm norm= sqrfNorm (f, alpha, shift);
286  TIMING_END_AND_PRINT (fac_alg_norm, "time to compute sqrf norm: ");
287  ASSERT (degree (norm, alpha) <= 0, "wrong norm computed");
288  TIMING_START (fac_alg_factor_norm);
289  bool save_sort= !isOn (SW_USE_NTL_SORT);
290  On (SW_USE_NTL_SORT);
291  CFFList normFactors= factorize (norm);
292  if (save_sort)
293    Off (SW_USE_NTL_SORT);
294  TIMING_END_AND_PRINT (fac_alg_factor_norm, "time to factor norm: ");
295  CFList factors;
296  if (normFactors.length() <= 2)
297  {
298    if (save_rat) Off(SW_RATIONAL);
299    return CFList (F);
300  }
301
302  normFactors.removeFirst();
303  CFFListIterator i= normFactors;
304  CanonicalForm buf;
305  bool shiftBuf= false;
306  if (!(normFactors.length() == 2 && degree (i.getItem().factor()) <= degree (f)))
307  {
308    TIMING_START (fac_alg_time_shift);
309    if (shift != 0)
310      buf= f (f.mvar() - shift*alpha, f.mvar());
311    else
312      buf= f;
313    shiftBuf= true;
314    TIMING_END_AND_PRINT (fac_alg_time_shift, "time to shift: ");
315  }
316  else
317    buf= f;
318  CanonicalForm factor;
319  int count= 0;
320  for (; i.hasItem(); i++)
321  {
322    ASSERT (i.getItem().exp() == 1, "norm not squarefree");
323    TIMING_START (fac_alg_gcd);
324    if (shiftBuf)
325      factor= gcd (buf, i.getItem().factor());
326    else
327    {
328      if (shift == 0)
329        factor= gcd (buf, i.getItem().factor());
330      else
331        factor= gcd (buf, i.getItem().factor() (f.mvar() + shift*alpha, f.mvar()));
332    }
333    buf /= factor;
334    if (shiftBuf)
335    {
336      if (shift != 0)
337        factor= factor (f.mvar() + shift*alpha, f.mvar());
338    }
339    TIMING_END_AND_PRINT (fac_alg_gcd, "time to recover factors: ");
340    factors.append (factor);
341    count++;
342    if (normFactors.length() - 1 == count)
343    {
344      if (shiftBuf)
345        factors.append (buf (f.mvar() + shift*alpha, f.mvar()));
346      else
347        factors.append (buf);
348      buf= 1;
349      break;
350    }
351  }
352  ASSERT (degree (buf) <= 0, "incomplete factorization");
353  if (save_rat) Off(SW_RATIONAL);
354  return factors;
355}*/
356
357CFFList
358AlgExtFactorize (const CanonicalForm& F, const Variable& alpha)
359{
360  ASSERT (F.isUnivariate(), "univariate input expected");
361  ASSERT (getCharacteristic() == 0, "characteristic 0 expected");
362
363
364  if (F.inCoeffDomain())
365    return CFFList (CFFactor (F, 1));
366
367  bool save_rat=!isOn (SW_RATIONAL);
368  On (SW_RATIONAL);
369  TIMING_START (fac_alg_sqrf);
370  CFFList sqrf= sqrFreeZ (F);
371  TIMING_END_AND_PRINT (fac_alg_sqrf, "time for sqrf factors in Q(a)[x]: ");
372  CFList factorsSqrf;
373  CFFList factors;
374  CFListIterator j;
375
376  CanonicalForm lcinv;
377  for (CFFListIterator i= sqrf; i.hasItem(); i++)
378  {
379    if (i.getItem().factor().inCoeffDomain()) continue;
380    TIMING_START (fac_alg_factor_sqrf);
381    factorsSqrf= AlgExtSqrfFactorize (i.getItem().factor(), alpha);
382    TIMING_END_AND_PRINT (fac_alg_factor_sqrf,
383                          "time to factor sqrf factors in Q(a)[x]: ");
384    for (j= factorsSqrf; j.hasItem(); j++)
385    {
386      lcinv= 1/Lc (j.getItem());
387      factors.append (CFFactor (j.getItem()*lcinv, i.getItem().exp()));
388    }
389  }
390
391  factors.insert (CFFactor (Lc(F), 1));
392  if (save_rat) Off(SW_RATIONAL);
393  return factors;
394}
395
Note: See TracBrowser for help on using the repository browser.