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