00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_SERIES_NAIVE__HPP
00014 #define __MMX_SERIES_NAIVE__HPP
00015 #include <algebramix/polynomial.hpp>
00016
00017 namespace mmx {
00018 #define TMPL template<typename C, typename V>
00019 #define Format format<C>
00020 #define Series series<C,V>
00021 #define Series_rep series_rep<C,V>
00022 #define Recursive_series_rep recursive_series_rep<C,V>
00023 #define Polynomial polynomial<C, typename series_polynomial_helper<C,V>::PV>
00024
00025 template<typename C>
00026 struct series_variant_helper;
00027
00028 template <typename C, typename V=typename series_variant_helper<C>::SV>
00029 class series_rep;
00030
00031 template <typename C, typename V=typename series_variant_helper<C>::SV>
00032 class series;
00033
00034 template <typename C, typename V=typename series_variant_helper<C>::SV>
00035 class recursive_series_rep;
00036
00037
00038
00039
00040
00041 struct series_naive {};
00042
00043 template<typename C>
00044 struct series_variant_helper {
00045 typedef series_naive SV;
00046 };
00047
00048 template<typename C, typename V>
00049 struct series_polynomial_helper {
00050 typedef typename Polynomial_variant(C) PV;
00051 };
00052
00053
00054
00055
00056
00057 struct series_defaults {};
00058
00059 template<typename V>
00060 struct implementation<series_defaults,V,series_naive> {
00061
00062 template<typename S>
00063 class global_variables {
00064
00065 static inline generic& dyn_name () {
00066 static generic name = "z";
00067 return name; }
00068
00069 static inline nat& dyn_output_order () {
00070 static nat n = 10;
00071 return n; }
00072
00073 static inline nat& dyn_cancel_order () {
00074 static nat n = 5;
00075 return n; }
00076
00077 static inline bool& dyn_formula_output () {
00078 static bool b = false;
00079 return b; }
00080
00081 public:
00082 static inline void set_variable_name (const generic& x) {
00083 dyn_name () = x; }
00084 static inline generic get_variable_name () {
00085 return dyn_name (); }
00086 static inline void set_output_order (const nat& x) {
00087 dyn_output_order () = x; }
00088 static inline nat get_output_order () {
00089 return dyn_output_order (); }
00090 static inline void set_cancel_order (const nat& x) {
00091 dyn_cancel_order () = x; }
00092 static inline nat get_cancel_order () {
00093 return dyn_cancel_order (); }
00094 static inline void set_formula_output (const bool& x) {
00095 dyn_formula_output () = x; }
00096 static inline bool get_formula_output () {
00097 return dyn_formula_output (); }
00098 };
00099 };
00100
00101
00102
00103
00104
00105 struct series_scalar_abstractions {};
00106
00107 template<typename U>
00108 struct implementation<series_scalar_abstractions,U,series_naive>:
00109 public implementation<series_defaults,U>
00110 {
00111
00112 template<typename Op,typename C,typename V,typename X>
00113 class binary_scalar_series_rep: public Series_rep {
00114 protected:
00115 Series f;
00116 X x;
00117 public:
00118 inline binary_scalar_series_rep (const Series& f2, const X& x2):
00119 Series_rep (CF(f2)), f(f2), x (x2) {}
00120 syntactic expression (const syntactic& z) const {
00121 return Op::op (flatten (f, z), flatten (x)); }
00122 virtual void Increase_order (nat l) {
00123 Series_rep::Increase_order (l);
00124 increase_order (f, l); }
00125 virtual C next () { return Op::op (f[this->n], x); }
00126 };
00127
00128 template<typename Op,typename C,typename V,typename X,typename Y>
00129 class ternary_scalar_series_rep: public Series_rep {
00130 protected:
00131 Series f;
00132 X x;
00133 Y y;
00134 public:
00135 inline ternary_scalar_series_rep (const Series& f2, const X& a, const Y& b):
00136 Series_rep (CF(f2)), f (f2), x (a), y (b) {}
00137 syntactic expression (const syntactic& z) const {
00138 return Op::op (flatten (f, z), flatten (x), flatten (y)); }
00139 virtual void Increase_order (nat l) {
00140 Series_rep::Increase_order (l);
00141 increase_order (f, l); }
00142 virtual C next () { return Op::op (f[this->n], x, y); }
00143 };
00144 };
00145
00146
00147
00148
00149
00150 struct series_map_as_abstractions {};
00151
00152 template<typename U>
00153 struct implementation<series_map_as_abstractions,U,series_naive>:
00154 public implementation<series_scalar_abstractions,U>
00155 {
00156
00157 template<typename Op, typename C, typename V, typename S, typename SV>
00158 class unary_map_as_series_rep: public Series_rep {
00159 protected:
00160 series<S,SV> f;
00161 public:
00162 inline unary_map_as_series_rep (const series<S,SV>& f2, const Format& fm):
00163 Series_rep (fm), f (f2) {}
00164 inline unary_map_as_series_rep (const series<S,SV>& f2):
00165 Series_rep (unary_map<Op> (CF(f2))), f (f2) {}
00166 syntactic expression (const syntactic& z) const {
00167 return Op::op (flatten (f, z)); }
00168 virtual void Increase_order (nat l) {
00169 Series_rep::Increase_order (l);
00170 increase_order (f, l); }
00171 virtual C next () {
00172 C r= get_sample (this->tfm ());
00173 Op::set_op (r, f[this->n]);
00174 return r; }
00175 };
00176
00177 template<typename Op, typename C, typename V, typename S, typename SV>
00178 static inline Series
00179 unary_map_as (const series<S,SV>& f, const Format& fm) {
00180 typedef unary_map_as_series_rep<Op,C,V,S,SV> Unary;
00181 return (Series_rep*) new Unary (f, fm);
00182 }
00183
00184 template<typename Op, typename C, typename V, typename S, typename SV>
00185 static inline Series
00186 unary_map_as (const series<S,SV>& f) {
00187 return unary_map_as<Op,C,V,S,SV> (f, unary_map<Op> (CF(f)));
00188 }
00189
00190 };
00191
00192
00193
00194
00195
00196 struct series_recursive_abstractions {};
00197
00198 template<typename U>
00199 struct implementation<series_recursive_abstractions,U,series_naive>:
00200 public implementation<series_map_as_abstractions,U>
00201 {
00202
00203 template<typename Op,typename C,
00204 typename V=typename series_variant_helper<C>::SV>
00205 class nullary_recursive_series_rep: public Recursive_series_rep {
00206 protected:
00207 C c;
00208 public:
00209 nullary_recursive_series_rep (const C& c2):
00210 Recursive_series_rep (get_format (c)), c (c2) {}
00211 syntactic expression (const syntactic& z) const {
00212 (void) z;
00213 return Op::op_init (flatten (c)); }
00214 Series initialize () {
00215 this->initial (0)= c;
00216 return Op::def (Series (this->me ())); }
00217 };
00218
00219 template<typename Op,typename C,typename V>
00220 class unary_recursive_series_rep: public Recursive_series_rep {
00221 protected:
00222 Series f;
00223 bool with_init;
00224 C c;
00225 public:
00226 unary_recursive_series_rep (const Series& f2):
00227 Recursive_series_rep (CF(f2)), f (f2), with_init (false) {}
00228 unary_recursive_series_rep (const Series& f2, const C& c2):
00229 Recursive_series_rep (CF(f2)), f (f2), with_init (true), c (c2) {}
00230 syntactic expression (const syntactic& z) const {
00231 return Op::op (flatten (f, z)); }
00232 virtual void Increase_order (nat l) {
00233 Recursive_series_rep::Increase_order (l);
00234 increase_order (f, l); }
00235 Series initialize () {
00236 if (with_init) this->initial (0)= c;
00237 else {
00238 ASSERT (Op::nr_init () <= 1, "wrong number of initial conditions");
00239 if (Op::nr_init () == 1)
00240 this->initial (0)= Op::op (f[0]);
00241 }
00242 return Op::def (Series (this->me ()), f); }
00243 };
00244
00245 template<typename Op,typename C,typename V,typename L>
00246 class unary_polynomial_recursive_series_rep: public Recursive_series_rep {
00247 protected:
00248 polynomial<L> P;
00249 bool with_init;
00250 C c;
00251
00252 public:
00253 unary_polynomial_recursive_series_rep (const polynomial<L>& P2):
00254 Recursive_series_rep (CF(P2)), P (P2), with_init (false) {}
00255 unary_polynomial_recursive_series_rep (const polynomial<L>& P2, const C& c2):
00256 Recursive_series_rep (CF(P2)), P (P2), with_init (true), c (c2) {}
00257 syntactic expression (const syntactic& z) const {
00258 return Op::op (flatten (P, z)); }
00259 virtual void Increase_order (nat l) {
00260 Recursive_series_rep::Increase_order (l);}
00261 Series initialize () {
00262 if (with_init) this->initial (0)= c;
00263 else {
00264 ASSERT (Op::nr_init () <= 1, "wrong number of initial conditions");
00265 if (Op::nr_init () == 1) {
00266 this->initial (0)= Op::template op<C> (P);
00267 }
00268 }
00269 return Op::def (Series (this->me ()), P); }
00270 };
00271
00272 template<typename Op,typename C,typename V>
00273 class binary_recursive_series_rep: public Recursive_series_rep {
00274 protected:
00275 Series f, g;
00276 public:
00277 binary_recursive_series_rep (const Series& f2, const Series& g2):
00278 Recursive_series_rep (CF(f2)), f(f2), g(g2) {}
00279 syntactic expression (const syntactic& z) const {
00280 return Op::op (flatten (f, z), flatten (g, z)); }
00281 virtual void Increase_order (nat l) {
00282 Recursive_series_rep::Increase_order (l);
00283 increase_order (f, l);
00284 increase_order (g, l); }
00285 Series initialize () {
00286 ASSERT (Op::nr_init () <= 1, "wrong number of initial conditions");
00287 if (Op::nr_init () == 1) this->initial (0)= Op::op (f[0], g[0]);
00288 return Op::def (Series (this->me ()), f, g);
00289 }
00290 };
00291
00292 template<typename Op,typename C,typename V,typename X>
00293 class binary_scalar_recursive_series_rep: public Recursive_series_rep {
00294 protected:
00295 Series f;
00296 bool with_init;
00297 X x;
00298 C c;
00299
00300 public:
00301 binary_scalar_recursive_series_rep (const Series& f2, const X& x2):
00302 Recursive_series_rep (CF(f2)), f (f2), with_init (false), x (x2) {}
00303 binary_scalar_recursive_series_rep (const Series& f2, const X& x2,
00304 const C& c2):
00305 Recursive_series_rep (CF(f2)), f (f2), with_init (true), x (x2), c (c2) {}
00306 syntactic expression (const syntactic& z) const {
00307 return Op::op (flatten (f, z), flatten (x)); }
00308 virtual void Increase_order (nat l) {
00309 Recursive_series_rep::Increase_order (l);
00310 increase_order (f, l); }
00311 Series initialize () {
00312 if (with_init) this->initial (0)= c;
00313 else {
00314 ASSERT (Op::nr_init () <= 1, "wrong number of initial conditions");
00315 if (Op::nr_init () == 1)
00316 this->initial (0)= Op::op (f[0], x);
00317 }
00318 return Op::def (Series (this->me ()), f, x); }
00319 };
00320
00321 };
00322
00323
00324
00325
00326
00327 struct series_abstractions {};
00328
00329 template<typename U>
00330 struct implementation<series_abstractions,U,series_naive>:
00331 public implementation<series_recursive_abstractions,U>
00332 {
00333
00334 template<typename Op,typename C,typename V>
00335 class unary_series_rep: public Series_rep {
00336 protected:
00337 Series f;
00338 public:
00339 inline unary_series_rep (const Series& f2):
00340 Series_rep (CF(f2)), f(f2) {}
00341 syntactic expression (const syntactic& z) const {
00342 return Op::op (flatten (f, z)); }
00343 virtual void Increase_order (nat l) {
00344 Series_rep::Increase_order (l);
00345 increase_order (f, l); }
00346 virtual C next () { return Op::op (f[this->n]); }
00347 };
00348
00349 template<typename Op,typename C,typename V>
00350 class binary_series_rep: public Series_rep {
00351 protected:
00352 Series f, g;
00353 public:
00354 inline binary_series_rep (const Series& f2, const Series& g2):
00355 Series_rep (CF(f2)), f(f2), g (g2) {}
00356 syntactic expression (const syntactic& z) const {
00357 return Op::op (flatten (f, z), flatten (g, z)); }
00358 virtual void Increase_order (nat l) {
00359 Series_rep::Increase_order (l);
00360 increase_order (f, l);
00361 increase_order (g, l); }
00362 virtual C next () { return Op::op (f[this->n], g[this->n]); }
00363 };
00364
00365 };
00366
00367
00368
00369
00370
00371 struct series_multiply {};
00372
00373 template<typename U>
00374 struct implementation<series_multiply,U,series_naive>:
00375 public implementation<series_abstractions,U>
00376 {
00377 typedef implementation<series_abstractions,U> Ser;
00378 TMPL
00379 class mul_series_rep: public Ser::template binary_series_rep<mul_op,C,V> {
00380 public:
00381 inline mul_series_rep (const Series& f, const Series& g):
00382 Ser::template binary_series_rep<mul_op,C,V> (f, g) {}
00383 C next () {
00384 nat i;
00385 C acc= this->f[0] * this->g[this->n];
00386 for (i=1; i<=this->n; i++)
00387 acc += this->f[i] * this->g[this->n-i];
00388 return acc; }
00389 };
00390
00391 TMPL
00392 class truncate_mul_series_rep:
00393 public Ser::template binary_series_rep<mul_op,C,V> {
00394 nat nf, ng;
00395 public:
00396 inline truncate_mul_series_rep (const Series& f, const Series& g,
00397 nat nf2, nat ng2):
00398 Ser::template binary_series_rep<mul_op,C,V> (f, g), nf (nf2), ng (ng2) {}
00399 C next () {
00400 nat i, n= this->n;
00401 C acc= this->zero ();
00402 for (i= n >= ng ? n - ng + 1 : 0; i < min (1 + n, nf); i++)
00403 acc += this->f[i] * this->g[n-i];
00404 return acc; }
00405 };
00406
00407 TMPL static inline Series
00408 ser_mul (const Series& f, const Series& g) {
00409 typedef mul_series_rep<C,V> Mul_rep;
00410 if (is_exact_zero (f) || is_exact_zero (g))
00411 return Series (CF(f));
00412 return (Series_rep*) new Mul_rep (f, g); }
00413
00414 TMPL static inline Series
00415 ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
00416 typedef truncate_mul_series_rep<C,V> Mul_rep;
00417 if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
00418 return Series (CF(f));
00419 return (Series_rep*) new Mul_rep (f, g, nf, ng); }
00420
00421 };
00422
00423
00424
00425
00426
00427 struct series_divide {};
00428
00429 template<typename U>
00430 struct implementation<series_divide,U,series_naive>:
00431 public implementation<series_multiply,U>
00432 {
00433
00434 TMPL static inline Series
00435 ser_rdiv_sc (const Series& f, const C& c) {
00436 typedef implementation<series_scalar_abstractions,U> Ser;
00437 typedef typename Ser::
00438 template binary_scalar_series_rep<rdiv_op,C,V,C> Div_sc_rep;
00439 if (is_exact_zero (f)) return Series (CF(f));
00440 return (Series_rep*) new Div_sc_rep (f, c); }
00441
00442 TMPL static inline Series
00443 ser_div (const Series& f, const Series& g) {
00444 typedef implementation<series_recursive_abstractions,U> Ser;
00445 typedef typename Ser::
00446 template binary_recursive_series_rep<div_op,C,V> Div_rep;
00447 if (is_exact_zero (f)) return Series (CF(f));
00448 return recursive (Series ((Series_rep*) new Div_rep (f, g))); }
00449
00450 TMPL static inline Series
00451 ser_rquo_sc (const Series& f, const C& c) {
00452 typedef implementation<series_scalar_abstractions,U> Ser;
00453 typedef typename Ser
00454 ::template binary_scalar_series_rep<rquo_op,C,V,C> Rquo_rep;
00455 if (is_exact_zero (f)) return Series (CF(f));
00456 return (Series_rep*) new Rquo_rep (f, c); }
00457
00458 TMPL static inline Series
00459 ser_quo (const Series& f, const Series& g) {
00460 typedef implementation<series_recursive_abstractions,U> Ser;
00461 typedef typename Ser
00462 ::template binary_recursive_series_rep<quo_op,C,V> Quo_rep;
00463 if (is_exact_zero (f)) return Series (CF(f));
00464 return recursive (Series ((Series_rep*) new Quo_rep (f, g))); }
00465
00466 TMPL static inline Series
00467 ser_rrem_sc (const Series& f, const C& c) {
00468 typedef implementation<series_scalar_abstractions,U> Ser;
00469 typedef typename Ser
00470 ::template binary_scalar_series_rep<rrem_op,C,V,C> Rrem_rep;
00471 if (is_exact_zero (f)) return Series (CF(f));
00472 return (Series_rep*) new Rrem_rep (f, c); }
00473
00474 };
00475
00476
00477
00478
00479
00480 struct series_compose {};
00481
00482 template<typename U>
00483 struct implementation<series_compose,U,series_naive>:
00484 public implementation<series_multiply,U>
00485 {
00486
00487 TMPL
00488 class compose_series_rep: public Series_rep {
00489 Series f;
00490 Series g;
00491 Series cumul;
00492 Series power;
00493 public:
00494 inline compose_series_rep (const Series& f2, const Series& g2):
00495 Series_rep (CF(f2)), f (f2), g (g2),
00496 cumul (CF(f2)), power (promote (1, CF(f2))) {}
00497 syntactic expression (const syntactic& z) const {
00498 return apply (GEN_COMPOSE, flatten (f, z), flatten (g, z)); }
00499 void Increase_order (nat l) {
00500 Series_rep::Increase_order (l);
00501 increase_order (f, l);
00502 increase_order (g, l);
00503 increase_order (cumul, l);
00504 increase_order (power, l); }
00505 C next () {
00506 cumul += f[this->n] * power;
00507 power *= g;
00508 return cumul[this->n];
00509 }
00510 };
00511
00512 TMPL static inline Series
00513 ser_compose (const Series& f, const Series& g) {
00514 typedef compose_series_rep<C,V> Compose_rep;
00515 if (is_exact_zero (f)) return Series (CF(f));
00516 ASSERT (g[0] == 0, "constant coefficient must be zero");
00517 return (Series_rep*) new Compose_rep (f, g); }
00518
00519 TMPL
00520 class reverse_series_rep: public Recursive_series_rep {
00521 protected:
00522 Series f;
00523 public:
00524 reverse_series_rep (const Series& f2):
00525 Recursive_series_rep (CF(f2)), f (f2) {}
00526 syntactic expression (const syntactic& z) const {
00527 return apply (GEN_REVERSE, flatten (f, z)); }
00528 virtual void Increase_order (nat l) {
00529 Recursive_series_rep::Increase_order (l);
00530 increase_order (f, l); }
00531 Series initialize () {
00532 ASSERT (f[0] == 0 && f[1] != 0, "invalid reversion");
00533 C a= f[1];
00534 Series z (this->one (), 1);
00535 Series s= square (rshiftz (this->me (), 1));
00536 Series c= compose (rshiftz (this->f, 2), this->me ());
00537 return (z - lshiftz (c * s, 2)) / a; }
00538 };
00539
00540 TMPL static inline Series
00541 ser_reverse (const Series& f) {
00542 typedef reverse_series_rep<C,V> Reverse_rep;
00543 Series_rep* rep= new Reverse_rep (f);
00544 return recursive (Series (rep)); }
00545
00546 };
00547
00548
00549
00550
00551
00552 struct ser_separable_root_op {
00553
00554
00555 private:
00556 template<typename S> static S
00557 binpow_no_tangent (const S& me, nat r) {
00558 typedef Scalar_type(S) C;
00559
00560 if (r <= 1) return S (0);
00561 if (r == 2) return lshiftz (square (rshiftz (me, 1)), 2);
00562
00563 nat h= r >> 1;
00564 C x= me[0];
00565 C x_h_1= binpow (x, h-1);
00566 S y= me - me[0];
00567 S z= rshiftz (binpow_no_tangent (me, h), 1);
00568 S t0= as<C> (h) * y;
00569 S t1= square (rshiftz (t0 * x_h_1, 1));
00570 S t2= x_h_1 * z * (t0 + x);
00571 S w= lshiftz (square (z) + t1, 2) + lshiftz (t2 + t2, 1);
00572 if (r & 1) {
00573 nat h2= h << 1;
00574 w= rshiftz (w, 1);
00575 w= lshiftz ((y + x) * w, 1)
00576 + lshiftz (as<C> (h2) * square (rshiftz (y * x_h_1, 1)) * x, 2);
00577 }
00578 return w; }
00579
00580 public:
00581 static generic name () { return "separable_root"; }
00582
00583 template<typename S> static inline S
00584 op (const S& x, nat r) { return separable_root (x, r); }
00585
00586 static inline syntactic
00587 op (const syntactic& x, const syntactic& r) {
00588 return apply ("separable_root", x, r); }
00589
00590 template<typename R, typename S> static inline void
00591 set_op (R& x, const S& y, nat r) { x= separable_root (y, r); }
00592 template<typename S, typename I> static inline S
00593 op_init (const S& x, nat r, const I& i) {
00594 return separable_root_init (x, r, i); }
00595
00596 static inline nat nr_init () { return 1; }
00597
00598 template<typename S> static inline S
00599 def (const S& me, const S& a, nat r) {
00600 typedef Scalar_type(S) C;
00601 C me0_r_1= binpow (me[0], r - 1);
00602 C me0_r= me[0] * me0_r_1;
00603 return (a - me0_r - binpow_no_tangent (me, r)) / (as<C> (r) * me0_r_1); }
00604
00605 };
00606
00607 struct series_separable_root {};
00608
00609 template<typename U>
00610 struct implementation<series_separable_root,U,series_naive> {
00611
00612 TMPL static inline Series
00613 sep_root (const Series& a, nat r) {
00614 typedef implementation<series_recursive_abstractions,U> Ser;
00615 typedef typename Ser::template binary_scalar_recursive_series_rep
00616 <ser_separable_root_op,C,V,nat> Root;
00617 ASSERT (C(r) != 0, "wrong argument");
00618 if (is_exact_zero (a)) return Series (CF(a));
00619 Series_rep* rep= new Root (a, r);
00620 return recursive (Series (rep)); }
00621
00622 TMPL static inline Series
00623 sep_root (const Series& a, nat r, const C& init) {
00624 typedef implementation<series_recursive_abstractions,U> Ser;
00625 typedef typename Ser::template binary_scalar_recursive_series_rep
00626 <ser_separable_root_op,C,V,nat> Root;
00627 ASSERT (C(r) != 0, "wrong argument");
00628 if (is_exact_zero (a)) return Series (CF(a));
00629 Series_rep* rep= new Root (a, r, init);
00630 return recursive (Series (rep)); }
00631
00632 };
00633
00634
00635
00636
00637
00638 struct series_pth_root_reg {};
00639 struct series_pth_root {};
00640
00641 template<typename U>
00642 struct implementation<series_pth_root,U,series_naive> {
00643
00644 TMPL static inline Series
00645 unsep_root (const Series& f) { ERROR ("not implemented"); }
00646
00647 };
00648
00649
00650
00651
00652
00653 struct ser_polynomial_regular_root_op {
00654
00655 public:
00656 static generic name () { return "polynomial_regular_root"; }
00657
00658 template<typename C> static inline C
00659 op (const polynomial<Lift_type (C)>& P) {
00660
00661 return polynomial_modular_root (as<polynomial<C> > (P)); }
00662
00663
00664
00665
00666
00667
00668 static inline syntactic
00669 op (const syntactic& P) {
00670 return apply ("polynomial_regular_root", P); }
00671
00672 template<typename R, typename L> static inline void
00673 set_op (R& x, const polynomial<L>& P) {
00674 x= polynomial_regular_root (P); }
00675
00676 template<typename C, typename V, typename L> static inline Series
00677 op_init (const polynomial<L>& P, const C& i) {
00678 return polynomial_regular_root_init (P, i); }
00679
00680 static inline nat nr_init () { return 1; }
00681
00682 template<typename S, typename L> static inline S
00683 def (const S& me, const polynomial<L>& P) {
00684 typedef Scalar_type (S) C;
00685 L y0 (* me[0]);
00686 polynomial<L> Q;
00687 polynomial<L> sqr (vec<L> (square (y0), - 2 * y0, 1));
00688 polynomial<L> R= rem (P, sqr, Q);
00689 S Py0 (R[0]);
00690 S dPy0 (R[1]);
00691 nat n = N(Q);
00692 vector<S> vs (S(0), n);
00693 for (nat i= 0; i<n; i++)
00694 vs[i] = S (Q[i]);
00695 polynomial<S> sQ (vs);
00696 return ((Py0 - me[0] * dPy0 + lshiftz (square (rshiftz (me))
00697 * eval (sQ, me), 2)) / (-dPy0));
00698 }
00699 };
00700
00701 struct series_polynomial_regular_root {};
00702
00703 template<typename U>
00704 struct implementation<series_polynomial_regular_root,U,series_naive> {
00705
00706 template<typename C, typename V, typename L> static inline Series
00707 pol_root (const polynomial<L>& P) {
00708 typedef implementation<series_recursive_abstractions,U> Ser;
00709 typedef typename Ser::template unary_polynomial_recursive_series_rep
00710 <ser_polynomial_regular_root_op,C,V,L> Root;
00711 ASSERT (N(P) > 1, "Polynomial must be non constant");
00712 Series_rep* rep= new Root (P);
00713 return recursive (Series (rep)); }
00714
00715 template<typename C, typename V, typename L> static inline Series
00716 pol_root (const polynomial<L>& P, const C& init) {
00717 typedef implementation<series_recursive_abstractions,U> Ser;
00718 typedef typename Ser::template unary_polynomial_recursive_series_rep
00719 <ser_polynomial_regular_root_op,C,V,L> Root;
00720 ASSERT (N(P) > 1, "Polynomial must be non constant");
00721 Series_rep* rep= new Root (P, init);
00722 return recursive (Series (rep)); }
00723
00724 };
00725
00726 #undef TMPL
00727 #undef Format
00728 #undef Series
00729 #undef Series_rep
00730 #undef Recursive_series_rep
00731 #undef Polynomial
00732 }
00733 #endif // __MMX_SERIES_NAIVE__HPP
