[dea3d2] | 1 | /*****************************************************************************\ |
---|
| 2 | * Computer Algebra System SINGULAR |
---|
| 3 | \*****************************************************************************/ |
---|
[058c1d] | 4 | /** @file facAbsFact.cc |
---|
[dea3d2] | 5 | * |
---|
| 6 | * @author Martin Lee |
---|
| 7 | * |
---|
| 8 | **/ |
---|
| 9 | /*****************************************************************************/ |
---|
| 10 | |
---|
[17a710] | 11 | #include "timing.h" |
---|
| 12 | #include "debug.h" |
---|
| 13 | |
---|
[5275c0] | 14 | #include "facAbsFact.h" |
---|
[17a710] | 15 | #include "facBivar.h" |
---|
| 16 | #include "facFqBivar.h" |
---|
[dea3d2] | 17 | #include "cf_reval.h" |
---|
| 18 | #include "cf_primes.h" |
---|
| 19 | #include "cf_algorithm.h" |
---|
| 20 | #ifdef HAVE_FLINT |
---|
| 21 | #include "FLINTconvert.h" |
---|
| 22 | #include <flint/fmpz_poly_factor.h> |
---|
| 23 | #endif |
---|
| 24 | #ifdef HAVE_NTL |
---|
| 25 | #include "NTLconvert.h" |
---|
| 26 | #include <NTL/LLL.h> |
---|
| 27 | #endif |
---|
| 28 | |
---|
| 29 | #ifdef HAVE_FLINT |
---|
[160ec6] | 30 | #ifdef HAVE_NTL |
---|
[dea3d2] | 31 | |
---|
[17a710] | 32 | TIMING_DEFINE_PRINT(fac_Qa_factorize) |
---|
| 33 | TIMING_DEFINE_PRINT(fac_evalpoint) |
---|
| 34 | |
---|
| 35 | void out_cf(const char *s1,const CanonicalForm &f,const char *s2); |
---|
| 36 | |
---|
[160ec6] | 37 | //TODO optimize choice of p -> choose p as large as possible (better than small p since factorization mod p does not require field extension, also less lifting) |
---|
[17a710] | 38 | int choosePoint (const CanonicalForm& F, int tdegF, CFArray& eval, bool rec) |
---|
[dea3d2] | 39 | { |
---|
| 40 | REvaluation E1 (1, 1, IntRandom (25)); |
---|
| 41 | REvaluation E2 (2, 2, IntRandom (25)); |
---|
[17a710] | 42 | if (rec) |
---|
| 43 | { |
---|
| 44 | E1.nextpoint(); |
---|
| 45 | E2.nextpoint(); |
---|
| 46 | } |
---|
| 47 | CanonicalForm f, f1, f2, Fp; |
---|
| 48 | int i, p; |
---|
[dea3d2] | 49 | eval=CFArray (2); |
---|
[17a710] | 50 | printf ("tdegF= %d\n", tdegF); |
---|
| 51 | out_cf ("F= ", F, "\n"); |
---|
| 52 | printf ("getCharacteristic()= %d\n", getCharacteristic()); |
---|
[dea3d2] | 53 | while (1) |
---|
| 54 | { |
---|
[17a710] | 55 | f1= E1(F); |
---|
| 56 | if (!f1.isZero() && factorize (f1).length() == 2) |
---|
[dea3d2] | 57 | { |
---|
| 58 | Off (SW_RATIONAL); |
---|
[17a710] | 59 | f= E2(f1); |
---|
| 60 | f2= E2 (F); |
---|
| 61 | out_cf ("f= ", f, "\n"); |
---|
| 62 | printf ("isOn (SW_RATIONAL)= %d\n", isOn (SW_RATIONAL)); |
---|
| 63 | printf ("cf_getSmallPrime (cf_getNumSmallPrimes())= %d\n", cf_getSmallPrime (cf_getNumSmallPrimes()-1)); |
---|
| 64 | if ((!f.isZero()) && (abs(f) > cf_getSmallPrime (cf_getNumSmallPrimes()-1))) |
---|
[dea3d2] | 65 | { |
---|
[17a710] | 66 | printf ("hier0\n"); |
---|
| 67 | for (i= cf_getNumPrimes()-1; i >= 0; i--) |
---|
[dea3d2] | 68 | { |
---|
| 69 | if (f % CanonicalForm (cf_getPrime (i)) == 0) |
---|
| 70 | { |
---|
[17a710] | 71 | p= cf_getPrime(i); |
---|
| 72 | Fp= mod (F,p); |
---|
| 73 | out_cf ("Fp0= ", Fp, "\n"); |
---|
| 74 | if (totaldegree (Fp) == tdegF && degree (mod (f2,p), 1) == degree (F,1) && degree (mod (f1, p),2) == degree (F,2)) |
---|
[dea3d2] | 75 | { |
---|
| 76 | eval[0]= E1[1]; |
---|
| 77 | eval[1]= E2[2]; |
---|
[17a710] | 78 | return p; |
---|
[dea3d2] | 79 | } |
---|
| 80 | } |
---|
| 81 | } |
---|
| 82 | } |
---|
| 83 | else if (!f.isZero()) |
---|
| 84 | { |
---|
[17a710] | 85 | printf ("hier\n"); |
---|
| 86 | for (i= cf_getNumSmallPrimes()-1; i >= 0; i--) |
---|
[dea3d2] | 87 | { |
---|
| 88 | if (f % CanonicalForm (cf_getSmallPrime (i)) == 0) |
---|
| 89 | { |
---|
[17a710] | 90 | p= cf_getSmallPrime (i); |
---|
| 91 | Fp= mod (F,p); |
---|
| 92 | out_cf ("Fp1= ", Fp, "\n"); |
---|
| 93 | if (totaldegree (Fp) == tdegF && degree (mod (f2, p),1) == degree (F,1) && degree (mod (f1,p),2) == degree (F,2)) |
---|
[dea3d2] | 94 | { |
---|
| 95 | eval[0]= E1[1]; |
---|
| 96 | eval[1]= E2[2]; |
---|
[17a710] | 97 | return p; |
---|
[dea3d2] | 98 | } |
---|
| 99 | } |
---|
| 100 | } |
---|
| 101 | } |
---|
| 102 | E2.nextpoint(); |
---|
| 103 | On (SW_RATIONAL); |
---|
[17a710] | 104 | out_cf ("E2= ", E2[2], "\n"); |
---|
[dea3d2] | 105 | } |
---|
| 106 | E1.nextpoint(); |
---|
[17a710] | 107 | out_cf ("E1= ", E1[1], "\n"); |
---|
[dea3d2] | 108 | } |
---|
| 109 | return 0; |
---|
| 110 | } |
---|
| 111 | |
---|
[17a710] | 112 | //TODO sowohl bzgl. x als auch y auswerten und minpoly berechnen |
---|
[efd410] | 113 | CFAFList absFactorizeMain (const CanonicalForm& G) |
---|
[dea3d2] | 114 | { |
---|
| 115 | //F is assumed to be bivariate, irreducible over Q, primitive wrt x and y, compressed |
---|
| 116 | |
---|
[17a710] | 117 | out_cf ("F= ", G, "\n"); |
---|
[dea3d2] | 118 | CanonicalForm F= bCommonDen (G)*G; |
---|
| 119 | Off (SW_RATIONAL); |
---|
| 120 | F /= icontent (F); |
---|
[17a710] | 121 | out_cf ("F after icontent= ", F, "\n"); |
---|
[dea3d2] | 122 | On (SW_RATIONAL); |
---|
| 123 | CFArray eval; |
---|
| 124 | int minTdeg, tdegF= totaldegree (F); |
---|
| 125 | CanonicalForm Fp, smallestFactor; |
---|
| 126 | int p; |
---|
[17a710] | 127 | CFFList factors; |
---|
| 128 | Variable alpha; |
---|
| 129 | bool rec= false; |
---|
| 130 | Variable x= Variable (1); |
---|
| 131 | Variable y= Variable (2); |
---|
| 132 | differentevalpoint: |
---|
[dea3d2] | 133 | while (1) |
---|
| 134 | { |
---|
[17a710] | 135 | TIMING_START (fac_evalpoint); |
---|
| 136 | p= choosePoint (F, tdegF, eval, rec); |
---|
| 137 | TIMING_END_AND_PRINT (fac_evalpoint, "time to find eval point: "); |
---|
[dea3d2] | 138 | |
---|
| 139 | setCharacteristic (p); |
---|
| 140 | Fp=F.mapinto(); |
---|
[17a710] | 141 | factors= factorize (Fp); |
---|
| 142 | |
---|
| 143 | for (CFFListIterator iter= factors; iter.hasItem(); iter++) |
---|
| 144 | { |
---|
| 145 | out_cf ("factors= ", iter.getItem().factor(), "\n"); |
---|
| 146 | printf ("exp= %d\n", iter.getItem().exp()); |
---|
| 147 | } |
---|
| 148 | printf ("p= %d\n", p); |
---|
| 149 | if (factors.getFirst().factor().inCoeffDomain()) |
---|
| 150 | factors.removeFirst(); |
---|
| 151 | printf ("factors.length()= %d\n", factors.length()); |
---|
| 152 | printf ("factors.getFirst().exp()= %d\n", factors.getFirst().exp()); |
---|
| 153 | if (factors.length() == 1 && factors.getFirst().exp() == 1) |
---|
| 154 | { |
---|
| 155 | if (absIrredTest (Fp)) //TODO absIrredTest mit shift, modular absIrredTest |
---|
| 156 | { |
---|
| 157 | printf ("irred after test\n"); |
---|
| 158 | printf ("absIrred\n"); |
---|
| 159 | setCharacteristic(0); |
---|
| 160 | alpha= rootOf (x); |
---|
| 161 | out_cf ("G= ", G, "\n"); |
---|
| 162 | out_cf ("getMipo (alpha)= ", getMipo (alpha), "\n"); |
---|
[dea3d2] | 163 | |
---|
[17a710] | 164 | return CFAFList (CFAFactor (G, getMipo (alpha), 1)); |
---|
| 165 | } |
---|
| 166 | else |
---|
| 167 | { |
---|
| 168 | setCharacteristic (0); |
---|
| 169 | if (modularIrredTestWithShift (F)) |
---|
| 170 | { |
---|
| 171 | printf ("irred after modular test\n"); |
---|
| 172 | alpha= rootOf (x); |
---|
| 173 | return CFAFList (CFAFactor (G, getMipo (alpha), 1)); |
---|
| 174 | } |
---|
| 175 | rec= true; |
---|
| 176 | continue; |
---|
| 177 | } |
---|
| 178 | } |
---|
[dea3d2] | 179 | CFFListIterator iter= factors; |
---|
| 180 | smallestFactor= iter.getItem().factor(); |
---|
[17a710] | 181 | out_cf ("smallestFactor before= ", smallestFactor, "\n"); |
---|
| 182 | while (smallestFactor.isUnivariate() && iter.hasItem()) |
---|
| 183 | { |
---|
| 184 | iter++; |
---|
| 185 | if (!iter.hasItem()) |
---|
| 186 | break; |
---|
| 187 | out_cf ("factors= ", iter.getItem().factor(), "\n"); |
---|
| 188 | printf ("exp= %d\n", iter.getItem().exp()); |
---|
| 189 | smallestFactor= iter.getItem().factor(); |
---|
| 190 | } |
---|
| 191 | //TODO univariate Faktoren rausschmeiÃen! |
---|
[dea3d2] | 192 | minTdeg= totaldegree (smallestFactor); |
---|
[17a710] | 193 | if (iter.hasItem()) |
---|
| 194 | iter++; |
---|
[dea3d2] | 195 | for (; iter.hasItem(); iter++) |
---|
| 196 | { |
---|
[17a710] | 197 | out_cf ("factors= ", iter.getItem().factor(), "\n"); |
---|
| 198 | printf ("exp= %d\n", iter.getItem().exp()); |
---|
| 199 | if (!iter.getItem().factor().isUnivariate() && (totaldegree (iter.getItem().factor()) < minTdeg)) |
---|
[dea3d2] | 200 | { |
---|
| 201 | smallestFactor= iter.getItem().factor(); |
---|
| 202 | minTdeg= totaldegree (smallestFactor); |
---|
| 203 | } |
---|
| 204 | } |
---|
| 205 | if (tdegF % minTdeg == 0) |
---|
| 206 | break; |
---|
[17a710] | 207 | setCharacteristic(0); |
---|
| 208 | rec=true; |
---|
[dea3d2] | 209 | } |
---|
| 210 | CanonicalForm Gp= Fp/smallestFactor; |
---|
[17a710] | 211 | out_cf ("Gp before= ", Gp, "\n"); |
---|
| 212 | out_cf ("smallestFactor= ", smallestFactor, "\n"); |
---|
| 213 | printf ("degree (Gp,1)= %d\n", degree (Gp, 1)); |
---|
| 214 | printf ("degree smallestFactor= %d\n", degree (smallestFactor, 1)); |
---|
| 215 | printf ("degree (Fp,1)= %d\n", degree (Fp,1)); |
---|
| 216 | printf ("degree (F,1)= %d\n", degree (F,1)); |
---|
| 217 | out_cf ("eval[1]= ", eval[1], "\n"); |
---|
| 218 | out_cf ("eval[0]= ", eval[0], "\n"); |
---|
| 219 | //printf ("Gp*smallestFactor==Fp ? %d\n", Gp*smallestFactor == Fp); |
---|
| 220 | Gp= Gp /Lc (Gp); |
---|
| 221 | |
---|
| 222 | CanonicalForm Gpy= Gp (eval[0].mapinto(), 1); |
---|
| 223 | CanonicalForm smallestFactorEvaly= smallestFactor (eval[0].mapinto(),1); |
---|
| 224 | CanonicalForm Gpx= Gp (eval[1].mapinto(), 2); |
---|
| 225 | out_cf ("Gp eval= ", Gp, "\n"); |
---|
| 226 | CanonicalForm smallestFactorEvalx= smallestFactor (eval[1].mapinto(),2); |
---|
| 227 | |
---|
| 228 | out_cf ("smallestFactorEvalx= ", smallestFactorEvalx, "\n"); |
---|
| 229 | out_cf ("gcd (Gpx, smallestFactorEvalx)= ", gcd (Gpx, smallestFactorEvalx), "\n"); |
---|
| 230 | bool xValid= !(Gpx.inCoeffDomain() || smallestFactorEvalx.inCoeffDomain() || !gcd (Gpx, smallestFactorEvalx).inCoeffDomain()); |
---|
| 231 | bool yValid= !(Gpy.inCoeffDomain() || smallestFactorEvaly.inCoeffDomain() || !gcd (Gpy, smallestFactorEvaly).inCoeffDomain()); |
---|
| 232 | if (!xValid && !yValid) |
---|
| 233 | { |
---|
| 234 | rec= true; |
---|
| 235 | setCharacteristic (0); |
---|
| 236 | printf ("goto1\n"); |
---|
| 237 | goto differentevalpoint; |
---|
| 238 | } |
---|
| 239 | |
---|
[dea3d2] | 240 | setCharacteristic (0); |
---|
| 241 | |
---|
[17a710] | 242 | CanonicalForm mipo; |
---|
| 243 | |
---|
| 244 | int loop, i; |
---|
| 245 | if (xValid && yValid) |
---|
| 246 | { |
---|
| 247 | loop= 3; |
---|
| 248 | i=1; |
---|
| 249 | } |
---|
| 250 | else if (xValid) |
---|
| 251 | { |
---|
| 252 | loop= 3; |
---|
| 253 | i=2; |
---|
| 254 | } |
---|
| 255 | else |
---|
| 256 | { |
---|
| 257 | loop= 2; |
---|
| 258 | i=1; |
---|
| 259 | } |
---|
| 260 | |
---|
| 261 | CFArray mipos= CFArray (loop-i); |
---|
| 262 | printf ("loop= %d\n", loop); |
---|
| 263 | printf ("xValid= %d\n", xValid); |
---|
| 264 | printf ("yValid= %d\n", yValid); |
---|
| 265 | |
---|
| 266 | for (; i < loop; i++) |
---|
| 267 | { |
---|
| 268 | CanonicalForm Fi= F(eval[i-1],i); |
---|
| 269 | //CanonicalForm Fx= F(eval[0],1); |
---|
| 270 | //CanonicalForm Fy= F(eval[1],2); |
---|
| 271 | |
---|
| 272 | int s= tdegF/minTdeg + 1; |
---|
| 273 | CanonicalForm bound= power (maxNorm (Fi), 2*(s-1)); |
---|
| 274 | bound *= power (CanonicalForm (s),s-1); |
---|
| 275 | bound *= power (CanonicalForm (2), ((s-1)*(s-1))/2); //possible int overflow |
---|
| 276 | |
---|
| 277 | CanonicalForm B = p; |
---|
| 278 | long k = 1L; |
---|
| 279 | while ( B < bound ) { |
---|
| 280 | B *= p; |
---|
| 281 | k++; |
---|
| 282 | } |
---|
| 283 | |
---|
| 284 | //TODO take floor (log2(k)) |
---|
| 285 | k= k+1; |
---|
| 286 | fmpz_poly_t FLINTFi; |
---|
| 287 | out_cf ("Fi= ", Fi, "\n"); |
---|
| 288 | convertFacCF2Fmpz_poly_t (FLINTFi, Fi); |
---|
| 289 | setCharacteristic (p); |
---|
| 290 | printf ("p= %d\n", p); |
---|
| 291 | nmod_poly_t FLINTFpi, FLINTGpi; |
---|
| 292 | if (i == 2) |
---|
| 293 | { |
---|
| 294 | convertFacCF2nmod_poly_t (FLINTFpi, smallestFactorEvalx/lc (smallestFactorEvalx)); |
---|
| 295 | convertFacCF2nmod_poly_t (FLINTGpi, Gpx/lc (Gpx)); |
---|
| 296 | } |
---|
| 297 | else |
---|
| 298 | { |
---|
| 299 | convertFacCF2nmod_poly_t (FLINTFpi, smallestFactorEvaly/lc (smallestFactorEvaly)); |
---|
| 300 | convertFacCF2nmod_poly_t (FLINTGpi, Gpy/lc (Gpy)); |
---|
| 301 | } |
---|
| 302 | nmod_poly_factor_t nmodFactors; |
---|
| 303 | nmod_poly_factor_init (nmodFactors); |
---|
| 304 | nmod_poly_factor_insert (nmodFactors, FLINTFpi, 1L); |
---|
| 305 | nmod_poly_factor_insert (nmodFactors, FLINTGpi, 1L); |
---|
| 306 | |
---|
| 307 | //out_cf ("Gpx= ", Gpx, "\n"); |
---|
| 308 | //out_cf ("smallestFactorEvalx= ", smallestFactorEvalx, "\n"); |
---|
| 309 | long * link= new long [2]; |
---|
| 310 | fmpz_poly_t *v= new fmpz_poly_t[2]; |
---|
| 311 | fmpz_poly_t *w= new fmpz_poly_t[2]; |
---|
| 312 | fmpz_poly_init(v[0]); |
---|
| 313 | fmpz_poly_init(v[1]); |
---|
| 314 | fmpz_poly_init(w[0]); |
---|
| 315 | fmpz_poly_init(w[1]); |
---|
| 316 | |
---|
| 317 | printf ("k= %ld\n", k); |
---|
| 318 | fmpz_poly_factor_t liftedFactors; |
---|
| 319 | fmpz_poly_factor_init (liftedFactors); |
---|
| 320 | _fmpz_poly_hensel_start_lift(liftedFactors, link, v, w, FLINTFi, nmodFactors, k); //lift factors up to p^k |
---|
| 321 | |
---|
| 322 | nmod_poly_factor_clear (nmodFactors); |
---|
| 323 | nmod_poly_clear (FLINTFpi); |
---|
| 324 | nmod_poly_clear (FLINTGpi); |
---|
| 325 | |
---|
| 326 | setCharacteristic(0); |
---|
| 327 | modpk pk= modpk (p,k); |
---|
| 328 | CanonicalForm liftedSmallestFactor= convertFmpz_poly_t2FacCF ((fmpz_poly_t &)liftedFactors->p[0],Variable (1)); |
---|
| 329 | |
---|
| 330 | CanonicalForm otherFactor= convertFmpz_poly_t2FacCF ((fmpz_poly_t &)liftedFactors->p[1],Variable (1)); |
---|
| 331 | CanonicalForm test= pk (otherFactor*liftedSmallestFactor); |
---|
| 332 | |
---|
| 333 | Off (SW_SYMMETRIC_FF); |
---|
| 334 | liftedSmallestFactor= pk (liftedSmallestFactor); |
---|
| 335 | if (i == 2) |
---|
| 336 | liftedSmallestFactor= liftedSmallestFactor (eval[0],1); |
---|
| 337 | else |
---|
| 338 | liftedSmallestFactor= liftedSmallestFactor (eval[1],1); |
---|
| 339 | |
---|
| 340 | On (SW_SYMMETRIC_FF); |
---|
| 341 | CFMatrix M= CFMatrix (s, s); |
---|
| 342 | M(s,s)= power (CanonicalForm (p), k); |
---|
| 343 | for (int j= 1; j < s; j++) |
---|
| 344 | { |
---|
| 345 | M (j,j)= 1; |
---|
| 346 | M (j+1,j)= -liftedSmallestFactor; |
---|
| 347 | } |
---|
| 348 | |
---|
| 349 | mat_ZZ NTLM= *convertFacCFMatrix2NTLmat_ZZ (M); |
---|
| 350 | |
---|
| 351 | ZZ det; |
---|
| 352 | |
---|
| 353 | transpose (NTLM, NTLM); |
---|
| 354 | (void) LLL (det, NTLM, 1L, 1L); //use floating point LLL ? |
---|
| 355 | transpose (NTLM, NTLM); |
---|
| 356 | M= *convertNTLmat_ZZ2FacCFMatrix (NTLM); |
---|
| 357 | |
---|
| 358 | mipo= 0; |
---|
| 359 | for (int j= 1; j <= s; j++) |
---|
| 360 | mipo += M (j,1)*power (x,s-j); |
---|
| 361 | |
---|
| 362 | CFFList mipoFactors= factorize (mipo); |
---|
| 363 | mipoFactors.removeFirst(); |
---|
| 364 | |
---|
| 365 | fmpz_poly_clear (v[0]); |
---|
| 366 | fmpz_poly_clear (v[1]); |
---|
| 367 | fmpz_poly_clear (w[0]); |
---|
| 368 | fmpz_poly_clear (w[1]); |
---|
| 369 | delete [] v; |
---|
| 370 | delete [] w; |
---|
| 371 | delete [] link; |
---|
| 372 | fmpz_poly_factor_clear (liftedFactors); |
---|
| 373 | |
---|
| 374 | if (mipoFactors.length() > 1 || |
---|
| 375 | (mipoFactors.length() == 1 && mipoFactors.getFirst().exp() > 1)) |
---|
| 376 | { |
---|
| 377 | if (i+1 >= loop && ((loop-i == 1) || (loop-i==2 && mipos[0].isZero()))) |
---|
| 378 | { |
---|
| 379 | rec=true; |
---|
| 380 | printf ("goto2\n"); |
---|
| 381 | goto differentevalpoint; |
---|
| 382 | //TODO check if mipoFactors has length 1 and multiplicity 1 - if not choose a new point! |
---|
| 383 | } |
---|
| 384 | } |
---|
| 385 | else |
---|
| 386 | mipos[loop-i-1]= mipo; |
---|
| 387 | } |
---|
| 388 | |
---|
| 389 | On (SW_RATIONAL); |
---|
| 390 | if (xValid && yValid && !mipos[0].isZero() && !mipos[1].isZero()) |
---|
| 391 | { |
---|
| 392 | if (maxNorm (mipos[0]) < maxNorm (mipos[1])) |
---|
| 393 | alpha= rootOf (mipos[0]); |
---|
| 394 | else |
---|
| 395 | alpha= rootOf (mipos[1]); |
---|
| 396 | } |
---|
| 397 | else |
---|
| 398 | alpha= rootOf (mipo); |
---|
| 399 | |
---|
| 400 | for (i= 0; i < mipos.size(); i++) |
---|
| 401 | { |
---|
| 402 | out_cf ("mipos= ", mipos [i], "\n"); |
---|
| 403 | out_cf ("maxNorm mipo= ", maxNorm (mipos[i]), "\n"); |
---|
| 404 | } |
---|
| 405 | |
---|
| 406 | CanonicalForm F1; |
---|
| 407 | CFFList QaF1Factors; |
---|
| 408 | int wrongMipo= 0; |
---|
| 409 | if (xValid && yValid) |
---|
| 410 | { |
---|
| 411 | if (degree (F,1) > minTdeg) |
---|
| 412 | F1= F (eval[1], 2); |
---|
| 413 | else |
---|
| 414 | F1= F (eval[0], 1); |
---|
| 415 | } |
---|
| 416 | else if (xValid) |
---|
| 417 | F1= F (eval[1], 2); |
---|
| 418 | else |
---|
| 419 | F1= F (eval[0], 1); |
---|
| 420 | |
---|
| 421 | QaF1Factors= factorize (F1, alpha); |
---|
| 422 | if (QaF1Factors.getFirst().factor().inCoeffDomain()) |
---|
| 423 | QaF1Factors.removeFirst(); |
---|
| 424 | out_cf ("mipo0= ", getMipo (alpha), "\n"); |
---|
| 425 | for (CFFListIterator iter= QaF1Factors; iter.hasItem(); iter++) |
---|
| 426 | { |
---|
| 427 | out_cf ("QaF1Factors= ", iter.getItem().factor(), "\n"); |
---|
| 428 | if (degree (iter.getItem().factor()) > minTdeg) |
---|
| 429 | wrongMipo++; |
---|
| 430 | } |
---|
| 431 | |
---|
| 432 | if (wrongMipo == QaF1Factors.length()) |
---|
| 433 | { |
---|
| 434 | if (xValid && yValid) |
---|
| 435 | { |
---|
| 436 | if (mipo==mipos[0]) |
---|
| 437 | alpha= rootOf (mipos[1]); |
---|
| 438 | else |
---|
| 439 | alpha= rootOf (mipos[0]); |
---|
| 440 | } |
---|
| 441 | |
---|
| 442 | wrongMipo= 0; |
---|
| 443 | out_cf ("mipo1= ", getMipo (alpha), "\n"); |
---|
| 444 | QaF1Factors= factorize (F1, alpha); |
---|
| 445 | if (QaF1Factors.getFirst().factor().inCoeffDomain()) |
---|
| 446 | QaF1Factors.removeFirst(); |
---|
| 447 | for (CFFListIterator iter= QaF1Factors; iter.hasItem(); iter++) |
---|
| 448 | { |
---|
| 449 | out_cf ("QaF1Factors= ", iter.getItem().factor(), "\n"); |
---|
| 450 | if (degree (iter.getItem().factor()) > minTdeg) |
---|
| 451 | wrongMipo++; |
---|
| 452 | } |
---|
| 453 | if (wrongMipo == QaF1Factors.length()) |
---|
| 454 | { |
---|
| 455 | rec= true; |
---|
| 456 | printf ("goto30\n"); |
---|
| 457 | goto differentevalpoint; |
---|
| 458 | } |
---|
| 459 | } |
---|
| 460 | |
---|
| 461 | CanonicalForm A= F; |
---|
| 462 | CanonicalForm Aeval= F1; |
---|
| 463 | |
---|
| 464 | out_cf ("F1= ", F1, "\n"); |
---|
| 465 | A *= bCommonDen (A); |
---|
| 466 | A= A (y + eval[1], y); //TODO find right evaluation point and swap if necessary |
---|
| 467 | |
---|
| 468 | out_cf ("A= ", A, "\n"); |
---|
| 469 | out_cf ("A[0]= ", A(0,y), "\n"); |
---|
| 470 | int liftBound= degree (A,y) + 1; |
---|
| 471 | |
---|
| 472 | modpk b= modpk(); |
---|
| 473 | |
---|
| 474 | //bool mipoHasDen= false; |
---|
| 475 | CanonicalForm den= 1; |
---|
| 476 | |
---|
| 477 | mipo= getMipo (alpha); |
---|
| 478 | |
---|
| 479 | CFList uniFactors; |
---|
| 480 | for (CFFListIterator iter=QaF1Factors; iter.hasItem(); iter++) |
---|
| 481 | { |
---|
| 482 | uniFactors.append (iter.getItem().factor()); |
---|
| 483 | out_cf ("uniFactors.getLast()= ", uniFactors.getLast(), "\n"); |
---|
| 484 | } |
---|
| 485 | |
---|
| 486 | |
---|
| 487 | A /= Lc (Aeval); |
---|
| 488 | //mipoHasDen= !bCommonDen(mipo).isOne(); |
---|
| 489 | //mipo *= bCommonDen (mipo); |
---|
| 490 | ZZX NTLmipo= convertFacCF2NTLZZX (mipo); |
---|
| 491 | ZZX NTLLcf= convertFacCF2NTLZZX (Lc (A*bCommonDen (A))); |
---|
| 492 | ZZ NTLf= resultant (NTLmipo, NTLLcf); |
---|
| 493 | ZZ NTLD= discriminant (NTLmipo); |
---|
| 494 | den= abs (convertZZ2CF (NTLD*NTLf)); |
---|
| 495 | |
---|
| 496 | // make factors elements of Z(a)[x] disable for modularDiophant |
---|
| 497 | CanonicalForm multiplier= 1; |
---|
| 498 | for (CFListIterator i= uniFactors; i.hasItem(); i++) |
---|
| 499 | { |
---|
| 500 | multiplier *= bCommonDen (i.getItem()); |
---|
| 501 | i.getItem()= i.getItem()*bCommonDen(i.getItem()); |
---|
| 502 | } |
---|
| 503 | A *= multiplier; |
---|
| 504 | A *= bCommonDen (A); |
---|
| 505 | |
---|
| 506 | Off (SW_RATIONAL); |
---|
| 507 | int ii= 0; |
---|
| 508 | CanonicalForm discMipo= convertZZ2CF (NTLD); |
---|
| 509 | findGoodPrime (F*discMipo,ii); |
---|
| 510 | findGoodPrime (Aeval*discMipo,ii); |
---|
| 511 | findGoodPrime (A*discMipo,ii); |
---|
| 512 | |
---|
| 513 | int pp=cf_getBigPrime(ii); |
---|
| 514 | b = coeffBound( A, pp, mipo ); |
---|
| 515 | modpk bb= coeffBound (Aeval, pp, mipo); |
---|
| 516 | if (bb.getk() > b.getk() ) b=bb; |
---|
| 517 | bb= coeffBound (F, pp, mipo); |
---|
| 518 | if (bb.getk() > b.getk() ) b=bb; |
---|
| 519 | |
---|
| 520 | ExtensionInfo dummy= ExtensionInfo (alpha, false); |
---|
| 521 | DegreePattern degs= DegreePattern (uniFactors); |
---|
| 522 | |
---|
| 523 | bool earlySuccess= false; |
---|
| 524 | CFList earlyFactors; |
---|
| 525 | TIMING_START (fac_bi_hensel_lift); |
---|
| 526 | uniFactors= henselLiftAndEarly |
---|
| 527 | (A, earlySuccess, earlyFactors, degs, liftBound, |
---|
| 528 | uniFactors, dummy, eval[1], b, den); |
---|
| 529 | TIMING_END_AND_PRINT (fac_bi_hensel_lift, |
---|
| 530 | "time for bivariate hensel lifting over Q: "); |
---|
| 531 | DEBOUTLN (cerr, "lifted factors= " << uniFactors); |
---|
| 532 | |
---|
| 533 | CanonicalForm MODl= power (y, liftBound); //TODO |
---|
| 534 | |
---|
[dea3d2] | 535 | On (SW_RATIONAL); |
---|
[17a710] | 536 | A *= bCommonDen (A); |
---|
| 537 | Off (SW_RATIONAL); |
---|
| 538 | |
---|
| 539 | printf ("earlyFactors.length()= %d\n", earlyFactors.length()); |
---|
| 540 | CFList biFactors; |
---|
| 541 | |
---|
| 542 | TIMING_START (fac_bi_factor_recombination); |
---|
| 543 | biFactors= factorRecombination (uniFactors, A, MODl, degs, 1, |
---|
| 544 | uniFactors.length()/2, b, den); |
---|
| 545 | TIMING_END_AND_PRINT (fac_bi_factor_recombination, |
---|
| 546 | "time for bivariate factor recombination over Q: "); |
---|
| 547 | |
---|
| 548 | On (SW_RATIONAL); |
---|
| 549 | |
---|
| 550 | if (earlySuccess) |
---|
| 551 | biFactors= Union (earlyFactors, biFactors); |
---|
| 552 | else if (!earlySuccess && degs.getLength() == 1) |
---|
| 553 | biFactors= earlyFactors; |
---|
| 554 | |
---|
| 555 | for (CFListIterator i= biFactors; i.hasItem(); i++) |
---|
| 556 | i.getItem()= i.getItem() (y - eval[1], y); //TODO |
---|
| 557 | |
---|
| 558 | bool swap= false; |
---|
| 559 | bool swap2= false; |
---|
| 560 | appendSwapDecompress (biFactors, CFList(), CFList(), swap, swap2, CFMap()); |
---|
| 561 | if (isOn (SW_RATIONAL)) |
---|
| 562 | normalize (biFactors); |
---|
| 563 | |
---|
| 564 | CFAFList result; |
---|
| 565 | bool found= false; |
---|
| 566 | |
---|
| 567 | out_cf ("mipo= ", mipo, "\n"); |
---|
| 568 | printf ("minTdeg= %d\n", minTdeg); |
---|
| 569 | for (CFListIterator iter= biFactors; iter.hasItem(); iter++) |
---|
| 570 | { |
---|
| 571 | out_cf ("biFactors= ", iter.getItem(), "\n"); |
---|
| 572 | printf ("totaldegree ()= %d\n", totaldegree (iter.getItem())); |
---|
| 573 | if (totaldegree (iter.getItem()) == minTdeg) |
---|
| 574 | { |
---|
| 575 | result= CFAFList (CFAFactor (iter.getItem(), getMipo (alpha), 1)); |
---|
| 576 | found= true; |
---|
| 577 | break; |
---|
| 578 | } |
---|
| 579 | } |
---|
[dea3d2] | 580 | |
---|
[17a710] | 581 | if (found) |
---|
| 582 | { |
---|
| 583 | printf ("thisexitexit\n\n"); |
---|
| 584 | return result; |
---|
| 585 | } |
---|
| 586 | |
---|
| 587 | /* A *= bCommonDen (A); |
---|
| 588 | A= A (y + evaluation, y); |
---|
| 589 | |
---|
| 590 | int liftBound= degree (A, y) + 1; |
---|
| 591 | |
---|
| 592 | modpk b= modpk(); |
---|
| 593 | bool mipoHasDen= false; |
---|
| 594 | CanonicalForm den= 1; |
---|
| 595 | |
---|
| 596 | if (!extension) |
---|
| 597 | { |
---|
| 598 | Off (SW_RATIONAL); |
---|
| 599 | int i= 0; |
---|
| 600 | findGoodPrime(F,i); |
---|
| 601 | findGoodPrime(Aeval,i); |
---|
| 602 | findGoodPrime(A,i); |
---|
| 603 | if (i >= cf_getNumBigPrimes()) |
---|
| 604 | printf ("out of primes\n"); //TODO exit |
---|
| 605 | |
---|
| 606 | int p=cf_getBigPrime(i); |
---|
| 607 | b = coeffBound( A, p ); |
---|
| 608 | modpk bb= coeffBound (Aeval, p); |
---|
| 609 | if (bb.getk() > b.getk() ) b=bb; |
---|
| 610 | bb= coeffBound (F, p); |
---|
| 611 | if (bb.getk() > b.getk() ) b=bb; |
---|
| 612 | } |
---|
| 613 | else |
---|
| 614 | { |
---|
| 615 | A /= Lc (Aeval); |
---|
| 616 | mipoHasDen= !bCommonDen(mipo).isOne(); |
---|
| 617 | mipo *= bCommonDen (mipo); |
---|
| 618 | ZZX NTLmipo= convertFacCF2NTLZZX (mipo); |
---|
| 619 | ZZX NTLLcf= convertFacCF2NTLZZX (Lc (A*bCommonDen (A))); |
---|
| 620 | ZZ NTLf= resultant (NTLmipo, NTLLcf); |
---|
| 621 | ZZ NTLD= discriminant (NTLmipo); |
---|
| 622 | den= abs (convertZZ2CF (NTLD*NTLf)); |
---|
| 623 | |
---|
| 624 | // make factors elements of Z(a)[x] disable for modularDiophant |
---|
| 625 | CanonicalForm multiplier= 1; |
---|
| 626 | for (CFListIterator i= uniFactors; i.hasItem(); i++) |
---|
| 627 | { |
---|
| 628 | multiplier *= bCommonDen (i.getItem()); |
---|
| 629 | i.getItem()= i.getItem()*bCommonDen(i.getItem()); |
---|
| 630 | } |
---|
| 631 | A *= multiplier; |
---|
| 632 | A *= bCommonDen (A); |
---|
| 633 | |
---|
| 634 | Off (SW_RATIONAL); |
---|
| 635 | int i= 0; |
---|
| 636 | CanonicalForm discMipo= convertZZ2CF (NTLD); |
---|
| 637 | findGoodPrime (F*discMipo,i); |
---|
| 638 | findGoodPrime (Aeval*discMipo,i); |
---|
| 639 | findGoodPrime (A*discMipo,i); |
---|
| 640 | |
---|
| 641 | int p=cf_getBigPrime(i); |
---|
| 642 | b = coeffBound( A, p, mipo ); |
---|
| 643 | modpk bb= coeffBound (Aeval, p, mipo); |
---|
| 644 | if (bb.getk() > b.getk() ) b=bb; |
---|
| 645 | bb= coeffBound (F, p, mipo); |
---|
| 646 | if (bb.getk() > b.getk() ) b=bb; |
---|
| 647 | } |
---|
| 648 | |
---|
| 649 | ExtensionInfo dummy= ExtensionInfo (false); |
---|
| 650 | if (extension) |
---|
| 651 | dummy= ExtensionInfo (v, false); |
---|
| 652 | bool earlySuccess= false; |
---|
| 653 | CFList earlyFactors; |
---|
| 654 | TIMING_START (fac_bi_hensel_lift); |
---|
| 655 | uniFactors= henselLiftAndEarly |
---|
| 656 | (A, earlySuccess, earlyFactors, degs, liftBound, |
---|
| 657 | uniFactors, dummy, evaluation, b, den); |
---|
| 658 | TIMING_END_AND_PRINT (fac_bi_hensel_lift, |
---|
| 659 | "time for bivariate hensel lifting over Q: "); |
---|
| 660 | DEBOUTLN (cerr, "lifted factors= " << uniFactors); |
---|
| 661 | |
---|
| 662 | CanonicalForm MODl= power (y, liftBound); |
---|
| 663 | |
---|
| 664 | if (mipoHasDen) |
---|
| 665 | { |
---|
| 666 | Variable vv; |
---|
| 667 | for (CFListIterator iter= uniFactors; iter.hasItem(); iter++) |
---|
| 668 | if (hasFirstAlgVar (iter.getItem(), vv)) |
---|
| 669 | break; |
---|
| 670 | for (CFListIterator iter= uniFactors; iter.hasItem(); iter++) |
---|
| 671 | iter.getItem()= replacevar (iter.getItem(), vv, v); |
---|
| 672 | } |
---|
| 673 | |
---|
| 674 | On (SW_RATIONAL); |
---|
| 675 | A *= bCommonDen (A); |
---|
| 676 | Off (SW_RATIONAL); |
---|
| 677 | |
---|
| 678 | TIMING_START (fac_bi_factor_recombination); |
---|
| 679 | factors= factorRecombination (uniFactors, A, MODl, degs, 1, |
---|
| 680 | uniFactors.length()/2, b, den); |
---|
| 681 | TIMING_END_AND_PRINT (fac_bi_factor_recombination, |
---|
| 682 | "time for bivariate factor recombination over Q: "); |
---|
| 683 | |
---|
| 684 | On (SW_RATIONAL); |
---|
| 685 | |
---|
| 686 | if (earlySuccess) |
---|
| 687 | factors= Union (earlyFactors, factors); |
---|
| 688 | else if (!earlySuccess && degs.getLength() == 1) |
---|
| 689 | factors= earlyFactors; |
---|
| 690 | |
---|
| 691 | for (CFListIterator i= factors; i.hasItem(); i++) |
---|
| 692 | i.getItem()= i.getItem() (y - evaluation, y); |
---|
| 693 | |
---|
| 694 | appendSwapDecompress (factors, conv (contentAxFactors), |
---|
| 695 | conv (contentAyFactors), swap, swap2, N); |
---|
| 696 | if (isOn (SW_RATIONAL)) |
---|
| 697 | normalize (factors);*/ |
---|
| 698 | |
---|
| 699 | TIMING_START (fac_Qa_factorize); |
---|
| 700 | CFFList QaFactors= factorize (F, alpha); //TODO lift these factors |
---|
| 701 | TIMING_END_AND_PRINT (fac_Qa_factorize, "time to factorize over Qa: "); |
---|
| 702 | |
---|
| 703 | /*mipo= getMipo (alpha); |
---|
| 704 | out_cf ("maxNorm (mipo)= ", maxNorm (mipo), "\n"); |
---|
[160ec6] | 705 | QaFactors.append (CFFactor (mipo, 1)); //last factor is the minimal polynomial that defines the extension |
---|
[17a710] | 706 | if (degree (mipo) < 3) |
---|
| 707 | printf ("scheissescheissescheissescheisse\n");*/ |
---|
| 708 | printf ("minTdeg= %d\n", minTdeg); |
---|
| 709 | //CFAFList result; |
---|
| 710 | //bool found= false; |
---|
| 711 | out_cf ("mipo= ", getMipo (alpha), "\n"); |
---|
| 712 | for (CFFListIterator iter= QaFactors; iter.hasItem(); iter++) |
---|
| 713 | { |
---|
| 714 | out_cf ("QaFactors= ", iter.getItem().factor(), "\n"); |
---|
| 715 | printf ("totaldegree ()= %d\n", totaldegree (iter.getItem().factor())); |
---|
| 716 | if (totaldegree (iter.getItem().factor()) == minTdeg) |
---|
| 717 | { |
---|
| 718 | result= CFAFList (CFAFactor (iter.getItem().factor(), getMipo (alpha), 1)); |
---|
| 719 | found= true; |
---|
| 720 | break; |
---|
| 721 | } |
---|
| 722 | } |
---|
| 723 | if (!found && xValid && yValid) |
---|
| 724 | { |
---|
| 725 | if (mipo == mipos [0]) |
---|
| 726 | mipo= mipos[1]; |
---|
| 727 | else |
---|
| 728 | mipo= mipos[0]; |
---|
| 729 | alpha= rootOf (mipo); |
---|
| 730 | |
---|
| 731 | QaFactors= factorize (F, alpha); |
---|
| 732 | |
---|
| 733 | for (CFFListIterator iter= QaFactors; iter.hasItem(); iter++) |
---|
| 734 | { |
---|
| 735 | out_cf ("QaFactors= ", iter.getItem().factor(), "\n"); |
---|
| 736 | printf ("totaldegree ()= %d\n", totaldegree (iter.getItem().factor())); |
---|
| 737 | if (totaldegree (iter.getItem().factor()) == minTdeg) |
---|
| 738 | { |
---|
| 739 | result= CFAFList (CFAFactor (iter.getItem().factor(), getMipo (alpha), 1)); |
---|
| 740 | found= true; |
---|
| 741 | break; |
---|
| 742 | } |
---|
| 743 | } |
---|
| 744 | if (!found) |
---|
| 745 | { |
---|
| 746 | rec= true; |
---|
| 747 | printf ("goto31\n"); |
---|
| 748 | goto differentevalpoint; |
---|
| 749 | } |
---|
| 750 | } |
---|
| 751 | else if (!found) |
---|
| 752 | { |
---|
| 753 | rec= true; |
---|
| 754 | printf ("goto32\n"); |
---|
| 755 | goto differentevalpoint; |
---|
| 756 | } |
---|
| 757 | |
---|
| 758 | return result; |
---|
[dea3d2] | 759 | } |
---|
| 760 | #endif |
---|
[160ec6] | 761 | #endif |
---|
[dea3d2] | 762 | |
---|
| 763 | |
---|