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