My Project
Loading...
Searching...
No Matches
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
19TIMING_DEFINE_PRINT(fac_solve);
20TIMING_DEFINE_PRINT(fac_modpk);
21TIMING_DEFINE_PRINT(fac_corrcoeff);
22TIMING_DEFINE_PRINT(fac_extgcd);
23
24static void
25extgcdrest ( 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
34static void
35solveF ( 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;
47 for ( j = 1; j <= r; j++ )
48 a[j] = mapinto( a[j] );
49}
50
51static CanonicalForm
52evalF ( 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
63static CanonicalForm
64derivAndEval ( 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 {
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
88static CanonicalForm
89checkCombination ( 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
105static CanonicalForm
106prodCombination ( 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
128static 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
140static 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
150static CFArray
151findCorrCoeffs ( 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
233static bool
234liftStep ( 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
325bool
326Hensel ( 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
static int Hensel(const CanonicalForm &UU, CFArray &G, const Evaluation &AA, const CFArray &LeadCoeffs)
Definition: cfEzgcd.cc:302
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: variable.h:33
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
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
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
#define A
Definition: sirandom.c:24
#define Q
Definition: sirandom.c:26
#define TIMING_DEFINE_PRINT(t)
Definition: timing.h:95
#define TIMING_END(t)
Definition: timing.h:93
#define TIMING_START(t)
Definition: timing.h:92