source: git/factory/facAlgExt.cc @ ae8a3a

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