[ad3c3ff] | 1 | /*****************************************************************************\ |
---|
[806c18] | 2 | * Computer Algebra System SINGULAR |
---|
[ad3c3ff] | 3 | \*****************************************************************************/ |
---|
| 4 | /** @file facHensel.cc |
---|
[806c18] | 5 | * |
---|
[0e2e23] | 6 | * This file implements functions to lift factors via Hensel lifting. |
---|
[806c18] | 7 | * |
---|
| 8 | * ABSTRACT: Hensel lifting is described in "Efficient Multivariate |
---|
| 9 | * Factorization over Finite Fields" by L. Bernardin & M. Monagon. Division with |
---|
| 10 | * remainder is described in "Fast Recursive Division" by C. Burnikel and |
---|
| 11 | * J. Ziegler. Karatsuba multiplication is described in "Modern Computer |
---|
| 12 | * Algebra" by J. von zur Gathen and J. Gerhard. |
---|
[ad3c3ff] | 13 | * |
---|
| 14 | * @author Martin Lee |
---|
| 15 | * |
---|
| 16 | **/ |
---|
| 17 | /*****************************************************************************/ |
---|
| 18 | |
---|
[e4fe2b] | 19 | #include "config.h" |
---|
| 20 | |
---|
[650f2d8] | 21 | #include "cf_assert.h" |
---|
[ad3c3ff] | 22 | #include "debug.h" |
---|
| 23 | #include "timing.h" |
---|
| 24 | |
---|
[81d96c] | 25 | #include "algext.h" |
---|
[ad3c3ff] | 26 | #include "facHensel.h" |
---|
[0e2e23] | 27 | #include "facMul.h" |
---|
[08daea] | 28 | #include "fac_util.h" |
---|
| 29 | #include "cf_algorithm.h" |
---|
[4a05ed] | 30 | #include "cf_primes.h" |
---|
[4a0a303] | 31 | #include "facBivar.h" |
---|
[ad3c3ff] | 32 | |
---|
| 33 | #ifdef HAVE_NTL |
---|
| 34 | #include <NTL/lzz_pEX.h> |
---|
| 35 | #include "NTLconvert.h" |
---|
| 36 | |
---|
[4a05ed] | 37 | static |
---|
| 38 | CFList productsNTL (const CFList& factors, const CanonicalForm& M) |
---|
| 39 | { |
---|
| 40 | zz_p::init (getCharacteristic()); |
---|
| 41 | zz_pX NTLMipo= convertFacCF2NTLzzpX (M); |
---|
| 42 | zz_pE::init (NTLMipo); |
---|
| 43 | zz_pEX prod; |
---|
| 44 | vec_zz_pEX v; |
---|
| 45 | v.SetLength (factors.length()); |
---|
| 46 | int j= 0; |
---|
| 47 | for (CFListIterator i= factors; i.hasItem(); i++, j++) |
---|
| 48 | { |
---|
| 49 | if (i.getItem().inCoeffDomain()) |
---|
| 50 | v[j]= to_zz_pEX (to_zz_pE (convertFacCF2NTLzzpX (i.getItem()))); |
---|
| 51 | else |
---|
| 52 | v[j]= convertFacCF2NTLzz_pEX (i.getItem(), NTLMipo); |
---|
| 53 | } |
---|
| 54 | CFList result; |
---|
| 55 | Variable x= Variable (1); |
---|
| 56 | for (int j= 0; j < factors.length(); j++) |
---|
| 57 | { |
---|
| 58 | int k= 0; |
---|
| 59 | set(prod); |
---|
| 60 | for (int i= 0; i < factors.length(); i++, k++) |
---|
| 61 | { |
---|
| 62 | if (k == j) |
---|
| 63 | continue; |
---|
| 64 | prod *= v[i]; |
---|
| 65 | } |
---|
| 66 | result.append (convertNTLzz_pEX2CF (prod, x, M.mvar())); |
---|
| 67 | } |
---|
| 68 | return result; |
---|
| 69 | } |
---|
| 70 | |
---|
| 71 | static |
---|
| 72 | void tryDiophantine (CFList& result, const CanonicalForm& F, |
---|
| 73 | const CFList& factors, const CanonicalForm& M, bool& fail) |
---|
| 74 | { |
---|
| 75 | ASSERT (M.isUnivariate(), "expected univariate poly"); |
---|
| 76 | |
---|
| 77 | CFList bufFactors= factors; |
---|
| 78 | bufFactors.removeFirst(); |
---|
| 79 | bufFactors.insert (factors.getFirst () (0,2)); |
---|
| 80 | CanonicalForm inv, leadingCoeff= Lc (F); |
---|
| 81 | CFListIterator i= bufFactors; |
---|
| 82 | if (bufFactors.getFirst().inCoeffDomain()) |
---|
| 83 | { |
---|
| 84 | if (i.hasItem()) |
---|
| 85 | i++; |
---|
| 86 | } |
---|
| 87 | for (; i.hasItem(); i++) |
---|
| 88 | { |
---|
| 89 | tryInvert (Lc (i.getItem()), M, inv ,fail); |
---|
| 90 | if (fail) |
---|
| 91 | return; |
---|
| 92 | i.getItem()= reduce (i.getItem()*inv, M); |
---|
| 93 | } |
---|
| 94 | bufFactors= productsNTL (bufFactors, M); |
---|
| 95 | |
---|
| 96 | CanonicalForm buf1, buf2, buf3, S, T; |
---|
| 97 | i= bufFactors; |
---|
| 98 | if (i.hasItem()) |
---|
| 99 | i++; |
---|
| 100 | buf1= bufFactors.getFirst(); |
---|
| 101 | buf2= i.getItem(); |
---|
| 102 | tryExtgcd (buf1, buf2, M, buf3, S, T, fail); |
---|
| 103 | if (fail) |
---|
| 104 | return; |
---|
| 105 | result.append (S); |
---|
| 106 | result.append (T); |
---|
| 107 | if (i.hasItem()) |
---|
| 108 | i++; |
---|
| 109 | for (; i.hasItem(); i++) |
---|
| 110 | { |
---|
| 111 | buf1= i.getItem(); |
---|
| 112 | tryExtgcd (buf3, buf1, M, buf3, S, T, fail); |
---|
| 113 | if (fail) |
---|
| 114 | return; |
---|
| 115 | CFListIterator k= factors; |
---|
| 116 | for (CFListIterator j= result; j.hasItem(); j++, k++) |
---|
| 117 | { |
---|
| 118 | j.getItem() *= S; |
---|
| 119 | j.getItem()= mod (j.getItem(), k.getItem()); |
---|
| 120 | j.getItem()= reduce (j.getItem(), M); |
---|
| 121 | } |
---|
| 122 | result.append (T); |
---|
| 123 | } |
---|
| 124 | } |
---|
| 125 | |
---|
[81d96c] | 126 | static |
---|
[4a05ed] | 127 | CFList mapinto (const CFList& L) |
---|
| 128 | { |
---|
| 129 | CFList result; |
---|
| 130 | for (CFListIterator i= L; i.hasItem(); i++) |
---|
| 131 | result.append (mapinto (i.getItem())); |
---|
| 132 | return result; |
---|
| 133 | } |
---|
| 134 | |
---|
[81d96c] | 135 | static |
---|
[4a05ed] | 136 | int mod (const CFList& L, const CanonicalForm& p) |
---|
| 137 | { |
---|
| 138 | for (CFListIterator i= L; i.hasItem(); i++) |
---|
| 139 | { |
---|
| 140 | if (mod (i.getItem(), p) == 0) |
---|
| 141 | return 0; |
---|
| 142 | } |
---|
| 143 | return 1; |
---|
| 144 | } |
---|
| 145 | |
---|
| 146 | |
---|
[81d96c] | 147 | static void |
---|
[4a05ed] | 148 | chineseRemainder (const CFList & x1, const CanonicalForm & q1, |
---|
| 149 | const CFList & x2, const CanonicalForm & q2, |
---|
| 150 | CFList & xnew, CanonicalForm & qnew) |
---|
| 151 | { |
---|
| 152 | ASSERT (x1.length() == x2.length(), "expected lists of equal length"); |
---|
| 153 | CanonicalForm tmp1, tmp2; |
---|
| 154 | CFListIterator j= x2; |
---|
| 155 | for (CFListIterator i= x1; i.hasItem() && j.hasItem(); i++, j++) |
---|
| 156 | { |
---|
| 157 | chineseRemainder (i.getItem(), q1, j.getItem(), q2, tmp1, tmp2); |
---|
| 158 | xnew.append (tmp1); |
---|
| 159 | } |
---|
| 160 | qnew= tmp2; |
---|
| 161 | } |
---|
| 162 | |
---|
[81d96c] | 163 | static |
---|
[4a05ed] | 164 | CFList Farey (const CFList& L, const CanonicalForm& q) |
---|
| 165 | { |
---|
| 166 | CFList result; |
---|
| 167 | for (CFListIterator i= L; i.hasItem(); i++) |
---|
| 168 | result.append (Farey (i.getItem(), q)); |
---|
| 169 | return result; |
---|
| 170 | } |
---|
| 171 | |
---|
[81d96c] | 172 | static |
---|
[4a05ed] | 173 | CFList replacevar (const CFList& L, const Variable& a, const Variable& b) |
---|
| 174 | { |
---|
| 175 | CFList result; |
---|
| 176 | for (CFListIterator i= L; i.hasItem(); i++) |
---|
| 177 | result.append (replacevar (i.getItem(), a, b)); |
---|
| 178 | return result; |
---|
| 179 | } |
---|
| 180 | |
---|
| 181 | CFList |
---|
| 182 | modularDiophant (const CanonicalForm& f, const CFList& factors, |
---|
| 183 | const CanonicalForm& M) |
---|
| 184 | { |
---|
| 185 | bool save_rat=!isOn (SW_RATIONAL); |
---|
| 186 | On (SW_RATIONAL); |
---|
| 187 | CanonicalForm F= f*bCommonDen (f); |
---|
| 188 | CFList products= factors; |
---|
| 189 | for (CFListIterator i= products; i.hasItem(); i++) |
---|
[9b2c8a] | 190 | { |
---|
| 191 | if (products.getFirst().level() == 1) |
---|
| 192 | i.getItem() /= Lc (i.getItem()); |
---|
[4a05ed] | 193 | i.getItem() *= bCommonDen (i.getItem()); |
---|
[9b2c8a] | 194 | } |
---|
[4a05ed] | 195 | if (products.getFirst().level() == 1) |
---|
| 196 | products.insert (Lc (F)); |
---|
| 197 | CanonicalForm bound= maxNorm (F); |
---|
| 198 | CFList leadingCoeffs; |
---|
| 199 | leadingCoeffs.append (lc (F)); |
---|
| 200 | CanonicalForm dummy; |
---|
| 201 | for (CFListIterator i= products; i.hasItem(); i++) |
---|
| 202 | { |
---|
| 203 | leadingCoeffs.append (lc (i.getItem())); |
---|
| 204 | dummy= maxNorm (i.getItem()); |
---|
| 205 | bound= (dummy > bound) ? dummy : bound; |
---|
| 206 | } |
---|
| 207 | bound *= maxNorm (Lc (F))*maxNorm (Lc(F))*bound; |
---|
| 208 | bound *= bound*bound; |
---|
| 209 | bound= power (bound, degree (M)); |
---|
| 210 | bound *= power (CanonicalForm (2),degree (f)); |
---|
| 211 | CanonicalForm bufBound= bound; |
---|
| 212 | int i = cf_getNumBigPrimes() - 1; |
---|
| 213 | int p; |
---|
| 214 | CFList resultModP, result, newResult; |
---|
| 215 | CanonicalForm q (0), newQ; |
---|
| 216 | bool fail= false; |
---|
| 217 | Variable a= M.mvar(); |
---|
| 218 | Variable b= Variable (2); |
---|
| 219 | setReduce (M.mvar(), false); |
---|
| 220 | CanonicalForm mipo= bCommonDen (M)*M; |
---|
| 221 | Off (SW_RATIONAL); |
---|
| 222 | CanonicalForm modMipo; |
---|
| 223 | leadingCoeffs.append (lc (mipo)); |
---|
| 224 | CFList tmp1, tmp2; |
---|
| 225 | bool equal= false; |
---|
| 226 | int count= 0; |
---|
| 227 | do |
---|
| 228 | { |
---|
| 229 | p = cf_getBigPrime( i ); |
---|
| 230 | i--; |
---|
| 231 | while ( i >= 0 && mod( leadingCoeffs, p ) == 0) |
---|
| 232 | { |
---|
| 233 | p = cf_getBigPrime( i ); |
---|
| 234 | i--; |
---|
| 235 | } |
---|
| 236 | |
---|
| 237 | ASSERT (i >= 0, "ran out of primes"); //sic |
---|
| 238 | |
---|
| 239 | setCharacteristic (p); |
---|
| 240 | modMipo= mapinto (mipo); |
---|
| 241 | modMipo /= lc (modMipo); |
---|
| 242 | resultModP= CFList(); |
---|
| 243 | tryDiophantine (resultModP, mapinto (F), mapinto (products), modMipo, fail); |
---|
| 244 | setCharacteristic (0); |
---|
| 245 | if (fail) |
---|
| 246 | { |
---|
| 247 | fail= false; |
---|
| 248 | continue; |
---|
| 249 | } |
---|
| 250 | |
---|
| 251 | if ( q.isZero() ) |
---|
| 252 | { |
---|
| 253 | result= replacevar (mapinto(resultModP), a, b); |
---|
| 254 | q= p; |
---|
| 255 | } |
---|
| 256 | else |
---|
| 257 | { |
---|
| 258 | result= replacevar (result, a, b); |
---|
| 259 | newResult= CFList(); |
---|
| 260 | chineseRemainder( result, q, replacevar (mapinto (resultModP), a, b), |
---|
| 261 | p, newResult, newQ ); |
---|
| 262 | q= newQ; |
---|
| 263 | result= newResult; |
---|
| 264 | if (newQ > bound) |
---|
| 265 | { |
---|
| 266 | count++; |
---|
| 267 | tmp1= replacevar (Farey (result, q), b, a); |
---|
| 268 | if (tmp2.isEmpty()) |
---|
| 269 | tmp2= tmp1; |
---|
| 270 | else |
---|
| 271 | { |
---|
| 272 | equal= true; |
---|
| 273 | CFListIterator k= tmp1; |
---|
| 274 | for (CFListIterator j= tmp2; j.hasItem(); j++, k++) |
---|
| 275 | { |
---|
| 276 | if (j.getItem() != k.getItem()) |
---|
| 277 | equal= false; |
---|
| 278 | } |
---|
| 279 | if (!equal) |
---|
| 280 | tmp2= tmp1; |
---|
| 281 | } |
---|
| 282 | if (count > 2) |
---|
| 283 | { |
---|
| 284 | bound *= bufBound; |
---|
| 285 | equal= false; |
---|
| 286 | count= 0; |
---|
| 287 | } |
---|
| 288 | } |
---|
| 289 | if (newQ > bound && equal) |
---|
| 290 | { |
---|
| 291 | On( SW_RATIONAL ); |
---|
| 292 | CFList bufResult= result; |
---|
[9b2c8a] | 293 | result= tmp2; |
---|
[4a05ed] | 294 | setReduce (M.mvar(), true); |
---|
| 295 | if (factors.getFirst().level() == 1) |
---|
| 296 | { |
---|
| 297 | result.removeFirst(); |
---|
| 298 | CFListIterator j= factors; |
---|
| 299 | CanonicalForm denf= bCommonDen (f); |
---|
| 300 | for (CFListIterator i= result; i.hasItem(); i++, j++) |
---|
| 301 | i.getItem() *= Lc (j.getItem())*denf; |
---|
| 302 | } |
---|
| 303 | if (factors.getFirst().level() != 1 && |
---|
| 304 | !bCommonDen (factors.getFirst()).isOne()) |
---|
| 305 | { |
---|
| 306 | CanonicalForm denFirst= bCommonDen (factors.getFirst()); |
---|
| 307 | for (CFListIterator i= result; i.hasItem(); i++) |
---|
| 308 | i.getItem() *= denFirst; |
---|
| 309 | } |
---|
| 310 | |
---|
| 311 | CanonicalForm test= 0; |
---|
| 312 | CFListIterator jj= factors; |
---|
| 313 | for (CFListIterator ii= result; ii.hasItem(); ii++, jj++) |
---|
| 314 | test += ii.getItem()*(f/jj.getItem()); |
---|
| 315 | if (!test.isOne()) |
---|
| 316 | { |
---|
| 317 | bound *= bufBound; |
---|
| 318 | equal= false; |
---|
| 319 | count= 0; |
---|
| 320 | setReduce (M.mvar(), false); |
---|
| 321 | result= bufResult; |
---|
| 322 | Off (SW_RATIONAL); |
---|
| 323 | } |
---|
| 324 | else |
---|
| 325 | break; |
---|
| 326 | } |
---|
| 327 | } |
---|
| 328 | } while (1); |
---|
| 329 | if (save_rat) Off(SW_RATIONAL); |
---|
| 330 | return result; |
---|
| 331 | } |
---|
| 332 | |
---|
[806c18] | 333 | void sortList (CFList& list, const Variable& x) |
---|
| 334 | { |
---|
[ad3c3ff] | 335 | int l= 1; |
---|
| 336 | int k= 1; |
---|
| 337 | CanonicalForm buf; |
---|
| 338 | CFListIterator m; |
---|
[806c18] | 339 | for (CFListIterator i= list; l <= list.length(); i++, l++) |
---|
[ad3c3ff] | 340 | { |
---|
[806c18] | 341 | for (CFListIterator j= list; k <= list.length() - l; k++) |
---|
[ad3c3ff] | 342 | { |
---|
| 343 | m= j; |
---|
| 344 | m++; |
---|
[806c18] | 345 | if (degree (j.getItem(), x) > degree (m.getItem(), x)) |
---|
[ad3c3ff] | 346 | { |
---|
| 347 | buf= m.getItem(); |
---|
| 348 | m.getItem()= j.getItem(); |
---|
| 349 | j.getItem()= buf; |
---|
| 350 | j++; |
---|
| 351 | j.getItem()= m.getItem(); |
---|
[806c18] | 352 | } |
---|
[ad3c3ff] | 353 | else |
---|
| 354 | j++; |
---|
| 355 | } |
---|
| 356 | k= 1; |
---|
| 357 | } |
---|
| 358 | } |
---|
| 359 | |
---|
[583cb9] | 360 | CFList |
---|
| 361 | diophantine (const CanonicalForm& F, const CFList& factors); |
---|
| 362 | |
---|
| 363 | CFList |
---|
| 364 | diophantineHensel (const CanonicalForm & F, const CFList& factors, |
---|
| 365 | const modpk& b) |
---|
| 366 | { |
---|
| 367 | int p= b.getp(); |
---|
| 368 | setCharacteristic (p); |
---|
| 369 | CFList recResult= diophantine (mapinto (F), mapinto (factors)); |
---|
| 370 | setCharacteristic (0); |
---|
| 371 | recResult= mapinto (recResult); |
---|
| 372 | CanonicalForm e= 1; |
---|
| 373 | CFList L; |
---|
| 374 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 375 | int k= 0; |
---|
| 376 | for (CFListIterator i= factors; i.hasItem(); i++, k++) |
---|
| 377 | { |
---|
| 378 | if (k == 0) |
---|
| 379 | bufFactors[k]= i.getItem() (0); |
---|
| 380 | else |
---|
| 381 | bufFactors [k]= i.getItem(); |
---|
| 382 | } |
---|
| 383 | CanonicalForm tmp, quot; |
---|
| 384 | for (k= 0; k < factors.length(); k++) //TODO compute b's faster |
---|
| 385 | { |
---|
| 386 | tmp= 1; |
---|
| 387 | for (int l= 0; l < factors.length(); l++) |
---|
| 388 | { |
---|
| 389 | if (l == k) |
---|
| 390 | continue; |
---|
| 391 | else |
---|
| 392 | { |
---|
| 393 | tmp= mulNTL (tmp, bufFactors[l]); |
---|
| 394 | } |
---|
| 395 | } |
---|
| 396 | L.append (tmp); |
---|
| 397 | } |
---|
| 398 | |
---|
| 399 | setCharacteristic (p); |
---|
| 400 | for (k= 0; k < factors.length(); k++) |
---|
| 401 | bufFactors [k]= bufFactors[k].mapinto(); |
---|
| 402 | setCharacteristic(0); |
---|
| 403 | |
---|
| 404 | CFListIterator j= L; |
---|
| 405 | for (CFListIterator i= recResult; i.hasItem(); i++, j++) |
---|
| 406 | e= b (e - mulNTL (i.getItem(),j.getItem(), b)); |
---|
| 407 | |
---|
| 408 | if (e.isZero()) |
---|
| 409 | return recResult; |
---|
| 410 | CanonicalForm coeffE; |
---|
| 411 | CFList s; |
---|
| 412 | CFList result= recResult; |
---|
[4e6d2a] | 413 | setCharacteristic (p); |
---|
| 414 | recResult= mapinto (recResult); |
---|
| 415 | setCharacteristic (0); |
---|
[583cb9] | 416 | CanonicalForm g; |
---|
| 417 | CanonicalForm modulus= p; |
---|
| 418 | int d= b.getk(); |
---|
| 419 | for (int i= 1; i < d; i++) |
---|
| 420 | { |
---|
| 421 | coeffE= div (e, modulus); |
---|
| 422 | setCharacteristic (p); |
---|
| 423 | coeffE= coeffE.mapinto(); |
---|
| 424 | setCharacteristic (0); |
---|
| 425 | if (!coeffE.isZero()) |
---|
| 426 | { |
---|
| 427 | CFListIterator k= result; |
---|
| 428 | CFListIterator l= L; |
---|
| 429 | int ii= 0; |
---|
| 430 | j= recResult; |
---|
| 431 | for (; j.hasItem(); j++, k++, l++, ii++) |
---|
| 432 | { |
---|
| 433 | setCharacteristic (p); |
---|
[4e6d2a] | 434 | g= mulNTL (coeffE, j.getItem()); |
---|
[c3ef86f] | 435 | g= modNTL (g, bufFactors[ii]); |
---|
[583cb9] | 436 | setCharacteristic (0); |
---|
| 437 | k.getItem() += g.mapinto()*modulus; |
---|
| 438 | e -= mulNTL (g.mapinto()*modulus, l.getItem(), b); |
---|
| 439 | e= b(e); |
---|
| 440 | } |
---|
| 441 | } |
---|
| 442 | modulus *= p; |
---|
| 443 | if (e.isZero()) |
---|
| 444 | break; |
---|
| 445 | } |
---|
| 446 | |
---|
| 447 | return result; |
---|
| 448 | } |
---|
| 449 | |
---|
[4a0a303] | 450 | CFList |
---|
| 451 | diophantineHenselQa (const CanonicalForm & F, const CanonicalForm& G, |
---|
| 452 | const CFList& factors, modpk& b, const Variable& alpha) |
---|
| 453 | { |
---|
| 454 | bool fail= false; |
---|
| 455 | CFList recResult; |
---|
| 456 | CanonicalForm modMipo, mipo; |
---|
[f9da5e] | 457 | //here SW_RATIONAL is off |
---|
| 458 | On (SW_RATIONAL); |
---|
[4a0a303] | 459 | mipo= getMipo (alpha); |
---|
[f9da5e] | 460 | bool mipoHasDen= false; |
---|
| 461 | if (!bCommonDen (mipo).isOne()) |
---|
| 462 | { |
---|
| 463 | mipo *= bCommonDen (mipo); |
---|
| 464 | mipoHasDen= true; |
---|
| 465 | } |
---|
| 466 | Off (SW_RATIONAL); |
---|
| 467 | int p= b.getp(); |
---|
| 468 | setCharacteristic (p); |
---|
[4a0a303] | 469 | setReduce (alpha, false); |
---|
| 470 | while (1) |
---|
| 471 | { |
---|
| 472 | setCharacteristic (p); |
---|
| 473 | modMipo= mapinto (mipo); |
---|
| 474 | modMipo /= lc (modMipo); |
---|
| 475 | tryDiophantine (recResult, mapinto (F), mapinto (factors), modMipo, fail); |
---|
| 476 | if (fail) |
---|
| 477 | { |
---|
| 478 | int i= 0; |
---|
| 479 | while (cf_getBigPrime (i) < p) |
---|
| 480 | i++; |
---|
| 481 | findGoodPrime (F, i); |
---|
| 482 | findGoodPrime (G, i); |
---|
| 483 | p=cf_getBigPrime(i); |
---|
| 484 | b = coeffBound( G, p, mipo ); |
---|
| 485 | modpk bb= coeffBound (F, p, mipo ); |
---|
| 486 | if (bb.getk() > b.getk() ) b=bb; |
---|
[1ddcde9] | 487 | fail= false; |
---|
[4a0a303] | 488 | } |
---|
| 489 | else |
---|
| 490 | break; |
---|
| 491 | } |
---|
| 492 | setCharacteristic (0); |
---|
| 493 | recResult= mapinto (recResult); |
---|
| 494 | setReduce (alpha, true); |
---|
| 495 | CanonicalForm e= 1; |
---|
| 496 | CFList L; |
---|
| 497 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 498 | int k= 0; |
---|
| 499 | for (CFListIterator i= factors; i.hasItem(); i++, k++) |
---|
| 500 | { |
---|
| 501 | if (k == 0) |
---|
| 502 | bufFactors[k]= i.getItem() (0); |
---|
| 503 | else |
---|
| 504 | bufFactors [k]= i.getItem(); |
---|
| 505 | } |
---|
[f9da5e] | 506 | CanonicalForm tmp; |
---|
| 507 | On (SW_RATIONAL); |
---|
[4a0a303] | 508 | for (k= 0; k < factors.length(); k++) //TODO compute b's faster |
---|
| 509 | { |
---|
| 510 | tmp= 1; |
---|
| 511 | for (int l= 0; l < factors.length(); l++) |
---|
| 512 | { |
---|
| 513 | if (l == k) |
---|
| 514 | continue; |
---|
| 515 | else |
---|
| 516 | tmp= mulNTL (tmp, bufFactors[l]); |
---|
| 517 | } |
---|
[f9da5e] | 518 | L.append (tmp*bCommonDen(tmp)); |
---|
| 519 | } |
---|
| 520 | |
---|
| 521 | Variable gamma; |
---|
| 522 | CanonicalForm den; |
---|
| 523 | if (mipoHasDen) |
---|
| 524 | { |
---|
| 525 | modMipo= getMipo (alpha); |
---|
| 526 | den= bCommonDen (modMipo); |
---|
| 527 | modMipo *= den; |
---|
| 528 | Off (SW_RATIONAL); |
---|
| 529 | setReduce (alpha, false); |
---|
| 530 | gamma= rootOf (b (modMipo*b.inverse (den))); |
---|
| 531 | setReduce (alpha, true); |
---|
[4a0a303] | 532 | } |
---|
| 533 | |
---|
| 534 | setCharacteristic (p); |
---|
[f9da5e] | 535 | Variable beta; |
---|
| 536 | Off (SW_RATIONAL); |
---|
| 537 | if (mipoHasDen) |
---|
| 538 | { |
---|
| 539 | setReduce (alpha, false); |
---|
| 540 | modMipo= modMipo.mapinto(); |
---|
| 541 | modMipo /= lc (modMipo); |
---|
| 542 | beta= rootOf (modMipo); |
---|
| 543 | setReduce (alpha, true); |
---|
| 544 | } |
---|
| 545 | |
---|
[4a0a303] | 546 | for (k= 0; k < factors.length(); k++) |
---|
[f9da5e] | 547 | { |
---|
| 548 | if (!mipoHasDen) |
---|
| 549 | bufFactors [k]= bufFactors[k].mapinto(); |
---|
| 550 | else |
---|
| 551 | { |
---|
| 552 | bufFactors [k]= bufFactors[k].mapinto(); |
---|
| 553 | bufFactors [k]= replacevar (bufFactors[k], alpha, beta); |
---|
| 554 | } |
---|
| 555 | } |
---|
[4a0a303] | 556 | setCharacteristic(0); |
---|
| 557 | |
---|
| 558 | CFListIterator j= L; |
---|
[f9da5e] | 559 | for (;j.hasItem(); j++) |
---|
| 560 | { |
---|
| 561 | if (mipoHasDen) |
---|
| 562 | j.getItem()= replacevar (b(j.getItem()*b.inverse(lc(j.getItem()))), |
---|
| 563 | alpha, gamma); |
---|
| 564 | else |
---|
| 565 | j.getItem()= b(j.getItem()*b.inverse(lc(j.getItem()))); |
---|
| 566 | } |
---|
| 567 | j= L; |
---|
[4a0a303] | 568 | for (CFListIterator i= recResult; i.hasItem(); i++, j++) |
---|
[f9da5e] | 569 | { |
---|
| 570 | if (mipoHasDen) |
---|
| 571 | e= b (e - mulNTL (replacevar (i.getItem(), alpha, gamma),j.getItem(), b)); |
---|
| 572 | else |
---|
| 573 | e= b (e - mulNTL (i.getItem(), j.getItem(), b)); |
---|
| 574 | } |
---|
[4a0a303] | 575 | |
---|
| 576 | if (e.isZero()) |
---|
[f9da5e] | 577 | { |
---|
| 578 | if (mipoHasDen) |
---|
| 579 | { |
---|
| 580 | for (CFListIterator i= recResult; i.hasItem(); i++) |
---|
| 581 | i.getItem()= replacevar (i.getItem(), alpha, gamma); |
---|
| 582 | } |
---|
[4a0a303] | 583 | return recResult; |
---|
[f9da5e] | 584 | } |
---|
[4a0a303] | 585 | CanonicalForm coeffE; |
---|
| 586 | CFList result= recResult; |
---|
[f9da5e] | 587 | if (mipoHasDen) |
---|
| 588 | { |
---|
| 589 | for (CFListIterator i= result; i.hasItem(); i++) |
---|
| 590 | i.getItem()= replacevar (i.getItem(), alpha, gamma); |
---|
| 591 | } |
---|
[4a0a303] | 592 | setCharacteristic (p); |
---|
| 593 | recResult= mapinto (recResult); |
---|
[f9da5e] | 594 | if (mipoHasDen) |
---|
| 595 | { |
---|
| 596 | for (CFListIterator i= recResult; i.hasItem(); i++) |
---|
| 597 | i.getItem()= replacevar (i.getItem(), alpha, beta); |
---|
| 598 | } |
---|
| 599 | |
---|
[4a0a303] | 600 | setCharacteristic (0); |
---|
| 601 | CanonicalForm g; |
---|
| 602 | CanonicalForm modulus= p; |
---|
| 603 | int d= b.getk(); |
---|
| 604 | for (int i= 1; i < d; i++) |
---|
| 605 | { |
---|
| 606 | coeffE= div (e, modulus); |
---|
| 607 | setCharacteristic (p); |
---|
| 608 | coeffE= coeffE.mapinto(); |
---|
[f9da5e] | 609 | if (mipoHasDen) |
---|
| 610 | coeffE= replacevar (coeffE, gamma, beta); |
---|
[4a0a303] | 611 | setCharacteristic (0); |
---|
| 612 | if (!coeffE.isZero()) |
---|
| 613 | { |
---|
| 614 | CFListIterator k= result; |
---|
| 615 | CFListIterator l= L; |
---|
| 616 | int ii= 0; |
---|
| 617 | j= recResult; |
---|
| 618 | for (; j.hasItem(); j++, k++, l++, ii++) |
---|
| 619 | { |
---|
| 620 | setCharacteristic (p); |
---|
| 621 | g= mulNTL (coeffE, j.getItem()); |
---|
| 622 | g= modNTL (g, bufFactors[ii]); |
---|
| 623 | setCharacteristic (0); |
---|
[f9da5e] | 624 | if (mipoHasDen) |
---|
| 625 | { |
---|
| 626 | k.getItem() += replacevar (g.mapinto()*modulus, beta, gamma); |
---|
| 627 | e -= mulNTL (replacevar (g.mapinto(), beta, gamma)*modulus, |
---|
| 628 | l.getItem(), b); |
---|
| 629 | } |
---|
| 630 | else |
---|
| 631 | { |
---|
| 632 | k.getItem() += g.mapinto()*modulus; |
---|
| 633 | e -= mulNTL (g.mapinto()*modulus, l.getItem(), b); |
---|
| 634 | } |
---|
[4a0a303] | 635 | e= b(e); |
---|
| 636 | } |
---|
| 637 | } |
---|
| 638 | modulus *= p; |
---|
| 639 | if (e.isZero()) |
---|
| 640 | break; |
---|
| 641 | } |
---|
| 642 | |
---|
| 643 | return result; |
---|
| 644 | } |
---|
| 645 | |
---|
[583cb9] | 646 | CFList |
---|
[69fdf90] | 647 | diophantine (const CanonicalForm& F, const CanonicalForm& G, |
---|
| 648 | const CFList& factors, modpk& b) |
---|
[ad3c3ff] | 649 | { |
---|
[4a05ed] | 650 | if (getCharacteristic() == 0) |
---|
| 651 | { |
---|
| 652 | Variable v; |
---|
| 653 | bool hasAlgVar= hasFirstAlgVar (F, v); |
---|
| 654 | for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++) |
---|
| 655 | hasAlgVar= hasFirstAlgVar (i.getItem(), v); |
---|
| 656 | if (hasAlgVar) |
---|
| 657 | { |
---|
[4a0a303] | 658 | if (b.getp() != 0) |
---|
| 659 | { |
---|
| 660 | CFList result= diophantineHenselQa (F, G, factors, b, v); |
---|
| 661 | return result; |
---|
| 662 | } |
---|
[4a05ed] | 663 | CFList result= modularDiophant (F, factors, getMipo (v)); |
---|
| 664 | return result; |
---|
| 665 | } |
---|
[4a0a303] | 666 | if (b.getp() != 0) |
---|
| 667 | return diophantineHensel (F, factors, b); |
---|
[4a05ed] | 668 | } |
---|
| 669 | |
---|
[ad3c3ff] | 670 | CanonicalForm buf1, buf2, buf3, S, T; |
---|
| 671 | CFListIterator i= factors; |
---|
| 672 | CFList result; |
---|
| 673 | if (i.hasItem()) |
---|
| 674 | i++; |
---|
[806c18] | 675 | buf1= F/factors.getFirst(); |
---|
[ad3c3ff] | 676 | buf2= divNTL (F, i.getItem()); |
---|
| 677 | buf3= extgcd (buf1, buf2, S, T); |
---|
| 678 | result.append (S); |
---|
| 679 | result.append (T); |
---|
| 680 | if (i.hasItem()) |
---|
| 681 | i++; |
---|
| 682 | for (; i.hasItem(); i++) |
---|
| 683 | { |
---|
| 684 | buf1= divNTL (F, i.getItem()); |
---|
| 685 | buf3= extgcd (buf3, buf1, S, T); |
---|
| 686 | CFListIterator k= factors; |
---|
| 687 | for (CFListIterator j= result; j.hasItem(); j++, k++) |
---|
| 688 | { |
---|
| 689 | j.getItem()= mulNTL (j.getItem(), S); |
---|
| 690 | j.getItem()= modNTL (j.getItem(), k.getItem()); |
---|
| 691 | } |
---|
| 692 | result.append (T); |
---|
| 693 | } |
---|
| 694 | return result; |
---|
| 695 | } |
---|
| 696 | |
---|
[583cb9] | 697 | CFList |
---|
| 698 | diophantine (const CanonicalForm& F, const CFList& factors) |
---|
| 699 | { |
---|
| 700 | modpk b= modpk(); |
---|
[69fdf90] | 701 | return diophantine (F, 1, factors, b); |
---|
[583cb9] | 702 | } |
---|
| 703 | |
---|
[806c18] | 704 | void |
---|
| 705 | henselStep12 (const CanonicalForm& F, const CFList& factors, |
---|
[ad3c3ff] | 706 | CFArray& bufFactors, const CFList& diophant, CFMatrix& M, |
---|
[a090c88] | 707 | CFArray& Pi, int j, const modpk& b) |
---|
[ad3c3ff] | 708 | { |
---|
| 709 | CanonicalForm E; |
---|
| 710 | CanonicalForm xToJ= power (F.mvar(), j); |
---|
| 711 | Variable x= F.mvar(); |
---|
| 712 | // compute the error |
---|
| 713 | if (j == 1) |
---|
[806c18] | 714 | E= F[j]; |
---|
[ad3c3ff] | 715 | else |
---|
| 716 | { |
---|
[806c18] | 717 | if (degree (Pi [factors.length() - 2], x) > 0) |
---|
[ad3c3ff] | 718 | E= F[j] - Pi [factors.length() - 2] [j]; |
---|
| 719 | else |
---|
[806c18] | 720 | E= F[j]; |
---|
[ad3c3ff] | 721 | } |
---|
| 722 | |
---|
[a090c88] | 723 | if (b.getp() != 0) |
---|
| 724 | E= b(E); |
---|
[ad3c3ff] | 725 | CFArray buf= CFArray (diophant.length()); |
---|
| 726 | bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1)); |
---|
| 727 | int k= 0; |
---|
| 728 | CanonicalForm remainder; |
---|
| 729 | // actual lifting |
---|
[806c18] | 730 | for (CFListIterator i= diophant; i.hasItem(); i++, k++) |
---|
[ad3c3ff] | 731 | { |
---|
| 732 | if (degree (bufFactors[k], x) > 0) |
---|
| 733 | { |
---|
| 734 | if (k > 0) |
---|
[a090c88] | 735 | remainder= modNTL (E, bufFactors[k] [0], b); |
---|
[ad3c3ff] | 736 | else |
---|
| 737 | remainder= E; |
---|
| 738 | } |
---|
| 739 | else |
---|
[a090c88] | 740 | remainder= modNTL (E, bufFactors[k], b); |
---|
[ad3c3ff] | 741 | |
---|
[a090c88] | 742 | buf[k]= mulNTL (i.getItem(), remainder, b); |
---|
[ad3c3ff] | 743 | if (degree (bufFactors[k], x) > 0) |
---|
[a090c88] | 744 | buf[k]= modNTL (buf[k], bufFactors[k] [0], b); |
---|
[806c18] | 745 | else |
---|
[a090c88] | 746 | buf[k]= modNTL (buf[k], bufFactors[k], b); |
---|
[ad3c3ff] | 747 | } |
---|
| 748 | for (k= 1; k < factors.length(); k++) |
---|
[a090c88] | 749 | { |
---|
[ad3c3ff] | 750 | bufFactors[k] += xToJ*buf[k]; |
---|
[a090c88] | 751 | if (b.getp() != 0) |
---|
| 752 | bufFactors[k]= b(bufFactors[k]); |
---|
| 753 | } |
---|
[ad3c3ff] | 754 | |
---|
| 755 | // update Pi [0] |
---|
| 756 | int degBuf0= degree (bufFactors[0], x); |
---|
| 757 | int degBuf1= degree (bufFactors[1], x); |
---|
| 758 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
[a090c88] | 759 | M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j], b); |
---|
[ad3c3ff] | 760 | CanonicalForm uIZeroJ; |
---|
| 761 | if (j == 1) |
---|
| 762 | { |
---|
| 763 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 764 | uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]), |
---|
[a090c88] | 765 | (bufFactors[1] [0] + buf[1]), b) - M(1, 1) - M(j + 1, 1); |
---|
[ad3c3ff] | 766 | else if (degBuf0 > 0) |
---|
[a090c88] | 767 | uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1], b); |
---|
[ad3c3ff] | 768 | else if (degBuf1 > 0) |
---|
[a090c88] | 769 | uIZeroJ= mulNTL (bufFactors[0], buf[1], b); |
---|
[ad3c3ff] | 770 | else |
---|
| 771 | uIZeroJ= 0; |
---|
[a090c88] | 772 | if (b.getp() != 0) |
---|
| 773 | uIZeroJ= b (uIZeroJ); |
---|
[806c18] | 774 | Pi [0] += xToJ*uIZeroJ; |
---|
[a090c88] | 775 | if (b.getp() != 0) |
---|
| 776 | Pi [0]= b (Pi[0]); |
---|
[806c18] | 777 | } |
---|
[ad3c3ff] | 778 | else |
---|
| 779 | { |
---|
| 780 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 781 | uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]), |
---|
[a090c88] | 782 | (bufFactors[1] [0] + buf[1]), b) - M(1, 1) - M(j + 1, 1); |
---|
[ad3c3ff] | 783 | else if (degBuf0 > 0) |
---|
[a090c88] | 784 | uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1], b); |
---|
[ad3c3ff] | 785 | else if (degBuf1 > 0) |
---|
[a090c88] | 786 | uIZeroJ= mulNTL (bufFactors[0], buf[1], b); |
---|
[ad3c3ff] | 787 | else |
---|
| 788 | uIZeroJ= 0; |
---|
[a090c88] | 789 | if (b.getp() != 0) |
---|
| 790 | uIZeroJ= b (uIZeroJ); |
---|
[806c18] | 791 | Pi [0] += xToJ*uIZeroJ; |
---|
[a090c88] | 792 | if (b.getp() != 0) |
---|
| 793 | Pi [0]= b (Pi[0]); |
---|
[ad3c3ff] | 794 | } |
---|
| 795 | CFArray tmp= CFArray (factors.length() - 1); |
---|
| 796 | for (k= 0; k < factors.length() - 1; k++) |
---|
| 797 | tmp[k]= 0; |
---|
| 798 | CFIterator one, two; |
---|
| 799 | one= bufFactors [0]; |
---|
| 800 | two= bufFactors [1]; |
---|
| 801 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 802 | { |
---|
| 803 | for (k= 1; k <= (int) ceil (j/2.0); k++) |
---|
[806c18] | 804 | { |
---|
[ad3c3ff] | 805 | if (k != j - k + 1) |
---|
| 806 | { |
---|
[e016ba] | 807 | if ((one.hasTerms() && one.exp() == j - k + 1) |
---|
| 808 | && (two.hasTerms() && two.exp() == j - k + 1)) |
---|
[806c18] | 809 | { |
---|
[e016ba] | 810 | tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), (bufFactors[1][k]+ |
---|
[a090c88] | 811 | two.coeff()), b) - M (k + 1, 1) - M (j - k + 2, 1); |
---|
[806c18] | 812 | one++; |
---|
| 813 | two++; |
---|
| 814 | } |
---|
[e368746] | 815 | else if (one.hasTerms() && one.exp() == j - k + 1) |
---|
[806c18] | 816 | { |
---|
[a090c88] | 817 | tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1][k], b) |
---|
| 818 | - M (k + 1, 1); |
---|
[806c18] | 819 | one++; |
---|
| 820 | } |
---|
[e368746] | 821 | else if (two.hasTerms() && two.exp() == j - k + 1) |
---|
[806c18] | 822 | { |
---|
[a090c88] | 823 | tmp[0] += mulNTL (bufFactors[0][k], (bufFactors[1][k]+two.coeff()), b) |
---|
| 824 | - M (k + 1, 1); |
---|
[806c18] | 825 | two++; |
---|
| 826 | } |
---|
[ad3c3ff] | 827 | } |
---|
| 828 | else |
---|
| 829 | { |
---|
[806c18] | 830 | tmp[0] += M (k + 1, 1); |
---|
[ad3c3ff] | 831 | } |
---|
| 832 | } |
---|
| 833 | } |
---|
[a090c88] | 834 | if (b.getp() != 0) |
---|
| 835 | tmp[0]= b (tmp[0]); |
---|
[ad3c3ff] | 836 | Pi [0] += tmp[0]*xToJ*F.mvar(); |
---|
| 837 | |
---|
[806c18] | 838 | // update Pi [l] |
---|
[ad3c3ff] | 839 | int degPi, degBuf; |
---|
[806c18] | 840 | for (int l= 1; l < factors.length() - 1; l++) |
---|
[ad3c3ff] | 841 | { |
---|
| 842 | degPi= degree (Pi [l - 1], x); |
---|
| 843 | degBuf= degree (bufFactors[l + 1], x); |
---|
| 844 | if (degPi > 0 && degBuf > 0) |
---|
[a090c88] | 845 | M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j], b); |
---|
[ad3c3ff] | 846 | if (j == 1) |
---|
| 847 | { |
---|
| 848 | if (degPi > 0 && degBuf > 0) |
---|
| 849 | Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j], |
---|
[a090c88] | 850 | bufFactors[l + 1] [0] + buf[l + 1], b) - M (j + 1, l +1) - |
---|
[ad3c3ff] | 851 | M (1, l + 1)); |
---|
| 852 | else if (degPi > 0) |
---|
[a090c88] | 853 | Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1], b)); |
---|
[ad3c3ff] | 854 | else if (degBuf > 0) |
---|
[a090c88] | 855 | Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1], b)); |
---|
[ad3c3ff] | 856 | } |
---|
| 857 | else |
---|
| 858 | { |
---|
| 859 | if (degPi > 0 && degBuf > 0) |
---|
| 860 | { |
---|
[a090c88] | 861 | uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0], b); |
---|
| 862 | uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1], b); |
---|
[ad3c3ff] | 863 | } |
---|
| 864 | else if (degPi > 0) |
---|
[a090c88] | 865 | uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1], b); |
---|
[ad3c3ff] | 866 | else if (degBuf > 0) |
---|
| 867 | { |
---|
[a090c88] | 868 | uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0], b); |
---|
| 869 | uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1], b); |
---|
[ad3c3ff] | 870 | } |
---|
| 871 | Pi[l] += xToJ*uIZeroJ; |
---|
| 872 | } |
---|
| 873 | one= bufFactors [l + 1]; |
---|
| 874 | two= Pi [l - 1]; |
---|
[e368746] | 875 | if (two.hasTerms() && two.exp() == j + 1) |
---|
[ad3c3ff] | 876 | { |
---|
| 877 | if (degBuf > 0 && degPi > 0) |
---|
| 878 | { |
---|
[a090c88] | 879 | tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0], b); |
---|
[806c18] | 880 | two++; |
---|
[ad3c3ff] | 881 | } |
---|
| 882 | else if (degPi > 0) |
---|
[806c18] | 883 | { |
---|
[a090c88] | 884 | tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1], b); |
---|
[806c18] | 885 | two++; |
---|
| 886 | } |
---|
[ad3c3ff] | 887 | } |
---|
[806c18] | 888 | if (degBuf > 0 && degPi > 0) |
---|
[ad3c3ff] | 889 | { |
---|
| 890 | for (k= 1; k <= (int) ceil (j/2.0); k++) |
---|
[806c18] | 891 | { |
---|
| 892 | if (k != j - k + 1) |
---|
| 893 | { |
---|
[c1b9927] | 894 | if ((one.hasTerms() && one.exp() == j - k + 1) && |
---|
[e368746] | 895 | (two.hasTerms() && two.exp() == j - k + 1)) |
---|
[806c18] | 896 | { |
---|
[e016ba] | 897 | tmp[l] += mulNTL ((bufFactors[l+1][k] + one.coeff()), (Pi[l-1][k] + |
---|
[a090c88] | 898 | two.coeff()),b) - M (k + 1, l + 1) - M (j - k + 2, l + 1); |
---|
[806c18] | 899 | one++; |
---|
| 900 | two++; |
---|
| 901 | } |
---|
[e368746] | 902 | else if (one.hasTerms() && one.exp() == j - k + 1) |
---|
[806c18] | 903 | { |
---|
[a090c88] | 904 | tmp[l] += mulNTL ((bufFactors[l+1][k]+one.coeff()), Pi[l-1][k], b) - |
---|
[806c18] | 905 | M (k + 1, l + 1); |
---|
| 906 | one++; |
---|
| 907 | } |
---|
[e368746] | 908 | else if (two.hasTerms() && two.exp() == j - k + 1) |
---|
[806c18] | 909 | { |
---|
[a090c88] | 910 | tmp[l] += mulNTL (bufFactors[l+1][k], (Pi[l-1][k] + two.coeff()), b) |
---|
| 911 | - M (k + 1, l + 1); |
---|
[806c18] | 912 | two++; |
---|
| 913 | } |
---|
| 914 | } |
---|
| 915 | else |
---|
| 916 | tmp[l] += M (k + 1, l + 1); |
---|
| 917 | } |
---|
| 918 | } |
---|
[a090c88] | 919 | if (b.getp() != 0) |
---|
| 920 | tmp[l]= b (tmp[l]); |
---|
[ad3c3ff] | 921 | Pi[l] += tmp[l]*xToJ*F.mvar(); |
---|
| 922 | } |
---|
| 923 | } |
---|
| 924 | |
---|
[806c18] | 925 | void |
---|
| 926 | henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi, |
---|
[69fdf90] | 927 | CFList& diophant, CFMatrix& M, modpk& b, bool sort) |
---|
[ad3c3ff] | 928 | { |
---|
[e368746] | 929 | if (sort) |
---|
| 930 | sortList (factors, Variable (1)); |
---|
[ad3c3ff] | 931 | Pi= CFArray (factors.length() - 1); |
---|
| 932 | CFListIterator j= factors; |
---|
[69fdf90] | 933 | diophant= diophantine (F[0], F, factors, b); |
---|
[f9da5e] | 934 | CanonicalForm bufF= F; |
---|
| 935 | if (getCharacteristic() == 0 && b.getp() != 0) |
---|
| 936 | { |
---|
| 937 | Variable v; |
---|
| 938 | bool hasAlgVar= hasFirstAlgVar (F, v); |
---|
| 939 | for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++) |
---|
| 940 | hasAlgVar= hasFirstAlgVar (i.getItem(), v); |
---|
| 941 | Variable w; |
---|
| 942 | bool hasAlgVar2= false; |
---|
| 943 | for (CFListIterator i= diophant; i.hasItem() && !hasAlgVar2; i++) |
---|
| 944 | hasAlgVar2= hasFirstAlgVar (i.getItem(), w); |
---|
| 945 | if (hasAlgVar && hasAlgVar2 && v!=w) |
---|
| 946 | { |
---|
| 947 | bufF= replacevar (bufF, v, w); |
---|
| 948 | for (CFListIterator i= factors; i.hasItem(); i++) |
---|
| 949 | i.getItem()= replacevar (i.getItem(), v, w); |
---|
| 950 | } |
---|
| 951 | } |
---|
| 952 | |
---|
[ad3c3ff] | 953 | DEBOUTLN (cerr, "diophant= " << diophant); |
---|
| 954 | j++; |
---|
[a090c88] | 955 | Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()), b); |
---|
[ad3c3ff] | 956 | M (1, 1)= Pi [0]; |
---|
| 957 | int i= 1; |
---|
| 958 | if (j.hasItem()) |
---|
| 959 | j++; |
---|
[c1b9927] | 960 | for (; j.hasItem(); j++, i++) |
---|
[ad3c3ff] | 961 | { |
---|
[a090c88] | 962 | Pi [i]= mulNTL (Pi [i - 1], j.getItem(), b); |
---|
[ad3c3ff] | 963 | M (1, i + 1)= Pi [i]; |
---|
| 964 | } |
---|
| 965 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 966 | i= 0; |
---|
| 967 | for (CFListIterator k= factors; k.hasItem(); i++, k++) |
---|
| 968 | { |
---|
| 969 | if (i == 0) |
---|
| 970 | bufFactors[i]= mod (k.getItem(), F.mvar()); |
---|
| 971 | else |
---|
| 972 | bufFactors[i]= k.getItem(); |
---|
| 973 | } |
---|
| 974 | for (i= 1; i < l; i++) |
---|
[f9da5e] | 975 | henselStep12 (bufF, factors, bufFactors, diophant, M, Pi, i, b); |
---|
[ad3c3ff] | 976 | |
---|
| 977 | CFListIterator k= factors; |
---|
| 978 | for (i= 0; i < factors.length (); i++, k++) |
---|
| 979 | k.getItem()= bufFactors[i]; |
---|
| 980 | factors.removeFirst(); |
---|
[69fdf90] | 981 | } |
---|
| 982 | |
---|
| 983 | void |
---|
| 984 | henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi, |
---|
| 985 | CFList& diophant, CFMatrix& M, bool sort) |
---|
| 986 | { |
---|
| 987 | modpk dummy= modpk(); |
---|
| 988 | henselLift12 (F, factors, l, Pi, diophant, M, dummy, sort); |
---|
[ad3c3ff] | 989 | } |
---|
| 990 | |
---|
[806c18] | 991 | void |
---|
[ad3c3ff] | 992 | henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int |
---|
[d9357b] | 993 | end, CFArray& Pi, const CFList& diophant, CFMatrix& M, |
---|
| 994 | const modpk& b) |
---|
[ad3c3ff] | 995 | { |
---|
| 996 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 997 | int i= 0; |
---|
| 998 | CanonicalForm xToStart= power (F.mvar(), start); |
---|
| 999 | for (CFListIterator k= factors; k.hasItem(); k++, i++) |
---|
| 1000 | { |
---|
| 1001 | if (i == 0) |
---|
| 1002 | bufFactors[i]= mod (k.getItem(), xToStart); |
---|
| 1003 | else |
---|
| 1004 | bufFactors[i]= k.getItem(); |
---|
| 1005 | } |
---|
| 1006 | for (i= start; i < end; i++) |
---|
[a090c88] | 1007 | henselStep12 (F, factors, bufFactors, diophant, M, Pi, i, b); |
---|
[806c18] | 1008 | |
---|
[ad3c3ff] | 1009 | CFListIterator k= factors; |
---|
| 1010 | for (i= 0; i < factors.length(); k++, i++) |
---|
| 1011 | k.getItem()= bufFactors [i]; |
---|
| 1012 | factors.removeFirst(); |
---|
[806c18] | 1013 | return; |
---|
| 1014 | } |
---|
[ad3c3ff] | 1015 | |
---|
| 1016 | CFList |
---|
[81d96c] | 1017 | biDiophantine (const CanonicalForm& F, const CFList& factors, int d) |
---|
[ad3c3ff] | 1018 | { |
---|
| 1019 | Variable y= F.mvar(); |
---|
| 1020 | CFList result; |
---|
| 1021 | if (y.level() == 1) |
---|
| 1022 | { |
---|
| 1023 | result= diophantine (F, factors); |
---|
| 1024 | return result; |
---|
| 1025 | } |
---|
| 1026 | else |
---|
| 1027 | { |
---|
| 1028 | CFList buf= factors; |
---|
| 1029 | for (CFListIterator i= buf; i.hasItem(); i++) |
---|
| 1030 | i.getItem()= mod (i.getItem(), y); |
---|
| 1031 | CanonicalForm A= mod (F, y); |
---|
| 1032 | int bufD= 1; |
---|
| 1033 | CFList recResult= biDiophantine (A, buf, bufD); |
---|
| 1034 | CanonicalForm e= 1; |
---|
| 1035 | CFList p; |
---|
| 1036 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 1037 | CanonicalForm yToD= power (y, d); |
---|
| 1038 | int k= 0; |
---|
| 1039 | for (CFListIterator i= factors; i.hasItem(); i++, k++) |
---|
| 1040 | { |
---|
| 1041 | bufFactors [k]= i.getItem(); |
---|
| 1042 | } |
---|
[21b8f4c] | 1043 | CanonicalForm b, quot; |
---|
[ad3c3ff] | 1044 | for (k= 0; k < factors.length(); k++) //TODO compute b's faster |
---|
| 1045 | { |
---|
| 1046 | b= 1; |
---|
[21b8f4c] | 1047 | if (fdivides (bufFactors[k], F, quot)) |
---|
| 1048 | b= quot; |
---|
[806c18] | 1049 | else |
---|
[ad3c3ff] | 1050 | { |
---|
[806c18] | 1051 | for (int l= 0; l < factors.length(); l++) |
---|
| 1052 | { |
---|
| 1053 | if (l == k) |
---|
| 1054 | continue; |
---|
| 1055 | else |
---|
| 1056 | { |
---|
| 1057 | b= mulMod2 (b, bufFactors[l], yToD); |
---|
| 1058 | } |
---|
| 1059 | } |
---|
[ad3c3ff] | 1060 | } |
---|
| 1061 | p.append (b); |
---|
[806c18] | 1062 | } |
---|
[ad3c3ff] | 1063 | |
---|
| 1064 | CFListIterator j= p; |
---|
| 1065 | for (CFListIterator i= recResult; i.hasItem(); i++, j++) |
---|
| 1066 | e -= i.getItem()*j.getItem(); |
---|
| 1067 | |
---|
| 1068 | if (e.isZero()) |
---|
[806c18] | 1069 | return recResult; |
---|
[ad3c3ff] | 1070 | CanonicalForm coeffE; |
---|
| 1071 | CFList s; |
---|
| 1072 | result= recResult; |
---|
| 1073 | CanonicalForm g; |
---|
| 1074 | for (int i= 1; i < d; i++) |
---|
| 1075 | { |
---|
| 1076 | if (degree (e, y) > 0) |
---|
| 1077 | coeffE= e[i]; |
---|
| 1078 | else |
---|
| 1079 | coeffE= 0; |
---|
| 1080 | if (!coeffE.isZero()) |
---|
| 1081 | { |
---|
[806c18] | 1082 | CFListIterator k= result; |
---|
| 1083 | CFListIterator l= p; |
---|
[ad3c3ff] | 1084 | int ii= 0; |
---|
[806c18] | 1085 | j= recResult; |
---|
| 1086 | for (; j.hasItem(); j++, k++, l++, ii++) |
---|
| 1087 | { |
---|
| 1088 | g= coeffE*j.getItem(); |
---|
[ad3c3ff] | 1089 | if (degree (bufFactors[ii], y) <= 0) |
---|
| 1090 | g= mod (g, bufFactors[ii]); |
---|
| 1091 | else |
---|
| 1092 | g= mod (g, bufFactors[ii][0]); |
---|
[806c18] | 1093 | k.getItem() += g*power (y, i); |
---|
| 1094 | e -= mulMod2 (g*power(y, i), l.getItem(), yToD); |
---|
| 1095 | DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " << |
---|
[ad3c3ff] | 1096 | mod (e, power (y, i + 1))); |
---|
[806c18] | 1097 | } |
---|
| 1098 | } |
---|
[ad3c3ff] | 1099 | if (e.isZero()) |
---|
| 1100 | break; |
---|
| 1101 | } |
---|
| 1102 | |
---|
| 1103 | DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y)); |
---|
| 1104 | |
---|
| 1105 | #ifdef DEBUGOUTPUT |
---|
| 1106 | CanonicalForm test= 0; |
---|
[806c18] | 1107 | j= p; |
---|
[ad3c3ff] | 1108 | for (CFListIterator i= result; i.hasItem(); i++, j++) |
---|
| 1109 | test += mod (i.getItem()*j.getItem(), power (y, d)); |
---|
| 1110 | DEBOUTLN (cerr, "test= " << test); |
---|
| 1111 | #endif |
---|
| 1112 | return result; |
---|
| 1113 | } |
---|
| 1114 | } |
---|
| 1115 | |
---|
| 1116 | CFList |
---|
[806c18] | 1117 | multiRecDiophantine (const CanonicalForm& F, const CFList& factors, |
---|
[81d96c] | 1118 | const CFList& recResult, const CFList& M, int d) |
---|
[ad3c3ff] | 1119 | { |
---|
| 1120 | Variable y= F.mvar(); |
---|
| 1121 | CFList result; |
---|
| 1122 | CFListIterator i; |
---|
| 1123 | CanonicalForm e= 1; |
---|
| 1124 | CFListIterator j= factors; |
---|
| 1125 | CFList p; |
---|
| 1126 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 1127 | CanonicalForm yToD= power (y, d); |
---|
| 1128 | int k= 0; |
---|
| 1129 | for (CFListIterator i= factors; i.hasItem(); i++, k++) |
---|
| 1130 | bufFactors [k]= i.getItem(); |
---|
[21b8f4c] | 1131 | CanonicalForm b, quot; |
---|
[ad3c3ff] | 1132 | CFList buf= M; |
---|
| 1133 | buf.removeLast(); |
---|
| 1134 | buf.append (yToD); |
---|
| 1135 | for (k= 0; k < factors.length(); k++) //TODO compute b's faster |
---|
| 1136 | { |
---|
| 1137 | b= 1; |
---|
[21b8f4c] | 1138 | if (fdivides (bufFactors[k], F, quot)) |
---|
| 1139 | b= quot; |
---|
[806c18] | 1140 | else |
---|
[ad3c3ff] | 1141 | { |
---|
| 1142 | for (int l= 0; l < factors.length(); l++) |
---|
| 1143 | { |
---|
[806c18] | 1144 | if (l == k) |
---|
| 1145 | continue; |
---|
| 1146 | else |
---|
| 1147 | { |
---|
| 1148 | b= mulMod (b, bufFactors[l], buf); |
---|
| 1149 | } |
---|
[ad3c3ff] | 1150 | } |
---|
| 1151 | } |
---|
| 1152 | p.append (b); |
---|
[806c18] | 1153 | } |
---|
[ad3c3ff] | 1154 | j= p; |
---|
| 1155 | for (CFListIterator i= recResult; i.hasItem(); i++, j++) |
---|
| 1156 | e -= mulMod (i.getItem(), j.getItem(), M); |
---|
| 1157 | |
---|
| 1158 | if (e.isZero()) |
---|
[806c18] | 1159 | return recResult; |
---|
[ad3c3ff] | 1160 | CanonicalForm coeffE; |
---|
| 1161 | CFList s; |
---|
| 1162 | result= recResult; |
---|
| 1163 | CanonicalForm g; |
---|
| 1164 | for (int i= 1; i < d; i++) |
---|
| 1165 | { |
---|
| 1166 | if (degree (e, y) > 0) |
---|
| 1167 | coeffE= e[i]; |
---|
| 1168 | else |
---|
| 1169 | coeffE= 0; |
---|
| 1170 | if (!coeffE.isZero()) |
---|
| 1171 | { |
---|
| 1172 | CFListIterator k= result; |
---|
| 1173 | CFListIterator l= p; |
---|
| 1174 | j= recResult; |
---|
| 1175 | int ii= 0; |
---|
| 1176 | CanonicalForm dummy; |
---|
| 1177 | for (; j.hasItem(); j++, k++, l++, ii++) |
---|
| 1178 | { |
---|
[806c18] | 1179 | g= mulMod (coeffE, j.getItem(), M); |
---|
| 1180 | if (degree (bufFactors[ii], y) <= 0) |
---|
| 1181 | divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy, |
---|
[ad3c3ff] | 1182 | g, M); |
---|
[806c18] | 1183 | else |
---|
| 1184 | divrem (g, bufFactors[ii][0], dummy, g, M); |
---|
| 1185 | k.getItem() += g*power (y, i); |
---|
| 1186 | e -= mulMod (g*power (y, i), l.getItem(), M); |
---|
| 1187 | } |
---|
[ad3c3ff] | 1188 | } |
---|
[806c18] | 1189 | |
---|
[ad3c3ff] | 1190 | if (e.isZero()) |
---|
| 1191 | break; |
---|
| 1192 | } |
---|
| 1193 | |
---|
| 1194 | #ifdef DEBUGOUTPUT |
---|
| 1195 | CanonicalForm test= 0; |
---|
[806c18] | 1196 | j= p; |
---|
[ad3c3ff] | 1197 | for (CFListIterator i= result; i.hasItem(); i++, j++) |
---|
| 1198 | test += mod (i.getItem()*j.getItem(), power (y, d)); |
---|
[27ab36] | 1199 | DEBOUTLN (cerr, "test in multiRecDiophantine= " << test); |
---|
[ad3c3ff] | 1200 | #endif |
---|
| 1201 | return result; |
---|
| 1202 | } |
---|
| 1203 | |
---|
[81d96c] | 1204 | static |
---|
[806c18] | 1205 | void |
---|
| 1206 | henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors, |
---|
| 1207 | const CFList& diophant, CFMatrix& M, CFArray& Pi, int j, |
---|
| 1208 | const CFList& MOD) |
---|
[ad3c3ff] | 1209 | { |
---|
[c4f4fd] | 1210 | CanonicalForm E; |
---|
| 1211 | CanonicalForm xToJ= power (F.mvar(), j); |
---|
| 1212 | Variable x= F.mvar(); |
---|
| 1213 | // compute the error |
---|
| 1214 | if (j == 1) |
---|
[ad3c3ff] | 1215 | { |
---|
[c4f4fd] | 1216 | E= F[j]; |
---|
[806c18] | 1217 | #ifdef DEBUGOUTPUT |
---|
[c4f4fd] | 1218 | CanonicalForm test= 1; |
---|
| 1219 | for (int i= 0; i < factors.length(); i++) |
---|
[ad3c3ff] | 1220 | { |
---|
[c4f4fd] | 1221 | if (i == 0) |
---|
| 1222 | test= mulMod (test, mod (bufFactors [i], xToJ), MOD); |
---|
| 1223 | else |
---|
| 1224 | test= mulMod (test, bufFactors[i], MOD); |
---|
[ad3c3ff] | 1225 | } |
---|
[c4f4fd] | 1226 | CanonicalForm test2= mod (F-test, xToJ); |
---|
[139b67] | 1227 | |
---|
[c4f4fd] | 1228 | test2= mod (test2, MOD); |
---|
[27ab36] | 1229 | DEBOUTLN (cerr, "test in henselStep= " << test2); |
---|
[c4f4fd] | 1230 | #endif |
---|
[806c18] | 1231 | } |
---|
[c4f4fd] | 1232 | else |
---|
[ad3c3ff] | 1233 | { |
---|
[806c18] | 1234 | #ifdef DEBUGOUTPUT |
---|
[c4f4fd] | 1235 | CanonicalForm test= 1; |
---|
| 1236 | for (int i= 0; i < factors.length(); i++) |
---|
[ad3c3ff] | 1237 | { |
---|
[c4f4fd] | 1238 | if (i == 0) |
---|
| 1239 | test *= mod (bufFactors [i], power (x, j)); |
---|
| 1240 | else |
---|
| 1241 | test *= bufFactors[i]; |
---|
[ad3c3ff] | 1242 | } |
---|
[c4f4fd] | 1243 | test= mod (test, power (x, j)); |
---|
| 1244 | test= mod (test, MOD); |
---|
| 1245 | CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1)); |
---|
[27ab36] | 1246 | DEBOUTLN (cerr, "test in henselStep= " << test2); |
---|
[c4f4fd] | 1247 | #endif |
---|
| 1248 | |
---|
| 1249 | if (degree (Pi [factors.length() - 2], x) > 0) |
---|
| 1250 | E= F[j] - Pi [factors.length() - 2] [j]; |
---|
| 1251 | else |
---|
| 1252 | E= F[j]; |
---|
[139b67] | 1253 | } |
---|
| 1254 | |
---|
| 1255 | CFArray buf= CFArray (diophant.length()); |
---|
| 1256 | bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1)); |
---|
| 1257 | int k= 0; |
---|
| 1258 | // actual lifting |
---|
| 1259 | CanonicalForm dummy, rest1; |
---|
[806c18] | 1260 | for (CFListIterator i= diophant; i.hasItem(); i++, k++) |
---|
[139b67] | 1261 | { |
---|
[ad3c3ff] | 1262 | if (degree (bufFactors[k], x) > 0) |
---|
| 1263 | { |
---|
| 1264 | if (k > 0) |
---|
| 1265 | divrem (E, bufFactors[k] [0], dummy, rest1, MOD); |
---|
| 1266 | else |
---|
| 1267 | rest1= E; |
---|
| 1268 | } |
---|
| 1269 | else |
---|
| 1270 | divrem (E, bufFactors[k], dummy, rest1, MOD); |
---|
| 1271 | |
---|
| 1272 | buf[k]= mulMod (i.getItem(), rest1, MOD); |
---|
| 1273 | |
---|
| 1274 | if (degree (bufFactors[k], x) > 0) |
---|
| 1275 | divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD); |
---|
[806c18] | 1276 | else |
---|
[ad3c3ff] | 1277 | divrem (buf[k], bufFactors[k], dummy, buf[k], MOD); |
---|
| 1278 | } |
---|
| 1279 | for (k= 1; k < factors.length(); k++) |
---|
| 1280 | bufFactors[k] += xToJ*buf[k]; |
---|
| 1281 | |
---|
| 1282 | // update Pi [0] |
---|
| 1283 | int degBuf0= degree (bufFactors[0], x); |
---|
| 1284 | int degBuf1= degree (bufFactors[1], x); |
---|
| 1285 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 1286 | M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD); |
---|
| 1287 | CanonicalForm uIZeroJ; |
---|
| 1288 | if (j == 1) |
---|
| 1289 | { |
---|
| 1290 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 1291 | uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]), |
---|
| 1292 | (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1); |
---|
| 1293 | else if (degBuf0 > 0) |
---|
| 1294 | uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD); |
---|
| 1295 | else if (degBuf1 > 0) |
---|
| 1296 | uIZeroJ= mulMod (bufFactors[0], buf[1], MOD); |
---|
| 1297 | else |
---|
| 1298 | uIZeroJ= 0; |
---|
[806c18] | 1299 | Pi [0] += xToJ*uIZeroJ; |
---|
| 1300 | } |
---|
[ad3c3ff] | 1301 | else |
---|
| 1302 | { |
---|
| 1303 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
[806c18] | 1304 | uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]), |
---|
[ad3c3ff] | 1305 | (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1); |
---|
| 1306 | else if (degBuf0 > 0) |
---|
| 1307 | uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD); |
---|
| 1308 | else if (degBuf1 > 0) |
---|
| 1309 | uIZeroJ= mulMod (bufFactors[0], buf[1], MOD); |
---|
| 1310 | else |
---|
| 1311 | uIZeroJ= 0; |
---|
[806c18] | 1312 | Pi [0] += xToJ*uIZeroJ; |
---|
[ad3c3ff] | 1313 | } |
---|
| 1314 | |
---|
| 1315 | CFArray tmp= CFArray (factors.length() - 1); |
---|
| 1316 | for (k= 0; k < factors.length() - 1; k++) |
---|
| 1317 | tmp[k]= 0; |
---|
| 1318 | CFIterator one, two; |
---|
| 1319 | one= bufFactors [0]; |
---|
| 1320 | two= bufFactors [1]; |
---|
| 1321 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 1322 | { |
---|
| 1323 | for (k= 1; k <= (int) ceil (j/2.0); k++) |
---|
[806c18] | 1324 | { |
---|
[ad3c3ff] | 1325 | if (k != j - k + 1) |
---|
| 1326 | { |
---|
[c1b9927] | 1327 | if ((one.hasTerms() && one.exp() == j - k + 1) && |
---|
[e368746] | 1328 | (two.hasTerms() && two.exp() == j - k + 1)) |
---|
[806c18] | 1329 | { |
---|
| 1330 | tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()), |
---|
| 1331 | (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) - |
---|
[ad3c3ff] | 1332 | M (j - k + 2, 1); |
---|
[806c18] | 1333 | one++; |
---|
| 1334 | two++; |
---|
| 1335 | } |
---|
[e368746] | 1336 | else if (one.hasTerms() && one.exp() == j - k + 1) |
---|
[806c18] | 1337 | { |
---|
| 1338 | tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()), |
---|
[ad3c3ff] | 1339 | bufFactors[1] [k], MOD) - M (k + 1, 1); |
---|
[806c18] | 1340 | one++; |
---|
| 1341 | } |
---|
[e368746] | 1342 | else if (two.hasTerms() && two.exp() == j - k + 1) |
---|
[806c18] | 1343 | { |
---|
| 1344 | tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] + |
---|
[ad3c3ff] | 1345 | two.coeff()), MOD) - M (k + 1, 1); |
---|
[806c18] | 1346 | two++; |
---|
| 1347 | } |
---|
[ad3c3ff] | 1348 | } |
---|
| 1349 | else |
---|
| 1350 | { |
---|
[806c18] | 1351 | tmp[0] += M (k + 1, 1); |
---|
[ad3c3ff] | 1352 | } |
---|
| 1353 | } |
---|
| 1354 | } |
---|
| 1355 | Pi [0] += tmp[0]*xToJ*F.mvar(); |
---|
| 1356 | |
---|
[806c18] | 1357 | // update Pi [l] |
---|
[ad3c3ff] | 1358 | int degPi, degBuf; |
---|
[806c18] | 1359 | for (int l= 1; l < factors.length() - 1; l++) |
---|
[ad3c3ff] | 1360 | { |
---|
| 1361 | degPi= degree (Pi [l - 1], x); |
---|
| 1362 | degBuf= degree (bufFactors[l + 1], x); |
---|
| 1363 | if (degPi > 0 && degBuf > 0) |
---|
[806c18] | 1364 | M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD); |
---|
[ad3c3ff] | 1365 | if (j == 1) |
---|
| 1366 | { |
---|
| 1367 | if (degPi > 0 && degBuf > 0) |
---|
[806c18] | 1368 | Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]), |
---|
| 1369 | (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)- |
---|
[ad3c3ff] | 1370 | M (1, l + 1)); |
---|
| 1371 | else if (degPi > 0) |
---|
| 1372 | Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD)); |
---|
| 1373 | else if (degBuf > 0) |
---|
| 1374 | Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD)); |
---|
| 1375 | } |
---|
| 1376 | else |
---|
| 1377 | { |
---|
| 1378 | if (degPi > 0 && degBuf > 0) |
---|
| 1379 | { |
---|
| 1380 | uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD); |
---|
[806c18] | 1381 | uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD); |
---|
[ad3c3ff] | 1382 | } |
---|
| 1383 | else if (degPi > 0) |
---|
| 1384 | uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD); |
---|
| 1385 | else if (degBuf > 0) |
---|
| 1386 | { |
---|
| 1387 | uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD); |
---|
| 1388 | uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD); |
---|
| 1389 | } |
---|
| 1390 | Pi[l] += xToJ*uIZeroJ; |
---|
| 1391 | } |
---|
| 1392 | one= bufFactors [l + 1]; |
---|
| 1393 | two= Pi [l - 1]; |
---|
[e368746] | 1394 | if (two.hasTerms() && two.exp() == j + 1) |
---|
[ad3c3ff] | 1395 | { |
---|
| 1396 | if (degBuf > 0 && degPi > 0) |
---|
| 1397 | { |
---|
[806c18] | 1398 | tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD); |
---|
| 1399 | two++; |
---|
[ad3c3ff] | 1400 | } |
---|
| 1401 | else if (degPi > 0) |
---|
[806c18] | 1402 | { |
---|
| 1403 | tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD); |
---|
| 1404 | two++; |
---|
| 1405 | } |
---|
[ad3c3ff] | 1406 | } |
---|
[806c18] | 1407 | if (degBuf > 0 && degPi > 0) |
---|
[ad3c3ff] | 1408 | { |
---|
| 1409 | for (k= 1; k <= (int) ceil (j/2.0); k++) |
---|
[806c18] | 1410 | { |
---|
| 1411 | if (k != j - k + 1) |
---|
| 1412 | { |
---|
[c1b9927] | 1413 | if ((one.hasTerms() && one.exp() == j - k + 1) && |
---|
[e368746] | 1414 | (two.hasTerms() && two.exp() == j - k + 1)) |
---|
[806c18] | 1415 | { |
---|
| 1416 | tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()), |
---|
| 1417 | (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) - |
---|
[ad3c3ff] | 1418 | M (j - k + 2, l + 1); |
---|
[806c18] | 1419 | one++; |
---|
| 1420 | two++; |
---|
| 1421 | } |
---|
[e368746] | 1422 | else if (one.hasTerms() && one.exp() == j - k + 1) |
---|
[806c18] | 1423 | { |
---|
| 1424 | tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()), |
---|
| 1425 | Pi[l - 1] [k], MOD) - M (k + 1, l + 1); |
---|
| 1426 | one++; |
---|
| 1427 | } |
---|
[e368746] | 1428 | else if (two.hasTerms() && two.exp() == j - k + 1) |
---|
[806c18] | 1429 | { |
---|
| 1430 | tmp[l] += mulMod (bufFactors[l + 1] [k], |
---|
[ad3c3ff] | 1431 | (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1); |
---|
[806c18] | 1432 | two++; |
---|
| 1433 | } |
---|
| 1434 | } |
---|
| 1435 | else |
---|
| 1436 | tmp[l] += M (k + 1, l + 1); |
---|
| 1437 | } |
---|
| 1438 | } |
---|
[ad3c3ff] | 1439 | Pi[l] += tmp[l]*xToJ*F.mvar(); |
---|
| 1440 | } |
---|
| 1441 | |
---|
| 1442 | return; |
---|
| 1443 | } |
---|
| 1444 | |
---|
[806c18] | 1445 | CFList |
---|
[81d96c] | 1446 | henselLift23 (const CFList& eval, const CFList& factors, int* l, CFList& |
---|
[ad3c3ff] | 1447 | diophant, CFArray& Pi, CFMatrix& M) |
---|
| 1448 | { |
---|
| 1449 | CFList buf= factors; |
---|
| 1450 | int k= 0; |
---|
| 1451 | int liftBoundBivar= l[k]; |
---|
| 1452 | diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar); |
---|
| 1453 | CFList MOD; |
---|
| 1454 | MOD.append (power (Variable (2), liftBoundBivar)); |
---|
| 1455 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 1456 | k= 0; |
---|
| 1457 | CFListIterator j= eval; |
---|
| 1458 | j++; |
---|
| 1459 | buf.removeFirst(); |
---|
| 1460 | buf.insert (LC (j.getItem(), 1)); |
---|
| 1461 | for (CFListIterator i= buf; i.hasItem(); i++, k++) |
---|
| 1462 | bufFactors[k]= i.getItem(); |
---|
| 1463 | Pi= CFArray (factors.length() - 1); |
---|
| 1464 | CFListIterator i= buf; |
---|
| 1465 | i++; |
---|
| 1466 | Variable y= j.getItem().mvar(); |
---|
| 1467 | Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD); |
---|
| 1468 | M (1, 1)= Pi [0]; |
---|
| 1469 | k= 1; |
---|
| 1470 | if (i.hasItem()) |
---|
| 1471 | i++; |
---|
[c1b9927] | 1472 | for (; i.hasItem(); i++, k++) |
---|
[ad3c3ff] | 1473 | { |
---|
| 1474 | Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD); |
---|
| 1475 | M (1, k + 1)= Pi [k]; |
---|
| 1476 | } |
---|
| 1477 | |
---|
| 1478 | for (int d= 1; d < l[1]; d++) |
---|
| 1479 | henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD); |
---|
| 1480 | CFList result; |
---|
| 1481 | for (k= 1; k < factors.length(); k++) |
---|
| 1482 | result.append (bufFactors[k]); |
---|
| 1483 | return result; |
---|
[806c18] | 1484 | } |
---|
[ad3c3ff] | 1485 | |
---|
[806c18] | 1486 | void |
---|
| 1487 | henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end, |
---|
| 1488 | CFArray& Pi, const CFList& diophant, CFMatrix& M, |
---|
[ad3c3ff] | 1489 | const CFList& MOD) |
---|
| 1490 | { |
---|
| 1491 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 1492 | int i= 0; |
---|
| 1493 | CanonicalForm xToStart= power (F.mvar(), start); |
---|
| 1494 | for (CFListIterator k= factors; k.hasItem(); k++, i++) |
---|
| 1495 | { |
---|
| 1496 | if (i == 0) |
---|
| 1497 | bufFactors[i]= mod (k.getItem(), xToStart); |
---|
| 1498 | else |
---|
| 1499 | bufFactors[i]= k.getItem(); |
---|
| 1500 | } |
---|
| 1501 | for (i= start; i < end; i++) |
---|
[806c18] | 1502 | henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD); |
---|
| 1503 | |
---|
[ad3c3ff] | 1504 | CFListIterator k= factors; |
---|
| 1505 | for (i= 0; i < factors.length(); k++, i++) |
---|
| 1506 | k.getItem()= bufFactors [i]; |
---|
| 1507 | factors.removeFirst(); |
---|
[806c18] | 1508 | return; |
---|
| 1509 | } |
---|
[ad3c3ff] | 1510 | |
---|
| 1511 | CFList |
---|
| 1512 | henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList& |
---|
[81d96c] | 1513 | diophant, CFArray& Pi, CFMatrix& M, int lOld, int lNew) |
---|
[ad3c3ff] | 1514 | { |
---|
| 1515 | diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld); |
---|
| 1516 | int k= 0; |
---|
| 1517 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 1518 | for (CFListIterator i= factors; i.hasItem(); i++, k++) |
---|
| 1519 | { |
---|
| 1520 | if (k == 0) |
---|
| 1521 | bufFactors[k]= LC (F.getLast(), 1); |
---|
| 1522 | else |
---|
| 1523 | bufFactors[k]= i.getItem(); |
---|
| 1524 | } |
---|
| 1525 | CFList buf= factors; |
---|
| 1526 | buf.removeFirst(); |
---|
| 1527 | buf.insert (LC (F.getLast(), 1)); |
---|
| 1528 | CFListIterator i= buf; |
---|
| 1529 | i++; |
---|
| 1530 | Variable y= F.getLast().mvar(); |
---|
| 1531 | Variable x= F.getFirst().mvar(); |
---|
| 1532 | CanonicalForm xToLOld= power (x, lOld); |
---|
| 1533 | Pi [0]= mod (Pi[0], xToLOld); |
---|
| 1534 | M (1, 1)= Pi [0]; |
---|
| 1535 | k= 1; |
---|
| 1536 | if (i.hasItem()) |
---|
| 1537 | i++; |
---|
[c1b9927] | 1538 | for (; i.hasItem(); i++, k++) |
---|
[ad3c3ff] | 1539 | { |
---|
| 1540 | Pi [k]= mod (Pi [k], xToLOld); |
---|
| 1541 | M (1, k + 1)= Pi [k]; |
---|
| 1542 | } |
---|
| 1543 | |
---|
| 1544 | for (int d= 1; d < lNew; d++) |
---|
| 1545 | henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD); |
---|
| 1546 | CFList result; |
---|
| 1547 | for (k= 1; k < factors.length(); k++) |
---|
| 1548 | result.append (bufFactors[k]); |
---|
| 1549 | return result; |
---|
| 1550 | } |
---|
| 1551 | |
---|
| 1552 | CFList |
---|
[81d96c] | 1553 | henselLift (const CFList& eval, const CFList& factors, int* l, int lLength, |
---|
| 1554 | bool sort) |
---|
[ad3c3ff] | 1555 | { |
---|
[806c18] | 1556 | CFList diophant; |
---|
| 1557 | CFList buf= factors; |
---|
[ad3c3ff] | 1558 | buf.insert (LC (eval.getFirst(), 1)); |
---|
[e368746] | 1559 | if (sort) |
---|
| 1560 | sortList (buf, Variable (1)); |
---|
[ad3c3ff] | 1561 | CFArray Pi; |
---|
| 1562 | CFMatrix M= CFMatrix (l[1], factors.length()); |
---|
| 1563 | CFList result= henselLift23 (eval, buf, l, diophant, Pi, M); |
---|
| 1564 | if (eval.length() == 2) |
---|
| 1565 | return result; |
---|
| 1566 | CFList MOD; |
---|
| 1567 | for (int i= 0; i < 2; i++) |
---|
| 1568 | MOD.append (power (Variable (i + 2), l[i])); |
---|
| 1569 | CFListIterator j= eval; |
---|
| 1570 | j++; |
---|
| 1571 | CFList bufEval; |
---|
| 1572 | bufEval.append (j.getItem()); |
---|
| 1573 | j++; |
---|
[806c18] | 1574 | |
---|
[ea88e0] | 1575 | for (int i= 2; i < lLength && j.hasItem(); i++, j++) |
---|
[ad3c3ff] | 1576 | { |
---|
| 1577 | result.insert (LC (bufEval.getFirst(), 1)); |
---|
| 1578 | bufEval.append (j.getItem()); |
---|
| 1579 | M= CFMatrix (l[i], factors.length()); |
---|
| 1580 | result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]); |
---|
| 1581 | MOD.append (power (Variable (i + 2), l[i])); |
---|
| 1582 | bufEval.removeFirst(); |
---|
| 1583 | } |
---|
[806c18] | 1584 | return result; |
---|
[ad3c3ff] | 1585 | } |
---|
| 1586 | |
---|
[81d96c] | 1587 | // nonmonic |
---|
| 1588 | |
---|
[08daea] | 1589 | void |
---|
[81d96c] | 1590 | nonMonicHenselStep12 (const CanonicalForm& F, const CFList& factors, |
---|
| 1591 | CFArray& bufFactors, const CFList& diophant, CFMatrix& M, |
---|
| 1592 | CFArray& Pi, int j, const CFArray& /*LCs*/) |
---|
[08daea] | 1593 | { |
---|
| 1594 | Variable x= F.mvar(); |
---|
| 1595 | CanonicalForm xToJ= power (x, j); |
---|
| 1596 | |
---|
| 1597 | CanonicalForm E; |
---|
| 1598 | // compute the error |
---|
| 1599 | if (degree (Pi [factors.length() - 2], x) > 0) |
---|
| 1600 | E= F[j] - Pi [factors.length() - 2] [j]; |
---|
| 1601 | else |
---|
| 1602 | E= F[j]; |
---|
| 1603 | |
---|
| 1604 | CFArray buf= CFArray (diophant.length()); |
---|
| 1605 | |
---|
| 1606 | int k= 0; |
---|
| 1607 | CanonicalForm remainder; |
---|
| 1608 | // actual lifting |
---|
| 1609 | for (CFListIterator i= diophant; i.hasItem(); i++, k++) |
---|
| 1610 | { |
---|
| 1611 | if (degree (bufFactors[k], x) > 0) |
---|
| 1612 | remainder= modNTL (E, bufFactors[k] [0]); |
---|
| 1613 | else |
---|
| 1614 | remainder= modNTL (E, bufFactors[k]); |
---|
| 1615 | buf[k]= mulNTL (i.getItem(), remainder); |
---|
| 1616 | if (degree (bufFactors[k], x) > 0) |
---|
| 1617 | buf[k]= modNTL (buf[k], bufFactors[k] [0]); |
---|
| 1618 | else |
---|
[c1b9927] | 1619 | buf[k]= modNTL (buf[k], bufFactors[k]); |
---|
[08daea] | 1620 | } |
---|
| 1621 | |
---|
| 1622 | for (k= 0; k < factors.length(); k++) |
---|
| 1623 | bufFactors[k] += xToJ*buf[k]; |
---|
| 1624 | |
---|
| 1625 | // update Pi [0] |
---|
| 1626 | int degBuf0= degree (bufFactors[0], x); |
---|
| 1627 | int degBuf1= degree (bufFactors[1], x); |
---|
| 1628 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 1629 | { |
---|
| 1630 | M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]); |
---|
| 1631 | if (j + 2 <= M.rows()) |
---|
| 1632 | M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]); |
---|
| 1633 | } |
---|
| 1634 | CanonicalForm uIZeroJ; |
---|
| 1635 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 1636 | uIZeroJ= mulNTL(bufFactors[0][0],buf[1])+mulNTL (bufFactors[1][0], buf[0]); |
---|
| 1637 | else if (degBuf0 > 0) |
---|
| 1638 | uIZeroJ= mulNTL (buf[0], bufFactors[1]); |
---|
| 1639 | else if (degBuf1 > 0) |
---|
| 1640 | uIZeroJ= mulNTL (bufFactors[0], buf [1]); |
---|
| 1641 | else |
---|
| 1642 | uIZeroJ= 0; |
---|
| 1643 | Pi [0] += xToJ*uIZeroJ; |
---|
| 1644 | |
---|
| 1645 | CFArray tmp= CFArray (factors.length() - 1); |
---|
| 1646 | for (k= 0; k < factors.length() - 1; k++) |
---|
| 1647 | tmp[k]= 0; |
---|
| 1648 | CFIterator one, two; |
---|
| 1649 | one= bufFactors [0]; |
---|
| 1650 | two= bufFactors [1]; |
---|
| 1651 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 1652 | { |
---|
[1b8e048] | 1653 | while (one.hasTerms() && one.exp() > j) one++; |
---|
| 1654 | while (two.hasTerms() && two.exp() > j) two++; |
---|
[08daea] | 1655 | for (k= 1; k <= (int) ceil (j/2.0); k++) |
---|
| 1656 | { |
---|
[e368746] | 1657 | if (one.hasTerms() && two.hasTerms()) |
---|
[08daea] | 1658 | { |
---|
[e368746] | 1659 | if (k != j - k + 1) |
---|
| 1660 | { |
---|
| 1661 | if ((one.hasTerms() && one.exp() == j - k + 1) && + |
---|
| 1662 | (two.hasTerms() && two.exp() == j - k + 1)) |
---|
[08daea] | 1663 | { |
---|
| 1664 | tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] + |
---|
| 1665 | two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1); |
---|
| 1666 | one++; |
---|
| 1667 | two++; |
---|
| 1668 | } |
---|
[e368746] | 1669 | else if (one.hasTerms() && one.exp() == j - k + 1) |
---|
[08daea] | 1670 | { |
---|
| 1671 | tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) - |
---|
| 1672 | M (k + 1, 1); |
---|
| 1673 | one++; |
---|
| 1674 | } |
---|
[e368746] | 1675 | else if (two.hasTerms() && two.exp() == j - k + 1) |
---|
[08daea] | 1676 | { |
---|
| 1677 | tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) - |
---|
| 1678 | M (k + 1, 1); |
---|
| 1679 | two++; |
---|
| 1680 | } |
---|
| 1681 | } |
---|
| 1682 | else |
---|
| 1683 | tmp[0] += M (k + 1, 1); |
---|
[e368746] | 1684 | } |
---|
[08daea] | 1685 | } |
---|
[1b8e048] | 1686 | } |
---|
[08daea] | 1687 | |
---|
[1b8e048] | 1688 | if (degBuf0 >= j + 1 && degBuf1 >= j + 1) |
---|
| 1689 | { |
---|
[08daea] | 1690 | if (j + 2 <= M.rows()) |
---|
[1b8e048] | 1691 | tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]), |
---|
| 1692 | (bufFactors [1] [j + 1] + bufFactors [1] [0])) |
---|
| 1693 | - M(1,1) - M (j + 2,1); |
---|
| 1694 | } |
---|
| 1695 | else if (degBuf0 >= j + 1) |
---|
| 1696 | { |
---|
| 1697 | if (degBuf1 > 0) |
---|
| 1698 | tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]); |
---|
| 1699 | else |
---|
| 1700 | tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1]); |
---|
| 1701 | } |
---|
| 1702 | else if (degBuf1 >= j + 1) |
---|
| 1703 | { |
---|
| 1704 | if (degBuf0 > 0) |
---|
| 1705 | tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]); |
---|
| 1706 | else |
---|
| 1707 | tmp[0] += mulNTL (bufFactors [0], bufFactors [1] [j + 1]); |
---|
[08daea] | 1708 | } |
---|
[1b8e048] | 1709 | |
---|
[08daea] | 1710 | Pi [0] += tmp[0]*xToJ*F.mvar(); |
---|
[327efa2] | 1711 | |
---|
| 1712 | int degPi, degBuf; |
---|
| 1713 | for (int l= 1; l < factors.length() - 1; l++) |
---|
| 1714 | { |
---|
| 1715 | degPi= degree (Pi [l - 1], x); |
---|
| 1716 | degBuf= degree (bufFactors[l + 1], x); |
---|
| 1717 | if (degPi > 0 && degBuf > 0) |
---|
| 1718 | { |
---|
| 1719 | M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]); |
---|
| 1720 | if (j + 2 <= M.rows()) |
---|
| 1721 | M (j + 2, l + 1)= mulNTL (Pi [l - 1][j + 1], bufFactors[l + 1] [j + 1]); |
---|
| 1722 | } |
---|
| 1723 | |
---|
| 1724 | if (degPi > 0 && degBuf > 0) |
---|
| 1725 | uIZeroJ= mulNTL (Pi[l -1] [0], buf[l + 1]) + |
---|
| 1726 | mulNTL (uIZeroJ, bufFactors[l+1] [0]); |
---|
| 1727 | else if (degPi > 0) |
---|
| 1728 | uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1]); |
---|
| 1729 | else if (degBuf > 0) |
---|
| 1730 | uIZeroJ= mulNTL (Pi[l - 1], buf[1]); |
---|
| 1731 | else |
---|
| 1732 | uIZeroJ= 0; |
---|
| 1733 | |
---|
| 1734 | Pi [l] += xToJ*uIZeroJ; |
---|
| 1735 | |
---|
| 1736 | one= bufFactors [l + 1]; |
---|
| 1737 | two= Pi [l - 1]; |
---|
| 1738 | if (degBuf > 0 && degPi > 0) |
---|
| 1739 | { |
---|
| 1740 | while (one.hasTerms() && one.exp() > j) one++; |
---|
| 1741 | while (two.hasTerms() && two.exp() > j) two++; |
---|
| 1742 | for (k= 1; k <= (int) ceil (j/2.0); k++) |
---|
| 1743 | { |
---|
| 1744 | if (k != j - k + 1) |
---|
| 1745 | { |
---|
| 1746 | if ((one.hasTerms() && one.exp() == j - k + 1) && |
---|
| 1747 | (two.hasTerms() && two.exp() == j - k + 1)) |
---|
| 1748 | { |
---|
| 1749 | tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), |
---|
| 1750 | (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1) - |
---|
| 1751 | M (j - k + 2, l + 1); |
---|
| 1752 | one++; |
---|
| 1753 | two++; |
---|
| 1754 | } |
---|
| 1755 | else if (one.hasTerms() && one.exp() == j - k + 1) |
---|
| 1756 | { |
---|
| 1757 | tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), |
---|
| 1758 | Pi[l - 1] [k]) - M (k + 1, l + 1); |
---|
| 1759 | one++; |
---|
| 1760 | } |
---|
| 1761 | else if (two.hasTerms() && two.exp() == j - k + 1) |
---|
| 1762 | { |
---|
| 1763 | tmp[l] += mulNTL (bufFactors[l + 1] [k], |
---|
| 1764 | (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1); |
---|
| 1765 | two++; |
---|
| 1766 | } |
---|
| 1767 | } |
---|
| 1768 | else |
---|
| 1769 | tmp[l] += M (k + 1, l + 1); |
---|
| 1770 | } |
---|
| 1771 | } |
---|
| 1772 | |
---|
| 1773 | if (degPi >= j + 1 && degBuf >= j + 1) |
---|
| 1774 | { |
---|
| 1775 | if (j + 2 <= M.rows()) |
---|
| 1776 | tmp [l] += mulNTL ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]), |
---|
| 1777 | (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0]) |
---|
| 1778 | ) - M(1,l+1) - M (j + 2,l+1); |
---|
| 1779 | } |
---|
| 1780 | else if (degPi >= j + 1) |
---|
| 1781 | { |
---|
| 1782 | if (degBuf > 0) |
---|
| 1783 | tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1] [0]); |
---|
| 1784 | else |
---|
| 1785 | tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1]); |
---|
| 1786 | } |
---|
| 1787 | else if (degBuf >= j + 1) |
---|
| 1788 | { |
---|
| 1789 | if (degPi > 0) |
---|
| 1790 | tmp[l] += mulNTL (Pi [l - 1] [0], bufFactors [l + 1] [j + 1]); |
---|
| 1791 | else |
---|
| 1792 | tmp[l] += mulNTL (Pi [l - 1], bufFactors [l + 1] [j + 1]); |
---|
| 1793 | } |
---|
| 1794 | |
---|
| 1795 | Pi[l] += tmp[l]*xToJ*F.mvar(); |
---|
| 1796 | } |
---|
[08daea] | 1797 | return; |
---|
| 1798 | } |
---|
| 1799 | |
---|
| 1800 | void |
---|
[81d96c] | 1801 | nonMonicHenselLift12 (const CanonicalForm& F, CFList& factors, int l, |
---|
| 1802 | CFArray& Pi, CFList& diophant, CFMatrix& M, |
---|
| 1803 | const CFArray& LCs, bool sort) |
---|
[08daea] | 1804 | { |
---|
| 1805 | if (sort) |
---|
| 1806 | sortList (factors, Variable (1)); |
---|
| 1807 | Pi= CFArray (factors.length() - 2); |
---|
| 1808 | CFList bufFactors2= factors; |
---|
| 1809 | bufFactors2.removeFirst(); |
---|
[327efa2] | 1810 | diophant= diophantine (F[0], bufFactors2); |
---|
[08daea] | 1811 | DEBOUTLN (cerr, "diophant= " << diophant); |
---|
| 1812 | |
---|
| 1813 | CFArray bufFactors= CFArray (bufFactors2.length()); |
---|
| 1814 | int i= 0; |
---|
| 1815 | for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++) |
---|
| 1816 | bufFactors[i]= replaceLc (k.getItem(), LCs [i]); |
---|
| 1817 | |
---|
| 1818 | Variable x= F.mvar(); |
---|
| 1819 | if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0) |
---|
| 1820 | { |
---|
| 1821 | M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]); |
---|
| 1822 | Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) + |
---|
| 1823 | mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x; |
---|
| 1824 | } |
---|
| 1825 | else if (degree (bufFactors[0], x) > 0) |
---|
| 1826 | { |
---|
| 1827 | M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]); |
---|
| 1828 | Pi [0]= M (1, 1) + |
---|
| 1829 | mulNTL (bufFactors [0] [1], bufFactors[1])*x; |
---|
| 1830 | } |
---|
| 1831 | else if (degree (bufFactors[1], x) > 0) |
---|
| 1832 | { |
---|
| 1833 | M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]); |
---|
| 1834 | Pi [0]= M (1, 1) + |
---|
| 1835 | mulNTL (bufFactors [0], bufFactors[1] [1])*x; |
---|
| 1836 | } |
---|
| 1837 | else |
---|
| 1838 | { |
---|
| 1839 | M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]); |
---|
| 1840 | Pi [0]= M (1, 1); |
---|
| 1841 | } |
---|
| 1842 | |
---|
[327efa2] | 1843 | for (i= 1; i < Pi.size(); i++) |
---|
| 1844 | { |
---|
| 1845 | if (degree (Pi[i-1], x) > 0 && degree (bufFactors [i+1], x) > 0) |
---|
| 1846 | { |
---|
| 1847 | M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors[i+1] [0]); |
---|
| 1848 | Pi [i]= M (1,i+1) + (mulNTL (Pi[i-1] [1], bufFactors[i+1] [0]) + |
---|
| 1849 | mulNTL (Pi[i-1] [0], bufFactors [i+1] [1]))*x; |
---|
| 1850 | } |
---|
| 1851 | else if (degree (Pi[i-1], x) > 0) |
---|
| 1852 | { |
---|
| 1853 | M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors [i+1]); |
---|
| 1854 | Pi [i]= M(1,i+1) + mulNTL (Pi[i-1] [1], bufFactors[i+1])*x; |
---|
| 1855 | } |
---|
| 1856 | else if (degree (bufFactors[i+1], x) > 0) |
---|
| 1857 | { |
---|
| 1858 | M (1,i+1)= mulNTL (Pi[i-1], bufFactors [i+1] [0]); |
---|
| 1859 | Pi [i]= M (1,i+1) + mulNTL (Pi[i-1], bufFactors[i+1] [1])*x; |
---|
| 1860 | } |
---|
| 1861 | else |
---|
| 1862 | { |
---|
| 1863 | M (1,i+1)= mulNTL (Pi [i-1], bufFactors [i+1]); |
---|
| 1864 | Pi [i]= M (1,i+1); |
---|
| 1865 | } |
---|
| 1866 | } |
---|
| 1867 | |
---|
[08daea] | 1868 | for (i= 1; i < l; i++) |
---|
[81d96c] | 1869 | nonMonicHenselStep12 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs); |
---|
[08daea] | 1870 | |
---|
| 1871 | factors= CFList(); |
---|
| 1872 | for (i= 0; i < bufFactors.size(); i++) |
---|
| 1873 | factors.append (bufFactors[i]); |
---|
| 1874 | return; |
---|
| 1875 | } |
---|
| 1876 | |
---|
[327efa2] | 1877 | |
---|
[c1b9927] | 1878 | /// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$ |
---|
[08daea] | 1879 | /// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$ |
---|
| 1880 | CFList |
---|
[c1b9927] | 1881 | diophantine (const CFList& recResult, const CFList& factors, |
---|
[e368746] | 1882 | const CFList& products, const CFList& M, const CanonicalForm& E, |
---|
| 1883 | bool& bad) |
---|
[08daea] | 1884 | { |
---|
| 1885 | if (M.isEmpty()) |
---|
| 1886 | { |
---|
| 1887 | CFList result; |
---|
| 1888 | CFListIterator j= factors; |
---|
| 1889 | CanonicalForm buf; |
---|
| 1890 | for (CFListIterator i= recResult; i.hasItem(); i++, j++) |
---|
| 1891 | { |
---|
| 1892 | ASSERT (E.isUnivariate() || E.inCoeffDomain(), |
---|
| 1893 | "constant or univariate poly expected"); |
---|
| 1894 | ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(), |
---|
| 1895 | "constant or univariate poly expected"); |
---|
| 1896 | ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(), |
---|
| 1897 | "constant or univariate poly expected"); |
---|
| 1898 | buf= mulNTL (E, i.getItem()); |
---|
| 1899 | result.append (modNTL (buf, j.getItem())); |
---|
| 1900 | } |
---|
| 1901 | return result; |
---|
| 1902 | } |
---|
| 1903 | Variable y= M.getLast().mvar(); |
---|
| 1904 | CFList bufFactors= factors; |
---|
| 1905 | for (CFListIterator i= bufFactors; i.hasItem(); i++) |
---|
| 1906 | i.getItem()= mod (i.getItem(), y); |
---|
| 1907 | CFList bufProducts= products; |
---|
| 1908 | for (CFListIterator i= bufProducts; i.hasItem(); i++) |
---|
| 1909 | i.getItem()= mod (i.getItem(), y); |
---|
| 1910 | CFList buf= M; |
---|
| 1911 | buf.removeLast(); |
---|
| 1912 | CanonicalForm bufE= mod (E, y); |
---|
| 1913 | CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf, |
---|
[e368746] | 1914 | bufE, bad); |
---|
| 1915 | |
---|
| 1916 | if (bad) |
---|
| 1917 | return CFList(); |
---|
[08daea] | 1918 | |
---|
| 1919 | CanonicalForm e= E; |
---|
| 1920 | CFListIterator j= products; |
---|
| 1921 | for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++) |
---|
[e368746] | 1922 | e -= i.getItem()*j.getItem(); |
---|
[08daea] | 1923 | |
---|
| 1924 | CFList result= recDiophantine; |
---|
| 1925 | int d= degree (M.getLast()); |
---|
| 1926 | CanonicalForm coeffE; |
---|
| 1927 | for (int i= 1; i < d; i++) |
---|
| 1928 | { |
---|
| 1929 | if (degree (e, y) > 0) |
---|
| 1930 | coeffE= e[i]; |
---|
| 1931 | else |
---|
| 1932 | coeffE= 0; |
---|
| 1933 | if (!coeffE.isZero()) |
---|
| 1934 | { |
---|
| 1935 | CFListIterator k= result; |
---|
| 1936 | recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf, |
---|
[e368746] | 1937 | coeffE, bad); |
---|
| 1938 | if (bad) |
---|
| 1939 | return CFList(); |
---|
[08daea] | 1940 | CFListIterator l= products; |
---|
| 1941 | for (j= recDiophantine; j.hasItem(); j++, k++, l++) |
---|
| 1942 | { |
---|
| 1943 | k.getItem() += j.getItem()*power (y, i); |
---|
[e368746] | 1944 | e -= j.getItem()*power (y, i)*l.getItem(); |
---|
[08daea] | 1945 | } |
---|
| 1946 | } |
---|
| 1947 | if (e.isZero()) |
---|
| 1948 | break; |
---|
| 1949 | } |
---|
[e368746] | 1950 | if (!e.isZero()) |
---|
[08daea] | 1951 | { |
---|
[e368746] | 1952 | bad= true; |
---|
| 1953 | return CFList(); |
---|
[08daea] | 1954 | } |
---|
| 1955 | return result; |
---|
| 1956 | } |
---|
| 1957 | |
---|
| 1958 | void |
---|
[81d96c] | 1959 | nonMonicHenselStep (const CanonicalForm& F, const CFList& factors, |
---|
| 1960 | CFArray& bufFactors, const CFList& diophant, CFMatrix& M, |
---|
| 1961 | CFArray& Pi, const CFList& products, int j, |
---|
| 1962 | const CFList& MOD, bool& noOneToOne) |
---|
[08daea] | 1963 | { |
---|
| 1964 | CanonicalForm E; |
---|
| 1965 | CanonicalForm xToJ= power (F.mvar(), j); |
---|
| 1966 | Variable x= F.mvar(); |
---|
| 1967 | |
---|
| 1968 | // compute the error |
---|
| 1969 | #ifdef DEBUGOUTPUT |
---|
| 1970 | CanonicalForm test= 1; |
---|
| 1971 | for (int i= 0; i < factors.length(); i++) |
---|
| 1972 | { |
---|
| 1973 | if (i == 0) |
---|
| 1974 | test *= mod (bufFactors [i], power (x, j)); |
---|
| 1975 | else |
---|
| 1976 | test *= bufFactors[i]; |
---|
| 1977 | } |
---|
| 1978 | test= mod (test, power (x, j)); |
---|
| 1979 | test= mod (test, MOD); |
---|
| 1980 | CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1)); |
---|
[27ab36] | 1981 | DEBOUTLN (cerr, "test in nonMonicHenselStep= " << test2); |
---|
[08daea] | 1982 | #endif |
---|
| 1983 | |
---|
| 1984 | if (degree (Pi [factors.length() - 2], x) > 0) |
---|
| 1985 | E= F[j] - Pi [factors.length() - 2] [j]; |
---|
| 1986 | else |
---|
| 1987 | E= F[j]; |
---|
| 1988 | |
---|
| 1989 | CFArray buf= CFArray (diophant.length()); |
---|
| 1990 | |
---|
| 1991 | // actual lifting |
---|
[e368746] | 1992 | CFList diophantine2= diophantine (diophant, factors, products, MOD, E, |
---|
| 1993 | noOneToOne); |
---|
| 1994 | |
---|
| 1995 | if (noOneToOne) |
---|
| 1996 | return; |
---|
[08daea] | 1997 | |
---|
| 1998 | int k= 0; |
---|
| 1999 | for (CFListIterator i= diophantine2; k < factors.length(); k++, i++) |
---|
| 2000 | { |
---|
| 2001 | buf[k]= i.getItem(); |
---|
| 2002 | bufFactors[k] += xToJ*i.getItem(); |
---|
| 2003 | } |
---|
| 2004 | |
---|
| 2005 | // update Pi [0] |
---|
| 2006 | int degBuf0= degree (bufFactors[0], x); |
---|
| 2007 | int degBuf1= degree (bufFactors[1], x); |
---|
| 2008 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 2009 | { |
---|
| 2010 | M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD); |
---|
| 2011 | if (j + 2 <= M.rows()) |
---|
| 2012 | M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD); |
---|
| 2013 | } |
---|
| 2014 | CanonicalForm uIZeroJ; |
---|
| 2015 | |
---|
[1b8e048] | 2016 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 2017 | uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) + |
---|
| 2018 | mulMod (bufFactors[1] [0], buf[0], MOD); |
---|
| 2019 | else if (degBuf0 > 0) |
---|
| 2020 | uIZeroJ= mulMod (buf[0], bufFactors[1], MOD); |
---|
| 2021 | else if (degBuf1 > 0) |
---|
| 2022 | uIZeroJ= mulMod (bufFactors[0], buf[1], MOD); |
---|
| 2023 | else |
---|
| 2024 | uIZeroJ= 0; |
---|
[c1b9927] | 2025 | Pi [0] += xToJ*uIZeroJ; |
---|
[08daea] | 2026 | |
---|
| 2027 | CFArray tmp= CFArray (factors.length() - 1); |
---|
| 2028 | for (k= 0; k < factors.length() - 1; k++) |
---|
| 2029 | tmp[k]= 0; |
---|
| 2030 | CFIterator one, two; |
---|
| 2031 | one= bufFactors [0]; |
---|
| 2032 | two= bufFactors [1]; |
---|
| 2033 | if (degBuf0 > 0 && degBuf1 > 0) |
---|
| 2034 | { |
---|
[1b8e048] | 2035 | while (one.hasTerms() && one.exp() > j) one++; |
---|
| 2036 | while (two.hasTerms() && two.exp() > j) two++; |
---|
[08daea] | 2037 | for (k= 1; k <= (int) ceil (j/2.0); k++) |
---|
| 2038 | { |
---|
| 2039 | if (k != j - k + 1) |
---|
| 2040 | { |
---|
[c1b9927] | 2041 | if ((one.hasTerms() && one.exp() == j - k + 1) && |
---|
[e368746] | 2042 | (two.hasTerms() && two.exp() == j - k + 1)) |
---|
[08daea] | 2043 | { |
---|
| 2044 | tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()), |
---|
| 2045 | (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) - |
---|
| 2046 | M (j - k + 2, 1); |
---|
| 2047 | one++; |
---|
| 2048 | two++; |
---|
| 2049 | } |
---|
[e368746] | 2050 | else if (one.hasTerms() && one.exp() == j - k + 1) |
---|
[08daea] | 2051 | { |
---|
| 2052 | tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()), |
---|
| 2053 | bufFactors[1] [k], MOD) - M (k + 1, 1); |
---|
| 2054 | one++; |
---|
| 2055 | } |
---|
[e368746] | 2056 | else if (two.hasTerms() && two.exp() == j - k + 1) |
---|
[08daea] | 2057 | { |
---|
| 2058 | tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] + |
---|
| 2059 | two.coeff()), MOD) - M (k + 1, 1); |
---|
| 2060 | two++; |
---|
| 2061 | } |
---|
| 2062 | } |
---|
| 2063 | else |
---|
| 2064 | { |
---|
| 2065 | tmp[0] += M (k + 1, 1); |
---|
| 2066 | } |
---|
| 2067 | } |
---|
[1b8e048] | 2068 | } |
---|
[08daea] | 2069 | |
---|
[1b8e048] | 2070 | if (degBuf0 >= j + 1 && degBuf1 >= j + 1) |
---|
| 2071 | { |
---|
[08daea] | 2072 | if (j + 2 <= M.rows()) |
---|
[1b8e048] | 2073 | tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]), |
---|
| 2074 | (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD) |
---|
| 2075 | - M(1,1) - M (j + 2,1); |
---|
| 2076 | } |
---|
| 2077 | else if (degBuf0 >= j + 1) |
---|
| 2078 | { |
---|
| 2079 | if (degBuf1 > 0) |
---|
| 2080 | tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD); |
---|
| 2081 | else |
---|
| 2082 | tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD); |
---|
| 2083 | } |
---|
| 2084 | else if (degBuf1 >= j + 1) |
---|
| 2085 | { |
---|
| 2086 | if (degBuf0 > 0) |
---|
| 2087 | tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD); |
---|
| 2088 | else |
---|
| 2089 | tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD); |
---|
[08daea] | 2090 | } |
---|
| 2091 | Pi [0] += tmp[0]*xToJ*F.mvar(); |
---|
| 2092 | |
---|
| 2093 | // update Pi [l] |
---|
| 2094 | int degPi, degBuf; |
---|
| 2095 | for (int l= 1; l < factors.length() - 1; l++) |
---|
| 2096 | { |
---|
| 2097 | degPi= degree (Pi [l - 1], x); |
---|
| 2098 | degBuf= degree (bufFactors[l + 1], x); |
---|
| 2099 | if (degPi > 0 && degBuf > 0) |
---|
| 2100 | { |
---|
[e368746] | 2101 | M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD); |
---|
| 2102 | if (j + 2 <= M.rows()) |
---|
| 2103 | M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1], |
---|
| 2104 | MOD); |
---|
[08daea] | 2105 | } |
---|
[e368746] | 2106 | |
---|
| 2107 | if (degPi > 0 && degBuf > 0) |
---|
| 2108 | uIZeroJ= mulMod (Pi[l -1] [0], buf[l + 1], MOD) + |
---|
| 2109 | mulMod (uIZeroJ, bufFactors[l+1] [0], MOD); |
---|
| 2110 | else if (degPi > 0) |
---|
| 2111 | uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD); |
---|
| 2112 | else if (degBuf > 0) |
---|
| 2113 | uIZeroJ= mulMod (Pi[l - 1], buf[1], MOD); |
---|
[08daea] | 2114 | else |
---|
[e368746] | 2115 | uIZeroJ= 0; |
---|
| 2116 | |
---|
| 2117 | Pi [l] += xToJ*uIZeroJ; |
---|
| 2118 | |
---|
[08daea] | 2119 | one= bufFactors [l + 1]; |
---|
| 2120 | two= Pi [l - 1]; |
---|
| 2121 | if (degBuf > 0 && degPi > 0) |
---|
| 2122 | { |
---|
[e368746] | 2123 | while (one.hasTerms() && one.exp() > j) one++; |
---|
| 2124 | while (two.hasTerms() && two.exp() > j) two++; |
---|
[08daea] | 2125 | for (k= 1; k <= (int) ceil (j/2.0); k++) |
---|
| 2126 | { |
---|
| 2127 | if (k != j - k + 1) |
---|
| 2128 | { |
---|
[e368746] | 2129 | if ((one.hasTerms() && one.exp() == j - k + 1) && |
---|
| 2130 | (two.hasTerms() && two.exp() == j - k + 1)) |
---|
[08daea] | 2131 | { |
---|
| 2132 | tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()), |
---|
| 2133 | (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) - |
---|
| 2134 | M (j - k + 2, l + 1); |
---|
| 2135 | one++; |
---|
| 2136 | two++; |
---|
| 2137 | } |
---|
[e368746] | 2138 | else if (one.hasTerms() && one.exp() == j - k + 1) |
---|
[08daea] | 2139 | { |
---|
| 2140 | tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()), |
---|
| 2141 | Pi[l - 1] [k], MOD) - M (k + 1, l + 1); |
---|
| 2142 | one++; |
---|
| 2143 | } |
---|
[e368746] | 2144 | else if (two.hasTerms() && two.exp() == j - k + 1) |
---|
[08daea] | 2145 | { |
---|
[c1b9927] | 2146 | tmp[l] += mulMod (bufFactors[l + 1] [k], |
---|
[08daea] | 2147 | (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1); |
---|
| 2148 | two++; |
---|
[e368746] | 2149 | } |
---|
[08daea] | 2150 | } |
---|
| 2151 | else |
---|
| 2152 | tmp[l] += M (k + 1, l + 1); |
---|
| 2153 | } |
---|
| 2154 | } |
---|
[e368746] | 2155 | |
---|
| 2156 | if (degPi >= j + 1 && degBuf >= j + 1) |
---|
| 2157 | { |
---|
| 2158 | if (j + 2 <= M.rows()) |
---|
| 2159 | tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]), |
---|
| 2160 | (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0]) |
---|
| 2161 | , MOD) - M(1,l+1) - M (j + 2,l+1); |
---|
| 2162 | } |
---|
| 2163 | else if (degPi >= j + 1) |
---|
| 2164 | { |
---|
| 2165 | if (degBuf > 0) |
---|
| 2166 | tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD); |
---|
| 2167 | else |
---|
| 2168 | tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD); |
---|
| 2169 | } |
---|
| 2170 | else if (degBuf >= j + 1) |
---|
| 2171 | { |
---|
| 2172 | if (degPi > 0) |
---|
| 2173 | tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD); |
---|
| 2174 | else |
---|
| 2175 | tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD); |
---|
| 2176 | } |
---|
| 2177 | |
---|
[08daea] | 2178 | Pi[l] += tmp[l]*xToJ*F.mvar(); |
---|
| 2179 | } |
---|
| 2180 | return; |
---|
| 2181 | } |
---|
| 2182 | |
---|
| 2183 | // wrt. Variable (1) |
---|
| 2184 | CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c) |
---|
| 2185 | { |
---|
| 2186 | if (degree (F, 1) <= 0) |
---|
| 2187 | return c; |
---|
| 2188 | else |
---|
| 2189 | { |
---|
| 2190 | CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1)); |
---|
| 2191 | result += (swapvar (c, Variable (F.level() + 1), Variable (1)) |
---|
| 2192 | - LC (result))*power (result.mvar(), degree (result)); |
---|
| 2193 | return swapvar (result, Variable (F.level() + 1), Variable (1)); |
---|
| 2194 | } |
---|
| 2195 | } |
---|
| 2196 | |
---|
| 2197 | CFList |
---|
[81d96c] | 2198 | nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList& |
---|
| 2199 | diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1, |
---|
| 2200 | const CFList& LCs2, bool& bad) |
---|
[08daea] | 2201 | { |
---|
| 2202 | CFList buf= factors; |
---|
| 2203 | int k= 0; |
---|
| 2204 | int liftBoundBivar= l[k]; |
---|
| 2205 | CFList bufbuf= factors; |
---|
| 2206 | Variable v= Variable (2); |
---|
| 2207 | |
---|
| 2208 | CFList MOD; |
---|
| 2209 | MOD.append (power (Variable (2), liftBoundBivar)); |
---|
| 2210 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 2211 | k= 0; |
---|
| 2212 | CFListIterator j= eval; |
---|
| 2213 | j++; |
---|
| 2214 | CFListIterator iter1= LCs1; |
---|
| 2215 | CFListIterator iter2= LCs2; |
---|
| 2216 | iter1++; |
---|
| 2217 | iter2++; |
---|
| 2218 | bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem()); |
---|
| 2219 | bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem()); |
---|
| 2220 | |
---|
| 2221 | CFListIterator i= buf; |
---|
| 2222 | i++; |
---|
| 2223 | Variable y= j.getItem().mvar(); |
---|
| 2224 | if (y.level() != 3) |
---|
| 2225 | y= Variable (3); |
---|
| 2226 | |
---|
| 2227 | Pi[0]= mod (Pi[0], power (v, liftBoundBivar)); |
---|
| 2228 | M (1, 1)= Pi[0]; |
---|
| 2229 | if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0) |
---|
| 2230 | Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) + |
---|
| 2231 | mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y; |
---|
| 2232 | else if (degree (bufFactors[0], y) > 0) |
---|
| 2233 | Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y; |
---|
| 2234 | else if (degree (bufFactors[1], y) > 0) |
---|
| 2235 | Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y; |
---|
| 2236 | |
---|
| 2237 | CFList products; |
---|
| 2238 | for (int i= 0; i < bufFactors.size(); i++) |
---|
| 2239 | { |
---|
| 2240 | if (degree (bufFactors[i], y) > 0) |
---|
[e368746] | 2241 | products.append (eval.getFirst()/bufFactors[i] [0]); |
---|
[08daea] | 2242 | else |
---|
[e368746] | 2243 | products.append (eval.getFirst()/bufFactors[i]); |
---|
[08daea] | 2244 | } |
---|
| 2245 | |
---|
| 2246 | for (int d= 1; d < l[1]; d++) |
---|
[e368746] | 2247 | { |
---|
[81d96c] | 2248 | nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products, |
---|
| 2249 | d, MOD, bad); |
---|
[e368746] | 2250 | if (bad) |
---|
| 2251 | return CFList(); |
---|
| 2252 | } |
---|
[08daea] | 2253 | CFList result; |
---|
| 2254 | for (k= 0; k < factors.length(); k++) |
---|
| 2255 | result.append (bufFactors[k]); |
---|
| 2256 | return result; |
---|
| 2257 | } |
---|
| 2258 | |
---|
| 2259 | |
---|
| 2260 | CFList |
---|
[81d96c] | 2261 | nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, |
---|
| 2262 | CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld, |
---|
| 2263 | int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad |
---|
| 2264 | ) |
---|
[08daea] | 2265 | { |
---|
| 2266 | int k= 0; |
---|
| 2267 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 2268 | bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast()); |
---|
| 2269 | bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast()); |
---|
| 2270 | CFList buf= factors; |
---|
| 2271 | Variable y= F.getLast().mvar(); |
---|
| 2272 | Variable x= F.getFirst().mvar(); |
---|
| 2273 | CanonicalForm xToLOld= power (x, lOld); |
---|
| 2274 | Pi [0]= mod (Pi[0], xToLOld); |
---|
| 2275 | M (1, 1)= Pi [0]; |
---|
| 2276 | |
---|
| 2277 | if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0) |
---|
[c1b9927] | 2278 | Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) + |
---|
[08daea] | 2279 | mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y; |
---|
| 2280 | else if (degree (bufFactors[0], y) > 0) |
---|
| 2281 | Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y; |
---|
| 2282 | else if (degree (bufFactors[1], y) > 0) |
---|
| 2283 | Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y; |
---|
| 2284 | |
---|
| 2285 | CFList products; |
---|
[21b8f4c] | 2286 | CanonicalForm quot; |
---|
[08daea] | 2287 | for (int i= 0; i < bufFactors.size(); i++) |
---|
| 2288 | { |
---|
| 2289 | if (degree (bufFactors[i], y) > 0) |
---|
| 2290 | { |
---|
[21b8f4c] | 2291 | if (!fdivides (bufFactors[i] [0], F.getFirst(), quot)) |
---|
[e368746] | 2292 | { |
---|
| 2293 | bad= true; |
---|
| 2294 | return CFList(); |
---|
| 2295 | } |
---|
[21b8f4c] | 2296 | products.append (quot); |
---|
[08daea] | 2297 | } |
---|
| 2298 | else |
---|
| 2299 | { |
---|
[21b8f4c] | 2300 | if (!fdivides (bufFactors[i], F.getFirst(), quot)) |
---|
[e368746] | 2301 | { |
---|
| 2302 | bad= true; |
---|
| 2303 | return CFList(); |
---|
| 2304 | } |
---|
[21b8f4c] | 2305 | products.append (quot); |
---|
[08daea] | 2306 | } |
---|
| 2307 | } |
---|
| 2308 | |
---|
| 2309 | for (int d= 1; d < lNew; d++) |
---|
[e368746] | 2310 | { |
---|
[81d96c] | 2311 | nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products, |
---|
| 2312 | d, MOD, bad); |
---|
[e368746] | 2313 | if (bad) |
---|
| 2314 | return CFList(); |
---|
| 2315 | } |
---|
[08daea] | 2316 | |
---|
| 2317 | CFList result; |
---|
| 2318 | for (k= 0; k < factors.length(); k++) |
---|
| 2319 | result.append (bufFactors[k]); |
---|
| 2320 | return result; |
---|
| 2321 | } |
---|
| 2322 | |
---|
| 2323 | CFList |
---|
[81d96c] | 2324 | nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int |
---|
| 2325 | lLength, bool sort, const CFList& LCs1, const CFList& LCs2, |
---|
| 2326 | const CFArray& Pi, const CFList& diophant, bool& bad) |
---|
[08daea] | 2327 | { |
---|
| 2328 | CFList bufDiophant= diophant; |
---|
| 2329 | CFList buf= factors; |
---|
| 2330 | if (sort) |
---|
| 2331 | sortList (buf, Variable (1)); |
---|
| 2332 | CFArray bufPi= Pi; |
---|
| 2333 | CFMatrix M= CFMatrix (l[1], factors.length()); |
---|
[81d96c] | 2334 | CFList result= |
---|
| 2335 | nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad); |
---|
[e368746] | 2336 | if (bad) |
---|
| 2337 | return CFList(); |
---|
| 2338 | |
---|
[08daea] | 2339 | if (eval.length() == 2) |
---|
| 2340 | return result; |
---|
| 2341 | CFList MOD; |
---|
| 2342 | for (int i= 0; i < 2; i++) |
---|
| 2343 | MOD.append (power (Variable (i + 2), l[i])); |
---|
| 2344 | CFListIterator j= eval; |
---|
| 2345 | j++; |
---|
| 2346 | CFList bufEval; |
---|
| 2347 | bufEval.append (j.getItem()); |
---|
| 2348 | j++; |
---|
| 2349 | CFListIterator jj= LCs1; |
---|
| 2350 | CFListIterator jjj= LCs2; |
---|
| 2351 | CFList bufLCs1, bufLCs2; |
---|
| 2352 | jj++, jjj++; |
---|
| 2353 | bufLCs1.append (jj.getItem()); |
---|
| 2354 | bufLCs2.append (jjj.getItem()); |
---|
| 2355 | jj++, jjj++; |
---|
| 2356 | |
---|
[ea88e0] | 2357 | for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++) |
---|
[08daea] | 2358 | { |
---|
| 2359 | bufEval.append (j.getItem()); |
---|
| 2360 | bufLCs1.append (jj.getItem()); |
---|
| 2361 | bufLCs2.append (jjj.getItem()); |
---|
| 2362 | M= CFMatrix (l[i], factors.length()); |
---|
[81d96c] | 2363 | result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, |
---|
| 2364 | l[i - 1], l[i], bufLCs1, bufLCs2, bad); |
---|
[e368746] | 2365 | if (bad) |
---|
| 2366 | return CFList(); |
---|
[08daea] | 2367 | MOD.append (power (Variable (i + 2), l[i])); |
---|
| 2368 | bufEval.removeFirst(); |
---|
| 2369 | bufLCs1.removeFirst(); |
---|
| 2370 | bufLCs2.removeFirst(); |
---|
| 2371 | } |
---|
| 2372 | return result; |
---|
| 2373 | } |
---|
| 2374 | |
---|
[e368746] | 2375 | CFList |
---|
[c1b9927] | 2376 | nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const |
---|
[e368746] | 2377 | CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound, |
---|
| 2378 | int bivarLiftBound, bool& bad) |
---|
| 2379 | { |
---|
| 2380 | CFList bufFactors2= factors; |
---|
| 2381 | |
---|
| 2382 | Variable y= Variable (2); |
---|
| 2383 | for (CFListIterator i= bufFactors2; i.hasItem(); i++) |
---|
| 2384 | i.getItem()= mod (i.getItem(), y); |
---|
| 2385 | |
---|
| 2386 | CanonicalForm bufF= F; |
---|
| 2387 | bufF= mod (bufF, y); |
---|
| 2388 | bufF= mod (bufF, Variable (3)); |
---|
| 2389 | |
---|
| 2390 | diophant= diophantine (bufF, bufFactors2); |
---|
| 2391 | |
---|
| 2392 | CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1); |
---|
| 2393 | |
---|
| 2394 | Pi= CFArray (bufFactors2.length() - 1); |
---|
| 2395 | |
---|
| 2396 | CFArray bufFactors= CFArray (bufFactors2.length()); |
---|
| 2397 | CFListIterator j= LCs; |
---|
| 2398 | int i= 0; |
---|
| 2399 | for (CFListIterator k= factors; k.hasItem(); j++, k++, i++) |
---|
| 2400 | bufFactors[i]= replaceLC (k.getItem(), j.getItem()); |
---|
| 2401 | |
---|
| 2402 | //initialise Pi |
---|
| 2403 | Variable v= Variable (3); |
---|
| 2404 | CanonicalForm yToL= power (y, bivarLiftBound); |
---|
| 2405 | if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0) |
---|
| 2406 | { |
---|
| 2407 | M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL); |
---|
[c1b9927] | 2408 | Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) + |
---|
[e368746] | 2409 | mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v; |
---|
| 2410 | } |
---|
| 2411 | else if (degree (bufFactors[0], v) > 0) |
---|
| 2412 | { |
---|
| 2413 | M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL); |
---|
| 2414 | Pi [0]= M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v; |
---|
| 2415 | } |
---|
| 2416 | else if (degree (bufFactors[1], v) > 0) |
---|
| 2417 | { |
---|
| 2418 | M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL); |
---|
| 2419 | Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v; |
---|
| 2420 | } |
---|
| 2421 | else |
---|
| 2422 | { |
---|
| 2423 | M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL); |
---|
| 2424 | Pi [0]= M (1,1); |
---|
| 2425 | } |
---|
| 2426 | |
---|
| 2427 | for (i= 1; i < Pi.size(); i++) |
---|
| 2428 | { |
---|
| 2429 | if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0) |
---|
| 2430 | { |
---|
| 2431 | M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL); |
---|
[c1b9927] | 2432 | Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) + |
---|
[e368746] | 2433 | mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v; |
---|
| 2434 | } |
---|
| 2435 | else if (degree (Pi[i-1], v) > 0) |
---|
| 2436 | { |
---|
| 2437 | M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL); |
---|
| 2438 | Pi [i]= M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v; |
---|
| 2439 | } |
---|
| 2440 | else if (degree (bufFactors[i+1], v) > 0) |
---|
| 2441 | { |
---|
| 2442 | M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL); |
---|
| 2443 | Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v; |
---|
| 2444 | } |
---|
| 2445 | else |
---|
| 2446 | { |
---|
| 2447 | M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL); |
---|
| 2448 | Pi [i]= M (1,i+1); |
---|
| 2449 | } |
---|
| 2450 | } |
---|
| 2451 | |
---|
| 2452 | CFList products; |
---|
| 2453 | bufF= mod (F, Variable (3)); |
---|
| 2454 | for (CFListIterator k= factors; k.hasItem(); k++) |
---|
| 2455 | products.append (bufF/k.getItem()); |
---|
| 2456 | |
---|
| 2457 | CFList MOD= CFList (power (v, liftBound)); |
---|
| 2458 | MOD.insert (yToL); |
---|
| 2459 | for (int d= 1; d < liftBound; d++) |
---|
| 2460 | { |
---|
[81d96c] | 2461 | nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d, |
---|
| 2462 | MOD, bad); |
---|
[e368746] | 2463 | if (bad) |
---|
| 2464 | return CFList(); |
---|
| 2465 | } |
---|
| 2466 | |
---|
| 2467 | CFList result; |
---|
| 2468 | for (i= 0; i < factors.length(); i++) |
---|
| 2469 | result.append (bufFactors[i]); |
---|
| 2470 | return result; |
---|
| 2471 | } |
---|
| 2472 | |
---|
| 2473 | CFList |
---|
| 2474 | nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs, |
---|
[81d96c] | 2475 | CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld, |
---|
[e368746] | 2476 | int& lNew, const CFList& MOD, bool& noOneToOne |
---|
| 2477 | ) |
---|
| 2478 | { |
---|
| 2479 | |
---|
| 2480 | int k= 0; |
---|
| 2481 | CFArray bufFactors= CFArray (factors.length()); |
---|
| 2482 | CFListIterator j= LCs; |
---|
| 2483 | for (CFListIterator i= factors; i.hasItem(); i++, j++, k++) |
---|
| 2484 | bufFactors [k]= replaceLC (i.getItem(), j.getItem()); |
---|
| 2485 | |
---|
| 2486 | Variable y= F.getLast().mvar(); |
---|
| 2487 | Variable x= F.getFirst().mvar(); |
---|
| 2488 | CanonicalForm xToLOld= power (x, lOld); |
---|
| 2489 | |
---|
| 2490 | Pi [0]= mod (Pi[0], xToLOld); |
---|
| 2491 | M (1, 1)= Pi [0]; |
---|
| 2492 | |
---|
| 2493 | if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0) |
---|
[c1b9927] | 2494 | Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) + |
---|
[e368746] | 2495 | mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y; |
---|
| 2496 | else if (degree (bufFactors[0], y) > 0) |
---|
| 2497 | Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y; |
---|
| 2498 | else if (degree (bufFactors[1], y) > 0) |
---|
| 2499 | Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y; |
---|
| 2500 | |
---|
| 2501 | for (int i= 1; i < Pi.size(); i++) |
---|
| 2502 | { |
---|
| 2503 | Pi [i]= mod (Pi [i], xToLOld); |
---|
| 2504 | M (1, i + 1)= Pi [i]; |
---|
| 2505 | |
---|
| 2506 | if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0) |
---|
| 2507 | Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) + |
---|
| 2508 | mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y; |
---|
| 2509 | else if (degree (Pi[i-1], y) > 0) |
---|
| 2510 | Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y; |
---|
| 2511 | else if (degree (bufFactors[i+1], y) > 0) |
---|
| 2512 | Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y; |
---|
| 2513 | } |
---|
| 2514 | |
---|
| 2515 | CFList products; |
---|
[21b8f4c] | 2516 | CanonicalForm quot, bufF= F.getFirst(); |
---|
[e368746] | 2517 | |
---|
| 2518 | for (int i= 0; i < bufFactors.size(); i++) |
---|
| 2519 | { |
---|
| 2520 | if (degree (bufFactors[i], y) > 0) |
---|
| 2521 | { |
---|
[21b8f4c] | 2522 | if (!fdivides (bufFactors[i] [0], bufF, quot)) |
---|
[e368746] | 2523 | { |
---|
| 2524 | noOneToOne= true; |
---|
| 2525 | return factors; |
---|
| 2526 | } |
---|
[21b8f4c] | 2527 | products.append (quot); |
---|
[e368746] | 2528 | } |
---|
| 2529 | else |
---|
| 2530 | { |
---|
[21b8f4c] | 2531 | if (!fdivides (bufFactors[i], bufF, quot)) |
---|
[e368746] | 2532 | { |
---|
| 2533 | noOneToOne= true; |
---|
| 2534 | return factors; |
---|
| 2535 | } |
---|
[21b8f4c] | 2536 | products.append (quot); |
---|
[e368746] | 2537 | } |
---|
| 2538 | } |
---|
| 2539 | |
---|
| 2540 | for (int d= 1; d < lNew; d++) |
---|
| 2541 | { |
---|
[81d96c] | 2542 | nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi, |
---|
| 2543 | products, d, MOD, noOneToOne); |
---|
[e368746] | 2544 | if (noOneToOne) |
---|
| 2545 | return CFList(); |
---|
| 2546 | } |
---|
| 2547 | |
---|
| 2548 | CFList result; |
---|
| 2549 | for (k= 0; k < factors.length(); k++) |
---|
| 2550 | result.append (bufFactors[k]); |
---|
| 2551 | return result; |
---|
| 2552 | } |
---|
| 2553 | |
---|
| 2554 | CFList |
---|
| 2555 | nonMonicHenselLift (const CFList& eval, const CFList& factors, |
---|
| 2556 | CFList* const& LCs, CFList& diophant, CFArray& Pi, |
---|
| 2557 | int* liftBound, int length, bool& noOneToOne |
---|
| 2558 | ) |
---|
| 2559 | { |
---|
| 2560 | CFList bufDiophant= diophant; |
---|
| 2561 | CFList buf= factors; |
---|
| 2562 | CFArray bufPi= Pi; |
---|
| 2563 | CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1); |
---|
| 2564 | int k= 0; |
---|
| 2565 | |
---|
| 2566 | CFList result= |
---|
| 2567 | nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi, |
---|
| 2568 | liftBound[1], liftBound[0], noOneToOne); |
---|
| 2569 | |
---|
| 2570 | if (noOneToOne) |
---|
| 2571 | return CFList(); |
---|
| 2572 | |
---|
| 2573 | if (eval.length() == 1) |
---|
| 2574 | return result; |
---|
| 2575 | |
---|
| 2576 | k++; |
---|
| 2577 | CFList MOD; |
---|
| 2578 | for (int i= 0; i < 2; i++) |
---|
| 2579 | MOD.append (power (Variable (i + 2), liftBound[i])); |
---|
| 2580 | |
---|
| 2581 | CFListIterator j= eval; |
---|
| 2582 | CFList bufEval; |
---|
| 2583 | bufEval.append (j.getItem()); |
---|
| 2584 | j++; |
---|
| 2585 | |
---|
| 2586 | for (int i= 2; i <= length && j.hasItem(); i++, j++, k++) |
---|
| 2587 | { |
---|
| 2588 | bufEval.append (j.getItem()); |
---|
| 2589 | M= CFMatrix (liftBound[i], factors.length() - 1); |
---|
| 2590 | result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M, |
---|
| 2591 | liftBound[i-1], liftBound[i], MOD, noOneToOne); |
---|
| 2592 | if (noOneToOne) |
---|
| 2593 | return result; |
---|
| 2594 | MOD.append (power (Variable (i + 2), liftBound[i])); |
---|
| 2595 | bufEval.removeFirst(); |
---|
| 2596 | } |
---|
| 2597 | |
---|
| 2598 | return result; |
---|
| 2599 | } |
---|
| 2600 | |
---|
[ad3c3ff] | 2601 | #endif |
---|
| 2602 | /* HAVE_NTL */ |
---|
| 2603 | |
---|