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