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