[6ead9d] | 1 | |
---|
[9c3d69] | 2 | #include "canonicalform.h" |
---|
[0c4a34b] | 3 | #ifdef HAVE_BIFAC |
---|
[6ead9d] | 4 | #ifndef NOSTREAMIO |
---|
| 5 | #include<fstream> |
---|
| 6 | #endif |
---|
| 7 | #include<math.h> |
---|
| 8 | #include <sys/timeb.h> |
---|
| 9 | |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | static void |
---|
| 13 | fillVarsRec ( const CanonicalForm & f, int * vars ) |
---|
| 14 | { |
---|
| 15 | int n; |
---|
| 16 | if ( (n = f.level()) > 0 ) |
---|
| 17 | { |
---|
| 18 | vars[n] = 1; |
---|
| 19 | CFIterator i; |
---|
| 20 | for ( i = f; i.hasTerms(); ++i ) |
---|
| 21 | fillVarsRec( i.coeff(), vars ); |
---|
| 22 | } |
---|
| 23 | } |
---|
| 24 | |
---|
| 25 | int ExtensionLevel(); |
---|
| 26 | void Reduce( bool); |
---|
| 27 | |
---|
| 28 | CanonicalForm MYGCD( const CanonicalForm& f, const CanonicalForm& g); |
---|
| 29 | |
---|
| 30 | |
---|
| 31 | CanonicalForm MyContent( const CanonicalForm& F, const Variable& x) |
---|
| 32 | { |
---|
| 33 | CanonicalForm r,t; |
---|
| 34 | CanonicalForm g=F; |
---|
| 35 | CanonicalForm one=1; |
---|
| 36 | |
---|
| 37 | if( F.isZero() ) return 0; |
---|
[806c18] | 38 | if( F.inBaseDomain() ) return F; |
---|
[6ead9d] | 39 | |
---|
[806c18] | 40 | if( level(F) < 0 ) return 1; |
---|
[6ead9d] | 41 | |
---|
| 42 | r = LC( F, x); |
---|
| 43 | |
---|
| 44 | g = g - power(x,degree(g,x))*r; |
---|
| 45 | |
---|
| 46 | while( g.isZero() != 1 && r!= 1 && r!=-1 ) |
---|
| 47 | { |
---|
| 48 | t = LC(g, x); |
---|
[806c18] | 49 | if( t == 1 || t == -1 ) return 1; |
---|
[6ead9d] | 50 | r = MYGCD( r, t); |
---|
[806c18] | 51 | if( r == 1 ) return 1; |
---|
| 52 | g = g - power(x,degree(g,x))*t; |
---|
[6ead9d] | 53 | } |
---|
| 54 | return r; |
---|
| 55 | } |
---|
| 56 | |
---|
| 57 | void CurrentExtension() |
---|
| 58 | { |
---|
| 59 | Variable x('x'); |
---|
| 60 | int i; |
---|
| 61 | #ifndef NOSTREAMIO |
---|
| 62 | cout << "Current Extension: "<<endl; |
---|
| 63 | #endif |
---|
| 64 | for (i = ExtensionLevel();i>0;i--) |
---|
| 65 | { |
---|
| 66 | Variable l(-i); |
---|
| 67 | #ifndef NOSTREAMIO |
---|
| 68 | cout << "Variable: "<<l<<" Level: "<<l.level()<<" Minimal Polynom: "<<getMipo(l,x)<<endl; |
---|
| 69 | #endif |
---|
| 70 | } |
---|
| 71 | } |
---|
| 72 | |
---|
| 73 | /*Liefert den ggt aller numerischen Koeffizienten einer Canonischen Form */ |
---|
| 74 | |
---|
| 75 | CanonicalForm MyNum(const CanonicalForm & a) |
---|
| 76 | { |
---|
| 77 | bool bruch = isOn(SW_RATIONAL); |
---|
| 78 | Off (SW_RATIONAL); |
---|
| 79 | |
---|
| 80 | CanonicalForm dummy =0; |
---|
| 81 | CanonicalForm dummy2; |
---|
| 82 | |
---|
| 83 | CFIterator F =a; |
---|
| 84 | |
---|
| 85 | for ( ; F.hasTerms(); F++) |
---|
| 86 | { |
---|
| 87 | if (F.coeff().inBaseDomain()) |
---|
| 88 | { |
---|
| 89 | dummy2 = F.coeff().num(); |
---|
| 90 | if (dummy == 0) |
---|
| 91 | { |
---|
| 92 | dummy = dummy2; |
---|
| 93 | } |
---|
| 94 | else |
---|
| 95 | { |
---|
| 96 | dummy = gcd(dummy, dummy2); |
---|
| 97 | } |
---|
| 98 | } |
---|
| 99 | else |
---|
| 100 | { |
---|
| 101 | dummy2 = MyNum(F.coeff()); |
---|
| 102 | if (dummy == 0) |
---|
| 103 | { |
---|
| 104 | dummy = dummy2; |
---|
| 105 | } |
---|
| 106 | else |
---|
| 107 | { |
---|
| 108 | dummy = gcd(dummy, dummy2); |
---|
| 109 | } |
---|
| 110 | } |
---|
| 111 | } |
---|
| 112 | if (bruch) |
---|
| 113 | On (SW_RATIONAL); |
---|
| 114 | else |
---|
| 115 | Off(SW_RATIONAL); |
---|
| 116 | return dummy; |
---|
| 117 | } |
---|
| 118 | |
---|
| 119 | /* Liefert den kgV aller Nenner der Koeffizenten einer Canonischen Form */ |
---|
| 120 | |
---|
| 121 | CanonicalForm MyDen(const CanonicalForm & a) |
---|
| 122 | { |
---|
| 123 | bool bruch = isOn(SW_RATIONAL); |
---|
| 124 | Off (SW_RATIONAL); |
---|
| 125 | |
---|
| 126 | CanonicalForm dummy(1); |
---|
| 127 | CanonicalForm dummy2; |
---|
| 128 | |
---|
| 129 | CFIterator F =a; |
---|
| 130 | |
---|
| 131 | for ( ; F.hasTerms(); F++) |
---|
| 132 | { |
---|
| 133 | if (F.coeff().inBaseDomain()) |
---|
| 134 | { |
---|
| 135 | dummy2 = gcd(dummy,F.coeff().den()); |
---|
| 136 | dummy = dummy * F.coeff().den()/dummy2; |
---|
| 137 | } |
---|
| 138 | else |
---|
| 139 | { |
---|
| 140 | dummy2 = MyDen(F.coeff()); |
---|
| 141 | dummy = dummy*dummy2/gcd(dummy,dummy2); |
---|
| 142 | } |
---|
| 143 | } |
---|
| 144 | if (bruch) |
---|
| 145 | On (SW_RATIONAL); |
---|
| 146 | else |
---|
| 147 | Off(SW_RATIONAL); |
---|
| 148 | return dummy; |
---|
| 149 | } |
---|
| 150 | |
---|
| 151 | /* Liefert die normierte Canonische Form a zurück, wenn LC(a) kein Nullteiler in Characteristic p ist */ |
---|
| 152 | /* sonst -1*/ |
---|
| 153 | CanonicalForm MyMonic(const CanonicalForm & a, const CanonicalForm & r, const int & l) |
---|
| 154 | { |
---|
| 155 | bool bruch = isOn(SW_RATIONAL); |
---|
| 156 | int zaehler; |
---|
| 157 | int Level = l; |
---|
| 158 | CanonicalForm dummy, dummy1, dummy2; |
---|
| 159 | CanonicalForm g = a; |
---|
| 160 | CanonicalForm p = r; |
---|
| 161 | |
---|
| 162 | On (SW_RATIONAL); |
---|
| 163 | |
---|
| 164 | if (Level == g.level()) |
---|
| 165 | { |
---|
| 166 | dummy = 1/g.LC(); |
---|
| 167 | } |
---|
| 168 | else |
---|
| 169 | { |
---|
| 170 | dummy = 1/g; |
---|
| 171 | } |
---|
| 172 | dummy1 = MyDen(dummy); |
---|
| 173 | dummy2 =dummy1; |
---|
| 174 | zaehler =1; |
---|
| 175 | |
---|
| 176 | Off (SW_RATIONAL); |
---|
| 177 | |
---|
| 178 | while ((mod(dummy2,p) != 1) && (mod(dummy2,p) !=0)) |
---|
| 179 | { |
---|
| 180 | dummy2 = dummy2*dummy1; |
---|
| 181 | dummy2 =mod(dummy2,p); |
---|
| 182 | zaehler++; |
---|
| 183 | } |
---|
| 184 | if (mod(dummy2,p).isZero()) |
---|
| 185 | { |
---|
| 186 | if (bruch) |
---|
| 187 | On (SW_RATIONAL); |
---|
| 188 | else |
---|
| 189 | Off(SW_RATIONAL); |
---|
| 190 | return -1; |
---|
| 191 | } |
---|
| 192 | else |
---|
| 193 | { |
---|
| 194 | zaehler--; |
---|
| 195 | dummy2 = power(dummy1,zaehler); |
---|
| 196 | dummy2 = mod(dummy2,p); |
---|
| 197 | dummy*= dummy1; |
---|
| 198 | dummy*= dummy2; |
---|
| 199 | g =mod(g*dummy,p); |
---|
| 200 | if (bruch) |
---|
| 201 | On (SW_RATIONAL); |
---|
| 202 | else |
---|
| 203 | Off(SW_RATIONAL); |
---|
| 204 | return g; |
---|
| 205 | } |
---|
| 206 | } |
---|
| 207 | |
---|
| 208 | |
---|
| 209 | |
---|
| 210 | /*Berechnet den ggT der Formen a und b in Characteristic p*/ |
---|
| 211 | |
---|
| 212 | CanonicalForm MyGCDlocal (const CanonicalForm & a, const CanonicalForm & b, const CanonicalForm & r, const int & l) |
---|
| 213 | { |
---|
| 214 | bool bruch = isOn(SW_RATIONAL); |
---|
| 215 | Off(SW_RATIONAL); |
---|
| 216 | CanonicalForm Rest, Result; |
---|
| 217 | |
---|
| 218 | CanonicalForm f=a; |
---|
| 219 | CanonicalForm g=b; |
---|
| 220 | CanonicalForm p=r; |
---|
| 221 | int Level =l; |
---|
| 222 | f = mod(f,p); |
---|
| 223 | g = mod(g,p); |
---|
| 224 | |
---|
| 225 | Rest=g; |
---|
| 226 | |
---|
| 227 | while (!Rest.isZero()) |
---|
| 228 | { |
---|
| 229 | g = MyMonic(g,p,Level); |
---|
| 230 | |
---|
| 231 | if (g == -1) |
---|
| 232 | { |
---|
| 233 | if (bruch) |
---|
| 234 | On(SW_RATIONAL); |
---|
| 235 | else |
---|
| 236 | Off(SW_RATIONAL); |
---|
| 237 | return -1; |
---|
| 238 | } |
---|
| 239 | else |
---|
| 240 | { |
---|
| 241 | Result =g; |
---|
| 242 | Rest = f%g; |
---|
| 243 | f = g; |
---|
| 244 | g = Rest; |
---|
| 245 | f =mod(f,p); |
---|
| 246 | g =mod(g,p); |
---|
| 247 | Rest =mod(Rest,p); |
---|
| 248 | } |
---|
| 249 | } |
---|
| 250 | if (bruch) |
---|
| 251 | On(SW_RATIONAL); |
---|
| 252 | else |
---|
| 253 | Off(SW_RATIONAL); |
---|
| 254 | return Result; |
---|
| 255 | } |
---|
| 256 | |
---|
| 257 | /* Chinese Remaindering für a mod m und b mod l */ |
---|
| 258 | |
---|
| 259 | CanonicalForm MyChiRem(const CanonicalForm & a,const CanonicalForm & m,const CanonicalForm & b,const CanonicalForm & l) |
---|
| 260 | { |
---|
| 261 | bool bruch = isOn(SW_RATIONAL); |
---|
| 262 | |
---|
| 263 | CanonicalForm u,v,Runner; |
---|
| 264 | CanonicalForm Result(0); |
---|
| 265 | CanonicalForm LeadTerm; |
---|
| 266 | CanonicalForm x1=a; |
---|
| 267 | CanonicalForm m1 = m; |
---|
| 268 | CanonicalForm x2=b; |
---|
| 269 | CanonicalForm m2 = l; |
---|
| 270 | |
---|
| 271 | while (!x1.isZero() || !x2.isZero()) |
---|
| 272 | { |
---|
| 273 | if (x1.degree() > x2.degree()) |
---|
| 274 | { |
---|
| 275 | LeadTerm = power(x1.mvar(),x1.degree()); |
---|
| 276 | u = x1.LC()*LeadTerm; |
---|
| 277 | v = 0; |
---|
| 278 | x1 = x1-u; |
---|
| 279 | } |
---|
| 280 | else |
---|
| 281 | { |
---|
| 282 | if (x1.degree() < x2.degree()) |
---|
| 283 | { |
---|
| 284 | u = 0; |
---|
| 285 | LeadTerm = power(x2.mvar(),x2.degree()); |
---|
| 286 | v = x2.LC()*LeadTerm; |
---|
| 287 | x2 = x2-v; |
---|
| 288 | } |
---|
| 289 | else |
---|
| 290 | { |
---|
| 291 | if (x1.degree() == x2.degree()) |
---|
| 292 | { |
---|
| 293 | LeadTerm = power(x2.mvar(),x2.degree()); |
---|
| 294 | u = x1.LC()*LeadTerm; |
---|
| 295 | v = x2.LC()*LeadTerm; |
---|
| 296 | x1 = x1-u; |
---|
| 297 | x2 = x2-v; |
---|
| 298 | } |
---|
| 299 | } |
---|
| 300 | } |
---|
| 301 | |
---|
| 302 | if (u.LC().inBaseDomain() && v.LC().inBaseDomain()) |
---|
| 303 | { |
---|
| 304 | Runner = u.LC(); |
---|
| 305 | Off(SW_RATIONAL); |
---|
| 306 | while(mod(Runner,m2) !=v.LC()) |
---|
| 307 | { |
---|
| 308 | Runner = Runner+m1; |
---|
| 309 | } |
---|
| 310 | |
---|
| 311 | Result = Result+Runner*LeadTerm; |
---|
| 312 | |
---|
| 313 | } |
---|
| 314 | else |
---|
| 315 | { |
---|
| 316 | Result = Result + MyChiRem(u.LC(),m1, v.LC(), m2)*LeadTerm; |
---|
| 317 | } |
---|
| 318 | } |
---|
| 319 | if (bruch) |
---|
| 320 | On(SW_RATIONAL); |
---|
| 321 | else |
---|
| 322 | Off(SW_RATIONAL); |
---|
| 323 | return Result; |
---|
| 324 | } |
---|
| 325 | |
---|
| 326 | /*Rational Rekonstruction für a mod b*/ |
---|
| 327 | |
---|
| 328 | CanonicalForm MyRatRed(const CanonicalForm & a,const CanonicalForm & b) |
---|
| 329 | { |
---|
| 330 | bool bruch = isOn(SW_RATIONAL); |
---|
| 331 | |
---|
| 332 | CanonicalForm f,dummy,dummy1,dummy2, Wurzel; |
---|
| 333 | CanonicalForm q,u0,u1,v0,v1; |
---|
| 334 | CanonicalForm Result(0); |
---|
| 335 | |
---|
| 336 | CFIterator F =a; |
---|
| 337 | |
---|
| 338 | for ( ; F.hasTerms(); F++) |
---|
| 339 | { |
---|
| 340 | if (F.coeff().inBaseDomain()) |
---|
| 341 | { |
---|
| 342 | Wurzel = sqrt(b); |
---|
| 343 | u0 =b; |
---|
| 344 | u1= F.coeff(); |
---|
| 345 | v1 = 1; |
---|
| 346 | v0 = 0; |
---|
| 347 | |
---|
| 348 | int i=0; |
---|
| 349 | |
---|
| 350 | while(!(u1<Wurzel)) |
---|
| 351 | { |
---|
| 352 | Off(SW_RATIONAL); |
---|
| 353 | q=u0/u1; |
---|
| 354 | dummy = u0-q*u1; |
---|
| 355 | u0=u1; |
---|
| 356 | u1=dummy; |
---|
| 357 | dummy = v0+q*v1; |
---|
| 358 | v0=v1; |
---|
| 359 | v1=dummy; |
---|
| 360 | i++; |
---|
| 361 | } |
---|
| 362 | f = -1; |
---|
| 363 | |
---|
| 364 | On(SW_RATIONAL); |
---|
| 365 | dummy2 = u1/v1; |
---|
| 366 | |
---|
| 367 | f = power(f,i)*dummy2; |
---|
| 368 | |
---|
| 369 | if (f.isZero()) |
---|
| 370 | { |
---|
| 371 | Off(SW_RATIONAL); |
---|
| 372 | |
---|
| 373 | if (!mod(F.coeff(),b).isZero()) |
---|
| 374 | { |
---|
| 375 | if (bruch) |
---|
| 376 | On(SW_RATIONAL); |
---|
| 377 | else |
---|
| 378 | Off(SW_RATIONAL); |
---|
| 379 | return -1; |
---|
| 380 | } |
---|
| 381 | |
---|
| 382 | } |
---|
| 383 | Result = Result+f*power(a.mvar(),F.exp()); |
---|
| 384 | |
---|
| 385 | } |
---|
| 386 | else |
---|
| 387 | { |
---|
| 388 | dummy1 = MyRatRed(F.coeff(),b); |
---|
| 389 | if (dummy1 == -1) |
---|
| 390 | { |
---|
| 391 | if (bruch) |
---|
| 392 | On(SW_RATIONAL); |
---|
| 393 | else |
---|
| 394 | Off(SW_RATIONAL); |
---|
| 395 | return -1; |
---|
| 396 | } |
---|
| 397 | else |
---|
| 398 | Result = Result + dummy1*power(a.mvar(),F.exp()); |
---|
| 399 | } |
---|
| 400 | } |
---|
| 401 | if (bruch) |
---|
| 402 | On(SW_RATIONAL); |
---|
| 403 | else |
---|
| 404 | Off(SW_RATIONAL); |
---|
| 405 | return Result; |
---|
| 406 | } |
---|
| 407 | |
---|
| 408 | /*Berechnet lokale ggT's der Formen a und b und liftet sie wieder hoch*/ |
---|
| 409 | |
---|
| 410 | CanonicalForm MyGCDmod( const CanonicalForm & a,const CanonicalForm & b) |
---|
| 411 | { |
---|
| 412 | bool bruch = isOn(SW_RATIONAL); |
---|
| 413 | // cout << "enter MyGCD mit a= "<<a<<endl; |
---|
| 414 | // cout << "und b= "<<b<<endl; |
---|
| 415 | CanonicalForm LeadA, LeadB; |
---|
| 416 | CanonicalForm Kandidat,Kandidat1; |
---|
| 417 | CanonicalForm f,g, Result; |
---|
| 418 | CanonicalForm NennerA =1; |
---|
| 419 | CanonicalForm NennerB=1; |
---|
| 420 | CanonicalForm ZahlerA=1; |
---|
| 421 | CanonicalForm ZahlerB=1; |
---|
| 422 | int treffer = 0; |
---|
| 423 | int Level; |
---|
| 424 | CanonicalForm Modulo; |
---|
| 425 | int i = 0; |
---|
| 426 | bool TryAgain = 1; |
---|
| 427 | int Primes[1228]; |
---|
| 428 | |
---|
| 429 | |
---|
| 430 | |
---|
| 431 | for (int i = 0;i <1228;i++) |
---|
| 432 | { |
---|
| 433 | Primes[i]=cf_getPrime(i+1); |
---|
| 434 | |
---|
| 435 | } |
---|
| 436 | Level=a.level(); |
---|
| 437 | |
---|
| 438 | if (a.degree() > b.degree()) |
---|
| 439 | { |
---|
| 440 | f = a; |
---|
| 441 | g = b; |
---|
| 442 | } |
---|
| 443 | else |
---|
| 444 | { |
---|
| 445 | g = a; |
---|
| 446 | f = b; |
---|
| 447 | } |
---|
| 448 | |
---|
| 449 | if (g.isZero()) |
---|
| 450 | { |
---|
| 451 | if (f.isZero()) return CanonicalForm(1); |
---|
| 452 | return f/f.LC(); |
---|
| 453 | } |
---|
| 454 | |
---|
| 455 | NennerA = MyDen(f); |
---|
| 456 | NennerB = MyDen(g); |
---|
| 457 | |
---|
| 458 | f = f*NennerA; |
---|
| 459 | g = g*NennerB; |
---|
| 460 | |
---|
| 461 | ZahlerA = MyNum(f); |
---|
| 462 | ZahlerB = MyNum(g); |
---|
| 463 | |
---|
| 464 | f=f/ZahlerA; |
---|
| 465 | g=g/ZahlerB; |
---|
| 466 | |
---|
| 467 | LeadA = f.LC(); |
---|
| 468 | while (!LeadA.inBaseDomain()) |
---|
| 469 | { |
---|
| 470 | LeadA =LeadA.LC(); |
---|
| 471 | } |
---|
| 472 | LeadB = g.LC(); |
---|
| 473 | while (!LeadB.inBaseDomain()) |
---|
| 474 | { |
---|
| 475 | LeadB =LeadB.LC(); |
---|
| 476 | } |
---|
| 477 | |
---|
| 478 | Off (SW_RATIONAL); |
---|
| 479 | |
---|
| 480 | |
---|
| 481 | while (TryAgain && i < 1228) |
---|
| 482 | { |
---|
| 483 | |
---|
| 484 | CanonicalForm p(Primes[i]); |
---|
| 485 | // cout << "p: "<<p<<endl; |
---|
| 486 | i++; |
---|
| 487 | if ( (mod(LeadA,p) != 0) && (mod(LeadB,p) != 0)) |
---|
| 488 | { |
---|
| 489 | Result = MyGCDlocal(f,g,p,Level); |
---|
| 490 | |
---|
| 491 | if (Result !=-1) |
---|
| 492 | { |
---|
| 493 | |
---|
| 494 | if (Result == 1) |
---|
| 495 | { |
---|
| 496 | if (bruch) |
---|
| 497 | On(SW_RATIONAL); |
---|
| 498 | else |
---|
| 499 | Off(SW_RATIONAL); |
---|
| 500 | return Result; |
---|
| 501 | } |
---|
| 502 | else |
---|
| 503 | { |
---|
| 504 | if (treffer == 0 || Kandidat.degree() > Result.degree()) |
---|
| 505 | { |
---|
| 506 | treffer = 1; |
---|
| 507 | Kandidat = Result; |
---|
| 508 | Modulo = p; |
---|
| 509 | } |
---|
| 510 | else |
---|
| 511 | { |
---|
| 512 | if (Kandidat.degree() == Result.degree()) |
---|
| 513 | { |
---|
| 514 | Kandidat = MyChiRem(Kandidat,Modulo,Result,p); |
---|
| 515 | Modulo = Modulo*p; |
---|
| 516 | treffer++; |
---|
| 517 | } |
---|
| 518 | } |
---|
| 519 | if (mod(treffer,4) ==1) |
---|
| 520 | { |
---|
| 521 | Kandidat1=MyRatRed(Kandidat, Modulo); |
---|
| 522 | |
---|
| 523 | |
---|
| 524 | if (Kandidat1 !=-1) |
---|
| 525 | { |
---|
| 526 | Off(SW_RATIONAL); |
---|
| 527 | if (mod(f,Kandidat1) == 0 && mod(g,Kandidat1) == 0) |
---|
| 528 | { |
---|
| 529 | break; |
---|
| 530 | } |
---|
| 531 | |
---|
| 532 | } |
---|
| 533 | } |
---|
| 534 | |
---|
| 535 | |
---|
| 536 | } |
---|
| 537 | } |
---|
| 538 | } |
---|
| 539 | else |
---|
| 540 | { |
---|
| 541 | |
---|
| 542 | } |
---|
| 543 | |
---|
| 544 | } |
---|
| 545 | |
---|
| 546 | if (bruch) |
---|
| 547 | On(SW_RATIONAL); |
---|
| 548 | else |
---|
| 549 | Off(SW_RATIONAL); |
---|
| 550 | return Kandidat1; |
---|
| 551 | } |
---|
| 552 | |
---|
| 553 | /* Berechnet die Norm eines Form h über zwei Körpererweiterungen und faktorisiert sie*/ |
---|
| 554 | |
---|
| 555 | CFFList FactorizeNorm (const CanonicalForm & h, const int & i ) |
---|
| 556 | { |
---|
| 557 | bool bruch = isOn(SW_RATIONAL); |
---|
| 558 | if (i ==0) |
---|
| 559 | { return factorize(h); |
---|
| 560 | } |
---|
| 561 | |
---|
| 562 | |
---|
| 563 | CanonicalForm g =h; |
---|
| 564 | |
---|
| 565 | Variable x =g.mvar(); |
---|
| 566 | |
---|
| 567 | int AnzExt = i; // Über welcher Erweiterung arbeite ich gerade ... |
---|
| 568 | Variable l(-AnzExt); //... und welche algebr. Variable gehört dazu ? |
---|
| 569 | Variable y('_'); |
---|
| 570 | |
---|
| 571 | |
---|
| 572 | CanonicalForm MiPo, Norm, NormAbl, Factor_Norm,dummy1, dummy2, Nenner,LeaC; |
---|
| 573 | |
---|
| 574 | CFFList Result; |
---|
| 575 | CFFList dummy; |
---|
| 576 | |
---|
| 577 | bool is = true; |
---|
| 578 | |
---|
| 579 | int k = 0; |
---|
| 580 | g = g(y,l); //die algeb. Variable wird durch eine Polynomvariable ersetzt |
---|
| 581 | MiPo = getMipo(l,y); |
---|
| 582 | |
---|
| 583 | Norm = resultant(MiPo,g,y); //norm von g als Funk. in x und y (l->y) bzgl y |
---|
| 584 | NormAbl = Norm.deriv(); |
---|
| 585 | // Off(SW_RATIONAL); |
---|
| 586 | is = !MyGCDmod(Norm,NormAbl).inBaseDomain(); //ist die Norm quadratfrei ? |
---|
| 587 | while (is) |
---|
| 588 | { |
---|
| 589 | k++; |
---|
| 590 | CanonicalForm t = g; |
---|
| 591 | t = t(x-k*y,x); //wenn nicht, wird g gestört und die neue Norm berechnet |
---|
| 592 | |
---|
| 593 | On(SW_RATIONAL); |
---|
| 594 | Norm = resultant(MiPo,t,y); //Problem tritt hier auf, bei AnzExt = 1 |
---|
| 595 | Off(SW_RATIONAL); |
---|
| 596 | NormAbl = Norm.deriv(); |
---|
| 597 | is = ! MyGCDmod(Norm,NormAbl).inBaseDomain(); |
---|
| 598 | //cout << "ggt der Norm: "<< MyGCDmod(Norm,NormAbl)<<endl; |
---|
| 599 | } |
---|
| 600 | AnzExt--; |
---|
| 601 | if (bruch) |
---|
| 602 | On(SW_RATIONAL); |
---|
| 603 | else |
---|
| 604 | Off(SW_RATIONAL); |
---|
| 605 | if (AnzExt == 0) //sind alle Erweiterungen abgearbeitet... |
---|
| 606 | { |
---|
| 607 | Result = factorize(Norm); //... wird die Norm Faktorisiert |
---|
| 608 | } |
---|
| 609 | else |
---|
| 610 | { |
---|
| 611 | Result = FactorizeNorm(Norm, AnzExt); //wenn nicht, kommt die nächste erweiterung dran |
---|
| 612 | } |
---|
| 613 | CFFListIterator J=Result; |
---|
| 614 | for ( ; J.hasItem(); J++) |
---|
[806c18] | 615 | { |
---|
[6ead9d] | 616 | Factor_Norm = J.getItem().factor(); |
---|
| 617 | Factor_Norm = Factor_Norm(x+k*l,x); // die Störungen werden rückgänig gemacht |
---|
| 618 | dummy.append(CFFactor(Factor_Norm)); |
---|
| 619 | } |
---|
| 620 | return dummy; |
---|
| 621 | } |
---|
| 622 | |
---|
| 623 | |
---|
| 624 | /* Bereitet die Form h vor, ruft FactorizeNorm(h) auf und rekonstruiert daraus die |
---|
| 625 | Faktoren von h */ |
---|
| 626 | |
---|
| 627 | CFFList MyFactorize(const CanonicalForm & h) |
---|
| 628 | { |
---|
| 629 | bool bruch = isOn(SW_RATIONAL); |
---|
| 630 | |
---|
| 631 | CanonicalForm g = h; |
---|
| 632 | CFFList FacNorm, Result; // Faktorisierung der Norm und das Ergebniss |
---|
| 633 | CanonicalForm Abl_g, LeaCoeff_g, normiert_g, g_quadrat ; //Ableitung, führender Koeff. und Normierung von g |
---|
| 634 | CanonicalForm Norm,NormAbl; |
---|
| 635 | CanonicalForm Factor_Norm, Fac; |
---|
| 636 | CanonicalForm dummy, g_origin; |
---|
| 637 | CanonicalForm Nenner,warte; |
---|
| 638 | |
---|
| 639 | Variable x = g.mvar(); |
---|
| 640 | |
---|
| 641 | int exp =0; |
---|
| 642 | int DegAlt, DegNeu; |
---|
| 643 | On(SW_RATIONAL); |
---|
| 644 | |
---|
| 645 | /* Initzialisierung, faktorisiert wird CF g */ |
---|
| 646 | g_origin = g; |
---|
| 647 | LeaCoeff_g = g.LC(); |
---|
| 648 | //g /=LeaCoeff_g; |
---|
| 649 | Nenner = MyDen(g); |
---|
| 650 | g *= power(Nenner, g.degree()); |
---|
| 651 | g *= power(LeaCoeff_g,g.degree()-1); |
---|
| 652 | g = g(x/(Nenner*LeaCoeff_g),x); |
---|
| 653 | Abl_g = g.deriv(); |
---|
| 654 | DegAlt = g.degree(); |
---|
| 655 | g_quadrat=g; |
---|
| 656 | g /= MyGCDmod(g,Abl_g); // g wird quadratfrei gemacht |
---|
| 657 | DegNeu = g.degree(); |
---|
| 658 | |
---|
| 659 | //g = g/LeaCoeff_g; // und normiert |
---|
| 660 | //CurrentExtension(); |
---|
| 661 | FacNorm = FactorizeNorm(g,ExtensionLevel()); |
---|
| 662 | CFFListIterator J=FacNorm; |
---|
| 663 | J.lastItem(); |
---|
| 664 | // g = g*MyDen(g); |
---|
| 665 | // |
---|
| 666 | g = h ; |
---|
| 667 | |
---|
| 668 | for ( ; J.hasItem(); J--) // Iteration über die Faktoren der Norm |
---|
| 669 | { |
---|
| 670 | Factor_Norm = J.getItem().factor(); |
---|
| 671 | |
---|
| 672 | Fac = MyGCDmod(g,Factor_Norm); //Ergebniss wird hochgeliftet |
---|
| 673 | |
---|
| 674 | Fac = Fac/Fac.LC(); // und normiert |
---|
| 675 | |
---|
| 676 | /* Ermittlung der Exponenten der einzelnen Faktoren */ |
---|
| 677 | |
---|
| 678 | exp = 1; // für den FaKtor mit Grad 0 |
---|
| 679 | dummy = g_quadrat; |
---|
| 680 | |
---|
| 681 | if (!Fac.inBaseDomain()) // echter Faktor ? |
---|
| 682 | { |
---|
| 683 | exp = 0; |
---|
| 684 | while ( 0==dummy%Fac && !dummy.inBaseDomain()) // Wie oft Teilt der Faktor das Polynom ? |
---|
| 685 | { |
---|
| 686 | dummy =dummy/Fac; |
---|
| 687 | exp++; |
---|
| 688 | } |
---|
| 689 | Fac = Fac(x*(Nenner*LeaCoeff_g),x); |
---|
| 690 | |
---|
| 691 | Fac /= Fac.LC(); |
---|
| 692 | } |
---|
| 693 | |
---|
| 694 | else |
---|
| 695 | { |
---|
| 696 | Fac *= LeaCoeff_g; |
---|
| 697 | g *= LeaCoeff_g; |
---|
| 698 | } |
---|
| 699 | |
---|
| 700 | g /=power(Fac,exp); |
---|
| 701 | |
---|
| 702 | Result.append(CFFactor( Fac, exp )); // Faktor wird an Result gehängt |
---|
| 703 | } |
---|
| 704 | if (bruch) |
---|
| 705 | On(SW_RATIONAL); |
---|
| 706 | else |
---|
| 707 | Off(SW_RATIONAL); |
---|
| 708 | return Result; // und zurückgegeben |
---|
| 709 | } |
---|
| 710 | |
---|
| 711 | CFFList AbsFactorize(const CanonicalForm & a) |
---|
| 712 | { |
---|
| 713 | CanonicalForm f = a; |
---|
| 714 | CanonicalForm Factor, NewFactor,dummy3, Nenner,LeadC; |
---|
| 715 | CanonicalForm Coeff=f.LC(); |
---|
| 716 | |
---|
| 717 | Variable x =a.mvar(); |
---|
| 718 | CFFList dummy, dummy2; |
---|
| 719 | CFFList result, empty; |
---|
| 720 | empty.append(CFFactor(1,1)); |
---|
| 721 | bool NewRoot = false; |
---|
| 722 | bool fertig; |
---|
| 723 | |
---|
| 724 | LeadC = f.LC(); |
---|
| 725 | f *= power(LeadC, f.degree()-1); |
---|
| 726 | Nenner = MyDen(f); |
---|
| 727 | f *= power(Nenner, f.degree()); |
---|
| 728 | f = f(x/(Nenner*LeadC), x); |
---|
| 729 | result = MyFactorize(f); |
---|
| 730 | |
---|
| 731 | CFFListIterator L = result; |
---|
| 732 | fertig = true; |
---|
| 733 | for(; L.hasItem();L++) |
---|
| 734 | { |
---|
| 735 | if (L.getItem().factor().degree() >1) |
---|
| 736 | { |
---|
| 737 | fertig = false; |
---|
| 738 | } |
---|
| 739 | } |
---|
| 740 | |
---|
| 741 | while(!fertig) |
---|
| 742 | { |
---|
| 743 | dummy = result; |
---|
| 744 | CFFListIterator J = dummy; |
---|
| 745 | result = empty; |
---|
| 746 | for ( ; J.hasItem(); J++) // Iteration über die Faktoren der Norm |
---|
| 747 | { |
---|
| 748 | Factor = J.getItem().factor(); |
---|
| 749 | |
---|
| 750 | if (Factor.degree() != 0 && Factor.degree() != 1 && !NewRoot) |
---|
| 751 | { |
---|
| 752 | Reduce(false); |
---|
| 753 | Variable u = rootOf(Factor); |
---|
| 754 | Reduce(true); |
---|
| 755 | NewRoot = true; |
---|
| 756 | result.append(CFFactor((x-u),1)); |
---|
| 757 | Factor /= (x-u); |
---|
| 758 | } |
---|
| 759 | |
---|
| 760 | |
---|
| 761 | if (Factor.degree() != 0 && Factor.degree() != 1 && NewRoot) |
---|
| 762 | { |
---|
| 763 | dummy2 = MyFactorize(Factor); |
---|
| 764 | |
---|
| 765 | CFFListIterator H = dummy2; |
---|
| 766 | for ( ; H.hasItem(); H++) // Iteration über die Faktoren der Norm |
---|
| 767 | { |
---|
| 768 | NewFactor = H.getItem().factor(); |
---|
| 769 | if (!NewFactor.inBaseDomain()) |
---|
| 770 | { |
---|
| 771 | result.append(CFFactor(NewFactor, H.getItem().exp()*J.getItem().exp())); |
---|
| 772 | } |
---|
| 773 | else |
---|
| 774 | { |
---|
| 775 | Coeff *=H.getItem().factor(); |
---|
| 776 | } |
---|
| 777 | } |
---|
| 778 | } |
---|
| 779 | if ( Factor.degree() == 0) |
---|
| 780 | { |
---|
| 781 | Coeff *=Factor; |
---|
| 782 | } |
---|
| 783 | if( Factor.degree() == 1) |
---|
| 784 | { |
---|
| 785 | result.append(CFFactor(Factor,J.getItem().exp())); |
---|
| 786 | } |
---|
| 787 | } |
---|
| 788 | NewRoot = false; |
---|
| 789 | CFFListIterator K = result; |
---|
| 790 | fertig = true; |
---|
| 791 | |
---|
| 792 | for(; K.hasItem();K++) |
---|
| 793 | { |
---|
| 794 | |
---|
| 795 | if (K.getItem().factor().degree() >1) |
---|
| 796 | { |
---|
| 797 | fertig = false; |
---|
| 798 | } |
---|
| 799 | } |
---|
| 800 | } |
---|
| 801 | CFFList result2; |
---|
| 802 | //result2.append(CFFactor(Coeff)); |
---|
| 803 | CFFListIterator K = result; |
---|
| 804 | for(; K.hasItem();K++) |
---|
| 805 | { |
---|
| 806 | dummy3 = K.getItem().factor(); |
---|
| 807 | if (dummy3.degree() == 0) |
---|
| 808 | { dummy3 *= Coeff; |
---|
| 809 | } |
---|
| 810 | else |
---|
| 811 | { |
---|
| 812 | dummy3 = dummy3(x*Nenner*LeadC,x); |
---|
| 813 | dummy3 /= dummy3.LC(); |
---|
| 814 | } |
---|
| 815 | |
---|
| 816 | result2.append(CFFactor(dummy3,K.getItem().exp())); |
---|
| 817 | } |
---|
| 818 | return result2; |
---|
| 819 | } |
---|
| 820 | |
---|
| 821 | |
---|
| 822 | // |
---|
| 823 | // |
---|
| 824 | CanonicalForm Bigcd( const CanonicalForm& f, const CanonicalForm& g) |
---|
| 825 | { |
---|
| 826 | |
---|
| 827 | |
---|
| 828 | if( f.level() < 0 ) return 1; |
---|
| 829 | if( g.level() < 0 ) return 1; |
---|
| 830 | |
---|
| 831 | CFArray A; |
---|
| 832 | |
---|
| 833 | int i=0; |
---|
| 834 | int r; |
---|
| 835 | |
---|
| 836 | Variable x = f.mvar(); |
---|
| 837 | Variable y = g.mvar(); |
---|
| 838 | |
---|
| 839 | // Wahl als Hauptvariable ? |
---|
| 840 | // |
---|
| 841 | if( x.level() >= y.level() ) x = y; |
---|
| 842 | |
---|
| 843 | CanonicalForm Cf, Cg, gamma, c,T; |
---|
| 844 | CanonicalForm F=f; |
---|
| 845 | CanonicalForm G=g; |
---|
| 846 | |
---|
| 847 | Cf = MyContent(f,x); //changed |
---|
| 848 | Cg = MyContent(g,x); //changed |
---|
| 849 | F = F/Cf; |
---|
| 850 | G = G/Cg; |
---|
| 851 | gamma = MyGCDmod( LC(F,x), LC(G,x) ); |
---|
| 852 | |
---|
| 853 | A = subResChain(F,G,x); |
---|
| 854 | |
---|
| 855 | c = MyGCDmod( Cf, Cg ); |
---|
| 856 | |
---|
| 857 | r = A.size(); |
---|
| 858 | |
---|
| 859 | while( A[i].isZero() ) i++; |
---|
| 860 | |
---|
| 861 | F = A[i]; |
---|
| 862 | |
---|
[806c18] | 863 | if( degree(F,x) == 0 ) |
---|
| 864 | if( c.level() < 0 ) return 1; else return c; |
---|
[6ead9d] | 865 | |
---|
| 866 | F = gamma*F/LC(F, x); |
---|
| 867 | F = F/content(F,x); |
---|
| 868 | |
---|
| 869 | F = c*F; |
---|
| 870 | |
---|
| 871 | c = F.LC(); |
---|
| 872 | |
---|
| 873 | while( c.level()>0 ) c = c.LC(); |
---|
| 874 | |
---|
| 875 | F=F/c; |
---|
| 876 | |
---|
| 877 | if( F.level() < 0 ) return 1; |
---|
| 878 | |
---|
| 879 | return F; |
---|
| 880 | } |
---|
| 881 | |
---|
| 882 | |
---|
| 883 | |
---|
| 884 | |
---|
| 885 | CanonicalForm MYGCD( const CanonicalForm& f, const CanonicalForm& g) |
---|
| 886 | { |
---|
| 887 | |
---|
| 888 | // FIX ME: CONSTANT FIELD |
---|
| 889 | // |
---|
| 890 | // |
---|
| 891 | // |
---|
| 892 | // |
---|
| 893 | // |
---|
| 894 | if( f.level() < 0 && g.level() < 0) return 1; |
---|
| 895 | if( (f.level() < 0 && g.level() > 0) || |
---|
| 896 | (f.level() > 0 && g.level() <0 ) ) return 1; |
---|
| 897 | |
---|
| 898 | int i; |
---|
| 899 | |
---|
| 900 | CFList L; |
---|
| 901 | for (i=1; i<= level(f); i++) |
---|
| 902 | if( f != f(0,i) ) |
---|
| 903 | L.append(i); |
---|
| 904 | |
---|
| 905 | int nvf = L.length(); |
---|
| 906 | |
---|
| 907 | for (i=1; i<= level(g); i++) |
---|
| 908 | if( g != g(0,i) ) |
---|
| 909 | L.append(i); |
---|
[806c18] | 910 | |
---|
[6ead9d] | 911 | int nvg = L.length(); |
---|
| 912 | |
---|
[806c18] | 913 | |
---|
[6ead9d] | 914 | |
---|
| 915 | if( f.level() < 0 && g.level() < 0 ) { ; |
---|
| 916 | return 1; } |
---|
[806c18] | 917 | |
---|
[6ead9d] | 918 | CFArray A; |
---|
[806c18] | 919 | |
---|
[6ead9d] | 920 | i=0; |
---|
| 921 | int r; |
---|
| 922 | |
---|
| 923 | Variable x = f.mvar(); |
---|
| 924 | Variable y = g.mvar(); |
---|
| 925 | |
---|
| 926 | // Wahl als Hauptvariable ? |
---|
| 927 | // |
---|
| 928 | if( x.level() >= y.level() ) x = y; |
---|
| 929 | |
---|
| 930 | CanonicalForm Cf, Cg, gamma, c,T; |
---|
| 931 | CanonicalForm F=f; |
---|
| 932 | CanonicalForm G=g; |
---|
| 933 | |
---|
| 934 | |
---|
| 935 | Cf = MyContent(f,x); |
---|
| 936 | Cg = MyContent(g,x); |
---|
| 937 | F = F/Cf; |
---|
| 938 | G = G/Cg; |
---|
| 939 | |
---|
| 940 | if( nvf <= 1 && nvg <=1 ) |
---|
| 941 | { |
---|
[806c18] | 942 | gamma = MyGCDmod( LC(F,x), LC(G,x) ); |
---|
[6ead9d] | 943 | c = MyGCDmod( Cf, Cg ); |
---|
| 944 | } |
---|
| 945 | else |
---|
| 946 | { |
---|
| 947 | gamma = MYGCD( LC(F,x), LC(G,x) ); |
---|
| 948 | c = MYGCD( Cf, Cg ); |
---|
| 949 | } |
---|
| 950 | A = subResChain(F,G,x); |
---|
| 951 | |
---|
| 952 | r = A.size(); |
---|
| 953 | |
---|
| 954 | while( A[i].isZero() ) i++; |
---|
| 955 | |
---|
| 956 | F = A[i]; |
---|
| 957 | |
---|
[806c18] | 958 | if( degree(F,x) == 0 ) |
---|
| 959 | if( c.level() < 0 ) return 1; else return c; |
---|
[6ead9d] | 960 | |
---|
| 961 | F = gamma*F/LC(F, x); |
---|
| 962 | F = F/MyContent(F,x); |
---|
| 963 | |
---|
| 964 | F = c*F; |
---|
| 965 | |
---|
| 966 | c = F.LC(); |
---|
| 967 | |
---|
| 968 | while( c.level()>0 ) c = c.LC(); |
---|
| 969 | |
---|
| 970 | F=F/c; |
---|
| 971 | |
---|
| 972 | //if( F.level() < 0 ) return 1; |
---|
| 973 | |
---|
| 974 | return F; |
---|
| 975 | } |
---|
| 976 | |
---|
| 977 | |
---|
| 978 | |
---|
| 979 | CFFList Mysqrfree_local( const CanonicalForm& F, const Variable& v) |
---|
| 980 | { |
---|
| 981 | int i=1; |
---|
| 982 | CanonicalForm f=F; |
---|
| 983 | CanonicalForm g, qi, fp, wp, temp1, temp2, temp3, temp4, pA; |
---|
| 984 | CFFList L; |
---|
| 985 | |
---|
| 986 | |
---|
| 987 | g = MyContent(f,v); |
---|
| 988 | |
---|
| 989 | if( g != 1 ) |
---|
| 990 | L.append( CFFactor(g,1) ); |
---|
| 991 | |
---|
| 992 | pA = f/g; |
---|
| 993 | |
---|
| 994 | fp = deriv( pA, v); |
---|
| 995 | |
---|
| 996 | temp1 = MYGCD( pA, fp ); |
---|
| 997 | |
---|
| 998 | if( temp1 == 1 ){ L.append( CFFactor(pA,1) ); return L; } |
---|
| 999 | else |
---|
| 1000 | { |
---|
| 1001 | temp2 = pA/temp1; |
---|
| 1002 | temp3 = fp/temp1; |
---|
| 1003 | wp = deriv(temp2,v); |
---|
| 1004 | temp4 = temp3 - wp; |
---|
| 1005 | |
---|
| 1006 | while( !temp4.isZero() ) |
---|
| 1007 | { |
---|
| 1008 | CanonicalForm qi = MYGCD( temp2, temp4); |
---|
[806c18] | 1009 | if( qi != 1 ) L.append( CFFactor( qi, i ) ); |
---|
[6ead9d] | 1010 | i++; |
---|
| 1011 | temp2 = temp2/qi; |
---|
| 1012 | temp3 = temp4/qi; |
---|
| 1013 | temp4 = temp3-deriv(temp2, v); |
---|
| 1014 | } |
---|
| 1015 | |
---|
| 1016 | if( temp2 != 1 ) L.append( CFFactor( temp2, i ) ); |
---|
| 1017 | |
---|
| 1018 | } |
---|
| 1019 | |
---|
| 1020 | return L; |
---|
| 1021 | } |
---|
| 1022 | |
---|
| 1023 | CFFList Mysqrfree( const CanonicalForm& F ) |
---|
| 1024 | { |
---|
| 1025 | CFFList L, M, V; |
---|
| 1026 | CFFList N; |
---|
| 1027 | CanonicalForm vars=getVars(F); |
---|
| 1028 | Variable v; |
---|
| 1029 | CanonicalForm s; |
---|
| 1030 | bool b; |
---|
| 1031 | |
---|
| 1032 | L.append( CFFactor(F,1) ); |
---|
| 1033 | |
---|
| 1034 | int n = F.level(); |
---|
| 1035 | int *vrs = new int[n+1]; |
---|
| 1036 | for ( int i = 0; i <= n; i++ ) vars[i] = 0; |
---|
| 1037 | for ( CFIterator I = F; I.hasTerms(); ++I ) fillVarsRec( I.coeff(), vrs ); |
---|
| 1038 | |
---|
[806c18] | 1039 | N.append( CFFactor(F,1) ); |
---|
[6ead9d] | 1040 | |
---|
| 1041 | int i = n+1; |
---|
| 1042 | |
---|
| 1043 | while( i >= 0 ) |
---|
| 1044 | { |
---|
[806c18] | 1045 | b = 0; |
---|
[6ead9d] | 1046 | |
---|
| 1047 | if( i == 0 ){ v = mvar(F); b=1 ;} |
---|
| 1048 | else |
---|
[806c18] | 1049 | if( vrs[i] != 0 ){ b=1; v= Variable(i);} |
---|
| 1050 | if( vrs[i] == 0 ) i--; |
---|
[6ead9d] | 1051 | |
---|
| 1052 | if( b ) |
---|
[806c18] | 1053 | { |
---|
[6ead9d] | 1054 | for( CFFListIterator J = L; J.hasItem(); J++ ) |
---|
| 1055 | { |
---|
| 1056 | M = Mysqrfree_local( J.getItem().factor() , v ); |
---|
| 1057 | |
---|
| 1058 | for( CFFListIterator K = M; K.hasItem(); K++ ) |
---|
| 1059 | { |
---|
| 1060 | if( K.getItem().factor().level() > 0 ) |
---|
| 1061 | { |
---|
| 1062 | N.append( CFFactor( K.getItem().factor(), K.getItem().exp()+J.getItem().exp()-1 )); } |
---|
| 1063 | } |
---|
| 1064 | N.removeFirst(); |
---|
| 1065 | } |
---|
| 1066 | if( N.length() == L.length() ) i -= 1; |
---|
[806c18] | 1067 | L=N; |
---|
[6ead9d] | 1068 | } |
---|
| 1069 | } |
---|
| 1070 | |
---|
| 1071 | return L; |
---|
| 1072 | |
---|
| 1073 | } |
---|
| 1074 | #endif |
---|