[10af64] | 1 | // -*- c++ -*- |
---|
| 2 | //***************************************************************************** |
---|
| 3 | /** @file cf_gcd_smallp.cc |
---|
| 4 | * |
---|
| 5 | * @author Martin Lee |
---|
| 6 | * @date 22.10.2009 |
---|
| 7 | * |
---|
| 8 | * This file implements the GCD of two polynomials over \f$ F_{p} \f$ , |
---|
| 9 | * \f$ F_{p}(\alpha ) \f$ or GF based on Alg. 7.2. as described in "Algorithms |
---|
| 10 | * for Computer Algebra" by Geddes, Czapor, Labahnn |
---|
| 11 | * |
---|
| 12 | * @par Copyright: |
---|
| 13 | * (c) by The SINGULAR Team, see LICENSE file |
---|
| 14 | * |
---|
| 15 | * @internal |
---|
| 16 | * @version \$Id$ |
---|
| 17 | * |
---|
| 18 | **/ |
---|
| 19 | //***************************************************************************** |
---|
| 20 | |
---|
| 21 | #include <config.h> |
---|
| 22 | |
---|
| 23 | #include "assert.h" |
---|
| 24 | #include "debug.h" |
---|
| 25 | #include "timing.h" |
---|
| 26 | |
---|
| 27 | #include "canonicalform.h" |
---|
| 28 | #include "cf_map.h" |
---|
[88f3644] | 29 | #include "cf_util.h" |
---|
[10af64] | 30 | #include "ftmpl_functions.h" |
---|
| 31 | #include "cf_random.h" |
---|
| 32 | |
---|
[c4f4fd] | 33 | // iinline helper functions: |
---|
[51615d6] | 34 | #include "cf_map_ext.h" |
---|
[9c115e1] | 35 | |
---|
[10af64] | 36 | #ifdef HAVE_NTL |
---|
| 37 | #include <NTL/ZZ_pEX.h> |
---|
[04dd0c] | 38 | #include <NTLconvert.h> |
---|
[10af64] | 39 | #endif |
---|
| 40 | |
---|
[911444] | 41 | #include "cf_gcd_smallp.h" |
---|
| 42 | |
---|
[10af64] | 43 | #ifdef HAVE_NTL |
---|
| 44 | |
---|
| 45 | TIMING_DEFINE_PRINT(gcd_recursion); |
---|
| 46 | TIMING_DEFINE_PRINT(newton_interpolation); |
---|
| 47 | |
---|
| 48 | /// compressing two polynomials F and G, M is used for compressing, |
---|
| 49 | /// N to reverse the compression |
---|
[04dd0c] | 50 | static inline |
---|
[10af64] | 51 | int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M, |
---|
[c4f4fd] | 52 | CFMap & N, bool& top_level) |
---|
[10af64] | 53 | { |
---|
| 54 | int n= tmax (F.level(), G.level()); |
---|
| 55 | int * degsf= new int [n + 1]; |
---|
| 56 | int * degsg= new int [n + 1]; |
---|
| 57 | |
---|
| 58 | for (int i = 0; i <= n; i++) |
---|
| 59 | degsf[i]= degsg[i]= 0; |
---|
| 60 | |
---|
| 61 | degsf= degrees (F, degsf); |
---|
| 62 | degsg= degrees (G, degsg); |
---|
| 63 | |
---|
| 64 | int both_non_zero= 0; |
---|
| 65 | int f_zero= 0; |
---|
| 66 | int g_zero= 0; |
---|
| 67 | int both_zero= 0; |
---|
| 68 | |
---|
[c4f4fd] | 69 | if (top_level) |
---|
[10af64] | 70 | { |
---|
| 71 | for (int i= 1; i <= n; i++) |
---|
| 72 | { |
---|
| 73 | if (degsf[i] != 0 && degsg[i] != 0) |
---|
| 74 | { |
---|
| 75 | both_non_zero++; |
---|
| 76 | continue; |
---|
| 77 | } |
---|
| 78 | if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) |
---|
| 79 | { |
---|
| 80 | f_zero++; |
---|
| 81 | continue; |
---|
| 82 | } |
---|
| 83 | if (degsg[i] == 0 && degsf[i] && i <= F.level()) |
---|
| 84 | { |
---|
| 85 | g_zero++; |
---|
| 86 | continue; |
---|
| 87 | } |
---|
| 88 | } |
---|
| 89 | |
---|
| 90 | if (both_non_zero == 0) return 0; |
---|
| 91 | |
---|
| 92 | // map Variables which do not occur in both polynomials to higher levels |
---|
| 93 | int k= 1; |
---|
| 94 | int l= 1; |
---|
| 95 | for (int i= 1; i <= n; i++) |
---|
| 96 | { |
---|
| 97 | if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level()) |
---|
| 98 | { |
---|
| 99 | if (k + both_non_zero != i) |
---|
| 100 | { |
---|
| 101 | M.newpair (Variable (i), Variable (k + both_non_zero)); |
---|
| 102 | N.newpair (Variable (k + both_non_zero), Variable (i)); |
---|
| 103 | } |
---|
| 104 | k++; |
---|
| 105 | } |
---|
| 106 | if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) |
---|
| 107 | { |
---|
| 108 | if (l + g_zero + both_non_zero != i) |
---|
| 109 | { |
---|
| 110 | M.newpair (Variable (i), Variable (l + g_zero + both_non_zero)); |
---|
| 111 | N.newpair (Variable (l + g_zero + both_non_zero), Variable (i)); |
---|
| 112 | } |
---|
| 113 | l++; |
---|
| 114 | } |
---|
| 115 | } |
---|
| 116 | |
---|
| 117 | // sort Variables x_{i} in increasing order of |
---|
| 118 | // min(deg_{x_{i}}(f),deg_{x_{i}}(g)) |
---|
| 119 | int m= tmin (F.level(), G.level()); |
---|
| 120 | int max_min_deg; |
---|
| 121 | k= both_non_zero; |
---|
| 122 | l= 0; |
---|
| 123 | int i= 1; |
---|
| 124 | while (k > 0) |
---|
| 125 | { |
---|
| 126 | max_min_deg= tmin (degsf[i], degsg[i]); |
---|
| 127 | while (max_min_deg == 0) |
---|
| 128 | { |
---|
| 129 | i++; |
---|
| 130 | max_min_deg= tmin (degsf[i], degsg[i]); |
---|
| 131 | } |
---|
| 132 | for (int j= i + 1; j <= m; j++) |
---|
| 133 | { |
---|
| 134 | if (tmin (degsf[j],degsg[j]) >= max_min_deg) |
---|
| 135 | { |
---|
| 136 | max_min_deg= tmin (degsf[j], degsg[j]); |
---|
| 137 | l= j; |
---|
| 138 | } |
---|
| 139 | } |
---|
| 140 | if (l != 0) |
---|
| 141 | { |
---|
| 142 | if (l != k) |
---|
| 143 | { |
---|
| 144 | M.newpair (Variable (l), Variable(k)); |
---|
| 145 | N.newpair (Variable (k), Variable(l)); |
---|
| 146 | degsf[l]= 0; |
---|
| 147 | degsg[l]= 0; |
---|
| 148 | l= 0; |
---|
| 149 | } |
---|
| 150 | else |
---|
| 151 | { |
---|
| 152 | degsf[l]= 0; |
---|
| 153 | degsg[l]= 0; |
---|
| 154 | l= 0; |
---|
| 155 | } |
---|
| 156 | } |
---|
| 157 | else if (l == 0) |
---|
| 158 | { |
---|
| 159 | if (i != k) |
---|
| 160 | { |
---|
| 161 | M.newpair (Variable (i), Variable (k)); |
---|
| 162 | N.newpair (Variable (k), Variable (i)); |
---|
| 163 | degsf[i]= 0; |
---|
| 164 | degsg[i]= 0; |
---|
| 165 | } |
---|
| 166 | else |
---|
| 167 | { |
---|
| 168 | degsf[i]= 0; |
---|
| 169 | degsg[i]= 0; |
---|
| 170 | } |
---|
| 171 | i++; |
---|
| 172 | } |
---|
| 173 | k--; |
---|
| 174 | } |
---|
| 175 | } |
---|
| 176 | else |
---|
| 177 | { |
---|
| 178 | //arrange Variables such that no gaps occur |
---|
| 179 | for (int i= 1; i <= n; i++) |
---|
| 180 | { |
---|
| 181 | if (degsf[i] == 0 && degsg[i] == 0) |
---|
| 182 | { |
---|
| 183 | both_zero++; |
---|
| 184 | continue; |
---|
| 185 | } |
---|
| 186 | else |
---|
| 187 | { |
---|
| 188 | if (both_zero != 0) |
---|
| 189 | { |
---|
| 190 | M.newpair (Variable (i), Variable (i - both_zero)); |
---|
| 191 | N.newpair (Variable (i - both_zero), Variable (i)); |
---|
| 192 | } |
---|
| 193 | } |
---|
| 194 | } |
---|
| 195 | } |
---|
| 196 | |
---|
| 197 | delete [] degsf; |
---|
| 198 | delete [] degsg; |
---|
| 199 | |
---|
[c4f4fd] | 200 | return 1; |
---|
[10af64] | 201 | } |
---|
| 202 | |
---|
| 203 | /// compute the content of F, where F is considered as an element of |
---|
| 204 | /// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$ |
---|
| 205 | static inline CanonicalForm |
---|
| 206 | uni_content (const CanonicalForm & F) |
---|
| 207 | { |
---|
| 208 | if (F.inBaseDomain()) |
---|
| 209 | return F.genOne(); |
---|
| 210 | if (F.level() == 1 && F.isUnivariate()) |
---|
| 211 | return F; |
---|
| 212 | if (F.level() != 1 && F.isUnivariate()) |
---|
| 213 | return F.genOne(); |
---|
[c4f4fd] | 214 | if (degree (F,1) == 0) return F.genOne(); |
---|
[10af64] | 215 | |
---|
| 216 | int l= F.level(); |
---|
| 217 | if (l == 2) |
---|
| 218 | return content(F); |
---|
| 219 | else |
---|
| 220 | { |
---|
| 221 | CanonicalForm pol, c = 0; |
---|
| 222 | CFIterator i = F; |
---|
| 223 | for (; i.hasTerms(); i++) |
---|
| 224 | { |
---|
| 225 | pol= i.coeff(); |
---|
| 226 | pol= uni_content (pol); |
---|
| 227 | c= gcd (c, pol); |
---|
| 228 | if (c.isOne()) |
---|
| 229 | return c; |
---|
| 230 | } |
---|
| 231 | return c; |
---|
| 232 | } |
---|
| 233 | } |
---|
| 234 | |
---|
| 235 | /// compute the leading coefficient of F, where F is considered as an element |
---|
| 236 | /// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on Mon (x_{2},\ldots ,x_{n}) |
---|
[c4f4fd] | 237 | /// is dp. |
---|
[10af64] | 238 | static inline |
---|
| 239 | CanonicalForm uni_lcoeff (const CanonicalForm& F) |
---|
| 240 | { |
---|
| 241 | if (F.level() <= 1) |
---|
| 242 | return F; |
---|
| 243 | else |
---|
| 244 | { |
---|
| 245 | Variable x= Variable (2); |
---|
| 246 | int deg= totaldegree (F, x, F.mvar()); |
---|
| 247 | for (CFIterator i= F; i.hasTerms(); i++) |
---|
| 248 | { |
---|
| 249 | if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg) |
---|
| 250 | return uni_lcoeff (i.coeff()); |
---|
| 251 | } |
---|
| 252 | } |
---|
| 253 | } |
---|
| 254 | |
---|
| 255 | /// Newton interpolation - Incremental algorithm. |
---|
| 256 | /// Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n, |
---|
| 257 | /// computes the interpolation polynomial assuming that |
---|
| 258 | /// the polynomials in u are the results of evaluating the variabe x |
---|
| 259 | /// of the unknown polynomial at the alpha values. |
---|
| 260 | /// This incremental version receives only the values of alpha_n and u_n and |
---|
| 261 | /// the previous interpolation polynomial for points 1 <= i <= n-1, and computes |
---|
| 262 | /// the polynomial interpolating in all the points. |
---|
| 263 | /// newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1}) |
---|
| 264 | static inline CanonicalForm |
---|
| 265 | newtonInterp(const CanonicalForm alpha, const CanonicalForm u, const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly, const Variable & x) |
---|
| 266 | { |
---|
| 267 | CanonicalForm interPoly; |
---|
| 268 | |
---|
| 269 | interPoly = oldInterPoly + ((u - oldInterPoly(alpha, x)) / newtonPoly(alpha, x)) * newtonPoly; |
---|
| 270 | return interPoly; |
---|
| 271 | } |
---|
| 272 | |
---|
| 273 | /// compute a random element a of \f$ F_{p}(\alpha ) \f$ , |
---|
| 274 | /// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns |
---|
| 275 | /// fail if there are no field elements left which have not been used before |
---|
| 276 | static inline CanonicalForm |
---|
| 277 | randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list, |
---|
| 278 | bool & fail) |
---|
| 279 | { |
---|
| 280 | fail= false; |
---|
| 281 | Variable x= F.mvar(); |
---|
| 282 | AlgExtRandomF genAlgExt (alpha); |
---|
| 283 | FFRandom genFF; |
---|
| 284 | CanonicalForm random, mipo; |
---|
| 285 | mipo= getMipo (alpha); |
---|
| 286 | int p= getCharacteristic (); |
---|
| 287 | int d= degree (mipo); |
---|
[c4f4fd] | 288 | int bound= ipower (p, d); |
---|
[10af64] | 289 | do |
---|
| 290 | { |
---|
| 291 | if (list.length() == bound) |
---|
| 292 | { |
---|
| 293 | fail= true; |
---|
| 294 | break; |
---|
| 295 | } |
---|
| 296 | if (list.length() < p) |
---|
| 297 | { |
---|
| 298 | random= genFF.generate(); |
---|
| 299 | while (find (list, random)) |
---|
| 300 | random= genFF.generate(); |
---|
| 301 | } |
---|
| 302 | else |
---|
| 303 | { |
---|
| 304 | random= genAlgExt.generate(); |
---|
| 305 | while (find (list, random)) |
---|
| 306 | random= genAlgExt.generate(); |
---|
| 307 | } |
---|
| 308 | if (F (random, x) == 0) |
---|
| 309 | { |
---|
| 310 | list.append (random); |
---|
| 311 | continue; |
---|
| 312 | } |
---|
| 313 | } while (find (list, random)); |
---|
| 314 | return random; |
---|
| 315 | } |
---|
| 316 | |
---|
| 317 | /// chooses a suitable extension of \f$ F_{p}(\alpha ) \f$ |
---|
| 318 | /// we do not extend \f$ F_{p}(\alpha ) \f$ itself but extend \f$ F_{p} \f$ , |
---|
| 319 | /// s.t. \f$ F_{p}(\alpha ) \subset F_{p}(\beta ) \f$ |
---|
| 320 | static inline |
---|
| 321 | void choose_extension (const int& d, const int& num_vars, Variable& beta) |
---|
| 322 | { |
---|
| 323 | int p= getCharacteristic(); |
---|
| 324 | ZZ NTLp= to_ZZ (p); |
---|
| 325 | ZZ_p::init (NTLp); |
---|
| 326 | ZZ_pX NTLirredpoly; |
---|
| 327 | //TODO: replace d by max_{i} (deg_x{i}(f)) |
---|
[c4f4fd] | 328 | int i= (int) (log ((double) ipower (d + 1, num_vars))/log ((double) p)); |
---|
[10af64] | 329 | int m= degree (getMipo (beta)); |
---|
| 330 | if (i <= 1) |
---|
| 331 | i= 2; |
---|
| 332 | BuildIrred (NTLirredpoly, i*m); |
---|
| 333 | CanonicalForm mipo= convertNTLZZpX2CF (NTLirredpoly, Variable(1)); |
---|
| 334 | beta= rootOf (mipo); |
---|
| 335 | } |
---|
| 336 | |
---|
| 337 | /// GCD of F and G over \f$ F_{p}(\alpha ) \f$ , |
---|
[c4f4fd] | 338 | /// l and top_level are only used internally, output is monic |
---|
[10af64] | 339 | /// based on Alg. 7.2. as described in "Algorithms for |
---|
| 340 | /// Computer Algebra" by Geddes, Czapor, Labahn |
---|
[911444] | 341 | CanonicalForm |
---|
[10af64] | 342 | GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, |
---|
[c4f4fd] | 343 | Variable & alpha, CFList& l, bool& top_level) |
---|
[10af64] | 344 | { |
---|
| 345 | CanonicalForm A= F; |
---|
| 346 | CanonicalForm B= G; |
---|
| 347 | if (F.isZero() && degree(G) > 0) return G/Lc(G); |
---|
| 348 | else if (G.isZero() && degree (F) > 0) return F/Lc(F); |
---|
| 349 | else if (F.isZero() && G.isZero()) return F.genOne(); |
---|
| 350 | if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne(); |
---|
| 351 | if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F); |
---|
| 352 | if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G); |
---|
| 353 | if (F == G) return F/Lc(F); |
---|
| 354 | |
---|
| 355 | CFMap M,N; |
---|
[c4f4fd] | 356 | int best_level= myCompress (A, B, M, N, top_level); |
---|
[10af64] | 357 | |
---|
| 358 | if (best_level == 0) return B.genOne(); |
---|
| 359 | |
---|
| 360 | A= M(A); |
---|
| 361 | B= M(B); |
---|
| 362 | |
---|
| 363 | Variable x= Variable(1); |
---|
| 364 | |
---|
| 365 | //univariate case |
---|
| 366 | if (A.isUnivariate() && B.isUnivariate()) |
---|
| 367 | return N (gcd(A,B)); |
---|
| 368 | |
---|
| 369 | CanonicalForm cA, cB; // content of A and B |
---|
| 370 | CanonicalForm ppA, ppB; // primitive part of A and B |
---|
| 371 | CanonicalForm gcdcAcB; |
---|
[c4f4fd] | 372 | |
---|
| 373 | cA = uni_content (A); |
---|
| 374 | cB = uni_content (B); |
---|
| 375 | |
---|
| 376 | gcdcAcB= gcd (cA, cB); |
---|
| 377 | |
---|
| 378 | ppA= A/cA; |
---|
| 379 | ppB= B/cB; |
---|
[10af64] | 380 | |
---|
| 381 | CanonicalForm lcA, lcB; // leading coefficients of A and B |
---|
| 382 | CanonicalForm gcdlcAlcB; |
---|
| 383 | |
---|
| 384 | lcA= uni_lcoeff (ppA); |
---|
| 385 | lcB= uni_lcoeff (ppB); |
---|
| 386 | |
---|
| 387 | if (fdivides (lcA, lcB)) { |
---|
| 388 | if (fdivides (A, B)) |
---|
| 389 | return F/Lc(F); |
---|
| 390 | } |
---|
| 391 | if (fdivides (lcB, lcA)) |
---|
| 392 | { |
---|
| 393 | if (fdivides (B, A)) |
---|
| 394 | return G/Lc(G); |
---|
| 395 | } |
---|
| 396 | |
---|
| 397 | gcdlcAlcB= gcd (lcA, lcB); |
---|
| 398 | |
---|
| 399 | int d= totaldegree (ppA, Variable(2), Variable (ppA.level())); |
---|
| 400 | |
---|
[c4f4fd] | 401 | if (d == 0) return N(gcdcAcB); |
---|
[10af64] | 402 | int d0= totaldegree (ppB, Variable(2), Variable (ppB.level())); |
---|
| 403 | if (d0 < d) |
---|
| 404 | d= d0; |
---|
[c4f4fd] | 405 | if (d == 0) return N(gcdcAcB); |
---|
[10af64] | 406 | |
---|
| 407 | CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; |
---|
| 408 | CanonicalForm newtonPoly; |
---|
| 409 | |
---|
| 410 | newtonPoly= 1; |
---|
| 411 | m= gcdlcAlcB; |
---|
| 412 | G_m= 0; |
---|
| 413 | H= 0; |
---|
| 414 | bool fail= false; |
---|
[c4f4fd] | 415 | top_level= false; |
---|
[10af64] | 416 | bool inextension= false; |
---|
| 417 | Variable V_buf= alpha; |
---|
| 418 | CanonicalForm prim_elem, im_prim_elem; |
---|
| 419 | CFList source, dest; |
---|
| 420 | do |
---|
| 421 | { |
---|
| 422 | random_element= randomElement (m, V_buf, l, fail); |
---|
| 423 | if (fail) |
---|
| 424 | { |
---|
| 425 | source= CFList(); |
---|
| 426 | dest= CFList(); |
---|
| 427 | int num_vars= tmin (getNumVars(A), getNumVars(B)); |
---|
| 428 | num_vars--; |
---|
[c4f4fd] | 429 | |
---|
[10af64] | 430 | choose_extension (d, num_vars, V_buf); |
---|
| 431 | bool prim_fail= false; |
---|
| 432 | Variable V_buf2; |
---|
| 433 | prim_elem= primitiveElement (alpha, V_buf2, prim_fail); |
---|
[c4f4fd] | 434 | |
---|
[10af64] | 435 | ASSERT (!prim_fail, "failure in integer factorizer"); |
---|
| 436 | if (prim_fail) |
---|
| 437 | ; //ERROR |
---|
| 438 | else |
---|
| 439 | im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf); |
---|
[c4f4fd] | 440 | |
---|
[10af64] | 441 | DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha)); |
---|
[04dd0c] | 442 | DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2)); |
---|
[10af64] | 443 | inextension= true; |
---|
| 444 | for (CFListIterator i= l; i.hasItem(); i++) |
---|
| 445 | i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, |
---|
| 446 | im_prim_elem, source, dest); |
---|
| 447 | m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); |
---|
| 448 | G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); |
---|
| 449 | newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, |
---|
| 450 | source, dest); |
---|
| 451 | ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest); |
---|
| 452 | ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest); |
---|
| 453 | gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, |
---|
| 454 | source, dest); |
---|
| 455 | |
---|
| 456 | fail= false; |
---|
| 457 | random_element= randomElement (m, V_buf, l, fail ); |
---|
[a5cc7b3] | 458 | DEBOUTLN (cerr, "fail= " << fail); |
---|
[10af64] | 459 | CFList list; |
---|
| 460 | TIMING_START (gcd_recursion); |
---|
| 461 | G_random_element= |
---|
| 462 | GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, |
---|
[c4f4fd] | 463 | list, top_level); |
---|
[10af64] | 464 | TIMING_END_AND_PRINT (gcd_recursion, |
---|
| 465 | "time for recursive call: "); |
---|
| 466 | DEBOUTLN (cerr, "G_random_element= " << G_random_element); |
---|
| 467 | } |
---|
| 468 | else |
---|
| 469 | { |
---|
| 470 | CFList list; |
---|
| 471 | TIMING_START (gcd_recursion); |
---|
| 472 | G_random_element= |
---|
| 473 | GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf, |
---|
[c4f4fd] | 474 | list, top_level); |
---|
[10af64] | 475 | TIMING_END_AND_PRINT (gcd_recursion, |
---|
| 476 | "time for recursive call: "); |
---|
| 477 | DEBOUTLN (cerr, "G_random_element= " << G_random_element); |
---|
| 478 | } |
---|
| 479 | |
---|
| 480 | d0= totaldegree (G_random_element, Variable(2), |
---|
| 481 | Variable (G_random_element.level())); |
---|
[c4f4fd] | 482 | if (d0 == 0) return N(gcdcAcB); |
---|
[10af64] | 483 | if (d0 > d) |
---|
| 484 | { |
---|
| 485 | if (!find (l, random_element)) |
---|
| 486 | l.append (random_element); |
---|
| 487 | continue; |
---|
| 488 | } |
---|
| 489 | |
---|
| 490 | G_random_element= |
---|
| 491 | (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element)) |
---|
| 492 | * G_random_element; |
---|
| 493 | |
---|
| 494 | d0= totaldegree (G_random_element, Variable(2), |
---|
| 495 | Variable(G_random_element.level())); |
---|
| 496 | if (d0 < d) |
---|
| 497 | { |
---|
| 498 | m= gcdlcAlcB; |
---|
| 499 | newtonPoly= 1; |
---|
| 500 | G_m= 0; |
---|
| 501 | d= d0; |
---|
| 502 | } |
---|
| 503 | |
---|
| 504 | TIMING_START (newton_interpolation); |
---|
| 505 | H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x); |
---|
| 506 | TIMING_END_AND_PRINT (newton_interpolation, |
---|
| 507 | "time for newton interpolation: "); |
---|
| 508 | |
---|
| 509 | //termination test |
---|
| 510 | if (uni_lcoeff (H) == gcdlcAlcB) |
---|
| 511 | { |
---|
| 512 | cH= uni_content (H); |
---|
| 513 | ppH= H/cH; |
---|
| 514 | if (inextension) |
---|
| 515 | { |
---|
| 516 | CFList u, v; |
---|
| 517 | //maybe it's better to test if ppH is an element of F(\alpha) before |
---|
| 518 | //mapping down |
---|
| 519 | DEBOUTLN (cerr, "ppH before mapDown= " << ppH); |
---|
| 520 | ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v); |
---|
| 521 | ppH /= Lc(ppH); |
---|
| 522 | DEBOUTLN (cerr, "ppH after mapDown= " << ppH); |
---|
[c4f4fd] | 523 | if (fdivides (ppH, A) && fdivides (ppH, B)) |
---|
[10af64] | 524 | return N(gcdcAcB*ppH); |
---|
| 525 | } |
---|
[c4f4fd] | 526 | else if (fdivides (ppH, A) && fdivides (ppH, B)) |
---|
| 527 | return N(gcdcAcB*(ppH/Lc(ppH))); |
---|
[10af64] | 528 | } |
---|
| 529 | |
---|
| 530 | G_m= H; |
---|
| 531 | newtonPoly= newtonPoly*(x - random_element); |
---|
| 532 | m= m*(x - random_element); |
---|
| 533 | if (!find (l, random_element)) |
---|
| 534 | l.append (random_element); |
---|
| 535 | } while (1); |
---|
| 536 | } |
---|
| 537 | |
---|
| 538 | /// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a |
---|
| 539 | /// univariate polynomial, returns fail if there are no field elements left |
---|
| 540 | /// which have not been used before |
---|
[04dd0c] | 541 | static inline |
---|
[10af64] | 542 | CanonicalForm |
---|
| 543 | GFRandomElement (const CanonicalForm& F, CFList& list, bool& fail) |
---|
| 544 | { |
---|
| 545 | fail= false; |
---|
| 546 | Variable x= F.mvar(); |
---|
| 547 | GFRandom genGF; |
---|
| 548 | CanonicalForm random; |
---|
| 549 | int p= getCharacteristic(); |
---|
| 550 | int d= getGFDegree(); |
---|
[c4f4fd] | 551 | int bound= ipower (p, d); |
---|
[10af64] | 552 | do |
---|
| 553 | { |
---|
| 554 | if (list.length() == bound) |
---|
| 555 | { |
---|
| 556 | fail= true; |
---|
| 557 | break; |
---|
| 558 | } |
---|
| 559 | if (list.length() < 1) |
---|
| 560 | random= 0; |
---|
| 561 | else |
---|
| 562 | { |
---|
| 563 | random= genGF.generate(); |
---|
| 564 | while (find (list, random)) |
---|
| 565 | random= genGF.generate(); |
---|
| 566 | } |
---|
| 567 | if (F (random, x) == 0) |
---|
| 568 | { |
---|
| 569 | list.append (random); |
---|
| 570 | continue; |
---|
| 571 | } |
---|
| 572 | } while (find (list, random)); |
---|
| 573 | return random; |
---|
| 574 | } |
---|
| 575 | |
---|
| 576 | /// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for |
---|
| 577 | /// Computer Algebra" by Geddes, Czapor, Labahn |
---|
| 578 | /// Usually this algorithm will be faster than GCD_Fp_extension since GF has |
---|
| 579 | /// faster field arithmetics, however it might fail if the input is large since |
---|
| 580 | /// the size of the base field is bounded by 2^16, output is monic |
---|
[911444] | 581 | CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G, |
---|
[c4f4fd] | 582 | CFList& l, bool& top_level) |
---|
[10af64] | 583 | { |
---|
| 584 | CanonicalForm A= F; |
---|
| 585 | CanonicalForm B= G; |
---|
| 586 | if (F.isZero() && degree(G) > 0) return G/Lc(G); |
---|
| 587 | else if (G.isZero() && degree (F) > 0) return F/Lc(F); |
---|
| 588 | else if (F.isZero() && G.isZero()) return F.genOne(); |
---|
| 589 | if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne(); |
---|
| 590 | if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F); |
---|
| 591 | if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G); |
---|
| 592 | if (F == G) return F/Lc(F); |
---|
| 593 | |
---|
| 594 | CFMap M,N; |
---|
[c4f4fd] | 595 | int best_level= myCompress (A, B, M, N, top_level); |
---|
[10af64] | 596 | |
---|
| 597 | if (best_level == 0) return B.genOne(); |
---|
| 598 | |
---|
| 599 | A= M(A); |
---|
| 600 | B= M(B); |
---|
| 601 | |
---|
| 602 | Variable x= Variable(1); |
---|
| 603 | |
---|
| 604 | //univariate case |
---|
| 605 | if (A.isUnivariate() && B.isUnivariate()) |
---|
| 606 | return N (gcd (A, B)); |
---|
| 607 | |
---|
| 608 | CanonicalForm cA, cB; // content of A and B |
---|
| 609 | CanonicalForm ppA, ppB; // primitive part of A and B |
---|
| 610 | CanonicalForm gcdcAcB; |
---|
| 611 | |
---|
| 612 | cA = uni_content (A); |
---|
| 613 | cB = uni_content (B); |
---|
| 614 | |
---|
| 615 | gcdcAcB= gcd (cA, cB); |
---|
| 616 | |
---|
| 617 | ppA= A/cA; |
---|
| 618 | ppB= B/cB; |
---|
| 619 | |
---|
| 620 | CanonicalForm lcA, lcB; // leading coefficients of A and B |
---|
| 621 | CanonicalForm gcdlcAlcB; |
---|
| 622 | |
---|
| 623 | lcA= uni_lcoeff (ppA); |
---|
| 624 | lcB= uni_lcoeff (ppB); |
---|
| 625 | |
---|
| 626 | if (fdivides (lcA, lcB)) |
---|
| 627 | { |
---|
| 628 | if (fdivides (A, B)) |
---|
| 629 | return F; |
---|
| 630 | } |
---|
| 631 | if (fdivides (lcB, lcA)) |
---|
| 632 | { |
---|
| 633 | if (fdivides (B, A)) |
---|
| 634 | return G; |
---|
| 635 | } |
---|
| 636 | |
---|
| 637 | gcdlcAlcB= gcd (lcA, lcB); |
---|
| 638 | |
---|
| 639 | int d= totaldegree (ppA, Variable(2), Variable (ppA.level())); |
---|
| 640 | if (d == 0) return N(gcdcAcB); |
---|
| 641 | int d0= totaldegree (ppB, Variable(2), Variable (ppB.level())); |
---|
| 642 | if (d0 < d) |
---|
| 643 | d= d0; |
---|
| 644 | if (d == 0) return N(gcdcAcB); |
---|
| 645 | |
---|
| 646 | CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; |
---|
| 647 | CanonicalForm newtonPoly; |
---|
| 648 | |
---|
| 649 | newtonPoly= 1; |
---|
| 650 | m= gcdlcAlcB; |
---|
| 651 | G_m= 0; |
---|
| 652 | H= 0; |
---|
| 653 | bool fail= false; |
---|
[c4f4fd] | 654 | top_level= false; |
---|
[10af64] | 655 | bool inextension= false; |
---|
| 656 | int p; |
---|
| 657 | int k= getGFDegree(); |
---|
| 658 | int kk; |
---|
[88f3644] | 659 | int expon; |
---|
[10af64] | 660 | char gf_name_buf= gf_name; |
---|
| 661 | do |
---|
| 662 | { |
---|
| 663 | random_element= GFRandomElement (m, l, fail); |
---|
| 664 | if (fail) |
---|
| 665 | { |
---|
| 666 | int num_vars= tmin (getNumVars(A), getNumVars(B)); |
---|
| 667 | num_vars--; |
---|
| 668 | p= getCharacteristic(); |
---|
[c4f4fd] | 669 | expon= (int) floor ((log ((double) ipower (d + 1, num_vars))/log ((double) p))); |
---|
[10af64] | 670 | if (expon < 2) |
---|
| 671 | expon= 2; |
---|
| 672 | kk= getGFDegree(); |
---|
[c4f4fd] | 673 | if (ipower (p, kk*expon) < (1 << 16)) |
---|
[10af64] | 674 | setCharacteristic (p, kk*(int)expon, 'b'); |
---|
| 675 | else |
---|
| 676 | { |
---|
[04dd0c] | 677 | expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk)); |
---|
[10af64] | 678 | ASSERT (expon >= 2, "not enough points in GCD_GF"); |
---|
| 679 | setCharacteristic (p, (int)(kk*expon), 'b'); |
---|
| 680 | } |
---|
| 681 | inextension= true; |
---|
| 682 | fail= false; |
---|
| 683 | for (CFListIterator i= l; i.hasItem(); i++) |
---|
| 684 | i.getItem()= GFMapUp (i.getItem(), kk); |
---|
| 685 | m= GFMapUp (m, kk); |
---|
| 686 | G_m= GFMapUp (G_m, kk); |
---|
| 687 | newtonPoly= GFMapUp (newtonPoly, kk); |
---|
| 688 | ppA= GFMapUp (ppA, kk); |
---|
| 689 | ppB= GFMapUp (ppB, kk); |
---|
| 690 | gcdlcAlcB= GFMapUp (gcdlcAlcB, kk); |
---|
| 691 | random_element= GFRandomElement (m, l, fail); |
---|
[a5cc7b3] | 692 | DEBOUTLN (cerr, "fail= " << fail); |
---|
[10af64] | 693 | CFList list; |
---|
| 694 | TIMING_START (gcd_recursion); |
---|
| 695 | G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), |
---|
[c4f4fd] | 696 | list, top_level); |
---|
[10af64] | 697 | TIMING_END_AND_PRINT (gcd_recursion, |
---|
| 698 | "time for recursive call: "); |
---|
| 699 | DEBOUTLN (cerr, "G_random_element= " << G_random_element); |
---|
| 700 | } |
---|
| 701 | else |
---|
| 702 | { |
---|
| 703 | CFList list; |
---|
| 704 | TIMING_START (gcd_recursion); |
---|
| 705 | G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), |
---|
[c4f4fd] | 706 | list, top_level); |
---|
[10af64] | 707 | TIMING_END_AND_PRINT (gcd_recursion, |
---|
| 708 | "time for recursive call: "); |
---|
| 709 | DEBOUTLN (cerr, "G_random_element= " << G_random_element); |
---|
| 710 | } |
---|
| 711 | |
---|
| 712 | d0= totaldegree (G_random_element, Variable(2), |
---|
| 713 | Variable (G_random_element.level())); |
---|
| 714 | if (d0 == 0) |
---|
| 715 | { |
---|
| 716 | if (inextension) |
---|
| 717 | { |
---|
| 718 | setCharacteristic (p, k, gf_name_buf); |
---|
| 719 | return N(gcdcAcB); |
---|
| 720 | } |
---|
| 721 | else |
---|
| 722 | return N(gcdcAcB); |
---|
| 723 | } |
---|
| 724 | if (d0 > d) |
---|
| 725 | { |
---|
| 726 | if (!find (l, random_element)) |
---|
| 727 | l.append (random_element); |
---|
| 728 | continue; |
---|
| 729 | } |
---|
| 730 | |
---|
| 731 | G_random_element= |
---|
| 732 | (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) * |
---|
| 733 | G_random_element; |
---|
| 734 | d0= totaldegree (G_random_element, Variable(2), |
---|
| 735 | Variable (G_random_element.level())); |
---|
| 736 | |
---|
| 737 | if (d0 < d) |
---|
| 738 | { |
---|
| 739 | m= gcdlcAlcB; |
---|
| 740 | newtonPoly= 1; |
---|
| 741 | G_m= 0; |
---|
| 742 | d= d0; |
---|
| 743 | } |
---|
| 744 | |
---|
| 745 | TIMING_START (newton_interpolation); |
---|
| 746 | H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x); |
---|
| 747 | TIMING_END_AND_PRINT (newton_interpolation, "time for newton interpolation: "); |
---|
| 748 | |
---|
| 749 | //termination test |
---|
| 750 | if (uni_lcoeff (H) == gcdlcAlcB) |
---|
| 751 | { |
---|
| 752 | cH= uni_content (H); |
---|
| 753 | ppH= H/cH; |
---|
| 754 | if (inextension) |
---|
| 755 | { |
---|
| 756 | if (fdivides(ppH, GFMapUp(A, k)) && fdivides(ppH, GFMapUp(B,k))) |
---|
| 757 | { |
---|
| 758 | DEBOUTLN (cerr, "ppH before mapDown= " << ppH); |
---|
| 759 | ppH= GFMapDown (ppH, k); |
---|
| 760 | DEBOUTLN (cerr, "ppH after mapDown= " << ppH); |
---|
| 761 | setCharacteristic (p, k, gf_name_buf); |
---|
| 762 | return N(gcdcAcB*ppH); |
---|
| 763 | } |
---|
| 764 | } |
---|
| 765 | else |
---|
| 766 | { |
---|
| 767 | if (fdivides (ppH, A) && fdivides (ppH, B)) |
---|
| 768 | return N(gcdcAcB*ppH); |
---|
| 769 | } |
---|
| 770 | } |
---|
| 771 | |
---|
| 772 | G_m= H; |
---|
| 773 | newtonPoly= newtonPoly*(x - random_element); |
---|
| 774 | m= m*(x - random_element); |
---|
| 775 | if (!find (l, random_element)) |
---|
| 776 | l.append (random_element); |
---|
| 777 | } while (1); |
---|
| 778 | } |
---|
| 779 | |
---|
| 780 | /// F is assumed to be an univariate polynomial in x, |
---|
| 781 | /// computes a random monic irreducible univariate polynomial of random |
---|
| 782 | /// degree < i in x which does not divide F |
---|
| 783 | CanonicalForm |
---|
[04dd0c] | 784 | randomIrredpoly (int i, const Variable & x) |
---|
[10af64] | 785 | { |
---|
| 786 | int p= getCharacteristic(); |
---|
| 787 | ZZ NTLp= to_ZZ (p); |
---|
| 788 | ZZ_p::init (NTLp); |
---|
| 789 | ZZ_pX NTLirredpoly; |
---|
| 790 | CanonicalForm CFirredpoly; |
---|
[04dd0c] | 791 | BuildIrred (NTLirredpoly, i + 1); |
---|
| 792 | CFirredpoly= convertNTLZZpX2CF (NTLirredpoly, x); |
---|
[10af64] | 793 | return CFirredpoly; |
---|
| 794 | } |
---|
| 795 | |
---|
[04dd0c] | 796 | static inline |
---|
[10af64] | 797 | CanonicalForm |
---|
| 798 | FpRandomElement (const CanonicalForm& F, CFList& list, bool& fail) |
---|
| 799 | { |
---|
| 800 | fail= false; |
---|
| 801 | Variable x= F.mvar(); |
---|
| 802 | FFRandom genFF; |
---|
| 803 | CanonicalForm random; |
---|
| 804 | int p= getCharacteristic(); |
---|
[88f3644] | 805 | int bound= p; |
---|
[10af64] | 806 | do |
---|
| 807 | { |
---|
| 808 | if (list.length() == bound) |
---|
| 809 | { |
---|
| 810 | fail= true; |
---|
| 811 | break; |
---|
| 812 | } |
---|
| 813 | if (list.length() < 1) |
---|
| 814 | random= 0; |
---|
| 815 | else |
---|
| 816 | { |
---|
| 817 | random= genFF.generate(); |
---|
| 818 | while (find (list, random)) |
---|
| 819 | random= genFF.generate(); |
---|
| 820 | } |
---|
| 821 | if (F (random, x) == 0) |
---|
| 822 | { |
---|
| 823 | list.append (random); |
---|
| 824 | continue; |
---|
| 825 | } |
---|
| 826 | } while (find (list, random)); |
---|
| 827 | return random; |
---|
| 828 | } |
---|
| 829 | |
---|
| 830 | CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, |
---|
[c4f4fd] | 831 | bool& top_level, CFList& l) |
---|
[10af64] | 832 | { |
---|
| 833 | CanonicalForm A= F; |
---|
| 834 | CanonicalForm B= G; |
---|
| 835 | if (F.isZero() && degree(G) > 0) return G/Lc(G); |
---|
| 836 | else if (G.isZero() && degree (F) > 0) return F/Lc(F); |
---|
| 837 | else if (F.isZero() && G.isZero()) return F.genOne(); |
---|
| 838 | if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne(); |
---|
| 839 | if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F); |
---|
| 840 | if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G); |
---|
| 841 | if (F == G) return F/Lc(F); |
---|
| 842 | |
---|
| 843 | CFMap M,N; |
---|
[c4f4fd] | 844 | int best_level= myCompress (A, B, M, N, top_level); |
---|
[10af64] | 845 | |
---|
| 846 | if (best_level == 0) return B.genOne(); |
---|
| 847 | |
---|
| 848 | A= M(A); |
---|
| 849 | B= M(B); |
---|
| 850 | |
---|
[c4f4fd] | 851 | Variable x= Variable (1); |
---|
| 852 | |
---|
[10af64] | 853 | //univariate case |
---|
| 854 | if (A.isUnivariate() && B.isUnivariate()) |
---|
| 855 | return N (gcd (A, B)); |
---|
| 856 | |
---|
| 857 | CanonicalForm cA, cB; // content of A and B |
---|
| 858 | CanonicalForm ppA, ppB; // primitive part of A and B |
---|
| 859 | CanonicalForm gcdcAcB; |
---|
[c4f4fd] | 860 | cA = uni_content (A); |
---|
| 861 | cB = uni_content (B); |
---|
| 862 | gcdcAcB= gcd (cA, cB); |
---|
| 863 | ppA= A/cA; |
---|
| 864 | ppB= B/cB; |
---|
[10af64] | 865 | |
---|
| 866 | CanonicalForm lcA, lcB; // leading coefficients of A and B |
---|
| 867 | CanonicalForm gcdlcAlcB; |
---|
| 868 | lcA= uni_lcoeff (ppA); |
---|
| 869 | lcB= uni_lcoeff (ppB); |
---|
| 870 | |
---|
| 871 | if (fdivides (lcA, lcB)) |
---|
| 872 | { |
---|
| 873 | if (fdivides (A, B)) |
---|
| 874 | return F/Lc(F); |
---|
| 875 | } |
---|
| 876 | if (fdivides (lcB, lcA)) |
---|
| 877 | { |
---|
| 878 | if (fdivides (B, A)) |
---|
| 879 | return G/Lc(G); |
---|
| 880 | } |
---|
[c4f4fd] | 881 | |
---|
[10af64] | 882 | gcdlcAlcB= gcd (lcA, lcB); |
---|
| 883 | |
---|
| 884 | int d= totaldegree (ppA, Variable (2), Variable (ppA.level())); |
---|
| 885 | int d0; |
---|
| 886 | |
---|
[c4f4fd] | 887 | if (d == 0) return N(gcdcAcB); |
---|
[10af64] | 888 | d0= totaldegree (ppB, Variable (2), Variable (ppB.level())); |
---|
| 889 | |
---|
| 890 | if (d0 < d) |
---|
| 891 | d= d0; |
---|
| 892 | |
---|
[c4f4fd] | 893 | if (d == 0) return N(gcdcAcB); |
---|
[10af64] | 894 | |
---|
| 895 | CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; |
---|
| 896 | CanonicalForm newtonPoly= 1; |
---|
| 897 | m= gcdlcAlcB; |
---|
| 898 | H= 0; |
---|
| 899 | G_m= 0; |
---|
| 900 | Variable alpha, V_buf; |
---|
[c4f4fd] | 901 | int p, k; |
---|
[10af64] | 902 | bool fail= false; |
---|
| 903 | bool inextension= false; |
---|
| 904 | bool inextensionextension= false; |
---|
[c4f4fd] | 905 | top_level= false; |
---|
| 906 | CanonicalForm CF_buf; |
---|
[10af64] | 907 | CFList source, dest; |
---|
[c4f4fd] | 908 | CanonicalForm gcdcheck; |
---|
[10af64] | 909 | do |
---|
| 910 | { |
---|
| 911 | if (inextension) |
---|
| 912 | random_element= randomElement (m, alpha, l, fail); |
---|
| 913 | else |
---|
| 914 | random_element= FpRandomElement (m, l, fail); |
---|
| 915 | |
---|
| 916 | if (!fail && !inextension) |
---|
| 917 | { |
---|
| 918 | CFList list; |
---|
| 919 | TIMING_START (gcd_recursion); |
---|
| 920 | G_random_element= |
---|
[c4f4fd] | 921 | GCD_small_p (ppA (random_element,x), ppB (random_element,x), top_level, |
---|
[10af64] | 922 | list); |
---|
| 923 | TIMING_END_AND_PRINT (gcd_recursion, |
---|
| 924 | "time for recursive call: "); |
---|
| 925 | DEBOUTLN (cerr, "G_random_element= " << G_random_element); |
---|
| 926 | } |
---|
| 927 | else if (!fail && inextension) |
---|
| 928 | { |
---|
| 929 | CFList list; |
---|
| 930 | TIMING_START (gcd_recursion); |
---|
| 931 | G_random_element= |
---|
| 932 | GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, |
---|
[c4f4fd] | 933 | list, top_level); |
---|
[10af64] | 934 | TIMING_END_AND_PRINT (gcd_recursion, |
---|
| 935 | "time for recursive call: "); |
---|
| 936 | DEBOUTLN (cerr, "G_random_element= " << G_random_element); |
---|
| 937 | } |
---|
| 938 | else if (fail && !inextension) |
---|
| 939 | { |
---|
| 940 | source= CFList(); |
---|
| 941 | dest= CFList(); |
---|
| 942 | CFList list; |
---|
[a5cc7b3] | 943 | CanonicalForm mipo; |
---|
[c4f4fd] | 944 | int deg= 3; |
---|
[a5cc7b3] | 945 | do { |
---|
[04dd0c] | 946 | mipo= randomIrredpoly (deg, x); |
---|
[a5cc7b3] | 947 | alpha= rootOf (mipo); |
---|
| 948 | inextension= true; |
---|
| 949 | fail= false; |
---|
| 950 | random_element= randomElement (m, alpha, l, fail); |
---|
| 951 | deg++; |
---|
[c4f4fd] | 952 | } while (fail); |
---|
[10af64] | 953 | list= CFList(); |
---|
| 954 | TIMING_START (gcd_recursion); |
---|
| 955 | G_random_element= |
---|
| 956 | GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, |
---|
[c4f4fd] | 957 | list, top_level); |
---|
[10af64] | 958 | TIMING_END_AND_PRINT (gcd_recursion, |
---|
| 959 | "time for recursive call: "); |
---|
| 960 | DEBOUTLN (cerr, "G_random_element= " << G_random_element); |
---|
| 961 | } |
---|
| 962 | else if (fail && inextension) |
---|
| 963 | { |
---|
| 964 | source= CFList(); |
---|
| 965 | dest= CFList(); |
---|
| 966 | int num_vars= tmin (getNumVars(A), getNumVars(B)); |
---|
| 967 | num_vars--; |
---|
| 968 | V_buf= alpha; |
---|
| 969 | choose_extension (d, num_vars, V_buf); |
---|
| 970 | bool prim_fail= false; |
---|
| 971 | Variable V_buf2; |
---|
| 972 | CanonicalForm prim_elem, im_prim_elem; |
---|
| 973 | prim_elem= primitiveElement (alpha, V_buf2, prim_fail); |
---|
[c4f4fd] | 974 | |
---|
[10af64] | 975 | ASSERT (!prim_fail, "failure in integer factorizer"); |
---|
| 976 | if (prim_fail) |
---|
| 977 | ; //ERROR |
---|
| 978 | else |
---|
| 979 | im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf); |
---|
| 980 | |
---|
| 981 | DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha)); |
---|
| 982 | DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2)); |
---|
| 983 | |
---|
| 984 | inextensionextension= true; |
---|
| 985 | for (CFListIterator i= l; i.hasItem(); i++) |
---|
| 986 | i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, |
---|
| 987 | im_prim_elem, source, dest); |
---|
| 988 | m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); |
---|
| 989 | G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); |
---|
| 990 | newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, |
---|
| 991 | source, dest); |
---|
| 992 | ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest); |
---|
| 993 | ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest); |
---|
| 994 | gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, |
---|
| 995 | source, dest); |
---|
| 996 | fail= false; |
---|
| 997 | random_element= randomElement (m, V_buf, l, fail ); |
---|
[a5cc7b3] | 998 | DEBOUTLN (cerr, "fail= " << fail); |
---|
[10af64] | 999 | CFList list; |
---|
| 1000 | TIMING_START (gcd_recursion); |
---|
| 1001 | G_random_element= |
---|
| 1002 | GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, |
---|
[c4f4fd] | 1003 | list, top_level); |
---|
[10af64] | 1004 | TIMING_END_AND_PRINT (gcd_recursion, |
---|
| 1005 | "time for recursive call: "); |
---|
| 1006 | DEBOUTLN (cerr, "G_random_element= " << G_random_element); |
---|
| 1007 | alpha= V_buf; |
---|
| 1008 | } |
---|
| 1009 | |
---|
| 1010 | d0= totaldegree (G_random_element, Variable(2), |
---|
| 1011 | Variable (G_random_element.level())); |
---|
| 1012 | |
---|
| 1013 | if (d0 == 0) |
---|
| 1014 | { |
---|
[c4f4fd] | 1015 | return N(gcdcAcB); |
---|
[10af64] | 1016 | } |
---|
| 1017 | if (d0 > d) |
---|
| 1018 | { |
---|
| 1019 | if (!find (l, random_element)) |
---|
| 1020 | l.append (random_element); |
---|
| 1021 | continue; |
---|
| 1022 | } |
---|
| 1023 | |
---|
| 1024 | G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element)) |
---|
| 1025 | *G_random_element; |
---|
| 1026 | |
---|
| 1027 | |
---|
| 1028 | d0= totaldegree (G_random_element, Variable(2), |
---|
| 1029 | Variable(G_random_element.level())); |
---|
| 1030 | |
---|
| 1031 | if (d0 < d) |
---|
| 1032 | { |
---|
| 1033 | m= gcdlcAlcB; |
---|
| 1034 | newtonPoly= 1; |
---|
| 1035 | G_m= 0; |
---|
| 1036 | d= d0; |
---|
| 1037 | } |
---|
[04dd0c] | 1038 | |
---|
[10af64] | 1039 | TIMING_START (newton_interpolation); |
---|
| 1040 | H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x); |
---|
| 1041 | TIMING_END_AND_PRINT (newton_interpolation, |
---|
| 1042 | "time for newton_interpolation: "); |
---|
| 1043 | |
---|
| 1044 | //termination test |
---|
| 1045 | if (uni_lcoeff (H) == gcdlcAlcB) |
---|
| 1046 | { |
---|
| 1047 | cH= uni_content (H); |
---|
| 1048 | ppH= H/cH; |
---|
| 1049 | ppH /= Lc (ppH); |
---|
| 1050 | DEBOUTLN (cerr, "ppH= " << ppH); |
---|
[c4f4fd] | 1051 | if (fdivides (ppH, A) && fdivides (ppH, B)) |
---|
[10af64] | 1052 | return N(gcdcAcB*ppH); |
---|
| 1053 | } |
---|
| 1054 | |
---|
| 1055 | G_m= H; |
---|
| 1056 | newtonPoly= newtonPoly*(x - random_element); |
---|
| 1057 | m= m*(x - random_element); |
---|
| 1058 | if (!find (l, random_element)) |
---|
| 1059 | l.append (random_element); |
---|
| 1060 | } while (1); |
---|
| 1061 | } |
---|
| 1062 | |
---|
| 1063 | #endif |
---|