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