My Project
fac_multivar.cc
Go to the documentation of this file.
1 /* emacs edit mode for this file is -*- C++ -*- */
2 /* $Id: fac_multivar.cc 14377 2011-09-01 13:40:30Z mlee $ */
3 
4 #include <config.h>
5 
6 #include "assert.h"
7 #include "debug.h"
8 #include "timing.h"
9 
10 #include "cf_defs.h"
11 #include "cf_algorithm.h"
12 #include "fac_multivar.h"
13 #include "fac_univar.h"
14 #include "cf_reval.h"
15 #include "cf_map.h"
16 #include "fac_util.h"
17 #include "cf_binom.h"
18 #include "cf_iter.h"
19 #include "cf_primes.h"
20 #include "fac_distrib.h"
21 #include "fac_multihensel.h"
22 #include "facBivar.h"
23 
24 #ifndef HAVE_NTL
25 
26 void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
27 void out_cff(CFFList &L);
28 
29 TIMING_DEFINE_PRINT(fac_content);
30 TIMING_DEFINE_PRINT(fac_findeval);
31 TIMING_DEFINE_PRINT(fac_distrib);
32 TIMING_DEFINE_PRINT(fac_hensel);
33 
34 static CFArray
35 conv_to_factor_array( const CFFList & L )
36 {
37  int n;
38  CFFListIterator I = L;
39  bool negate = false;
40 
41  if ( ! I.hasItem() )
42  n = 0;
43  else if ( I.getItem().factor().inBaseDomain() ) {
44  negate = I.getItem().factor().sign() < 0;
45  I++;
46  n = L.length();
47  }
48  else
49  n = L.length() + 1;
50  CFFListIterator J = I;
51  while ( J.hasItem() ) {
52  n += J.getItem().exp() - 1;
53  J++;
54  }
55  CFArray result( 1, n-1 );
56  int i, j, k;
57  i = 1;
58  while ( I.hasItem() ) {
59  k = I.getItem().exp();
60  for ( j = 1; j <= k; j++ ) {
61  result[i] = I.getItem().factor();
62  i++;
63  }
64  I++;
65  }
66  if ( negate )
67  result[1] = -result[1];
68  return result;
69 }
70 
71 static modpk
72 coeffBound_old ( const CanonicalForm & f, int p )
73 {
74  int * degs = degrees( f );
75  int M = 0, i, k = f.level();
76  for ( i = 1; i <= k; i++ )
77  M += degs[i];
78  CanonicalForm b = 2 * maxNorm( f ) * power( CanonicalForm( 3 ), M );
79  CanonicalForm B = p;
80  k = 1;
81  while ( B < b ) {
82  B *= p;
83  k++;
84  }
85  return modpk( p, k );
86 }
87 
88 // static bool
89 // nonDivisors ( CanonicalForm omega, CanonicalForm delta, const CFArray & F, CFArray & d )
90 // {
91 // DEBOUTLN( cerr, "nondivisors omega = " << omega );
92 // DEBOUTLN( cerr, "nondivisors delta = " << delta );
93 // DEBOUTLN( cerr, "nondivisors F = " << F );
94 // CanonicalForm q, r;
95 // int k = F.size();
96 // d = CFArray( 0, k );
97 // d[0] = delta * omega;
98 // for ( int i = 1; i <= k; i++ ) {
99 // q = abs(F[i]);
100 // for ( int j = i-1; j >= 0; j-- ) {
101 // r = d[j];
102 // do {
103 // r = gcd( r, q );
104 // q = q / r;
105 // } while ( r != 1 );
106 // if ( q == 1 )
107 // return false;
108 // }
109 // d[i] = q;
110 // }
111 // return true;
112 // }
113 
114 static void
115 findEvaluation ( const CanonicalForm & U, const CanonicalForm & V, const CanonicalForm & omega, const CFFList & F, Evaluation & A, CanonicalForm & U0, CanonicalForm & delta, CFArray & D, int r )
116 {
117  DEBINCLEVEL( cerr, "findEvaluation" );
118  CanonicalForm Vn;
119  CFFListIterator I;
120  int j;
121  bool found = false;
122  CFArray FF = CFArray( 1, F.length() );
123  if ( r > 0 )
124  A.nextpoint();
125  while ( ! found )
126  {
127  Vn = A( V );
128  if ( Vn != 0 )
129  {
130  U0 = A( U );
131  DEBOUTLN( cerr, "U0 = " << U0 );
132  if ( isSqrFree( U0 ) )
133  {
134  delta = content( U0 );
135  DEBOUTLN( cerr, "content( U0 ) = " << delta );
136  for ( I = F, j = 1; I.hasItem(); I++, j++ )
137  FF[j] = A( I.getItem().factor() );
138  found = nonDivisors( omega, delta, FF, D );
139  }
140  else
141  {
142  DEBOUTLN( cerr, "not sqrfree : " << sqrFree( U0 ) );
143  }
144  }
145  if ( ! found )
146  A.nextpoint();
147  }
148  DEBDECLEVEL( cerr, "findEvaluation" );
149 }
150 
151 
152 static int prime_number=0;
153 void find_good_prime(const CanonicalForm &f, int &start)
154 {
155  if (! f.inBaseDomain() )
156  {
157  CFIterator i = f;
158  for(;;)
159  {
160  if ( i.hasTerms() )
161  {
162  find_good_prime(i.coeff(),start);
163  if (start==cf_getNumSmallPrimes()) return;
164  if((i.exp()!=0) && ((i.exp() % cf_getSmallPrime(start))==0))
165  {
166  start++;
167  if (start==cf_getNumSmallPrimes()) return;
168  i=f;
169  }
170  else i++;
171  }
172  else break;
173  }
174  }
175  else
176  {
177  if (f.inZ())
178  {
179  if (start==cf_getNumSmallPrimes()) return;
180  while((!f.isZero()) && (mod(f,cf_getSmallPrime(start))==0))
181  {
182  start++;
183  if (start==cf_getNumSmallPrimes()) return;
184  }
185  }
186 /* should not happen!
187  else if (f.inQ())
188  {
189  while((f.den()!=0) && (mod(f.den(),cf_getSmallPrime(start))==0))
190  {
191  start++;
192  }
193  while((f.num()!=0) && (mod(f.num(),cf_getSmallPrime(start))==0))
194  {
195  start++;
196  }
197  }
198  else
199  cout <<"??"<< f <<"\n";
200 */
201  }
202 }
203 static CFArray ZFactorizeMulti ( const CanonicalForm & arg )
204 {
205  prime_number=0;
206  bool is_rat=isOn(SW_RATIONAL);
207  Off(SW_RATIONAL);
208  DEBINCLEVEL( cerr, "ZFactorizeMulti" );
209  CFMap M;
210  CanonicalForm UU, U = compress( arg, M );
211  CanonicalForm delta, omega, V = LC( U, 1 );
212  int t = U.level();
213  CFFList F = factorize( V );
214  CFFListIterator I, J;
215  CFArray G, lcG, D;
216  int i, j, r, maxdeg;
217  REvaluation A( 2, t, IntRandom( 50 ) );
218  CanonicalForm U0;
219  CanonicalForm ft, ut, gt, d;
220  modpk b;
221  bool negate = false;
222 
223  DEBOUTLN( cerr, "-----------------------------------------------------" );
224  DEBOUTLN( cerr, "trying to factorize U = " << U );
225  DEBOUTLN( cerr, "U is a polynomial of level = " << arg.level() );
226  DEBOUTLN( cerr, "U will be factorized with respect to variable " << Variable(1) );
227  DEBOUTLN( cerr, "the leading coefficient of U with respect to that variable is " << V );
228  DEBOUTLN( cerr, "which is factorized as " << F );
229 
230  maxdeg = 0;
231  for ( i = 2; i <= t; i++ )
232  {
233  j = U.degree( Variable( i ) );
234  if ( j > maxdeg ) maxdeg = j;
235  }
236 
237  if ( F.getFirst().factor().inCoeffDomain() )
238  {
239  omega = F.getFirst().factor();
240  F.removeFirst();
241  if ( omega < 0 )
242  {
243  negate = true;
244  omega = -omega;
245  U = -U;
246  }
247  }
248  else
249  omega = 1;
250 
251  bool goodeval = false;
252  r = 0;
253 
254 // for ( i = 0; i < 10*t; i++ )
255 // A.nextpoint();
256 
257  while ( ! goodeval )
258  {
259  TIMING_START(fac_findeval);
260  findEvaluation( U, V, omega, F, A, U0, delta, D, r );
261  TIMING_END(fac_findeval);
262  DEBOUTLN( cerr, "the evaluation point to reduce to an univariate problem is " << A );
263  DEBOUTLN( cerr, "corresponding delta = " << delta );
264  DEBOUTLN( cerr, " omega = " << omega );
265  DEBOUTLN( cerr, " D = " << D );
266  DEBOUTLN( cerr, "now factorize the univariate polynomial " << U0 );
267  G = conv_to_factor_array( factorize( U0, false ) );
268  DEBOUTLN( cerr, "which factorizes into " << G );
269  {
270  int i=prime_number;
271  find_good_prime(arg,i);
272  find_good_prime(U0,i);
273  find_good_prime(U,i);
274  int p;
275  if (i==cf_getNumSmallPrimes()) p=0;
276  else p=cf_getSmallPrime(i);
277  //printf("found:p=%d (%d)\n",p,i);
278  if (p==0)
279  {
280  return conv_to_factor_array(CFFactor(arg,1));
281  //printf("out of primes - switch to non-NTL\n");
282  }
283  else if (((i==0)||(i!=prime_number)))
284  {
285  b = coeffBound_old( U, p );
286  prime_number=i;
287  }
288  else prime_number++;
289  // p!=0:
290  modpk bb=coeffBound_old(U0,p);
291  if (bb.getk() > b.getk() ) b=bb;
292  bb=coeffBound_old(arg,p);
293  if (bb.getk() > b.getk() ) b=bb;
294  }
295  if ( getZFacModulus().getpk() > b.getpk() )
296  b = getZFacModulus();
297  //printf("p=%d, k=%d\n",b.getp(),b.getk());
298  DEBOUTLN( cerr, "the coefficient bound of the factors of U is " << b.getpk() );
299 
300  r = G.size();
301  lcG = CFArray( 1, r );
302  UU = U;
303  DEBOUTLN( cerr, "now trying to distribute the leading coefficients ..." );
304  TIMING_START(fac_distrib);
305  goodeval = distributeLeadingCoeffs( UU, G, lcG, F, D, delta, omega, A, r );
306  TIMING_END(fac_distrib);
307 #ifdef DEBUGOUTPUT
308  if ( goodeval )
309  {
310  DEBOUTLN( cerr, "the univariate factors after distribution are " << G );
311  DEBOUTLN( cerr, "the distributed leading coeffs are " << lcG );
312  DEBOUTLN( cerr, "U may have changed and is now " << UU );
313  DEBOUTLN( cerr, "which has leading coefficient " << LC( UU, Variable(1) ) );
314 
315  if ( LC( UU, Variable(1) ) != prod( lcG ) || A(UU) != prod( G ) )
316  {
317  DEBOUTLN( cerr, "!!! distribution was not correct !!!" );
318  DEBOUTLN( cerr, "product of leading coeffs is " << prod( lcG ) );
319  DEBOUTLN( cerr, "product of univariate factors is " << prod( G ) );
320  DEBOUTLN( cerr, "the new U is evaluated as " << A(UU) );
321  }
322  else
323  DEBOUTLN( cerr, "leading coeffs correct" );
324  }
325  else
326  {
327  DEBOUTLN( cerr, "we have found a bad evaluation point" );
328  }
329 #endif
330  if ( goodeval )
331  {
332  TIMING_START(fac_hensel);
333  goodeval = Hensel( UU, G, lcG, A, b, Variable(1) );
334  TIMING_END(fac_hensel);
335  }
336  }
337  for ( i = 1; i <= r; i++ )
338  {
339  G[i] /= icontent( G[i] );
340  G[i] = M(G[i]);
341  }
342  // negate noch beachten !
343  if ( negate )
344  G[1] = -G[1];
345  DEBDECLEVEL( cerr, "ZFactorMulti" );
346  if(is_rat) On(SW_RATIONAL);
347  return G;
348 }
349 
350 CFFList ZFactorizeMultivariate ( const CanonicalForm & f, bool issqrfree )
351 {
352  CFFList G, F, R;
353  CFArray GG;
355  CFMap M;
356  CanonicalForm g, cont;
357  Variable v1, vm;
358  int k, m, n;
359 
360  v1 = Variable(1);
361  if ( issqrfree )
362  F = CFFactor( f, 1 );
363  else
364  F = sqrFree( f );
365 
366  for ( i = F; i.hasItem(); i++ )
367  {
368  if ( i.getItem().factor().inCoeffDomain() )
369  {
370  if ( ! i.getItem().factor().isOne() )
371  R.append( CFFactor( i.getItem().factor(), i.getItem().exp() ) );
372  }
373  else
374  {
375  TIMING_START(fac_content);
376  g = compress( i.getItem().factor(), M );
377  // now after compress g contains Variable(1)
378  vm = g.mvar();
379  g = swapvar( g, v1, vm );
380  cont = content( g );
381  g = swapvar( g / cont, v1, vm );
382  cont = swapvar( cont, v1, vm );
383  n = i.getItem().exp();
384  TIMING_END(fac_content);
385  DEBOUTLN( cerr, "now after content ..." );
386  if ( g.isUnivariate() )
387  {
388  G = factorize( g, true );
389  for ( j = G; j.hasItem(); j++ )
390  if ( ! j.getItem().factor().isOne() )
391  R.append( CFFactor( M( j.getItem().factor() ), n ) );
392  }
393  else
394  {
395  GG = ZFactorizeMulti( g );
396  m = GG.max();
397  for ( k = GG.min(); k <= m; k++ )
398  if ( ! GG[k].isOne() )
399  R.append( CFFactor( M( GG[k] ), n ) );
400  }
401  G = factorize( cont, true );
402  for ( j = G; j.hasItem(); j++ )
403  if ( ! j.getItem().factor().isOne() )
404  R.append( CFFactor( M( j.getItem().factor() ), n ) );
405  }
406  }
407  return R;
408 }
409 #endif
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
CanonicalForm icontent(const CanonicalForm &f)
CanonicalForm icontent ( const CanonicalForm & f )
Definition: cf_gcd.cc:74
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Array< CanonicalForm > CFArray
CanonicalForm swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
Factor< CanonicalForm > CFFactor
CanonicalForm content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:603
int * degrees(const CanonicalForm &f, int *degs=0)
int * degrees ( const CanonicalForm & f, int * degs )
Definition: cf_ops.cc:493
CanonicalForm LC(const CanonicalForm &f)
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
int p
Definition: cfModGcd.cc:4080
g
Definition: cfModGcd.cc:4092
const CanonicalForm & GG
Definition: cfModGcd.cc:4078
CanonicalForm b
Definition: cfModGcd.cc:4105
CanonicalForm maxNorm(const CanonicalForm &f)
CanonicalForm maxNorm ( const CanonicalForm & f )
declarations of higher level algorithms.
CFFList FACTORY_PUBLIC sqrFree(const CanonicalForm &f, bool sort=false)
squarefree factorization
Definition: cf_factor.cc:946
CFFList FACTORY_PUBLIC factorize(const CanonicalForm &f, bool issqrfree=false)
factorization over or
Definition: cf_factor.cc:405
factory switches.
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:30
void out_cff(CFFList &L)
Definition: cf_factor.cc:202
Iterators for CanonicalForm's.
CanonicalForm compress(const CanonicalForm &f, CFMap &m)
CanonicalForm compress ( const CanonicalForm & f, CFMap & m )
Definition: cf_map.cc:210
map polynomials
int cf_getNumSmallPrimes()
Definition: cf_primes.cc:34
int cf_getSmallPrime(int i)
Definition: cf_primes.cc:28
access to prime tables
generate random evaluation points
FILE * f
Definition: checklibs.c:9
void out_cf(const char *s1, const CanonicalForm &f, const char *s2)
Definition: cf_factor.cc:99
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
class CFMap
Definition: cf_map.h:85
factory's main class
Definition: canonicalform.h:86
int degree() const
Returns -1 for the zero polynomial and 0 if CO is in a base domain.
int level() const
level() returns the level of CO.
class to evaluate a polynomial at points
Definition: cf_eval.h:32
generate random integers
Definition: cf_random.h:56
T & getItem() const
Definition: ftmpl_list.cc:431
T getFirst() const
Definition: ftmpl_list.cc:279
void removeFirst()
Definition: ftmpl_list.cc:287
int length() const
Definition: ftmpl_list.cc:273
class to generate random evaluation points
Definition: cf_reval.h:26
factory's class for variables
Definition: factory.h:134
class to do operations mod p^k for int's p and k
Definition: fac_util.h:23
int getk() const
Definition: fac_util.h:36
functions to print debug output
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45
return result
Definition: facAbsBiFact.cc:75
TIMING_START(fac_alg_resultant)
return modpk(p, k)
b *CanonicalForm B
Definition: facBivar.cc:52
bivariate factorization over Q(a)
bool found
Definition: facFactorize.cc:55
int j
Definition: facHensel.cc:110
fq_nmod_poly_t prod
Definition: facHensel.cc:100
bool nonDivisors(CanonicalForm omega, CanonicalForm delta, const CFArray &F, CFArray &d)
bool distributeLeadingCoeffs(CanonicalForm &U, CFArray &G, CFArray &lcG, const CFFList &F, const CFArray &D, CanonicalForm &delta, CanonicalForm &omega, const Evaluation &A, int r)
bool Hensel(const CanonicalForm &U, CFArray &G, const CFArray &lcG, const Evaluation &A, const modpk &bound, const Variable &)
CFFList ZFactorizeMultivariate(const CanonicalForm &f, bool issqrfree)
bool isSqrFree(const CanonicalForm &f)
modpk getZFacModulus()
operations mod p^k and some other useful functions for factorization
#define D(A)
Definition: gentable.cc:131
STATIC_VAR TreeM * G
Definition: janet.cc:31
bool delta(X x, Y y, D d)
Definition: TestSuite.h:160
#define R
Definition: sirandom.c:27
#define A
Definition: sirandom.c:24
#define M
Definition: sirandom.c:25
#define TIMING_DEFINE_PRINT(t)
Definition: timing.h:95
#define TIMING_END(t)
Definition: timing.h:93