source: git/factory/facAlgExt.cc @ 386b3d

spielwiese
Last change on this file since 386b3d was 4f6d99, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: better recovery of factors in univariate factorization over Q(a)
  • Property mode set to 100644
File size: 6.2 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#include "config.h"
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
50// i is an integer such that Norm (F (x-i*alpha)) is squarefree
51CanonicalForm sqrfNorm (const CanonicalForm& F, const Variable& alpha, int& i)
52{
53  Variable x= Variable (F.level() + 1);
54  Variable y= F.mvar();
55  CanonicalForm g= F (x, alpha);
56  CanonicalForm mipo= getMipo (alpha);
57  mipo= mipo (x, alpha);
58  mipo *= bCommonDen (mipo);
59
60  int degg= degree (g);
61  int degmipo= degree (mipo);
62  CanonicalForm norm;
63  TIMING_START (fac_alg_resultant);
64  if (degg >= 8 || degmipo >= 8)
65    norm= resultantZ (g, mipo, x);
66  else
67    norm= resultant (g, mipo, x);
68  TIMING_END_AND_PRINT (fac_alg_resultant, "time to compute resultant0: ");
69
70  i= 0;
71  int k;
72  if (degree (gcd (deriv (norm, y), norm)) <= 0)
73    return norm;
74  i= 1;
75  do
76  {
77    k= 1;
78    while (k < 3)
79    {
80      if (k == 1)
81      {
82        g= F (y - i*alpha, y);
83        g *= bCommonDen (g);
84        TIMING_START (fac_alg_resultant);
85        if (degg >= 8 || degmipo >= 8)
86          norm= resultantZ (g (x, alpha), mipo, x);
87        else
88          norm= resultant (g (x, alpha), mipo, x);
89        TIMING_END_AND_PRINT (fac_alg_resultant,"time to compute resultant1: ");
90      }
91      else
92      {
93        g= F (y + i*alpha, y);
94        g *= bCommonDen (g);
95        TIMING_START (fac_alg_resultant);
96        if (degg >= 8 || degmipo >= 8)
97          norm= resultantZ (g (x, alpha), mipo, x);
98        else
99          norm= resultant (g (x, alpha), mipo, x);
100        TIMING_END_AND_PRINT (fac_alg_resultant,"time to compute resultant2: ");
101      }
102      if (degree (gcd (deriv (norm, y), norm)) <= 0)
103      {
104        if (k == 2)
105          i= -i;
106        return norm;
107      }
108      k++;
109    }
110    i++;
111  } while (1);
112}
113
114CFList
115AlgExtSqrfFactorize (const CanonicalForm& F, const Variable& alpha)
116{
117  ASSERT (F.isUnivariate(), "univariate input expected");
118  ASSERT (getCharacteristic() == 0, "characteristic 0 expected");
119
120  bool save_rat=!isOn (SW_RATIONAL);
121  On (SW_RATIONAL);
122  CanonicalForm f= F*bCommonDen (F);
123  int shift;
124  TIMING_START (fac_alg_norm);
125  CanonicalForm norm= sqrfNorm (f, alpha, shift);
126  TIMING_END_AND_PRINT (fac_alg_norm, "time to compute sqrf norm: ");
127  ASSERT (degree (norm, alpha) <= 0, "wrong norm computed");
128  TIMING_START (fac_alg_factor_norm);
129  bool save_sort= !isOn (SW_USE_NTL_SORT);
130  On (SW_USE_NTL_SORT);
131  CFFList normFactors= factorize (norm);
132  if (save_sort)
133    Off (SW_USE_NTL_SORT);
134  TIMING_END_AND_PRINT (fac_alg_factor_norm, "time to factor norm: ");
135  CFList factors;
136  if (normFactors.length() <= 2)
137  {
138    if (save_rat) Off(SW_RATIONAL);
139    return CFList (F);
140  }
141
142  normFactors.removeFirst();
143  CFFListIterator i= normFactors;
144  CanonicalForm buf;
145  bool shiftBuf= false;
146  if (!(normFactors.length() == 2 && degree (i.getItem().factor()) <= degree (f)))
147  {
148    TIMING_START (fac_alg_time_shift);
149    if (shift != 0)
150      buf= f (f.mvar() - shift*alpha, f.mvar());
151    else
152      buf= f;
153    shiftBuf= true;
154    TIMING_END_AND_PRINT (fac_alg_time_shift, "time to shift: ");
155  }
156  else
157    buf= f;
158  CanonicalForm factor;
159  int count= 0;
160  for (; i.hasItem(); i++)
161  {
162    ASSERT (i.getItem().exp() == 1, "norm not squarefree");
163    TIMING_START (fac_alg_gcd);
164    if (shiftBuf)
165      factor= gcd (buf, i.getItem().factor());
166    else
167    {
168      if (shift == 0)
169        factor= gcd (buf, i.getItem().factor());
170      else
171        factor= gcd (buf, i.getItem().factor() (f.mvar() + shift*alpha, f.mvar()));
172    }
173    buf /= factor;
174    if (shiftBuf)
175    {
176      if (shift != 0)
177        factor= factor (f.mvar() + shift*alpha, f.mvar());
178    }
179    TIMING_END_AND_PRINT (fac_alg_gcd, "time to recover factors: ");
180    factors.append (factor);
181    count++;
182    if (normFactors.length() - 1 == count)
183    {
184      if (shiftBuf)
185        factors.append (buf (f.mvar() + shift*alpha, f.mvar()));
186      else
187        factors.append (buf);
188      buf= 1;
189      break;
190    }
191  }
192  ASSERT (degree (buf) <= 0, "incomplete factorization");
193  if (save_rat) Off(SW_RATIONAL);
194  return factors;
195}
196
197CFFList
198AlgExtFactorize (const CanonicalForm& F, const Variable& alpha)
199{
200  ASSERT (F.isUnivariate(), "univariate input expected");
201  ASSERT (getCharacteristic() == 0, "characteristic 0 expected");
202
203
204  if (F.inCoeffDomain())
205    return CFFList (CFFactor (F, 1));
206
207  bool save_rat=!isOn (SW_RATIONAL);
208  On (SW_RATIONAL);
209  TIMING_START (fac_alg_sqrf);
210  CFFList sqrf= sqrFreeZ (F);
211  TIMING_END_AND_PRINT (fac_alg_sqrf, "time for sqrf factors in Q(a)[x]: ");
212  CFList factorsSqrf;
213  CFFList factors;
214  CFListIterator j;
215
216  CanonicalForm lcinv;
217  for (CFFListIterator i= sqrf; i.hasItem(); i++)
218  {
219    if (i.getItem().factor().inCoeffDomain()) continue;
220    TIMING_START (fac_alg_factor_sqrf);
221    factorsSqrf= AlgExtSqrfFactorize (i.getItem().factor(), alpha);
222    TIMING_END_AND_PRINT (fac_alg_factor_sqrf,
223                          "time to factor sqrf factors in Q(a)[x]: ");
224    for (j= factorsSqrf; j.hasItem(); j++)
225    {
226      lcinv= 1/Lc (j.getItem());
227      factors.append (CFFactor (j.getItem()*lcinv, i.getItem().exp()));
228    }
229  }
230
231  factors.insert (CFFactor (Lc(F), 1));
232  if (save_rat) Off(SW_RATIONAL);
233  return factors;
234}
235
Note: See TracBrowser for help on using the repository browser.