My Project
fac_multihensel.cc
Go to the documentation of this file.
1 /* emacs edit mode for this file is -*- C++ -*- */
2 /* $Id: fac_multihensel.cc 12231 2009-11-02 10:12:22Z hannes $ */
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_eval.h"
12 #include "cf_binom.h"
13 #include "fac_util.h"
14 #include "fac_iterfor.h"
15 #include "cf_iter.h"
16 
17 #ifndef HAVE_NTL
18 
19 TIMING_DEFINE_PRINT(fac_solve);
20 TIMING_DEFINE_PRINT(fac_modpk);
21 TIMING_DEFINE_PRINT(fac_corrcoeff);
22 TIMING_DEFINE_PRINT(fac_extgcd);
23 
24 static void
25 extgcdrest ( const CanonicalForm & a, const CanonicalForm & b, const CanonicalForm & s, const CanonicalForm & t, const CanonicalForm & c, CanonicalForm & S, CanonicalForm & T, const modpk & /*pk*/ )
26 {
27  CanonicalForm sigma = s * c, tau = t * c;
28 // divremainder( sigma, b, T, S, pk );
29 // T = pk( tau + T * a );
30  divrem( sigma, b, T, S );
31  T = tau + T * a;
32 }
33 
34 static void
35 solveF ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, const modpk & pk, int r, CFArray & a )
36 {
37  setCharacteristic( pk.getp(), pk.getk() );
38  CanonicalForm g, bb, b = mapinto( C );
39  int j;
40  for ( j = 1; j < r; j++ )
41  {
42  extgcdrest( mapinto( P[j] ), mapinto( Q[j] ), mapinto( S[j] ), mapinto( T[j] ), b, bb, a[j], pk );
43  b = bb;
44  }
45  a[r] = b;
46  setCharacteristic( 0 );
47  for ( j = 1; j <= r; j++ )
48  a[j] = mapinto( a[j] );
49 }
50 
51 static CanonicalForm
52 evalF ( const CFArray & P, const CFArray & Q, const CFArray & A, int r )
53 {
54  CanonicalForm pprod = 1, sum = 0;
55  for ( int i = 1; i <= r; i++ )
56  {
57  sum += pprod * A[i] * Q[i];
58  pprod *= P[i];
59  }
60  return sum;
61 }
62 
63 static CanonicalForm
64 derivAndEval ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a )
65 {
66  if ( n == 0 )
67  return f( a, x );
68  else if ( f.degree( x ) < n )
69  return 0;
70  else
71  {
72  CFIterator i;
73  CanonicalForm sum = 0, fact;
74  int min, j;
75  Variable v = Variable( f.level() + 1 );
76  for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ )
77  {
78  fact = 1;
79  min = i.exp() - n;
80  for ( j = i.exp(); j > min; j-- )
81  fact *= j;
82  sum += fact * i.coeff() * power( v, min );
83  }
84  return sum( a, v );
85  }
86 }
87 
88 static CanonicalForm
89 checkCombination ( const CanonicalForm & W, const Evaluation & a, const IteratedFor & e, int k )
90 {
91  CanonicalForm dW = W;
92  int i, j;
93  for ( i = k-1; i >= 2 && ! dW.isZero(); i-- )
94  dW = derivAndEval( dW, e[i], Variable( i ), a[i] );
95  if ( ! dW.isZero() ) {
96  CanonicalForm fact = 1;
97  for ( i = 2; i < k; i++ )
98  for ( j = 2; j <= e[i]; j++ )
99  fact *= j;
100  dW /= fact;
101  }
102  return dW;
103 }
104 
105 static CanonicalForm
106 prodCombination ( const Evaluation & a, const IteratedFor & e, int k )
107 {
108  CanonicalForm p = 1;
109  for ( int i = k-1; i > 1; i-- )
110  p *= binomialpower( Variable(i), -a[i], e[i] );
111  return p;
112 }
113 
114 //static CanonicalForm check_dummy( const CFArray &a, const CFArray & P, const CFArray & Q )
115 //{
116 // int i, r = a.size();
117 // CanonicalForm res, prod;
118 // res = 0;
119 // prod = 1;
120 // for ( i = 1; i <= r; i++ )
121 // {
122 // res += prod * a[i] * Q[i];
123 // prod *= P[i];
124 // }
125 // return res;
126 //}
127 
128 static bool check_e( const IteratedFor & e, int k, int m, int * n )
129 {
130  int sum = 0;
131  for ( int i = 2; i < k; i++ )
132  {
133  sum += e[i];
134  if ( e[i] > n[i] )
135  return false;
136  }
137  return sum == m+1;
138 }
139 
140 static CanonicalForm modDeltak ( const CanonicalForm & f, const Evaluation & A, int k, int * n )
141 {
143  for ( int i = 2; i < k; i++ )
144  {
145  result.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) );
146  }
147  return result;
148 }
149 
150 static CFArray
151 findCorrCoeffs ( const CFArray & P, const CFArray & Q, const CFArray & P0, const CFArray & Q0, const CFArray & S, const CFArray & T, const CanonicalForm & C, const Evaluation & I, const modpk & pk, int r, int k, int h, int * n )
152 {
153  DEBINCLEVEL( cerr, "findCorrCoeffs" );
154  int i, m;
155  CFArray A(1,r), a(1,r);
156  CanonicalForm C0, Dm, g, prodcomb;
157 
158  TIMING_START(fac_solve);
159  C0 = pk( I( C, 2, k-1 ), true );
160  solveF( P0, Q0, S, T, 1, pk, r, a );
161  TIMING_END(fac_solve);
162 
163  DEBOUTLN( cerr, "trying to find correction coefficients for " << C );
164  DEBOUTLN( cerr, "which evaluates to " << C0 );
165 
166  for ( i = 1; i <= r; i++ )
167  A[i] = remainder( pk( a[i] * C0 ), P0[i], pk );
168  DEBOUTLN( cerr, "the first approximation of the correction coefficients is " << A );
169 #ifdef DEBUGOUTPUT
170  if ( check_dummy( A, P, Q ) - C != 0 )
171  {
172  DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not" );
173  DEBOUTLN( cerr, "fulfill equation F(A)" );
174  DEBOUTLN( cerr, "corresponding P " << P );
175  DEBOUTLN( cerr, " Q " << Q );
176  }
177 #endif
178  for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ )
179  {
180  Dm = pk( evalF( P, Q, A, r ) - C );
181  if ( Dm != 0 )
182  {
183  if ( k == 2 )
184  {
185  TIMING_START(fac_solve);
186  solveF( P0, Q0, S, T, Dm, pk, r, a );
187  TIMING_END(fac_solve);
188  for ( i = 1; i <= r; i++ )
189  A[i] -= a[i];
190  }
191  else
192  {
193  IteratedFor e( 2, k-1, m+1 );
194  while ( e.iterations_left() )
195  {
196  if ( check_e( e, k, m, n ) )
197  {
198  g = pk( checkCombination( Dm, I, e, k ) );
199  if ( ! g.isZero() && ! (g.mvar() > Variable(1)) )
200  {
201  prodcomb = prodCombination( I, e, k );
202 // Dm = Dm - g * prodcomb;
203  TIMING_START(fac_solve);
204  solveF( P0, Q0, S, T, g, pk, r, a );
205  TIMING_END(fac_solve);
206  for ( i = 1; i <= r; i++ )
207  {
208 // A[i] -= a[i] * prodcomb;
209  A[i] = pk( A[i] - a[i] * prodcomb );
210  }
211  }
212  }
213  e++;
214  }
215  }
216  }
217  DEBOUTLN( cerr, "the correction coefficients at step " << m );
218  DEBOUTLN( cerr, "are now " << A );
219 #ifdef DEBUGOUTPUT
220  if ( check_dummy( A, P, Q ) - C != 0 ) {
221  DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not" );
222  DEBOUTLN( cerr, "fulfill equation F(A)" );
223  DEBOUTLN( cerr, "corresponding P " << P );
224  DEBOUTLN( cerr, " Q " << Q );
225  }
226 #endif
227  }
228  DEBDECLEVEL( cerr, "findCorrCoeffs" );
229  return A;
230 }
231 
232 
233 static bool
234 liftStep ( CFArray & P, int k, int r, int t, const modpk & b, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h )
235 {
236  DEBINCLEVEL( cerr, "liftStep" );
237  CFArray K( 1, r ), Q( 1, r ), Q0( 1, r ), P0( 1, r ), S( 1, r ), T( 1, r ), alpha( 1, r );
238  CanonicalForm Rm, C, D, xa = Variable(k) - A[k];
239  CanonicalForm xa1 = xa, xa2 = xa*xa;
240  CanonicalForm dummy;
241  int i, m;
242 
243  DEBOUTLN( cerr, "we are now performing the liftstep to reach " << Variable(k) );
244  DEBOUTLN( cerr, "the factors so far are " << P );
245  DEBOUTLN( cerr, "modulus p^k= " << b.getpk() << "=" << b.getp() <<"^"<< b.getk() );
246 
247  for ( i = 1; i <= r; i++ )
248  {
249  Variable vm = Variable( t + 1 );
250  Variable v1 = Variable(1);
251  K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), swapvar( A( lcG[i], k+1, t ), v1, vm ) ), v1, vm );
252  P[i] = A( K[i], k, t );
253  }
254  DEBOUTLN( cerr, "lift K = " << K );
255 
256 // d = degree( Uk, Variable( k ) );
257 
258  TIMING_START(fac_extgcd);
259  Q[r] = 1;
260  for ( i = r; i > 1; i-- )
261  {
262  Q[i-1] = Q[i] * P[i];
263  P0[i] = A( P[i], 2, k-1 );
264  Q0[i] = A( Q[i], 2, k-1 );
265  extgcd( P0[i], Q0[i], S[i], T[i], b );
266  }
267  P0[1] = A( P[1], 2, k-1 );
268  Q0[1] = A( Q[1], 2, k-1 );
269  extgcd( P0[1], Q0[1], S[1], T[1], b );
270  TIMING_END(fac_extgcd);
271 
272  for ( m = 1; m <= n[k]+1; m++ )
273  {
274  TIMING_START(fac_modpk);
275  Rm = modDeltak( prod( K ) - Uk, A, k, n );
276  TIMING_END(fac_modpk);
277 #ifdef DEBUGOUTPUT
278  if ( mod( Rm, xa1 ) != 0 )
279  {
280  DEBOUTLN( cerr, "something seems not to be ok with Rm which is " << Rm );
281  DEBOUTLN( cerr, "and should reduce to zero modulo " << xa1 );
282  }
283 #endif
284  if ( mod( Rm, xa2 ) != 0 )
285  {
286  C = derivAndEval( Rm, m, Variable( k ), A[k] );
287  D = 1;
288  for ( i = 2; i <= m; i++ ) D *= i;
289  C /= D;
290 
291  TIMING_START(fac_corrcoeff);
292  alpha = findCorrCoeffs( P, Q, P0, Q0, S, T, C, A, b, r, k, h, n ); // -> h berechnen
293  TIMING_END(fac_corrcoeff);
294 // #ifdef DEBUGOUTPUT
295 // dummy = check_dummy( alpha, P, Q );
296 // if ( b(modDeltak( dummy, A, k, n )) != b(modDeltak( C, A, k, n )) )
297 // {
298 // DEBOUTLN( cerr, "lift fault" );
299 // DEBDECLEVEL( cerr, "liftStep" );
300 // return false;
301 // }
302 // #endif
303  for ( i = 1; i <= r; i++ )
304  K[i] = b(K[i] - alpha[i] * xa1);
305  DEBOUTLN( cerr, "the corrected K's are now " << K );
306  }
307  xa1 = xa2;
308  xa2 *= xa;
309  }
310  for ( i = 1; i <= r; i++ )
311  P[i] = K[i];
312  if ( prod( K ) - Uk != 0 )
313  {
314  DEBOUTLN( cerr, "the liftstep produced the wrong result" );
315  DEBOUTLN( cerr, "the product of the factors calculated so far is " << prod(K) );
316  DEBOUTLN( cerr, "and the Uk that should have been reached is " << Uk );
317  DEBDECLEVEL( cerr, "liftStep" );
318  return false;
319  }
320  DEBOUTLN( cerr, "the lift seems ok so far" );
321  DEBDECLEVEL( cerr, "liftStep" );
322  return true; // check for divisibility
323 }
324 
325 bool
326 Hensel ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const modpk & bound, const Variable & /*x*/ )
327 {
328  DEBINCLEVEL( cerr, "Hensel" );
329  int k, i, h, t = A.max();
330  bool goodeval = true;
331  CFArray Uk( A.min(), A.max() );
332  int * n = new int[t+1];
333 
334  Uk[t] = U;
335  for ( k = t-1; k > 1; k-- )
336  {
337  Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
338  n[k] = degree( Uk[k], Variable( k ) );
339  }
340  for ( k = A.min(); goodeval && (k <= t); k++ )
341  {
342  h = totaldegree( Uk[k], Variable(A.min()), Variable(k-1) );
343  for ( i = A.min(); i <= k; i++ )
344  n[i] = degree( Uk[k], Variable(i) );
345  goodeval = liftStep( G, k, G.max(), t, bound, A, lcG, Uk[k], n, h );
346  DEBOUTLN( cerr, "Factorization so far: " << G );
347  }
348  DEBDECLEVEL( cerr, "Hensel" );
349  delete[] n;
350  return goodeval;
351 }
352 #endif
void divrem(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
CanonicalForm mapinto(const CanonicalForm &f)
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
int degree(const CanonicalForm &f)
void FACTORY_PUBLIC setCharacteristic(int c)
Definition: cf_char.cc:28
CanonicalForm FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
int totaldegree(const CanonicalForm &f)
int totaldegree ( const CanonicalForm & f )
Definition: cf_ops.cc:523
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
g
Definition: cfModGcd.cc:4090
CanonicalForm b
Definition: cfModGcd.cc:4103
void tau(int **points, int sizePoints, int k)
CanonicalForm extgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
Definition: cfUnivarGcd.cc:174
CanonicalForm binomialpower(const Variable &, const CanonicalForm &, int)
factory switches.
evaluate polynomials at points
Iterators for CanonicalForm's.
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
FILE * f
Definition: checklibs.c:9
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
factory's main class
Definition: canonicalform.h:86
CF_NO_INLINE bool isZero() const
class to evaluate a polynomial at points
Definition: cf_eval.h:32
bool iterations_left() const
Definition: fac_iterfor.h:41
factory's class for variables
Definition: factory.h:127
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
int getp() const
Definition: fac_util.h:35
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
Variable alpha
Definition: facAbsBiFact.cc:51
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
TIMING_START(fac_alg_resultant)
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
fq_nmod_poly_t prod
Definition: facHensel.cc:100
bool Hensel(const CanonicalForm &U, CFArray &G, const CFArray &lcG, const Evaluation &A, const modpk &bound, const Variable &)
CanonicalForm remainder(const CanonicalForm &f, const CanonicalForm &g, const modpk &pk)
Definition: fac_util.cc:115
CanonicalForm replaceLc(const CanonicalForm &f, const CanonicalForm &c)
Definition: fac_util.cc:90
operations mod p^k and some other useful functions for factorization
static int min(int a, int b)
Definition: fast_mult.cc:268
#define D(A)
Definition: gentable.cc:131
STATIC_VAR int64 * Q0
Definition: hilb.cc:58
STATIC_VAR jList * T
Definition: janet.cc:30
STATIC_VAR TreeM * G
Definition: janet.cc:31
STATIC_VAR Poly * h
Definition: janet.cc:971
STATIC_VAR jList * Q
Definition: janet.cc:30
#define A
Definition: sirandom.c:24
#define TIMING_DEFINE_PRINT(t)
Definition: timing.h:95
#define TIMING_END(t)
Definition: timing.h:93