Changeset 8eb6da in git
- Timestamp:
- Nov 15, 2010, 10:29:11 AM (13 years ago)
- Branches:
- (u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
- Children:
- 7993adcf5a444486c5b07eacda51e7dfc6d67afc
- Parents:
- 8b2ef529c04c7b3517f5001ff6c69e582bb9fd09
- Location:
- libfac/charset
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
libfac/charset/alg_factor.cc
r8b2ef5 r8eb6da 117 117 118 118 return resultant(fz,gz,v); 119 120 CanonicalForm h, beta, help, F, G;121 int delta;122 123 DEBOUTLN( CERR, "resultante: called f= ", f);124 DEBOUTLN( CERR, "resultante: called g= ", g);125 DEBOUTLN( CERR, "resultante: called v= ", v);126 if ( f.mvar() < v || g.mvar() < v ){127 DEBOUTMSG(CERR, "resultante: f.mvar() < v || g.mvar() < v");128 return 1;129 }130 131 if ( f.degree( v ) < 1 || g.degree( v ) < 1 ){132 DEBOUTMSG(CERR, "resultante: f.degree( v ) < 1 || g.degree( v ) < 1");133 // If deg(F,v) == 0 , then resultante(F,G,v) = F^n, where n=deg(G,v)134 if ( f.degree( v ) < 1 ) return power(f,degree(g,v));135 else return power(g,degree(f,v));136 }137 138 if ( f.degree( v ) >= g.degree( v ) ) { F = f; G = g; }139 else { G = f; F = g; }140 141 h = CanonicalForm(1);142 while ( G != G.genZero() ) {143 delta= degree(F,v) -degree(G,v);144 beta = power(CanonicalForm(-1), delta+1) * LC(F,v)* power(h, delta);145 h= (h * power(LC(G,v), delta)) / power(h, delta);146 help= G;147 G= mypsr(F,G,v);148 G= G/beta;149 F=help;150 }151 if ( degree(F,v) != 0 ) F= CanonicalForm(0);152 return F;153 119 } 154 120 … … 530 496 if (degree(i.getItem().factor()) > 0 ) 531 497 { 532 // undo linear transformation!!!! and then gcd!533 //CERR << "algcd(" << g << "," << fnew << ",as" << as << ")" << "\n";534 //out_cf("g=",g,"\n");535 //out_cf("fnew=",fnew,"\n");536 //h= algcd(g,fnew, as, oldord);537 //if (as.length() >1)538 // h= algcd(g,fnew, as, oldord);539 //else540 498 h=alg_gcd(g,fnew,as); 541 //out_cf(" -> algcd=",algcd(g,fnew, as, oldord),"\n");542 //out_cf(" -> alg_gcd=",alg_gcd(g,fnew,as),"\n");543 //CERR << "algcd result:" << h << "\n";544 499 DEBOUTLN(CERR, " alg_factor: h= ", h); 545 500 DEBOUTLN(CERR, " alg_factor: oldord= ", oldord); 546 501 if ( degree(h) > 0 ) 547 502 { //otherwise it's a constant 548 //CanonicalForm c=LC(h,g.mvar());549 //out_cf("new lc h:",c,"\n");550 //h= divide(h,c,as);551 //out_cf("new factor h/c:",h,"\n");552 503 g= divide(g, h,as); 553 504 DEBOUTLN(CERR, "alg_factor: g/h= ", g); … … 736 687 Varlist gcdord= Union(ord,newuord); gcdord.append(f.mvar()); 737 688 // This is for now. we need alg_sqrfree implemented! 738 //CERR << "algcd(" << f << "," << f.deriv() << " as:" << Astar <<"\n";739 //CanonicalForm Fgcd= algcd(f,f.deriv(),Astar,gcdord);740 689 CanonicalForm Fgcd; 741 //if (Astar.length() >1)742 // Fgcd= algcd(f,f.deriv(),Astar,gcdord);743 //else744 690 Fgcd= alg_gcd(f,f.deriv(),Astar); 745 //out_cf("algcd:",algcd(f,f.deriv(),Astar,gcdord),"\n");746 //out_cf("alg_gcd:",alg_gcd(f,f.deriv(),Astar),"\n");747 // CERR << "algcd result:" << Fgcd << "\n";748 691 if ( Fgcd == 0 ) DEBOUTMSG(CERR, "WARNING: p'th root ?"); 749 692 if (( degree(Fgcd, f.mvar()) > 0) && (!(f.deriv().isZero())) ){ -
libfac/charset/algfactor.cc
r8b2ef5 r8eb6da 182 182 return result; 183 183 } 184 185 CFFList186 factoras( const CanonicalForm & f, const CFList & as, int & success ){187 Variable vf=f.mvar();188 CFListIterator i;189 CFFListIterator jj;190 CFList reduceresult;191 CFFList result;192 193 DEBINCLEVEL(CERR, "factoras");194 DEBOUTLN(CERR, "factoras called with f= ", f);195 DEBOUTLN(CERR, " content(f)= ", content(f));196 DEBOUTLN(CERR, " as= ", as);197 DEBOUTLN(CERR, "factoras: cls(vf)= ", cls(vf));198 DEBOUTLN(CERR, "factoras: cls(as.getLast())= ", cls(as.getLast()));199 DEBOUTLN(CERR, "factoras: degree(f,vf)= ", degree(f,vf));200 201 // F1: [Test trivial cases]202 // 1) first trivial cases:203 if ( (cls(vf) <= cls(as.getLast())) ||204 degree(f,vf)<=1205 // ||( (as.length()==1) && (degree(f,vf)==3) && (degree(as.getFirst()==2)) )206 ){207 success=1;208 DEBDECLEVEL(CERR,"factoras");209 return CFFList(CFFactor(f,1));210 }211 212 213 // 2) List of variables:214 // 2a) Setup list of those polys in AS having degree(AS[i], AS[i].mvar()) > 1215 // 2b) Setup variableordering216 CFList Astar;217 Variable x;218 CanonicalForm elem;219 Varlist ord, uord,oldord;220 for ( int ii=1; ii< level(vf) ; ii++ ) { uord.append(Variable(ii)); }221 oldord= uord; oldord.append(vf);222 223 for ( i=as; i.hasItem(); i++ ){224 elem= i.getItem();225 x= elem.mvar();226 if ( degree(elem,x) > 1){ // otherwise it's not an extension227 //if ( degree(f,x) > 0 ){ // does it occure in f? RICHTIG?228 Astar.append(elem);229 ord.append(x);230 //}231 }232 }233 uord= Difference(uord,ord);234 DEBOUTLN(CERR, "Astar is: ", Astar);235 DEBOUTLN(CERR, "ord is: ", ord);236 DEBOUTLN(CERR, "uord is: ", uord);237 238 // 3) second trivial cases239 if ( Astar.length() == 0 ){240 success=1;241 DEBDECLEVEL(CERR,"factoras");242 return CFFList(CFFactor(f,1));243 }244 245 // 4) Try to obtain a partial factorization using prop2 and prop3246 // Not yet implemented247 248 // 5) choose (d_1, ..., d_s)249 FFRandom gen;250 if ( getCharacteristic() == 0 ) IntRandom gen(100);251 REvaluation D(1, ord.length(), gen);252 //D.nextpoint(); // first we use 0 as evaluation point253 254 F2: //[Characteristic set computation]255 // 6) Compute f^* (<- fstar)256 CanonicalForm substhin=CanonicalForm(vf), substback=CanonicalForm(vf);257 CanonicalForm monom;258 DEBOUTLN(CERR, "substhin= ", substhin);259 DEBOUTLN(CERR, "substback= ", substback);260 int j=1;261 for ( VarlistIterator jjj=ord; jjj.hasItem(); jjj++){262 if ( getCharacteristic() > 0 ) monom= jjj.getItem()*D[j];263 else monom= mapinto(jjj.getItem()*mapinto(D[j]));264 j++;265 substhin-= monom; substback+= monom;266 }267 DEBOUTLN(CERR, "substhin= ", substhin);268 DEBOUTLN(CERR, "substback= ", substback);269 CanonicalForm fstar=f(substhin,vf);270 DEBOUTLN(CERR, "fstar= ", fstar);271 272 // 7) Set up Variable ordering from ord,vf to vf,ord ! i.e. vf is the variable273 // with lowest level.274 Varlist nord=uord;275 nord.append(vf); nord= Union(ord,nord);276 DEBOUTLN(CERR, " nord= ", nord);277 CFList Astarnord= Astar; Astarnord.insert(fstar);278 DEBOUTLN(CERR, " Astarnord= ", Astarnord);279 Astarnord= reorder(nord,Astarnord);280 DEBOUTLN(CERR, "original Astar= ", Astar);281 DEBOUTLN(CERR, "reorderd Astar= ", Astarnord);282 DEBOUTLN(CERR, " f= ", f);283 DEBOUTLN(CERR, " fstar= ", fstar);284 285 // 8) Compute Charset Cstar of Astar \cup { fstar } wrt. ordering {vf, ord}286 PremForm Remembern;287 CFList Cstar= MCharSetN(Astarnord, Remembern);288 DEBOUTLN(CERR, " Cstar= ", Cstar );289 DEBOUTLN(CERR, " Factors removed= ", Remembern.FS2 );290 DEBOUTLN(CERR, " Possible Factors considered:= ", Remembern.FS1 );291 DEBOUTLN(CERR, " Reorderd Cstar= ", reorder(nord,Cstar));292 293 // 9) Compute Delta: the set of all irr. factors (over K_0) of initials of294 // Cstar295 CFList iniset= initalset1(Cstar);296 DEBOUTLN(CERR, "Set of initials: ", iniset);297 CFFList Delta;298 CFFList temp;299 for ( i=iniset; i.hasItem(); i++){300 Off(SW_RATIONAL);301 temp= Factorize(i.getItem());302 On(SW_RATIONAL);303 temp.removeFirst();304 Delta= Union(temp,Delta);305 }306 DEBOUTLN(CERR, "Delta= ", Delta);307 308 // 10) Compute Psi: the irreduzible factors (over K_0) of the first polynomial309 // in Cstar, which are not factors of any polynomial in Delta310 Off(SW_RATIONAL);311 CFFList Psi=Factorize(Cstar.getFirst());312 On(SW_RATIONAL);313 Psi.removeFirst();314 Psi= myDifference(Psi,Delta);315 DEBOUTLN(CERR, "Psi= ", Psi);316 317 // F3: [Test quasilinearity]318 // If Cstar is quasilinear -> F4319 // else if Psi.length() > 1 then Delta<- Delta \cup Psi; Psi<- \emptyset -> F4320 // else D.nextpoint() -> F2321 if ( isquasilinear(Cstar) ) {322 DEBOUTLN(CERR, "Cstar is quasilinear; going to F4, Cstar= ", Cstar);323 goto F4;324 }325 else if ( Psi.length() > 1 ){326 Delta= Union(Psi,Delta);327 Psi= CFFList();328 DEBOUTMSG(CERR, "Psi.length() == 1; going to F4");329 goto F4;330 }331 else{332 D.nextpoint();333 DEBOUTMSG(CERR, "Choosing next evaluation point. going to F2");334 goto F2;335 }336 337 F4: //[GCD Computation]338 CanonicalForm g=f;339 Delta= reorder(nord,Delta);340 Psi= reorder(nord,Psi);341 DEBOUTLN(CERR, "Reordered: Delta= ", Delta);342 DEBOUTLN(CERR, " Psi= ", Psi);343 CanonicalForm fp;344 345 DEBOUTMSG(CERR, "Testing Psi: this gives irreducible Factors!");346 for (jj=Psi; jj.hasItem(); jj++){347 if ( degree(g,vf) == 1 ) { // g is linear348 break;349 }350 fp= jj.getItem().factor();351 DEBOUT(CERR, "Calculating fp= gcd(", g);352 DEBOUT(CERR, ",", fp(substback,vf));353 DEBOUT(CERR, ") over K_r wrt ", vf);354 fp=alg_gcd(g,fp(substback,vf), as);355 //fp= algcd(g,fp(substback,vf), as, oldord);356 DEBOUTLN(CERR, " = ", fp);357 if ( degree(fp,vf) > 0 ){ //otherwise it's a constant358 g= divide(g, fp,as);359 DEBOUTLN(CERR, "f/fp= ", g);360 result.append(CFFactor(fp,1));361 }362 }363 364 DEBOUTMSG(CERR, "Testing Delta: this gives Factors (not nec. irreduzible!)");365 for (jj=Delta; jj.hasItem(); jj++){366 if ( degree(g,vf) <= 1 ) { // g is linear (or a constant..)367 break;368 }369 fp= jj.getItem().factor();370 DEBOUT(CERR, "Calculating fp= gcd(", g);371 DEBOUT(CERR, ",", fp(substback,vf));372 DEBOUT(CERR, ") over K_r wrt ", vf);373 fp= alg_gcd(g,fp(substback,vf), as);374 //fp= algcd(g,fp(substback,vf), as, oldord);375 DEBOUTLN(CERR, " = ", fp);376 if ( degree(fp,vf) > 0 ){ //otherwise it's a constant377 g= divide(g, fp,as);378 DEBOUTLN(CERR, "f/fp= ", g);379 // reduceresult.append(fp); // a facctor but not nec. irreduzible380 result.append(CFFactor(fp,1));381 success=0;382 }383 }384 385 // if (reduceresult.length() > 0){ //we have factors;386 // for ( CFListIterator I=reduceresult; I.hasItem(); I++ ){387 // // try to obtain an irreducible factorization388 // result= Union( result,389 // factoras(I.getItem(), Astar, success) );390 // }391 // }392 393 if ( result.length()==0 ){394 if ( isquasilinear(Cstar) ){ // Cstar quasilinear => f is irreduzible395 success=1;396 DEBDECLEVEL(CERR,"factoras");397 return CFFList(CFFactor(f,1));398 }399 else {400 DEBOUTMSG(CERR, "Going to F2!");401 D.nextpoint(); // choose a new set of int's.402 goto F2;403 }404 }405 // result.length() > 0 ==> we have factors!406 success=1;407 result.append(CFFactor(g,1)); //append the rest408 DEBDECLEVEL(CERR,"factoras");409 return result;410 }411 412 // for now: success=-1 no success, success=1 factored poly,413 // success = 0 factored poly partially414 CFFList415 cfactor(const CanonicalForm & f, const CFList & as, int success ){416 Off(SW_RATIONAL);417 CFFList Output, output, Factors=Factorize(f); On(SW_RATIONAL);418 int csuccess=0;419 Factors.removeFirst();420 421 if ( as.length() == 0 ){ success=1; return Factors;}422 if ( cls(f) <= cls(as.getLast()) ) { success=1; return Factors;}423 424 success=1;425 for ( CFFListIterator i=Factors; i.hasItem(); i++ ){426 CFFList CERR=factoras(i.getItem().factor(),as,csuccess);427 success= min(success,csuccess);428 for ( CFFListIterator j=CERR; j.hasItem(); j++)429 Output = myappend(Output,CFFactor(j.getItem().factor(),j.getItem().exp()*i.getItem().exp()));430 }431 return Output;432 }433 434 // Gegenbeispiel: char(k)=0 (original maple von wang):435 // f:= x^2*z^2-6909*x*z^2-3537*z^2-14344*x^2;436 // as:= [x^3-9977*x^2+15570*x-961]437 // ord:= [x,z]438 // cfactor(f,as,ord)=439 // -1/153600624478379*(-153600624478379*z^2+96784194432*x^2-246879593179080*x+440 // 13599716437688)*(x^2-6909*x-3537)441 //442 // maple:443 // alias(alpha=RootOf(x^3-9977*x^2+15570*x-961,x));444 // factor(f,alpha); // faktorisierung in char=0445 // x^2*z^2-6909*x*z^2-3537*z^2-14344*x^2446 -
libfac/charset/algfactor.h
r8b2ef5 r8eb6da 12 12 13 13 CanonicalForm algcd(const CanonicalForm & f, const CanonicalForm & g, const CFList & as, const Varlist & order); 14 CFFList factoras( const CanonicalForm & f, const CFList & as, int & success );15 CFFList cfactor(const CanonicalForm & f, const CFList & as, int success );16 14 /*BEGINPUBLIC*/ 17 15 int hasVar(const CanonicalForm &f, const Variable &v);
Note: See TracChangeset
for help on using the changeset viewer.