source: git/factory/abs_fac.cc @ 6ead9d

spielwiese
Last change on this file since 6ead9d was 6ead9d, checked in by Hans Schönemann <hannes@…>, 18 years ago
*Hannnes: bivar stuff git-svn-id: file:///usr/local/Singular/svn/trunk@8173 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 19.6 KB
Line 
1
2#ifdef HAVE_BIFAC
3#include "factory.h"
4#ifndef NOSTREAMIO
5#include<fstream>
6#endif
7#include<math.h>
8#include <sys/timeb.h>
9
10
11
12static void
13fillVarsRec ( 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
25int ExtensionLevel();
26void Reduce( bool);
27
28CanonicalForm MYGCD( const CanonicalForm& f, const CanonicalForm& g);
29
30
31CanonicalForm 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
57void 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
75CanonicalForm 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
121CanonicalForm 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*/
153CanonicalForm 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
212CanonicalForm 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
259CanonicalForm 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 }
319if (bruch)
320       On(SW_RATIONAL);
321      else
322       Off(SW_RATIONAL);
323return Result;
324}
325
326/*Rational Rekonstruction für a mod b*/
327
328CanonicalForm 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  }
401if (bruch)
402       On(SW_RATIONAL);
403      else
404       Off(SW_RATIONAL);
405return Result;
406}
407
408/*Berechnet lokale ggT's der Formen a und b und liftet sie wieder hoch*/
409
410CanonicalForm 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;
415CanonicalForm LeadA, LeadB;
416CanonicalForm Kandidat,Kandidat1;
417CanonicalForm f,g, Result;
418CanonicalForm NennerA =1;
419CanonicalForm NennerB=1;
420CanonicalForm ZahlerA=1;
421CanonicalForm ZahlerB=1;
422int treffer = 0;
423int Level;
424CanonicalForm Modulo;
425int i = 0;
426bool TryAgain = 1;
427int Primes[1228];
428
429
430
431for (int i = 0;i <1228;i++)
432{
433 Primes[i]=cf_getPrime(i+1);
434
435}
436Level=a.level();
437
438if (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
555CFFList 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
625Faktoren von h */
626
627CFFList 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
711CFFList AbsFactorize(const CanonicalForm  & a)
712{
713CanonicalForm f = a;
714CanonicalForm Factor, NewFactor,dummy3, Nenner,LeadC;
715CanonicalForm Coeff=f.LC();
716
717Variable x =a.mvar();
718CFFList dummy, dummy2;
719CFFList result, empty;
720empty.append(CFFactor(1,1));
721bool NewRoot = false;
722bool fertig;
723
724LeadC = f.LC();
725f *= power(LeadC, f.degree()-1);
726Nenner = MyDen(f);
727f *= power(Nenner, f.degree());
728f = f(x/(Nenner*LeadC), x);
729result = MyFactorize(f);
730
731CFFListIterator 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
741while(!fertig)
742{
743 dummy = result;
744 CFFListIterator J = dummy;
745 result = empty;
746for ( ; 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 }
801CFFList result2;
802//result2.append(CFFactor(Coeff));
803CFFListIterator K = result;
804for(; 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     }
818return result2;
819}
820
821
822//
823//
824CanonicalForm 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
885CanonicalForm 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
979CFFList 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
1023CFFList 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
Note: See TracBrowser for help on using the repository browser.