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