1 | #include <stdio.h> |
---|
2 | #include <gmp.h> |
---|
3 | #include <math.h> |
---|
4 | #include "AEQ.h" |
---|
5 | |
---|
6 | |
---|
7 | |
---|
8 | using namespace std; |
---|
9 | |
---|
10 | //Konstruktoren |
---|
11 | |
---|
12 | Q_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 | |
---|
21 | Q_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 | |
---|
36 | Q_poly::~Q_poly() |
---|
37 | { |
---|
38 | delete[] coef; |
---|
39 | } |
---|
40 | |
---|
41 | */ |
---|
42 | |
---|
43 | // KÃŒrzen -- MACHT NOCH MIST! |
---|
44 | Q_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 |
---|
81 | Q_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 |
---|
97 | Q_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 | |
---|
131 | Q_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 |
---|
137 | Q_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 | |
---|
158 | Q_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 |
---|
164 | Q_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 |
---|
203 | Q_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 | |
---|
210 | Q_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 | |
---|
221 | Q_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 |
---|
227 | Q_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 | |
---|
247 | Q_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 |
---|
254 | Q_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 |
---|
263 | Q_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 |
---|
272 | Q_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 | |
---|
286 | Q_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 | |
---|
294 | Q_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 | |
---|
310 | Q_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 | |
---|
325 | Q_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 | |
---|
334 | Q_poly Q_poly::Q_poly_neg() |
---|
335 | { |
---|
336 | mpz_neg(denom,denom); |
---|
337 | } |
---|
338 | |
---|
339 | // Naive Multiplikation |
---|
340 | Q_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 | |
---|
382 | Q_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 !!! |
---|
388 | Q_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 | |
---|
438 | Q_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 | |
---|
448 | Q_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 |
---|
454 | Q_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 |
---|
473 | Q_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 | |
---|
503 | Q_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 | |
---|
551 | Q_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 |
---|
558 | Q_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 | |
---|
615 | Q_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 |
---|
624 | Q_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 |
---|
631 | Q_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 |
---|
641 | Q_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 |
---|
653 | void 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 | |
---|
668 | void 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 |
---|
687 | Q_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 |
---|
700 | Q_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 |
---|
708 | Q_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 |
---|
717 | Q_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 | |
---|
725 | int 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 |
---|
745 | int 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 |
---|
756 | int Q_poly::is_one() |
---|
757 | { |
---|
758 | if (deg==0) |
---|
759 | { |
---|
760 | if (mpz_cmp(coef[0],denom)==0) { return 1; } |
---|
761 | else { return 0; } |
---|
762 | } |
---|
763 | else { return 0; } |
---|
764 | } |
---|
765 | |
---|
766 | int 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 | |
---|
776 | Q_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 |
---|
807 | Q_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 | |
---|
878 | void 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 |
---|
898 | void 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 | |
---|