[493c477] | 1 | /* emacs edit mode for this file is -*- C++ -*- */ |
---|
[341696] | 2 | /* $Id$ */ |
---|
[ba1fde] | 3 | |
---|
[e4fe2b] | 4 | #include "config.h" |
---|
[173e86] | 5 | |
---|
[650f2d8] | 6 | #include "cf_assert.h" |
---|
[043814] | 7 | #include "debug.h" |
---|
[ba1fde] | 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 | |
---|
[e76d7a6] | 18 | TIMING_DEFINE_PRINT(fac_solve) |
---|
| 19 | TIMING_DEFINE_PRINT(fac_modpk) |
---|
| 20 | TIMING_DEFINE_PRINT(fac_corrcoeff) |
---|
| 21 | TIMING_DEFINE_PRINT(fac_extgcd) |
---|
[ba1fde] | 22 | |
---|
| 23 | static void |
---|
[5b8726d] | 24 | extgcdrest ( const CanonicalForm & a, const CanonicalForm & b, const CanonicalForm & s, const CanonicalForm & t, const CanonicalForm & c, CanonicalForm & S, CanonicalForm & T, const modpk & /*pk*/ ) |
---|
[ba1fde] | 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; |
---|
[01e8874] | 39 | for ( j = 1; j < r; j++ ) |
---|
| 40 | { |
---|
[157108] | 41 | extgcdrest( mapinto( P[j] ), mapinto( Q[j] ), mapinto( S[j] ), mapinto( T[j] ), b, bb, a[j], pk ); |
---|
| 42 | b = bb; |
---|
[ba1fde] | 43 | } |
---|
| 44 | a[r] = b; |
---|
| 45 | setCharacteristic( 0 ); |
---|
| 46 | for ( j = 1; j <= r; j++ ) |
---|
[157108] | 47 | a[j] = mapinto( a[j] ); |
---|
[ba1fde] | 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; |
---|
[01e8874] | 54 | for ( int i = 1; i <= r; i++ ) |
---|
| 55 | { |
---|
[157108] | 56 | sum += pprod * A[i] * Q[i]; |
---|
| 57 | pprod *= P[i]; |
---|
[ba1fde] | 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 ) |
---|
[157108] | 66 | return f( a, x ); |
---|
[ba1fde] | 67 | else if ( f.degree( x ) < n ) |
---|
[157108] | 68 | return 0; |
---|
[01e8874] | 69 | else |
---|
| 70 | { |
---|
[157108] | 71 | CFIterator i; |
---|
| 72 | CanonicalForm sum = 0, fact; |
---|
| 73 | int min, j; |
---|
| 74 | Variable v = Variable( f.level() + 1 ); |
---|
[01e8874] | 75 | for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ ) |
---|
| 76 | { |
---|
[157108] | 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 ); |
---|
[ba1fde] | 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-- ) |
---|
[157108] | 93 | dW = derivAndEval( dW, e[i], Variable( i ), a[i] ); |
---|
[ba1fde] | 94 | if ( ! dW.isZero() ) { |
---|
[157108] | 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; |
---|
[ba1fde] | 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-- ) |
---|
[157108] | 109 | p *= binomialpower( Variable(i), -a[i], e[i] ); |
---|
[ba1fde] | 110 | return p; |
---|
| 111 | } |
---|
| 112 | |
---|
[01e8874] | 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 | //} |
---|
[ba1fde] | 126 | |
---|
[01e8874] | 127 | static bool check_e( const IteratedFor & e, int k, int m, int * n ) |
---|
[ba1fde] | 128 | { |
---|
| 129 | int sum = 0; |
---|
[01e8874] | 130 | for ( int i = 2; i < k; i++ ) |
---|
| 131 | { |
---|
[157108] | 132 | sum += e[i]; |
---|
| 133 | if ( e[i] > n[i] ) |
---|
| 134 | return false; |
---|
[ba1fde] | 135 | } |
---|
| 136 | return sum == m+1; |
---|
| 137 | } |
---|
| 138 | |
---|
[01e8874] | 139 | static CanonicalForm modDeltak ( const CanonicalForm & f, const Evaluation & A, int k, int * n ) |
---|
[ba1fde] | 140 | { |
---|
| 141 | CanonicalForm result = f; |
---|
[01e8874] | 142 | for ( int i = 2; i < k; i++ ) |
---|
| 143 | { |
---|
[157108] | 144 | result.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) ); |
---|
[ba1fde] | 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 | |
---|
[d81ed62] | 162 | DEBOUTLN( cerr, "trying to find correction coefficients for " << C ); |
---|
| 163 | DEBOUTLN( cerr, "which evaluates to " << C0 ); |
---|
[ba1fde] | 164 | |
---|
| 165 | for ( i = 1; i <= r; i++ ) |
---|
[157108] | 166 | A[i] = remainder( pk( a[i] * C0 ), P0[i], pk ); |
---|
[d81ed62] | 167 | DEBOUTLN( cerr, "the first approximation of the correction coefficients is " << A ); |
---|
[f3a82f4] | 168 | /*#ifdef DEBUGOUTPUT |
---|
[01e8874] | 169 | if ( check_dummy( A, P, Q ) - C != 0 ) |
---|
| 170 | { |
---|
[157108] | 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 ); |
---|
[ba1fde] | 175 | } |
---|
[f3a82f4] | 176 | #endif*/ |
---|
[01e8874] | 177 | for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ ) |
---|
| 178 | { |
---|
[157108] | 179 | Dm = pk( evalF( P, Q, A, r ) - C ); |
---|
[01e8874] | 180 | if ( Dm != 0 ) |
---|
| 181 | { |
---|
| 182 | if ( k == 2 ) |
---|
| 183 | { |
---|
[157108] | 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 | } |
---|
[01e8874] | 190 | else |
---|
| 191 | { |
---|
[157108] | 192 | IteratedFor e( 2, k-1, m+1 ); |
---|
[01e8874] | 193 | while ( e.iterations_left() ) |
---|
| 194 | { |
---|
| 195 | if ( check_e( e, k, m, n ) ) |
---|
| 196 | { |
---|
[157108] | 197 | g = pk( checkCombination( Dm, I, e, k ) ); |
---|
[01e8874] | 198 | if ( ! g.isZero() && ! (g.mvar() > Variable(1)) ) |
---|
| 199 | { |
---|
[157108] | 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); |
---|
[01e8874] | 205 | for ( i = 1; i <= r; i++ ) |
---|
| 206 | { |
---|
[157108] | 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 ); |
---|
[f3a82f4] | 218 | /*#ifdef DEBUGOUTPUT |
---|
[ba1fde] | 219 | if ( check_dummy( A, P, Q ) - C != 0 ) { |
---|
[157108] | 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 ); |
---|
[ba1fde] | 224 | } |
---|
[f3a82f4] | 225 | #endif*/ |
---|
[ba1fde] | 226 | } |
---|
| 227 | DEBDECLEVEL( cerr, "findCorrCoeffs" ); |
---|
| 228 | return A; |
---|
| 229 | } |
---|
[043814] | 230 | |
---|
[ba1fde] | 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 | |
---|
[d81ed62] | 242 | DEBOUTLN( cerr, "we are now performing the liftstep to reach " << Variable(k) ); |
---|
| 243 | DEBOUTLN( cerr, "the factors so far are " << P ); |
---|
[cae0b6] | 244 | DEBOUTLN( cerr, "modulus p^k= " << b.getpk() << "=" << b.getp() <<"^"<< b.getk() ); |
---|
[ba1fde] | 245 | |
---|
[01e8874] | 246 | for ( i = 1; i <= r; i++ ) |
---|
| 247 | { |
---|
[157108] | 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 ); |
---|
[ba1fde] | 252 | } |
---|
[d81ed62] | 253 | DEBOUTLN( cerr, "lift K = " << K ); |
---|
[ba1fde] | 254 | |
---|
| 255 | // d = degree( Uk, Variable( k ) ); |
---|
| 256 | |
---|
| 257 | TIMING_START(fac_extgcd); |
---|
| 258 | Q[r] = 1; |
---|
[01e8874] | 259 | for ( i = r; i > 1; i-- ) |
---|
| 260 | { |
---|
[157108] | 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 ); |
---|
[ba1fde] | 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 | |
---|
[01e8874] | 271 | for ( m = 1; m <= n[k]+1; m++ ) |
---|
| 272 | { |
---|
[157108] | 273 | TIMING_START(fac_modpk); |
---|
| 274 | Rm = modDeltak( prod( K ) - Uk, A, k, n ); |
---|
| 275 | TIMING_END(fac_modpk); |
---|
[ba1fde] | 276 | #ifdef DEBUGOUTPUT |
---|
[01e8874] | 277 | if ( mod( Rm, xa1 ) != 0 ) |
---|
| 278 | { |
---|
[157108] | 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 | } |
---|
[ba1fde] | 282 | #endif |
---|
[01e8874] | 283 | if ( mod( Rm, xa2 ) != 0 ) |
---|
| 284 | { |
---|
[157108] | 285 | C = derivAndEval( Rm, m, Variable( k ), A[k] ); |
---|
| 286 | D = 1; |
---|
| 287 | for ( i = 2; i <= m; i++ ) D *= i; |
---|
| 288 | C /= D; |
---|
[ba1fde] | 289 | |
---|
[157108] | 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); |
---|
[ba1fde] | 293 | // #ifdef DEBUGOUTPUT |
---|
[157108] | 294 | // dummy = check_dummy( alpha, P, Q ); |
---|
[01e8874] | 295 | // if ( b(modDeltak( dummy, A, k, n )) != b(modDeltak( C, A, k, n )) ) |
---|
| 296 | // { |
---|
[157108] | 297 | // DEBOUTLN( cerr, "lift fault" ); |
---|
| 298 | // DEBDECLEVEL( cerr, "liftStep" ); |
---|
| 299 | // return false; |
---|
| 300 | // } |
---|
[ba1fde] | 301 | // #endif |
---|
[157108] | 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; |
---|
[ba1fde] | 308 | } |
---|
| 309 | for ( i = 1; i <= r; i++ ) |
---|
[157108] | 310 | P[i] = K[i]; |
---|
[01e8874] | 311 | if ( prod( K ) - Uk != 0 ) |
---|
| 312 | { |
---|
[157108] | 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; |
---|
[ba1fde] | 318 | } |
---|
[d81ed62] | 319 | DEBOUTLN( cerr, "the lift seems ok so far" ); |
---|
[ba1fde] | 320 | DEBDECLEVEL( cerr, "liftStep" ); |
---|
| 321 | return true; // check for divisibility |
---|
| 322 | } |
---|
| 323 | |
---|
| 324 | bool |
---|
[5b8726d] | 325 | Hensel ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const modpk & bound, const Variable & /*x*/ ) |
---|
[ba1fde] | 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; |
---|
[01e8874] | 334 | for ( k = t-1; k > 1; k-- ) |
---|
| 335 | { |
---|
[157108] | 336 | Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) ); |
---|
| 337 | n[k] = degree( Uk[k], Variable( k ) ); |
---|
[ba1fde] | 338 | } |
---|
[01e8874] | 339 | for ( k = A.min(); goodeval && (k <= t); k++ ) |
---|
| 340 | { |
---|
[157108] | 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 ); |
---|
[ba1fde] | 346 | } |
---|
| 347 | DEBDECLEVEL( cerr, "Hensel" ); |
---|
[157108] | 348 | delete[] n; |
---|
[ba1fde] | 349 | return goodeval; |
---|
| 350 | } |
---|