[493c477] | 1 | /* emacs edit mode for this file is -*- C++ -*- */ |
---|
[341696] | 2 | /* $Id$ */ |
---|
[a2c0302] | 3 | |
---|
[e4fe2b] | 4 | #include "config.h" |
---|
[173e86] | 5 | |
---|
[2dd068] | 6 | #include <math.h> |
---|
[57f13b] | 7 | |
---|
[650f2d8] | 8 | #include "cf_assert.h" |
---|
[57f13b] | 9 | #include "debug.h" |
---|
| 10 | #include "timing.h" |
---|
| 11 | |
---|
[2dd068] | 12 | #include "cf_defs.h" |
---|
[f58e95] | 13 | #include "cf_algorithm.h" |
---|
[2dd068] | 14 | #include "fac_util.h" |
---|
| 15 | #include "fac_univar.h" |
---|
| 16 | #include "fac_cantzass.h" |
---|
| 17 | #include "fac_berlekamp.h" |
---|
| 18 | #include "cf_iter.h" |
---|
| 19 | #include "cf_primes.h" |
---|
| 20 | #include "fac_sqrfree.h" |
---|
| 21 | |
---|
[e76d7a6] | 22 | TIMING_DEFINE_PRINT(fac_choosePrimes) |
---|
| 23 | TIMING_DEFINE_PRINT(fac_facModPrimes) |
---|
| 24 | TIMING_DEFINE_PRINT(fac_liftFactors) |
---|
| 25 | TIMING_DEFINE_PRINT(fac_combineFactors) |
---|
[57f13b] | 26 | |
---|
[a2c0302] | 27 | |
---|
| 28 | const int max_fp_fac = 3; |
---|
[2dd068] | 29 | |
---|
| 30 | static modpk theModulus; |
---|
| 31 | |
---|
[57f13b] | 32 | #ifdef DEBUGOUTPUT |
---|
[a2c0302] | 33 | #define DEBOUTHPRINT(stream, msg, hg) \ |
---|
[346edc8] | 34 | {std::stream << deb_level_msg << msg, std::stream.flush(); hprint( hg ); std::stream << std::endl;} |
---|
[a2c0302] | 35 | |
---|
[2dd068] | 36 | static void |
---|
| 37 | hprint ( int * a ) |
---|
| 38 | { |
---|
| 39 | int n = a[0]; |
---|
[346edc8] | 40 | std::cerr << "( " << n << ": "; |
---|
[2dd068] | 41 | int i = 1; |
---|
| 42 | while ( i < n ) { |
---|
[806c18] | 43 | if ( a[i] != 0 ) |
---|
[346edc8] | 44 | std::cerr << i << " "; |
---|
[806c18] | 45 | i++; |
---|
[2dd068] | 46 | } |
---|
[346edc8] | 47 | std::cerr << ")"; |
---|
[2dd068] | 48 | } |
---|
[a2c0302] | 49 | #else /* DEBUGOUTPUT */ |
---|
| 50 | #define DEBOUTHPRINT(stream, msg, hg) |
---|
| 51 | #endif /* DEBUGOUTPUT */ |
---|
[2dd068] | 52 | |
---|
| 53 | static void |
---|
| 54 | hgroup ( int * a ) |
---|
| 55 | { |
---|
| 56 | int n = a[0]; |
---|
| 57 | int i, j, k; |
---|
| 58 | for ( i = 1; i < n; i++ ) |
---|
[806c18] | 59 | if ( a[i] != 0 ) |
---|
| 60 | for ( j = 1; j <= i; j++ ) |
---|
| 61 | if ( a[j] != 0 ) |
---|
| 62 | for ( k = i; k < n; k += j ) |
---|
| 63 | a[k] = 1; |
---|
[2dd068] | 64 | } |
---|
| 65 | |
---|
| 66 | static void |
---|
| 67 | hintersect( int * a, const int * const b ) |
---|
| 68 | { |
---|
| 69 | int i, n, na = a[0], nb = b[0]; |
---|
| 70 | if ( nb < na ) |
---|
[806c18] | 71 | n = nb; |
---|
[2dd068] | 72 | else |
---|
[806c18] | 73 | n = na; |
---|
[2dd068] | 74 | for ( i = 1; i < n; i++ ) |
---|
[806c18] | 75 | if ( b[i] == 0 ) |
---|
| 76 | a[i] = 0; |
---|
[2dd068] | 77 | a[0] = n; |
---|
| 78 | } |
---|
| 79 | |
---|
| 80 | /* |
---|
| 81 | static int |
---|
| 82 | hcount ( const int * const a ) |
---|
| 83 | { |
---|
| 84 | int n = a[0], sum = 0, i; |
---|
| 85 | for ( i = 1; i < n; i++ ) |
---|
[806c18] | 86 | if ( a[i] != 0 ) sum++; |
---|
[2dd068] | 87 | return sum; |
---|
| 88 | } |
---|
| 89 | */ |
---|
| 90 | |
---|
| 91 | static void |
---|
| 92 | initHG ( int * a, const CFFList & F ) |
---|
| 93 | { |
---|
| 94 | ListIterator<CFFactor> i; |
---|
| 95 | |
---|
| 96 | int n = a[0], k; |
---|
| 97 | for ( int j = 1; j < n; j++ ) a[j] = 0; |
---|
| 98 | for ( i = F; i.hasItem(); i++ ) |
---|
[806c18] | 99 | if ( (k = i.getItem().factor().degree()) < n ) |
---|
[c1b9927] | 100 | { |
---|
[806c18] | 101 | if ( k == -1 ) { |
---|
| 102 | STICKYWARN( k == -1, "there occured an error. factory was not able to factorize\n" |
---|
| 103 | "correctly mod p. Please send the example which caused\n" |
---|
| 104 | "this error to the authors. Nonetheless we will go on with the\n" |
---|
| 105 | "calculations hoping the result will be correct. Thank you." ); |
---|
| 106 | } |
---|
| 107 | else if ( k != 0 ) |
---|
| 108 | a[k] = 1; |
---|
[c1b9927] | 109 | } |
---|
[2dd068] | 110 | } |
---|
| 111 | |
---|
| 112 | static void |
---|
| 113 | initHG ( int * a, const Array<CanonicalForm> & F ) |
---|
| 114 | { |
---|
| 115 | int i, n = a[0], m = F.size(), k; |
---|
| 116 | for ( i = 1; i < n; i++ ) a[i] = 0; |
---|
| 117 | for ( i = 1; i < m; i++ ) |
---|
[806c18] | 118 | if ( (k = F[i].degree()) < n ) |
---|
[c1b9927] | 119 | { |
---|
| 120 | if ( k == -1 ) |
---|
| 121 | { |
---|
[806c18] | 122 | STICKYWARN( k == -1, "there occured an error. factory was not able to factorize\n" |
---|
| 123 | "correctly mod p. Please send the example which caused\n" |
---|
| 124 | "this error to the authors. Nonetheless we will go on with the\n" |
---|
| 125 | "calculations hoping the result will be correct. Thank you." ); |
---|
| 126 | } |
---|
| 127 | else if ( k != 0 ) |
---|
| 128 | a[k] = 1; |
---|
[c1b9927] | 129 | } |
---|
[2dd068] | 130 | } |
---|
| 131 | |
---|
[01e8874] | 132 | static int cmpFactor( const CFFactor & a, const CFFactor & b ) |
---|
[2dd068] | 133 | { |
---|
| 134 | CFFactor A( a ), B( b ); |
---|
| 135 | return degree( A.factor() ) > degree( B.factor() ); |
---|
| 136 | } |
---|
| 137 | |
---|
[01e8874] | 138 | //static double cf2double ( const CanonicalForm & f ) |
---|
| 139 | //{ |
---|
| 140 | // CanonicalForm a = f, q, r; |
---|
| 141 | // double m = 1, res = 0; |
---|
| 142 | // if ( a.sign() < 0 ) a = -a; |
---|
| 143 | // while ( ! a.isZero() ) { |
---|
[806c18] | 144 | // divrem( a, 10, q, r ); |
---|
| 145 | // res += m * (double)(r.intval()); |
---|
| 146 | // m *= 10; |
---|
| 147 | // a = q; |
---|
[01e8874] | 148 | // } |
---|
| 149 | // if ( f.sign() < 0 ) res = -res; |
---|
| 150 | // return res; |
---|
| 151 | //} |
---|
[2dd068] | 152 | |
---|
[e1e8a9] | 153 | //{{{ static int kBound ( const CanonicalForm & f, int p ) |
---|
| 154 | //{{{ docu |
---|
| 155 | // |
---|
[19f74fe] | 156 | // kBound() - return bound of coefficients of factors of f. |
---|
| 157 | // |
---|
| 158 | // The bound is returned as an integer k such that p^k is larger |
---|
| 159 | // than all coefficients of all possible factors of f. f should |
---|
| 160 | // be an univariate polynomial over Z. |
---|
[e1e8a9] | 161 | // |
---|
[d7eca75] | 162 | // For a discussion of the formula, see the article Mignotte - |
---|
| 163 | // 'Some Usefull Bounds' in Buchberger, Collins, Loos (eds.) - |
---|
| 164 | // 'Computer Algebra: Symbolic and Algebraic Computation', 2nd |
---|
| 165 | // ed. |
---|
| 166 | // |
---|
[db84b4] | 167 | // Use by ZFactorizeUnivariate(). |
---|
| 168 | // |
---|
[e1e8a9] | 169 | //}}} |
---|
[2dd068] | 170 | static int |
---|
| 171 | kBound ( const CanonicalForm & f, int p ) |
---|
| 172 | { |
---|
[80532d7] | 173 | return (int)(f.degree() + (double)(ilog2( euclideanNorm(f)+1 ) + 1) / (double)ilog2(p)) + 1; |
---|
[2dd068] | 174 | } |
---|
[e1e8a9] | 175 | //}}} |
---|
[2dd068] | 176 | |
---|
| 177 | modpk |
---|
| 178 | getZFacModulus() |
---|
| 179 | { |
---|
| 180 | return theModulus; |
---|
| 181 | } |
---|
| 182 | |
---|
| 183 | static bool |
---|
| 184 | liftDegreeFactRec( CFArray & theFactors, CanonicalForm & F, const CanonicalForm & recip_lf, const CanonicalForm & f, const modpk & pk, int i, int d, CFFList & ZF, int exp ) |
---|
| 185 | { |
---|
| 186 | if ( i >= theFactors.size() ) |
---|
[806c18] | 187 | return false; |
---|
[2dd068] | 188 | else if ( degree( f ) + degree( theFactors[i] ) == d ) { |
---|
[806c18] | 189 | DEBOUTLN( cerr, "ldfr (f) = " << f ); |
---|
| 190 | DEBOUTLN( cerr, "ldfr (g) = " << theFactors[i] ); |
---|
| 191 | CanonicalForm g = pp( pk( recip_lf * f * theFactors[i] ) ); |
---|
| 192 | DEBOUTLN( cerr, "ldfr (pk(f*g)) = " << g ); |
---|
| 193 | CanonicalForm gg, hh; |
---|
| 194 | DEBOUTLN( cerr, "F = " << F ); |
---|
| 195 | DEBOUTLN( cerr, "g = " << g ); |
---|
| 196 | if ( divremt( F, g, gg, hh ) && hh.isZero() ) { |
---|
| 197 | ZF.append( CFFactor( g, exp ) ); |
---|
| 198 | F = gg; |
---|
| 199 | theFactors[i] = 1; |
---|
| 200 | return true; |
---|
| 201 | } |
---|
| 202 | else { |
---|
| 203 | return liftDegreeFactRec( theFactors, F, recip_lf, f, pk, i+1, d, ZF, exp ); |
---|
| 204 | } |
---|
[2dd068] | 205 | } |
---|
| 206 | else if ( degree( f ) + degree( theFactors[i] ) > d ) |
---|
[806c18] | 207 | return false; |
---|
[2dd068] | 208 | else { |
---|
[806c18] | 209 | bool ok = liftDegreeFactRec( theFactors, F, recip_lf, pk( recip_lf * f * theFactors[i] ), pk, i+1, d, ZF, exp ); |
---|
| 210 | if ( ok ) |
---|
| 211 | theFactors[i] = 1; |
---|
| 212 | else |
---|
| 213 | ok = liftDegreeFactRec( theFactors, F, recip_lf, f, pk, i+1, d, ZF, exp ); |
---|
| 214 | return ok; |
---|
[2dd068] | 215 | } |
---|
| 216 | } |
---|
| 217 | |
---|
[a2c0302] | 218 | |
---|
[aadd442] | 219 | static int choosePrimes ( int * p, const CanonicalForm & f ) |
---|
[2dd068] | 220 | { |
---|
| 221 | int ptr = 0; |
---|
| 222 | int i = 0; |
---|
| 223 | int maxp = cf_getNumPrimes(); |
---|
| 224 | int prime; |
---|
| 225 | |
---|
[aadd442] | 226 | while ( ptr < maxp && i < max_fp_fac ) |
---|
| 227 | { |
---|
[806c18] | 228 | prime = cf_getPrime( ptr ); |
---|
| 229 | if ( mod( lc( f ), prime ) != 0 ) |
---|
| 230 | { |
---|
| 231 | setCharacteristic( prime ); |
---|
| 232 | if ( isSqrFree( mapinto( f ) ) ) |
---|
| 233 | { |
---|
| 234 | p[i] = prime; |
---|
| 235 | i++; |
---|
| 236 | } |
---|
| 237 | setCharacteristic( 0 ); |
---|
| 238 | } |
---|
| 239 | ptr++; |
---|
[2dd068] | 240 | } |
---|
[a2c0302] | 241 | return ( i == max_fp_fac ); |
---|
[2dd068] | 242 | } |
---|
| 243 | |
---|
| 244 | |
---|
| 245 | static int |
---|
[5b8726d] | 246 | UnivariateQuadraticLift ( const CanonicalForm &F, const CanonicalForm & G, const CanonicalForm &H, const modpk & pk, const CanonicalForm & /*Gamma*/, CanonicalForm & gk, CanonicalForm & hk ) |
---|
[2dd068] | 247 | { |
---|
| 248 | CanonicalForm lf, f, gamma; |
---|
| 249 | CanonicalForm a, b, aa, bb, c, g, h, g1, h1, e, modulus, tmp, q, r; |
---|
| 250 | int i, j, save; |
---|
| 251 | int p = pk.getp(), k = pk.getk(); |
---|
[f5581d] | 252 | int no_iter = (int)(log( (double)k )/log(2.0)+2); |
---|
[2dd068] | 253 | int * kvals = new int[no_iter]; |
---|
| 254 | |
---|
[a2c0302] | 255 | DEBOUTLN( cerr, "quadratic lift called with p = " << p << " and k = " << k ); |
---|
[2dd068] | 256 | for ( j = 0, i = k; i > 1; i = ( i+1 ) / 2, j++ ) kvals[j] = i; |
---|
| 257 | kvals[j] = 1; |
---|
| 258 | |
---|
| 259 | save = getCharacteristic(); setCharacteristic( 0 ); |
---|
| 260 | |
---|
| 261 | lf = lc( F ); |
---|
| 262 | f = lf * F; |
---|
| 263 | { |
---|
[806c18] | 264 | setCharacteristic( p ); |
---|
| 265 | g1 = mapinto( lf ) / lc( G ) * G; |
---|
| 266 | h1 = mapinto( lf ) / lc( H ) * H; |
---|
| 267 | (void)extgcd( g1, h1, a, b ); |
---|
| 268 | setCharacteristic( 0 ); |
---|
[2dd068] | 269 | } |
---|
| 270 | a = mapinto( a ); b = mapinto( b ); |
---|
| 271 | g = mapinto( g1 ); h = mapinto( h1 ); |
---|
| 272 | g = replaceLc( g, lf ); h = replaceLc( h, lf ); |
---|
| 273 | e = f - g * h; |
---|
| 274 | modulus = p; |
---|
| 275 | i = 1; |
---|
| 276 | |
---|
| 277 | while ( ! e.isZero() && j > 0 ) { |
---|
[806c18] | 278 | c = e / modulus; |
---|
| 279 | { |
---|
| 280 | j--; |
---|
| 281 | setCharacteristic( p, kvals[j+1] ); |
---|
| 282 | DEBOUTLN( cerr, "lifting from p^" << kvals[j+1] << " to p^" << kvals[j] ); |
---|
| 283 | c = mapinto( c ); |
---|
| 284 | DEBOUTLN( cerr, " !!! g = " << mapinto( g ) ); |
---|
| 285 | g1 = mapinto( lf ) / mapinto( lc( g ) ) * mapinto( g ); |
---|
| 286 | h1 = mapinto( lf ) / mapinto( lc( h ) ) * mapinto( h ); |
---|
| 287 | // (void)extgcd( g1, h1, a, b ); |
---|
| 288 | // DEBOUTLN( cerr, " a = " << aa ); |
---|
| 289 | // DEBOUTLN( cerr, " b = " << bb ); |
---|
| 290 | a = mapinto( a ); b = mapinto( b ); |
---|
| 291 | a += ( ( 1 - a * g1 ) * a ) % h1; |
---|
| 292 | b += ( ( 1 - b * h1 ) * b ) % g1; |
---|
| 293 | DEBOUTLN( cerr, " a = " << a ); |
---|
| 294 | DEBOUTLN( cerr, " b = " << b ); |
---|
| 295 | divrem( a * c, h1, q, r ); |
---|
| 296 | tmp = b * c + q * g1; |
---|
| 297 | setCharacteristic( 0 ); |
---|
| 298 | } |
---|
| 299 | a = mapinto( a ); b = mapinto( b ); |
---|
| 300 | g += mapinto( tmp ) * modulus; |
---|
| 301 | h += mapinto( r ) * modulus; |
---|
| 302 | // g = replaceLc( g, lf ); h = replaceLc( h, lf ); |
---|
| 303 | e = f - g * h; |
---|
| 304 | modulus = power( CanonicalForm(p), kvals[j] ); |
---|
[82db91a] | 305 | if ( mod( f - g * h, modulus ) != 0 ) { |
---|
[806c18] | 306 | DEBOUTLN( cerr, "error at lift stage " << i ); |
---|
[82db91a] | 307 | } |
---|
[806c18] | 308 | i++; |
---|
[2dd068] | 309 | } |
---|
| 310 | if ( e.isZero() ) { |
---|
[806c18] | 311 | tmp = content( g ); |
---|
| 312 | gk = g / tmp; hk = h / ( lf / tmp ); |
---|
[2dd068] | 313 | } |
---|
| 314 | else { |
---|
[806c18] | 315 | gk = pk(g); hk = pk(h); |
---|
[2dd068] | 316 | } |
---|
| 317 | setCharacteristic( save ); |
---|
| 318 | return e.isZero(); |
---|
| 319 | } |
---|
| 320 | |
---|
| 321 | static int |
---|
[5b8726d] | 322 | UnivariateLinearLift ( const CanonicalForm &F, const CanonicalForm & G, const CanonicalForm &H, const modpk & pk, const CanonicalForm & /*Gamma*/, CanonicalForm & gk, CanonicalForm & hk ) |
---|
[2dd068] | 323 | { |
---|
| 324 | CanonicalForm lf, f, gamma; |
---|
| 325 | CanonicalForm a, b, c, g, h, g1, h1, e, modulus, tmp, q, r; |
---|
| 326 | int i, save; |
---|
| 327 | int p = pk.getp(), k = pk.getk(); |
---|
| 328 | save = getCharacteristic(); setCharacteristic( 0 ); |
---|
| 329 | |
---|
| 330 | lf = lc( F ); |
---|
| 331 | f = lf * F; |
---|
| 332 | { |
---|
[806c18] | 333 | setCharacteristic( p ); |
---|
| 334 | g1 = mapinto( lf ) / lc( G ) * G; |
---|
| 335 | h1 = mapinto( lf ) / lc( H ) * H; |
---|
| 336 | (void)extgcd( g1, h1, a, b ); |
---|
| 337 | setCharacteristic( 0 ); |
---|
[2dd068] | 338 | } |
---|
| 339 | g = mapinto( g1 ); h = mapinto( h1 ); |
---|
| 340 | g = replaceLc( g, lf ); h = replaceLc( h, lf ); |
---|
| 341 | e = f - g * h; |
---|
| 342 | modulus = p; |
---|
| 343 | i = 1; |
---|
| 344 | |
---|
| 345 | while ( ! e.isZero() && i <= k ) { |
---|
[806c18] | 346 | c = e / modulus; |
---|
| 347 | { |
---|
| 348 | setCharacteristic( p ); |
---|
| 349 | c = mapinto( c ); |
---|
| 350 | divrem( a * c, h1, q, r ); |
---|
| 351 | tmp = b * c + q * g1; |
---|
| 352 | setCharacteristic( 0 ); |
---|
| 353 | } |
---|
| 354 | g += mapinto( tmp ) * modulus; |
---|
| 355 | h += mapinto( r ) * modulus; |
---|
| 356 | // g = replaceLc( g, lf ); h = replaceLc( h, lf ); |
---|
| 357 | e = f - g * h; |
---|
| 358 | modulus *= p; |
---|
| 359 | ASSERT( mod( f - g * h, modulus ) == 0, "error at lift stage" ); |
---|
| 360 | i++; |
---|
[2dd068] | 361 | } |
---|
| 362 | if ( e.isZero() ) { |
---|
[806c18] | 363 | tmp = content( g ); |
---|
| 364 | gk = g / tmp; hk = h / ( lf / tmp ); |
---|
[2dd068] | 365 | } |
---|
| 366 | else { |
---|
[806c18] | 367 | gk = pk(g); hk = pk(h); |
---|
[2dd068] | 368 | } |
---|
| 369 | setCharacteristic( save ); |
---|
| 370 | // return e.isZero(); |
---|
| 371 | return (F-gk*hk).isZero(); |
---|
| 372 | } |
---|
| 373 | |
---|
[a2c0302] | 374 | |
---|
[2dd068] | 375 | CFFList |
---|
| 376 | ZFactorizeUnivariate( const CanonicalForm& ff, bool issqrfree ) |
---|
| 377 | { |
---|
| 378 | bool symmsave = isOn( SW_SYMMETRIC_FF ); |
---|
| 379 | CanonicalForm cont = content( ff ); |
---|
| 380 | CanonicalForm lf, recip_lf, fp, f, g = ff / cont, dummy1, dummy2; |
---|
| 381 | int i, k, exp, n; |
---|
| 382 | bool ok; |
---|
[a2c0302] | 383 | CFFList H, F[max_fp_fac]; |
---|
[2dd068] | 384 | CFFList ZF; |
---|
[a2c0302] | 385 | int * p = new int [max_fp_fac]; |
---|
[2dd068] | 386 | int * D = 0; |
---|
| 387 | int * Dh = 0; |
---|
| 388 | ListIterator<CFFactor> J, I; |
---|
[a2c0302] | 389 | |
---|
| 390 | DEBINCLEVEL( cerr, "ZFactorizeUnivariate" ); |
---|
[2dd068] | 391 | On( SW_SYMMETRIC_FF ); |
---|
[a2c0302] | 392 | |
---|
| 393 | // get squarefree decomposition of f |
---|
[2dd068] | 394 | if ( issqrfree ) |
---|
[806c18] | 395 | H.append( CFFactor( g, 1 ) ); |
---|
[2dd068] | 396 | else |
---|
[806c18] | 397 | H = sqrFree( g ); |
---|
[a2c0302] | 398 | |
---|
| 399 | DEBOUTLN( cerr, "H = " << H ); |
---|
| 400 | |
---|
| 401 | // cycle through squarefree factors of f |
---|
[2dd068] | 402 | for ( J = H; J.hasItem(); ++J ) { |
---|
[806c18] | 403 | f = J.getItem().factor(); |
---|
| 404 | if ( f.inCoeffDomain() ) continue; |
---|
| 405 | n = f.degree() / 2 + 1; |
---|
| 406 | delete [] D; |
---|
| 407 | delete [] Dh; |
---|
| 408 | D = new int [n]; D[0] = n; |
---|
| 409 | Dh = new int [n]; Dh[0] = n; |
---|
| 410 | exp = J.getItem().exp(); |
---|
| 411 | |
---|
| 412 | // choose primes to factor f |
---|
| 413 | TIMING_START(fac_choosePrimes); |
---|
| 414 | ok = choosePrimes( p, f ); |
---|
| 415 | TIMING_END_AND_PRINT(fac_choosePrimes, "time to choose the primes: "); |
---|
| 416 | if ( ! ok ) { |
---|
| 417 | DEBOUTLN( cerr, "warning: no good prime found to factorize " << f ); |
---|
| 418 | STICKYWARN( ok, "there occured an error. We went out of primes p\n" |
---|
| 419 | "to factorize mod p. Please send the example which caused\n" |
---|
| 420 | "this error to the authors. Nonetheless we will go on with the\n" |
---|
| 421 | "calculations hoping the result will be correct. Thank you."); |
---|
| 422 | ZF.append( CFFactor( f, exp ) ); |
---|
| 423 | continue; |
---|
| 424 | } |
---|
| 425 | |
---|
| 426 | // factorize f modulo certain primes |
---|
| 427 | TIMING_START(fac_facModPrimes); |
---|
| 428 | for ( i = 0; i < max_fp_fac; i++ ) { |
---|
| 429 | setCharacteristic( p[i] ); |
---|
| 430 | fp = mapinto( f ); |
---|
| 431 | F[i] = FpFactorizeUnivariateCZ( fp, true, 0, Variable(), Variable() ); |
---|
| 432 | // if ( p[i] < 23 && fp.degree() < 10 ) |
---|
| 433 | // F[i] = FpFactorizeUnivariateB( fp, true ); |
---|
| 434 | // else |
---|
| 435 | // F[i] = FpFactorizeUnivariateCZ( fp, true, 0, Variable, Variable() ); |
---|
| 436 | DEBOUTLN( cerr, "F[i] = " << F[i] << ", p = " << p[i] ); |
---|
| 437 | } |
---|
| 438 | TIMING_END_AND_PRINT(fac_facModPrimes, "time to factorize mod primes: "); |
---|
| 439 | setCharacteristic( 0 ); |
---|
| 440 | |
---|
| 441 | // do some strange things with the D's |
---|
| 442 | initHG( D, F[0] ); |
---|
| 443 | hgroup( D ); |
---|
| 444 | DEBOUTHPRINT( cerr, "D = ", D ); |
---|
| 445 | for ( i = 1; i < max_fp_fac; i++ ) { |
---|
| 446 | initHG( Dh, F[i] ); |
---|
| 447 | hgroup( Dh ); |
---|
| 448 | DEBOUTHPRINT( cerr, "Dh = ", Dh ); |
---|
| 449 | hintersect( D, Dh ); |
---|
| 450 | DEBOUTHPRINT( cerr, "D = ", D ); |
---|
| 451 | } |
---|
| 452 | |
---|
| 453 | // look which p gives the shortest factorization of f mod p |
---|
| 454 | // j: index of that p in p[] |
---|
| 455 | int min, j; |
---|
| 456 | min = F[0].length(), j = 0; |
---|
| 457 | for ( i = 1; i < max_fp_fac; i++ ) { |
---|
| 458 | if ( min >= F[i].length() ) { |
---|
| 459 | j = i; min = F[i].length(); |
---|
| 460 | } |
---|
| 461 | } |
---|
| 462 | k = kBound( f, p[j] ); |
---|
| 463 | CFArray theFactors( F[j].length() ); |
---|
| 464 | // pk = power( CanonicalForm( p[j] ), k ); |
---|
| 465 | // pkhalf = pk / 2; |
---|
| 466 | modpk pk( p[j], k ); |
---|
| 467 | DEBOUTLN( cerr, "coeff bound = " << pk.getpk() ); |
---|
| 468 | theModulus = pk; |
---|
| 469 | setCharacteristic( p[j] ); |
---|
| 470 | fp = mapinto( f ); |
---|
| 471 | F[j].sort( cmpFactor ); |
---|
| 472 | I = F[j]; i = 0; |
---|
| 473 | TIMING_START(fac_liftFactors); |
---|
| 474 | while ( I.hasItem() ) { |
---|
| 475 | DEBOUTLN( cerr, "factor to lift = " << I.getItem().factor() ); |
---|
| 476 | if ( isOn( SW_FAC_QUADRATICLIFT ) ) |
---|
| 477 | ok = UnivariateQuadraticLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 ); |
---|
| 478 | else |
---|
| 479 | ok = UnivariateLinearLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 ); |
---|
| 480 | if ( ok ) { |
---|
| 481 | // should be done in a more efficient way |
---|
| 482 | DEBOUTLN( cerr, "dummy1 = " << dummy1 ); |
---|
| 483 | DEBOUTLN( cerr, "dummy2 = " << dummy2 ); |
---|
| 484 | f = dummy2; |
---|
| 485 | fp /= I.getItem().factor(); |
---|
| 486 | ZF.append( CFFactor( dummy1, exp ) ); |
---|
| 487 | I.remove( 0 ); |
---|
| 488 | I = F[j]; |
---|
| 489 | i = 0; |
---|
| 490 | DEBOUTLN( cerr, "F[j] = " << F[j] ); |
---|
| 491 | } |
---|
| 492 | else { |
---|
| 493 | DEBOUTLN( cerr, "i = " << i ); |
---|
| 494 | DEBOUTLN( cerr, "dummy1 = " << dummy1 ); |
---|
| 495 | setCharacteristic( 0 ); |
---|
| 496 | // theFactors[i] = pk( dummy1 * pk.inverse( lc( dummy1 ) ) ); |
---|
| 497 | theFactors[i] = pk( dummy1 ); |
---|
| 498 | setCharacteristic( p[j] ); |
---|
| 499 | i++; |
---|
| 500 | I++; |
---|
| 501 | } |
---|
| 502 | } |
---|
| 503 | TIMING_END_AND_PRINT(fac_liftFactors, "time to lift the factors: "); |
---|
| 504 | DEBOUTLN( cerr, "ZF = " << ZF ); |
---|
| 505 | initHG( Dh, theFactors ); |
---|
| 506 | hgroup( Dh ); |
---|
| 507 | DEBOUTHPRINT( cerr, "Dh = ", Dh ); |
---|
| 508 | hintersect( D, Dh ); |
---|
| 509 | setCharacteristic( 0 ); |
---|
| 510 | for ( int l = i; l < F[j].length(); l++ ) |
---|
| 511 | theFactors[l] = 1; |
---|
| 512 | DEBOUTLN( cerr, "theFactors = " << theFactors ); |
---|
| 513 | DEBOUTLN( cerr, "f = " << f ); |
---|
| 514 | DEBOUTLN( cerr, "p = " << pk.getp() << ", k = " << pk.getk() ); |
---|
| 515 | DEBOUTHPRINT( cerr, "D = ", D ); |
---|
| 516 | lf = lc( f ); |
---|
| 517 | (void)bextgcd( pk.getpk(), lf, dummy1, recip_lf ); |
---|
| 518 | DEBOUTLN( cerr, "recip_lf = " << recip_lf ); |
---|
| 519 | TIMING_START(fac_combineFactors); |
---|
| 520 | for ( i = 1; i < D[0]; i++ ) |
---|
| 521 | if ( D[i] != 0 ) |
---|
| 522 | while ( liftDegreeFactRec( theFactors, f, recip_lf, lf, pk, 0, i, ZF, exp ) ); |
---|
| 523 | if ( degree( f ) > 0 ) |
---|
| 524 | ZF.append( CFFactor( f, exp ) ); |
---|
| 525 | TIMING_END_AND_PRINT(fac_combineFactors, "time to combine the factors: "); |
---|
[2dd068] | 526 | } |
---|
[a2c0302] | 527 | |
---|
| 528 | // brush up our result |
---|
[2dd068] | 529 | if ( ZF.getFirst().factor().inCoeffDomain() ) |
---|
[806c18] | 530 | ZF.removeFirst(); |
---|
[af4c56] | 531 | if ( lc( ff ).sign() < 0 ) |
---|
[806c18] | 532 | ZF.insert( CFFactor( -cont, 1 ) ); |
---|
[0af18dd] | 533 | else |
---|
[806c18] | 534 | ZF.insert( CFFactor( cont, 1 ) ); |
---|
[a2c0302] | 535 | delete [] D; |
---|
| 536 | delete [] Dh; |
---|
[2dd068] | 537 | if ( ! symmsave ) |
---|
[806c18] | 538 | Off( SW_SYMMETRIC_FF ); |
---|
[a2c0302] | 539 | |
---|
| 540 | DEBDECLEVEL( cerr, "ZFactorizeUnivariate" ); |
---|
[2dd068] | 541 | return ZF; |
---|
| 542 | } |
---|