Changeset f7a4e9 in git
 Timestamp:
 Apr 14, 2012, 1:39:51 PM (11 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
 Children:
 c244d04150a9fd45331938d2c9fb6f89985acb05
 Parents:
 a9a6dcbcbc4f0be99516f14eae42b1776235ca9f
 gitauthor:
 Martin Lee <martinlee84@web.de>20120414 13:39:51+02:00
 gitcommitter:
 Martin Lee <martinlee84@web.de>20120507 14:14:17+02:00
 Location:
 factory
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

factory/cf_gcd.cc
ra9a6dcb rf7a4e9 757 757 else if (!fc_and_gc_Univariate) 758 758 { 759 if ( 760 isOn(SW_USE_CHINREM_GCD) 761 && (isPurePoly_m(fc)) && (isPurePoly_m(gc)) 762 ) 763 { 764 #if 0 765 if ( p1 == fc.level() ) 766 fc = chinrem_gcd( fc, gc ); 767 else 768 { 769 fc = replacevar( fc, Variable(p1), Variable(mp) ); 770 gc = replacevar( gc, Variable(p1), Variable(mp) ); 771 fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) ); 772 } 773 #else 774 fc = chinrem_gcd( fc, gc); 775 #endif 776 } 777 else if ( isOn( SW_USE_EZGCD ) ) 778 { 779 if ( pe == 1 ) 759 if ( isOn( SW_USE_EZGCD ) ) 760 { 761 fc= ezgcd (fc, gc); 762 /*if ( pe == 1 ) 780 763 fc = ezgcd( fc, gc ); 781 764 else if ( pe > 0 )// no variable at position 1 … … 791 774 gc = swapvar( gc, Variable(pe), Variable(1) ); 792 775 fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) ); 776 }*/ 777 } 778 else if ( 779 isOn(SW_USE_CHINREM_GCD) 780 && (isPurePoly_m(fc)) && (isPurePoly_m(gc)) 781 ) 782 { 783 #if 0 784 if ( p1 == fc.level() ) 785 fc = chinrem_gcd( fc, gc ); 786 else 787 { 788 fc = replacevar( fc, Variable(p1), Variable(mp) ); 789 gc = replacevar( gc, Variable(p1), Variable(mp) ); 790 fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) ); 793 791 } 792 #else 793 fc = chinrem_gcd( fc, gc); 794 #endif 794 795 } 795 796 else 
factory/fac_ezgcd.cc
ra9a6dcb rf7a4e9 16 16 #include "fac_distrib.h" 17 17 #include "templates/ftmpl_functions.h" 18 #include "cf_map.h" 19 #include "facHensel.h" 20 21 static 22 int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M, 23 CFMap & N, int& both_non_zero) 24 { 25 int n= tmax (F.level(), G.level()); 26 int * degsf= new int [n + 1]; 27 int * degsg= new int [n + 1]; 28 29 for (int i = 0; i <= n; i++) 30 degsf[i]= degsg[i]= 0; 31 32 degsf= degrees (F, degsf); 33 degsg= degrees (G, degsg); 34 35 both_non_zero= 0; 36 int f_zero= 0; 37 int g_zero= 0; 38 39 for (int i= 1; i <= n; i++) 40 { 41 if (degsf[i] != 0 && degsg[i] != 0) 42 { 43 both_non_zero++; 44 continue; 45 } 46 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 47 { 48 f_zero++; 49 continue; 50 } 51 if (degsg[i] == 0 && degsf[i] && i <= F.level()) 52 { 53 g_zero++; 54 continue; 55 } 56 } 57 58 if (both_non_zero == 0) 59 { 60 delete [] degsf; 61 delete [] degsg; 62 return 0; 63 } 64 65 // map Variables which do not occur in both polynomials to higher levels 66 int k= 1; 67 int l= 1; 68 for (int i= 1; i <= n; i++) 69 { 70 if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level()) 71 { 72 if (k + both_non_zero != i) 73 { 74 M.newpair (Variable (i), Variable (k + both_non_zero)); 75 N.newpair (Variable (k + both_non_zero), Variable (i)); 76 } 77 k++; 78 } 79 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 80 { 81 if (l + g_zero + both_non_zero != i) 82 { 83 M.newpair (Variable (i), Variable (l + g_zero + both_non_zero)); 84 N.newpair (Variable (l + g_zero + both_non_zero), Variable (i)); 85 } 86 l++; 87 } 88 } 89 90 // sort Variables x_{i} in decreasing order of 91 // min(deg_{x_{i}}(f),deg_{x_{i}}(g)) 92 int m= tmin (F.level(), G.level()); 93 int max_min_deg; 94 k= both_non_zero; 95 l= 0; 96 int i= 1; 97 while (k > 0) 98 { 99 max_min_deg= tmin (degsf[i], degsg[i]); 100 while (max_min_deg == 0) 101 { 102 i++; 103 max_min_deg= tmin (degsf[i], degsg[i]); 104 } 105 for (int j= i + 1; j <= m; j++) 106 { 107 if ((tmin (degsf[j],degsg[j]) < max_min_deg) && 108 (tmin (degsf[j], degsg[j]) != 0)) 109 { 110 max_min_deg= tmin (degsf[j], degsg[j]); 111 l= j; 112 } 113 } 114 115 if (l != 0) 116 { 117 if (l != k) 118 { 119 M.newpair (Variable (l), Variable(k)); 120 N.newpair (Variable (k), Variable(l)); 121 degsf[l]= 0; 122 degsg[l]= 0; 123 l= 0; 124 } 125 else 126 { 127 degsf[l]= 0; 128 degsg[l]= 0; 129 l= 0; 130 } 131 } 132 else if (l == 0) 133 { 134 if (i != k) 135 { 136 M.newpair (Variable (i), Variable (k)); 137 N.newpair (Variable (k), Variable (i)); 138 degsf[i]= 0; 139 degsg[i]= 0; 140 } 141 else 142 { 143 degsf[i]= 0; 144 degsg[i]= 0; 145 } 146 i++; 147 } 148 k; 149 } 150 151 delete [] degsf; 152 delete [] degsg; 153 154 return both_non_zero; 155 } 156 157 static inline 158 CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval, 159 const CFList& evaluation) 160 { 161 CanonicalForm A= F; 162 int k= 2; 163 for (CFListIterator i= evaluation; i.hasItem(); i++, k++) 164 A= A (Variable (k) + i.getItem(), k); 165 166 CanonicalForm buf= A; 167 Feval= CFList(); 168 Feval.append (buf); 169 for (k= evaluation.length() + 1; k > 2; k) 170 { 171 buf= mod (buf, Variable (k)); 172 Feval.insert (buf); 173 } 174 return A; 175 } 176 177 static inline 178 CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation) 179 { 180 int l= evaluation.length() + 1; 181 CanonicalForm result= F; 182 CFListIterator j= evaluation; 183 for (int i= 2; i < l + 1; i++, j++) 184 { 185 if (F.level() < i) 186 continue; 187 result= result (Variable (i)  j.getItem(), i); 188 } 189 return result; 190 } 191 192 static inline 193 Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M, 194 CFMap & N, const Evaluation& A) 195 { 196 int n= F.level(); 197 int * degsf= new int [n + 1]; 198 199 for (int i = 0; i <= n; i++) 200 degsf[i]= 0; 201 202 degsf= degrees (F, degsf); 203 204 Evaluation result= Evaluation (A.min(), A.max()); 205 ASSERT (A.min() == 2, "expected A.min() == 2"); 206 ASSERT (A.max() == n, "expected A.max() == n"); 207 int max_deg; 208 int k= n; 209 int l= 1; 210 int i= 2; 211 int pos= 2; 212 while (k > 1) 213 { 214 max_deg= degsf [i]; 215 while (max_deg == 0) 216 { 217 i++; 218 max_deg= degsf [i]; 219 } 220 l= i; 221 for (int j= i + 1; j <= n; j++) 222 { 223 if (degsf[j] > max_deg) 224 { 225 max_deg= degsf[j]; 226 l= j; 227 } 228 } 229 230 if (l <= n) 231 { 232 if (l != pos) 233 { 234 result.setValue (pos, A [l]); 235 M.newpair (Variable (l), Variable (pos)); 236 N.newpair (Variable (pos), Variable (l)); 237 degsf[l]= 0; 238 l= 2; 239 if (k == 2 && n == 3) 240 { 241 result.setValue (l, A [pos]); 242 M.newpair (Variable (pos), Variable (l)); 243 N.newpair (Variable (l), Variable (pos)); 244 degsf[pos]= 0; 245 } 246 } 247 else 248 { 249 result.setValue (l, A [l]); 250 degsf [l]= 0; 251 } 252 } 253 pos++; 254 k; 255 l= 2; 256 } 257 258 delete [] degsf; 259 260 return result; 261 } 262 263 static inline 264 int Hensel (const CanonicalForm & UU, CFArray & G, const Evaluation & AA, 265 const CFArray& LeadCoeffs ) 266 { 267 CFList factors; 268 factors.append (G[1]); 269 factors.append (G[2]); 270 271 CFMap NN, MM; 272 Evaluation A= optimize4Lift (UU, MM, NN, AA); 273 274 CanonicalForm U= MM (UU); 275 CFArray LCs= CFArray (1,2); 276 LCs [1]= MM (LeadCoeffs [1]); 277 LCs [2]= MM (LeadCoeffs [2]); 278 279 CFList evaluation; 280 for (int i= A.min(); i <= A.max(); i++) 281 evaluation.append (A [i]); 282 CFList UEval; 283 CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation); 284 285 if (size (shiftedU)/getNumVars (U) > 500) 286 return 1; 287 288 CFArray shiftedLCs= CFArray (2); 289 CFList shiftedLCsEval1, shiftedLCsEval2; 290 shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation); 291 shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation); 292 factors.insert (1); 293 int liftBound= degree (UEval.getLast(), 2) + 1; 294 CFArray Pi; 295 CFMatrix M= CFMatrix (liftBound, factors.length()  1); 296 CFList diophant; 297 CFArray lcs= CFArray (2); 298 lcs [0]= shiftedLCsEval1.getFirst(); 299 lcs [1]= shiftedLCsEval2.getFirst(); 300 nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M, 301 lcs, false); 302 303 for (CFListIterator i= factors; i.hasItem(); i++) 304 { 305 if (!fdivides (i.getItem(), UEval.getFirst())) 306 return 0; 307 } 308 309 int * liftBounds; 310 bool noOneToOne= false; 311 if (U.level() > 2) 312 { 313 liftBounds= new int [U.level()  1]; /* index: 0.. U.level()2 */ 314 liftBounds[0]= liftBound; 315 for (int i= 1; i < U.level()  1; i++) 316 liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1; 317 factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level()  1, 318 false, shiftedLCsEval1, shiftedLCsEval2, Pi, 319 diophant, noOneToOne); 320 delete [] liftBounds; 321 if (noOneToOne) 322 return 0; 323 } 324 G[1]= factors.getFirst(); 325 G[2]= factors.getLast(); 326 G[1]= myReverseShift (G[1], evaluation); 327 G[2]= myReverseShift (G[2], evaluation); 328 G[1]= NN (G[1]); 329 G[2]= NN (G[2]); 330 return 1; 331 } 18 332 19 333 static void findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG ); … … 37 351 ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal ) 38 352 { 39 CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD ;353 CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD, cand, result; 40 354 CFArray DD( 1, 2 ), lcDD( 1, 2 ); 41 355 int degF, degG, delta, t; … … 45 359 modpk bound; 46 360 361 F= FF; 362 G= GG; 363 364 CFMap M,N; 365 int smallestDegLev; 366 int best_level= compress4EZGCD (F, G, M, N, smallestDegLev); 367 368 if (best_level == 0) return G.genOne(); 369 370 F= M (F); 371 G= M (G); 372 373 47 374 DEBINCLEVEL( cerr, "ezgcd" ); 48 375 DEBOUTLN( cerr, "FF = " << FF ); 49 376 DEBOUTLN( cerr, "GG = " << GG ); 50 f = content( F F, x ); g = content( GG, x ); d = gcd( f, g );377 f = content( F, x ); g = content( G, x ); d = gcd( f, g ); 51 378 DEBOUTLN( cerr, "f = " << f ); 52 379 DEBOUTLN( cerr, "g = " << g ); 53 F = FF / f; G = GG /g;380 F /= f; G /= g; 54 381 if ( F.isUnivariate() && G.isUnivariate() ) 55 382 { … … 57 384 if(F.mvar()==G.mvar()) 58 385 d*=gcd(F,G); 59 return d;386 return N (d); 60 387 } 61 388 else if ( gcd_test_one( F, G, false ) ) 62 389 { 63 390 DEBDECLEVEL( cerr, "ezgcd" ); 64 return d;391 return N (d); 65 392 } 66 393 lcF = LC( F, x ); lcG = LC( G, x ); … … 69 396 degF = degree( F, x ); degG = degree( G, x ); 70 397 t = tmax( F.level(), G.level() ); 71 bound = findBound( F, G, lcF, lcG, degF, degG );398 //bound = findBound( F, G, lcF, lcG, degF, degG ); 72 399 if ( ! internal ) 73 400 b = REvaluation( 2, t, IntRandom( 25 ) ); 74 while ( ! gcdfound ) { 401 while ( ! gcdfound ) 402 { 75 403 /// > A2 76 404 DEBOUTLN( cerr, "search for evaluation, delta = " << delta ); … … 87 415 { 88 416 DEBDECLEVEL( cerr, "ezgcd" ); 89 return d;417 return N (d); 90 418 } 91 419 /// > A4 … … 98 426 { 99 427 DEBDECLEVEL( cerr, "ezgcd" ); 100 return d;428 return N (d); 101 429 } 102 430 if ( dd /*degree( Dbt )*/ == delta ) … … 110 438 DEBOUTLN( cerr, "now after A4, delta = " << delta ); 111 439 /// > A5 112 if ( degF <= degG && delta == degF && fdivides( F, G ))440 if (delta == degF) 113 441 { 442 if (degF <= degG && fdivides (F, G)) 443 { 114 444 DEBDECLEVEL( cerr, "ezgcd" ); 115 return d*F; 116 } 117 if ( degG < degF && delta == degG && fdivides( G, F ) ) 445 return N (d*F); 446 } 447 else 448 delta; 449 } 450 else if (delta == degG) 118 451 { 452 if (degG <= degF && fdivides( G, F )) 453 { 119 454 DEBDECLEVEL( cerr, "ezgcd" ); 120 return d*G; 121 } 122 if ( delta != degF && delta != degG ) { 455 return N (d*G); 456 } 457 else 458 delta; 459 } 460 if ( delta != degF && delta != degG ) 461 { 462 /// > A6 123 463 bool B_is_F; 124 /// > A6 125 CanonicalForm xxx; 464 CanonicalForm xxx1, xxx2; 465 CanonicalForm buf; 466 DD[1] = Fb / Db; 467 buf= Gb/Db; 468 xxx1 = gcd( DD[1], Db ); 469 xxx2 = gcd( buf, Db ); 470 if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) && 471 (size (F) <= size (G))) 472  (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain())) 473 { 474 B = F; 475 DD[2] = Db; 476 lcDD[1] = lcF; 477 lcDD[2] = lcD; 478 B_is_F = true; 479 } 480 else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) && 481 (size (G) < size (F))) 482  (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain())) 483 { 484 DD[1] = buf; 485 B = G; 486 DD[2] = Db; 487 lcDD[1] = lcG; 488 lcDD[2] = lcD; 489 B_is_F = false; 490 } 491 else 492 { 493 Off(SW_USE_EZGCD); 494 result=gcd( F, G ); 495 On(SW_USE_EZGCD); 496 return N (d*result); 497 } 498 /*CanonicalForm xxx; 126 499 //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) { 127 500 DD[1] = Fb / Db; … … 165 538 #endif 166 539 } 167 } 540 }*/ 168 541 /// > A7 169 542 DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) ); 170 543 DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) ); 171 bound = enlargeBound( B, DD[2], DD[1], bound );544 //bound = enlargeBound( B, DD[2], DD[1], bound ); 172 545 DEBOUTLN( cerr, "(hensel) B = " << B ); 173 546 DEBOUTLN( cerr, "(hensel) lcB = " << LC( B, Variable(1) ) ); … … 175 548 DEBOUTLN( cerr, "(hensel) DD = " << DD ); 176 549 DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD ); 177 gcdfound = Hensel( B, DD, lcDD, b, bound, x);550 gcdfound= Hensel (B*lcD, DD, b, lcDD); 178 551 DEBOUTLN( cerr, "(hensel finished) DD = " << DD ); 552 553 if (gcdfound == 1) 554 { 555 Off (SW_USE_EZGCD); 556 result= gcd (F,G); 557 On (SW_USE_EZGCD_P); 558 return N (d*result); 559 } 179 560 180 561 if (gcdfound) 181 562 { 182 CanonicalForm xxx=content(DD[2],Variable(1)); 183 CanonicalForm cand=DD[2] / xxx; //content(DD[2],Variable(1)); 184 #if 0 185 gcdfound= fdivides(cand,G) && fdivides(cand,F); 186 #else 187 if (B_is_F /*B==F*lcF*/) 188 { 189 DEBOUTLN( cerr, "(test) G: "<<G<<" % gcd:"<<cand<<" > " << G%cand ); 190 gcdfound= fdivides(cand,G); 191 } 192 else 193 { 194 DEBOUTLN( cerr, "(test) F: "<<F<<" % gcd:"<<cand<<" > " << F%cand); 195 gcdfound= fdivides(cand,F); 196 } 197 #endif 563 cand = DD[2] / content( DD[2], Variable(1) ); 564 gcdfound = fdivides( cand, G ) && fdivides ( cand, F ); 198 565 } 199 566 /// > A8 (gcdfound) 200 567 } 201 delta ++;568 delta; 202 569 } 203 570 /// > A9 204 571 DEBDECLEVEL( cerr, "ezgcd" ); 205 return d * DD[2] / content(DD[2],Variable(1));572 return N (d*cand); 206 573 } 207 574
Note: See TracChangeset
for help on using the changeset viewer.