[e4fe2b] | 1 | #include "config.h" |
---|
| 2 | |
---|
[9c3d69] | 3 | #include "canonicalform.h" |
---|
[e4fe2b] | 4 | |
---|
[6ead9d] | 5 | #ifdef HAVE_BIFAC |
---|
[e4fe2b] | 6 | |
---|
| 7 | # include "lgs.h" |
---|
[0c4a34b] | 8 | #include "bifacConfig.h" |
---|
[6ead9d] | 9 | |
---|
[806c18] | 10 | #define BIFAC_BASIS_OF_G_CHECK 1 |
---|
[e4fe2b] | 11 | |
---|
[6ead9d] | 12 | void Reduce( bool ); |
---|
| 13 | CanonicalForm Bigcd( const CanonicalForm& f, const CanonicalForm& g); |
---|
| 14 | |
---|
| 15 | CanonicalForm MyContent( const CanonicalForm& F, const Variable& x); |
---|
| 16 | CFFList Mysqrfree( const CanonicalForm& F ); |
---|
| 17 | |
---|
| 18 | |
---|
| 19 | CanonicalForm MyGCDmod( const CanonicalForm & a,const CanonicalForm & b); |
---|
| 20 | CFFList RelFactorize(const CanonicalForm & h); |
---|
| 21 | |
---|
| 22 | //====== global definitions =================== |
---|
| 23 | Variable x( 'x' ); |
---|
| 24 | Variable y( 'y' ); |
---|
| 25 | Variable z( 'z' ); |
---|
| 26 | Variable e( 'e' ); |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | /////////////////////////////////////////////////////// |
---|
| 30 | // Class for storing polynomials as vectors. |
---|
| 31 | // Enables fast access to a certain degree. |
---|
| 32 | /////////////////////////////////////////////////////// |
---|
| 33 | |
---|
| 34 | //================================================== |
---|
[806c18] | 35 | class PolyVector |
---|
[6ead9d] | 36 | //================================================== |
---|
| 37 | { |
---|
| 38 | public: |
---|
| 39 | PolyVector ( void ){ |
---|
| 40 | m = -1; |
---|
| 41 | } |
---|
| 42 | virtual ~PolyVector( void ){ |
---|
| 43 | if( m!= -1) delete[] value; |
---|
| 44 | } |
---|
| 45 | void init (CanonicalForm f){ |
---|
| 46 | if( f.level()<0 ) |
---|
| 47 | { |
---|
| 48 | m = 0; |
---|
| 49 | n = 0; |
---|
| 50 | value = new CanonicalForm[1]; |
---|
| 51 | value[0] = f; |
---|
| 52 | } |
---|
| 53 | else |
---|
| 54 | { |
---|
| 55 | m = degree(f,x); |
---|
| 56 | n = degree(f,y); |
---|
| 57 | ASSERT( m>0 || n>0, "Input is not a polynomial"); |
---|
| 58 | int correction = 1; // univariate polynomials |
---|
| 59 | if( n==0) correction = n+1; |
---|
[806c18] | 60 | |
---|
[6ead9d] | 61 | value = new CanonicalForm[m*(n+1)+n+1]; |
---|
| 62 | for(int i=0; i<=m*(n+1)+n; i++) value[i]=0; |
---|
| 63 | |
---|
[806c18] | 64 | |
---|
| 65 | for ( CFIterator i = f; i.hasTerms(); i++ ) { |
---|
| 66 | for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ){ |
---|
| 67 | if( i.coeff().mvar().level()< 0 ){ |
---|
| 68 | value[ 0*(n+1) + i.exp()*correction ] = j.coeff();} |
---|
| 69 | else{ |
---|
| 70 | value[ j.exp()*(n+1) + i.exp()*correction ] = j.coeff();}}} |
---|
[6ead9d] | 71 | } |
---|
| 72 | } |
---|
| 73 | |
---|
| 74 | void push(int mm, int nn, CanonicalForm v){ |
---|
| 75 | ASSERT( 0<=mm<=m && 0<=nn<=n, "Wrong Index in PolyVector"); |
---|
| 76 | value[mm*(n+1)+nn] = v; |
---|
| 77 | } |
---|
| 78 | CanonicalForm get(int mm, int nn){ |
---|
| 79 | |
---|
| 80 | if( 0<=mm && mm<=m && 0<=nn && nn<=n ) |
---|
| 81 | return value[mm*(n+1)+nn]; |
---|
| 82 | else |
---|
| 83 | return 0; |
---|
| 84 | } |
---|
| 85 | #ifndef NOSTREAMIO |
---|
[181148] | 86 | friend OSTREAM & operator<< ( OSTREAM & s, const PolyVector& V ){ |
---|
[6ead9d] | 87 | for (int i=0;i<=V.m;i++) |
---|
| 88 | { |
---|
| 89 | s << "["; |
---|
| 90 | for (int j=0;j<=V.n;j++) |
---|
[806c18] | 91 | s << V.value[i*(V.n+1)+j] << ", "; |
---|
[6ead9d] | 92 | s << "]\n"; |
---|
| 93 | } |
---|
| 94 | return s; |
---|
| 95 | } |
---|
| 96 | #endif /* NOSTREAMIO */ |
---|
| 97 | |
---|
| 98 | |
---|
| 99 | private: |
---|
| 100 | int m; // Degree in x |
---|
| 101 | int n; // Degree in y |
---|
| 102 | CanonicalForm* value; // Value: index = m*(n+1)+n |
---|
| 103 | }; |
---|
| 104 | ////////// END of PolyVector /////////////////////////// |
---|
| 105 | |
---|
| 106 | |
---|
| 107 | |
---|
| 108 | |
---|
| 109 | ///////////////////////////////////////////////////////// |
---|
| 110 | // |
---|
| 111 | // Default class declarations |
---|
| 112 | // |
---|
| 113 | ///////////////////////////////////////////////////////// |
---|
| 114 | |
---|
| 115 | |
---|
| 116 | |
---|
| 117 | //--<>--------------------------------- |
---|
| 118 | BIFAC::BIFAC( void )// KONSTRUKTOR |
---|
| 119 | //--<>--------------------------------- |
---|
| 120 | { |
---|
| 121 | } |
---|
| 122 | |
---|
| 123 | //--<>--------------------------------- |
---|
| 124 | BIFAC::~BIFAC( void )// DESTRUKTOR |
---|
| 125 | //--<>--------------------------------- |
---|
| 126 | { |
---|
[806c18] | 127 | } |
---|
[6ead9d] | 128 | |
---|
| 129 | |
---|
| 130 | ///////////////////////////////////////////////////////// |
---|
| 131 | // |
---|
| 132 | // Auxiliary functions |
---|
| 133 | // |
---|
| 134 | ///////////////////////////////////////////////////////// |
---|
| 135 | |
---|
| 136 | // //--<>--------------------------------- |
---|
| 137 | // void BIFAC::matrix_drucken( CFMatrix M ) |
---|
| 138 | // //--<>--------------------------------- |
---|
| 139 | // { |
---|
| 140 | // int i,j; |
---|
| 141 | // char* name="matrix.ppm"; |
---|
| 142 | |
---|
| 143 | // // === Datei löschen === |
---|
[806c18] | 144 | |
---|
| 145 | // ofstream* aus = new ofstream(name, ios::out); |
---|
[6ead9d] | 146 | // delete aus; |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | // // === Jetzt immer nur anhängen === |
---|
| 150 | |
---|
[806c18] | 151 | // aus = new ofstream(name, ios::app); |
---|
[6ead9d] | 152 | // *aus << "// Zeilen Spalten\n" |
---|
| 153 | // << "// x-Koord. y-Koord. Wert\n"; |
---|
| 154 | |
---|
| 155 | // *aus << M.rows() << " " << M.columns() << endl; |
---|
| 156 | |
---|
| 157 | |
---|
| 158 | // // === Noch nicht bearbeitet Teile === |
---|
| 159 | |
---|
| 160 | // for( i=0; i<M.rows(); i++) |
---|
| 161 | // for( j=0; j<M.columns(); j++) |
---|
| 162 | // *aus << i << " " << j << " " << M(i+1,j+1) << endl;; |
---|
| 163 | // delete aus; |
---|
[806c18] | 164 | // } |
---|
[6ead9d] | 165 | |
---|
| 166 | //======================================================= |
---|
| 167 | void BIFAC::passedTime() |
---|
| 168 | //======================================================= |
---|
| 169 | { |
---|
[806c18] | 170 | ; |
---|
[6ead9d] | 171 | } |
---|
| 172 | |
---|
| 173 | |
---|
| 174 | //======================================================= |
---|
| 175 | long int BIFAC::anz_terme( CanonicalForm & f ) |
---|
| 176 | //======================================================= |
---|
| 177 | { |
---|
| 178 | long int z=0; |
---|
| 179 | |
---|
| 180 | for ( CFIterator i = f; i.hasTerms(); i++ ) |
---|
[806c18] | 181 | for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) |
---|
[6ead9d] | 182 | z++; |
---|
| 183 | return( z ); |
---|
| 184 | } |
---|
[806c18] | 185 | |
---|
[6ead9d] | 186 | //======================================================= |
---|
| 187 | void BIFAC::biGanzMachen( CanonicalForm & f ) |
---|
| 188 | //======================================================= |
---|
| 189 | { |
---|
| 190 | CanonicalForm ggT; |
---|
| 191 | bool init = false; |
---|
| 192 | Off( SW_RATIONAL ); |
---|
| 193 | |
---|
| 194 | for ( CFIterator i = f; i.hasTerms(); i++ ) |
---|
[806c18] | 195 | for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) |
---|
[6ead9d] | 196 | { |
---|
| 197 | if( !init ) |
---|
[806c18] | 198 | { |
---|
| 199 | ggT = j.coeff(); |
---|
| 200 | init = true; |
---|
[6ead9d] | 201 | } |
---|
[806c18] | 202 | else |
---|
| 203 | ggT = gcd(j.coeff(), ggT); |
---|
[6ead9d] | 204 | } |
---|
| 205 | f /= ggT; |
---|
| 206 | On( SW_RATIONAL ); |
---|
| 207 | } |
---|
| 208 | |
---|
| 209 | //======================================================= |
---|
[806c18] | 210 | void BIFAC::biNormieren( CanonicalForm & f ) |
---|
[6ead9d] | 211 | //======================================================= |
---|
| 212 | { |
---|
| 213 | if ( getCharacteristic() == 0 ) |
---|
| 214 | { |
---|
| 215 | for ( CFIterator i = f; i.hasTerms(); i++ ) |
---|
[806c18] | 216 | for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) |
---|
| 217 | if( j.coeff().den() != 1 ) |
---|
| 218 | { |
---|
| 219 | f *= j.coeff().den(); |
---|
| 220 | biNormieren( f ); |
---|
| 221 | } |
---|
[6ead9d] | 222 | biGanzMachen( f ); |
---|
[806c18] | 223 | } |
---|
[6ead9d] | 224 | else |
---|
| 225 | { |
---|
| 226 | f /= LC(f); |
---|
| 227 | } |
---|
| 228 | } |
---|
| 229 | |
---|
| 230 | |
---|
| 231 | //======================================================= |
---|
| 232 | // * Convert the basis vectors of G into polynomials |
---|
| 233 | // * Validate the solutions |
---|
| 234 | //======================================================= |
---|
| 235 | CFList BIFAC::matrix2basis(CFMatrix A, int dim, int m, int n, CanonicalForm f) |
---|
| 236 | //======================================================= |
---|
| 237 | { |
---|
| 238 | Variable x('x'), y('y'); |
---|
| 239 | int i,j,k; |
---|
| 240 | CanonicalForm g,h,ff; |
---|
| 241 | CFList Lg, Lh; |
---|
| 242 | |
---|
| 243 | // === Construction of the 'g's ==== |
---|
| 244 | for(k=1; k<=dim; k++) |
---|
| 245 | { |
---|
| 246 | g=0; |
---|
| 247 | for(i=0; i<=m-1; i++) |
---|
| 248 | for(j=0; j<=n; j++) |
---|
[806c18] | 249 | g += A(k, i*(n+1)+j+1)* power(x,i) * power(y,j); |
---|
[6ead9d] | 250 | Lg.append(g); |
---|
| 251 | } |
---|
| 252 | |
---|
| 253 | /////////// START VALIDATION //////////////////////////////////// |
---|
| 254 | if (BIFAC_BASIS_OF_G_CHECK) |
---|
| 255 | { |
---|
| 256 | |
---|
| 257 | // === Construction of the 'h's ==== |
---|
| 258 | for(k=1; k<=dim; k++) |
---|
| 259 | { |
---|
| 260 | h=0; |
---|
| 261 | for(i=0; i<=m; i++) |
---|
[806c18] | 262 | for(j=0; j<n; j++) |
---|
| 263 | h += A(k, i*n+j+1 +m*(n+1))* power(x,i) * power(y,j); |
---|
[6ead9d] | 264 | Lh.append(h); |
---|
| 265 | } |
---|
[806c18] | 266 | |
---|
[6ead9d] | 267 | // === Is the solution correct? === |
---|
| 268 | CFListIterator itg=Lg; |
---|
| 269 | CFListIterator ith=Lh; |
---|
| 270 | for( ; itg.hasItem(); itg++, ith++) |
---|
| 271 | { |
---|
| 272 | g = itg.getItem(); |
---|
| 273 | h = ith.getItem(); |
---|
| 274 | ff = f*(deriv(g,y)-deriv(h,x)) +h*deriv(f,x) -g*deriv(f,y); |
---|
[806c18] | 275 | if( !ff.isZero()) { |
---|
[0c4a34b] | 276 | #ifndef NOSTREAMIO |
---|
[806c18] | 277 | AUSGABE_ERR("* Falsche Polynome!"); |
---|
| 278 | exit (1); |
---|
[0c4a34b] | 279 | #else |
---|
| 280 | printf("wrong polys\n"); |
---|
| 281 | break; |
---|
| 282 | #endif |
---|
[6ead9d] | 283 | } |
---|
[806c18] | 284 | } |
---|
[6ead9d] | 285 | } |
---|
| 286 | /////////// END VALIDATION //////////////////////////////////// |
---|
| 287 | |
---|
| 288 | return (Lg); |
---|
| 289 | } |
---|
| 290 | |
---|
| 291 | ///////////////////////////////////////////////////////// |
---|
| 292 | // |
---|
| 293 | // Main functions |
---|
| 294 | // |
---|
| 295 | ///////////////////////////////////////////////////////// |
---|
| 296 | |
---|
| 297 | //======================================================= |
---|
| 298 | // * Create the matrix belonging to G |
---|
| 299 | // * Compute a basis of the kernel |
---|
| 300 | //======================================================= |
---|
| 301 | CFList BIFAC::basisOfG(CanonicalForm f) |
---|
| 302 | //======================================================= |
---|
| 303 | { |
---|
| 304 | |
---|
| 305 | |
---|
| 306 | int m = degree(f,x); |
---|
| 307 | int n = degree(f,y); |
---|
| 308 | int r,s, ii,jj; |
---|
[806c18] | 309 | |
---|
[6ead9d] | 310 | |
---|
| 311 | // ======= Creation of the system of linear equations for G ============= |
---|
| 312 | int rows = 4*m*n; |
---|
| 313 | int columns = m*(n+1) + (m+1)*n; |
---|
| 314 | |
---|
| 315 | CFMatrix M(rows, columns); // Remember: The first index is (1,1) -- not (0,0)! |
---|
[806c18] | 316 | |
---|
[6ead9d] | 317 | for ( CFIterator i = f; i.hasTerms(); i++ ) // All coeffizients of y |
---|
| 318 | { |
---|
| 319 | for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) // All coeffizients of x |
---|
| 320 | { |
---|
| 321 | r = j.exp(); // x^r |
---|
| 322 | s = i.exp(); // y^s |
---|
| 323 | |
---|
| 324 | // Now we regard g_{ii,jj) |
---|
| 325 | for( ii=0; ii<m; ii++) |
---|
[806c18] | 326 | for( jj=0; jj<=n; jj++) |
---|
| 327 | { |
---|
| 328 | if( s>= 1) M( (r+ii)*2*n +(jj+s-1)+1, ii*(n+1)+jj +1) += -j.coeff() * s; |
---|
| 329 | if( jj>= 1) M( (r+ii)*2*n +(jj+s-1)+1, ii*(n+1)+jj +1) += j.coeff() * jj; |
---|
| 330 | } |
---|
[6ead9d] | 331 | |
---|
| 332 | // Now we regard h_{ii,jj} |
---|
| 333 | for( ii=0; ii<=m; ii++) |
---|
[806c18] | 334 | for( jj=0; jj<n; jj++) |
---|
| 335 | { |
---|
| 336 | if( r>= 1) M( (r+ii-1)*2*n +(jj+s)+1, (ii*n)+jj +m*(n+1) +1) += j.coeff() * r; |
---|
| 337 | if( ii>= 1) M( (r+ii-1)*2*n +(jj+s)+1, (ii*n) +jj +m*(n+1) +1) += -j.coeff() * ii; |
---|
| 338 | } |
---|
[6ead9d] | 339 | } |
---|
| 340 | } |
---|
| 341 | // ========= Solving the system of linear equations for G ============= |
---|
| 342 | |
---|
| 343 | // matrix_drucken(M); // ********************************** |
---|
| 344 | |
---|
| 345 | LGS L(rows,columns); |
---|
| 346 | |
---|
| 347 | CFMatrix Z(1,columns); |
---|
| 348 | for( ii=1; ii<=rows; ii++) |
---|
| 349 | { |
---|
| 350 | for( jj=1; jj<=columns; jj++) |
---|
| 351 | Z(1,jj) = M(ii,jj); // Copy the ii-th row |
---|
| 352 | L.new_row(Z); |
---|
| 353 | } |
---|
| 354 | |
---|
| 355 | if( L.corank() == 1 ) |
---|
| 356 | { |
---|
| 357 | CFList Lg; |
---|
| 358 | Lg.append(f); |
---|
| 359 | return(Lg); |
---|
| 360 | } |
---|
| 361 | // L.print(); |
---|
| 362 | CFMatrix basis = L.GetKernelBasis(); |
---|
| 363 | |
---|
| 364 | // ============= TEST AUF KORREKTHEIT /start) ==== |
---|
| 365 | CanonicalForm tmp; |
---|
| 366 | for(int k=1; k<= L.corank(); k++) |
---|
| 367 | for(int i=1; i<=rows; i++) |
---|
| 368 | { |
---|
| 369 | tmp =0; |
---|
| 370 | for(int j=1; j<=columns; j++) |
---|
[806c18] | 371 | tmp += M(i,j) * basis(k,j); |
---|
[6ead9d] | 372 | if( tmp!= 0) { |
---|
[806c18] | 373 | exit(17); |
---|
[6ead9d] | 374 | } |
---|
| 375 | } |
---|
| 376 | // ============= TEST AUF KORREKTHEIT (ende) ==== |
---|
| 377 | return ( matrix2basis( basis, L.corank(), m,n,f ) ); |
---|
| 378 | } |
---|
| 379 | |
---|
| 380 | //======================================================= |
---|
[806c18] | 381 | // Compute a r x r - matrix A=(a_ij) for |
---|
[6ead9d] | 382 | // gg_i = SUM a_ij * g_j * f_x (mod f) |
---|
[806c18] | 383 | // Return a list consisting of |
---|
[6ead9d] | 384 | // r x (r+1) Matrix A |
---|
[806c18] | 385 | // the last columns contains only the indices of the |
---|
[6ead9d] | 386 | // first r linear independent lines |
---|
| 387 | // REMARK: this is used by BIFAC::createEg but NOT by createEgUni!! |
---|
| 388 | //======================================================= |
---|
| 389 | CFMatrix BIFAC::createA (CFList G, CanonicalForm f) |
---|
| 390 | //======================================================= |
---|
| 391 | { |
---|
| 392 | // === Declarations === |
---|
| 393 | int m,n; |
---|
| 394 | int i,j,e; |
---|
| 395 | int r = G.length(); // number of factors |
---|
| 396 | |
---|
[806c18] | 397 | LGS L(r,r,true); |
---|
| 398 | // LGS L(r,r); |
---|
[6ead9d] | 399 | CFMatrix Z(1,r); |
---|
| 400 | CFMatrix A(r,r+2); // the last two column contain the bi-degree |
---|
| 401 | |
---|
| 402 | CanonicalForm fx = deriv(f,x); |
---|
| 403 | PolyVector* gifx = new PolyVector[r]; |
---|
| 404 | |
---|
| 405 | // === Convert polynomials into vectors === |
---|
| 406 | i=0; |
---|
| 407 | CanonicalForm q; |
---|
| 408 | |
---|
[806c18] | 409 | for( CFListIterator it=G; it.hasItem(); it++, i++){ |
---|
[6ead9d] | 410 | |
---|
| 411 | gifx[i].init( (it.getItem()*fx)%f ); |
---|
| 412 | } |
---|
| 413 | |
---|
| 414 | // === Search linear independent lines === |
---|
| 415 | |
---|
| 416 | e=1; // row number of A |
---|
[806c18] | 417 | n=0; // |
---|
[6ead9d] | 418 | m=0; // |
---|
| 419 | while (L.rank() != r ) |
---|
| 420 | { |
---|
| 421 | for(j=0;j<r;j++) |
---|
| 422 | Z(1,j+1) = gifx[j].get(m,n); |
---|
| 423 | if( L.new_row(Z,0) ) // linear independent row? |
---|
| 424 | { |
---|
| 425 | ASSERT( e<=r, "Wrong index in matrix"); |
---|
| 426 | A(e,r+1) = m; // Degree in x |
---|
| 427 | A(e,r+2) = n; // Degree in y |
---|
| 428 | e++; |
---|
| 429 | } |
---|
| 430 | if (m>n) n++; |
---|
| 431 | else { m++; n=0; } |
---|
| 432 | } |
---|
| 433 | L.print(); |
---|
| 434 | |
---|
| 435 | L.inverse(A); |
---|
| 436 | |
---|
| 437 | // === Clean up == |
---|
| 438 | delete[] gifx; |
---|
| 439 | |
---|
| 440 | return A; |
---|
| 441 | } |
---|
| 442 | |
---|
| 443 | //======================================================= |
---|
| 444 | CanonicalForm BIFAC::create_g (CFList G) |
---|
| 445 | //======================================================= |
---|
| 446 | { |
---|
| 447 | CanonicalForm g = 0; |
---|
| 448 | int i = 0; |
---|
| 449 | |
---|
| 450 | int r = G.length(); // number of factors |
---|
| 451 | float SS = 10*( r*(r-1) / ( 2*( (100- (float) EgSeparable)/100)) ); |
---|
| 452 | int S = (int) SS +1; |
---|
| 453 | |
---|
| 454 | IntRandom RANDOM(S); |
---|
| 455 | |
---|
| 456 | int* rand_coeff1 = new int[r]; |
---|
| 457 | |
---|
| 458 | |
---|
| 459 | // static for debugging |
---|
| 460 | // rand_coeff1[0] = 12; rand_coeff1[1] = 91; rand_coeff1[2] = 42; |
---|
| 461 | |
---|
| 462 | for( CFListIterator it=G; it.hasItem(); it++, i++) |
---|
| 463 | { |
---|
| 464 | rand_coeff1[i] = RANDOM.generate().intval(); |
---|
| 465 | |
---|
| 466 | g += rand_coeff1[i] * it.getItem(); |
---|
| 467 | } |
---|
[806c18] | 468 | |
---|
[6ead9d] | 469 | delete[] rand_coeff1; |
---|
| 470 | |
---|
| 471 | return g; |
---|
| 472 | } |
---|
| 473 | |
---|
| 474 | ///////////////////////////////////////////////////////////// |
---|
| 475 | // This subroutine creates the polynomials Eg(x) and g |
---|
| 476 | // by using the 'bivariate' methode'. |
---|
| 477 | // REMARK: There is a 'univariate methode' as well |
---|
| 478 | // which ought to be faster! |
---|
| 479 | //////////////////////////////////////////////////////////// |
---|
| 480 | //======================================================= |
---|
| 481 | CFList BIFAC::createEg (CFList G, CanonicalForm f) |
---|
| 482 | //======================================================= |
---|
| 483 | { |
---|
| 484 | |
---|
| 485 | |
---|
| 486 | CFMatrix NEU = createA(G,f); |
---|
| 487 | // passedTime(); |
---|
| 488 | |
---|
| 489 | bool suitable1 = false; // Is Eg by chance unsuitable? |
---|
[806c18] | 490 | bool suitable2 = false; // Is on of g*g_i or g_i*f_x zero? |
---|
[6ead9d] | 491 | |
---|
| 492 | // === (0) Preparation === |
---|
| 493 | CanonicalForm g; |
---|
| 494 | CanonicalForm Eg; |
---|
| 495 | CanonicalForm fx = deriv(f,x); |
---|
| 496 | |
---|
| 497 | int i,j,e; |
---|
| 498 | int r = G.length(); // number of factors |
---|
| 499 | // float SS = ( r*(r-1) / ( 2*( (100- (float) EgSeparable)/100)) ); |
---|
| 500 | // int S = (int) SS +1; |
---|
| 501 | |
---|
| 502 | // IntRandom RANDOM(S); |
---|
| 503 | // int* rand_coeff = new int[r]; |
---|
| 504 | CFMatrix A(r,r); |
---|
| 505 | CanonicalForm* gi = new CanonicalForm[r]; |
---|
| 506 | CanonicalForm* ggi = new CanonicalForm[r]; |
---|
| 507 | PolyVector* v_ggi = new PolyVector [r]; |
---|
| 508 | |
---|
| 509 | |
---|
| 510 | |
---|
| 511 | i=0; |
---|
| 512 | for( CFListIterator it=G; it.hasItem(); it++, i++) |
---|
| 513 | gi[i] = it.getItem(); |
---|
| 514 | |
---|
| 515 | while ( !suitable1 ) |
---|
| 516 | { |
---|
| 517 | |
---|
[806c18] | 518 | suitable2 = false; |
---|
[6ead9d] | 519 | // === (1) Creating g === |
---|
| 520 | while ( !suitable2 ) |
---|
| 521 | { |
---|
| 522 | // i=0; |
---|
| 523 | // g=0; |
---|
| 524 | // for( CFListIterator it=G; it.hasItem(); it++, i++) |
---|
| 525 | // { |
---|
[806c18] | 526 | // gi[i] = it.getItem(); |
---|
| 527 | // rand_coeff[i] = RANDOM.generate().intval(); |
---|
| 528 | // g += rand_coeff[i] * it.getItem(); |
---|
[6ead9d] | 529 | // } |
---|
| 530 | g = create_g( G ); |
---|
[806c18] | 531 | |
---|
[6ead9d] | 532 | // === (2) Computing g_i * g === |
---|
[806c18] | 533 | // |
---|
| 534 | for(i=0; i<r; i++){ |
---|
[6ead9d] | 535 | |
---|
[806c18] | 536 | ggi[i] = (g*gi[i])%f; // seite 10 |
---|
[6ead9d] | 537 | } |
---|
[806c18] | 538 | |
---|
[6ead9d] | 539 | // === Check if all polynomials are <> 0 === |
---|
| 540 | suitable2 = true; // It should be fine, but ... |
---|
[806c18] | 541 | if( g.isZero() ) |
---|
| 542 | suitable2 = false; |
---|
[6ead9d] | 543 | // else |
---|
[806c18] | 544 | // for(i=0; i<r; i++) |
---|
| 545 | // if( ggi[i].isZero() ) |
---|
| 546 | // suitable2 = false; |
---|
[6ead9d] | 547 | |
---|
| 548 | } // end of Žwhile ( !suitable2 )Ž |
---|
[806c18] | 549 | |
---|
[6ead9d] | 550 | // === (3) Computing Eg(x) === |
---|
[806c18] | 551 | |
---|
[6ead9d] | 552 | for(i=0;i<r;i++) // Get Polynomials as vectors |
---|
| 553 | v_ggi[i].init(ggi[i]); |
---|
| 554 | |
---|
| 555 | // Matrix A |
---|
[806c18] | 556 | for(i=1; i<=r; i++) |
---|
[6ead9d] | 557 | for( j=1; j<=r; j++) |
---|
| 558 | { |
---|
[806c18] | 559 | A(i,j) = 0; |
---|
| 560 | for( e=1; e<=r; e++) |
---|
| 561 | { |
---|
[6ead9d] | 562 | |
---|
| 563 | |
---|
| 564 | A(i,j) += ( NEU(j,e ) * v_ggi[i-1].get(NEU(e,r+1).intval(),(NEU(e,r+2).intval() ))); |
---|
| 565 | |
---|
| 566 | |
---|
| 567 | // |
---|
| 568 | |
---|
[806c18] | 569 | } |
---|
[6ead9d] | 570 | } |
---|
| 571 | |
---|
| 572 | for(j=1; j<=r; j++) |
---|
| 573 | A(j,j) -= x; |
---|
| 574 | Eg = determinant(A,r); |
---|
| 575 | // exit(1); |
---|
| 576 | // === (4) Is Eg(x) suitable? === |
---|
| 577 | if( MyGCDmod(Eg, deriv(Eg,x)) == 1 ) |
---|
| 578 | suitable1 = true; |
---|
| 579 | else |
---|
| 580 | { |
---|
| 581 | } |
---|
| 582 | } // end of Žwhile ( !suitable1 )Ž |
---|
| 583 | |
---|
| 584 | // Delete trash |
---|
| 585 | |
---|
| 586 | |
---|
| 587 | |
---|
| 588 | delete[] v_ggi; |
---|
| 589 | delete[] gi; |
---|
| 590 | delete[] ggi; |
---|
| 591 | // delete[] rand_coeff; |
---|
| 592 | |
---|
| 593 | CFList LL; |
---|
| 594 | LL.append(Eg); |
---|
| 595 | LL.append(g); |
---|
| 596 | return (LL); |
---|
| 597 | } |
---|
| 598 | // ///////////////////////////////////////////////////////////// |
---|
| 599 | // // It is possible to take univariate polynomials |
---|
| 600 | // // with y:=c for a suitable c. |
---|
| 601 | // // c is suitable iff gcd( f(x,c), f_x(x,c)) = 1. |
---|
| 602 | // //////////////////////////////////////////////////////////// |
---|
| 603 | // |
---|
| 604 | //======================================================= |
---|
| 605 | CFList BIFAC::createEgUni (CFList G, CanonicalForm f) |
---|
| 606 | //======================================================= |
---|
| 607 | { |
---|
| 608 | |
---|
| 609 | int i,ii,k; |
---|
| 610 | CanonicalForm ff, ffx,g, gg, Eg; |
---|
| 611 | |
---|
[806c18] | 612 | |
---|
[6ead9d] | 613 | bool suitable1 = false; // Is Eg unsuitable? |
---|
[806c18] | 614 | bool suitable2 = false; // Is on of g*g_i or g_i*f_x zero? |
---|
[6ead9d] | 615 | bool suitable3 = false; // Is 'konst' unsuitable? |
---|
| 616 | |
---|
| 617 | // ======================== |
---|
| 618 | // = (0) Preparation = |
---|
| 619 | // ======================== |
---|
| 620 | int konst = 0; |
---|
| 621 | CanonicalForm fx = deriv(f,x); |
---|
| 622 | int m = degree(f,x); |
---|
| 623 | int r = G.length(); // number of factors |
---|
| 624 | int S = (int) ((float) ( r*(r-1) / ( 2*( (100- (float) EgSeparable)/100)) )+1); |
---|
| 625 | |
---|
| 626 | |
---|
| 627 | int* rand_coeff = new int[r]; |
---|
| 628 | CanonicalForm* gi = new CanonicalForm[r]; |
---|
| 629 | CanonicalForm* ggi = new CanonicalForm[r]; |
---|
| 630 | |
---|
| 631 | CFMatrix A(r,r); // We have to find the matrix A, |
---|
| 632 | CFMatrix Z(1,r); // `Vector` for data transportation |
---|
| 633 | CFMatrix AA(m,r); // but first we generate AA. |
---|
[806c18] | 634 | CFMatrix AI(r,r+1); // |
---|
| 635 | LGS L(r,r,true); |
---|
[6ead9d] | 636 | IntRandom RANDOM(S); |
---|
| 637 | |
---|
| 638 | |
---|
| 639 | // ========================================================== |
---|
| 640 | // = (1) Find a suitable constant to make bivariate = |
---|
| 641 | // = polynomials univariate. Try the following numbers = |
---|
| 642 | // = 0, 1, -1, 2, -2, 3,... = |
---|
| 643 | // ========================================================== |
---|
| 644 | |
---|
| 645 | while ( !suitable3 ) |
---|
| 646 | { |
---|
| 647 | ff = f(konst,'y'); |
---|
| 648 | ffx = fx(konst,'y'); |
---|
[806c18] | 649 | |
---|
[6ead9d] | 650 | if( gcd(ff, ffx) == 1) |
---|
| 651 | suitable3 = true; |
---|
| 652 | else |
---|
| 653 | { |
---|
| 654 | konst *= -1; |
---|
| 655 | if( konst >= 0 ) |
---|
[806c18] | 656 | konst++; |
---|
[6ead9d] | 657 | } |
---|
| 658 | } |
---|
| 659 | |
---|
| 660 | |
---|
| 661 | // =============================================== |
---|
| 662 | // = (2) Make g_i univariate = |
---|
| 663 | // =============================================== |
---|
| 664 | i=0; |
---|
| 665 | for( CFListIterator it=G; it.hasItem(); it++, i++) |
---|
| 666 | { |
---|
[806c18] | 667 | gi[i] = it.getItem()(konst,'y'); |
---|
[6ead9d] | 668 | } |
---|
| 669 | |
---|
| 670 | // =============================================== |
---|
| 671 | // = (3) Compute the matrices 'AA' and 'AI' = |
---|
| 672 | // =============================================== |
---|
[806c18] | 673 | |
---|
[6ead9d] | 674 | |
---|
| 675 | for( i=0; i<r; i++) // First store all coeffizients in AA. |
---|
| 676 | { |
---|
| 677 | ggi[i] = (gi[i]*ffx)%ff; // now we have degree < m. |
---|
| 678 | //biNormieren(ggi[i]); |
---|
| 679 | for ( CFIterator j = ggi[i]; j.hasTerms(); j++ ) |
---|
[806c18] | 680 | AA( j.exp()+1, i+1) = j.coeff(); |
---|
[6ead9d] | 681 | } |
---|
| 682 | |
---|
| 683 | |
---|
| 684 | // Now find the lin. indep. rows. |
---|
| 685 | i = 1; |
---|
| 686 | ii = 1; // row number of A |
---|
| 687 | while (L.rank() != r ) |
---|
| 688 | { |
---|
| 689 | ASSERT( i<=m, "Too few linear independent rows!"); |
---|
| 690 | |
---|
[806c18] | 691 | for (k=1; k<=r; k++) |
---|
[6ead9d] | 692 | Z(1,k) = AA(i,k); |
---|
| 693 | if( L.new_row(Z,0) ) // linear independent row? |
---|
| 694 | { |
---|
| 695 | ASSERT( ii<=r, "Wrong index in matrix"); |
---|
| 696 | AI(ii,r+1) = i; // Degree in x |
---|
| 697 | ii++; |
---|
| 698 | } |
---|
| 699 | i++; |
---|
| 700 | L.print(); |
---|
| 701 | } |
---|
| 702 | L.inverse(AI); |
---|
| 703 | |
---|
| 704 | |
---|
| 705 | // ============================================== |
---|
| 706 | // = (4) Big loop to find a suitable 'Eg(x) = |
---|
| 707 | // ============================================== |
---|
[806c18] | 708 | |
---|
[6ead9d] | 709 | while ( !suitable1 ) // Is Eg(x) suitable? -> Check at the end of this procedure! |
---|
| 710 | { |
---|
| 711 | suitable2 = false; // In case we need a second loop |
---|
| 712 | |
---|
| 713 | // ================================================ |
---|
| 714 | // = (4a) Find a suitable 'g' = |
---|
| 715 | // ================================================ |
---|
| 716 | // rand_coeff[0] = 0; |
---|
| 717 | // rand_coeff[1] = 4; |
---|
[806c18] | 718 | |
---|
[6ead9d] | 719 | |
---|
| 720 | while ( !suitable2 ) |
---|
| 721 | { |
---|
| 722 | // === (i) Creating g === |
---|
| 723 | i=0; |
---|
| 724 | g=0; |
---|
| 725 | for( CFListIterator it=G; it.hasItem(); it++, i++) |
---|
| 726 | { |
---|
[806c18] | 727 | rand_coeff[i] = RANDOM.generate().intval(); |
---|
| 728 | g += rand_coeff[i] * it.getItem(); |
---|
[6ead9d] | 729 | } |
---|
| 730 | gg = g(konst,'y'); // univariate! |
---|
| 731 | for(i=0; i<r; i++) ggi[i] = (gi[i]*gg)%ff; // !! Redefinition of ggi !! |
---|
| 732 | |
---|
| 733 | // === (ii) Check if all polynomials are <> 0 === |
---|
| 734 | suitable2 = true; // It should be fine, but ... |
---|
[806c18] | 735 | if( gg.isZero() ) |
---|
| 736 | suitable2 = false; |
---|
[6ead9d] | 737 | // else |
---|
[806c18] | 738 | // for(i=0; i<r; i++) |
---|
| 739 | // if( ggi[i].isZero() ) |
---|
| 740 | // suitable2 = false; |
---|
[6ead9d] | 741 | } // end of Žwhile ( !suitable2 )Ž |
---|
[806c18] | 742 | |
---|
[6ead9d] | 743 | // createRg(g,f); |
---|
| 744 | |
---|
| 745 | // =============================================== |
---|
| 746 | // = (b) Compute matrix 'A' = |
---|
| 747 | // =============================================== |
---|
[806c18] | 748 | for(i=1; i<=r; i++) |
---|
[6ead9d] | 749 | { |
---|
[806c18] | 750 | for( ii=1; ii<=m; ii++) |
---|
| 751 | AA (ii,1) = 0; // !! Redefinition of AA !! |
---|
[6ead9d] | 752 | for ( CFIterator j = ggi[i-1]; j.hasTerms(); j++ ) |
---|
[806c18] | 753 | AA( j.exp()+1, 1) = j.coeff(); |
---|
[6ead9d] | 754 | |
---|
| 755 | for( ii=1; ii<=r; ii++) |
---|
| 756 | { |
---|
[806c18] | 757 | A(i,ii) = 0; |
---|
| 758 | for( k=1; k<=r; k++) |
---|
| 759 | A(i,ii) += ( AI(ii,k ) * AA( AI(k, r+1 ).intval(),1) ); |
---|
[6ead9d] | 760 | } |
---|
| 761 | } |
---|
| 762 | for(i=1; i<=r; i++) |
---|
| 763 | A(i,i) -= x; |
---|
| 764 | |
---|
| 765 | // =============================================== |
---|
| 766 | // = (c) Compute Eg(x) and check it = |
---|
| 767 | // =============================================== |
---|
[806c18] | 768 | |
---|
[6ead9d] | 769 | Eg = determinant(A,r); |
---|
| 770 | if( gcd(Eg, deriv(Eg,x)) == 1 ) |
---|
| 771 | { |
---|
| 772 | suitable1 = true; |
---|
| 773 | } |
---|
| 774 | } // end of Žwhile ( !suitable1 )Ž |
---|
| 775 | |
---|
[806c18] | 776 | |
---|
[6ead9d] | 777 | // ============================================== |
---|
| 778 | // = (5) Prepare for leaving = |
---|
| 779 | // ============================================== |
---|
| 780 | |
---|
| 781 | delete[] gi; |
---|
| 782 | delete[] ggi; |
---|
| 783 | delete[] rand_coeff; |
---|
[806c18] | 784 | |
---|
[6ead9d] | 785 | CFList LL; |
---|
| 786 | LL.append(Eg); |
---|
| 787 | LL.append(g); |
---|
| 788 | |
---|
| 789 | return (LL); |
---|
| 790 | } |
---|
| 791 | ///////////////////////////////////////////////////////////// |
---|
| 792 | // This subroutine creates the polynomials Rg(x) |
---|
| 793 | // which can be used instead of Eg(x). |
---|
| 794 | // No basis of G is neccessary, only one element |
---|
| 795 | //////////////////////////////////////////////////////////// |
---|
| 796 | //======================================================= |
---|
| 797 | CFList BIFAC::createRg (CFList G, CanonicalForm f) |
---|
| 798 | //======================================================= |
---|
| 799 | { |
---|
| 800 | |
---|
| 801 | // cerr << "* Was ist wenn g versagt???? -> Ausbauen\n"; |
---|
| 802 | |
---|
| 803 | CanonicalForm fx = deriv(f,x); |
---|
| 804 | CanonicalForm Rg; |
---|
| 805 | CanonicalForm g = create_g(G); |
---|
| 806 | |
---|
| 807 | |
---|
| 808 | // =============================================== |
---|
| 809 | // = (1) Find a suitable constant = |
---|
| 810 | // =============================================== |
---|
| 811 | |
---|
[806c18] | 812 | CanonicalForm alpha=1; |
---|
[6ead9d] | 813 | |
---|
[806c18] | 814 | while( resultant( f, fx, x)(alpha) == 0 ) |
---|
[6ead9d] | 815 | { |
---|
[806c18] | 816 | //while( resultant( f, fx, x)(alpha).inCoeffDomain() != true ) |
---|
[6ead9d] | 817 | //alpha +=1; |
---|
| 818 | } |
---|
| 819 | |
---|
| 820 | |
---|
| 821 | // =============================================== |
---|
| 822 | // = (2) Find a suitable constant = |
---|
| 823 | // =============================================== |
---|
[806c18] | 824 | |
---|
[6ead9d] | 825 | Rg = resultant( f(alpha,y), g(alpha,y)-z*fx(alpha,y), x); |
---|
[806c18] | 826 | |
---|
[6ead9d] | 827 | |
---|
| 828 | CFList LL; |
---|
| 829 | LL.append(Rg(x,z)); |
---|
| 830 | LL.append(g); |
---|
| 831 | return (LL); |
---|
| 832 | } |
---|
| 833 | ///////////////////////////////////////////////////////// |
---|
| 834 | // Compute the absolute and rational factorization of |
---|
| 835 | // the univariate polynomial 'ff^grad'. |
---|
| 836 | //======================================================= |
---|
| 837 | void BIFAC::unifac (CanonicalForm ff, int grad) |
---|
| 838 | //======================================================= |
---|
| 839 | { |
---|
| 840 | |
---|
| 841 | CFFList factorsUni; |
---|
| 842 | CFFList factorsAbs; |
---|
| 843 | CanonicalForm tmp; |
---|
| 844 | |
---|
[806c18] | 845 | factorsUni = AbsFactorize(ff); |
---|
[6ead9d] | 846 | |
---|
| 847 | for( CFFListIterator l=factorsUni; l.hasItem(); l++) |
---|
| 848 | if( ! l.getItem().factor().inBaseDomain() ) |
---|
| 849 | { |
---|
| 850 | gl_RL.append( CFFactor( l.getItem().factor(),l.getItem().exp()*grad) ); |
---|
| 851 | } |
---|
| 852 | |
---|
| 853 | |
---|
| 854 | } |
---|
| 855 | |
---|
| 856 | |
---|
| 857 | /////////////////////////////////////////////////////// |
---|
| 858 | // Compute the rational factor of f belonging to phi |
---|
| 859 | //======================================================= |
---|
| 860 | CanonicalForm BIFAC::RationalFactor (CanonicalForm phi, CanonicalForm ff, \ |
---|
[806c18] | 861 | CanonicalForm fx, CanonicalForm g) |
---|
[6ead9d] | 862 | //======================================================= |
---|
| 863 | { |
---|
| 864 | |
---|
| 865 | CanonicalForm h,hh; |
---|
| 866 | // CanonicalForm fx = deriv(f,x); |
---|
| 867 | |
---|
| 868 | for ( CFIterator it = phi; it.hasTerms(); it++ ) |
---|
| 869 | h += it.coeff() * power(fx,phi.degree()-it.exp())*power(g,it.exp()); |
---|
| 870 | |
---|
| 871 | |
---|
| 872 | hh = Bigcd(ff, h); |
---|
[806c18] | 873 | |
---|
[6ead9d] | 874 | return(hh); |
---|
| 875 | } |
---|
| 876 | //======================================================= |
---|
| 877 | void BIFAC::RationalFactorizationOnly (CFFList Phis, CanonicalForm f0, CanonicalForm g) |
---|
| 878 | //======================================================= |
---|
| 879 | { |
---|
| 880 | CanonicalForm h,ff; |
---|
| 881 | CanonicalForm fx = deriv(f0,x); |
---|
| 882 | |
---|
| 883 | for( CFFListIterator i=Phis; i.hasItem(); i++) |
---|
| 884 | { |
---|
| 885 | ASSERT( i.getItem().exp() == 1 , "Wrong factor of Eg"); // degree must be 1 |
---|
| 886 | CanonicalForm phi = i.getItem().factor(); |
---|
[806c18] | 887 | |
---|
[6ead9d] | 888 | if( ! phi.inBaseDomain()) |
---|
| 889 | { |
---|
| 890 | h = RationalFactor(phi,f0,fx,g); |
---|
| 891 | gl_RL.append( CFFactor(h,exponent )); |
---|
| 892 | ff = f0; |
---|
| 893 | f0 /= h; |
---|
| 894 | ASSERT( f0*h==ff, "Wrong factor found"); |
---|
| 895 | } |
---|
| 896 | } |
---|
| 897 | } |
---|
| 898 | |
---|
| 899 | //======================================================= |
---|
| 900 | CFList BIFAC::getAbsoluteFactors (CanonicalForm f1, CanonicalForm phi) |
---|
| 901 | //======================================================= |
---|
| 902 | { |
---|
| 903 | CanonicalForm fac; |
---|
| 904 | CanonicalForm root; |
---|
| 905 | CFList AbsFac; |
---|
| 906 | |
---|
| 907 | CFFList Fac = factorize(phi,e); |
---|
| 908 | for( CFFListIterator i=Fac; i.hasItem(); i++) |
---|
| 909 | { |
---|
| 910 | fac = i.getItem().factor(); |
---|
| 911 | if( taildegree(fac) > 0 ) // case: phi = a * x |
---|
| 912 | root = 0; |
---|
| 913 | else |
---|
| 914 | root = -tailcoeff(fac)/lc(fac); |
---|
[806c18] | 915 | |
---|
| 916 | |
---|
[6ead9d] | 917 | AbsFac.append( f1(root,e) ); |
---|
| 918 | AbsFac.append( i.getItem().exp() * exponent); |
---|
| 919 | AbsFac.append( phi ); // Polynomial of the field extension |
---|
| 920 | } |
---|
| 921 | return AbsFac; |
---|
| 922 | } |
---|
| 923 | //======================================================= |
---|
| 924 | void BIFAC::AbsoluteFactorization (CFFList Phis, CanonicalForm ff, CanonicalForm g) |
---|
| 925 | //======================================================= |
---|
| 926 | { |
---|
| 927 | |
---|
| 928 | int ii; |
---|
| 929 | if( getCharacteristic() == 0 ) |
---|
| 930 | { |
---|
| 931 | //cerr << "* Charcteristic 0 is not yet implemented! => Aborting!\n"; |
---|
| 932 | exit(1); |
---|
| 933 | } |
---|
| 934 | |
---|
| 935 | |
---|
| 936 | CFList AbsFac; |
---|
| 937 | CanonicalForm phi; |
---|
| 938 | CanonicalForm h, h_abs, h_res, h_rat; |
---|
| 939 | CanonicalForm fx = deriv(ff,x); |
---|
[806c18] | 940 | |
---|
[6ead9d] | 941 | |
---|
| 942 | for( CFFListIterator i=Phis; i.hasItem(); i++) |
---|
| 943 | { |
---|
| 944 | ASSERT( i.getItem().exp() == 1 , "Wrong factor of Eg"); // degree must be 1 |
---|
| 945 | phi = i.getItem().factor(); |
---|
[806c18] | 946 | |
---|
[6ead9d] | 947 | if( ! phi.inBaseDomain()) |
---|
| 948 | { |
---|
| 949 | |
---|
| 950 | // === Case 1: phi has degree 1 === |
---|
| 951 | if( phi.degree() == 1 ) |
---|
| 952 | { |
---|
[806c18] | 953 | if( taildegree(phi) > 0 ) // case: phi = a * x |
---|
| 954 | h = gcd( ff,g ); |
---|
| 955 | else // case: phi = a * x + c |
---|
| 956 | { |
---|
| 957 | h = gcd( ff, g+tailcoeff(phi)/lc(phi)*fx); |
---|
| 958 | } |
---|
| 959 | |
---|
| 960 | //biNormieren( h ); |
---|
| 961 | gl_AL.append(h); // Factor of degree 1 |
---|
| 962 | gl_AL.append(exponent); // Multiplicity (exponent) |
---|
| 963 | gl_AL.append(0); // No field extension |
---|
[6ead9d] | 964 | } else |
---|
| 965 | { |
---|
[806c18] | 966 | // === Case 2: phi has degree > 1 === |
---|
| 967 | e=rootOf(phi, 'e'); |
---|
| 968 | h = gcd( ff, g-e*fx); |
---|
| 969 | //biNormieren( h ); |
---|
| 970 | |
---|
| 971 | AbsFac = getAbsoluteFactors(h, phi); |
---|
| 972 | for( CFListIterator l=AbsFac; l.hasItem(); l++) |
---|
| 973 | gl_AL.append( l.getItem() ); |
---|
| 974 | |
---|
| 975 | |
---|
| 976 | // === (1) Get the rational factor by multi- === |
---|
| 977 | // === plication of the absolute factor. === |
---|
| 978 | h_abs=1; |
---|
| 979 | ii = 0; |
---|
| 980 | |
---|
| 981 | for( CFListIterator l=AbsFac; l.hasItem(); l++) |
---|
| 982 | { |
---|
| 983 | ii++; |
---|
| 984 | if (ii%3 == 1 ) |
---|
| 985 | h_abs *= l.getItem(); |
---|
| 986 | } |
---|
| 987 | //biNormieren( h_abs ); |
---|
| 988 | |
---|
| 989 | |
---|
| 990 | // === (2) Compute the rational factor === |
---|
| 991 | // === by using the resultant. === |
---|
| 992 | h_res = resultant(phi(z,x), h(z,e), z); |
---|
| 993 | //biNormieren( h_res ); |
---|
| 994 | |
---|
| 995 | |
---|
| 996 | // === (3) Compute the rational factor by ignoring === |
---|
| 997 | // === all knowledge of absolute factors. === |
---|
| 998 | h_rat = RationalFactor(phi, ff,fx, g); |
---|
| 999 | //biNormieren( h_rat ); |
---|
| 1000 | |
---|
| 1001 | ASSERT( (h_abs == h_res) && (h_res == h_rat), "Wrong rational factor ?!?"); |
---|
| 1002 | h = h_abs; |
---|
[6ead9d] | 1003 | } |
---|
| 1004 | // End of absolute factorization. |
---|
| 1005 | gl_RL.append(CFFactor( h,exponent )); // Save the rational factor |
---|
| 1006 | ff/=h; |
---|
| 1007 | } |
---|
| 1008 | } |
---|
| 1009 | } |
---|
| 1010 | |
---|
| 1011 | |
---|
[806c18] | 1012 | //====================================================== |
---|
[6ead9d] | 1013 | // Factorization of a squarefree bivariate polynomial |
---|
| 1014 | // in which every factor appears only once. |
---|
| 1015 | // Do we need a complete factorization ('absolute' is true) |
---|
| 1016 | // or only a rational factorization ('absolute' false)? |
---|
[806c18] | 1017 | //====================================================== |
---|
[6ead9d] | 1018 | void BIFAC::bifacSqrFree(CanonicalForm ff) |
---|
| 1019 | //======================================================= |
---|
| 1020 | { |
---|
| 1021 | |
---|
| 1022 | int anz=0; // number of factors without field elements |
---|
| 1023 | |
---|
| 1024 | CFList G = basisOfG(ff); |
---|
| 1025 | |
---|
| 1026 | CFList LL; |
---|
| 1027 | CanonicalForm Eg,g; |
---|
| 1028 | |
---|
| 1029 | |
---|
| 1030 | |
---|
| 1031 | // Case 1: There is only one rational & absolute factor === |
---|
| 1032 | if( G.length() == 1 ){ // There is only one |
---|
| 1033 | gl_RL.append( CFFactor(ff, exponent)); // rational factor |
---|
| 1034 | gl_AL.append( ff ); |
---|
| 1035 | gl_AL.append( exponent ); |
---|
| 1036 | gl_AL.append( 0 ); |
---|
| 1037 | } |
---|
| 1038 | else // Case 2: There is more than one absolute factor === |
---|
| 1039 | { |
---|
| 1040 | // LL = createEg(G,ff); |
---|
| 1041 | // LL = createEgUni(G,ff); // Hier ist noch ein FEHLER !!!! |
---|
[806c18] | 1042 | |
---|
[6ead9d] | 1043 | LL = createRg( G, ff); // viel langsamer als EgUni |
---|
[806c18] | 1044 | |
---|
| 1045 | |
---|
[6ead9d] | 1046 | Eg = LL.getFirst(); |
---|
[806c18] | 1047 | Eg = Eg/LC(Eg); |
---|
| 1048 | |
---|
[6ead9d] | 1049 | g = LL.getLast(); |
---|
| 1050 | |
---|
| 1051 | // g = G.getFirst(); |
---|
| 1052 | |
---|
| 1053 | |
---|
[806c18] | 1054 | CFFList PHI = AbsFactorize( Eg ); |
---|
| 1055 | |
---|
| 1056 | CFFListIterator J=PHI; |
---|
| 1057 | CanonicalForm Eg2=1; |
---|
| 1058 | for ( ; J.hasItem(); J++) |
---|
| 1059 | { Eg2 = Eg2 * J.getItem().factor(); } |
---|
[6ead9d] | 1060 | |
---|
| 1061 | // === Is Eg(x) irreducible ? === |
---|
| 1062 | anz=0; |
---|
| 1063 | |
---|
[806c18] | 1064 | // PHI = AbsFactorize( Eg) ; |
---|
| 1065 | // |
---|
| 1066 | |
---|
| 1067 | for( CFFListIterator i=PHI; i.hasItem(); i++) { |
---|
[6ead9d] | 1068 | if( !i.getItem().factor().inBaseDomain()) |
---|
[806c18] | 1069 | anz++; |
---|
| 1070 | } |
---|
[6ead9d] | 1071 | |
---|
| 1072 | /* if( absolute ) // Only for a absolute factorization |
---|
| 1073 | AbsoluteFactorization( PHI,ff, g); |
---|
[806c18] | 1074 | else // only for a rational factorization |
---|
[6ead9d] | 1075 | { */ |
---|
| 1076 | if( anz==1 ){ ; |
---|
[806c18] | 1077 | gl_RL.append( CFFactor(ff,exponent));} |
---|
| 1078 | else |
---|
| 1079 | RationalFactorizationOnly( PHI,ff, g); |
---|
[6ead9d] | 1080 | /* } */ |
---|
| 1081 | } |
---|
| 1082 | } |
---|
| 1083 | |
---|
| 1084 | ///////////////////////////////////////////// |
---|
| 1085 | // Main procedure for the factorization |
---|
| 1086 | // of the bivariate polynomial 'f'. |
---|
| 1087 | // REMARK: 'f' might be univariate, too. |
---|
| 1088 | //--<>--------------------------------- |
---|
| 1089 | void BIFAC::bifacMain(CanonicalForm f) |
---|
| 1090 | //--<>--------------------------------- |
---|
| 1091 | { |
---|
| 1092 | |
---|
| 1093 | |
---|
| 1094 | CanonicalForm ff, ggT; |
---|
| 1095 | |
---|
| 1096 | // =============================================== |
---|
| 1097 | // = (1) Trivial case: Input is a constant = |
---|
| 1098 | // =============================================== |
---|
| 1099 | |
---|
| 1100 | if( f.inBaseDomain() ) |
---|
| 1101 | { |
---|
| 1102 | gl_AL.append(f); // store polynomial |
---|
| 1103 | gl_AL.append(1); // store exponent |
---|
| 1104 | gl_AL.append(0); // store ŽpolynomialŽ for field extension |
---|
| 1105 | |
---|
| 1106 | gl_RL.append( CFFactor(f,1) ); // store polynomial |
---|
| 1107 | return; |
---|
| 1108 | } |
---|
| 1109 | |
---|
| 1110 | // =============================================== |
---|
| 1111 | // = STEP: Squarefree decomposition = |
---|
| 1112 | // =============================================== |
---|
| 1113 | |
---|
[806c18] | 1114 | |
---|
| 1115 | CFFList Q =Mysqrfree(f); |
---|
| 1116 | // |
---|
| 1117 | // cout << Q << endl; |
---|
[6ead9d] | 1118 | // |
---|
| 1119 | |
---|
| 1120 | |
---|
| 1121 | |
---|
| 1122 | // ========================================================= |
---|
| 1123 | // = STEP: Factorization of the squarefree decomposition = |
---|
| 1124 | // ========================================================= |
---|
| 1125 | |
---|
| 1126 | |
---|
| 1127 | for( CFFListIterator i=Q; i.hasItem(); i++) |
---|
| 1128 | { |
---|
| 1129 | |
---|
[806c18] | 1130 | if( i.getItem().factor().level() < 0 ) ; |
---|
| 1131 | else |
---|
| 1132 | { |
---|
[6ead9d] | 1133 | if( ( degree(i.getItem().factor(),x) == 0 || degree( i.getItem().factor(),y) == 0) ) { |
---|
| 1134 | // case: univariate |
---|
| 1135 | unifac( i.getItem().factor(), i.getItem().exp() ); } |
---|
| 1136 | else // case: bivariate |
---|
| 1137 | { |
---|
[806c18] | 1138 | exponent = i.getItem().exp(); // global variable |
---|
| 1139 | CanonicalForm dumm = i.getItem().factor(); |
---|
| 1140 | dumm = dumm.LC(); |
---|
| 1141 | if( dumm.level() > 0 ){ dumm = 1; } |
---|
| 1142 | bifacSqrFree(i.getItem().factor()/dumm ); |
---|
[6ead9d] | 1143 | } |
---|
| 1144 | }} |
---|
| 1145 | |
---|
| 1146 | |
---|
| 1147 | } |
---|
| 1148 | |
---|
| 1149 | |
---|
| 1150 | /////////////////////////////////////////////////////// |
---|
| 1151 | // Find the least prime so that the factorization |
---|
| 1152 | // works. |
---|
| 1153 | /////////////////////////////////////////////////////// |
---|
| 1154 | |
---|
| 1155 | //======================================================= |
---|
| 1156 | int BIFAC::findCharacteristic(CanonicalForm f) |
---|
| 1157 | //======================================================= |
---|
| 1158 | { |
---|
| 1159 | int min = (2*degree(f,'x')-1)*degree(f,'y'); |
---|
| 1160 | int nr=0; |
---|
| 1161 | |
---|
| 1162 | if( min >= 32003 ) return ( 32003 ); // this is the maximum |
---|
[806c18] | 1163 | |
---|
[6ead9d] | 1164 | // Find the smallest poosible prime |
---|
| 1165 | while ( cf_getPrime(nr) < min) { nr++; } |
---|
| 1166 | return ( cf_getPrime(nr) ); |
---|
| 1167 | } |
---|
| 1168 | |
---|
| 1169 | ///////////////////////////////////////////////////////// |
---|
| 1170 | // |
---|
| 1171 | // PUBLIC functions |
---|
| 1172 | // |
---|
| 1173 | ///////////////////////////////////////////////////////// |
---|
| 1174 | |
---|
| 1175 | // convert the result of the factorization from |
---|
| 1176 | // the intern storage type into the public one. |
---|
| 1177 | // Also, check the correctness of the solution |
---|
| 1178 | // and, if neccessary, change the characteristic. |
---|
| 1179 | //--<>--------------------------------- |
---|
| 1180 | void BIFAC::convertResult(CanonicalForm & f, int ch, int sw) |
---|
| 1181 | //--<>--------------------------------- |
---|
| 1182 | { |
---|
| 1183 | |
---|
| 1184 | CanonicalForm ff = 1; |
---|
| 1185 | CanonicalForm c; |
---|
| 1186 | |
---|
| 1187 | CFFList aL; |
---|
| 1188 | |
---|
| 1189 | //cout << gl_RL<<endl; |
---|
| 1190 | |
---|
[806c18] | 1191 | if( sw ) |
---|
| 1192 | { |
---|
| 1193 | Variable W('W'); |
---|
| 1194 | for( CFFListIterator i=gl_RL; i.hasItem(); i++) |
---|
| 1195 | { |
---|
| 1196 | c = i.getItem().factor(); |
---|
| 1197 | c = c(W,y); |
---|
| 1198 | c = c(y,x); |
---|
| 1199 | c = c(x,W); |
---|
| 1200 | aL.append( CFFactor( c, i.getItem().exp() )); |
---|
| 1201 | } |
---|
[6ead9d] | 1202 | |
---|
[806c18] | 1203 | f = f(W,y); f=f(y,x); f=f(x,W); |
---|
| 1204 | } |
---|
| 1205 | else aL = gl_RL; |
---|
[6ead9d] | 1206 | |
---|
[806c18] | 1207 | gl_RL = aL; |
---|
[6ead9d] | 1208 | |
---|
[806c18] | 1209 | //cout << aL; |
---|
[6ead9d] | 1210 | |
---|
| 1211 | |
---|
| 1212 | |
---|
| 1213 | // ========== OUTPUT ===================== |
---|
| 1214 | /* for( CFFListIterator i=aL; i.hasItem(); i++) |
---|
| 1215 | { |
---|
| 1216 | cout << "(" << i.getItem().factor() << ")"; |
---|
| 1217 | if( i.getItem().exp() != 1 ) |
---|
| 1218 | cout << "^" << i.getItem().exp(); |
---|
| 1219 | cout << " * "; |
---|
| 1220 | } */ |
---|
| 1221 | |
---|
| 1222 | |
---|
| 1223 | // cout << "\n* Test auf Korrektheit ..."; |
---|
[806c18] | 1224 | |
---|
| 1225 | |
---|
[6ead9d] | 1226 | for( CFFListIterator i=aL; i.hasItem(); i++) |
---|
| 1227 | { |
---|
| 1228 | ff *= power(i.getItem().factor(), i.getItem().exp() ); |
---|
| 1229 | // cout << " ff = " << ff |
---|
[806c18] | 1230 | // << "\n a^b = " << i.getItem().factor() << " ^ " << i.getItem().exp() << endl; |
---|
[6ead9d] | 1231 | } |
---|
| 1232 | c = f.LC()/ff.LC(); |
---|
| 1233 | |
---|
| 1234 | ff *= c; |
---|
| 1235 | |
---|
| 1236 | |
---|
| 1237 | // cout << "\n\nOriginal f = " << f << "\n\nff = " << ff |
---|
[806c18] | 1238 | // << "\n\nDiff = " << f-ff << endl << "Quot "<< f/ff <<endl; |
---|
[6ead9d] | 1239 | // cout << "degree 0: " << c << endl; |
---|
[806c18] | 1240 | |
---|
| 1241 | |
---|
[6ead9d] | 1242 | #ifndef NOSTREAMIO |
---|
[806c18] | 1243 | if( f != ff ) cout << "\n\nOriginal f = " << f << "\n\nff = " << ff |
---|
| 1244 | << "\n\nDiff = " << f-ff << endl << "Quot "<< f/ff <<endl; |
---|
[6ead9d] | 1245 | #endif |
---|
| 1246 | ASSERT( f==ff, "Wrong rational factorization. Abborting!"); |
---|
| 1247 | // cout << " [OK]\n"; |
---|
[806c18] | 1248 | |
---|
[6ead9d] | 1249 | } |
---|
| 1250 | //--<>--------------------------------- |
---|
| 1251 | void BIFAC::bifac(CanonicalForm f, bool abs) |
---|
| 1252 | //--<>--------------------------------- |
---|
| 1253 | { |
---|
| 1254 | absolute = 1; // global variables |
---|
| 1255 | CFList factors; |
---|
| 1256 | int ch = getCharacteristic(); |
---|
| 1257 | int ch2; |
---|
| 1258 | |
---|
| 1259 | |
---|
| 1260 | ASSERT( ch==0 && !isOn(SW_RATIONAL), "Integer numbers not allowed" ); |
---|
[806c18] | 1261 | |
---|
[6ead9d] | 1262 | |
---|
| 1263 | // === Check the characteristic === |
---|
[806c18] | 1264 | if( ch != 0 ) |
---|
[6ead9d] | 1265 | { |
---|
| 1266 | ch2 = findCharacteristic(f); |
---|
| 1267 | if( ch < ch2 ) |
---|
| 1268 | { |
---|
| 1269 | // setCharacteristic( ch2 ); |
---|
| 1270 | f = mapinto(f); |
---|
| 1271 | |
---|
| 1272 | // PROVISORISCH |
---|
| 1273 | //cerr << "\n Characteristic is too small!" |
---|
[806c18] | 1274 | // << "\n The result might be wrong!\n\n"; |
---|
[6ead9d] | 1275 | exit(1); |
---|
| 1276 | |
---|
| 1277 | } else ; |
---|
| 1278 | } |
---|
| 1279 | |
---|
[806c18] | 1280 | Variable W('W'); |
---|
| 1281 | CanonicalForm l; |
---|
| 1282 | int sw = 0; |
---|
[6ead9d] | 1283 | |
---|
[806c18] | 1284 | if( degree(f,x) < degree(f,y) ) { |
---|
| 1285 | f = f(W,x); f = f(x,y); f=f(y,W); |
---|
| 1286 | sw = 1; |
---|
| 1287 | } |
---|
| 1288 | l = f.LC(); |
---|
[6ead9d] | 1289 | |
---|
[806c18] | 1290 | if( l.level()<0 ) { f = f/f.LC(); gl_RL.append( CFFactor(l,1) ); } |
---|
[6ead9d] | 1291 | |
---|
| 1292 | |
---|
| 1293 | bifacMain(f); // start the computation |
---|
| 1294 | |
---|
| 1295 | convertResult(f,ch, sw) ; // and convert the result |
---|
| 1296 | } |
---|
| 1297 | |
---|
| 1298 | // ============== end of 'bifac.cc' ================== |
---|
[e4fe2b] | 1299 | #endif /* #ifdef HAVE_BIFAC */ |
---|