- Timestamp:
- Apr 29, 2014, 5:54:00 PM (10 years ago)
- Branches:
- (u'spielwiese', '2a584933abf2a2d3082034c7586d38bb6de1a30a')
- Children:
- ea017ffdd4140ca784627c9a17bd45f8978c3818
- Parents:
- 24819474a270d2480ed2d9f37a85ab3f5eafee1c
- git-author:
- Martin Lee <martinlee84@web.de>2014-04-29 17:54:00+02:00
- git-committer:
- Martin Lee <martinlee84@web.de>2014-05-12 14:35:01+02:00
- Location:
- factory
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
factory/facAlgFunc.cc
r248194 r5355082 3 3 //////////////////////////////////////////////////////////// 4 4 5 // FACTORY - Includes 6 #include <factory.h> 7 // Factor - Includes 8 #include <tmpl_inst.h> 9 #include <helpstuff.h> 10 // Charset - Includes 11 #include "csutil.h" 12 #include "charset.h" 13 #include "reorder.h" 14 #include "algfactor.h" 15 // some CC's need this: 16 #include "alg_factor.h" 5 #include "config.h" 6 7 #include "cf_assert.h" 8 #include "debug.h" 9 10 #include "algext.h" 11 #include "cf_random.h" 12 #include "cf_generator.h" 13 #include "cf_irred.h" 14 #include "cf_iter.h" 15 #include "cf_util.h" 16 #include "cf_algorithm.h" 17 #include "cf_map.h" 18 #include "cfModResultant.h" 19 #include "cfCharSets.h" 20 #include "facAlgFunc.h" 17 21 18 22 void out_cf(const char *s1,const CanonicalForm &f,const char *s2); 19 23 20 #ifdef ALGFACTORDEBUG 21 # define DEBUGOUTPUT 22 #else 23 # undef DEBUGOUTPUT 24 #endif 25 26 #include <libfac/factor/debug.h> 24 static 25 CanonicalForm 26 Prem (const CanonicalForm& F, const CanonicalForm& G) 27 { 28 CanonicalForm f, g, l, test, lu, lv, t, retvalue; 29 int degF, degG, levelF, levelG; 30 bool reord; 31 Variable v, vg= G.mvar(); 32 33 if ( (levelF= F.level()) < (levelG= G.level())) 34 return F; 35 else 36 { 37 if ( levelF == levelG ) 38 { 39 f= F; 40 g= G; 41 reord= false; 42 v= F.mvar(); 43 } 44 else 45 { 46 v= Variable (levelF + 1); 47 f= swapvar (F, vg, v); 48 g= swapvar (G, vg, v); 49 reord= true; 50 } 51 degG= degree (g, v ); 52 degF= degree (f, v ); 53 if (degG <= degF) 54 { 55 l= LC (g); 56 g= g - l*power (v, degG); 57 } 58 else 59 l= 1; 60 while ((degG <= degF) && (!f.isZero())) 61 { 62 test= gcd (l, LC(f)); 63 lu= l / test; 64 lv= LC(f) / test; 65 66 t= power (v, degF - degG)*g*lv; 67 68 if (degF == 0) 69 f= 0; 70 else 71 f= f - LC (f)*power (v, degF); 72 73 f= lu*f - t; 74 degF= degree (f, v); 75 } 76 77 if (reord) 78 retvalue= swapvar (f, vg, v); 79 else 80 retvalue= f; 81 82 return retvalue; 83 } 84 } 85 86 static 87 CanonicalForm 88 Prem( const CanonicalForm &f, const CFList &L) 89 { 90 CanonicalForm rem= f; 91 CFListIterator i= L; 92 93 for (i.lastItem(); i.hasItem(); i--) 94 rem= normalize (Prem (rem, i.getItem())); 95 96 return rem; 97 } 98 99 static 100 CanonicalForm 101 Prem (const CanonicalForm &f, const CFList &L, const CFList &as) 102 { 103 CanonicalForm rem = f; 104 CFListIterator i = L; 105 for ( i.lastItem(); i.hasItem(); i-- ) 106 rem = Prem( rem, i.getItem() ); 107 return normalize (rem); //TODO better normalize in case as.length() == 1 && as.getFirst().level() == 1 ??? 108 } 109 110 CFFList 111 myappend( const CFFList & Inputlist, const CFFactor & TheFactor) 112 { 113 CFFList Outputlist ; 114 CFFactor copy; 115 CFFListIterator i; 116 int exp=0; 117 118 for ( i=Inputlist ; i.hasItem() ; i++ ) 119 { 120 copy = i.getItem(); 121 if ( copy.factor() == TheFactor.factor() ) 122 exp += copy.exp(); 123 else 124 Outputlist.append(copy); 125 } 126 Outputlist.append( CFFactor(TheFactor.factor(), exp + TheFactor.exp())); 127 return Outputlist; 128 } 129 130 CFFList 131 myUnion(const CFFList & Inputlist1,const CFFList & Inputlist2) 132 { 133 CFFList Outputlist; 134 CFFListIterator i; 135 136 for ( i=Inputlist1 ; i.hasItem() ; i++ ) 137 Outputlist = myappend(Outputlist, i.getItem() ); 138 for ( i=Inputlist2 ; i.hasItem() ; i++ ) 139 Outputlist = myappend(Outputlist, i.getItem() ); 140 141 return Outputlist; 142 } 27 143 28 144 /////////////////////////////////////////////////////////////// … … 34 150 { 35 151 FFRandom gen; 36 if (degree (Extension) < 0)37 factoryError("libfac: evaluate: Extension not inFF() or inGF() !"); 152 /*if (degree (Extension) < 0) 153 factoryError("libfac: evaluate: Extension not inFF() or inGF() !");*/ 38 154 return find_irreducible( degree_of_Extension, gen, Variable(1) ); 39 155 } … … 64 180 { 65 181 current++; 182 } 183 184 int hasAlgVar(const CanonicalForm &f, const Variable &v) 185 { 186 if (f.inBaseDomain()) return 0; 187 if (f.inCoeffDomain()) 188 { 189 if (f.mvar()==v) return 1; 190 return hasAlgVar(f.LC(),v); 191 } 192 if (f.inPolyDomain()) 193 { 194 if (hasAlgVar(f.LC(),v)) return 1; 195 for( CFIterator i=f; i.hasTerms(); i++) 196 { 197 if (hasAlgVar(i.coeff(),v)) return 1; 198 } 199 } 200 return 0; 201 } 202 203 int hasVar(const CanonicalForm &f, const Variable &v) 204 { 205 if (f.inBaseDomain()) return 0; 206 if (f.inCoeffDomain()) 207 { 208 if (f.mvar()==v) return 1; 209 return hasAlgVar(f.LC(),v); 210 } 211 if (f.inPolyDomain()) 212 { 213 if (f.mvar()==v) return 1; 214 if (hasVar(f.LC(),v)) return 1; 215 for( CFIterator i=f; i.hasTerms(); i++) 216 { 217 if (hasVar(i.coeff(),v)) return 1; 218 } 219 } 220 return 0; 221 } 222 223 int hasAlgVar(const CanonicalForm &f) 224 { 225 if (f.inBaseDomain()) return 0; 226 if (f.inCoeffDomain()) 227 { 228 if (f.level()!=0) 229 return 1; 230 return hasAlgVar(f.LC()); 231 } 232 if (f.inPolyDomain()) 233 { 234 if (hasAlgVar(f.LC())) return 1; 235 for( CFIterator i=f; i.hasTerms(); i++) 236 { 237 if (hasAlgVar(i.coeff())) return 1; 238 } 239 } 240 return 0; 241 } 242 243 /// pseudo division of f and g wrt. x s.t. multiplier*f=q*g+r 244 void 245 psqr (const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, 246 CanonicalForm & r, CanonicalForm& multiplier, const Variable& x) 247 { 248 ASSERT( x.level() > 0, "type error: polynomial variable expected" ); 249 ASSERT( ! g.isZero(), "math error: division by zero" ); 250 251 // swap variables such that x's level is larger or equal 252 // than both f's and g's levels. 253 Variable X; 254 if (f.level() > g.level()) 255 X= f.mvar(); 256 else 257 X= g.mvar(); 258 if (X.level() < x.level()) 259 X= x; 260 CanonicalForm F= swapvar (f, x, X); 261 CanonicalForm G= swapvar (g, x, X); 262 263 // now, we have to calculate the pseudo remainder of F and G 264 // w.r.t. X 265 int fDegree= degree (F, X); 266 int gDegree= degree (G, X); 267 if (fDegree < 0 || fDegree < gDegree) 268 { 269 q= 0; 270 r= f; 271 } 272 else 273 { 274 CanonicalForm LCG= LC (G, X); 275 multiplier= power (LCG, fDegree - gDegree + 1); 276 divrem (multiplier*F, G, q, r); 277 q= swapvar (q, x, X); 278 r= swapvar (r, x, X); 279 } 280 } 281 282 static CanonicalForm 283 Sprem ( const CanonicalForm &f, const CanonicalForm &g, CanonicalForm & m, CanonicalForm & q ) 284 { 285 CanonicalForm ff, gg, l, test, retvalue; 286 int df, dg,n; 287 bool reord; 288 Variable vf, vg, v; 289 290 if ( (vf = f.mvar()) < (vg = g.mvar()) ) 291 { 292 m=CanonicalForm(0); q=CanonicalForm(0); 293 return f; 294 } 295 else 296 { 297 if ( vf == vg ) 298 { 299 ff = f; gg = g; 300 reord = false; 301 v = vg; // == x 302 } 303 else 304 { 305 v = Variable(level(f.mvar()) + 1); 306 ff = swapvar(f,vg,v); // == r 307 gg = swapvar(g,vg,v); // == v 308 reord=true; 309 } 310 dg = degree( gg, v ); // == dv 311 df = degree( ff, v ); // == dr 312 if (dg <= df) {l=LC(gg); gg = gg -LC(gg)*power(v,dg);} 313 else { l = 1; } 314 n= 0; 315 while ( ( dg <= df ) && ( !ff.isZero()) ) 316 { 317 test= power(v,df-dg) * gg * LC(ff); 318 if ( df == 0 ){ff= ff.genZero();} 319 else {ff= ff - LC(ff)*power(v,df);} 320 ff = l*ff-test; 321 df= degree(ff,v); 322 n++; 323 } 324 if ( reord ) 325 { 326 retvalue= swapvar( ff, vg, v ); 327 } 328 else 329 { 330 retvalue= ff; 331 } 332 m= power(l,n); 333 if ( fdivides(g,m*f-retvalue) ) 334 q= (m*f-retvalue)/g; 335 else 336 q= CanonicalForm(0); 337 return retvalue; 338 } 339 } 340 341 CanonicalForm 342 divide( const CanonicalForm & ff, const CanonicalForm & f, const CFList & as) 343 { 344 CanonicalForm r,m,q; 345 346 if (f.inCoeffDomain()) 347 { 348 bool isRat=isOn(SW_RATIONAL); 349 if (getCharacteristic() == 0) 350 On(SW_RATIONAL); 351 q=ff/f; 352 if (!isRat && getCharacteristic() == 0) 353 Off(SW_RATIONAL); 354 } 355 else 356 r= Sprem(ff,f,m,q); 357 358 r= Prem(q,as); 359 return r; 360 } 361 362 CanonicalForm 363 alg_content (const CanonicalForm& f, const CFList& as) 364 { 365 if (!f.inCoeffDomain()) 366 { 367 CFIterator i= f; 368 CanonicalForm result= abs (i.coeff()); 369 i++; 370 while (i.hasTerms() && !result.isOne()) 371 { 372 result= alg_gcd (i.coeff(), result, as); 373 i++; 374 } 375 return result; 376 } 377 378 return abs (f); 379 } 380 381 CanonicalForm alg_gcd(const CanonicalForm & fff, const CanonicalForm &ggg, 382 const CFList &as) 383 { 384 if (fff.inCoeffDomain() || ggg.inCoeffDomain()) 385 return 1; 386 CanonicalForm f=fff; 387 CanonicalForm g=ggg; 388 f=Prem(f,as); 389 g=Prem(g,as); 390 if ( f.isZero() ) 391 { 392 if ( g.lc().sign() < 0 ) return -g; 393 else return g; 394 } 395 else if ( g.isZero() ) 396 { 397 if ( f.lc().sign() < 0 ) return -f; 398 else return f; 399 } 400 401 int v= as.getLast().level(); 402 if (f.level() <= v || g.level() <= v) 403 return 1; 404 405 CanonicalForm res; 406 407 // does as appear in f and g ? 408 bool has_alg_var=false; 409 for ( CFListIterator j=as;j.hasItem(); j++ ) 410 { 411 Variable v=j.getItem().mvar(); 412 if (hasVar (f, v)) 413 has_alg_var=true; 414 if (hasVar (g, v)) 415 has_alg_var=true; 416 } 417 if (!has_alg_var) 418 { 419 /*if ((hasAlgVar(f)) 420 || (hasAlgVar(g))) 421 { 422 Varlist ord; 423 for ( CFListIterator j=as;j.hasItem(); j++ ) 424 ord.append(j.getItem().mvar()); 425 res=algcd(f,g,as,ord); 426 } 427 else*/ 428 if (!hasAlgVar (f) && !hasAlgVar (g)) 429 res=gcd(f,g); 430 431 return res; 432 } 433 434 int mvf=f.level(); 435 int mvg=g.level(); 436 if (mvg > mvf) 437 { 438 CanonicalForm tmp=f; f=g; g=tmp; 439 int tmp2=mvf; mvf=mvg; mvg=tmp2; 440 } 441 if (g.inBaseDomain() || f.inBaseDomain()) 442 return CanonicalForm(1); 443 444 CanonicalForm c_f= alg_content (f, as); 445 446 if (mvf != mvg) 447 { 448 res= alg_gcd (g, c_f, as); 449 return res; 450 } 451 Variable x= f.mvar(); 452 453 // now: mvf==mvg, f.level()==g.level() 454 CanonicalForm c_g= alg_content (g, as); 455 456 int delta= degree (f) - degree (g); 457 458 f= divide (f, c_f, as); 459 g= divide (g, c_g, as); 460 461 // gcd of contents 462 CanonicalForm c_gcd= alg_gcd (c_f, c_g, as); 463 CanonicalForm tmp; 464 465 if (delta < 0) 466 { 467 tmp= f; 468 f= g; 469 g= tmp; 470 delta= -delta; 471 } 472 473 CanonicalForm r= 1; 474 475 while (degree (g, x) > 0) 476 { 477 r= Prem (f, g, as); 478 r= Prem (r, as); 479 if (!r.isZero()) 480 { 481 r= divide (r, alg_content (r,as), as); 482 r /= vcontent (r,Variable (v+1)); 483 } 484 f= g; 485 g= r; 486 } 487 488 if (degree (g, x) == 0) 489 return c_gcd; 490 491 c_f= alg_content (f, as); 492 493 f= divide (f, c_f, as); 494 495 f *= c_gcd; 496 f /= vcontent (f, Variable (v+1)); 497 498 return f; 66 499 } 67 500 … … 112 545 CFFListIterator i; 113 546 114 DEBOUTLN(CERR, "sqrf_norm_sub: f= ", f);115 DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);116 547 myrandom.reset(); s=myrandom.item(); g=f; 117 548 R= CanonicalForm(0); 118 DEBOUTLN(CERR, "sqrf_norm_sub: myrandom s= ", s);119 549 120 550 // Norm, resultante taken with respect to y 121 551 while ( !sqfreetest ) 122 552 { 123 DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);124 553 R = resultante(Palpha, g, y); R= R* bCommonDen(R); 125 DEBOUTLN(CERR, "sqrf_norm_sub: R= ", R);126 554 R /= content (R); 127 555 // sqfree check ; R is a polynomial in K[x] … … 129 557 { 130 558 temp= gcd(R, R.deriv(vf)); 131 DEBOUTLN(CERR, "sqrf_norm_sub: temp= ", temp);132 559 if (degree(temp,vf) != 0 || temp == temp.genZero() ) 133 560 sqfreetest= 0; 134 561 else 135 562 sqfreetest= 1; 136 DEBOUTLN(CERR, "sqrf_norm_sub: sqfreetest= ", sqfreetest);137 563 } 138 564 else 139 565 { 140 DEBOUTMSG(CERR, "Starting SqrFreeTest(R)!");141 // Look at SqrFreeTest!142 // (z+a^5+w)^4 with z<w<a should not give sqfreetest=1 !143 // for now we use this workaround with Factorize...144 // ...but it should go away soon!!!!145 566 Variable X; 146 567 if (hasFirstAlgVar(R,X)) … … 148 569 else 149 570 testlist= factorize(R); 150 DEBOUTLN(CERR, "testlist= ", testlist); 571 151 572 if (testlist.getFirst().factor().inCoeffDomain()) 152 573 testlist.removeFirst(); … … 160 581 } 161 582 } 162 DEBOUTLN(CERR, "SqrFreeTest(R)= ", sqfreetest);163 583 } 164 584 if ( ! sqfreetest ) 165 585 { 166 586 myrandom.next(); 167 DEBOUTLN(CERR, "sqrf_norm_sub generated new myrandom item: ", myrandom.item());168 587 if ( getCharacteristic() == 0 ) 169 588 t= CanonicalForm(mapinto(myrandom.item())); … … 171 590 t= CanonicalForm(myrandom.item()); 172 591 s= t; 173 DEBOUTLN(CERR, "sqrf_norm_sub: testing s= ", s);174 592 g= f(f.mvar()-t*Palpha.mvar(), f.mvar()); 175 DEBOUTLN(CERR, " gives g= ", g);176 593 } 177 594 } … … 188 605 CFFListIterator i; 189 606 190 DEBOUTLN(CERR, "sqrf_norm_sub: f= ", f);191 DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);192 607 myrandom.reset(); s=myrandom.item(); g=f; 193 608 R= CanonicalForm(0); 194 DEBOUTLN(CERR, "sqrf_norm_sub: myrandom s= ", s);195 609 196 610 // Norm, resultante taken with respect to y 197 611 while ( !sqfreetest ) 198 612 { 199 DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);200 613 R = resultante(Palpha, g, y); R= R* bCommonDen(R); 201 DEBOUTLN(CERR, "sqrf_norm_sub: R= ", R);202 614 R /= content (R); 203 615 // sqfree check ; R is a polynomial in K[x] … … 205 617 { 206 618 temp= gcd(R, R.deriv(vf)); 207 DEBOUTLN(CERR, "sqrf_norm_sub: temp= ", temp);208 619 if (degree(temp,vf) != 0 || temp == temp.genZero() ) 209 620 sqfreetest= 0; 210 621 else 211 622 sqfreetest= 1; 212 DEBOUTLN(CERR, "sqrf_norm_sub: sqfreetest= ", sqfreetest);213 623 } 214 624 else 215 625 { 216 DEBOUTMSG(CERR, "Starting SqrFreeTest(R)!");217 // Look at SqrFreeTest!218 // (z+a^5+w)^4 with z<w<a should not give sqfreetest=1 !219 // for now we use this workaround with Factorize...220 // ...but it should go away soon!!!!221 626 Variable X; 222 627 if (hasFirstAlgVar(R,X)) … … 224 629 else 225 630 testlist= factorize(R); 226 DEBOUTLN(CERR, "testlist= ", testlist);227 631 if (testlist.getFirst().factor().inCoeffDomain()) 228 632 testlist.removeFirst(); … … 236 640 } 237 641 } 238 DEBOUTLN(CERR, "SqrFreeTest(R)= ", sqfreetest);239 642 } 240 643 if ( ! sqfreetest ) 241 644 { 242 645 myrandom.next(); 243 DEBOUTLN(CERR, "sqrf_norm_sub generated new myrandom item: ", myrandom.item());244 646 if ( getCharacteristic() == 0 ) 245 647 t= CanonicalForm(mapinto(myrandom.item())); … … 247 649 t= CanonicalForm(myrandom.item()); 248 650 s= t; 249 DEBOUTLN(CERR, "sqrf_norm_sub: testing s= ", s);250 651 g= f(f.mvar()-t*Palpha.mvar(), f.mvar()); 251 DEBOUTLN(CERR, " gives g= ", g);252 652 } 253 653 } … … 259 659 CanonicalForm & R) 260 660 { 261 DEBOUTLN(CERR, "sqrf_norm: f= ", f);262 DEBOUTLN(CERR, "sqrf_norm: Palpha= ", PPalpha);263 661 if ( getCharacteristic() == 0 ) 264 662 { 265 663 IntGenerator myrandom; 266 DEBOUTMSG(CERR, "sqrf_norm: no extension, char=0");267 664 sqrf_norm_sub(f,PPalpha, myrandom, s,g,R); 268 DEBOUTLN(CERR, "sqrf_norm: f= ", f);269 DEBOUTLN(CERR, "sqrf_norm: Palpha= ", PPalpha);270 DEBOUTLN(CERR, "sqrf_norm: s= ", s);271 DEBOUTLN(CERR, "sqrf_norm: g= ", g);272 DEBOUTLN(CERR, "sqrf_norm: R= ", R);273 665 } 274 666 else if ( degree(Extension) > 0 ) // working over Extensions 275 667 { 276 DEBOUTLN(CERR, "sqrf_norm: degree of extension is ", degree(Extension));277 668 AlgExtGenerator myrandom(Extension); 278 669 sqrf_agnorm_sub(f,PPalpha, myrandom, s,g,R); … … 281 672 { 282 673 FFGenerator myrandom; 283 DEBOUTMSG(CERR, "sqrf_norm: degree of extension is 0");284 674 sqrf_norm_sub(f,PPalpha, myrandom, s,g,R); 285 675 } … … 347 737 l= l+1; 348 738 if ( igcd(k,i.getItem()) == 1){ 349 DEBOUTLN(CERR, "getextension: gcd == 1, l=",l);350 739 if ( l==length ){ setCharacteristic(charac); return k; } 351 740 } 352 else { DEBOUTMSG(CERR, "getextension: Next iteration");break; }741 else { break; } 353 742 } 354 743 k= k+1; l=0; 355 744 } 356 745 while ( 1 ); 357 }358 359 360 /// pseudo division of f and g wrt. x s.t. multiplier*f=q*g+r361 /// but only if the leading coefficient of g is of level lower than coeffLevel362 void363 psqr (const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q,364 CanonicalForm & r, CanonicalForm& multiplier, const Variable& x,365 int coeffLevel)366 {367 ASSERT( x.level() > 0, "type error: polynomial variable expected" );368 ASSERT( ! g.isZero(), "math error: division by zero" );369 370 // swap variables such that x's level is larger or equal371 // than both f's and g's levels.372 Variable X;373 if (f.level() > g.level())374 X= f.mvar();375 else376 X= g.mvar();377 if (X.level() < x.level())378 X= x;379 CanonicalForm F= swapvar (f, x, X);380 CanonicalForm G= swapvar (g, x, X);381 382 // now, we have to calculate the pseudo remainder of F and G383 // w.r.t. X384 int fDegree= degree (F, X);385 int gDegree= degree (G, X);386 if (fDegree < 0 || fDegree < gDegree)387 {388 q= 0;389 r= f;390 }391 else392 {393 CanonicalForm LCG= LC (G, X);394 if (LCG.level() < coeffLevel)395 {396 multiplier= power (LCG, fDegree - gDegree + 1);397 divrem (multiplier*F, G, q, r);398 q= swapvar (q, x, X);399 r= swapvar (r, x, X);400 }401 else402 {403 q= 0;404 r= f;405 }406 }407 }408 409 /// pseudo division of f and g wrt. x s.t. multiplier*f=q*g+r410 void411 psqr (const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q,412 CanonicalForm & r, CanonicalForm& multiplier, const Variable& x)413 {414 ASSERT( x.level() > 0, "type error: polynomial variable expected" );415 ASSERT( ! g.isZero(), "math error: division by zero" );416 417 // swap variables such that x's level is larger or equal418 // than both f's and g's levels.419 Variable X;420 if (f.level() > g.level())421 X= f.mvar();422 else423 X= g.mvar();424 if (X.level() < x.level())425 X= x;426 CanonicalForm F= swapvar (f, x, X);427 CanonicalForm G= swapvar (g, x, X);428 429 // now, we have to calculate the pseudo remainder of F and G430 // w.r.t. X431 int fDegree= degree (F, X);432 int gDegree= degree (G, X);433 if (fDegree < 0 || fDegree < gDegree)434 {435 q= 0;436 r= f;437 }438 else439 {440 CanonicalForm LCG= LC (G, X);441 multiplier= power (LCG, fDegree - gDegree + 1);442 divrem (multiplier*F, G, q, r);443 q= swapvar (q, x, X);444 r= swapvar (r, x, X);445 }446 746 } 447 747 … … 568 868 bool isRat= isOn (SW_RATIONAL); 569 869 570 DEBOUTLN(CERR, "simpleextension: Astar= ", Astar);571 DEBOUTLN(CERR, "simpleextension: R= ", R);572 DEBOUTLN(CERR, "simpleextension: Extension= ", Extension);573 870 CFListIterator j; 574 871 if (Astar.length() == 1) … … 654 951 } 655 952 656 DEBOUTLN(CERR, "simpleextension: g= ", g);657 DEBOUTLN(CERR, "simpleextension: s= ", s);658 DEBOUTLN(CERR, "simpleextension: R= ", R);659 953 Returnlist.append (ra); 660 954 if (isFunctionField) … … 680 974 return alg_lc(f.LC()); 681 975 } 682 //assert(f.inCoeffDomain()); 976 683 977 return f; 684 978 } … … 851 1145 852 1146 // the heart of the algorithm: the one from Trager 853 #ifndef DEBUGOUTPUT854 static CFFList855 alg_factor( const CanonicalForm & F, const CFList & Astar, const Variable & vminpoly, const Varlist /*& oldord*/, const CFList & as, bool isFunctionField)856 #else857 1147 static CFFList 858 1148 alg_factor( const CanonicalForm & F, const CFList & Astar, const Variable & vminpoly, const Varlist & oldord, const CFList & as, bool isFunctionField) 859 #endif860 1149 { 861 1150 bool isRat= isOn (SW_RATIONAL); … … 864 1153 CFList substlist, backSubsts; 865 1154 866 DEBINCLEVEL(CERR,"alg_factor");867 DEBOUTLN(CERR, "alg_factor: f= ", f);868 869 1155 substlist= simpleextension (backSubsts, Astar, vminpoly, isFunctionField, Rstar); 870 DEBOUTLN(CERR, "alg_factor: substlist= ", substlist);871 DEBOUTLN(CERR, "alg_factor: minpoly Rstar= ", Rstar);872 DEBOUTLN(CERR, "alg_factor: vminpoly= ", vminpoly);873 1156 874 1157 f= subst (f, Astar, substlist, Rstar, isFunctionField); … … 902 1185 // after here we are over an extension of a function field 903 1186 904 905 1187 // make quasi monic 906 1188 CFList Rstarlist= CFList (Rstar); … … 927 1209 928 1210 sqrf_norm(f, Rstar, vminpoly, s, g, R ); 929 //out_cf("sqrf_norm R:",R,"\n"); 930 //out_cf("sqrf_norm s:",s,"\n"); 931 //out_cf("sqrf_norm g:",g,"\n"); 932 DEBOUTLN(CERR, "alg_factor: g= ", g); 933 DEBOUTLN(CERR, "alg_factor: s= ", s); 934 DEBOUTLN(CERR, "alg_factor: R= ", R); 1211 935 1212 Variable X; 936 1213 if (hasFirstAlgVar(R,X)) 937 1214 { 938 1215 // factorize R over alg.extension with X 939 //CERR << "alg: "<< X << " mipo=" << getMipo(X,Variable('X')) <<"\n";940 DEBOUTLN(CERR, "alg_factor: factorize( ", R);941 1216 Factorlist = factorize( R, X ); 942 1217 } … … 944 1219 { 945 1220 // factor R over k 946 DEBOUTLN(CERR, "alg_factor: Factorize( ", R);947 1221 Factorlist = factorize(R); 948 1222 } 949 1223 950 DEBOUTLN(CERR, "alg_factor: Factorize(R)= ", Factorlist);951 1224 if ( !Factorlist.getFirst().factor().inCoeffDomain() ) 952 1225 Factorlist.insert(CFFactor(1,1)); … … 963 1236 { 964 1237 g= f; 965 DEBOUTLN(CERR, "alg_factor: g= ", g);966 1238 for ( CFFListIterator i=Factorlist; i.hasItem(); i++) 967 1239 { … … 974 1246 fnew= reduce (fnew, Rstar); 975 1247 } 976 977 DEBOUTLN(CERR, "alg_factor: fnew= ", fnew);978 1248 979 1249 h= alg_gcd (g, fnew, Rstarlist); … … 1018 1288 LL=L; 1019 1289 } 1020 //CFFListIterator i=LL; 1021 //for(;i.hasItem(); i++ ) 1022 // out_cf("end alg_f:",i.getItem().factor(),"\n"); 1023 //printf("end alg_factor\n"); 1024 DEBOUTLN(CERR, "alg_factor: L= ", LL); 1025 DEBDECLEVEL(CERR,"alg_factor"); 1290 1026 1291 if (!isRat && getCharacteristic() == 0) 1027 1292 Off (SW_RATIONAL); … … 1225 1490 1226 1491 // factor F with minimal polys given in asnew 1227 PremForm rem;1228 1492 asnew.append (F); 1229 asnew= mcharset (asnew, rem); // TODO use modCharSet1493 asnew= charSetViaModCharSet (asnew); // TODO use modCharSet 1230 1494 1231 1495 F= asnew.getLast(); … … 1269 1533 for (CFFListIterator k= tmp; k.hasItem(); k++) 1270 1534 { 1271 rem= PremForm();1272 1535 transform= transBack; 1273 1536 CanonicalForm factor= k.getItem().factor(); 1274 1537 factor= M (factor); 1275 1538 transform.append (factor); 1276 transform= mcharset (transform, rem); //TODO use modCharSet1539 transform= charSetViaModCharSet (transform); //TODO use modCharSet 1277 1540 for (i= transform; i.hasItem(); i++) 1278 1541 { … … 1341 1604 1342 1605 success=1; 1343 DEBINCLEVEL(CERR, "newfactoras");1344 DEBOUTLN(CERR, "newfactoras called with f= ", f);1345 DEBOUTLN(CERR, " content(f)= ", content(f));1346 DEBOUTLN(CERR, " as= ", as);1347 DEBOUTLN(CERR, "newfactoras: cls(vf)= ", cls(vf));1348 DEBOUTLN(CERR, "newfactoras: cls(as.getLast())= ", cls(as.getLast()));1349 DEBOUTLN(CERR, "newfactoras: degree(f,vf)= ", degree(f,vf));1350 1606 1351 1607 // F1: [Test trivial cases] 1352 1608 // 1) first trivial cases: 1353 if ( cls(vf) <= cls(as.getLast()))1609 if (vf.level() <= as.getLast().level()) 1354 1610 { 1355 1611 // ||( (as.length()==1) && (degree(f,vf)==3) && (degree(as.getFirst()==2)) ) 1356 DEBDECLEVEL(CERR,"newfactoras");1357 1612 if (!isRat && getCharacteristic() == 0) 1358 1613 Off (SW_RATIONAL); … … 1379 1634 } 1380 1635 uord= Difference(uord,ord); 1381 DEBOUTLN(CERR, "Astar is: ", Astar);1382 DEBOUTLN(CERR, "ord is: ", ord);1383 DEBOUTLN(CERR, "uord is: ", uord);1384 1636 1385 1637 // 3) second trivial cases: we already proved irr. of f over no extensions 1386 1638 if ( Astar.length() == 0 ){ 1387 DEBDECLEVEL(CERR,"newfactoras");1388 1639 if (!isRat && getCharacteristic() == 0) 1389 1640 Off (SW_RATIONAL); … … 1399 1650 // polynomials, we don't have transzendentals. 1400 1651 Varlist newuord=Var_is_in_AS(uord,Astar); 1401 DEBOUTLN(CERR, "newuord is: ", newuord);1402 1652 1403 1653 CFFList Factorlist; … … 1410 1660 Fgcd= alg_gcd(f,f.deriv(),Astar); 1411 1661 1412 if ( Fgcd == 0 ) {DEBOUTMSG(CERR, "WARNING: p'th root ?");}1413 1414 1662 bool derivZero= f.deriv().isZero(); 1415 1663 if (isFunctionField && (degree (Fgcd, f.mvar()) > 0) && !derivZero) 1416 1664 { 1417 DEBOUTLN(CERR, "Nontrivial GCD found of ", f);1418 1665 CanonicalForm Ggcd= divide(f, Fgcd,Astar); 1419 1666 if (getCharacteristic() == 0) … … 1425 1672 return result; 1426 1673 } 1427 DEBOUTLN(CERR, " split into ", Fgcd); 1428 DEBOUTLN(CERR, " and ", Ggcd); 1674 1429 1675 Fgcd= pp(Fgcd); Ggcd= pp(Ggcd); 1430 DEBDECLEVEL(CERR,"newfactoras");1431 1676 if (!isRat && getCharacteristic() == 0) 1432 1677 Off (SW_RATIONAL); … … 1441 1686 for (i=Astar; i.hasItem(); i++){degreelist.append(degree(i.getItem()));} 1442 1687 int extdeg= getextension(degreelist, degree(f)); 1443 DEBOUTLN(CERR, "Extension needed of degree ", extdeg);1444 1688 1445 1689 // Now the real stuff! 1446 1690 if ( newuord.length() == 0 ) // no transzendentals 1447 1691 { 1448 DEBOUTMSG(CERR, "No transzendentals!");1449 1692 if ( extdeg > 1 ){ 1450 1693 CanonicalForm MIPO= generate_mipo( extdeg, vminpoly); 1451 DEBOUTLN(CERR, "Minpoly produced ", MIPO);1452 1694 vminpoly= rootOf(MIPO); 1453 1695 } 1454 1696 Factorlist= alg_factor(f, Astar, vminpoly, oldord, as, isFunctionField); 1455 DEBDECLEVEL(CERR,"newfactoras");1456 1697 return Factorlist; 1457 1698 } 1458 1699 else if ( inseperable(Astar) > 0 || derivZero) // Look if extensions are seperable 1459 1700 { 1460 DEBOUTMSG(CERR, "Inseperable extensions! Using Endler!");1461 1701 Factorlist= SteelTrager(f, Astar, newuord); 1462 DEBOUTLN(CERR, "Endler gives: ", Factorlist);1463 1702 return Factorlist; 1464 1703 } 1465 1704 else{ // we are on the save side: Use trager 1466 DEBOUTMSG(CERR, "Only seperable extensions!");1467 1705 if (extdeg > 1 ){ 1468 1706 CanonicalForm MIPO=generate_mipo(extdeg, vminpoly ); 1469 1707 vminpoly= rootOf(MIPO); 1470 DEBOUTLN(CERR, "Minpoly generated: ", MIPO);1471 DEBOUTLN(CERR, "vminpoly= ", vminpoly);1472 DEBOUTLN(CERR, "degree(vminpoly)= ", degree(vminpoly));1473 1708 } 1474 1709 Factorlist= alg_factor(f, Astar, vminpoly, oldord, as, isFunctionField); 1475 DEBDECLEVEL(CERR,"newfactoras");1476 1710 return Factorlist; 1477 1711 } 1478 1712 } 1479 1713 else{ // char=0 apply trager directly 1480 DEBOUTMSG(CERR, "Char=0! Apply Trager!");1481 1714 Variable vminpoly; 1482 1715 Factorlist= alg_factor(f, Astar, vminpoly, oldord, as, isFunctionField); 1483 DEBDECLEVEL(CERR,"newfactoras");1484 1716 if (!isRat && getCharacteristic() == 0) 1485 1717 Off (SW_RATIONAL); … … 1487 1719 } 1488 1720 1489 DEBDECLEVEL(CERR,"newfactoras");1490 1721 return CFFList(CFFactor(f,1)); 1491 1722 } … … 1508 1739 return Factors; 1509 1740 } 1510 if ( cls(f) <= cls(as.getLast()))1741 if (f.level() <= as.getLast().level()) 1511 1742 { 1512 1743 if (!isRat && getCharacteristic() == 0) -
factory/facAlgFunc.h
r248194 r5355082 4 4 //////////////////////////////////////////////////////////// 5 5 6 #ifndef INCL_NEW_ALGFACTOR_H7 #define INCL_NEW_ALGFACTOR_H6 #ifndef FAC_ALG_FUNC_H 7 #define FAC_ALG_FUNC_H 8 8 9 #include <factory.h> 10 #include <tmpl_inst.h> // for typedef's 9 #include "canonicalform.h" 11 10 12 11 // missing class: IntGenerator: … … 26 25 }; 27 26 27 CanonicalForm alg_gcd(const CanonicalForm &, const CanonicalForm &, const CFList &); 28 /*BEGINPUBLIC*/ 28 29 CFFList newfactoras( const CanonicalForm & f, const CFList & as, int &success); 29 /*BEGINPUBLIC*/30 30 CFFList newcfactor(const CanonicalForm & f, const CFList & as, int & success ); 31 31 /*ENDPUBLIC*/ -
factory/libfac/charset/csutil.cc
r248194 r5355082 11 11 extern void out_cf(const char *s1,const CanonicalForm &f,const char *s2); 12 12 extern CanonicalForm alg_lc(const CanonicalForm &f); 13 14 int hasAlgVar(const CanonicalForm &f, const Variable &v)15 {16 if (f.inBaseDomain()) return 0;17 if (f.inCoeffDomain())18 {19 if (f.mvar()==v) return 1;20 return hasAlgVar(f.LC(),v);21 }22 if (f.inPolyDomain())23 {24 if (hasAlgVar(f.LC(),v)) return 1;25 for( CFIterator i=f; i.hasTerms(); i++)26 {27 if (hasAlgVar(i.coeff(),v)) return 1;28 }29 }30 return 0;31 }32 33 int hasVar(const CanonicalForm &f, const Variable &v)34 {35 if (f.inBaseDomain()) return 0;36 if (f.inCoeffDomain())37 {38 if (f.mvar()==v) return 1;39 return hasAlgVar(f.LC(),v);40 }41 if (f.inPolyDomain())42 {43 if (f.mvar()==v) return 1;44 if (hasVar(f.LC(),v)) return 1;45 for( CFIterator i=f; i.hasTerms(); i++)46 {47 if (hasVar(i.coeff(),v)) return 1;48 }49 }50 return 0;51 }52 53 int hasAlgVar(const CanonicalForm &f)54 {55 if (f.inBaseDomain()) return 0;56 if (f.inCoeffDomain())57 {58 if (f.level()!=0)59 {60 //CERR << "hasAlgVar:" << f.mvar() <<"\n";61 return 1;62 }63 return hasAlgVar(f.LC());64 }65 if (f.inPolyDomain())66 {67 if (hasAlgVar(f.LC())) return 1;68 for( CFIterator i=f; i.hasTerms(); i++)69 {70 if (hasAlgVar(i.coeff())) return 1;71 }72 }73 return 0;74 }75 76 13 77 14 static bool … … 278 215 } 279 216 280 CanonicalForm281 divide( const CanonicalForm & ff, const CanonicalForm & f, const CFList & as)282 {283 CanonicalForm r,m,q;284 285 //out_cf("divide f=",ff,"\n");286 //out_cf("divide g=",f,"\n");287 if (f.inCoeffDomain())288 {289 bool isRat=isOn(SW_RATIONAL);290 if (getCharacteristic() == 0)291 On(SW_RATIONAL);292 q=ff/f;293 if (!isRat && getCharacteristic() == 0)294 Off(SW_RATIONAL);295 }296 else297 r= Sprem(ff,f,m,q); //result in q, ignore r,m298 //CERR << "r= " << r << " , m= " << m << " , q= " << q << "\n";299 r= Prem(q,as);300 //CERR << "r= " << r << "\n";301 //out_cf(" ->",r,"\n");302 return r;303 }304 305 217 // This function allows as to be empty; in that case, it is equivalent 306 218 // to the previous version (treating no variables as algebraic). … … 870 782 #endif 871 783 872 CanonicalForm 873 alg_content (const CanonicalForm& f, const CFList& as) 874 { 875 if (!f.inCoeffDomain()) 876 { 877 CFIterator i= f; 878 CanonicalForm result= abs (i.coeff()); 879 i++; 880 while (i.hasTerms() && !result.isOne()) 881 { 882 result= alg_gcd (i.coeff(), result, as); 883 i++; 884 } 885 return result; 886 } 887 888 return abs (f); 889 } 890 891 CanonicalForm alg_gcd(const CanonicalForm & fff, const CanonicalForm &ggg, 892 const CFList &as) 893 { 894 if (fff.inCoeffDomain() || ggg.inCoeffDomain()) 895 return 1; 896 CanonicalForm f=fff; 897 CanonicalForm g=ggg; 898 f=Prem(f,as); 899 g=Prem(g,as); 900 if ( f.isZero() ) 901 { 902 if ( g.lc().sign() < 0 ) return -g; 903 else return g; 904 } 905 else if ( g.isZero() ) 906 { 907 if ( f.lc().sign() < 0 ) return -f; 908 else return f; 909 } 910 //out_cf("alg_gcd fff(",fff," \n "); 911 //out_cf("ggg",ggg,")\n"); 912 int v= as.getLast().level(); 913 if (f.level() <= v || g.level() <= v) 914 return 1; 915 916 CanonicalForm res; 917 // does as appear in f and g ? 918 bool has_alg_var=false; 919 for ( CFListIterator j=as;j.hasItem(); j++ ) 920 { 921 Variable v=j.getItem().mvar(); 922 if (hasVar(f,v)) {has_alg_var=true; /*break;*/} 923 if (hasVar(g,v)) {has_alg_var=true; /*break;*/} 924 //out_cf("as:",j.getItem(),"\n"); 925 } 926 if (!has_alg_var) 927 { 928 if ((hasAlgVar(f)) 929 || (hasAlgVar(g))) 930 { 931 Varlist ord; 932 for ( CFListIterator j=as;j.hasItem(); j++ ) 933 ord.append(j.getItem().mvar()); 934 res=algcd(f,g,as,ord); 935 } 936 else 937 res=gcd(f,g); 938 //out_cf("gcd=",res,"\n"); 939 //out_cf("of f=",fff," , "); 940 //out_cf("and g=",ggg,"\n"); 941 942 return res; 943 } 944 945 int mvf=f.level(); 946 int mvg=g.level(); 947 if (mvg > mvf) 948 { 949 CanonicalForm tmp=f; f=g; g=tmp; 950 int tmp2=mvf; mvf=mvg; mvg=tmp2; 951 } 952 if (g.inBaseDomain() || f.inBaseDomain()) 953 { 954 //printf("const\n"); 955 //out_cf("of f=",fff," , "); 956 //out_cf("and g=",ggg,"\n"); 957 return CanonicalForm(1); 958 } 959 960 // gcd of all coefficients: 961 CanonicalForm c_f= alg_content (f, as); 962 963 //printf("f.mvar=%d (%d), g.mvar=%d (%d)\n",f.level(),mvf,g.level(),mvg); 964 if (mvf != mvg) // => mvf > mvg 965 { 966 res= alg_gcd (g, c_f, as); 967 return res; 968 } 969 Variable x= f.mvar(); 970 // now: mvf==mvg, f.level()==g.level() 971 // content of g 972 CanonicalForm c_g= alg_content (g, as); 973 974 int delta= degree (f) - degree (g); 975 976 //f/=c_gcd; 977 //g/=c_gcd; 978 f= divide (f, c_f, as); 979 g= divide (g, c_g, as); 980 981 // gcd of contents 982 CanonicalForm c_gcd= alg_gcd (c_f, c_g, as); 983 CanonicalForm tmp; 984 985 if (delta < 0) 986 { 987 tmp= f; 988 f= g; 989 g= tmp; 990 delta= -delta; 991 } 992 993 CanonicalForm r=1; 994 995 while (degree (g, x) > 0) 996 { 997 r= Prem (f, g, as); 998 r= Prem (r, as); 999 if (!r.isZero()) 1000 { 1001 r= divide (r, alg_content (r,as), as); 1002 r /= vcontent (r,Variable (v+1)); 1003 } 1004 f= g; 1005 g= r; 1006 } 1007 1008 if (degree (g, x) == 0) 1009 return c_gcd; 1010 1011 c_f= alg_content (f, as); 1012 1013 f= divide (f, c_f, as); 1014 1015 f *= c_gcd; 1016 f /= vcontent (f, Variable (v+1)); 1017 1018 return f; 1019 } 1020 784 785 -
factory/libfac/charset/csutil.h
r248194 r5355082 52 52 CanonicalForm Prem( const CanonicalForm &f, const CFList &L ); 53 53 CFList Prem( const CFList &AS, const CFList &L ); 54 CanonicalForm alg_gcd(const CanonicalForm &, const CanonicalForm &, const CFList &);55 54 /*ENDPUBLIC*/ 56 CanonicalForm divide( const CanonicalForm & ff, const CanonicalForm & f, const CFList & as);57 55 CFList remsetb( const CFList & ps, const CFList & as); 58 56 CanonicalForm lowestRank( const CFList & F ); -
factory/libfac/factor/helpstuff.cc
r248194 r5355082 21 21 } 22 22 23 // Now some procedures used in SqrFree and in Factor24 ///////////////////////////////////////////////////////////////25 ///////////////////////////////////////////////////////////////26 // We have to include a version of <CFFList>.append(CFFactor)//27 // and Union( CFFList, CFFList) //28 // because we have to look for multiplicities in SqrFree. //29 // e.g.: SqrFree( f^3 ) with char <> 3 //30 ///////////////////////////////////////////////////////////////31 CFFList32 myappend( const CFFList & Inputlist, const CFFactor & TheFactor)33 {34 CFFList Outputlist ;35 CFFactor copy;36 CFFListIterator i;37 int exp=0;38 39 for ( i=Inputlist ; i.hasItem() ; i++ )40 {41 copy = i.getItem();42 if ( copy.factor() == TheFactor.factor() )43 exp += copy.exp();44 else45 Outputlist.append(copy);46 }47 Outputlist.append( CFFactor(TheFactor.factor(), exp + TheFactor.exp()));48 return Outputlist;49 }50 51 CFFList52 myUnion(const CFFList & Inputlist1,const CFFList & Inputlist2)53 {54 CFFList Outputlist;55 CFFListIterator i;56 57 for ( i=Inputlist1 ; i.hasItem() ; i++ )58 Outputlist = myappend(Outputlist, i.getItem() );59 for ( i=Inputlist2 ; i.hasItem() ; i++ )60 Outputlist = myappend(Outputlist, i.getItem() );61 62 return Outputlist;63 }64 -
factory/libfac/factor/helpstuff.h
r248194 r5355082 9 9 // Now some procedures used in SqrFree and in Factor 10 10 /////////////////////////////////////////////////////////////// 11 CFFList myappend( const CFFList & Inputlist, const CFFactor & TheFactor) ;12 CFFList myUnion(const CFFList & Inputlist1,const CFFList & Inputlist2);13 11 inline int min ( const int a, const int b ){ 14 12 return (a<=b ? a:b);
Note: See TracChangeset
for help on using the changeset viewer.