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