source: git/factory/abs_fac.cc @ 2db7ae

spielwiese
Last change on this file since 2db7ae was 2db7ae, checked in by Hans Schoenemann <hannes@…>, 13 years ago
disabled cf_hnf: does not compile with NTL with NTL_STD_CXX set
  • Property mode set to 100644
File size: 19.6 KB
Line 
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
11static void
12fillVarsRec ( 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
24int ExtensionLevel();
25void Reduce( bool);
26
27CanonicalForm MYGCD( const CanonicalForm& f, const CanonicalForm& g);
28
29
30CanonicalForm 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
56void 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
74CanonicalForm 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
120CanonicalForm 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*/
152CanonicalForm 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
211CanonicalForm 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
258CanonicalForm 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 }
318if (bruch)
319       On(SW_RATIONAL);
320      else
321       Off(SW_RATIONAL);
322return Result;
323}
324
325/*Rational Rekonstruction für a mod b*/
326
327CanonicalForm 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  }
400if (bruch)
401       On(SW_RATIONAL);
402      else
403       Off(SW_RATIONAL);
404return Result;
405}
406
407/*Berechnet lokale ggT's der Formen a und b und liftet sie wieder hoch*/
408
409CanonicalForm 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;
414CanonicalForm LeadA, LeadB;
415CanonicalForm Kandidat,Kandidat1;
416CanonicalForm f,g, Result;
417CanonicalForm NennerA =1;
418CanonicalForm NennerB=1;
419CanonicalForm ZahlerA=1;
420CanonicalForm ZahlerB=1;
421int treffer = 0;
422int Level;
423CanonicalForm Modulo;
424int i = 0;
425bool TryAgain = 1;
426int Primes[1228];
427
428
429
430for (int i = 0;i <1228;i++)
431{
432 Primes[i]=cf_getPrime(i+1);
433
434}
435Level=a.level();
436
437if (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
554CFFList 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
624Faktoren von h */
625
626CFFList 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
710CFFList AbsFactorize(const CanonicalForm  & a)
711{
712CanonicalForm f = a;
713CanonicalForm Factor, NewFactor,dummy3, Nenner,LeadC;
714CanonicalForm Coeff=f.LC();
715
716Variable x =a.mvar();
717CFFList dummy, dummy2;
718CFFList result, empty;
719empty.append(CFFactor(1,1));
720bool NewRoot = false;
721bool fertig;
722
723LeadC = f.LC();
724f *= power(LeadC, f.degree()-1);
725Nenner = MyDen(f);
726f *= power(Nenner, f.degree());
727f = f(x/(Nenner*LeadC), x);
728result = MyFactorize(f);
729
730CFFListIterator 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
740while(!fertig)
741{
742 dummy = result;
743 CFFListIterator J = dummy;
744 result = empty;
745for ( ; 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 }
800CFFList result2;
801//result2.append(CFFactor(Coeff));
802CFFListIterator K = result;
803for(; 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     }
817return result2;
818}
819
820
821//
822//
823CanonicalForm 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
884CanonicalForm 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
978CFFList 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
1022CFFList 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
Note: See TracBrowser for help on using the repository browser.