00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_SERIES_CARRY_NAIVE_HPP
00014 #define __MMX_SERIES_CARRY_NAIVE_HPP
00015 #include <algebramix/p_expansion.hpp>
00016 #include <algebramix/series_naive.hpp>
00017 #include <algebramix/root_modular.hpp>
00018 #include <algebramix/series_vector.hpp>
00019 #include <basix/literal.hpp>
00020 #include <basix/compound.hpp>
00021 #include <basix/table.hpp>
00022
00023 namespace mmx {
00024 #define Series series<M,V>
00025 #define Series_rep series_rep<M,V>
00026 #define Recursive_series_rep recursive_series_rep<M,V>
00027 #define C typename M::C
00028 #define Modulus typename M::modulus
00029 #define TMPL template<typename M, typename V>
00030 #define TMPLW template<typename M, typename V, typename W>
00031 #define Vector vector<M,W>
00032
00033
00034
00035
00036
00037 struct series_carry_naive {};
00038
00039 template<typename V,typename W>
00040 struct implementation<V,W,series_carry_naive>:
00041 public implementation<V,W,series_naive> {};
00042
00043 template<typename M>
00044 struct series_polynomial_helper<M,series_carry_naive> {
00045 typedef typename P_expansion_variant (M) PV;
00046 };
00047
00048 template<typename M>
00049 struct series_carry_variant_helper {
00050 typedef series_carry_naive SV;
00051 };
00052
00053 #define Series_carry_variant(C) series_carry_variant_helper<C>::SV
00054
00055
00056
00057
00058
00059 template<typename U>
00060 struct implementation<series_defaults,U,series_carry_naive> {
00061 typedef implementation<series_defaults,U,series_naive> Ser;
00062
00063 template<typename S>
00064 class global_variables :
00065 public Ser::template global_variables<S> {
00066
00067 static inline generic& dyn_name () {
00068 static generic name = "p";
00069 return name; }
00070
00071 public:
00072 static inline void set_variable_name (const generic& x) {
00073 dyn_name () = x; }
00074 static inline generic get_variable_name () {
00075 return dyn_name (); }
00076 };
00077 };
00078
00079
00080
00081
00082
00083 template<typename U>
00084 struct implementation<series_scalar_abstractions,U,series_carry_naive>:
00085 public implementation<series_defaults,U>
00086 {
00087
00088 template<typename Op,typename M,typename V,typename X>
00089 class binary_scalar_series_rep: public Series_rep {
00090 protected:
00091 Series f;
00092 M x;
00093 C c;
00094 public:
00095 inline binary_scalar_series_rep (const Series& f2, const X& x2):
00096 Series_rep (CF(f2)), f(f2), x (as<M> (x2)), c (0) {}
00097 syntactic expression (const syntactic& z) const {
00098 return Op::op (flatten (f, z), flatten (x)); }
00099 virtual void Increase_order (nat l) {
00100 Series_rep::Increase_order (l);
00101 increase_order (f, l); }
00102 virtual M next () {
00103 return Op::op_mod (f[this->n].rep, x.rep, M::get_modulus (), c); }
00104 };
00105
00106 };
00107
00108
00109
00110
00111
00112 template<typename U>
00113 struct implementation<series_abstractions,U,series_carry_naive>:
00114 public implementation<series_recursive_abstractions,U>
00115 {
00116
00117 template<typename Op,typename M,typename V>
00118 class unary_series_rep: public Series_rep {
00119 protected:
00120 const Series f;
00121 C carry;
00122 public:
00123 inline unary_series_rep (const Series& f2):
00124 Series_rep (CF(f2)), f(f2), carry (C(0)) {}
00125 syntactic expression (const syntactic& z) const {
00126 return Op::op (flatten (f, z)); }
00127 virtual void Increase_order (nat l) {
00128 Series_rep::Increase_order (l);
00129 increase_order (f, l); }
00130 virtual M next () {
00131 return Op::op_mod (f[this->n].rep, M::get_modulus (), carry); }
00132 };
00133
00134 template<typename Op,typename M,typename V>
00135 class binary_series_rep: public Series_rep {
00136 protected:
00137 Series f, g;
00138 C carry;
00139 public:
00140 inline binary_series_rep (const Series& f2, const Series& g2):
00141 Series_rep (CF(f2)), f(f2), g (g2), carry (0) {}
00142 syntactic expression (const syntactic& z) const {
00143 return Op::op (flatten (f, z), flatten (g, z)); }
00144 virtual void Increase_order (nat l) {
00145 Series_rep::Increase_order (l);
00146 increase_order (f, l);
00147 increase_order (g, l); }
00148 virtual M next () {
00149 return Op::op_mod (f[this->n].rep, g[this->n].rep,
00150 M::get_modulus (), carry); }
00151 };
00152
00153 };
00154
00155
00156
00157
00158
00159 template<typename U>
00160 struct implementation<series_multiply,U,series_carry_naive>:
00161 public implementation<series_abstractions,U>
00162 {
00163 typedef implementation<series_abstractions,U> Ser;
00164
00165 TMPL
00166 class mul_series_rep: public Ser::template binary_series_rep<mul_op,M,V> {
00167 typedef Lift_type (M) L;
00168 L carry;
00169 public:
00170 inline mul_series_rep (const Series& f, const Series& g):
00171 Ser::template binary_series_rep<mul_op,M,V> (f, g), carry (0) {}
00172 M next () {
00173 Modulus p= M::get_modulus ();
00174 L acc= carry;
00175 for (nat i= 0; i <= this->n; i++)
00176 acc += L ((this->f[i]).rep) * L((this->g[this->n-i]).rep);
00177 return rem (acc, * p, carry); }
00178 };
00179
00180 TMPL
00181 class truncate_mul_series_rep:
00182 public Ser::template binary_series_rep<mul_op,M,V> {
00183 nat nf, ng;
00184 typedef Lift_type (M) L;
00185 L carry;
00186 public:
00187 inline truncate_mul_series_rep (const Series& f, const Series& g,
00188 nat nf2, nat ng2):
00189 Ser::template binary_series_rep<mul_op,M,V> (f, g), nf (nf2), ng (ng2),
00190 carry (0) {}
00191 M next () {
00192 nat i, n= this->n;
00193 Modulus p= M::get_modulus ();
00194 L acc= carry;
00195 for (i= n >= ng ? n - ng + 1 : 0; i < min (1 + n, nf); i++)
00196 acc += L ((this->f[i]).rep) * L((this->g[n-i]).rep);
00197 return rem (acc, * p, carry); }
00198 };
00199
00200 TMPL static inline Series
00201 ser_mul (const Series& f, const Series& g) {
00202 typedef mul_series_rep<M,V> Mul_rep;
00203 if (is_exact_zero (f) || is_exact_zero (g))
00204 return Series (CF(f));
00205 return (Series_rep*) new Mul_rep (f, g); }
00206
00207 TMPL static inline Series
00208 ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
00209 typedef truncate_mul_series_rep<M,V> Mul_rep;
00210 if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
00211 return Series (CF(f));
00212 return (Series_rep*) new Mul_rep (f, g, nf, ng); }
00213
00214 };
00215
00216
00217
00218
00219
00220 template<typename V>
00221 struct series_carry_lift {};
00222
00223 template<typename U,typename V,typename W>
00224 struct implementation<V,U,series_carry_lift<W> >:
00225 public implementation<V,U,W> {};
00226
00227 template<typename U,typename W>
00228 struct implementation<series_multiply,U,series_carry_lift<W> >:
00229 public implementation<series_abstractions,U>
00230 {
00231 typedef implementation<series_abstractions,U> Ser;
00232
00233 template<typename M,typename V>
00234 class mul_series_rep: public Ser::template binary_series_rep<mul_op,M,V> {
00235 typedef C L;
00236 series<L> H;
00237 L c;
00238 public:
00239 inline mul_series_rep (const Series& f, const Series& g):
00240 Ser::template binary_series_rep<mul_op,M,V> (f, g), c(0) {
00241 H= as<series<L> > (f) * as<series<L> > (g); }
00242 M next () {
00243 c += H[this->n];
00244 return M (as<C> (rem (c, * M::get_modulus (), c)), true); }
00245 };
00246
00247 template<typename M,typename V>
00248 class truncate_mul_series_rep:
00249 public Ser::template binary_series_rep<mul_op,M,V> {
00250 typedef C L;
00251 series<L> H;
00252 L c;
00253 public:
00254 inline truncate_mul_series_rep (const Series& f, const Series& g,
00255 nat nf, nat ng):
00256 Ser::template binary_series_rep<mul_op,M,V> (f, g), c(0) {
00257 H= truncate_mul (as<series<L> > (f), as<series<L> > (g), nf, ng); }
00258 M next () {
00259 c += H[this->n];
00260 return M (as<C> (rem (c, * M::get_modulus (), c)), true); }
00261 };
00262
00263 template<typename M,typename V> static inline Series
00264 ser_mul (const Series& f, const Series& g) {
00265 typedef mul_series_rep<M,V> Mul_rep;
00266 if (is_exact_zero (f) || is_exact_zero (g))
00267 return Series (CF(f));
00268 return (Series_rep*) new Mul_rep (f, g); }
00269
00270 TMPL static inline Series
00271 ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
00272 typedef truncate_mul_series_rep<M,V> Mul_rep;
00273 if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
00274 return Series (CF(f));
00275 return (Series_rep*) new Mul_rep (f, g, nf, ng); }
00276
00277 };
00278
00279
00280
00281
00282
00283 template<typename U>
00284 struct implementation<series_divide,U,series_carry_naive>:
00285 public implementation<series_recursive_abstractions,U>
00286 {
00287
00288 template<typename M,typename V,typename X>
00289 class carry_mul_sc_series_rep: public Series_rep {
00290 protected:
00291 Series f;
00292 M x;
00293 public:
00294 inline carry_mul_sc_series_rep (const Series& f2, const X& x2):
00295 Series_rep (CF(f2)), f(f2), x (as<M> (x2)) {}
00296 syntactic expression (const syntactic& z) const {
00297 return apply (syntactic ("carry_mul"), flatten (f, z), flatten (x)); }
00298 virtual void Increase_order (nat l) {
00299 Series_rep::Increase_order (l);
00300 increase_order (f, l); }
00301 virtual M next () {
00302 if (this->n == 0) return M (0);
00303 C unused= 0;
00304 C ret= 0;
00305 mul_mod (unused, *f[this->n-1], *x, M::get_modulus (), ret);
00306 return M(ret); }
00307 };
00308
00309 template<typename M,typename V> static inline Series
00310 carry_mul_series (const Series& f, const M& c) {
00311 return (Series_rep*) new carry_mul_sc_series_rep<M,V,M> (f, c);
00312 }
00313
00314 template<typename M,typename V,typename X>
00315 class rdiv_sc_series_rep: public Recursive_series_rep {
00316 protected:
00317 Series f;
00318 M c;
00319 public:
00320 rdiv_sc_series_rep (const Series& f2, const X& c2):
00321 Recursive_series_rep (CF(f2)), f(f2), c(as<M> (c2)) {}
00322 syntactic expression (const syntactic& z) const {
00323 return flatten (f, z) / flatten (c); }
00324 virtual void Increase_order (nat l) {
00325 Recursive_series_rep::Increase_order (l);
00326 increase_order (f, l); }
00327 Series initialize () {
00328 typedef typename Series_variant(C) SV;
00329 ASSERT (c != 0, "division by zero");
00330 M u= invert (c);
00331 this->initial (0)= u * f[0];
00332 Series tmp= f - carry_mul_series (this -> me (), c);
00333 return as<Series> (u * as<series<M,SV> > (tmp));
00334 }
00335 };
00336
00337 TMPL static inline Series
00338 ser_rdiv_sc (const Series& f, const M& c) {
00339 typedef rdiv_sc_series_rep<M,V,M> Div_sc_rep;
00340 if (is_exact_zero (f)) return Series (CF(f));
00341 return recursive (Series ((Series_rep*) new Div_sc_rep (f, c))); }
00342
00343 template<typename M,typename V>
00344 class div_series_rep: public Recursive_series_rep {
00345 protected:
00346 Series f, g;
00347
00348 public:
00349 div_series_rep (const Series& f2, const Series& g2):
00350 Recursive_series_rep (CF(f2)), f(f2), g(g2) {}
00351 syntactic expression (const syntactic& z) const {
00352 return flatten (f, z) / flatten (g, z); }
00353 virtual void Increase_order (nat l) {
00354 Recursive_series_rep::Increase_order (l);
00355 increase_order (f, l);
00356 increase_order (g, l); }
00357 Series initialize () {
00358 typedef typename Series_variant(C) SV;
00359 ASSERT (g[0] != 0, "division by zero");
00360 M u= invert (g[0]);
00361 this->initial (0)= u * f[0];
00362 Series tmp= f - lshiftz (this -> me () * rshiftz (g))
00363 - carry_mul_series (this -> me (), g[0]);
00364 return as<Series> (u * as<series<M,SV> > (tmp)); }
00365 };
00366
00367 TMPL static inline Series
00368 ser_div (const Series& f, const Series& g) {
00369 typedef div_series_rep<M,V> Div_rep;
00370 if (is_exact_zero (f)) return Series (CF(f));
00371 return recursive (Series ((Series_rep*) new Div_rep (f, g))); }
00372
00373 TMPL static inline Series
00374 ser_rquo_sc (const Series& f, const M& c) {
00375 return ser_rdiv_sc (f, c); }
00376
00377 TMPL static inline Series
00378 ser_quo (const Series& f, const Series& g) {
00379 return ser_div (f, g); }
00380
00381 TMPL static inline Series
00382 ser_rrem_sc (const Series& f, const M& c) {
00383 return Series (CF(f)); }
00384
00385 };
00386
00387
00388
00389
00390
00391 struct ser_carry_separable_root_op {
00392
00393
00394 template<typename S, typename I> static inline S
00395 binpow_no_tangent (const S& me, const I& r) {
00396 typedef Scalar_type(S) M;
00397
00398
00399 if (r <= 1) return S (0);
00400 if (r == 2) return lshiftz (square (rshiftz (me, 1)), 2);
00401
00402 I h= r >> 1;
00403 S s_h (coefficients (as<default_p_expansion (M)> (integer (h))));
00404 M x= me[0];
00405 S x_h_1= binpow (S(x), h-1);
00406 S y= me - me[0];
00407 S z= rshiftz (binpow_no_tangent (me, h), 1);
00408 S t0= s_h * y;
00409 S t1= square (rshiftz (t0 * x_h_1, 1));
00410 S t2= x_h_1 * z * (t0 + x);
00411 S w= lshiftz (square (z) + t1, 2) + lshiftz (t2 + t2, 1);
00412 if ((r & 1) == 1) {
00413 I h2= h << 1;
00414 S s_h2 (coefficients (as<default_p_expansion (M)> (integer (h2))));
00415 w= lshiftz (me * rshiftz (w, 1), 1)
00416 + lshiftz (s_h2 * square (rshiftz (y * x_h_1, 1)) * x, 2);
00417 }
00418 return w; }
00419
00420 static generic name () { return "separable_root"; }
00421
00422 template<typename S> static inline S
00423 op (const S& x, nat r) { return separable_root (x, r); }
00424
00425 static inline syntactic
00426 op (const syntactic& x, const syntactic& r) {
00427 return apply ("separable_root", x, r); }
00428
00429 template<typename R, typename S> static inline void
00430 set_op (R& x, const S& y, nat r) { x= separable_root (y, r); }
00431
00432 template<typename S, typename I> static inline S
00433 op_init (const S& x, nat r, const I& i) {
00434 return separable_root_init (x, r, i); }
00435
00436 static inline nat nr_init () { return 1; }
00437
00438 template<typename S> static inline S
00439 def (const S& me, const S& a, nat r) {
00440 typedef Scalar_type(S) M;
00441 if (r == 2) {
00442 return ((a - S (me[0]) * me[0] - binpow_no_tangent (me, r))
00443 / M (r)) / me[0];
00444 }
00445 S ser_r (coefficients (as<default_p_expansion (M)> (integer (r))));
00446 S me_0_r_1= binpow (S (me[0]), r - 1);
00447 S me_0_r= me[0] * me_0_r_1;
00448 return (a - me_0_r - binpow_no_tangent (me, r))
00449 / (ser_r * me_0_r_1); }
00450
00451 };
00452
00453 template<typename U>
00454 struct implementation<series_separable_root,U,series_carry_naive> {
00455
00456 TMPL static inline Series
00457 sep_root (const Series& a, nat r) {
00458 ASSERT (M(r) != 0, "wrong argument");
00459 if (is_exact_zero (a)) return Series (CF(a));
00460 return binary_scalar_recursive_series
00461 <ser_carry_separable_root_op> (a, r); }
00462
00463 TMPL static inline Series
00464 sep_root (const Series& a, nat r, const M& init) {
00465 ASSERT (M(r) != 0, "wrong argument");
00466 if (is_exact_zero (a)) return Series (CF(a));
00467 return binary_scalar_recursive_series
00468 <ser_carry_separable_root_op> (a, r, init); }
00469
00470 };
00471
00472
00473
00474
00475
00476 TMPL inline Series
00477 integer_to_series_carry (integer a) {
00478 if (a == integer (0)) return Series ();
00479 if (a > integer (0))
00480 return Series (coefficients (as<default_p_expansion (M)> (a)));
00481 return -Series (coefficients (as<default_p_expansion (M)> (-a)));
00482 }
00483
00484
00485
00486
00487 struct series_slp_polynomial_regular_root {};
00488
00489 template<typename U>
00490 struct implementation<series_slp_polynomial_regular_root,U,series_naive> {
00491
00492 template<typename M, typename V, typename L>
00493 class slp_polynomial_regular_root_series_rep: public Recursive_series_rep {
00494 public:
00495
00496 generic f, x;
00497 M y0;
00498 Series rme, sq_rme;
00499 table<Series,generic> ht_ser;
00500 table<vector<L>,generic> ht_tau;
00501 table<L,generic> ht_eval;
00502
00503 static inline Series
00504 rec_add (const Series& ep, vector<L>& tp,
00505 const Series& eq, const vector<L>& tq) {
00506 tp= tp + tq;
00507 return ep + eq;
00508 }
00509
00510 static inline Series
00511 rec_minus (const Series& ep, vector<L>& tp) {
00512 tp= -tp;
00513 return -ep;
00514 }
00515
00516 static inline Series
00517 rec_minus (const Series& ep, vector<L>& tp,
00518 const Series& eq, const vector<L>& tq) {
00519 tp= tp - tq;
00520 return ep - eq;
00521 }
00522
00523 static Series
00524 rec_prod (const Series& ep, vector<L>& tp,
00525 const Series& eq, const vector<L>& tq,
00526 const Series& rme, const Series& sq_rme) {
00527 ASSERT (N(tp) == 2 && N (tq) == 2, "Wrong size of tau");
00528 const L& ap= tp[0], bp= tp[1];
00529 const L& aq= tq[0], bq= tq[1];
00530 Series c1= (integer_to_series_carry<M,V> (bp) * rshiftz (eq, 2)
00531 + integer_to_series_carry<M,V> (bq) * rshiftz (ep, 2));
00532 Series c2= (integer_to_series_carry<M,V> (ap) * rshiftz (eq, 2)
00533 + integer_to_series_carry<M,V> (aq) * rshiftz (ep, 2));
00534 tp= vec<L> (ap * aq, ap * bq + bp * aq);
00535 return (ep * eq + lshiftz (c1 * lshiftz (rme) + c2
00536 + integer_to_series_carry<M,V> (bp * bq)
00537 * sq_rme, 2));
00538 }
00539
00540 static inline Series
00541 rec_square (const Series& ep, vector<L>& tp,
00542 const Series& rme, const Series sq_rme) {
00543 ASSERT (N(tp) == 2, "Wrong size of tau");
00544 const L& ap= tp[0], bp= tp[1];
00545 Series taup = (integer_to_series_carry<M,V> (2 * ap) +
00546 integer_to_series_carry<M,V> (2 * bp) * lshiftz (rme));
00547 tp= vec<L> (square (ap), 2 * ap * bp);
00548 return (square (ep) + lshiftz (taup * rshiftz (ep, 2), 2) +
00549 lshiftz (integer_to_series_carry<M,V> (square (bp)) * sq_rme, 2));
00550 }
00551
00552
00553 static Series
00554 _eval2 (const generic& e, const Series& a) {
00555 if (is<literal> (e))
00556 return a;
00557 if (is<L> (e))
00558 return integer_to_series_carry<M,V> (as<L> (e));
00559 if (is<compound> (e)) {
00560 if (N(e) == 2) {
00561 if (e[0] == GEN_MINUS)
00562 return - _eval2 (e[1], a);
00563 if (e[0] == GEN_SQUARE)
00564 return square (_eval2 (e[1], a));
00565 }
00566 if (N(e) == 3) {
00567 if (e[0] == GEN_MINUS)
00568 return _eval2 (e[1], a) - _eval2 (e[2], a);
00569 if (e[0] == GEN_PLUS)
00570 return _eval2 (e[1], a) + _eval2 (e[2], a);
00571 if (e[0] == GEN_TIMES)
00572 return _eval2 (e[1], a) * _eval2 (e[2], a);
00573 }
00574 }
00575 ERROR ("bug, unexpected type with generic"); }
00576
00577 static L
00578 _eval (const generic& e, const L& a, table<L,generic>& ht) {
00579 if (is<literal> (e))
00580 return a;
00581 if (is<L> (e))
00582 return as<L> (e);
00583 if (contains (ht, e))
00584 return ht[e];
00585 L l;
00586 if (is<compound> (e)) {
00587 if (N(e) == 2) {
00588 if (e[0] == GEN_MINUS)
00589 l = - _eval (e[1], a, ht);
00590 if (e[0] == GEN_SQUARE)
00591 l = square (_eval (e[1], a, ht));
00592 }
00593 if (N(e) == 3) {
00594 if (e[0] == GEN_MINUS)
00595 l = _eval (e[1], a, ht) - _eval (e[2], a, ht);
00596 if (e[0] == GEN_PLUS)
00597 l = _eval (e[1], a, ht) + _eval (e[2], a, ht);
00598 if (e[0] == GEN_TIMES)
00599 l = _eval (e[1], a, ht) * _eval (e[2], a, ht);
00600 }
00601 ht[e] = l;
00602 return l; }
00603 ERROR ("bug, unexpected type with generic"); }
00604
00605 static generic
00606 _derive (const generic& e, const generic& x,
00607 table<generic, generic>& ht_derive) {
00608 if (e == x) return as<generic> (L (1));
00609 if (is<L> (e))
00610 return as<generic> (as<L> (0));
00611 if (contains (ht_derive, e)) {
00612 return ht_derive [e];
00613 }
00614 generic der;
00615 if (is<compound> (e)) {
00616 if (N(e) == 2) {
00617 if (e[0] == GEN_MINUS)
00618 der = gen (GEN_MINUS, _derive (e[1], x, ht_derive));
00619 else if (e[0] == GEN_SQUARE) {
00620 der = gen (GEN_TIMES, as<generic> (L (2)),
00621 gen (GEN_TIMES, _derive (e[1], x, ht_derive), e[1]));
00622 }
00623 }
00624 if (N(e) == 3) {
00625 if (e[0] == GEN_MINUS)
00626 der = gen (GEN_MINUS, _derive (e[1], x, ht_derive),
00627 _derive (e[2], x, ht_derive));
00628 else if (e[0] == GEN_PLUS)
00629 der = gen (GEN_PLUS, _derive (e[1], x, ht_derive),
00630 _derive (e[2], x, ht_derive));
00631 else if (e[0] == GEN_TIMES)
00632 der = gen (GEN_PLUS,
00633 gen (GEN_TIMES, _derive (e[1], x, ht_derive), e[2]),
00634 gen (GEN_TIMES, e[1], _derive (e[2], x, ht_derive)));
00635 }
00636 ht_derive [e] = der;
00637 return der;
00638 }
00639 ERROR ("bug, unexpected type with generic"); }
00640
00641 Series
00642 _eps (const generic& e, const generic& x, vector<L>& tau) {
00643 if (contains (ht_ser, e)) {
00644 tau = ht_tau[e];
00645 return ht_ser[e];
00646 }
00647 Series ser;
00648 if (is<L> (e)) {
00649 tau= vec<L> (as<L> (e), L (0));
00650 ser = Series ();
00651 ht_tau[e] = tau;
00652 ht_ser[e] = ser;
00653 return ser;
00654 }
00655 if (e == x) {
00656 tau= vec<L> (L (* y0), L(1));
00657 ser = Series ();
00658 ht_tau[e] = tau;
00659 ht_ser[e] = ser;
00660 return ser;
00661 }
00662 if (is<compound> (e)) {
00663 if (N(e) == 2) {
00664 Series s= _eps (e[1], x, tau);
00665 if (e[0] == GEN_MINUS)
00666 ser = rec_minus (s, tau);
00667 if (e[0] == GEN_SQUARE)
00668 ser = rec_square (s, tau, rme, sq_rme);
00669 }
00670 if (N(e) == 3) {
00671 vector<L> tau2 (L (0), 2);
00672 Series s1= _eps (e[1], x, tau), s2= _eps (e[2], x, tau2);
00673 if (e[0] == GEN_MINUS)
00674 ser = rec_minus (s1, tau, s2, tau2);
00675 if (e[0] == GEN_PLUS)
00676 ser = rec_add (s1, tau, s2, tau2);
00677 if (e[0] == GEN_TIMES)
00678 ser = rec_prod (s1, tau, s2, tau2, rme, sq_rme);
00679 }
00680 ht_tau[e] = tau;
00681 ht_ser[e] = ser;
00682 return ser; }
00683 ERROR ("bug, unexpected type with generic"); }
00684
00685 static void
00686 increase_order_generic (generic ff, nat l) {
00687 if (is<Series> (ff))
00688 increase_order (as<Series> (ff), l);
00689 if (is<compound> (ff)) {
00690 if (N(ff) == 2)
00691 increase_order_generic (ff[1], l);
00692 if (N(ff) == 3) {
00693 increase_order_generic (ff[1], l);
00694 increase_order_generic (ff[2], l);
00695 } } }
00696
00697 public:
00698 slp_polynomial_regular_root_series_rep (const generic& f2, const generic& x2,
00699 const M& y02):
00700 Recursive_series_rep (get_format (y02)), f (f2), x (x2), y0 (y02) {
00701 rme= rshiftz (this->me ());
00702 sq_rme= square (rme);
00703 }
00704 syntactic expression (const syntactic& z) const {
00705 return z; }
00706
00707 virtual void Increase_order (nat l) {
00708 Recursive_series_rep::Increase_order (l);
00709 increase_order_generic (f, l); }
00710 Series initialize () {
00711 ht_ser = table<Series,generic> ();
00712 ht_tau = table<vector<L>,generic> ();
00713 table<L,generic> ht_eval = table<L,generic> ();
00714 L Ly0 (* y0);
00715 this->initial (0) = y0;
00716 table<generic, generic> ht_derive =
00717 table<generic, generic> ();
00718 generic der = _derive (f, x, ht_derive);
00719 Series evder= integer_to_series_carry<M,V> ( _eval (der, Ly0, ht_eval));
00720 ht_eval = table<L,generic> ();
00721 L deriv = _eval (der, Ly0, ht_eval);
00722 ht_eval = table<L,generic> ();
00723 Series cst= integer_to_series_carry<M,V> (Ly0 * deriv
00724 - _eval (f, Ly0, ht_eval));
00725
00726
00727
00728
00729
00730 vector<L> tau (L (0), 2);
00731 return (cst - _eps (f, x, tau)) / evder;
00732 }
00733 };
00734
00735 template<typename M, typename V, typename L> static inline Series
00736 slp_pol_root (const generic& f, const generic& x, const M& y0) {
00737 return recursive (Series ((Series_rep*)
00738 new slp_polynomial_regular_root_series_rep<M,V,L>
00739 (f, x, y0)));
00740 }
00741
00742 };
00743
00744 template<typename M, typename V, typename L> static inline Series
00745 slp_polynomial_regular_root (const generic& f, const generic& x, const M& y0) {
00746 typedef implementation<series_slp_polynomial_regular_root,V> Ser;
00747 return Ser::template slp_pol_root<M,V,L> (f, x, y0);
00748 }
00749
00750
00751
00752
00753
00754
00755 struct ser_carry_polynomial_regular_root_op {
00756
00757 static generic name () { return "polynomial_regular_root"; }
00758
00759 template<typename M> static inline M
00760 op (const polynomial<Lift_type (M)>& P) {
00761 typedef Lift_type (M) L;
00762 L p (* M::get_modulus());
00763 vector<M> cP (M(0), N(P));
00764 for (nat i= 0; i < N(P); i++)
00765 cP[i]= M (rem (P[i], p));
00766 return polynomial_modular_root (polynomial<M> (cP)); }
00767
00768 template<typename M> static inline M
00769 op (const polynomial<M>& P, const M& init) {
00770 return init; }
00771
00772 template<typename M, typename V, typename L> static inline Series
00773 op (const polynomial<L>& P) {
00774 return polynomial_regular_root<M,V,L> (P); }
00775
00776 template<typename M, typename V, typename L> static inline Series
00777 op (const polynomial<L>& P, const M& init) {
00778 return polynomial_regular_root<M,V,L> (P,init); }
00779
00780 static inline syntactic
00781 op (const syntactic& P) {
00782 return apply ("polynomial_regular_root", P); }
00783
00784 template<typename R, typename L> static inline void
00785 set_op (R& x, const polynomial<L>& P) {
00786 x= polynomial_regular_root (P); }
00787
00788 template<typename M, typename V, typename L> static inline Series
00789 op_init (const polynomial<L>& P, const M& i) {
00790 return polynomial_regular_root_init<M,V,L> (P, i); }
00791
00792 static inline nat nr_init () { return 1; }
00793
00794 template<typename M, typename V, typename L> static inline Series
00795 def (const Series& me, const polynomial<L>& P) {
00796 L y0 (* me[0]);
00797 polynomial<L> Q;
00798 polynomial<L> sqr (vec<L> (square (y0), - 2 * y0, 1));
00799 polynomial<L> R= rem (P, sqr, Q);
00800 Series Py0= integer_to_series_carry<M, V> (R[0]);
00801 Series dPy0= integer_to_series_carry<M, V> (R[1]);
00802 nat n = N(Q);
00803 vector<Series> vs (Series(), n);
00804 for (nat i= 0; i<n; i++)
00805 vs[i] = integer_to_series_carry<M, V> (Q[i]);
00806 polynomial<Series> sQ (vs);
00807 return ((Py0 - me[0] * dPy0 + lshiftz (square (rshiftz (me))
00808 * eval (sQ, me), 2)) / (-dPy0));
00809 }
00810 };
00811
00812 template<typename U>
00813 struct implementation<series_polynomial_regular_root,U,series_carry_naive> {
00814
00815 template<typename M, typename V, typename L> static inline Series
00816 pol_root (const polynomial<L>& P) {
00817 ASSERT (N(P) > 1, "Polynomial must be non constant");
00818 return unary_polynomial_recursive_series
00819 <ser_carry_polynomial_regular_root_op,M,V,L> (P); }
00820
00821 template<typename M, typename V, typename L> static inline Series
00822 pol_root (const polynomial<L>& P, const M& init) {
00823 ASSERT (N(P) > 1, "Polynomial must be non constant");
00824 return unary_polynomial_recursive_series
00825 <ser_carry_polynomial_regular_root_op,M,V,L> (P, init); }
00826 };
00827
00828
00829 #undef Series
00830 #undef Series_rep
00831 #undef Recursive_series_rep
00832 #undef C
00833 #undef Modulus
00834 #undef TMPL
00835 #undef TMPLW
00836 #undef Vector
00837
00838 }
00839 #endif // __MMX_SERIES_CARRY_NAIVE_HPP
