1 | |
---|
2 | #include "canonicalform.h" |
---|
3 | #ifdef HAVE_BIFAC |
---|
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; |
---|
38 | if( F.inBaseDomain() ) return F; |
---|
39 | |
---|
40 | if( level(F) < 0 ) return 1; |
---|
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); |
---|
49 | if( t == 1 || t == -1 ) return 1; |
---|
50 | r = MYGCD( r, t); |
---|
51 | if( r == 1 ) return 1; |
---|
52 | g = g - power(x,degree(g,x))*t; |
---|
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++) |
---|
615 | { |
---|
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 | |
---|
863 | if( degree(F,x) == 0 ) |
---|
864 | if( c.level() < 0 ) return 1; else return c; |
---|
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); |
---|
910 | |
---|
911 | int nvg = L.length(); |
---|
912 | |
---|
913 | |
---|
914 | |
---|
915 | if( f.level() < 0 && g.level() < 0 ) { ; |
---|
916 | return 1; } |
---|
917 | |
---|
918 | CFArray A; |
---|
919 | |
---|
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 | { |
---|
942 | gamma = MyGCDmod( LC(F,x), LC(G,x) ); |
---|
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 | |
---|
958 | if( degree(F,x) == 0 ) |
---|
959 | if( c.level() < 0 ) return 1; else return c; |
---|
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); |
---|
1009 | if( qi != 1 ) L.append( CFFactor( qi, i ) ); |
---|
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 | |
---|
1039 | N.append( CFFactor(F,1) ); |
---|
1040 | |
---|
1041 | int i = n+1; |
---|
1042 | |
---|
1043 | while( i >= 0 ) |
---|
1044 | { |
---|
1045 | b = 0; |
---|
1046 | |
---|
1047 | if( i == 0 ){ v = mvar(F); b=1 ;} |
---|
1048 | else |
---|
1049 | if( vrs[i] != 0 ){ b=1; v= Variable(i);} |
---|
1050 | if( vrs[i] == 0 ) i--; |
---|
1051 | |
---|
1052 | if( b ) |
---|
1053 | { |
---|
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; |
---|
1067 | L=N; |
---|
1068 | } |
---|
1069 | } |
---|
1070 | |
---|
1071 | return L; |
---|
1072 | |
---|
1073 | } |
---|
1074 | #endif |
---|