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