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