source: git/libpolys/coeffs/AEp.cc @ 0acf3e

jengelh-datetimespielwiese
Last change on this file since 0acf3e was 4154bb, checked in by Oleksandr Motsak <motsak@…>, 9 years ago
Fix "if (deg=0)" in new exp. coeffs (see issue #374)
  • Property mode set to 100755
File size: 20.5 KB
Line 
1#include <stdio.h>
2#include <gmp.h>
3#include <math.h>
4#include "AEp.h"
5
6
7
8
9using namespace std;
10
11//Konstruktoren
12
13p_poly::p_poly()
14{
15    //coef = reinterpret_cast<mpz_t*> (omAlloc(100*sizeof(mpz_t)));
16    deg=-1;
17    mod=2;
18    mpz_init_set_ui(coef[0],0);
19}
20
21
22
23p_poly::p_poly(int n,int p, mpz_t *a)
24{
25
26    deg=n;
27    mod=p;
28
29    for ( int i=0;i<=n;i++)
30    {
31        mpz_mod_ui(a[i],a[i],mod);
32        mpz_init_set(coef[i], a[i]);
33    }
34
35}
36
37/*
38//Destruktor
39
40p_poly::~p_poly()
41{
42        delete[] coef;
43}
44
45*/
46
47//Reduktion modulo p
48
49p_poly p_poly::p_poly_reduce(p_poly f,int p)
50{
51    if (f.is_zero()==0)
52    {
53
54        for (int i=f.deg;i>=0;i--)
55        {
56            mpz_mod_ui(coef[i],f.coef[i],p);
57        }
58    }
59    //Hier nötige Grad Korrektur
60    int k=deg;
61    while(mpz_sgn(coef[k])==0 && k>=0)
62    {deg--;k--;}
63}
64
65
66// Arithmetik
67
68
69//Additionen
70
71//Standard - Addition
72p_poly p_poly::p_poly_add(const p_poly a, const p_poly b)
73{
74    if (a.deg >=b.deg)
75    {
76
77        deg=a.deg;
78        mod=a.mod;
79
80        for ( int i=0;i<=b.deg;i++)
81        {
82            mpz_add(coef[i],a.coef[i],b.coef[i]);
83        }
84
85        for ( int i=b.deg+1;i<=a.deg;i++)
86        {
87            mpz_init_set(coef[i],a.coef[i]);
88        }
89
90        //Hier nötige Grad Korrektur
91        int k=deg;
92        while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
93        {deg--;k--;}
94
95
96    }
97
98    else {p_poly_add(b,a);  }
99
100}
101
102//Überschreibende Addition
103
104p_poly p_poly::p_poly_add_to(const p_poly g)
105{
106    this->p_poly_add(*this,g);
107}
108
109//Addition einer Konstanten
110p_poly p_poly::p_poly_add_const(p_poly f,const mpz_t a)
111{
112    if (f.is_zero()==1 && mpz_divisible_ui_p(a,f.mod)==0)
113        p_poly_set(a,f.mod);
114
115    else if (mpz_divisible_ui_p(a,f.mod)==0)
116    {
117        p_poly_set(f);
118        mpz_add(coef[0],coef[0],a);
119        /*/Hier nötige Grad Korrektur
120        int k=deg;
121        while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
122        {deg--;k--;}
123        */
124    }
125
126}
127
128
129//To Variante Addition einer Konstanten
130
131p_poly p_poly::p_poly_add_const_to(const mpz_t a)
132{
133    this->p_poly_add_const(*this,a);
134}
135
136//Monom Addition
137p_poly p_poly::p_poly_add_mon(const p_poly f, mpz_t a, int i)
138{
139    p_poly_set(f);
140    if (i<=deg  && is_zero()==0)
141    {
142        mpz_add(coef[i],coef[i],a);
143    }
144    else if (is_zero()==1 && mpz_divisible_ui_p(a,f.mod)==0)
145    {
146        deg=i;
147        for(int j=0;j<=i;j++)
148        {
149            mpz_init_set_ui(coef[j],0);
150        }
151        mpz_t temp;
152        mpz_init_set_ui(temp,0);
153        mpz_mod_ui(temp,a,mod);
154        mpz_add(coef[i],coef[i],temp);
155
156    }
157    else if(i>deg && mpz_divisible_ui_p(a,f.mod)==0)
158    {
159        deg=i;
160        for(int j=i;j>deg;j--)
161        {
162            mpz_init_set_ui(coef[j],0);
163        }
164        mpz_t temp;
165        mpz_init_set_ui(temp,0);
166        mpz_mod_ui(temp,a,mod);
167        mpz_add(coef[i],coef[i],temp);
168
169    }
170    /*/Hier nötige Grad Korrektur
171    int k=deg;
172    while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
173    {deg--;k--;}
174    */
175
176}
177
178//To Variante Monomaddition
179p_poly p_poly::p_poly_add_mon_to(mpz_t a, int i)
180{
181
182    if (i<=deg  && is_zero()==0)
183    {
184        mpz_add(coef[i],coef[i],a);
185    }
186    else if (is_zero()==1 && mpz_divisible_ui_p(a,mod)==0)
187    {
188        deg=i;
189        for(int j=0;j<=i;j++)
190        {
191            mpz_init_set_ui(coef[j],0);
192        }
193        mpz_t temp;
194        mpz_init_set_ui(temp,0);
195        mpz_mod_ui(temp,a,mod);
196        mpz_add(coef[i],coef[i],temp);
197
198    }
199    else if(i>deg && mpz_divisible_ui_p(a,mod)==0)
200    {
201        deg=i;
202        for(int j=i;j>deg;j--)
203        {
204            mpz_init_set_ui(coef[j],0);
205        }
206        mpz_t temp;
207        mpz_init_set_ui(temp,0);
208        mpz_mod_ui(temp,a,mod);
209        mpz_add(coef[i],coef[i],temp);
210
211    }
212    /*Hier nötige Grad Korrektur
213    int k=deg;
214    while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
215    {deg--;k--;} */
216}
217
218//Subtraktionen
219
220p_poly p_poly::p_poly_sub(const p_poly a, const p_poly b)
221{
222    if (a.deg >=b.deg)
223    {
224
225        deg=a.deg;
226        mod=a.mod;
227
228        for ( int i=0;i<=b.deg;i++)
229        {
230            mpz_sub(coef[i],a.coef[i],b.coef[i]);
231        }
232
233        for ( int i=b.deg+1;i<=a.deg;i++)
234        {
235            mpz_init_set(coef[i],a.coef[i]);
236        }
237
238        //Hier nötige Grad Korrektur
239        int k=deg;
240        while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
241        {deg--;k--;}
242
243    }
244
245    else {p_poly_sub(b,a);p_poly_neg();  }
246
247}
248
249
250//Überschreibende Subtraktion
251
252p_poly p_poly::p_poly_sub_to(const p_poly b)
253{
254    this->p_poly_sub(*this,b);
255}
256
257//Subtraktion einer Konstanten
258p_poly p_poly::p_poly_sub_const(p_poly f,const mpz_t a)
259{
260    if (f.is_zero()==1)
261    {
262        p_poly_set(a,f.mod);
263        p_poly_neg();
264    }
265    else
266    {
267        p_poly_set(f);
268        mpz_sub(coef[0],coef[0],a);
269    }
270    /*/*Hier nötige Grad Korrektur
271    int k=deg;
272    while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
273    {deg--;k--;} */
274
275}
276
277
278//To Variante Subtraktion einer Konstanten
279
280p_poly p_poly::p_poly_sub_const_to(const mpz_t a)
281{
282    this->p_poly_sub_const(*this,a);
283}
284
285
286//Monom Subtraktion
287p_poly p_poly::p_poly_sub_mon(const p_poly f , mpz_t a, int i)
288{
289    mpz_t temp;
290    mpz_neg(temp,a);
291    p_poly_add_mon(f,temp,i);
292}
293
294//To Variante Monomaddition
295p_poly p_poly::p_poly_sub_mon_to(mpz_t a, int i)
296{
297    mpz_t temp;
298    mpz_neg(temp,a);
299    p_poly_add_mon_to(temp,i);
300}
301
302
303//Multiplikationen
304
305//Multiplikation mit Monom
306p_poly p_poly::p_poly_mon_mult( p_poly f, int n)
307{
308    if (f.is_zero()==1)
309    {p_poly_set_zero();}
310    else
311    {
312        deg=f.deg+n;
313        mod=f.mod;
314        for (int i=deg;i>=n;i--)
315        {
316            mpz_init_set(coef[i],f.coef[i-n]);
317        }
318        for (int i=n-1;i>=0;i--)
319        {
320            mpz_init_set_ui(coef[i],0);
321        }
322
323        /*/Hier nötige Grad Korrektur
324    int k=deg;
325    while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
326    {deg--;k--;} */
327    }
328}
329
330p_poly p_poly::p_poly_mon_mult_to(const int n)
331{
332    this->p_poly_mon_mult(*this,n);
333}
334
335
336//Multiplikation mit Skalar
337
338p_poly p_poly::p_poly_scalar_mult(const p_poly g, const mpz_t n)
339{
340    if (mpz_divisible_ui_p(n,g.mod)!=0)
341        p_poly_set_zero();
342    else
343    {
344        deg=g.deg;
345        mod=g.mod;
346
347        mpz_t temp;
348        mpz_init_set_ui(temp,0);
349        for(int i=0;i<=deg;i++)
350        {
351            mpz_mul(temp,n,g.coef[i]);
352            mpz_init_set(coef[i],temp);
353        }
354    }
355}
356
357
358
359p_poly p_poly::p_poly_scalar_mult(const mpz_t n, const p_poly g)
360{
361    if (mpz_divisible_ui_p(n,g.mod)!=0)
362        p_poly_set_zero();
363    else
364    {
365        deg=g.deg;
366        mod=g.mod;
367
368
369        mpz_t temp;
370        mpz_init_set_ui(temp,0);
371        for(int i=0;i<=deg;i++)
372        {
373            mpz_mul(temp,n,g.coef[i]);
374            mpz_init_set(coef[i],temp);
375        }
376        /*/Hier nötige Grad Korrektur
377        int k=deg;
378        while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
379        {deg--;k--;} */
380    }
381}
382
383
384p_poly p_poly::p_poly_scalar_mult_to(const mpz_t n)
385{
386    this->p_poly_scalar_mult(*this,n);
387}
388
389
390
391// Negation
392
393p_poly p_poly::p_poly_neg()
394{
395    for (int i=0;i<=deg;i++)
396    {
397        mpz_neg(coef[i],coef[i]);
398    }
399}
400
401// Naive Multiplikation
402p_poly p_poly::p_poly_mult_n(p_poly a,p_poly b)
403{
404    //Reduktion mod p
405    a.p_poly_reduce(a,a.mod);
406    b.p_poly_reduce(b,b.mod);
407
408    if (a.is_zero()==1 || b.is_zero()==1)
409        p_poly_set_zero();
410    else
411    {
412        mpz_t temp;
413        mpz_init_set_ui(temp,0);
414        deg = a.deg + b.deg;
415
416        // Kopien atemp und btemp von a bzw. b, mit Nullen aufgefÃŒllt
417        p_poly atemp, btemp;
418        atemp.p_poly_set(a);
419        btemp.p_poly_set(b);
420        for(int i=a.deg+1;i<=deg;i++)
421        {
422            mpz_init_set_ui(atemp.coef[i],0);
423        }
424        for(int i=b.deg+1;i<=deg;i++)
425        {
426            mpz_init_set_ui(btemp.coef[i],0);
427        }
428        atemp.deg = deg;
429        btemp.deg = deg;
430
431        // Multiplikationsalgorithmus
432        for (int k=0; k<=deg; k++)
433        {
434            mpz_init_set_ui(coef[k],0);        // k-ter Koeffizient zunÀchst 0
435            for (int i=0; i<=k; i++)        // dann schrittweise Summe der a[i]*b[k-i]/
436            {
437                mpz_mul(temp,atemp.coef[i],btemp.coef[k-i]);
438                mpz_add(coef[k],coef[k],temp);
439            }
440        }
441
442        /*/Hier nötige Grad Korrektur
443        int k=deg;
444        while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
445        {deg--;k--;} */
446
447
448    }
449}
450
451
452//Überschreibende Multiplikation
453
454p_poly p_poly::p_poly_mult_n_to(const p_poly g)
455{
456    this->p_poly_mult_n(*this,g);
457
458}
459
460// Karatsuba-Multiplikation (Weimerskirch/Paar Alg. 1), ACHTUNG VORLÄUFIGE VERSION, macht noch Fehler beim Grad und ist unelegant !!!
461p_poly p_poly::p_poly_mult_ka( p_poly A,  p_poly B)
462{
463
464    if (A.is_zero()==1 || B.is_zero()==1)
465    {
466        p_poly_set_zero();
467    }
468    else
469    {
470        //Reduktion mod p
471        A.p_poly_reduce(A,A.mod);
472        B.p_poly_reduce(B,B.mod);
473
474        // Größeren Grad feststellen
475        int n;
476        if(A.deg>=B.deg){n=A.deg+1;}
477        else{n=B.deg+1;}
478        // n auf nÀchste 2er-Potenz setzen (VORLÄUFIG!)
479        n = static_cast<int>(ceil(log(n)/log(2)));
480        n = static_cast<int>(pow(2,n));
481
482        if (n==1)
483        {
484            mpz_t AB;
485            mpz_mul(AB,A.coef[0],B.coef[0]);
486            p_poly_set(AB,A.mod);
487            this->p_poly_reduce(*this,A.mod);
488        }
489        else
490        {
491            // p_polynome A und B aufspalten
492            p_poly Au, Al, Bu, Bl;
493            Au.p_poly_mon_div(A,n/2);
494            Al.p_poly_mon_div_rem(A,n/2);
495            Bu.p_poly_mon_div(B,n/2);
496            Bl.p_poly_mon_div_rem(B,n/2);
497            p_poly Alu,Blu;
498            Alu.p_poly_add(Al,Au);
499            Blu.p_poly_add(Bl,Bu);
500
501            // Teile rekursiv multiplizieren
502            p_poly D0, D1, D01;
503            D0.p_poly_mult_ka(Al,Bl);
504            D1.p_poly_mult_ka(Au,Bu);
505            D01.p_poly_mult_ka(Alu,Blu);
506
507            // Ergebnis zusammensetzen
508            p_poly temp;
509            D01.p_poly_sub_to(D0);
510            D01.p_poly_sub_to(D1);
511            D01.p_poly_mon_mult_to(n/2);
512            D1.p_poly_mon_mult_to(n);
513            D1.p_poly_add_to(D01);
514            D1.p_poly_add_to(D0);
515            p_poly_set(D1);
516
517        }
518
519        //Hier nötige Grad Korrektur
520        int k=deg;
521        while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
522        {deg--;k--;}
523    }
524}
525
526
527
528//Skalare Divisionen
529
530p_poly p_poly::p_poly_scalar_div( const p_poly g, const mpz_t n)
531{
532
533    if (mpz_divisible_ui_p(n,g.mod)==0) // ÃŒberprÃŒft invertierbarkeit
534    {
535        deg=g.deg;
536        mod=g.mod;
537
538
539        mpz_t temp;
540        mpz_t p;
541        mpz_init_set_ui(temp,0);
542        mpz_init_set_ui(p,mod);
543        for(int i=0;i<=deg;i++)
544        {
545            mpz_invert(temp,n,p);
546            mpz_mul(temp,g.coef[i],temp);
547            mpz_init_set(coef[i],temp);
548        }
549
550        /*/Hier nötige Grad Korrektur
551        int k=deg;
552        while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
553        {deg--;k--;} */
554    }
555}
556
557
558p_poly p_poly::p_poly_scalar_div_to(const mpz_t n)
559{
560    this->p_poly_scalar_div(*this,n);
561}
562
563// Division durch Monom - Quotient
564p_poly p_poly::p_poly_mon_div(const p_poly f, const int n)
565{
566    if (f.deg<n)
567    {
568        p_poly_set_zero();
569    }
570    else
571    {
572        deg=f.deg-n;
573        mod=f.mod;
574
575        for (int i=0;i<=f.deg-n;i++)
576        {
577            mpz_init_set(coef[i],f.coef[n+i]);
578        }
579    }
580}
581
582// Division durch Monom - Rest
583p_poly p_poly::p_poly_mon_div_rem(const p_poly f, const int n)
584{
585    if (f.deg<n)
586    {
587        p_poly_set(f);
588    }
589    else
590    {
591        deg=n-1;
592        mod=f.mod;
593
594
595        for (int i=0;i<=n-1;i++)
596        {
597            mpz_init_set(coef[i],f.coef[i]);
598        }
599    }
600}
601
602
603
604
605//Euklidische Division nach Cohen Algo 3.1.1 (degA muss größer gleich deg B sein)!!
606
607p_poly p_poly::p_poly_div_rem( p_poly A,  p_poly B)
608{
609
610    if (B.is_zero()==0)
611    {
612        //Reduktion mod p
613        A.p_poly_reduce(A,A.mod);
614        B.p_poly_reduce(B,B.mod);
615
616        p_poly R;
617        p_poly temp;
618        R.p_poly_set(A);
619        mpz_t a;
620        mpz_t u;
621        mpz_t p;
622        mpz_init_set_ui(p,A.mod);
623        mpz_init_set_ui(a,0);
624        mpz_init_set_ui(u,0);
625        int i;
626
627        mpz_invert(u,B.coef[B.deg],p); // Inverse von lc(B)
628
629
630        while (R.deg>=B.deg)
631        {
632
633            mpz_mul(a,R.coef[R.deg],u);
634            i=R.deg-B.deg;
635
636            temp.p_poly_mon_mult(B,i);
637            temp.p_poly_scalar_mult_to(a);
638
639            R.p_poly_sub_to(temp);
640
641        }
642        p_poly_set(R);
643
644        /*/Hier nötige Grad Korrektur
645    int k=deg;
646    while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
647    {deg--;k--;}*/
648    }
649}
650//To Variante von Algo 3.1.1 im Cohen
651
652p_poly p_poly::p_poly_div_rem_to(const p_poly B)
653{
654    this->p_poly_div_rem(*this,B);
655
656
657}
658
659
660
661//Exakte Division nach Cohen 3.1.1
662p_poly p_poly::p_poly_div(p_poly &Q, p_poly &R, p_poly A,  p_poly B)
663{
664    if (B.is_zero()==0)
665    {
666        //Reduktion mod p
667        A.p_poly_reduce(A,A.mod);
668        B.p_poly_reduce(B,B.mod);
669
670        //Initialisierungen
671        p_poly temp;
672        R.p_poly_set(A);
673        Q.p_poly_set_zero();
674        Q.mod=A.mod;
675
676        mpz_t a;
677        mpz_t b;
678        mpz_t p;
679        mpz_init_set_ui(p,A.mod);
680        mpz_init_set_ui(a,0);
681        mpz_init_set_ui(b,0);
682        int i;
683        mpz_invert(b,B.coef[B.deg],p);
684
685        //Algorithmus TO DO: evtl hier mit auch den Rest ausgeben
686        while (R.deg>=B.deg)
687        {
688
689            mpz_mul(a,R.coef[R.deg],b);
690            i=R.deg-B.deg;
691
692            temp.p_poly_mon_mult(B,i);
693            temp.p_poly_scalar_mult_to(a);
694
695            R.p_poly_sub_to(temp);
696
697            Q.p_poly_add_mon_to(a,i);
698
699            R.p_poly_reduce(R,R.mod);
700            Q.p_poly_reduce(Q,Q.mod);
701        }
702
703        /*/Hier nötige Grad Korrektur Q
704    int k=Q.deg;
705    while(mpz_divisible_ui_p(Q.coef[k],Q.mod)!=0 && k>=0)
706    {Q.deg--;k--;}*/
707
708        /*/Hier nötige Grad Korrektur R
709    k=R.deg;
710    while(mpz_divisible_ui_p(R.coef[k],R.mod)!=0 && k>=0)
711    {R.deg--;k--;} */
712    }
713}
714
715
716//To Varainte der exakten Division
717
718p_poly p_poly::p_poly_div_to(p_poly &Q,p_poly &R, p_poly B)
719{
720    this->p_poly_div(Q ,R,*this,B);
721}
722
723
724// Kombinationen
725
726// a := a*b + c
727p_poly p_poly::p_poly_multadd_to(const p_poly b, const p_poly c)
728{
729    p_poly_mult_n_to(b);
730    p_poly_add_to(c);
731}
732
733//a=a*b-c
734p_poly p_poly::p_poly_multsub_to(const p_poly b, const p_poly c)
735{
736    p_poly_mult_n_to(b);
737    p_poly_sub_to(c);
738}
739
740
741
742/*
743// a := (a+b)* c
744p_poly p_poly::poly_addmult_to(const p_poly b, const p_poly c)
745{
746        p_poly a(deg,coef);
747        a.poly_add_to(b);
748        a.poly_mult_n_to(c);
749        poly_set(a);
750}
751*/
752
753
754
755//Sonstiges
756void p_poly::p_poly_horner(mpz_t erg, const mpz_t u)
757{
758    if (is_zero()==0)
759    {
760        mpz_init_set(erg,coef[deg]);
761        for (int i=deg;i>=1;i--)
762        {
763            mpz_mul(erg,erg,u);
764            mpz_add(erg,erg,coef[i-1]);
765        }
766
767        //Reduktion
768        mpz_mod_ui(erg,erg,mod);
769    }
770    else
771    {
772        mpz_set_ui(erg,0);
773    }
774}
775
776// p_polynom in p_polynom einsetzen(Horner-Schema) KRITISCHE EINGABE x^2, x^2  ....
777
778void p_poly::p_poly_horner_p_poly( p_poly A, p_poly B)
779{
780    //Reduktion mod p
781    A.p_poly_reduce(A,A.mod);
782    B.p_poly_reduce(B,B.mod);
783
784    p_poly_set(A.coef[A.deg],A.mod);
785    for (int i=A.deg;i>=1;i--)
786    {
787        p_poly_mult_n_to(B);
788        p_poly_add_const_to(A.coef[i-1]);
789    }
790    /*/Hier nötige Grad Korrektur
791    int i=deg;
792    while(mpz_divisible_ui_p(coef[i],mod)!=0 && i>=0)
793    {deg--;i--;} */
794}
795
796
797
798//Hilfsfunktionen
799
800
801//setze p_polynom auf p_polynom b
802p_poly p_poly::p_poly_set(const p_poly b)
803{
804    deg=b.deg;
805    mod=b.mod;
806
807
808    for(int i=0;i<=deg;i++)
809    {
810        mpz_init_set(coef[i],b.coef[i]); // Hier wird init set dringendst benötigt
811    }
812
813}
814
815// setze p_polynom auf konstantes p_polynom b
816p_poly p_poly::p_poly_set(const mpz_t b,int p)
817{
818    deg=0;
819    mod=p;
820
821    if (mpz_divisible_ui_p(b,mod)!=0)
822        p_poly_set_zero();
823    else
824    {
825        mpz_t temp;
826        mpz_init_set(temp,b);
827        mpz_mod_ui(temp,temp,p);
828        mpz_init_set(coef[0],b);
829    }
830}
831
832
833//setze p_polynom auf Nullpolynom
834p_poly p_poly::p_poly_set_zero()
835{
836    deg = -1;
837}
838
839
840//Vergleiche ob 2 p_polynome gleich return 1 falls ja sont 0
841
842int p_poly::is_equal(const p_poly g)
843{
844    if (deg!=g.deg)
845        return 0;
846    else
847
848        for (int i=deg;i>=0; i--)
849        {
850        if (mpz_cmp(coef[i],g.coef[i])!=0)
851        {return 0;}
852    }
853    return 1;
854}
855
856//ÜberprÃŒft ob das p_polynom 0 ist
857
858int p_poly::is_zero()
859{
860    if (deg<0)
861        return 1;
862    else
863        return 0;
864
865}
866
867int p_poly::is_one()
868{
869    if (deg==0)
870    {
871        if (mpz_cmp_ui(coef[0],1)==0) { return 1; }
872        else { return 0; }
873    }
874    else { return 0; }
875}
876
877int p_poly::is_monic()
878{
879    if (mpz_cmpabs_ui(coef[deg],1)==0)
880        return 1;
881    else
882        return 0;
883}
884
885// klassischer GGT nach Cohen 3.2.1
886
887p_poly p_poly::p_poly_gcd( p_poly A,  p_poly B)
888{
889
890    //Reduktion mod p
891    A.p_poly_reduce(A,A.mod);
892    B.p_poly_reduce(B,B.mod);
893
894    if (A.deg<B.deg)
895        p_poly_gcd(B,A);
896    else if (B.is_zero()==1)
897        p_poly_set(A);
898    else
899    {
900        p_poly App;
901        p_poly Bpp;
902        p_poly R;
903        App.p_poly_set(A);
904        Bpp.p_poly_set(B);
905
906        while (Bpp.is_zero()==0)
907        {
908            R.p_poly_div_rem(App,Bpp);
909            App.p_poly_set(Bpp);
910            Bpp.p_poly_set(R);
911        }
912        p_poly_set(App);
913    }
914
915}
916
917//Nach nach Fieker 2.12 Symbolisches Rechnen (2012)
918// gibt g=s*A+t*B aus
919p_poly p_poly::p_poly_extgcd(p_poly &s,p_poly &t,p_poly &g, p_poly A, p_poly B)
920{
921
922    //Reduktion mod p
923    A.p_poly_reduce(A,A.mod);
924    B.p_poly_reduce(B,B.mod);
925
926
927    if (A.deg<B.deg)
928        p_poly_extgcd(t,s,g,B,A);
929    else if (B.is_zero()==1)
930    {
931        g.p_poly_set(A);
932        t.p_poly_set_zero();
933
934        mpz_t temp;
935        mpz_init_set_ui(temp,1);
936        s.p_poly_set(temp,A.mod);
937    }
938
939    else
940    {
941        mpz_t temp;
942        mpz_init_set_ui(temp,1);
943
944        p_poly R1;
945        R1.p_poly_set(A);
946        p_poly R2;
947        R2.p_poly_set(B);
948        p_poly R3;
949        R3.mod=A.mod;
950
951
952        p_poly S1;
953        S1.p_poly_set(temp,A.mod);
954        p_poly S2;
955        S2.p_poly_set_zero();
956        S2.mod=A.mod;
957        p_poly S3;
958        S3.mod=A.mod;
959
960        p_poly T1;
961        T1.p_poly_set_zero();
962        T1.mod=A.mod;
963        p_poly T2;
964        T2.p_poly_set(temp,A.mod);
965        p_poly T3;
966        T3.mod=A.mod;
967
968        p_poly Q;
969
970        while (R2.is_zero()!=1)
971        {
972            p_poly_div(Q,R3,R1,R2);
973
974            S3.p_poly_mult_n(Q,S2);
975            S3.p_poly_neg();
976            S3.p_poly_add_to(S1);
977
978            T3.p_poly_mult_n(Q,T2);
979            T3.p_poly_neg();
980            T3.p_poly_add_to(T1);
981
982            R1.p_poly_set(R2);
983            R2.p_poly_set(R3);
984
985            S1.p_poly_set(S2);
986            S2.p_poly_set(S3);
987
988            T1.p_poly_set(T2);
989            T2.p_poly_set(T3);
990        }
991        t.p_poly_set(T1);
992        s.p_poly_set(S1);
993        g.p_poly_set(R1);
994    }
995}
996
997
998//Ein & Ausgabe
999
1000//Eingabe
1001
1002void p_poly::p_poly_insert()
1003{
1004#if 0
1005    cout << "Bitte geben Sie ein p_polynom ein! ZunÀchst den Grad: " << endl;
1006    cin >> deg;
1007    cout << "Jetzt den modul: " << endl;
1008    cin >> mod;
1009
1010    for ( int i=0; i<=deg;i++)
1011    {
1012        mpz_init_set_ui(coef[i],0);
1013        printf("Geben Sie nun f[%i] ein:",i);
1014        mpz_inp_str(coef[i],stdin, 10);
1015    }
1016    //Reduktion
1017    this->p_poly_reduce(*this,mod);
1018#endif
1019}
1020
1021
1022//Ausgabe
1023void p_poly::p_poly_print()
1024{
1025#if 0
1026
1027    //Reduktion
1028    this->p_poly_reduce(*this,mod);
1029
1030
1031    if (is_zero()==1)
1032        cout << "0" << "\n" <<endl;
1033    else
1034    {
1035        for (int i=deg;i>=1;i--)
1036        {
1037            mpz_out_str(stdout,10, coef[i]);
1038            printf("X%i+",i);
1039        }
1040        mpz_out_str(stdout,10, coef[0]);
1041        printf("\n");
1042    }
1043#endif
1044}
1045
Note: See TracBrowser for help on using the repository browser.