source: git/factory/facAlgExt.cc @ 80b8fe

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