AEp.cc
Go to the documentation of this file.
1 #include "misc/auxiliary.h"
2 #include "omalloc/omalloc.h"
3 
4 #ifdef SINGULAR_4_2
5 #include "AEp.h"
6 
7 #include <stdio.h>
8 #include <math.h>
9 
10 
11 using namespace std;
12 
13 //Konstruktoren
14 
15 p_poly::p_poly()
16 {
17  //coef = reinterpret_cast<mpz_t*> (omAlloc(100*sizeof(mpz_t)));
18  deg=-1;
19  mod=2;
20  mpz_init_set_ui(coef[0],0);
21 }
22 
23 
24 
25 p_poly::p_poly(int n,int p, mpz_t *a)
26 {
27 
28  deg=n;
29  mod=p;
30 
31  for ( int i=0;i<=n;i++)
32  {
33  mpz_mod_ui(a[i],a[i],mod);
34  mpz_init_set(coef[i], a[i]);
35  }
36 
37 }
38 
39 /*
40 //Destruktor
41 
42 p_poly::~p_poly()
43 {
44  delete[] coef;
45 }
46 
47 */
48 
49 //Reduktion modulo p
50 
51 void p_poly::p_poly_reduce(p_poly f,int p)
52 {
53  if (f.is_zero()==0)
54  {
55 
56  for (int i=f.deg;i>=0;i--)
57  {
58  mpz_mod_ui(coef[i],f.coef[i],p);
59  }
60  }
61  //Hier nötige Grad Korrektur
62  int k=deg;
63  while(mpz_sgn(coef[k])==0 && k>=0)
64  {deg--;k--;}
65 }
66 
67 
68 // Arithmetik
69 
70 
71 //Additionen
72 
73 //Standard - Addition
74 void p_poly::p_poly_add(const p_poly a, const p_poly b)
75 {
76  if (a.deg >=b.deg)
77  {
78 
79  deg=a.deg;
80  mod=a.mod;
81 
82  for ( int i=0;i<=b.deg;i++)
83  {
84  mpz_add(coef[i],a.coef[i],b.coef[i]);
85  }
86 
87  for ( int i=b.deg+1;i<=a.deg;i++)
88  {
89  mpz_init_set(coef[i],a.coef[i]);
90  }
91 
92  //Hier nötige Grad Korrektur
93  int k=deg;
94  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
95  {deg--;k--;}
96 
97 
98  }
99 
100  else {p_poly_add(b,a); }
101 
102 }
103 
104 //Überschreibende Addition
105 
106 void p_poly::p_poly_add_to(const p_poly g)
107 {
108  this->p_poly_add(*this,g);
109 }
110 
111 //Addition einer Konstanten
112 void p_poly::p_poly_add_const(p_poly f,const mpz_t a)
113 {
114  if (f.is_zero()==1 && mpz_divisible_ui_p(a,f.mod)==0)
115  p_poly_set(a,f.mod);
116 
117  else if (mpz_divisible_ui_p(a,f.mod)==0)
118  {
119  p_poly_set(f);
120  mpz_add(coef[0],coef[0],a);
121  /*/Hier nötige Grad Korrektur
122  int k=deg;
123  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
124  {deg--;k--;}
125  */
126  }
127 
128 }
129 
130 
131 //To Variante Addition einer Konstanten
132 
133 void p_poly::p_poly_add_const_to(const mpz_t a)
134 {
135  this->p_poly_add_const(*this,a);
136 }
137 
138 //Monom Addition
139 void p_poly::p_poly_add_mon(const p_poly f, mpz_t a, int i)
140 {
141  p_poly_set(f);
142  if (i<=deg && is_zero()==0)
143  {
144  mpz_add(coef[i],coef[i],a);
145  }
146  else if (is_zero()==1 && mpz_divisible_ui_p(a,f.mod)==0)
147  {
148  deg=i;
149  for(int j=0;j<=i;j++)
150  {
151  mpz_init_set_ui(coef[j],0);
152  }
153  mpz_t temp;
154  mpz_init_set_ui(temp,0);
155  mpz_mod_ui(temp,a,mod);
156  mpz_add(coef[i],coef[i],temp);
157 
158  }
159  else if(i>deg && mpz_divisible_ui_p(a,f.mod)==0)
160  {
161  deg=i;
162  for(int j=i;j>deg;j--)
163  {
164  mpz_init_set_ui(coef[j],0);
165  }
166  mpz_t temp;
167  mpz_init_set_ui(temp,0);
168  mpz_mod_ui(temp,a,mod);
169  mpz_add(coef[i],coef[i],temp);
170 
171  }
172  /*/Hier nötige Grad Korrektur
173  int k=deg;
174  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
175  {deg--;k--;}
176  */
177 
178 }
179 
180 //To Variante Monomaddition
181 void p_poly::p_poly_add_mon_to(mpz_t a, int i)
182 {
183 
184  if (i<=deg && is_zero()==0)
185  {
186  mpz_add(coef[i],coef[i],a);
187  }
188  else if (is_zero()==1 && mpz_divisible_ui_p(a,mod)==0)
189  {
190  deg=i;
191  for(int j=0;j<=i;j++)
192  {
193  mpz_init_set_ui(coef[j],0);
194  }
195  mpz_t temp;
196  mpz_init_set_ui(temp,0);
197  mpz_mod_ui(temp,a,mod);
198  mpz_add(coef[i],coef[i],temp);
199 
200  }
201  else if(i>deg && mpz_divisible_ui_p(a,mod)==0)
202  {
203  deg=i;
204  for(int j=i;j>deg;j--)
205  {
206  mpz_init_set_ui(coef[j],0);
207  }
208  mpz_t temp;
209  mpz_init_set_ui(temp,0);
210  mpz_mod_ui(temp,a,mod);
211  mpz_add(coef[i],coef[i],temp);
212 
213  }
214  /*Hier nötige Grad Korrektur
215  int k=deg;
216  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
217  {deg--;k--;} */
218 }
219 
220 //Subtraktionen
221 
222 void p_poly::p_poly_sub(const p_poly a, const p_poly b)
223 {
224  if (a.deg >=b.deg)
225  {
226 
227  deg=a.deg;
228  mod=a.mod;
229 
230  for ( int i=0;i<=b.deg;i++)
231  {
232  mpz_sub(coef[i],a.coef[i],b.coef[i]);
233  }
234 
235  for ( int i=b.deg+1;i<=a.deg;i++)
236  {
237  mpz_init_set(coef[i],a.coef[i]);
238  }
239 
240  //Hier nötige Grad Korrektur
241  int k=deg;
242  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
243  {deg--;k--;}
244 
245  }
246 
247  else {p_poly_sub(b,a);p_poly_neg(); }
248 
249 }
250 
251 
252 //Überschreibende Subtraktion
253 
254 void p_poly::p_poly_sub_to(const p_poly b)
255 {
256  this->p_poly_sub(*this,b);
257 }
258 
259 //Subtraktion einer Konstanten
260 void p_poly::p_poly_sub_const(p_poly f,const mpz_t a)
261 {
262  if (f.is_zero()==1)
263  {
264  p_poly_set(a,f.mod);
265  p_poly_neg();
266  }
267  else
268  {
269  p_poly_set(f);
270  mpz_sub(coef[0],coef[0],a);
271  }
272  /*/*Hier nötige Grad Korrektur
273  int k=deg;
274  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
275  {deg--;k--;} */
276 
277 }
278 
279 
280 //To Variante Subtraktion einer Konstanten
281 
282 void p_poly::p_poly_sub_const_to(const mpz_t a)
283 {
284  this->p_poly_sub_const(*this,a);
285 }
286 
287 
288 //Monom Subtraktion
289 void p_poly::p_poly_sub_mon(const p_poly f , mpz_t a, int i)
290 {
291  mpz_t temp;
292  mpz_neg(temp,a);
293  p_poly_add_mon(f,temp,i);
294 }
295 
296 //To Variante Monomaddition
297 void p_poly::p_poly_sub_mon_to(mpz_t a, int i)
298 {
299  mpz_t temp;
300  mpz_neg(temp,a);
301  p_poly_add_mon_to(temp,i);
302 }
303 
304 
305 //Multiplikationen
306 
307 //Multiplikation mit Monom
308 void p_poly::p_poly_mon_mult( p_poly f, int n)
309 {
310  if (f.is_zero()==1)
311  {p_poly_set_zero();}
312  else
313  {
314  deg=f.deg+n;
315  mod=f.mod;
316  for (int i=deg;i>=n;i--)
317  {
318  mpz_init_set(coef[i],f.coef[i-n]);
319  }
320  for (int i=n-1;i>=0;i--)
321  {
322  mpz_init_set_ui(coef[i],0);
323  }
324 
325  /*/Hier nötige Grad Korrektur
326  int k=deg;
327  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
328  {deg--;k--;} */
329  }
330 }
331 
332 void p_poly::p_poly_mon_mult_to(const int n)
333 {
334  this->p_poly_mon_mult(*this,n);
335 }
336 
337 
338 //Multiplikation mit Skalar
339 
340 void p_poly::p_poly_scalar_mult(const p_poly g, const mpz_t n)
341 {
342  if (mpz_divisible_ui_p(n,g.mod)!=0)
343  p_poly_set_zero();
344  else
345  {
346  deg=g.deg;
347  mod=g.mod;
348 
349  mpz_t temp;
350  mpz_init_set_ui(temp,0);
351  for(int i=0;i<=deg;i++)
352  {
353  mpz_mul(temp,n,g.coef[i]);
354  mpz_init_set(coef[i],temp);
355  }
356  }
357 }
358 
359 
360 
361 void p_poly::p_poly_scalar_mult(const mpz_t n, const p_poly g)
362 {
363  if (mpz_divisible_ui_p(n,g.mod)!=0)
364  p_poly_set_zero();
365  else
366  {
367  deg=g.deg;
368  mod=g.mod;
369 
370 
371  mpz_t temp;
372  mpz_init_set_ui(temp,0);
373  for(int i=0;i<=deg;i++)
374  {
375  mpz_mul(temp,n,g.coef[i]);
376  mpz_init_set(coef[i],temp);
377  }
378  /*/Hier nötige Grad Korrektur
379  int k=deg;
380  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
381  {deg--;k--;} */
382  }
383 }
384 
385 
386 void p_poly::p_poly_scalar_mult_to(const mpz_t n)
387 {
388  this->p_poly_scalar_mult(*this,n);
389 }
390 
391 
392 
393 // Negation
394 
395 void p_poly::p_poly_neg()
396 {
397  for (int i=0;i<=deg;i++)
398  {
399  mpz_neg(coef[i],coef[i]);
400  }
401 }
402 
403 // Naive Multiplikation
404 void p_poly::p_poly_mult_n(p_poly a,p_poly b)
405 {
406  //Reduktion mod p
407  a.p_poly_reduce(a,a.mod);
408  b.p_poly_reduce(b,b.mod);
409 
410  if (a.is_zero()==1 || b.is_zero()==1)
411  p_poly_set_zero();
412  else
413  {
414  mpz_t temp;
415  mpz_init_set_ui(temp,0);
416  deg = a.deg + b.deg;
417 
418  // Kopien atemp und btemp von a bzw. b, mit Nullen aufgefüllt
419  p_poly atemp, btemp;
420  atemp.p_poly_set(a);
421  btemp.p_poly_set(b);
422  for(int i=a.deg+1;i<=deg;i++)
423  {
424  mpz_init_set_ui(atemp.coef[i],0);
425  }
426  for(int i=b.deg+1;i<=deg;i++)
427  {
428  mpz_init_set_ui(btemp.coef[i],0);
429  }
430  atemp.deg = deg;
431  btemp.deg = deg;
432 
433  // Multiplikationsalgorithmus
434  for (int k=0; k<=deg; k++)
435  {
436  mpz_init_set_ui(coef[k],0); // k-ter Koeffizient zunächst 0
437  for (int i=0; i<=k; i++) // dann schrittweise Summe der a[i]*b[k-i]/
438  {
439  mpz_mul(temp,atemp.coef[i],btemp.coef[k-i]);
440  mpz_add(coef[k],coef[k],temp);
441  }
442  }
443 
444  /*/Hier nötige Grad Korrektur
445  int k=deg;
446  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
447  {deg--;k--;} */
448 
449 
450  }
451 }
452 
453 
454 //Überschreibende Multiplikation
455 
456 void p_poly::p_poly_mult_n_to(const p_poly g)
457 {
458  this->p_poly_mult_n(*this,g);
459 
460 }
461 
462 // Karatsuba-Multiplikation (Weimerskirch/Paar Alg. 1), ACHTUNG VORLÄUFIGE VERSION, macht noch Fehler beim Grad und ist unelegant !!!
463 void p_poly::p_poly_mult_ka( p_poly A, p_poly B)
464 {
465 
466  if (A.is_zero()==1 || B.is_zero()==1)
467  {
468  p_poly_set_zero();
469  }
470  else
471  {
472  //Reduktion mod p
473  A.p_poly_reduce(A,A.mod);
474  B.p_poly_reduce(B,B.mod);
475 
476  // Größeren Grad feststellen
477  int n;
478  if(A.deg>=B.deg){n=A.deg+1;}
479  else{n=B.deg+1;}
480  // n auf nächste 2er-Potenz setzen (VORLÄUFIG!)
481  n = static_cast<int>(ceil(log(n)/log(2)));
482  n = static_cast<int>(pow(2,n));
483 
484  if (n==1)
485  {
486  mpz_t AB;
487  mpz_mul(AB,A.coef[0],B.coef[0]);
488  p_poly_set(AB,A.mod);
489  this->p_poly_reduce(*this,A.mod);
490  }
491  else
492  {
493  // p_polynome A und B aufspalten
494  p_poly Au, Al, Bu, Bl;
495  Au.p_poly_mon_div(A,n/2);
496  Al.p_poly_mon_div_rem(A,n/2);
497  Bu.p_poly_mon_div(B,n/2);
498  Bl.p_poly_mon_div_rem(B,n/2);
499  p_poly Alu,Blu;
500  Alu.p_poly_add(Al,Au);
501  Blu.p_poly_add(Bl,Bu);
502 
503  // Teile rekursiv multiplizieren
504  p_poly D0, D1, D01;
505  D0.p_poly_mult_ka(Al,Bl);
506  D1.p_poly_mult_ka(Au,Bu);
507  D01.p_poly_mult_ka(Alu,Blu);
508 
509  // Ergebnis zusammensetzen
510  p_poly temp;
511  D01.p_poly_sub_to(D0);
512  D01.p_poly_sub_to(D1);
513  D01.p_poly_mon_mult_to(n/2);
514  D1.p_poly_mon_mult_to(n);
515  D1.p_poly_add_to(D01);
516  D1.p_poly_add_to(D0);
517  p_poly_set(D1);
518 
519  }
520 
521  //Hier nötige Grad Korrektur
522  int k=deg;
523  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
524  {deg--;k--;}
525  }
526 }
527 
528 
529 
530 //Skalare Divisionen
531 
532 void p_poly::p_poly_scalar_div( const p_poly g, const mpz_t n)
533 {
534 
535  if (mpz_divisible_ui_p(n,g.mod)==0) // überprüft invertierbarkeit
536  {
537  deg=g.deg;
538  mod=g.mod;
539 
540 
541  mpz_t temp;
542  mpz_t p;
543  mpz_init_set_ui(temp,0);
544  mpz_init_set_ui(p,mod);
545  for(int i=0;i<=deg;i++)
546  {
547  mpz_invert(temp,n,p);
548  mpz_mul(temp,g.coef[i],temp);
549  mpz_init_set(coef[i],temp);
550  }
551 
552  /*/Hier nötige Grad Korrektur
553  int k=deg;
554  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
555  {deg--;k--;} */
556  }
557 }
558 
559 
560 void p_poly::p_poly_scalar_div_to(const mpz_t n)
561 {
562  this->p_poly_scalar_div(*this,n);
563 }
564 
565 // Division durch Monom - Quotient
566 void p_poly::p_poly_mon_div(const p_poly f, const int n)
567 {
568  if (f.deg<n)
569  {
570  p_poly_set_zero();
571  }
572  else
573  {
574  deg=f.deg-n;
575  mod=f.mod;
576 
577  for (int i=0;i<=f.deg-n;i++)
578  {
579  mpz_init_set(coef[i],f.coef[n+i]);
580  }
581  }
582 }
583 
584 // Division durch Monom - Rest
585 void p_poly::p_poly_mon_div_rem(const p_poly f, const int n)
586 {
587  if (f.deg<n)
588  {
589  p_poly_set(f);
590  }
591  else
592  {
593  deg=n-1;
594  mod=f.mod;
595 
596 
597  for (int i=0;i<=n-1;i++)
598  {
599  mpz_init_set(coef[i],f.coef[i]);
600  }
601  }
602 }
603 
604 
605 
606 
607 //Euklidische Division nach Cohen Algo 3.1.1 (degA muss größer gleich deg B sein)!!
608 
609 void p_poly::p_poly_div_rem( p_poly A, p_poly B)
610 {
611 
612  if (B.is_zero()==0)
613  {
614  //Reduktion mod p
615  A.p_poly_reduce(A,A.mod);
616  B.p_poly_reduce(B,B.mod);
617 
618  p_poly R;
619  p_poly temp;
620  R.p_poly_set(A);
621  mpz_t a;
622  mpz_t u;
623  mpz_t p;
624  mpz_init_set_ui(p,A.mod);
625  mpz_init_set_ui(a,0);
626  mpz_init_set_ui(u,0);
627  int i;
628 
629  mpz_invert(u,B.coef[B.deg],p); // Inverse von lc(B)
630 
631 
632  while (R.deg>=B.deg)
633  {
634 
635  mpz_mul(a,R.coef[R.deg],u);
636  i=R.deg-B.deg;
637 
638  temp.p_poly_mon_mult(B,i);
639  temp.p_poly_scalar_mult_to(a);
640 
641  R.p_poly_sub_to(temp);
642 
643  }
644  p_poly_set(R);
645 
646  /*/Hier nötige Grad Korrektur
647  int k=deg;
648  while(mpz_divisible_ui_p(coef[k],mod)!=0 && k>=0)
649  {deg--;k--;}*/
650  }
651 }
652 //To Variante von Algo 3.1.1 im Cohen
653 
654 void p_poly::p_poly_div_rem_to(const p_poly B)
655 {
656  this->p_poly_div_rem(*this,B);
657 
658 
659 }
660 
661 
662 
663 //Exakte Division nach Cohen 3.1.1
664 void p_poly::p_poly_div(p_poly &Q, p_poly &R, p_poly A, p_poly B)
665 {
666  if (B.is_zero()==0)
667  {
668  //Reduktion mod p
669  A.p_poly_reduce(A,A.mod);
670  B.p_poly_reduce(B,B.mod);
671 
672  //Initialisierungen
673  p_poly temp;
674  R.p_poly_set(A);
675  Q.p_poly_set_zero();
676  Q.mod=A.mod;
677 
678  mpz_t a;
679  mpz_t b;
680  mpz_t p;
681  mpz_init_set_ui(p,A.mod);
682  mpz_init_set_ui(a,0);
683  mpz_init_set_ui(b,0);
684  int i;
685  mpz_invert(b,B.coef[B.deg],p);
686 
687  //Algorithmus TO DO: evtl hier mit auch den Rest ausgeben
688  while (R.deg>=B.deg)
689  {
690 
691  mpz_mul(a,R.coef[R.deg],b);
692  i=R.deg-B.deg;
693 
694  temp.p_poly_mon_mult(B,i);
695  temp.p_poly_scalar_mult_to(a);
696 
697  R.p_poly_sub_to(temp);
698 
699  Q.p_poly_add_mon_to(a,i);
700 
701  R.p_poly_reduce(R,R.mod);
702  Q.p_poly_reduce(Q,Q.mod);
703  }
704 
705  /*/Hier nötige Grad Korrektur Q
706  int k=Q.deg;
707  while(mpz_divisible_ui_p(Q.coef[k],Q.mod)!=0 && k>=0)
708  {Q.deg--;k--;}*/
709 
710  /*/Hier nötige Grad Korrektur R
711  k=R.deg;
712  while(mpz_divisible_ui_p(R.coef[k],R.mod)!=0 && k>=0)
713  {R.deg--;k--;} */
714  }
715 }
716 
717 
718 //To Varainte der exakten Division
719 
720 void p_poly::p_poly_div_to(p_poly &Q,p_poly &R, p_poly B)
721 {
722  this->p_poly_div(Q ,R,*this,B);
723 }
724 
725 
726 // Kombinationen
727 
728 // a := a*b + c
729 void p_poly::p_poly_multadd_to(const p_poly b, const p_poly c)
730 {
731  p_poly_mult_n_to(b);
732  p_poly_add_to(c);
733 }
734 
735 //a=a*b-c
736 void p_poly::p_poly_multsub_to(const p_poly b, const p_poly c)
737 {
738  p_poly_mult_n_to(b);
739  p_poly_sub_to(c);
740 }
741 
742 
743 
744 /*
745 // a := (a+b)* c
746 void p_poly::poly_addmult_to(const p_poly b, const p_poly c)
747 {
748  p_poly a(deg,coef);
749  a.poly_add_to(b);
750  a.poly_mult_n_to(c);
751  poly_set(a);
752 }
753 */
754 
755 
756 
757 //Sonstiges
758 void p_poly::p_poly_horner(mpz_t erg, const mpz_t u)
759 {
760  if (is_zero()==0)
761  {
762  mpz_init_set(erg,coef[deg]);
763  for (int i=deg;i>=1;i--)
764  {
765  mpz_mul(erg,erg,u);
766  mpz_add(erg,erg,coef[i-1]);
767  }
768 
769  //Reduktion
770  mpz_mod_ui(erg,erg,mod);
771  }
772  else
773  {
774  mpz_set_ui(erg,0);
775  }
776 }
777 
778 // p_polynom in p_polynom einsetzen(Horner-Schema) KRITISCHE EINGABE x^2, x^2 ....
779 
780 void p_poly::p_poly_horner_p_poly( p_poly A, p_poly B)
781 {
782  //Reduktion mod p
783  A.p_poly_reduce(A,A.mod);
784  B.p_poly_reduce(B,B.mod);
785 
786  p_poly_set(A.coef[A.deg],A.mod);
787  for (int i=A.deg;i>=1;i--)
788  {
789  p_poly_mult_n_to(B);
790  p_poly_add_const_to(A.coef[i-1]);
791  }
792  /*/Hier nötige Grad Korrektur
793  int i=deg;
794  while(mpz_divisible_ui_p(coef[i],mod)!=0 && i>=0)
795  {deg--;i--;} */
796 }
797 
798 
799 
800 //Hilfsfunktionen
801 
802 
803 //setze p_polynom auf p_polynom b
804 void p_poly::p_poly_set(const p_poly b)
805 {
806  deg=b.deg;
807  mod=b.mod;
808 
809 
810  for(int i=0;i<=deg;i++)
811  {
812  mpz_init_set(coef[i],b.coef[i]); // Hier wird init set dringendst benötigt
813  }
814 
815 }
816 
817 // setze p_polynom auf konstantes p_polynom b
818 void p_poly::p_poly_set(const mpz_t b,int p)
819 {
820  deg=0;
821  mod=p;
822 
823  if (mpz_divisible_ui_p(b,mod)!=0)
824  p_poly_set_zero();
825  else
826  {
827  mpz_t temp;
828  mpz_init_set(temp,b);
829  mpz_mod_ui(temp,temp,p);
830  mpz_init_set(coef[0],b);
831  }
832 }
833 
834 
835 //setze p_polynom auf Nullpolynom
836 void p_poly::p_poly_set_zero()
837 {
838  deg = -1;
839 }
840 
841 
842 //Vergleiche ob 2 p_polynome gleich return 1 falls ja sont 0
843 
844 int p_poly::is_equal(const p_poly g) const
845 {
846  if (deg!=g.deg)
847  return 0;
848  else
849 
850  for (int i=deg;i>=0; i--)
851  {
852  if (mpz_cmp(coef[i],g.coef[i])!=0)
853  {return 0;}
854  }
855  return 1;
856 }
857 
858 //Überprüft ob das p_polynom 0 ist
859 
860 int p_poly::is_zero() const
861 {
862  if (deg<0)
863  return 1;
864  else
865  return 0;
866 
867 }
868 
869 int p_poly::is_one() const
870 {
871  if (deg==0)
872  {
873  if (mpz_cmp_ui(coef[0],1)==0) { return 1; }
874  else { return 0; }
875  }
876  else { return 0; }
877 }
878 
879 int p_poly::is_monic() const
880 {
881  if (mpz_cmpabs_ui(coef[deg],1)==0)
882  return 1;
883  else
884  return 0;
885 }
886 
887 // klassischer GGT nach Cohen 3.2.1
888 
889 void p_poly::p_poly_gcd( p_poly A, p_poly B)
890 {
891 
892  //Reduktion mod p
893  A.p_poly_reduce(A,A.mod);
894  B.p_poly_reduce(B,B.mod);
895 
896  if (A.deg<B.deg)
897  p_poly_gcd(B,A);
898  else if (B.is_zero()==1)
899  p_poly_set(A);
900  else
901  {
902  p_poly App;
903  p_poly Bpp;
904  p_poly R;
905  App.p_poly_set(A);
906  Bpp.p_poly_set(B);
907 
908  while (Bpp.is_zero()==0)
909  {
910  R.p_poly_div_rem(App,Bpp);
911  App.p_poly_set(Bpp);
912  Bpp.p_poly_set(R);
913  }
914  p_poly_set(App);
915  }
916 
917 }
918 
919 //Nach nach Fieker 2.12 Symbolisches Rechnen (2012)
920 // gibt g=s*A+t*B aus
921 void p_poly::p_poly_extgcd(p_poly &s,p_poly &t,p_poly &g, p_poly A, p_poly B)
922 {
923 
924  //Reduktion mod p
925  A.p_poly_reduce(A,A.mod);
926  B.p_poly_reduce(B,B.mod);
927 
928 
929  if (A.deg<B.deg)
930  p_poly_extgcd(t,s,g,B,A);
931  else if (B.is_zero()==1)
932  {
933  g.p_poly_set(A);
934  t.p_poly_set_zero();
935 
936  mpz_t temp;
937  mpz_init_set_ui(temp,1);
938  s.p_poly_set(temp,A.mod);
939  }
940 
941  else
942  {
943  mpz_t temp;
944  mpz_init_set_ui(temp,1);
945 
946  p_poly R1;
947  R1.p_poly_set(A);
948  p_poly R2;
949  R2.p_poly_set(B);
950  p_poly R3;
951  R3.mod=A.mod;
952 
953 
954  p_poly S1;
955  S1.p_poly_set(temp,A.mod);
956  p_poly S2;
957  S2.p_poly_set_zero();
958  S2.mod=A.mod;
959  p_poly S3;
960  S3.mod=A.mod;
961 
962  p_poly T1;
963  T1.p_poly_set_zero();
964  T1.mod=A.mod;
965  p_poly T2;
966  T2.p_poly_set(temp,A.mod);
967  p_poly T3;
968  T3.mod=A.mod;
969 
970  p_poly Q;
971 
972  while (R2.is_zero()!=1)
973  {
974  p_poly_div(Q,R3,R1,R2);
975 
976  S3.p_poly_mult_n(Q,S2);
977  S3.p_poly_neg();
978  S3.p_poly_add_to(S1);
979 
980  T3.p_poly_mult_n(Q,T2);
981  T3.p_poly_neg();
982  T3.p_poly_add_to(T1);
983 
984  R1.p_poly_set(R2);
985  R2.p_poly_set(R3);
986 
987  S1.p_poly_set(S2);
988  S2.p_poly_set(S3);
989 
990  T1.p_poly_set(T2);
991  T2.p_poly_set(T3);
992  }
993  t.p_poly_set(T1);
994  s.p_poly_set(S1);
995  g.p_poly_set(R1);
996  }
997 }
998 
999 
1000 //Ein & Ausgabe
1001 
1002 //Eingabe
1003 
1004 void p_poly::p_poly_insert()
1005 {
1006 #if 0
1007  cout << "Bitte geben Sie ein p_polynom ein! Zunächst den Grad: " << endl;
1008  cin >> deg;
1009  cout << "Jetzt den modul: " << endl;
1010  cin >> mod;
1011 
1012  for ( int i=0; i<=deg;i++)
1013  {
1014  mpz_init_set_ui(coef[i],0);
1015  printf("Geben Sie nun f[%i] ein:",i);
1016  mpz_inp_str(coef[i],stdin, 10);
1017  }
1018  //Reduktion
1019  this->p_poly_reduce(*this,mod);
1020 #endif
1021 }
1022 
1023 
1024 //Ausgabe
1025 void p_poly::p_poly_print()
1026 {
1027 #if 0
1028 
1029  //Reduktion
1030  this->p_poly_reduce(*this,mod);
1031 
1032 
1033  if (is_zero()==1)
1034  cout << "0" << "\n" <<endl;
1035  else
1036  {
1037  for (int i=deg;i>=1;i--)
1038  {
1039  mpz_out_str(stdout,10, coef[i]);
1040  printf("X%i+",i);
1041  }
1042  mpz_out_str(stdout,10, coef[0]);
1043  printf("\n");
1044  }
1045 #endif
1046 }
1047 
1048 #endif
void T1(ideal h)
Definition: cohomo.cc:2473
const CanonicalForm int s
Definition: facAbsFact.cc:55
int j
Definition: facHensel.cc:105
CF_NO_INLINE CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Definition: cf_inline.cc:564
gmp_float log(const gmp_float &a)
Definition: mpr_complex.cc:344
STL namespace.
g
Definition: cfModGcd.cc:4031
int k
Definition: cfEzgcd.cc:92
#define Q
Definition: sirandom.c:25
CanonicalForm b
Definition: cfModGcd.cc:4044
const signed long ceil(const ampf< Precision > &x)
Definition: amp.h:789
#define A
Definition: sirandom.c:23
All the auxiliary stuff.
void T2(ideal h)
Definition: cohomo.cc:2754
FILE * f
Definition: checklibs.c:9
int i
Definition: cfEzgcd.cc:125
CanonicalForm & mod(const CanonicalForm &)
b *CanonicalForm B
Definition: facBivar.cc:52
#define R
Definition: sirandom.c:26
int p
Definition: cfModGcd.cc:4019
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:414