source: git/libpolys/coeffs/AE.cc @ 5b5409b

spielwiese
Last change on this file since 5b5409b was 4154bb, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Fix "if (deg=0)" in new exp. coeffs (see issue #374)
  • Property mode set to 100755
File size: 25.7 KB
Line 
1#include <gmp.h>
2#include <math.h>
3#include "AE.h"
4
5
6
7using namespace std;
8
9//Konstruktoren
10
11int_poly::int_poly()
12{
13    //coef = reinterpret_cast<mpz_t*> (omAlloc(100*sizeof(mpz_t)));
14    deg=-1;
15    mpz_init_set_ui(coef[0],0);
16}
17
18
19
20int_poly::int_poly(int n, mpz_t *a)
21{
22    //coef = reinterpret_cast<mpz_t*> (omAlloc(100*sizeof(mpz_t)));
23    deg=n;
24
25    for ( int i=0;i<=n;i++)
26    {
27        mpz_init_set(coef[i], a[i]);
28    }
29
30}
31
32/*
33//Destruktor
34
35int_poly::~int_poly()
36{
37        delete[] coef;
38}
39
40*/
41
42
43
44
45// Arithmetik
46
47
48//Additionen
49
50//Standard - Addition
51int_poly int_poly::poly_add(const int_poly a, const int_poly b)
52{
53    if (a.deg >=b.deg)
54    {
55
56        deg=a.deg;
57
58        for ( int i=0;i<=b.deg;i++)
59        {
60            mpz_add(coef[i],a.coef[i],b.coef[i]);
61        }
62
63        for ( int i=b.deg+1;i<=a.deg;i++)
64        {
65            mpz_init_set(coef[i],a.coef[i]);
66        }
67        //Hier nötige Grad Korrektur
68        int i;
69        i=deg;
70        while(mpz_sgn(coef[i])==0 && i>=0)
71        {deg--;i--;}
72    }
73
74    else {poly_add(b,a);  }
75
76}
77
78//Überschreibende Addition
79
80int_poly int_poly::poly_add_to(const int_poly g)
81{
82    this->poly_add(*this,g);
83}
84
85//Addition einer Konstanten
86int_poly int_poly::poly_add_const(int_poly f,const mpz_t a)
87{
88    if (f.is_zero()==1)
89        poly_set(a);
90    else
91    {
92        poly_set(f);
93        mpz_add(coef[0],coef[0],a);
94        // Grad Korrektur
95        if (deg==0 && mpz_sgn(coef[0])==0)
96            poly_set_zero();
97    }
98
99}
100
101
102//To Variante Addition einer Konstanten
103
104int_poly int_poly::poly_add_const_to(const mpz_t a)
105{
106    this->poly_add_const(*this,a);
107}
108
109//Monom Addition
110int_poly int_poly::poly_add_mon(int_poly f, mpz_t a, int i)
111{
112    poly_set(f);
113
114    if (i<=deg  && is_zero()==0)
115    {
116        mpz_add(coef[i],coef[i],a);
117        // Grad Korrektur
118        if (deg==i && mpz_sgn(coef[i])==0)
119            deg--;
120    }
121    else if (is_zero()==1)
122    {
123        deg=i;
124        for(int j=0;j<=i;j++)
125        {
126            mpz_init_set_ui(coef[j],0);
127        }
128        mpz_add(coef[i],coef[i],a);
129
130    }
131    else if (i>deg)
132    {
133        for(int j=i;j>deg;j--)
134        {
135            mpz_init_set_ui(coef[j],0);
136        }
137        mpz_add(coef[i],coef[i],a);
138        deg=i;
139    }
140}
141
142//To Variante Monomaddition
143int_poly int_poly::poly_add_mon_to(mpz_t a, int i)
144{
145    if (i<=deg  && is_zero()==0)
146    {
147        mpz_add(coef[i],coef[i],a);
148    }
149    else if (is_zero()==1)
150    {
151        deg=i;
152        for(int j=0;j<=i;j++)
153        {
154            mpz_init_set_ui(coef[j],0);
155        }
156        mpz_add(coef[i],coef[i],a);
157
158    }
159    else if (i>deg)
160    {
161        for(int j=i;j>deg;j--)
162        {
163            mpz_init_set_ui(coef[j],0);
164        }
165        mpz_add(coef[i],coef[i],a);
166        deg=i;
167    }
168    //Hier nötige Grad Korrektur
169    int k=deg;
170    while(mpz_sgn(coef[k])==0 && k>=0)
171    {deg--;k--;}
172
173}
174
175//Subtraktionen
176
177int_poly int_poly::poly_sub(const int_poly a, const int_poly b)
178{
179    int_poly temp;
180    temp.poly_set(b);
181    temp.poly_neg();
182    poly_add(a,temp);
183
184    // Grad Korrektur
185    int i;
186    i=deg;
187    while(mpz_sgn(coef[i])==0 && i>=0)
188    {deg--;i--;}
189
190}
191
192
193//Überschreibende Subtraktion
194
195int_poly int_poly::poly_sub_to(const int_poly b)
196{
197    this->poly_sub(*this,b);
198}
199
200//Subtraktion einer Konstanten
201int_poly int_poly::poly_sub_const(int_poly f,const mpz_t a)
202{
203    if (f.is_zero()==1)
204    {
205        poly_set(a);
206        poly_neg();
207    }
208    else
209    {
210        poly_set(f);
211        mpz_sub(coef[0],coef[0],a);
212
213    }
214    //Hier nötige Grad Korrektur
215    int i=deg;
216    while(mpz_sgn(coef[i])==0 && i>=0)
217    {deg--;i--;}
218
219}
220
221
222//To Variante Subtraktion einer Konstanten
223
224int_poly int_poly::poly_sub_const_to(const mpz_t a)
225{
226    this->poly_sub_const(*this,a);
227}
228
229
230//Monom Subtraktion
231int_poly int_poly::poly_sub_mon(const int_poly f , mpz_t a, int i)
232{
233    poly_set(f);
234
235    if (i<=deg && is_zero()!=1)
236    {
237        mpz_sub(coef[i],coef[i],a);
238        // Grad Korrektur
239        if (deg==i && mpz_sgn(coef[i])==0)
240            deg--;
241    }
242    else if (is_zero()==1)
243    {
244        for(int j=0;j<=i;j++)
245        {
246            mpz_init_set_ui(coef[j],0);
247        }
248        mpz_sub(coef[i],coef[i],a);
249        deg=i;
250    }
251    else
252    {
253        for(int j=i;j>deg;j--)
254        {
255            mpz_init_set_ui(coef[j],0);
256        }
257        mpz_sub(coef[i],coef[i],a);
258        deg=i;
259    }
260}
261
262//To Variante Monomaddition
263int_poly int_poly::poly_sub_mon_to(mpz_t a, int i)
264{
265
266    if (i<=deg && is_zero()!=1)
267    {
268        mpz_sub(coef[i],coef[i],a);
269        // Grad Korrektur
270        if (deg==i && mpz_sgn(coef[i])==0)
271            deg--;
272    }
273    else if (is_zero()==1)
274    {
275        for(int j=0;j<=i;j++)
276        {
277            mpz_init_set_ui(coef[j],0);
278        }
279        mpz_sub(coef[i],coef[i],a);
280        deg=i;
281    }
282    else
283    {
284        for(int j=i;j>deg;j--)
285        {
286            mpz_init_set_ui(coef[j],0);
287        }
288        mpz_sub(coef[i],coef[i],a);
289        deg=i;
290    }
291}
292
293
294//Multiplikationen
295
296//Multiplikation mit Monom
297int_poly int_poly::poly_mon_mult(const int_poly f, int n)
298{
299    if (f,is_zero()==1)
300    {
301        poly_set_zero();
302    }
303    else
304    {
305        deg=f.deg+n;
306        for (int i=deg;i>=n;i--)
307        {
308            mpz_init_set(coef[i],f.coef[i-n]);
309        }
310        for (int i=n-1;i>=0;i--)
311        {
312            mpz_init_set_ui(coef[i],0);
313        }
314    }
315}
316
317int_poly int_poly::poly_mon_mult_to(const int n)
318{
319    this->poly_mon_mult(*this,n);
320}
321
322
323//Multiplikation mit Skalar
324
325int_poly int_poly::poly_scalar_mult(const int_poly g, const mpz_t n)
326{
327    if (mpz_sgn(n)==0)
328        poly_set_zero();
329    else
330    {
331        deg=g.deg;
332        mpz_t temp;
333        mpz_init_set_ui(temp,0);
334        for(int i=0;i<=deg;i++)
335        {
336            mpz_mul(temp,n,g.coef[i]);
337            mpz_init_set(coef[i],temp);
338        }
339    }
340}
341
342
343
344int_poly int_poly::poly_scalar_mult(const mpz_t n, const int_poly g)
345{
346    if (mpz_sgn(n)==0)
347        poly_set_zero();
348    else
349    {
350        deg=g.deg;
351        mpz_t temp;
352        mpz_init_set_ui(temp,0);
353        for(int i=0;i<=deg;i++)
354        {
355            mpz_mul(temp,n,g.coef[i]);
356            mpz_init_set(coef[i],temp);
357        }
358    }
359}
360
361
362int_poly int_poly::poly_scalar_mult_to(const mpz_t n)
363{
364    this->poly_scalar_mult(*this,n);
365}
366
367
368
369// Negation
370
371int_poly int_poly::poly_neg()
372{
373    for (int i=0;i<=deg;i++)
374    {
375        mpz_neg(coef[i],coef[i]);
376    }
377}
378
379// Naive Multiplikation
380int_poly int_poly::poly_mult_n(int_poly a,int_poly b)
381{
382
383    if (a.is_zero()==1 || b.is_zero()==1)
384    {
385        poly_set_zero();
386    }
387    else
388    {
389        mpz_t temp;
390        mpz_init_set_ui(temp,0);
391        deg = a.deg + b.deg;
392
393        // Kopien atemp und btemp von a bzw. b, mit Nullen aufgefÃŒllt
394        int_poly atemp, btemp;
395        atemp.poly_set(a);
396        btemp.poly_set(b);
397        for(int i=a.deg+1;i<=deg;i++)
398        {
399            mpz_init_set_ui(atemp.coef[i],0);
400        }
401        for(int i=b.deg+1;i<=deg;i++)
402        {
403            mpz_init_set_ui(btemp.coef[i],0);
404        }
405        atemp.deg = deg;
406        btemp.deg = deg;
407
408        // Multiplikationsalgorithmus
409        for (int k=0; k<=deg; k++)
410        {
411            mpz_init_set_ui(coef[k],0);        // k-ter Koeffizient zunÀchst 0
412            for (int i=0; i<=k; i++)        // dann schrittweise Summe der a[i]*b[k-i]/
413            {
414                mpz_mul(temp,atemp.coef[i],btemp.coef[k-i]);
415                mpz_add(coef[k],coef[k],temp);
416            }
417        }
418
419    }
420}
421
422//Überschreibende Multiplikation
423
424int_poly int_poly::poly_mult_n_to(const int_poly g)
425{
426    this->poly_mult_n(*this,g);
427
428}
429
430// Karatsuba-Multiplikation (Weimerskirch/Paar Alg. 1), ACHTUNG VORLÄUFIGE VERSION, macht noch Fehler beim Grad und ist unelegant !!!
431int_poly int_poly::poly_mult_ka( int_poly A,  int_poly B)
432{
433
434    if (A.is_zero()==1 || B.is_zero()==1)
435    {
436        poly_set_zero();
437    }
438    else
439    {
440        // Größeren Grad feststellen
441        int n;
442        if(A.deg>=B.deg){n=A.deg+1;}
443        else{n=B.deg+1;}
444        // n auf nÀchste 2er-Potenz setzen (VORLÄUFIG!)
445        n = static_cast<int>(ceil(log(n)/log(2)));
446        n = static_cast<int>(pow(2,n));
447
448        if (n==1)
449        {
450            mpz_t AB;
451            mpz_mul(AB,A.coef[0],B.coef[0]);
452            poly_set(AB);
453        }
454        else
455        {
456            // int_polynome A und B aufspalten
457            int_poly Au, Al, Bu, Bl;
458            Au.poly_mon_div(A,n/2);
459            Al.poly_mon_div_rem(A,n/2);
460            Bu.poly_mon_div(B,n/2);
461            Bl.poly_mon_div_rem(B,n/2);
462            int_poly Alu,Blu;
463            Alu.poly_add(Al,Au);
464            Blu.poly_add(Bl,Bu);
465
466            // Teile rekursiv multiplizieren
467            int_poly D0, D1, D01;
468            D0.poly_mult_ka(Al,Bl);
469            D1.poly_mult_ka(Au,Bu);
470            D01.poly_mult_ka(Alu,Blu);
471
472            // Ergebnis zusammensetzen
473            int_poly temp;
474            D01.poly_sub_to(D0);
475            D01.poly_sub_to(D1);
476            D01.poly_mon_mult_to(n/2);
477            D1.poly_mon_mult_to(n);
478            D1.poly_add_to(D01);
479            D1.poly_add_to(D0);
480            poly_set(D1);
481        }
482    }
483}
484
485
486
487//Skalare Divisionen
488
489int_poly int_poly::poly_scalar_div( const int_poly g, const mpz_t n)
490{
491    deg=g.deg;
492    mpz_t temp;
493    mpz_init_set_ui(temp,0);
494    for(int i=0;i<=deg;i++)
495    {
496        mpz_divexact(temp,g.coef[i],n);
497        mpz_init_set(coef[i],temp);
498    }
499}
500
501
502int_poly int_poly::poly_scalar_div_to(const mpz_t n)
503{
504    this->poly_scalar_div(*this,n);
505}
506
507// Division durch Monom -  results in Quotient without remainder
508int_poly int_poly::poly_mon_div(const int_poly f, const int n)
509{
510    if (f.deg<n)
511    {
512        poly_set_zero();
513    }
514    else
515    {
516        deg=f.deg-n;
517        for (int i=0;i<=f.deg-n;i++)
518        {
519            mpz_init_set(coef[i],f.coef[n+i]);
520        }
521    }
522}
523
524// Division durch Monom - Rest
525int_poly int_poly::poly_mon_div_rem(const int_poly f, const int n)
526{
527    if (f.deg<n)
528    {
529        poly_set(f);
530    }
531    else
532    {
533        deg=n-1;
534        for (int i=0;i<=n-1;i++)
535        {
536            mpz_init_set(coef[i],f.coef[i]);
537        }
538    }
539}
540
541
542
543
544//Exakte Division nach Cohen 3.1.1 (works only if B!=0)
545int_poly int_poly::poly_div(int_poly &Q,int_poly &R, int_poly A,  int_poly B)
546{
547    if (B.is_zero()==0)
548    {
549        //Initialisierungen
550        int_poly temp;
551        R.poly_set(A);
552        Q.poly_set_zero();
553        mpz_t a;
554        mpz_init_set_ui(a,0);
555        int i;
556
557        //Algorithmus TO DO: evtl hier mit auch den Rest ausgeben
558        while (R.deg>=B.deg)
559        {
560            mpz_divexact(a,R.coef[R.deg],B.coef[B.deg]);
561            i=R.deg-B.deg;
562
563            temp.poly_mon_mult(B,i);
564            temp.poly_scalar_mult_to(a);
565
566            R.poly_sub_to(temp);
567            Q.poly_add_mon_to(a,i);
568        }
569        poly_set(Q);
570    }
571}
572
573
574//To Varainte der exakten Division
575
576int_poly int_poly::poly_div_to(int_poly &Q,int_poly &R,const int_poly B)
577{
578    this->poly_div( Q, R,*this,B);
579}
580
581// pseudo Division nach Cohen 3.1.2 (geht eleganter)
582
583int_poly int_poly::poly_pseudodiv_rem( int_poly A,  int_poly B)
584{
585
586    if (B.is_zero()==0)
587    {
588        int_poly temp;
589        int_poly R;
590        R.poly_set(A);
591        int e=A.deg-B.deg+1;
592        while (R.deg>=B.deg)
593        {
594
595            temp.poly_mon_mult(B,R.deg-B.deg);
596            temp.poly_scalar_mult_to(R.coef[R.deg]);
597            R.poly_scalar_mult_to(B.coef[B.deg]);
598            R.poly_sub_to(temp);
599            e--;
600
601        }
602        mpz_t q;
603        mpz_init_set(q,B.coef[B.deg]);
604        mpz_pow_ui(q,q,e);
605        R.poly_scalar_mult_to(q);
606        poly_set(R);
607    }
608}
609
610//To Variante Algo 3.1.2 nach Cohen
611
612int_poly int_poly::poly_pseudodiv_rem_to(const int_poly B)
613{
614    this->poly_pseudodiv_rem(*this,B);
615}
616
617
618//Pseudodivision nach Kaplan, M. Computeralgebra 4.6 welche q^e*A=Q*B+R
619//berechnet und entsprechendes in Q und R hineinschreibt
620
621int_poly int_poly::poly_pseudodiv(int_poly &Q, int_poly &R, int_poly A,  int_poly B)
622{
623
624    if (B.is_zero()==0)
625    {
626        // Initialisierungen: Vergiss zunÀchst die Hauptnenner von A und B (--> R bzw. Bint)
627        int_poly temp;
628        R.poly_set(A);
629
630
631        int k = A.deg - B.deg;
632
633        //Initialisiere Q mit 0en
634        Q.deg=k;
635        for (int i=0;i<=k;i++)
636        {
637            mpz_init_set_ui(Q.coef[i],0);
638        }
639
640
641        // Algorithmus
642        while (k>=0)
643        {
644
645            mpz_set(Q.coef[k],R.coef[R.deg]);
646
647            temp.poly_mon_mult(B,k);
648            temp.poly_scalar_mult_to(R.coef[R.deg]);
649
650            R.poly_scalar_mult_to(B.coef[B.deg]);
651            R.poly_sub_to(temp);
652
653            k=R.deg-B.deg;
654        }
655
656        int delta;
657        delta=0;
658        mpz_t dummy;
659        mpz_init_set_ui(dummy,0);
660
661        for (int i=0;i<=A.deg-B.deg;i++)
662        {
663            if (mpz_cmp_ui(Q.coef[i],0)!=0)
664            {
665                mpz_pow_ui(dummy,B.coef[B.deg],delta);
666                mpz_mul(Q.coef[i],Q.coef[i],dummy);
667                delta++;
668            }
669
670        }
671
672    }
673
674
675}
676
677
678//To Variante Algo 3.1.2 nach Cohen
679
680int_poly int_poly::poly_pseudodiv_to(int_poly &Q, int_poly &R, int_poly B)
681{
682    this->poly_pseudodiv(Q, R,*this,B);
683}
684
685// Kombinationen
686
687// a := a*b + c
688int_poly int_poly::poly_multadd_to(const int_poly b, const int_poly c)
689{
690    poly_mult_n_to(b);
691    poly_add_to(c);
692}
693
694//a=a*b-c
695int_poly int_poly::poly_multsub_to(const int_poly b, const int_poly c)
696{
697    poly_mult_n_to(b);
698    poly_sub_to(c);
699}
700
701
702
703/*
704// a := (a+b)* c
705int_poly int_poly::poly_addmult_to(const int_poly b, const int_poly c)
706{
707        int_poly a(deg,coef);
708        a.poly_add_to(b);
709        a.poly_mult_n_to(c);
710        poly_set(a);
711}
712*/
713
714// Eigenschaften
715
716// Content (Cohen 3.2.7), schreibt Ergebnis ins Argument cont
717void int_poly::poly_cont(mpz_t& cont)
718{
719    if (is_zero()==1)
720    {
721        mpz_init_set_ui(cont,0);
722    }
723    else
724    {
725        mpz_t temp;
726        int i=1;
727        mpz_init_set(temp,coef[0]);
728        while (mpz_cmp_ui(temp,1)!=0 && i<=deg)
729        {
730            mpz_gcd(temp,temp,coef[i]);
731            i++;
732        }
733        mpz_init_set(cont,temp);
734    }
735}
736
737
738// Primitive Part (Cohen 3.2.7) unter Verwendung von mpz_divexact
739// da wir wissen,dass der Inhalt alle Elemente teilt
740//ÜBERSCHREIBT DAS int_polyNOM WELCHES DEN BEFEHL AUFRUFT!!!!
741
742
743int_poly int_poly::poly_pp(int_poly f)
744{
745    mpz_t cont;
746    f.poly_cont(cont); // cont ist auf den Inhalt von f gesetzt.
747
748    if (mpz_cmp_ui(cont,1)==0)
749        poly_set(f);
750    else
751    {
752        deg=f.deg;
753        for (int i=0;i<=deg;i++)
754        {
755            mpz_init_set_ui(coef[i],0);
756            mpz_divexact(coef[i],f.coef[i],cont);
757        }
758
759    }
760}
761
762
763
764//Sonstiges
765void int_poly::poly_horner(mpz_t erg, const mpz_t u)
766{
767    mpz_init_set(erg,coef[deg]);
768    for (int i=deg;i>=1;i--)
769    {
770        mpz_mul(erg,erg,u);
771        mpz_add(erg,erg,coef[i-1]);
772    }
773}
774
775// int_polynom in int_polynom einsetzen(Horner-Schema) KRITISCHE EINGABE x^2, x^2  ....
776
777void int_poly::poly_horner_int_poly(const int_poly A,const int_poly B)
778{
779    poly_set(A.coef[A.deg]);
780    for (int i=A.deg;i>=1;i--)
781    {
782        poly_mult_n_to(B);
783        poly_add_const_to(A.coef[i-1]);
784    }
785}
786
787
788
789//Hilfsfunktionen
790
791
792//setze int_polynom auf int_polynom b
793int_poly int_poly::poly_set(const int_poly b)
794{
795    deg=b.deg;
796    for(int i=0;i<=deg;i++)
797    {
798        mpz_init_set(coef[i],b.coef[i]); // Hier wird init set dringendst benötigt
799    }
800
801}
802
803// setze int_polynom auf konstantes int_polynom b
804int_poly int_poly::poly_set(const mpz_t b)
805{
806    deg=0;
807    mpz_init_set(coef[0],b);
808}
809
810
811//setze int_polynom auf Nullint_polynom
812int_poly int_poly::poly_set_zero()
813{
814    deg = -1;
815}
816
817
818//Vergleiche ob 2 int_polynome gleich return 1 falls ja sont 0
819
820int int_poly::is_equal(const int_poly g)
821{
822    if (deg!=g.deg)
823        return 0;
824    else
825
826        for (int i=deg;i>=0; i--)
827        {
828        if (mpz_cmp(coef[i],g.coef[i])!=0)
829        {return 0;}
830    }
831    return 1;
832}
833
834//ÜberprÃŒft ob das int_polynom 0 ist
835
836int int_poly::is_zero()
837{
838    if (deg<0)
839        return 1;
840    else
841        return 0;
842
843}
844
845int int_poly::is_one()
846{
847    if (deg==0)
848    {
849        if (mpz_cmpabs_ui(coef[0],1)==0) { return 1; }
850        else { return 0; }
851    }
852    else { return 0; }
853}
854
855int int_poly::is_monic()
856{
857    if (mpz_cmpabs_ui(coef[deg],1)==0)
858        return 1;
859    else
860        return 0;
861}
862
863// klassischer GGT nach Cohen 3.2.1
864
865int_poly int_poly::poly_gcd( int_poly A,  int_poly B)
866{
867    if (A.deg<B.deg)
868        poly_gcd(B,A);
869    else if (B.is_zero()==1)
870        poly_set(A);
871    else if (B.is_monic()==0)
872    {
873        //cout << "Subresultanten GGT wird benutzt"<<endl;
874        poly_subgcd(A,B);
875    }
876    else
877    {
878        int_poly Q;
879        int_poly App;
880        int_poly Bpp;
881        int_poly R;
882        App.poly_set(A);
883        Bpp.poly_set(B);
884
885        while (Bpp.is_zero()==0)
886        {
887            R.poly_div(Q,R,App,Bpp);
888            App.poly_set(Bpp);
889            Bpp.poly_set(R);
890        }
891        poly_set(App);
892    }
893
894}
895
896// GGT nach Cohen, Algorithmus 3.2.10 (Primitive int_polynomial GCD) TO DO: Optimierung bzgl. Mehrfachberechnung)
897// Bpp ist das B in den Schritten ab 2
898
899
900int_poly int_poly::poly_ppgcd(int_poly A,int_poly B)
901{
902    if(A.deg<B.deg)
903    {
904        poly_ppgcd(B,A);
905
906    }
907    else if(B.is_zero()==1)
908    {
909        poly_set(A);
910
911    }
912    else
913    {
914        //Initialisierung und Reduktionen
915        mpz_t a;
916        mpz_init_set_ui(a,0);
917        mpz_t b;
918        mpz_init_set_ui(b,0);
919        mpz_t d;
920        mpz_init_set_ui(d,0);
921        A.poly_cont(a);
922        B.poly_cont(b);
923        mpz_gcd(d,a,b);
924
925        int_poly App;
926        int_poly Bpp;
927        int_poly R;
928
929        //Erster Schritt im Algo
930        App.poly_pp(A);
931        Bpp.poly_pp(B);
932        R.poly_pseudodiv_rem(App,Bpp);
933
934        //Algo
935
936        while (R.deg>0)
937        {
938            App.poly_set(Bpp);
939            Bpp.poly_pp(R);
940            R.poly_pseudodiv_rem(App,Bpp);
941        }
942
943        if (R.is_zero()==0)
944        {
945            deg=0;
946            mpz_init_set(coef[0],d);
947        }
948        else
949        {
950            poly_set(Bpp);
951            poly_scalar_mult_to(d);
952        }
953    }
954}
955// To Variante ppgcd
956
957
958int_poly int_poly::poly_ppgcd_to(int_poly B)
959{
960    this->poly_ppgcd(*this,B);
961}
962
963
964
965// GGT nach Cohen, Algorithmus 3.3.1 (Subresultant int_polynomial GCD) TO DO: Optimierung bzgl. Mehrfachberechnung)
966// Bpp ist das B in den Schritten ab 2
967int_poly int_poly::poly_subgcd(int_poly A, int_poly B)
968{
969    //Initialisierung und Reduktionen
970    if(A.deg<B.deg)
971    {
972        poly_subgcd(B,A);
973
974    }
975    else if(B.is_zero()==1)
976    {
977        poly_set(A);
978
979    }
980    else
981    {
982        mpz_t a;
983        mpz_init_set_ui(a,0);
984        mpz_t b;
985        mpz_init_set_ui(b,0);
986        mpz_t d;
987        mpz_init_set_ui(d,0);
988        mpz_t h;
989        mpz_init_set_ui(h,1);
990        mpz_t g;
991        mpz_init_set_ui(g,1);
992        mpz_t temp1;
993        mpz_init_set_ui(temp1,0);
994        mpz_t temp2;
995        mpz_init_set_ui(temp2,0);
996
997        A.poly_cont(a);
998        B.poly_cont(b);
999        mpz_gcd(d,a,b);
1000
1001        int_poly App;
1002        int_poly Bpp;
1003        int_poly R;
1004
1005        //Erster Schritt im Algo
1006        int delta;
1007        App.poly_pp(A);
1008        Bpp.poly_pp(B);
1009        R.poly_pseudodiv_rem(App,Bpp);
1010
1011        //Algo
1012
1013        while (R.deg>0)
1014        {
1015            delta =App.deg-Bpp.deg;
1016
1017            mpz_pow_ui(temp1,h,delta);
1018            mpz_mul(temp2,temp1,g);
1019            App.poly_set(Bpp);
1020            Bpp.poly_pp(R);
1021            Bpp.poly_scalar_div_to(temp2);
1022
1023            mpz_set(g,App.coef[App.deg]);
1024            mpz_pow_ui(temp1,h,1-delta);
1025            mpz_pow_ui(temp2,g,delta);
1026            mpz_mul(h,temp1,temp2);
1027
1028
1029            App.poly_set(Bpp);
1030            Bpp.poly_pp(R);
1031            R.poly_pseudodiv_rem(App,Bpp);
1032
1033        }
1034
1035        if (R.is_zero()==0)
1036        {
1037            deg=0;
1038            mpz_init_set(coef[0],d);
1039        }
1040        else
1041        {
1042            poly_set(Bpp);
1043            poly_cont(temp1);
1044            poly_scalar_mult_to(d);
1045            poly_scalar_div_to(temp1);
1046        }
1047    }
1048}
1049
1050// To Varianta Subresultanten
1051
1052int_poly int_poly::poly_subgcd_to(int_poly B)
1053{
1054    this->poly_subgcd(*this,B);
1055}
1056
1057
1058//Extended Subresultant GCD; see Kaplan, M. Computeralgebra, chapter 4.6
1059//returns g=r*A+t*B IT WORKS DONT TOUCH IT!!!!!!!!
1060int_poly int_poly::poly_extsubgcd(int_poly& r, int_poly& t,int_poly &g,int_poly A,int_poly B)
1061{
1062    if (A.deg<B.deg)
1063        poly_extsubgcd(t,r,g,B,A);
1064    else if (B.is_zero()==1)
1065    {
1066        g.poly_set(A);
1067        t.poly_set_zero();
1068
1069        mpz_t temp;
1070        mpz_init_set_ui(temp,1);
1071        r.poly_set(temp);
1072    }
1073
1074    else
1075    {
1076        //Initialization (temp will be -1 in the algorithm)
1077        mpz_t temp;
1078        mpz_t temp2;
1079        mpz_t psi;
1080        int alpha;
1081        int delta;
1082        int delta2;
1083        mpz_t base; //will be always (-1)^ ...
1084        mpz_t base2; //will be always (-1)^ .../LK ...
1085        mpz_t base3;
1086        mpz_init_set_ui(temp,1);
1087        mpz_init_set_ui(base,1);
1088        mpz_init_set_ui(base2,1);
1089        mpz_init_set_ui(base3,1);
1090        mpz_init_set_ui(psi,1);
1091        delta=A.deg-B.deg;
1092        delta2=delta;
1093        alpha=delta2+1;
1094
1095        int_poly dummy; // is needed in the main algorithm for -q*r_i resp. -q*t_i
1096        dummy.poly_set_zero();
1097
1098        int_poly dummy2; // is needed in the main algorithm for LK * r_i resp LK* t_i
1099        dummy2.poly_set_zero();
1100
1101        int_poly f1;
1102        int_poly f2;
1103        int_poly f3;
1104        int_poly f;
1105
1106        int_poly q;
1107
1108        int_poly r1;
1109        int_poly r2;
1110        int_poly r3;
1111
1112        int_poly t1;
1113        int_poly t2;
1114        int_poly t3;
1115
1116        f1.poly_set(A);
1117        f2.poly_set(B);
1118        f.poly_set_zero();
1119
1120        r1.poly_set(temp);
1121        r2.poly_set_zero();
1122
1123        t1.poly_set_zero();
1124        t2.poly_set(temp);
1125
1126        mpz_set_si(temp,-1);
1127
1128        //Calculating first step
1129        mpz_init_set_ui(temp2,0);
1130        mpz_pow_ui(temp2,f2.coef[f2.deg],alpha);
1131        f.poly_scalar_mult(f1,temp2);
1132
1133
1134        A.poly_div(q,f3,f,f2);
1135        mpz_pow_ui(base,temp,alpha);
1136        f3.poly_scalar_mult_to(base);
1137
1138
1139        r3.poly_set(base);
1140        mpz_pow_ui(base2,f2.coef[f2.deg],alpha);
1141        r3.poly_scalar_mult_to(base2);
1142
1143
1144        mpz_pow_ui(base,temp,delta);
1145        t3.poly_set(base);
1146        t3.poly_mult_n_to(q);
1147
1148
1149
1150        //Correcting the polynomials and constants
1151
1152        f1.poly_set(f2);
1153        f2.poly_set(f3);
1154
1155        r1.poly_set(r2);
1156        r2.poly_set(r3);
1157
1158        t1.poly_set(t2);
1159        t2.poly_set(t3);
1160
1161        delta=delta2;
1162        delta2=f1.deg-f2.deg;
1163        alpha=delta2+1;
1164
1165        //Main Algorithm
1166        while (f2.is_zero()==0)
1167        {
1168
1169
1170            //Step 1.1+1.2: base and base 2 will be psi^ ... and LK^...
1171
1172            mpz_pow_ui(temp2,f2.coef[f2.deg],alpha);
1173            f.poly_scalar_mult(f1,temp2);
1174            A.poly_div(q,f3,f,f2);
1175
1176
1177            mpz_pow_ui(base,psi,delta);
1178            mpz_pow_ui(base2,f1.coef[f1.deg],delta);
1179
1180
1181            mpz_mul(base2,base2,psi);
1182            mpz_divexact(psi,base2,base);
1183
1184            //Step 1.3
1185
1186            mpz_pow_ui(base,temp,alpha);
1187            mpz_pow_ui(base2,psi,delta2);
1188            mpz_mul(base2,base2,f1.coef[f1.deg]);
1189            f3.poly_scalar_div_to(base2);
1190            f3.poly_scalar_mult_to(base);
1191
1192
1193            //Step 1.4
1194
1195            mpz_pow_ui(base3,f2.coef[f2.deg],alpha);
1196
1197            //computing r_i
1198            dummy.poly_mult_n(q,r2);
1199            dummy2.poly_scalar_mult(r1,base3);
1200            r3.poly_sub(dummy2,dummy);
1201            r3.poly_scalar_mult_to(base);
1202            r3.poly_scalar_div_to(base2);
1203
1204            //computing t_i
1205            dummy.poly_mult_n(q,t2);
1206            dummy2.poly_scalar_mult(t1,base3);
1207            t3.poly_sub(dummy2,dummy);
1208            t3.poly_scalar_mult_to(base);
1209            t3.poly_scalar_div_to(base2);
1210
1211            //Correcting the polynomials and constants
1212
1213            f1.poly_set(f2);
1214            f2.poly_set(f3);
1215
1216            r1.poly_set(r2);
1217            r2.poly_set(r3);
1218
1219            t1.poly_set(t2);
1220            t2.poly_set(t3);
1221
1222            delta=delta2;
1223            delta2=f1.deg-f2.deg;
1224            alpha=delta2+1;
1225
1226        }
1227
1228        //Computing result
1229        g.poly_set(f1);
1230        r.poly_set(r1);
1231        t.poly_set(t1);
1232
1233    }
1234
1235
1236}
1237
1238//Ein & Ausgabe
1239
1240//Eingabe
1241
1242void int_poly::poly_insert()
1243{
1244#if 0
1245    cout << "Bitte geben Sie ein int_polynom ein! ZunÀchst den Grad: " << endl;
1246    cin >> deg;
1247
1248
1249    for ( int i=0; i<=deg;i++)
1250    {
1251        mpz_init_set_ui(coef[i],0);
1252        printf("Geben Sie nun f[%i] ein:",i);
1253        mpz_inp_str(coef[i],stdin, 10);
1254    }
1255#endif
1256}
1257
1258
1259//Ausgabe
1260void int_poly::poly_print()
1261{
1262#if 0
1263    if (is_zero()==1)
1264        cout << "0" << "\n" <<endl;
1265    else
1266    {
1267        for (int i=deg;i>=1;i--)
1268        {
1269            mpz_out_str(stdout,10, coef[i]);
1270            printf("X%i+",i);
1271        }
1272        mpz_out_str(stdout,10, coef[0]);
1273        printf("\n");
1274    }
1275#endif
1276}
1277
Note: See TracBrowser for help on using the repository browser.