00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_POLYNOMIAL_HPP
00014 #define __MMX_POLYNOMIAL_HPP
00015 #include <algebramix/polynomial_naive.hpp>
00016
00017 namespace mmx {
00018 #define Polynomial_variant(C) polynomial_variant_helper<C>::PV
00019 #define Polynomial_multiply(C) polynomial_variant_helper<C>::PV
00020 #define TMPL_DEF \
00021 template<typename C, typename V=typename Polynomial_variant(C) >
00022 #define TMPL template<typename C, typename V>
00023 #define TMPLK template<typename C, typename V, typename K>
00024 #define Format format<C>
00025 #define Vector vector<C>
00026 #define Polynomial polynomial<C,V>
00027 #define Polynomial_rep polynomial_rep<C,V>
00028 #define Abs_polynomial polynomial<Abs_type(C),V>
00029 #define Real_polynomial polynomial<Real_type(C),V>
00030 #define Center_polynomial polynomial<Center_type(C),V>
00031 #define Radius_polynomial polynomial<Radius_type(C),V>
00032 #define Lifted_polynomial polynomial<Lift_type(C)>
00033 #define Projected_polynomial polynomial<Project_type(C)>
00034 #define Reconstructed_polynomial polynomial<Reconstruct_type(C)>
00035 TMPL class polynomial_rep;
00036 TMPL class polynomial;
00037 TMPL inline int deg (const Polynomial& P);
00038 TMPL inline int degree (const Polynomial& P);
00039 TMPL inline nat N (const Polynomial& P);
00040 TMPL inline C* seg (Polynomial& P);
00041 TMPL inline const C* seg (const Polynomial& P);
00042
00043
00044
00045
00046
00047 TMPL_DEF
00048 class polynomial_rep REP_STRUCT_1(C) {
00049 public:
00050 typedef implementation<vector_linear,V> Vec;
00051 typedef implementation<polynomial_linear,V> Pol;
00052 private:
00053 C* a;
00054 nat n;
00055 nat l;
00056 public:
00057 inline void normalize () {
00058 while ((n != 0) && is_exact_zero (a[n-1])) n--; }
00059 inline polynomial_rep (C* a2, nat n2, const Format& fm):
00060 Format (fm), a(a2), n(n2), l(n2) { normalize (); }
00061 inline polynomial_rep (C* a2, nat n2, nat l2, const Format& fm):
00062 Format (fm), a(a2), n(n2), l(l2) { normalize (); }
00063 inline ~polynomial_rep () {
00064 mmx_delete<C> (a, l); }
00065 void extend (nat n2) {
00066 nat l2= aligned_size<C,V> (n2);
00067 C* b= mmx_formatted_new<C> (l2, this->tfm ());
00068 Pol::copy (b, a, n);
00069 Pol::clear (b+n, l2-n);
00070 mmx_delete<C> (a, l);
00071 a= b;
00072 n= n2;
00073 l= l2;
00074 }
00075 void resize (nat n2) {
00076 if (n < n2) extend (n2);
00077 n= n2;
00078 }
00079
00080 friend class Polynomial;
00081 friend int deg LESSGTR (const Polynomial& P);
00082 friend int degree LESSGTR (const Polynomial& P);
00083 friend nat N LESSGTR (const Polynomial& P);
00084 friend C* seg LESSGTR (Polynomial& P);
00085 friend const C* seg LESSGTR (const Polynomial& P);
00086 };
00087
00088 TMPL_DEF
00089 class polynomial {
00090 INDIRECT_PROTO_2 (polynomial, polynomial_rep, C, V)
00091 public:
00092 typedef implementation<vector_linear,V> Vec;
00093 typedef implementation<polynomial_linear,V> Pol;
00094 typedef typename Pol::template global_variables<Polynomial> S;
00095 public:
00096 static inline generic get_variable_name () {
00097 return S::get_variable_name (); }
00098 static inline void set_variable_name (const generic& x) {
00099 S::set_variable_name (x); }
00100
00101
00102 inline polynomial () {
00103 nat l= aligned_size<C,V> (0);
00104 rep= new Polynomial_rep (mmx_new<C> (l), 0, l, Format (no_format ()));
00105 }
00106
00107 inline polynomial (const Format& fm) {
00108 nat l= aligned_size<C,V> (0);
00109 rep= new Polynomial_rep (mmx_new<C> (l), 0, l, fm);
00110 }
00111 template<typename T> inline polynomial (const T& x) {
00112 nat l= aligned_size<C,V> (1);
00113 C c= as<C> (x);
00114 C* a= mmx_formatted_new<C> (l, get_format (c));
00115 a[0]= c;
00116 rep= new Polynomial_rep (a, 1, l, get_format (c));
00117 rep->normalize ();
00118 }
00119
00120 template<typename T> inline
00121 polynomial (const polynomial<T,V>& P) {
00122 Format fm= as<Format > (CF(P));
00123 nat l= aligned_size<C,V> (N(P));
00124 C* a= mmx_formatted_new<C> (l, fm);
00125 Pol::cast (a, seg (P), N(P));
00126 rep= new Polynomial_rep (a, N(P), l, fm);
00127 rep->normalize ();
00128 }
00129
00130 template<typename T> inline
00131 polynomial (const polynomial<T,V>& P, const Format& fm) {
00132 nat l= aligned_size<C,V> (N(P));
00133 C* a= mmx_formatted_new<C> (l, fm);
00134 Pol::cast (a, seg (P), N(P));
00135 rep= new Polynomial_rep (a, N(P), l, fm);
00136 rep->normalize ();
00137 }
00138 polynomial (const C& x, nat n) {
00139 nat l= aligned_size<C,V> (n+1);
00140 C* a= mmx_formatted_new<C> (l, get_format (x));
00141 Pol::clear (a, n);
00142 a[n]= x;
00143 rep= new Polynomial_rep (a, n+1, l, get_format (x));
00144 rep->normalize ();
00145 }
00146 inline polynomial (C* a, nat n, const Format& fm) {
00147 rep= new Polynomial_rep (a, n, fm);
00148 }
00149 inline polynomial (C* a, nat n, nat l, const Format& fm) {
00150 rep= new Polynomial_rep (a, n, l, fm);
00151 }
00152 polynomial (const vector<C>& v) {
00153 nat l= aligned_size<C,V> (N(v));
00154 C* a= mmx_formatted_new<C> (l, CF(v));
00155 Pol::copy (a, seg (v), N(v));
00156 rep= new Polynomial_rep (a, N(v), l, CF(v));
00157 rep->normalize ();
00158 }
00159 polynomial (const iterator<C>& it) {
00160 nat i, l=0;
00161 nat s= aligned_size<C,V> (0);
00162 C* a= mmx_formatted_new<C> (s, CF(it));
00163 rep= new Polynomial_rep (a, 0, s, CF(it));
00164 for (i=0; busy (it); i++, ++it) {
00165 if (i >= l) {
00166 l= max ((nat) 1, l<<1);
00167 rep->extend (l);
00168 }
00169 rep->a[i]= *it;
00170 }
00171 rep->normalize ();
00172 }
00173 inline friend Polynomial copy (const Polynomial& P) {
00174 nat n= N(P), l= aligned_size<C,V> (n);
00175 C* r= mmx_formatted_new<C> (l, CF(P));
00176 Pol::copy (r, seg (P), n);
00177 return Polynomial (r, n, l, CF(P));
00178 }
00179 inline C operator[] (nat i) const {
00180 return i<rep->n? rep->a[i]: promote (0, CF(*this)); }
00181 Polynomial& operator += (const Polynomial& P);
00182 Polynomial& operator -= (const Polynomial& P);
00183 Polynomial& operator *= (const C& c);
00184 Polynomial& operator /= (const C& c);
00185 };
00186
00187 INDIRECT_IMPL_2 (polynomial, polynomial_rep, typename C, C, typename V, V)
00188 DEFINE_UNARY_FORMAT_2 (polynomial)
00189
00190 TMPL inline nat N (const Polynomial& P) { return P->n; }
00191 TMPL inline int deg (const Polynomial& P) { return ((int) P->n) - 1; }
00192 TMPL inline int degree (const Polynomial& P) { return ((int) P->n) - 1; }
00193 TMPL inline C* seg (Polynomial& P) { return P->a; }
00194 TMPL inline const C* seg (const Polynomial& P) { return P->a; }
00195 TMPL inline Format CF (const Polynomial& P) { return P->tfm (); }
00196
00197 TMPL inline vector<C>
00198 coefficients (const Polynomial& P) {
00199 vector<C> v (get_sample (CF(P)), N(P));
00200 for (nat i= 0; i < N(P); i++) v[i]= P[i];
00201 return v;
00202 }
00203
00204 TMPL inline vector<C>
00205 coefficients (const Polynomial& P, nat b, nat e) {
00206 if (b >= e) return vector<C> ();
00207 vector<C> v (get_sample (CF(P)), e - b);
00208 for (nat i= 0; i < e - b; i++) v[i]= P[b+i];
00209 return v;
00210 }
00211
00212 STYPE_TO_TYPE(TMPL,scalar_type,Polynomial,C);
00213 STYPE_TO_TYPE(TMPL,monomial_type,Polynomial,nat);
00214 UNARY_RETURN_TYPE(TMPL,Re_op,Polynomial,Real_polynomial);
00215 UNARY_RETURN_TYPE(TMPL,abs_op,Polynomial,Abs_polynomial);
00216 UNARY_RETURN_TYPE(TMPL,center_op,Polynomial,Center_polynomial);
00217 UNARY_RETURN_TYPE(TMPL,radius_op,Polynomial,Radius_polynomial);
00218 UNARY_RETURN_TYPE(TMPL,lift_op,Polynomial,Lifted_polynomial);
00219 UNARY_RETURN_TYPE(TMPL,project_op,Polynomial,Projected_polynomial);
00220 UNARY_RETURN_TYPE(TMPL,reconstruct_op,Polynomial,Reconstructed_polynomial);
00221
00222 TMPL
00223 struct fast_helper<Polynomial > {
00224 typedef Fast_type(C) FC;
00225 typedef polynomial<FC> fast_type;
00226 typedef typename Polynomial_variant(FC) FV;
00227 static inline fast_type dd (const Polynomial& p) {
00228 format<FC> fm= fast (CF(p));
00229 nat l= aligned_size<FC,FV> (N(p));
00230 FC* r= mmx_formatted_new<FC> (l, fm);
00231 for (nat i=0; i<N(p); i++) r[i]= fast (p[i]);
00232 return fast_type (r, N(p), l, fm); }
00233 static inline Polynomial uu (const fast_type& p) {
00234 Format fm= slow<C> (CF(p));
00235 nat l= aligned_size<C,V> (N(p));
00236 C* r= mmx_formatted_new<C> (l, fm);
00237 for (nat i=0; i<N(p); i++) r[i]= slow<C> (p[i]);
00238 return Polynomial (r, N(p), l, fm); }
00239 };
00240
00241 template<typename T, typename F, typename TV, typename FV>
00242 struct as_helper<polynomial<T,TV>,polynomial<F,FV> > {
00243 static inline polynomial<T,TV>
00244 cv (const polynomial<F,FV>& p) {
00245 format<T> fm= as<format<T> > (CF(p));
00246 nat l= aligned_size<T,TV> (N(p));
00247 T* c= mmx_formatted_new<T> (l, fm);
00248 for (nat i=0; i<N(p); i++) c[i]= as<T> (p[i]);
00249
00250
00251 polynomial<T,TV> r= polynomial<T,TV> (c, N(p), l, fm);
00252 return r;
00253 }
00254 };
00255
00256 template<typename T, typename F, typename TV, typename FV> inline void
00257 set_as (polynomial<T,TV>& r, const polynomial<F,FV>& p) {
00258 r.secure ();
00259 inside (r) -> resize (N(p));
00260 T* c= seg (r);
00261 for (nat i=0; i<N(p); i++) set_as (c[i], p[i]);
00262 inside (r) -> normalize ();
00263 }
00264
00265 template<typename C, typename V, typename T> inline void
00266 set_as (Polynomial& r, const T& x) {
00267 r.secure ();
00268 inside (r) -> resize (1);
00269 set_as (seg(r)[0], x);
00270 inside (r) -> normalize ();
00271 }
00272
00273 template<typename Op, typename C, typename V> nat
00274 unary_hash (const Polynomial& p) {
00275 register nat i, h= 642531, n= N(p);
00276 for (i=0; i<n; i++)
00277 h= (h<<1) ^ (h<<5) ^ (h>>27) ^ Op::op (p[i]);
00278 return h;
00279 }
00280
00281
00282
00283
00284
00285 TMPL_DEF
00286 class polynomial_iterator_rep: public iterator_rep<C> {
00287 Polynomial P;
00288 nat i;
00289 protected:
00290 bool is_busy () { return ((int) i) <= deg (P); }
00291 void advance () { i++; }
00292 C current () { return P[i]; }
00293 public:
00294 polynomial_iterator_rep (const Polynomial& P2):
00295 iterator_rep<C> (CF(P2)), P(P2), i(0) {}
00296 };
00297
00298 TMPL iterator<C>
00299 iterate (const Polynomial& P) {
00300 return iterator<C> (new polynomial_iterator_rep<C,V> (P));
00301 }
00302
00303
00304
00305
00306
00307 TMPL inline generic var (const Polynomial& P) {
00308 (void) P; return Polynomial::get_variable_name (); }
00309
00310 TMPL inline void set_variable_name (const Polynomial& P, const generic& x) {
00311 (void) P; return Polynomial::set_variable_name (x); }
00312
00313 TMPL syntactic
00314 flatten (const Polynomial& P, const syntactic& v) {
00315 int i, len= deg(P);
00316 syntactic s= 0;
00317 for (i=len; i>=0; i--)
00318 s= s + flatten (P[i]) * pow (v, i);
00319 return s;
00320 }
00321
00322 TMPL syntactic
00323 flatten (const Polynomial& P) {
00324 return flatten (P, as_syntactic (var (P)));
00325 }
00326
00327 TMPL
00328 struct binary_polynomial_helper: public void_binary_helper<Polynomial > {
00329 static inline string short_type_name () {
00330 return "P" * Short_type_name (C); }
00331 static inline generic full_type_name () {
00332 return gen ("Polynomial", Full_type_name (C)); }
00333 };
00334
00335 TMPL
00336 struct binary_helper<Polynomial >: public binary_polynomial_helper<C,V> {
00337 static inline generic disassemble (const Polynomial& p) {
00338 return binary_disassemble<vector<C> > (coefficients (p)); }
00339 static inline Polynomial assemble (const generic& v) {
00340 return Polynomial (binary_assemble<vector<C> > (v)); }
00341 static inline void write (const port& out, const Polynomial& p) {
00342 binary_write<Format > (out, CF(p));
00343 binary_write<nat> (out, N(p));
00344 for (nat i=0; i<N(p); i++)
00345 binary_write<C> (out, p[i]); }
00346 static inline Polynomial read (const port& in) {
00347 Format fm= binary_read<Format > (in);
00348 nat n= binary_read<nat> (in);
00349 nat l= aligned_size<C,V> (n);
00350 C* r= mmx_formatted_new<C> (l, fm);
00351 for (nat i=0; i<n; i++)
00352 r[i]= binary_read<C> (in);
00353 return Polynomial (r, n, l, fm); }
00354 };
00355
00356
00357
00358
00359
00360 TMPL int
00361 val (const Polynomial& P) {
00362 for (nat i=0; i<N(P); i++)
00363 if (P[i] != 0) return (int) i;
00364 return (int) (((nat) (-1)) >> 1);
00365 }
00366
00367 TMPL int
00368 sign (const Polynomial& P) {
00369 for (nat i=0; i<N(P); i++) {
00370 int s= sign (P[i]);
00371 if (s != 0) return s;
00372 }
00373 return 0;
00374 }
00375
00376 TMPL int
00377 compare (const Polynomial& P1, const Polynomial& P2) {
00378 nat n1= N(P1), n2= N(P2), n= max (n1, n2);
00379 for (nat i=0; i<n; i++) {
00380 int s= sign (P1[i] - P2[i]);
00381 if (s != 0) return s;
00382 }
00383 return 0;
00384 }
00385
00386 template<typename Op,typename C,typename V> inline bool
00387 binary_test (const Polynomial& P1, const Polynomial& P2) {
00388 nat n1= N(P1), n2= N(P2), n= max (n1, n2);
00389 for (nat i=0; i<n; i++)
00390 if (Op::not_op (P1[i], P2[i])) return false;
00391 return true;
00392 }
00393
00394 TRUE_IDENTITY_OP_SUGAR(TMPL,Polynomial)
00395 EXACT_IDENTITY_OP_SUGAR(TMPL,Polynomial)
00396 COMPARE_SUGAR(TMPL,Polynomial)
00397 EQUAL_SCALAR_SUGAR(TMPL,Polynomial,C)
00398 COMPARE_SCALAR_SUGAR(TMPL,Polynomial,C)
00399 EQUAL_INT_SUGAR(TMPL,Polynomial)
00400 COMPARE_INT_SUGAR(TMPL,Polynomial)
00401 EQUAL_SCALAR_SUGAR_BIS(TMPL,Polynomial,C)
00402 COMPARE_SCALAR_SUGAR_BIS(TMPL,Polynomial,C)
00403
00404
00405
00406
00407
00408 TMPL Polynomial&
00409 Polynomial::operator += (const Polynomial& P) {
00410 typedef implementation<polynomial_linear,V> Pol;
00411 secure ();
00412 nat n1= rep->n, n2= N(P);
00413 if (n1 < n2) rep -> extend (n2);
00414 Pol::add (rep->a, seg (P), n2);
00415 if (n1 == n2) rep -> normalize ();
00416 return *this;
00417 }
00418
00419 TMPL Polynomial&
00420 Polynomial::operator -= (const Polynomial& P) {
00421 typedef implementation<polynomial_linear,V> Pol;
00422 secure ();
00423 nat n1= rep->n, n2= N(P);
00424 if (n1 < n2) rep->extend (n2);
00425 Pol::sub (rep->a, seg (P), n2);
00426 if (n1 == n2) rep->normalize ();
00427 return *this;
00428 }
00429
00430 TMPL Polynomial&
00431 Polynomial::operator *= (const C& c) {
00432 typedef implementation<polynomial_linear,V> Pol;
00433 secure ();
00434 Pol::mul_sc (rep->a, c, rep->n);
00435 rep->normalize ();
00436 return *this;
00437 }
00438
00439 TMPL Polynomial&
00440 Polynomial::operator /= (const C& c) {
00441 typedef implementation<polynomial_linear,V> Pol;
00442 secure ();
00443 Pol::div_sc (rep->a, c, rep->n);
00444 rep->normalize ();
00445 return *this;
00446 }
00447
00448
00449
00450
00451
00452 TMPL Polynomial
00453 operator - (const Polynomial& P) {
00454 typedef implementation<polynomial_linear,V> Pol;
00455 nat n= N(P);
00456 nat l= aligned_size<C,V> (n);
00457 C* r= mmx_formatted_new<C> (l, CF(P));
00458 Pol::neg (r, seg (P), n);
00459 return Polynomial (r, n, l, CF(P));
00460 }
00461
00462 TMPL Polynomial
00463 operator + (const Polynomial& P1, const Polynomial& P2) {
00464 typedef implementation<polynomial_linear,V> Pol;
00465 nat m= min (N(P1), N(P2));
00466 nat n= max (N(P1), N(P2));
00467 nat l= aligned_size<C,V> (n);
00468 C* r= mmx_formatted_new<C> (l, CF(P1));
00469 Pol::add (r, seg (P1), seg (P2), m);
00470 if (N(P1) > m) Pol::copy (r+m, seg(P1)+m, n-m);
00471 if (N(P2) > m) Pol::copy (r+m, seg(P2)+m, n-m);
00472 return Polynomial (r, n, l, CF(P1));
00473 }
00474
00475 TMPL Polynomial operator + (const Polynomial& P, const C& c) {
00476 return P + Polynomial (c); }
00477 TMPL Polynomial operator + (const C& c, const Polynomial& P) {
00478 return Polynomial (c) + P; }
00479
00480 TMPL Polynomial
00481 big_add (const vector<Polynomial >& a) {
00482 typedef implementation<polynomial_linear,V> Pol;
00483 nat i, k= N(a), n=0;
00484 nat l= aligned_size<C,V> (n);
00485 for (i=0; i<k; i++) n= max (N(a[i]), n);
00486 C* r= mmx_formatted_new<C> (l, CF (get_sample (CF (a))));
00487 Pol::set (r, 0, n);
00488 for (i=0; i<k; i++) Pol::add (r, seg (a[i]), N(a[i]));
00489 return Polynomial (r, n, l, Format (CF(a)));
00490 }
00491
00492 TMPL Polynomial
00493 operator - (const Polynomial& P1, const Polynomial& P2) {
00494 typedef implementation<polynomial_linear,V> Pol;
00495 nat m= min (N(P1), N(P2));
00496 nat n= max (N(P1), N(P2));
00497 nat l= aligned_size<C,V> (n);
00498 C* r= mmx_formatted_new<C> (l, CF(P1));
00499 Pol::sub (r, seg (P1), seg (P2), m);
00500 if (N(P1) > m) Pol::copy (r+m, seg(P1)+m, n-m);
00501 if (N(P2) > m) Pol::neg (r+m, seg(P2)+m, n-m);
00502 return Polynomial (r, n, l, CF(P1));
00503 }
00504
00505 TMPL Polynomial operator - (const Polynomial& P, const C& c) {
00506 return P - Polynomial (c); }
00507 TMPL Polynomial operator - (const C& c, const Polynomial& P) {
00508 return Polynomial (c) - P; }
00509
00510 TMPL Polynomial
00511 operator * (const Polynomial& P1, const Polynomial& P2) {
00512 typedef implementation<polynomial_multiply,V> Pol;
00513 nat n1= N(P1), n2= N(P2);
00514 if (n1 == 0 || n2 == 0) return Polynomial (CF(P1));
00515 nat l= aligned_size<C,V> (n1+n2-1);
00516 C* r= mmx_formatted_new<C> (l, CF(P1));
00517 Pol::mul (r, seg (P1), seg (P2), n1, n2);
00518 return Polynomial (r, n1+n2-1, l, CF(P1));
00519 }
00520
00523 TMPL Polynomial
00524 tmul (int d2, const Polynomial& P1, const Polynomial& P2) {
00525 typedef implementation<polynomial_multiply,V> Pol;
00526 nat n2 = max (0, d2 + 1), n1= N(P1), n= N(P2);
00527 if (n1 == 0 || n2 == 0) return Polynomial (CF(P1));
00528 ASSERT (n < n1 + n2, "bad dimension in tmul");
00529 nat l2= aligned_size<C,V> (n2);
00530 C* r= mmx_formatted_new<C> (l2, CF(P1));
00531 if (n != n1 + n2 - 1) {
00532 nat l= aligned_size<C,V> (n1+n2-1);
00533 C* s2= mmx_formatted_new<C> (l, CF(P1));
00534 Pol::copy (s2, seg (P2), n);
00535 Pol::clear (s2 + n, n1+n2-n-1);
00536 Pol::tmul (r, seg (P1), s2, n1, n2);
00537 mmx_delete<C> (s2, l);
00538 }
00539 else
00540 Pol::tmul (r, seg (P1), seg (P2), n1, n2);
00541 return Polynomial (r, n2, l2, CF(P1));
00542 }
00543
00544 TMPL Polynomial
00545 square (const Polynomial& P) {
00546 typedef implementation<polynomial_multiply,V> Pol;
00547 nat n= N(P);
00548 if (n == 0) return P;
00549 nat l= aligned_size<C,V> (2*n-1);
00550 C* r= mmx_formatted_new<C> (l, CF(P));
00551 Pol::square (r, seg (P), n);
00552 return Polynomial (r, 2*n-1, l, CF(P));
00553 }
00554
00555 TMPL Polynomial
00556 operator * (const Polynomial& P, const C& c) {
00557 typedef implementation<polynomial_linear,V> Pol;
00558 nat n= N(P);
00559 nat l= aligned_size<C,V> (n);
00560 C* r= mmx_formatted_new<C> (l, CF(P));
00561 Pol::mul_sc (r, seg (P), c, n);
00562 return Polynomial (r, n, l, CF(P));
00563 }
00564
00565 TMPL Polynomial
00566 operator * (const C& c, const Polynomial& P) {
00567 typedef implementation<polynomial_linear,V> Pol;
00568 nat n= N(P);
00569 nat l= aligned_size<C,V> (n);
00570 C* r= mmx_formatted_new<C> (l, CF(P));
00571 Pol::mul_sc (r, seg (P), c, n);
00572 return Polynomial (r, n, l, CF(P));
00573 }
00574
00575 TMPL Polynomial
00576 big_mul (const vector<Polynomial >& a) {
00577 typedef implementation<vector_linear,V> Vec;
00578 return Vec::template vec_unary_big_dicho<mul_op> (seg (a), N(a));
00579 }
00580
00581 TMPL Polynomial
00582 operator / (const Polynomial& P, const C& c) {
00583 typedef implementation<polynomial_linear,V> Pol;
00584 nat n= N(P);
00585 nat l= aligned_size<C,V> (n);
00586 C* r= mmx_formatted_new<C> (l, CF(P));
00587 Pol::div_sc (r, seg (P), c, n);
00588 return Polynomial (r, n, l, CF(P));
00589 }
00590
00591 TMPL Polynomial
00592 operator / (const Polynomial& P1, const Polynomial& P2) {
00593 typedef implementation<polynomial_exact_divide,V> Pol;
00594 nat n1= N(P1), n2= N(P2);
00595 ASSERT (n2 != 0, "division by zero");
00596 if (n1 < n2) return Polynomial (CF(P1));
00597 nat lq= aligned_size<C,V> (n1-n2+1);
00598 C* q= mmx_formatted_new<C> (lq, CF(P1));
00599 Pol::div (q, seg (P1), seg (P2), n1, n2);
00600 return Polynomial (q, n1-n2+1, lq, CF(P1));
00601 }
00602
00603 TMPL Polynomial
00604 skew_div (const Polynomial& P, const C& c, bool left) {
00605 typedef implementation<polynomial_linear,V> Pol;
00606 if (left) {
00607 nat n= N(P);
00608 nat l= aligned_size<C,V> (n);
00609 C* r= mmx_formatted_new<C> (n, CF(P));
00610 Pol::div_sc (r, seg (P), c, n);
00611 return Polynomial (r, n, l, CF(P));
00612 }
00613 else return P / c;
00614 }
00615
00616 TMPL Polynomial
00617 quo (const Polynomial& P, const C& c) {
00618 typedef implementation<polynomial_linear,V> Pol;
00619 nat n= N(P);
00620 nat l= aligned_size<C,V> (n);
00621 C* r= mmx_formatted_new<C> (l, CF(P));
00622 Pol::quo_sc (r, seg (P), c, n);
00623 return Polynomial (r, n, l, CF(P));
00624 }
00625
00626 TMPL Polynomial
00627 rem (const Polynomial& P, const C& c) {
00628 typedef implementation<polynomial_linear,V> Pol;
00629 nat n= N(P);
00630 nat l= aligned_size<C,V> (n);
00631 C* r= mmx_formatted_new<C> (l, CF(P));
00632 Pol::rem_sc (r, seg (P), c, n);
00633 return Polynomial (r, n, l, CF(P));
00634 }
00635
00636 TMPL Polynomial
00637 quo (const Polynomial& P1, const Polynomial& P2) {
00638 typedef implementation<polynomial_divide,V> Pol;
00639 nat n1= N(P1), n2= N(P2);
00640 if (n1 < n2 || n2 == 0) return Polynomial (CF(P1));
00641 nat lq= aligned_size<C,V> (n1-n2+1);
00642 nat lr= aligned_size<C,V> (n1);
00643 C* q= mmx_formatted_new<C> (lq, CF(P1));
00644 C* r= mmx_formatted_new<C> (lr, CF(P1));
00645 Pol::copy (r, seg (P1), n1);
00646 Pol::quo_rem (q, r, seg (P2), n1, n2);
00647 mmx_delete<C> (r, lr);
00648 return Polynomial (q, n1-n2+1, lq, CF(P1));
00649 }
00650
00652 TMPL Polynomial
00653 tquo (int d1, const Polynomial& P1, const Polynomial& P2) {
00654 typedef implementation<polynomial_divide,V> Pol;
00655 nat n1= max (0, d1+1), n= N(P1), n2= N(P2);
00656 ASSERT (n <= n1-n2+1, "bad dimension in tquo");
00657 if (n1 < n2 || n2 == 0) return Polynomial (CF(P1));
00658 nat l= aligned_size<C,V> (n1);
00659 C* q= mmx_formatted_new<C> (l, CF(P1));
00660 C* r= mmx_formatted_new<C> (l, CF(P1));
00661 Pol::clear (r, n1);
00662 Pol::copy (r+n2-1, seg (P1), n);
00663 Pol::tquo_rem (q, r, seg (P2), n1, n2);
00664 mmx_delete<C> (r, l);
00665 return Polynomial (q, n1, l, CF(P1));
00666 }
00667
00668 TMPL Polynomial
00669 rem (const Polynomial& P1, const Polynomial& P2) {
00670 typedef implementation<polynomial_divide,V> Pol;
00671 nat n1= N(P1), n2= N(P2);
00672 if (n1 < n2 || n2 == 0) return P1;
00673 nat lq= aligned_size<C,V> (n1-n2+1);
00674 nat lr= aligned_size<C,V> (n1);
00675 C* q= mmx_formatted_new<C> (lq, CF(P1));
00676 C* r= mmx_formatted_new<C> (lr, CF(P1));
00677 Pol::copy (r, seg (P1), n1);
00678 Pol::quo_rem (q, r, seg (P2), n1, n2);
00679 mmx_delete<C> (q, lq);
00680 return Polynomial (r, n2 - 1, lr, CF(P1));
00681 }
00682
00683 TMPL bool
00684 divides (const Polynomial& P1, const Polynomial& P2) {
00685 return rem (P2, P1) == 0;
00686 }
00687
00689 TMPL Polynomial
00690 rem (const Polynomial& P1, const Polynomial& P2, Polynomial& Q) {
00691 typedef implementation<polynomial_divide,V> Pol;
00692 nat n1= N(P1), n2= N(P2);
00693 if (n1 < n2 || n2 == 0) return P1;
00694 nat lq= aligned_size<C,V> (n1-n2+1);
00695 nat lr= aligned_size<C,V> (n1);
00696 C* q= mmx_formatted_new<C> (lq, CF(P1));
00697 C* r= mmx_formatted_new<C> (lr, CF(P1));
00698 Pol::copy (r, seg (P1), n1);
00699 Pol::quo_rem (q, r, seg (P2), n1, n2);
00700 Q= Polynomial (q, n1-n2+1, lq, CF(P1));
00701 return Polynomial (r, n2 - 1, lr, CF(P1));
00702 }
00703
00705 TMPL Polynomial
00706 trem (int d1, const Polynomial& P1, const Polynomial& P2) {
00707 typedef implementation<polynomial_divide,V> Pol;
00708 nat n1= max (0, d1+1), n= N(P1), n2= N(P2);
00709 ASSERT (n <= n2-1, "bad dimension in trem");
00710 if (n1 < n2 || n2 == 0) return P1;
00711 nat l= aligned_size<C,V> (n1);
00712 C* q= mmx_formatted_new<C> (l, CF(P1));
00713 C* r= mmx_formatted_new<C> (l, CF(P1));
00714 Pol::copy (r, seg (P1), n);
00715 Pol::clear (r+n, n1-n);
00716 Pol::tquo_rem (q, r, seg (P2), n1, n2);
00717 mmx_delete<C> (r, l);
00718 return Polynomial (q, n1, l, CF(P1));
00719 }
00720
00723 TMPL inline int
00724 pexponent (const Polynomial& P1, const Polynomial& P2) {
00725 return max (0, deg (P1) - deg (P2) + 1);
00726 }
00727
00728 TMPL Polynomial
00729 pquo (const Polynomial& P1, const Polynomial& P2) {
00730 typedef implementation<polynomial_divide,V> Pol;
00731 nat n1= N(P1), n2= N(P2);
00732 if (n1 < n2 || n2 == 0) return Polynomial (CF(P1));
00733 nat lq= aligned_size<C,V> (n1-n2+1);
00734 nat lr= aligned_size<C,V> (n1);
00735 C* q= mmx_formatted_new<C> (lq, CF(P1));
00736 C* r= mmx_formatted_new<C> (lr, CF(P1));
00737 Pol::copy (r, seg (P1), n1);
00738 Pol::pquo_rem (q, r, seg (P2), n1, n2);
00739 mmx_delete<C> (r, lr);
00740 return Polynomial (q, n1-n2+1, lq, CF(P1));
00741 }
00742
00743 TMPL Polynomial
00744 prem (const Polynomial& P1, const Polynomial& P2) {
00745 typedef implementation<polynomial_divide,V> Pol;
00746 nat n1= N(P1), n2= N(P2);
00747 if (n1 < n2 || n2 == 0) return P1;
00748 nat lq= aligned_size<C,V> (n1-n2+1);
00749 nat lr= aligned_size<C,V> (n1);
00750 C* q= mmx_formatted_new<C> (lq, CF(P1));
00751 C* r= mmx_formatted_new<C> (lr, CF(P1));
00752 Pol::copy (r, seg (P1), n1);
00753 Pol::pquo_rem (q, r, seg (P2), n1, n2);
00754 mmx_delete<C> (q, lq);
00755 return Polynomial (r, n2 - 1, lr, CF(P1));
00756 }
00757
00759 TMPL Polynomial
00760 prem (const Polynomial& P1, const Polynomial& P2, Polynomial& Q) {
00761 typedef implementation<polynomial_divide,V> Pol;
00762 nat n1= N(P1), n2= N(P2);
00763 if (n1 < n2 || n2 == 0) return Polynomial (CF(P1));
00764 nat lq= aligned_size<C,V> (n1-n2+1);
00765 nat lr= aligned_size<C,V> (n1);
00766 C* q= mmx_formatted_new<C> (lq, CF(P1));
00767 C* r= mmx_formatted_new<C> (lr, CF(P1));
00768 Pol::copy (r, seg (P1), n1);
00769 Pol::pquo_rem (q, r, seg (P2), n1, n2);
00770 Q= Polynomial (q, n1-n2+1, lq, CF(P1));
00771 return Polynomial (r, n2 - 1, lr, CF(P1));
00772 }
00773
00774 TMPL inline Polynomial operator / (const C& c, const Polynomial& P) {
00775 return Polynomial (c) / P; }
00776
00777 ARITH_SCALAR_INT_SUGAR(TMPL,Polynomial)
00778
00779
00780
00781
00782
00783
00784
00785
00786 TMPL inline C
00787 contents (const Polynomial& P) {
00788 typedef implementation<polynomial_linear,V> Pol;
00789 return Pol::big_gcd (seg (P), N(P));
00790 }
00791
00792 TMPL inline Polynomial
00793 primitive_part (const Polynomial& P, C& c) {
00794 typedef implementation<polynomial_linear,V> Pol;
00795 c= contents (P);
00796 nat n= N(P);
00797 nat l= aligned_size<C,V> (n);
00798 C* r= mmx_formatted_new<C> (l, CF(P));
00799 Pol::div_sc (r, seg (P), c, n);
00800 return Polynomial (r, n, l, CF(P));
00801 }
00802
00803 TMPL inline Polynomial
00804 primitive_part (const Polynomial& P) {
00805 C c;
00806 return primitive_part (P, c);
00807 }
00808
00809 TMPL inline Polynomial
00810 gcd (const Polynomial& P1, const Polynomial& P2,
00811 Polynomial& U1, Polynomial& U2) {
00812 typedef implementation<polynomial_gcd,V> Pol;
00813 C one= promote (1, CF(P1));
00814 U1= one; U2= one;
00815 return Pol::gcd (P1, P2, U1, U2);
00816 }
00817
00818 TMPL inline Polynomial
00819 gcd (const Polynomial& P1, const Polynomial& P2, Polynomial& U1) {
00820 Polynomial U2(0); U1= promote (1, CF(P1));
00821 return gcd (P1, P2, U1, U2);
00822 }
00823
00824 TMPL inline Polynomial
00825 gcd (const Polynomial& P1, const Polynomial& P2) {
00826 typedef implementation<polynomial_gcd,V> Pol;
00827 return Pol::gcd (P1, P2);
00828 }
00829
00830 TMPL static vector<Polynomial>
00831 gcd (const Polynomial& p, const vector<Polynomial>& q) {
00832 typedef implementation<polynomial_evaluate,V> Pol;
00833 return Pol::multi_gcd (p, q);
00834 }
00835
00836 TMPL inline Polynomial
00837 lcm (const Polynomial& P1, const Polynomial& P2) {
00838 typedef implementation<polynomial_gcd,V> Pol;
00839 return Pol::lcm (P1, P2);
00840 }
00841
00842 TMPL inline Polynomial
00843 invert_modulo (const Polynomial& P, const Polynomial& Q) {
00844 typedef implementation<polynomial_gcd,V> Pol;
00845 return Pol::invert_mod (P, Q);
00846 }
00847
00848 TMPL inline bool
00849 reconstruct (const Polynomial& P, const Polynomial& Q, nat k,
00850 Polynomial& Num, Polynomial &Den) {
00851 typedef implementation<polynomial_gcd,V> Pol;
00852 if (deg (Q) <= 0) {
00853 if (deg (P) < (int) k) {
00854 Num= P;
00855 Den= Polynomial (promote (1, CF(P)));
00856 return true;
00857 }
00858 else return false;
00859 }
00860 nat n= deg (Q);
00861 ASSERT (N (P) <= n+1 && k <= n && k > 0, "invalid argument");
00862 Pol::reconstruct (P, Q, k, Num, Den);
00863 return deg (gcd (Num, Den)) == 1;
00864 }
00865
00866 TMPL inline void
00867 pade (const Polynomial& P, nat n, nat k, Polynomial& Num, Polynomial &Den) {
00868 typedef implementation<polynomial_gcd,V> Pol;
00869 ASSERT (N (P) <= n+1 && k <= n && k > 0, "invalid argument");
00870 Polynomial Q (promote (1, CF(P)), n);
00871 Pol::reconstruct (P, Q, k, Num, Den);
00872 }
00873
00874 template<typename C,typename V,typename W> inline void
00875 minimal_polynomial_bis (Polynomial& p, const vector<C,W>& v) {
00876
00877
00878 nat n= N(v); nat k= n >> 1;
00879 ASSERT ((n & 1) == 0, "even size expected");
00880 Polynomial h (v), s, t;
00881 pade (h, n, k, s, t);
00882 t= reverse (t);
00883 p= (deg(t) < 1 + deg(s)) ? lshiftz (t, 1 + deg (s) - deg (t)) : t;
00884 }
00885
00886 template<typename C,typename W> inline
00887 polynomial<C,typename Polynomial_variant(C) >
00888 minimal_polynomial (const vector<C,W>& v) {
00889 polynomial<C,typename Polynomial_variant(C) > P;
00890 minimal_polynomial_bis (P, v);
00891 return P;
00892 }
00893
00894
00895
00896
00897
00898
00899
00900
00901 TMPL inline vector<Polynomial>
00902 subresultants (const Polynomial& P, const Polynomial& Q,
00903 vector<Polynomial>& co_P, vector<Polynomial>& co_Q) {
00904 typedef implementation<polynomial_subresultant,V> Pol;
00905 int n= deg (P), m= deg (Q); nat l= max (min (n, m), 0);
00906 C one= promote (1, CF(P));
00907 Polynomial d; vector<Polynomial> res (Polynomial (one), l);
00908 co_P= vector<Polynomial> (Polynomial (one), l);
00909 co_Q= vector<Polynomial> (Polynomial (one), l);
00910 Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
00911 d, d, d, d, d, d, 0);
00912 return res;
00913 }
00914
00915 TMPL inline vector<Polynomial>
00916 subresultants (const Polynomial& P, const Polynomial& Q,
00917 vector<Polynomial>& co_P) {
00918 typedef implementation<polynomial_subresultant,V> Pol;
00919 int n= deg (P), m= deg (Q); nat l= max (min (n, m), 0);
00920 Polynomial d;
00921 C zero= promote (0, CF(P)), one= promote (1, CF(P));
00922 vector<Polynomial> res (Polynomial (one), l), co_Q (Polynomial (zero), 0);
00923 co_P= vector<Polynomial> (Polynomial (one), l);
00924 Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
00925 d, d, d, d, d, d, 0);
00926 return res;
00927 }
00928
00929 TMPL inline vector<Polynomial>
00930 subresultants (const Polynomial& P, const Polynomial& Q) {
00931 typedef implementation<polynomial_subresultant,V> Pol;
00932 int n= deg (P), m= deg (Q); nat l= max (min (n, m), 0);
00933 Polynomial d;
00934 C zero= promote (0, CF(P)), one= promote (1, CF(P));
00935 vector<Polynomial>
00936 res (Polynomial(one), l),
00937 co_P (Polynomial (zero), 0),
00938 co_Q (Polynomial (zero), 0);
00939 Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
00940 d, d, d, d, d, d, 0);
00941 return res;
00942 }
00943
00944 TMPL inline Polynomial
00945 subresultant (const Polynomial& P, const Polynomial& Q, int k,
00946 Polynomial& coP, Polynomial& coQ) {
00947 typedef implementation<polynomial_subresultant,V> Pol;
00948 int n= deg (P), m= deg (Q);
00949 nat l= max (min (n, m), 0);
00950 ASSERT (k < l, "index out of range");
00951 C zero= promote (0, CF(P)), one= promote (1, CF(P));
00952 Polynomial d;
00953 vector<Polynomial>
00954 res (Polynomial (zero), l),
00955 co_P (Polynomial (zero), l),
00956 co_Q (Polynomial (zero), l);
00957 res [k]= Polynomial (one);
00958 co_P[k]= Polynomial (one);
00959 co_Q[k]= Polynomial (one);
00960 Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
00961 d, d, d, d, d, d, 0);
00962 coP= co_P[k]; coQ= co_Q[k];
00963 return res[k];
00964 }
00965
00966 TMPL inline Polynomial
00967 subresultant (const Polynomial& P, const Polynomial& Q, int k,
00968 Polynomial& coP) {
00969 typedef implementation<polynomial_subresultant,V> Pol;
00970 int n= deg (P), m= deg (Q);
00971 nat l= max (min (n, m), 0);
00972 ASSERT (k < l, "index out of range");
00973 Polynomial d;
00974 C zero= promote (0, CF(P)), one= promote (1, CF(P));
00975 vector<Polynomial>
00976 res (Polynomial (zero), l),
00977 co_P (Polynomial (zero), l),
00978 co_Q (Polynomial (zero), 0);
00979 res [k]= Polynomial (one);
00980 co_P[k]= Polynomial (one);
00981 Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
00982 d, d, d, d, d, d, 0);
00983 coP= co_P[k];
00984 return res[k];
00985 }
00986
00987 TMPL inline Polynomial
00988 subresultant (const Polynomial& P, const Polynomial& Q, int k) {
00989 typedef implementation<polynomial_subresultant,V> Pol;
00990 if (k < 0) return promote (0, P);
00991 int n= deg (P), m= deg (Q);
00992 nat l= max (min (n, m), 0);
00993 ASSERT ((nat) k < l, "index out of range");
00994 Polynomial d;
00995 C zero= promote (0, CF(P)), one= promote (1, CF(P));
00996 vector<Polynomial>
00997 res (Polynomial (zero), l),
00998 co_P (Polynomial (zero), 0),
00999 co_Q (Polynomial (zero), 0);
01000 res [k]= Polynomial (one);
01001 Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
01002 d, d, d, d, d, d, 0);
01003 return res[k];
01004 }
01005
01006 TMPL inline C
01007 resultant (const Polynomial& P, const Polynomial& Q) {
01008 typedef implementation<polynomial_subresultant,V> Pol;
01009 int n= degree (P), m= degree (Q);
01010 if (n < 0 || m < 0) return promote (0, CF(P));
01011 if (m == 0) return binpow (Q[0], n);
01012 if (n == 0) {
01013 C r= binpow (P[0], m);
01014 return (n & 1) ? -r : r;
01015 }
01016 Polynomial d;
01017 C zero= promote (0, CF(P)), one= promote (1, CF(P));
01018 vector<Polynomial>
01019 res (Polynomial (one), 1),
01020 co_P (Polynomial (zero), 0),
01021 co_Q (Polynomial (zero), 0);
01022 Pol::subresultant_sequence (P, Q, res, co_P, co_Q,
01023 d, d, d, d, d, d, 0);
01024 return res[0][0];
01025 }
01026
01027 TMPL inline C
01028 discriminant (const Polynomial& P) {
01029 typedef implementation<polynomial_subresultant,V> Pol;
01030 return resultant (derive (P), P);
01031 }
01032
01033
01034
01035
01036
01037 TMPL Polynomial
01038 derive (const Polynomial& P) {
01039 typedef implementation<polynomial_linear,V> Pol;
01040 nat n= N(P);
01041 if (n <= 1) return promote (0, P);
01042 nat l= aligned_size<C,V> (n-1);
01043 C* r= mmx_formatted_new<C> (l, CF(P));
01044 Pol::derive (r, seg (P), n);
01045 return Polynomial (r, n-1, l, CF(P));
01046 }
01047
01048 TMPL Polynomial
01049 derive (const Polynomial& P, const nat& order) {
01050 typedef implementation<polynomial_linear,V> Pol;
01051 nat n= N(P);
01052 if (n <= order) return promote (0, P);
01053 nat l= aligned_size<C,V> (n-order);
01054 C* r= mmx_formatted_new<C> (l, CF(P));
01055 Pol::derive (r, seg (P), n, order);
01056 return Polynomial (r, n-order, l, CF(P));
01057 }
01058
01059 TMPL Polynomial
01060 xderive (const Polynomial& P) {
01061 typedef implementation<polynomial_linear,V> Pol;
01062 nat n= N(P);
01063 nat l= aligned_size<C,V> (n);
01064 C* r= mmx_formatted_new<C> (l, CF(P));
01065 Pol::xderive (r, seg (P), n);
01066 return Polynomial (r, n, l, CF(P));
01067 }
01068
01069 TMPL Polynomial
01070 integrate (const Polynomial& P) {
01071 typedef implementation<polynomial_linear,V> Pol;
01072 nat n= N(P);
01073 nat l= aligned_size<C,V> (n+1);
01074 C* r= mmx_formatted_new<C> (l, CF(P));
01075 Pol::integrate (r, seg (P), n);
01076 return Polynomial (r, n+1, l, CF(P));
01077 }
01078
01079 template<typename C,typename V,typename K> Polynomial
01080 compose (const Polynomial& P1, const polynomial<K,V>& P2) {
01081 typedef implementation<polynomial_compose,V> Pol;
01082 nat n1= N(P1), n2= N(P2);
01083 if (n1 <= 1) return P1;
01084 if (n2 == 0) return P1[0];
01085 nat n= (n1-1) * (n2-1) + 1;
01086 nat l= aligned_size<C,V> (n);
01087 C* r= mmx_formatted_new<C> (l, CF(P1));
01088 Pol::compose (r, seg (P1), seg (P2), n1, n2);
01089 return Polynomial (r, n, l, CF(P1));
01090 }
01091
01092 TMPL Polynomial
01093 shift (const Polynomial& P, const C& sh) {
01094 typedef implementation<polynomial_compose,V> Pol;
01095 nat n= N(P);
01096 if (n <= 1 || sh == 0) return P;
01097 nat l= aligned_size<C,V> (n);
01098 C* r= mmx_formatted_new<C> (l, CF(P));
01099 Pol::shift (r, seg (P), sh, n);
01100 return Polynomial (r, n, l, CF(P));
01101 }
01102
01103 TMPL Polynomial
01104 shift (const Polynomial& P, int i) {
01105 return shift (P, promote (i, CF(P)));
01106 }
01107
01108 template<typename C,typename V,typename K> Polynomial
01109 q_difference (const Polynomial& P, const K& q) {
01110 typedef implementation<polynomial_linear,V> Pol;
01111 nat n= N(P);
01112 if (n <= 1) return P;
01113 nat l= aligned_size<C,V> (n);
01114 C* r= mmx_formatted_new<C> (l, CF(P));
01115 Pol::q_difference (r, seg (P), q, n);
01116 return Polynomial (r, n, l, CF(P));
01117 }
01118
01119 TMPL Polynomial
01120 dilate (const Polynomial& P, nat p) {
01121 typedef implementation<polynomial_linear,V> Pol;
01122 if (p == 1) return P;
01123 nat n= N(P);
01124 if (n <= 1) return P;
01125 nat k= (n-1)*p + 1;
01126 nat l= aligned_size<C,V> (k);
01127 C* r= mmx_formatted_new<C> (l, CF(P));
01128 Pol::dilate (r, seg (P), p, n);
01129 return Polynomial (r, k, l, CF(P));
01130 }
01131
01132
01133
01134
01135
01136 TMPL inline vector<Polynomial>
01137 rem (const Polynomial& p, const vector<Polynomial>& q) {
01138 typedef implementation<polynomial_evaluate,V> Pol;
01139 return Pol::multi_rem (p, q);
01140 }
01141
01142 TMPL inline vector<Polynomial>
01143 prem (const Polynomial& p, const vector<Polynomial>& q) {
01144 typedef implementation<polynomial_evaluate,V> Pol;
01145 return Pol::multi_prem (p, q);
01146 }
01147
01148 TMPL inline Polynomial
01149 annulator_bis (const Vector& x) {
01150 typedef implementation<polynomial_evaluate,V> Pol;
01151 return Pol::template annulator<Polynomial> (x);
01152 }
01153
01154 template<typename C> inline polynomial<C,typename Polynomial_variant(C) >
01155 annulator (const Vector& x) {
01156 return annulator_bis<C,typename Polynomial_variant(C) > (x);
01157 }
01158
01159 BINARY_RETURN_TYPE(TMPL,evaluate_op,Polynomial,C,C);
01160
01161 TMPL inline C
01162 evaluate (const Polynomial& p, const C& x) {
01163 typedef implementation<polynomial_evaluate,V> Pol;
01164 return Pol::evaluate (seg (p), x, N(p));
01165 }
01166
01167 TMPL inline bool
01168 is_evaluable (const Polynomial& p, const C& x, C& y) {
01169 y= evaluate (p, x);
01170 return true;
01171 }
01172
01173 TMPL inline C
01174 eval (const Polynomial& p, const C& x) {
01175 return evaluate (p, x);
01176 }
01177
01178 TMPL inline Polynomial
01179 tevaluate_bis (const C& v, const C& x, nat n) {
01180 typedef implementation<polynomial_evaluate,V> Pol;
01181 nat l= aligned_size<C,V> (n);
01182 C* buf= mmx_formatted_new<C> (l, get_format (v));
01183 Pol::tevaluate (v, buf, x, n);
01184 return Polynomial (buf, n, l, get_format (v));
01185 }
01186
01187 template<typename C> inline polynomial<C,typename Polynomial_variant(C) >
01188 tevaluate (const C& v, const C& x, nat n) {
01189 return tevaluate_bis<C,typename Polynomial_variant(C) > (v, x, n);
01190 }
01191
01192 TMPL inline Vector
01193 evaluate (const Polynomial& p, const Vector& x) {
01194 typedef implementation<polynomial_evaluate,V> Pol;
01195 return Pol::evaluate (p, x);
01196 }
01197
01198 TMPL inline Vector
01199 eval (const Polynomial& p, const Vector& x) {
01200 return evaluate (p, x);
01201 }
01202
01203 TMPL inline Polynomial
01204 tevaluate_bis (const Vector& v, const Vector& x, nat l) {
01205 typedef implementation<polynomial_evaluate,V> Pol;
01206 return Pol::template tevaluate<Polynomial> (v, x, l);
01207 }
01208
01209 template<typename C> inline polynomial<C,typename Polynomial_variant(C) >
01210 tevaluate (const Vector& v, const Vector& x, nat l) {
01211 return tevaluate_bis<C,typename Polynomial_variant(C) > (v, x, l);
01212 }
01213
01214 TMPL inline Polynomial
01215 interpolate_bis (const Vector& v, const Vector& x) {
01216 typedef implementation<polynomial_evaluate,V> Pol;
01217 return Pol::template interpolate<Polynomial> (v, x);
01218 }
01219
01220 template<typename C> inline polynomial<C,typename Polynomial_variant(C) >
01221 interpolate (const Vector& v, const Vector& x) {
01222 return interpolate_bis<C,typename Polynomial_variant(C) > (v, x);
01223 }
01224
01225 TMPL Vector
01226 tinterpolate (const Polynomial& p, const Vector& x) {
01227 typedef implementation<polynomial_evaluate,V> Pol;
01228 return Pol::tinterpolate (p, x);
01229 }
01230
01231 TMPL vector<Polynomial >
01232 expand (const Polynomial& p, const Vector& v, const vector<nat>& mu) {
01233 typedef implementation<polynomial_evaluate,V> Pol;
01234 nat k= N(v);
01235 ASSERT (N(mu) == k, "dimensions don't match");
01236 C** r= mmx_new<C*> (k);
01237 for (nat i=0; i<k; i++)
01238 r[i]= mmx_formatted_new<C> (aligned_size<C,V> (mu[i]), CF(p));
01239 Pol::expand (r, seg (p), seg (v), seg (mu), N (p), k);
01240 nat l= default_aligned_size<Polynomial > (k);
01241 Polynomial* ret= mmx_formatted_new<Polynomial > (l, get_format (p));
01242 for (nat i=0; i<k; i++)
01243 ret[i]= Polynomial (r[i], mu[i], aligned_size<C,V> (mu[i]), CF(p));
01244 mmx_delete<C*> (r, k);
01245 return vector<Polynomial > (ret, k, l);
01246 }
01247
01248
01249
01250
01251
01252 TMPL Polynomial
01253 range (const Polynomial& P, nat start, nat end) {
01254 typedef implementation<polynomial_linear,V> Pol;
01255 nat l= aligned_size<C,V> (end-start);
01256 C* r= mmx_formatted_new<C> (l, CF(P));
01257 for (nat i=start; i<end; i++) r[i-start]= P[i];
01258 return Polynomial (r, end-start, l, CF(P));
01259 }
01260
01261 TMPL Polynomial
01262 extract_mod (const Polynomial& P, nat k, nat p) {
01263 typedef implementation<polynomial_linear,V> Pol;
01264 nat n= (N(P) - k + p - 1) / p;
01265 nat l= aligned_size<C,V> (n);
01266 C* r= mmx_formatted_new<C> (l, CF(P));
01267 for (nat i=0; i<n; i++) r[i]= P[i*p+k];
01268 return Polynomial (r, n, l, CF(P));
01269 }
01270
01271 TMPL Polynomial
01272 lshiftz (const Polynomial& P, const int& shift) {
01273 typedef implementation<polynomial_linear,V> Pol;
01274 if (shift == 0) return P;
01275 else if (shift <= -((int) N(P))) return promote (0, P);
01276 else if (shift < 0) return range (P, (nat) (-shift), N(P));
01277 else {
01278 nat n= N(P) + shift;
01279 nat l= aligned_size<C,V> (n);
01280 C* r= mmx_formatted_new<C> (l, CF(P));
01281 Pol::copy (r+shift, seg (P), N(P));
01282 return Polynomial (r, n, l, CF(P));
01283 }
01284 }
01285
01286 TMPL Polynomial
01287 reverse (const Polynomial& P) {
01288 typedef implementation<polynomial_linear,V> Pol;
01289 typedef implementation<vector_linear,V> Vec;
01290 nat n= N(P);
01291 nat l= aligned_size<C,V> (n);
01292 C* r= mmx_formatted_new<C> (l, CF(P));
01293 Vec::vec_reverse (r, seg (P), n);
01294 return Polynomial (r, n, l, CF(P));
01295 }
01296
01297 TMPL Polynomial
01298 graeffe (const Polynomial& P) {
01299 typedef implementation<polynomial_graeffe,V> Pol;
01300 nat n= N(P);
01301 nat l= aligned_size<C,V> (n);
01302 C* r= mmx_formatted_new<C> (l, CF(P));
01303 Pol::graeffe (r, seg (P), n);
01304 return Polynomial (r, n, l, CF(P));
01305 }
01306
01307 TMPL Polynomial
01308 invert_lo (const Polynomial& P, nat m) {
01309 typedef implementation<polynomial_divide,V> Pol;
01310 nat n= N(P);
01311 nat l= aligned_size<C,V> (m);
01312 C* r= mmx_formatted_new<C> (l, CF(P));
01313 if (n >= m)
01314 Pol::invert_lo (r, seg (P), m);
01315 else {
01316 C* t= mmx_formatted_new<C> (l, CF(P));
01317 Pol::copy (t, seg (P), n);
01318 Pol::clear (t + n, m - n);
01319 Pol::invert_lo (r, t, m);
01320 }
01321 return Polynomial (r, m, l, CF(P));
01322 }
01323
01324 TMPL inline Polynomial
01325 invert_lo (const Polynomial& P) {
01326 return invert_lo (P, N(P));
01327 }
01328
01329 TMPL Polynomial
01330 invert_hi (const Polynomial& P) {
01331 typedef implementation<polynomial_divide,V> Pol;
01332 nat n= N(P);
01333 nat l= aligned_size<C,V> (n);
01334 C* r= mmx_formatted_new<C> (l, CF(P));
01335 Pol::invert_hi (r, seg (P), n);
01336 return Polynomial (r, n, l, CF(P));
01337 }
01338
01339
01340
01341
01342
01343 template<typename Op, typename C, typename V>
01344 polynomial<Unary_return_type(Op,C),V>
01345 unary_map (const polynomial<C,V>& p) {
01346 typedef implementation<vector_linear,V> Vec;
01347 typedef Unary_return_type(Op,C) T;
01348 nat n= N(p);
01349 nat l= aligned_size<T,V> (n);
01350 format<T> fm= unary_map<Op> (CF(p));
01351 T* r= mmx_formatted_new<T> (l, fm);
01352 Vec::template vec_unary<Op> (r, seg (p), n);
01353 return polynomial<T,V> (r, n, l, fm);
01354 }
01355
01356 template<typename Op, typename C, typename V, typename X>
01357 polynomial<Binary_return_type(Op,C,X),V>
01358 binary_map_scalar (const polynomial<C,V>& p, const X& x) {
01359 typedef implementation<vector_linear,V> Vec;
01360 typedef Binary_return_type(Op,C,X) T;
01361 nat n= N(p);
01362 nat l= aligned_size<T,V> (n);
01363 format<T> fm= binary_map_scalar<Op> (CF(p), x);
01364 T* r= mmx_formatted_new<T> (l, fm);
01365 Vec::template vec_binary_scalar<Op> (r, seg (p), x, n);
01366 return polynomial<T,V> (r, n, l, fm);
01367 }
01368
01369 template<typename Op, typename C, typename V> Unary_return_type(Op,C)
01370 big (const polynomial<C,V>& v) {
01371 typedef implementation<vector_linear,V> Vec;
01372 return Vec::template vec_unary_big<Op> (seg (v), N(v));
01373 }
01374
01375 TMPL Polynomial copy (const Polynomial& P) {
01376 return unary_map<id_op> (P); }
01377 TMPL Polynomial duplicate (const Polynomial& P) {
01378 return unary_map<duplicate_op> (P); }
01379
01380
01381
01382
01383
01384 TMPL
01385 struct lift_helper<Polynomial> {
01386 template<typename T,typename W> static inline void
01387 set_op (polynomial<T,W>& w, const Polynomial& v) {
01388 typedef implementation<vector_linear,V> Vec;
01389 format<T> fm= CF(w);
01390 nat n= N(v);
01391 nat l= aligned_size<T,W> (n);
01392 T* r= mmx_formatted_new<T> (l, fm);
01393 Vec::template vec_unary<lift_op> (r, seg (v), n);
01394 w= polynomial<T,W> (r, n, l, fm); }
01395 static inline Lifted_polynomial
01396 op (const Polynomial& v) {
01397 Lifted_polynomial w (lift (CF(v)));
01398 set_op (w, v);
01399 return w; }
01400 };
01401
01402 TMPL Reconstructed_polynomial reconstruct (const Polynomial& p) {
01403 return as<Reconstructed_polynomial> (unary_map<reconstruct_op> (p)); }
01404
01405 template<typename C, typename V, typename W> bool
01406 is_reconstructible (const Polynomial& p,
01407 polynomial<Reconstruct_type(C),W>& q) {
01408 vector<Reconstruct_type(C)> v (CF(q));
01409 if (!is_reconstructible (coefficients (p), v)) return false;
01410 q= polynomial<Reconstruct_type(C),W> (v);
01411 return true; }
01412
01413
01414
01415
01416
01417 template<typename S1, typename D, typename Fun> polynomial<D>
01418 map (const Fun& fun, const polynomial<S1>& p1, const format<D>& fm) {
01419 nat n= N(p1);
01420 nat l= default_aligned_size<D> (n);
01421 D* r= mmx_formatted_new<D> (l, fm);
01422 for (nat i=0; i<n; i++) r[i]= fun (p1[i]);
01423 return polynomial<D> (r, n, l, fm);
01424 }
01425
01426
01427
01428
01429
01430 TMPL inline bool is_finite (const Polynomial& p) {
01431 return big<and_is_finite_op> (p); }
01432 TMPL inline bool is_nan (const Polynomial& p) {
01433 return big<or_is_nan_op> (p); }
01434 TMPL inline bool is_infinite (const Polynomial& p) {
01435 return !is_nan (p) && big<or_is_infinite_op> (p); }
01436 TMPL inline bool is_fuzz (const Polynomial& p) {
01437 return !is_nan (p) && big<or_is_fuzz_op> (p); }
01438 TMPL inline bool is_reliable (const Polynomial& p) {
01439 return is_reliable (promote (0, CF(p))); }
01440
01441 TMPL inline Polynomial change_precision (const Polynomial& P, xnat p) {
01442 return binary_map_scalar<change_precision_op> (P, p); }
01443 TMPL inline xnat precision (const Polynomial& p) {
01444 return big<min_precision_op> (p); }
01445
01446 TMPL inline xint exponent (const Polynomial& p) {
01447 return big<max_exponent_op> (p); }
01448 TMPL inline double magnitude (const Polynomial& p) {
01449 return big<max_magnitude_op> (p); }
01450
01451
01452
01453
01454
01455 TMPL Abs_polynomial abs (const Polynomial& p) { return unary_map<abs_op> (p); }
01456 TMPL Real_polynomial Re (const Polynomial& p) { return unary_map<Re_op> (p); }
01457 TMPL Real_polynomial Im (const Polynomial& p) { return unary_map<Im_op> (p); }
01458 TMPL Polynomial conj (const Polynomial& p) { return unary_map<conj_op> (p); }
01459
01460 TMPL Center_polynomial center (const Polynomial& p) {
01461 return unary_map<center_op> (p); }
01462 TMPL Radius_polynomial radius (const Polynomial& p) {
01463 return unary_map<radius_op> (p); }
01464 TMPL Center_polynomial lower (const Polynomial& p) {
01465 return unary_map<lower_op> (p); }
01466 TMPL Center_polynomial upper (const Polynomial& p) {
01467 return unary_map<upper_op> (p); }
01468 TMPL inline Polynomial sharpen (const Polynomial& p) {
01469 return unary_map<sharpen_op> (p); }
01470 TMPLK inline Polynomial blur (const Polynomial& p, const K& x) {
01471 return binary_map_scalar<blur_op> (p, x); }
01472
01473 #undef TMPL_DEF
01474 #undef TMPL
01475 #undef TMPLK
01476 #undef Format
01477 #undef Vector
01478 #undef Polynomial
01479 #undef Polynomial_rep
01480 #undef Abs_polynomial
01481 #undef Real_polynomial
01482 #undef Center_polynomial
01483 #undef Radius_polynomial
01484 #undef Lifted_polynomial
01485 #undef Projected_polynomial
01486 #undef Reconstructed_polynomial
01487 }
01488 #endif // __MMX_POLYNOMIAL_HPP
