source: git/libpolys/coeffs/AEQ.cc @ 45cc512

spielwiese
Last change on this file since 45cc512 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: 18.4 KB
RevLine 
[0bfec5]1#include <stdio.h>
2#include <gmp.h>
3#include <math.h>
4#include "AEQ.h"
5
6
7
8using namespace std;
9
10//Konstruktoren
11
12Q_poly::Q_poly()
13{
14    deg=-1;
15    mpz_init_set_ui(denom,1);
16    mpz_init_set_ui(coef[0],0);
17}
18
19
20
21Q_poly::Q_poly(int n,mpz_t d, mpz_t *a)
22{
23    deg=n;
24
25    mpz_init_set(denom,d);
26
27    for ( int i=0;i<=n;i++)
28    {
29        mpz_init_set(coef[i], a[i]);
30    }
31}
32
33/*
34//Destruktor
35
36Q_poly::~Q_poly()
37{
38        delete[] coef;
39}
40
41*/
42
43// KÃŒrzen  --  MACHT NOCH MIST!
44Q_poly Q_poly::Q_poly_reduce()
45{
46    if (is_zero()==1)
47    {
48        mpz_init_set_ui(denom,1);
49    }
50    else
51    {
52        mpz_t d;
53        mpz_init_set(d,denom);
54        int i=0;
55        while (mpz_cmp_ui(d,1)!=0 && i<=deg)
56        {
57            mpz_gcd(d,d,coef[i]);
58            i++;
59        }
60        if (mpz_sgn(denom)==-1)
61        {
62            mpz_neg(d,d);
63        }
64        if (mpz_cmp_ui(d,1)!=0)
65        {
66            mpz_div(denom,denom,d);
67            for (int j=0;j<=deg;j++)
68            {
69                mpz_div(coef[j],coef[j],d);
70            }
71        }
72    }
73    // Grad-Korrektur
74    int j;
75    j=deg;
76    while(mpz_sgn(coef[j])==0 && j>=0)
77    {deg--;j--;}
78}
79
80// Koeffizienten mit b erweitern
81Q_poly Q_poly::Q_poly_extend(mpz_t b)
82{
83        mpz_mul(denom,denom,b);
84        for (int i=0;i<=deg;i++)
85        {
86                mpz_mul(coef[i],coef[i],b);
87        }
88}
89
90
91// Arithmetik
92
93
94//Additionen
95
96//Standard - Addition
97Q_poly Q_poly::Q_poly_add(const Q_poly a, const Q_poly b)
98{
99    if (a.deg >= b.deg)
100    {
101        deg=a.deg;
102        mpz_t atemp, btemp;
103        mpz_init_set_ui(atemp,0);
104        mpz_init_set_ui(btemp,0);
105
106        for (int i=0;i<=b.deg;i++)
107        {
108                mpz_mul(atemp,a.coef[i],b.denom);
109                mpz_mul(btemp,b.coef[i],a.denom);
110                mpz_add(coef[i],atemp,btemp);
111        }
112
113        for ( int i=b.deg+1;i<=a.deg;i++)
114        {
115            mpz_mul(coef[i],a.coef[i],b.denom);
116        }
117        mpz_mul(denom,a.denom,b.denom);
118
119        // Grad-Korrektur
120        int i=deg;
121        while(mpz_sgn(coef[i])==0 && i>=0)
122        {deg--;i--;}
123    }
124
125    else {Q_poly_add(b,a);}
126
127}
128
129//Überschreibende Addition
130
131Q_poly Q_poly::Q_poly_add_to(const Q_poly g)
132{
133    this->Q_poly_add(*this,g);
134}
135
136//Addition einer Konstanten
137Q_poly Q_poly::Q_poly_add_const(Q_poly f, const mpz_t a)
138{
139    if (f.is_zero()==1)
140    {
141        Q_poly_set(a,f.denom);
142    }
143    else
144    {
145        Q_poly_set(f);
146        mpz_t atemp;
147        mpz_mul(atemp,a,f.denom);
148        mpz_add(coef[0],coef[0],atemp);
149        // Grad Korrektur
150        if (deg==0 && mpz_sgn(coef[0])==0)
151            Q_poly_set_zero();
152    }
153}
154
155
156//To Variante Addition einer Konstanten
157
158Q_poly Q_poly::Q_poly_add_const_to(const mpz_t a)
159{
160    this->Q_poly_add_const(*this,a);
161}
162
163//Monom Addition
164Q_poly Q_poly::Q_poly_add_mon(const Q_poly f, mpz_t a, int i)
165{
166    Q_poly_set(f);
167    if (i<=deg  && is_zero()==0)
168    {
169        mpz_t atemp;
170        mpz_init_set_ui(atemp,0);
171        mpz_mul(atemp,a,f.denom);
172        mpz_add(coef[i],coef[i],atemp);
173
174        // Grad Korrektur
175
176        if (deg==i && mpz_sgn(coef[i])==0)
177        {deg--;}
178    }
179    else if (is_zero()==1)
180    {
181        deg=i;
182        for(int j=0;j<i;j++)
183        {
184            mpz_init_set_ui(coef[j],0);
185        }
186        mpz_init_set(coef[i],a);
187        mpz_init_set_ui(denom,1);
188    }
189    else if(i>deg)
190    {
191        deg=i;
192        for(int j=i-1;j>deg;j--)
193        {
194            mpz_init_set_ui(coef[j],0);
195        }
196        mpz_t atemp;
197        mpz_mul(atemp,a,f.denom);
198        mpz_init_set(coef[i],atemp);
199    }
200}
201
202//To Variante Monomaddition
203Q_poly Q_poly::Q_poly_add_mon_to(mpz_t a, int i)
204{
205    this->Q_poly_add_mon(*this,a,i);
206}
207
208//Subtraktionen
209
210Q_poly Q_poly::Q_poly_sub(const Q_poly a, const Q_poly b)
211{
212    Q_poly temp;
213    temp.Q_poly_set(b);
214    temp.Q_poly_neg();
215    Q_poly_add(a,temp);
216}
217
218
219//Überschreibende Subtraktion
220
221Q_poly Q_poly::Q_poly_sub_to(const Q_poly b)
222{
223    this->Q_poly_sub(*this,b);
224}
225
226//Subtraktion einer Konstanten
227Q_poly Q_poly::Q_poly_sub_const(Q_poly f,const mpz_t a)
228{
229    if (f.is_zero()==1)
230    {
231        Q_poly_set(a);
232        Q_poly_neg();
233    }
234    else
235    {
236        Q_poly_set(f);
237        mpz_t atemp;
238        mpz_init_set_ui(atemp,1);
239        mpz_mul(atemp,a,f.denom);
240        mpz_sub(coef[0],coef[0],atemp);
241    }
242}
243
244
245//To Variante Subtraktion einer Konstanten
246
247Q_poly Q_poly::Q_poly_sub_const_to(const mpz_t a)
248{
249    this->Q_poly_sub_const(*this,a);
250}
251
252
253//Monom Subtraktion
254Q_poly Q_poly::Q_poly_sub_mon(const Q_poly f , mpz_t a, int i)
255{
256    mpz_t temp;
257    mpz_init_set_ui(temp,0);
258    mpz_neg(temp,a);
259    Q_poly_add_mon(f,temp,i);
260}
261
262//To Variante Monomsubtraktion
263Q_poly Q_poly::Q_poly_sub_mon_to(mpz_t a, int i)
264{
265    this->Q_poly_sub_mon(*this,a,i);
266}
267
268
269//Multiplikationen
270
271//Multiplikation mit Monom
272Q_poly Q_poly::Q_poly_mon_mult(const Q_poly f, int n)
273{
274    deg=f.deg+n;
275    mpz_init_set(denom,f.denom);
276    for (int i=deg;i>=n;i--)
277    {
278        mpz_init_set(coef[i],f.coef[i-n]);
279    }
280    for (int i=n-1;i>=0;i--)
281    {
282        mpz_init_set_ui(coef[i],0);
283    }
284}
285
286Q_poly Q_poly::Q_poly_mon_mult_to(const int n)
287{
288    this->Q_poly_mon_mult(*this,n);
289}
290
291
292//Multiplikation mit Skalar
293
294Q_poly Q_poly::Q_poly_scalar_mult(const Q_poly g, const mpz_t n)
295{
296    deg=g.deg;
297    mpz_init_set(denom,g.denom);
298
299    mpz_t temp;
300    mpz_init_set_ui(temp,0);
301    for(int i=0;i<=deg;i++)
302    {
303        mpz_mul(temp,n,g.coef[i]);
304        mpz_init_set(coef[i],temp);
305    }
306}
307
308
309
310Q_poly Q_poly::Q_poly_scalar_mult(const mpz_t n, const Q_poly g)
311{
312    deg=g.deg;
313    mpz_init_set(denom,g.denom);
314
315    mpz_t temp;
316    mpz_init_set_ui(temp,0);
317    for(int i=0;i<=deg;i++)
318    {
319        mpz_mul(temp,n,g.coef[i]);
320        mpz_init_set(coef[i],temp);
321    }
322}
323
324
325Q_poly Q_poly::Q_poly_scalar_mult_to(const mpz_t n)
326{
327    this->Q_poly_scalar_mult(*this,n);
328}
329
330
331
332// Negation
333
334Q_poly Q_poly::Q_poly_neg()
335{
336    mpz_neg(denom,denom);
337}
338
339// Naive Multiplikation
340Q_poly Q_poly::Q_poly_mult_n(Q_poly a,Q_poly b)
341{
342
343    if (a.is_zero()==1 || b.is_zero()==1)
344        Q_poly_set_zero();
345    else
346    {
347        mpz_t temp;
348        mpz_init_set_ui(temp,0);
349        deg = a.deg + b.deg;
350
351        // Kopien atemp und btemp von a bzw. b, mit Nullen aufgefÃŒllt
352        Q_poly atemp, btemp;
353        atemp.Q_poly_set(a);
354        btemp.Q_poly_set(b);
355        for(int i=a.deg+1;i<=deg;i++)
356        {
357            mpz_init_set_ui(atemp.coef[i],0);
358        }
359        for(int i=b.deg+1;i<=deg;i++)
360        {
361            mpz_init_set_ui(btemp.coef[i],0);
362        }
363        atemp.deg = deg;
364        btemp.deg = deg;
365
366        // Multiplikationsalgorithmus
367        for (int k=0; k<=deg; k++)
368        {
369            mpz_init_set_ui(coef[k],0);        // k-ter Koeffizient zunÀchst 0
370            for (int i=0; i<=k; i++)        // dann schrittweise Summe der a[i]*b[k-i]/
371            {
372                mpz_mul(temp,atemp.coef[i],btemp.coef[k-i]);
373                mpz_add(coef[k],coef[k],temp);
374            }
375        }
376        mpz_mul(denom,a.denom,b.denom);
377    }
378}
379
380//Überschreibende Multiplikation
381
382Q_poly Q_poly::Q_poly_mult_n_to(const Q_poly g)
383{
384    this->Q_poly_mult_n(*this,g);
385}
386
387// Karatsuba-Multiplikation (Weimerskirch/Paar Alg. 1), ACHTUNG VORLÄUFIGE VERSION, macht noch Fehler beim Grad und ist unelegant !!!
388Q_poly Q_poly::Q_poly_mult_ka(const Q_poly A, const Q_poly B)
389{
390    // Größeren Grad feststellen
391    int n;
392    if(A.deg>=B.deg){n=A.deg+1;}
393    else{n=B.deg+1;}
394    // n auf nÀchste 2er-Potenz setzen (VORLÄUFIG!)
395    n = static_cast<int>(ceil(log(n)/log(2)));
396    n = static_cast<int>(pow(2,n));
397
398    if (n==1)
399    {
400        mpz_t AB;
401        mpz_mul(AB,A.coef[0],B.coef[0]);
402        Q_poly_set(AB,A.denom);
403    }
404    else
405    {
406        // Q_polynome A und B aufspalten
407        Q_poly Au, Al, Bu, Bl;
408        Au.Q_poly_mon_div(A,n/2);
409        Al.Q_poly_mon_div_rem(A,n/2);
410        Bu.Q_poly_mon_div(B,n/2);
411        Bl.Q_poly_mon_div_rem(B,n/2);
412        Q_poly Alu,Blu;
413        Alu.Q_poly_add(Al,Au);
414        Blu.Q_poly_add(Bl,Bu);
415
416        // Teile rekursiv multiplizieren
417        Q_poly D0, D1, D01;
418        D0.Q_poly_mult_ka(Al,Bl);
419        D1.Q_poly_mult_ka(Au,Bu);
420        D01.Q_poly_mult_ka(Alu,Blu);
421
422        // Ergebnis zusammensetzen
423        Q_poly temp;
424        D01.Q_poly_sub_to(D0);
425        D01.Q_poly_sub_to(D1);
426        D01.Q_poly_mon_mult_to(n/2);
427        D1.Q_poly_mon_mult_to(n);
428        D1.Q_poly_add_to(D01);
429        D1.Q_poly_add_to(D0);
430        Q_poly_set(D1);
431    }
432}
433
434
435
436//Skalare Divisionen
437
438Q_poly Q_poly::Q_poly_scalar_div(const Q_poly g, const mpz_t n)
439{
440    if (mpz_sgn(n)!=0) // ÃŒberprÃŒft Teilung durch 0
441    {
442        Q_poly_set(g);
443        mpz_mul(denom,g.denom,n);
444    }
445}
446
447
448Q_poly Q_poly::Q_poly_scalar_div_to(const mpz_t n)
449{
450    this->Q_poly_scalar_div(*this,n);
451}
452
453// Division durch Monom - Quotient
454Q_poly Q_poly::Q_poly_mon_div(const Q_poly f, const int n)
455{
456    if (f.deg<n)
457    {
458        Q_poly_set_zero();
459    }
460    else
461    {
462        deg=f.deg-n;
463        mpz_init_set(denom,f.denom);
464
465        for (int i=0;i<=f.deg-n;i++)
466        {
467            mpz_init_set(coef[i],f.coef[n+i]);
468        }
469    }
470}
471
472// Division durch Monom - Rest
473Q_poly Q_poly::Q_poly_mon_div_rem(const Q_poly f, const int n)
474{
475    if (f.deg<n)
476    {
477        Q_poly_set(f);
478    }
479    else
480    {
481        // Grad-Korrektur ist inklusive
482        deg=n-1;
483        int j=deg;
484        while(mpz_sgn(f.coef[j])==0 && j>=0)
485        {
486            deg--;
487            j--;
488            mpz_init_set_ui(coef[j],0);
489        }
490        for (int i=j;i>=0;i--)
491        {
492            mpz_init_set(coef[i],f.coef[i]);
493        }
494        mpz_init_set(denom,f.denom);
495    }
496}
497
498
499
500
501// Euklidische Division nach Cohen Algo 3.1.1 (degA muss größer gleich deg B sein)!!
502
503Q_poly Q_poly::Q_poly_div_rem(const Q_poly A, const Q_poly B)
504{
505
506    // Initialisierungen: Vergiss zunÀchst die Hauptnenner von A und B (--> R bzw. Bint)
507    Q_poly temp, Bint;
508    Q_poly_set(A);
509    mpz_init_set_ui(denom,1);
510    Bint.Q_poly_set(B);
511    mpz_init_set_ui(Bint.denom,1);
512    int e = A.deg - B.deg + 1;
513
514    // Algorithmus
515    while (deg>=B.deg)
516    {
517        temp.Q_poly_mon_mult(Bint,deg-B.deg);
518        temp.Q_poly_scalar_mult_to(coef[deg]);
519
520        Q_poly_scalar_mult_to(B.coef[B.deg]);
521        Q_poly_sub_to(temp);
522
523        e--;
524    }
525
526    // Terminierung
527    mpz_t d,q;
528    mpz_init_set(d,B.coef[B.deg]);
529    if (e>0)
530    {
531        mpz_pow_ui(q,d,e);
532        Q_poly_scalar_mult_to(q);
533    }
534    else if (e<0)
535    {
536        mpz_pow_ui(q,d,-e);
537        Q_poly_scalar_div_to(q);
538    }
539
540    mpz_pow_ui(d,d,A.deg-B.deg+1);
541    mpz_mul(denom,denom,d);
542
543    // Hauptnenner von A und B berÃŒcksichtigen
544    mpz_mul(denom,denom,A.denom);
545    Q_poly_scalar_mult_to(B.denom);
546}
547
548
549//To Variante von Algo 3.1.1 im Cohen
550
551Q_poly Q_poly::Q_poly_div_rem_to(const Q_poly B)
552{
553    this->Q_poly_div_rem(*this,B);
554}
555
556
557// Division nach Cohen 3.1.2 (gibt R und Q aus) --> FÃŒhrt Pseudo-Division durch, korrigiert den Faktor aber im Nenner
558Q_poly Q_poly::Q_poly_div(Q_poly &Q, Q_poly &R, const Q_poly A, const Q_poly B)
559{
560
561    // Initialisierungen: Vergiss zunÀchst die Hauptnenner von A und B (--> R bzw. Bint)
562    Q_poly temp, Bint;
563    R.Q_poly_set(A);
564    mpz_init_set_ui(R.denom,1);
565    Q.Q_poly_set_zero();
566    Bint.Q_poly_set(B);
567    mpz_init_set_ui(Bint.denom,1);
568    int e = A.deg - B.deg + 1;
569
570    // Algorithmus
571    while (R.deg>=B.deg)
572    {
573        temp.Q_poly_mon_mult(Bint,R.deg-B.deg);
574        temp.Q_poly_scalar_mult_to(R.coef[R.deg]);
575
576        Q.Q_poly_scalar_mult_to(B.coef[B.deg]);
577        Q.Q_poly_add_mon_to(R.coef[R.deg],R.deg-B.deg);
578
579        R.Q_poly_scalar_mult_to(B.coef[B.deg]);
580        R.Q_poly_sub_to(temp);
581
582        e--;
583    }
584
585    // Terminierung
586    mpz_t d,q;
587    mpz_init_set(d,B.coef[B.deg]);
588    if (e>0)
589    {
590        mpz_pow_ui(q,d,e);
591        R.Q_poly_scalar_mult_to(q);
592        Q.Q_poly_scalar_mult_to(q);
593    }
594    else if (e<0)
595    {
596        mpz_pow_ui(q,d,-e);
597        R.Q_poly_scalar_div_to(q);
598        Q.Q_poly_scalar_div_to(q);
599    }
600
601    mpz_pow_ui(d,d,A.deg-B.deg+1);
602    mpz_mul(R.denom,R.denom,d);
603    mpz_mul(Q.denom,Q.denom,d);
604
605    // Hauptnenner von A und B berÃŒcksichtigen
606    mpz_mul(R.denom,R.denom,A.denom);
607    mpz_mul(Q.denom,Q.denom,A.denom);
608    R.Q_poly_scalar_mult_to(B.denom);
609    Q.Q_poly_scalar_mult_to(B.denom);
610}
611
612
613//To Variante der exakten Division
614
615Q_poly Q_poly::Q_poly_div_to(Q_poly &Q,Q_poly &R,const Q_poly B)
616{
617    this->Q_poly_div(Q,R,*this,B);
618}
619
620
621// Kombinationen
622
623// a := a*b + c
624Q_poly Q_poly::Q_poly_multadd_to(const Q_poly b, const Q_poly c)
625{
626    Q_poly_mult_n_to(b);
627    Q_poly_add_to(c);
628}
629
630//a=a*b-c
631Q_poly Q_poly::Q_poly_multsub_to(const Q_poly b, const Q_poly c)
632{
633    Q_poly_mult_n_to(b);
634    Q_poly_sub_to(c);
635}
636
637
638
639/*
640// a := (a+b)* c
641Q_poly Q_poly::poly_addmult_to(const Q_poly b, const Q_poly c)
642{
643        Q_poly a(deg,coef);
644        a.poly_add_to(b);
645        a.poly_mult_n_to(c);
646        poly_set(a);
647}
648*/
649
650
651
652//Sonstiges
653void Q_poly::Q_poly_horner(mpz_t erg, const mpz_t u)
654{
655    mpz_init_set(erg,coef[deg]);
656    for (int i=deg;i>=1;i--)
657    {
658        mpz_mul(erg,erg,u);
659        mpz_add(erg,erg,coef[i-1]);
660    }
661
662// erg noch durch denom dividieren
663
664}
665
666// Q_polynom in Q_polynom einsetzen(Horner-Schema) KRITISCHE EINGABE x^2, x^2  ....
667
668void Q_poly::Q_poly_horner_Q_poly(const Q_poly A,const Q_poly B)
669{
670    Q_poly_set(A.coef[A.deg],A.denom);
671    for (int i=A.deg;i>=1;i--)
672    {
673        Q_poly_mult_n_to(B);
674        Q_poly_add_const_to(A.coef[i-1]);
675    }
676
677  // Nenner anpassen
678
679}
680
681
682
683//Hilfsfunktionen
684
685
686//setze Q_polynom auf Q_polynom b
687Q_poly Q_poly::Q_poly_set(const Q_poly b)
688{
689    deg=b.deg;
690    mpz_init_set(denom,b.denom);
691
692    for(int i=0;i<=deg;i++)
693    {
694        mpz_init_set(coef[i],b.coef[i]); // Hier wird init set dringendst benötigt
695    }
696}
697
698
699// setze Q_polynom auf konstantes Q_polynom b/d
700Q_poly Q_poly::Q_poly_set(const mpz_t b, const mpz_t d)
701{
702    deg=0;
703    mpz_init_set(denom,d);
704    mpz_init_set(coef[0],b);
705}
706
707// setze Q_polynom auf konstantes Z_polynom b
708Q_poly Q_poly::Q_poly_set(const mpz_t b)
709{
710    deg=0;
711    mpz_init_set_ui(denom,1);
712    mpz_init_set(coef[0],b);
713}
714
715
716//setze Q_polynom auf Nullpolynom
717Q_poly Q_poly::Q_poly_set_zero()
718{
719    deg = -1;
720}
721
722
723// Vergleiche ob zwei Q_polynome gleich --> return 1 falls ja sont 0
724
725int Q_poly::is_equal(Q_poly &g)
726{
727    if (deg!=g.deg)
728    {
729      return 0;
730    }
731    else
732    {
733      g.Q_poly_reduce();
734      Q_poly_reduce();
735      for (int i=deg;i>=0; i--)
736      {
737        if (mpz_cmp(coef[i],g.coef[i])!=0)
738          {return 0;}
739      }
740      return 1;
741    }
742}
743
744//ÜberprÃŒft ob das Q_polynom 0 ist
745int Q_poly::is_zero()
746{
747    if (deg<0)
748        return 1;
749    else
750        return 0;
751
752}
753
754
755//ÜberprÃŒft ob das Q_polynom 1 ist
756int Q_poly::is_one()
757{
[4154bb]758    if (deg==0)
[0bfec5]759    {
760        if (mpz_cmp(coef[0],denom)==0) { return 1; }
761        else { return 0; }
762    }
763    else { return 0; }
764}
765
766int Q_poly::is_monic()
767{
768    if (mpz_cmp(coef[deg],denom)==0)
769        return 1;
770    else
771        return 0;
772}
773
774// klassischer GGT nach Cohen 3.2.1
775
776Q_poly Q_poly::Q_poly_gcd(Q_poly A, Q_poly B)
777{
778
779    if (A.deg<B.deg)
780        Q_poly_gcd(B,A);
781    else if (B.is_zero()==1)
782        Q_poly_set(A);
783    else
784    {
785        Q_poly App;
786        Q_poly Bpp;
787        Q_poly R;
788        App.Q_poly_set(A);
789        Bpp.Q_poly_set(B);
790        mpz_init_set_ui(App.denom,1);
791        mpz_init_set_ui(Bpp.denom,1);
792
793        while (Bpp.is_zero()==0)
794        {
795            R.Q_poly_div_rem(App,Bpp);
796            App.Q_poly_set(Bpp);
797            Bpp.Q_poly_set(R);
798        }
799        mpz_init_set(App.denom,App.coef[App.deg]);
800        Q_poly_set(App);
801    }
802}
803
804
805// Nach nach Fieker 2.12 Symbolisches Rechnen (2012) MACHT PROBLEME
806// gibt g=s*A+t*B aus
807Q_poly Q_poly::Q_poly_extgcd(Q_poly &s,Q_poly &t,Q_poly &g, Q_poly A, Q_poly B)
808{
809    if (A.deg<B.deg)
810        Q_poly_extgcd(t,s,g,B,A);
811    else if (B.is_zero()==1)
812    {
813        g.Q_poly_set(A);
814        t.Q_poly_set_zero();
815
816        mpz_t temp;
817        mpz_init_set_ui(temp,1);
818        s.Q_poly_set(temp,A.denom);
819    }
820
821    else
822    {
823        mpz_t temp;
824        mpz_init_set_ui(temp,1);
825
826        Q_poly R1;
827        R1.Q_poly_set(A);
828        Q_poly R2;
829        R2.Q_poly_set(B);
830        Q_poly R3;
831
832        Q_poly S1;
833        S1.Q_poly_set(temp,A.denom);
834        Q_poly S2;
835        S2.Q_poly_set_zero();
836        Q_poly S3;
837
838        Q_poly T1;
839        T1.Q_poly_set_zero();
840        Q_poly T2;
841        T2.Q_poly_set(temp,A.denom);
842        Q_poly T3;
843
844        Q_poly Q;
845
846        while (R2.is_zero()!=1)
847        {
848            Q_poly_div(Q,R3,R1,R2);
849
850            S3.Q_poly_mult_n(Q,S2);
851            S3.Q_poly_neg();
852            S3.Q_poly_add_to(S1);
853
854            T3.Q_poly_mult_n(Q,T2);
855            T3.Q_poly_neg();
856            T3.Q_poly_add_to(T1);
857
858            R1.Q_poly_set(R2);
859            R2.Q_poly_set(R3);
860
861            S1.Q_poly_set(S2);
862            S2.Q_poly_set(S3);
863
864            T1.Q_poly_set(T2);
865            T2.Q_poly_set(T3);
866        }
867        t.Q_poly_set(T1);
868        s.Q_poly_set(S1);
869        g.Q_poly_set(R1);
870    }
871}
872
873
874//Ein & Ausgabe
875
876//Eingabe
877
878void Q_poly::Q_poly_insert()
879{
880#if 0
881    cout << "Bitte geben Sie ein Q_polynom ein! ZunÀchst den Grad: " << endl;
882    cin >> deg;
883    mpz_init_set_ui(denom,1);
884    cout << "Jetzt den Hauptnenner: " << endl;
885    mpz_inp_str(denom,stdin, 10);
886
887    for ( int i=0; i<=deg;i++)
888    {
889        mpz_init_set_ui(coef[i],0);
890        printf("Geben Sie nun f[%i] ein:",i);
891        mpz_inp_str(coef[i],stdin, 10);
892    }
893#endif
894}
895
896
897//Ausgabe
898void Q_poly::Q_poly_print()
899{
900#if 0
901    if (is_zero()==1)
902        cout << "0" << "\n" <<endl;
903    else
904    {
905        printf("(");
906        for (int i=deg;i>=1;i--)
907        {
908            mpz_out_str(stdout,10,coef[i]);
909            printf("X%i+",i);
910        }
911        mpz_out_str(stdout,10,coef[0]);
912        printf(")/");
913        mpz_out_str(stdout,10,denom);
914        printf("\n");
915    }
916#endif
917}
918
Note: See TracBrowser for help on using the repository browser.