source: git/factory/facAlgExt.cc @ c1b52b

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