00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_SERIES_CARRY_BLOCKS_HPP
00014 #define __MMX_SERIES_CARRY_BLOCKS_HPP
00015 #include <algebramix/series_carry.hpp>
00016
00017 namespace mmx {
00018 #define Series series<M,V>
00019 #define Series_rep series_rep<M,V>
00020 #define C typename M::C
00021 #define Modulus typename M::modulus
00022 #define TMPL template<typename M,typename V>
00023 #define Monoblock_transformer series_carry_monoblock_transformer<M,V,s,BV>
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 template<typename V,
00034 nat s,
00035 typename BV,
00036 nat t=s>
00037 struct series_carry_monoblock {};
00038
00039 template<typename V,typename BV,typename U,typename W,nat s,nat t>
00040 struct implementation<U,W,series_carry_monoblock<V,s,BV,t> >:
00041 public implementation<U,W,V> {};
00042
00043
00044
00045
00046
00047 template<typename V,
00048 nat s,
00049 typename BV,
00050 nat t=s>
00051 struct series_carry_blocks {};
00052
00053 template<typename V,typename BV,typename U,typename W,nat s,nat t>
00054 struct implementation<U,W,series_carry_blocks<V,s,BV,t> >:
00055 public implementation<U,W,V> {};
00056
00057
00058
00059
00060
00061 #define Monoblock_modular modular<modulus<Lift_type(M)>, \
00062 modular_global_series_carry_monoblock \
00063 <M,s,BV> >
00064
00065 #define Monoblock_series_rep series_rep<Monoblock_modular,BV>
00066 #define Monoblock_series series <Monoblock_modular,BV>
00067
00068 template<typename M,nat s,typename BV>
00069 struct modular_global_series_carry_monoblock {
00070 template<typename P>
00071 class modulus_storage {
00072 static inline P& dyn_modulus () {
00073 static P modulus = P ();
00074 return modulus; }
00075 public:
00076 static inline void set_modulus (const P& p) { dyn_modulus () = p; }
00077 static inline P get_modulus () { return dyn_modulus (); }
00078 };
00079 };
00080
00081 template<typename M, typename V, nat s, typename BV>
00082 class series_carry_monoblock_transformer {
00083
00084 class to_monoblock_series_rep: public Monoblock_series_rep {
00085 public:
00086 typedef Lift_type(M) L;
00087 typedef typename Base_transformer_unsigned(L,C) Baser;
00088
00089 protected:
00090 Series f;
00091 C* buf;
00092 Baser* baser;
00093
00094 public:
00095 to_monoblock_series_rep (const Series& f2):
00096 Monoblock_series_rep (format<Monoblock_modular > ()), f(f2) {
00097 Monoblock_modular::set_modulus (binpow (L(* M::get_modulus ()), s));
00098 buf= mmx_new<C> (s);
00099 baser= mmx_new<Baser> (1, * M::get_modulus ()); }
00100
00101 ~to_monoblock_series_rep () {
00102 mmx_delete<Baser> (baser, 1);
00103 mmx_delete<C> (buf, s); }
00104
00105 syntactic expression (const syntactic& z) const {
00106 return apply ("to_monoblock", flatten (f, z), flatten (s)); }
00107
00108 void Increase_order (nat l) {
00109 Monoblock_series_rep::Increase_order (l);
00110 increase_order (f, l * s); }
00111
00112 Monoblock_modular next () {
00113 nat start = s * this->n;
00114 for (nat i= 0; i < s; i++) buf[i]= * f[start + i];
00115 L e; inverse_base (e, buf, s, * baser);
00116 return Monoblock_modular (e, true); }
00117 };
00118
00119 class from_monoblock_series_rep: public Series_rep {
00120 public:
00121 typedef Lift_type(M) L;
00122 typedef typename Base_transformer_unsigned(L,C) Baser;
00123
00124 protected:
00125 Monoblock_series F;
00126 C* buf;
00127 Baser* baser;
00128 public:
00129
00130 from_monoblock_series_rep (const Monoblock_series& F2):
00131 Series_rep (format<M> ()), F(F2) {
00132 buf= mmx_new<C> (s);
00133 baser= mmx_new<Baser> (1, * M::get_modulus ()); }
00134
00135 ~from_monoblock_series_rep () {
00136 mmx_delete<C> (buf, s);
00137 mmx_delete<Baser> (baser, 1); }
00138
00139 syntactic expression (const syntactic& z) const {
00140 return apply ("from_monoblock", flatten (F, z), flatten (s)); }
00141
00142 void Increase_order (nat l) {
00143 Series_rep::Increase_order (l);
00144 increase_order (F, (l + s - 1) / s); }
00145
00146 M next () {
00147 nat q= this->n / s, r= this->n % s;
00148 if (r == 0) {
00149 nat m= direct_base (buf, s, * F[q], * baser);
00150 for (nat i= m; i < s; i++) buf[i]= 0;
00151 }
00152 return buf[r];
00153 }
00154 };
00155
00156 public:
00157
00158 inline series_carry_monoblock_transformer () {}
00159
00160 inline Monoblock_series
00161 to_monoblock (const Series& f) const {
00162 if (is_exact_zero (f))
00163 return Monoblock_series (Monoblock_modular (0) );
00164 return (Monoblock_series_rep*) new to_monoblock_series_rep (f);
00165 }
00166
00167 inline Series
00168 from_monoblock (const Monoblock_series& f) const {
00169 if (is_exact_zero (f))
00170 return Series (C(0) );
00171 return (Series_rep*) new from_monoblock_series_rep (f);
00172 }
00173
00174 };
00175
00176
00177
00178 template<typename M, typename V, nat s, typename BV> Monoblock_series
00179 to_monoblock (const Series& f,
00180 const Monoblock_transformer& blocker) {
00181 return blocker.to_monoblock (f);
00182 }
00183
00184 template<typename M, typename V, nat s, typename BV> Series
00185 from_monoblock (const Monoblock_series& f,
00186 const Monoblock_transformer& blocker) {
00187 return blocker.from_monoblock (f);
00188 }
00189
00190
00191
00192
00193
00194 template<typename Op, typename M, typename V, nat s, typename BV, nat t>
00195 class binary_monoblock_series_rep: public Series_rep {
00196 Monoblock_transformer blocker;
00197 Series f, g;
00198 bool h1_init;
00199 Series h0, h1;
00200
00201 public:
00202 binary_monoblock_series_rep (const Series& f2, const Series& g2):
00203 Series_rep (CF(f2)), f(f2), g(g2), h1_init (false) {
00204 h0= Op::op (f, g); }
00205 syntactic expression (const syntactic& z) const {
00206 return Op::op (flatten (f, z), flatten (g, z)); }
00207 virtual void Increase_order (nat l) {
00208 Series_rep::Increase_order (l);
00209 if (l <= t)
00210 increase_order (h0, l);
00211 else {
00212 if (! h1_init) {
00213 h1= from_monoblock (Op::op (to_monoblock (f, blocker),
00214 to_monoblock (g, blocker)),
00215 blocker);
00216 h1_init= true;
00217 }
00218 increase_order (h1, l); } }
00219 M next () {
00220 return this->n < t ? h0[this->n] : h1[this->n]; }
00221 };
00222
00223 template<typename Op, typename M, typename V, nat s, typename BV, nat t>
00224 inline Series
00225 binary_monoblock_series (const Series& f, const Series& g) {
00226 typedef binary_monoblock_series_rep<Op,M,V,s,BV,t> Op_rep;
00227 return (Series_rep*) new Op_rep (f, g); }
00228
00229
00230
00231
00232
00233 template<typename M, typename V, nat s, typename BV, nat t>
00234 class truncate_mul_monoblock_series_rep: public Series_rep {
00235 Monoblock_transformer blocker;
00236 Series f, g;
00237 nat nf, ng;
00238 bool h1_init;
00239 Series h0, h1;
00240
00241 public:
00242 truncate_mul_monoblock_series_rep (const Series& f2, const Series& g2,
00243 nat nf2, nat ng2):
00244 Series_rep (CF(f2)), f(f2), g(g2), nf (nf2), ng (ng2), h1_init (false) {
00245 h0= truncate_mul (f, g, nf, ng); }
00246 syntactic expression (const syntactic& z) const {
00247 return apply ("truncate_mul", flatten (f, z), flatten (g, z),
00248 flatten (nf), flatten (ng)); }
00249 virtual void Increase_order (nat l) {
00250 Series_rep::Increase_order (l);
00251 if (l <= t)
00252 increase_order (h0, l);
00253 else {
00254 if (! h1_init) {
00255 Series f0= piecewise (f, Series (CF(f)), nf);
00256 Series g0= piecewise (g, Series (CF(g)), ng);
00257 h1= from_monoblock (truncate_mul (to_monoblock (f0, blocker),
00258 to_monoblock (g0, blocker),
00259 (nf + s - 1) / s,
00260 (ng + s - 1) / s),
00261 blocker);
00262 h1_init= true;
00263 }
00264 increase_order (h1, l); } }
00265 M next () {
00266 return this->n < t ? h0[this->n] : h1[this->n]; }
00267 };
00268
00269 template<typename M, typename V, nat s, typename BV, nat t>
00270 inline Series
00271 truncate_mul_monoblock_series (const Series& f, const Series& g,
00272 nat nf, nat ng) {
00273 typedef truncate_mul_monoblock_series_rep<M,V,s,BV,t> Mul_rep;
00274 return (Series_rep*) new Mul_rep (f, g, nf, ng); }
00275
00276
00277
00278
00279
00280 template<typename Op, typename M, typename V, nat s,
00281 typename BV, nat t, typename L>
00282 class unary_polynomial_recursive_monoblock_series_rep: public Series_rep {
00283 Monoblock_transformer blocker;
00284 polynomial<L> P;
00285 bool h1_init;
00286 Series h0, h1;
00287
00288 public:
00289 unary_polynomial_recursive_monoblock_series_rep (const polynomial<L>& P2):
00290 Series_rep (CF(P2)), P(P2), h1_init (false) {
00291 h0= Op::template op<M,V,L> (P); }
00292 unary_polynomial_recursive_monoblock_series_rep (const polynomial<L>& P2,
00293 const M& init):
00294 Series_rep (CF(P2)), P(P2), h1_init (false) {
00295 h0= Op::op (P, init); }
00296 syntactic expression (const syntactic& z) const {
00297 return apply (Op::name (), flatten (P, z)); }
00298 virtual void Increase_order (nat l) {
00299 Series_rep::Increase_order (l);
00300 if (l <= t)
00301 increase_order (h0, l);
00302 else {
00303 if (! h1_init) {
00304 Monoblock_modular::set_modulus (binpow (L(* M::get_modulus ()), s));
00305 Monoblock_modular init= as<Lift_type(M)> (truncate (h0, s));
00306 h1= from_monoblock (Op::template op_init<Monoblock_modular,BV,L>
00307 (P, init), blocker);
00308 h1_init= true;
00309 }
00310 increase_order (h1, l); } }
00311 M next () {
00312 return this->n < t ? h0[this->n] : h1[this->n]; }
00313 };
00314
00315 template<typename Op, typename M, typename V, nat s,
00316 typename BV, nat t, typename L> inline Series
00317 unary_polynomial_recursive_monoblock_series (const polynomial<L>& P) {
00318 typedef unary_polynomial_recursive_monoblock_series_rep<Op,M,V,s,BV,t,L>
00319 Op_rep;
00320 return (Series_rep*) new Op_rep (P); }
00321
00322 template<typename Op, typename M, typename V, nat s,
00323 typename BV, nat t, typename L> inline Series
00324 unary_polynomial_recursive_monoblock_series (const polynomial<L>& P,
00325 const M& init) {
00326 typedef unary_polynomial_recursive_monoblock_series_rep<Op,M,V,s,BV,t,L>
00327 Op_rep;
00328 return (Series_rep*) new Op_rep (P, init); }
00329
00330
00331 template<typename Op, typename M, typename V, nat s,
00332 typename BV, nat t, typename X>
00333 class binary_scalar_recursive_monoblock_series_rep: public Series_rep {
00334 Monoblock_transformer blocker;
00335 Series f;
00336 X x;
00337 bool h1_init;
00338 Series h0, h1;
00339
00340 public:
00341 binary_scalar_recursive_monoblock_series_rep (const Series& f2, const X& x2):
00342 Series_rep (CF(f2)), f(f2), x(x2), h1_init (false) {
00343 h0= Op::op (f, x); }
00344 binary_scalar_recursive_monoblock_series_rep (const Series& f2, const X& x2,
00345 const M& init):
00346 Series_rep (CF(f2)), f(f2), x(x2), h1_init (false) {
00347 h0= Op::op (f, x, init); }
00348 syntactic expression (const syntactic& z) const {
00349 return apply (Op::name (), flatten (f, z), flatten (x)); }
00350 virtual void Increase_order (nat l) {
00351 Series_rep::Increase_order (l);
00352 if (l <= t)
00353 increase_order (h0, l);
00354 else {
00355 if (! h1_init) {
00356 Monoblock_modular init= as<Lift_type(M)> (truncate (h0, s));
00357 h1= from_monoblock (Op::op_init (to_monoblock (f, blocker), x, init),
00358 blocker);
00359 h1_init= true;
00360 }
00361 increase_order (h1, l); } }
00362 M next () {
00363 return this->n < t ? h0[this->n] : h1[this->n]; }
00364 };
00365
00366 template<typename Op, typename M, typename V, nat s,
00367 typename BV, nat t, typename X> inline Series
00368 binary_scalar_recursive_monoblock_series (const Series& f, const X& x) {
00369 typedef binary_scalar_recursive_monoblock_series_rep<Op,M,V,s,BV,t,X> Op_rep;
00370 return (Series_rep*) new Op_rep (f, x); }
00371
00372 template<typename Op, typename M, typename V, nat s,
00373 typename BV, nat t, typename X> inline Series
00374 binary_scalar_recursive_monoblock_series (const Series& f, const X& x,
00375 const M& init) {
00376 typedef binary_scalar_recursive_monoblock_series_rep<Op,M,V,s,BV,t,X> Op_rep;
00377 return (Series_rep*) new Op_rep (f, x, init); }
00378
00379
00380 #undef Monoblock_series_rep
00381 #undef Monoblock_series
00382
00383
00384
00385
00386
00387 template<typename U,typename W,nat s,typename BV,nat t>
00388 struct implementation<series_multiply,U,series_carry_monoblock<W,s,BV,t> >:
00389 public implementation<series_abstractions,U>
00390 {
00391
00392 TMPL static inline Series
00393 ser_mul (const Series& f, const Series& g) {
00394 if (is_exact_zero (f) || is_exact_zero (g))
00395 return Series (CF(f));
00396 return as<Series> (binary_monoblock_series<mul_op,M,W,s,BV,t>
00397 (as<series<M,W> > (f), as<series<M,W> > (g))); }
00398
00399 TMPL static inline Series
00400 ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
00401 if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
00402 return Series (CF(f));
00403 return as<Series> (binary_monoblock_series<mul_op,M,W,s,BV,t>
00404 (as<series<M,W> > (f), as<series<M,W> > (g))); }
00405 };
00406
00407
00408
00409
00410
00411 template<typename U,typename W,nat s,typename BV,nat t>
00412 struct implementation<series_divide,U,series_carry_monoblock<W,s,BV,t> >:
00413 public implementation<series_multiply,U>
00414 {
00415
00416 TMPL static inline Series
00417 ser_rdiv_sc (const Series& f, const M& c) {
00418 typedef implementation<series_divide,W> Ser;
00419 return Ser::ser_rdiv_sc (f, c); }
00420
00421 TMPL static inline Series
00422 ser_div (const Series& f, const Series& g) {
00423 if (is_exact_zero (f))
00424 return Series (CF(f));
00425 return as<Series> (binary_monoblock_series<div_op,M,W,s,BV,t>
00426 (as<series<M,W> > (f), as<series<M,W> > (g))); }
00427
00428 TMPL static inline Series
00429 ser_rquo_sc (const Series& f, const M& c) {
00430 return ser_rdiv_sc (f, c); }
00431
00432 TMPL static inline Series
00433 ser_quo (const Series& f, const Series& g) {
00434 return ser_div (f, g); }
00435
00436 TMPL static inline Series
00437 ser_rrem_sc (const Series& f, const M& c) {
00438 return Series (CF(f)); }
00439
00440 };
00441
00442
00443
00444
00445
00446 template<typename U,typename W,nat s,typename BV,nat t>
00447 struct implementation<series_separable_root,U,
00448 series_carry_monoblock<W,s,BV,t> > {
00449
00450 TMPL static inline Series
00451 sep_root (const Series& f, nat r) {
00452 ASSERT (M(r) != 0, "not implemented");
00453 return as<Series>
00454 (binary_scalar_recursive_monoblock_series
00455 <ser_separable_root_op,M,W,s,BV,t,nat> (as<series<M,W> > (f), r)); }
00456
00457 TMPL static inline Series
00458 sep_root_init (const Series& f, nat r, const M& init) {
00459 ASSERT (M(r) != 0, "not implemented");
00460 return as<Series>
00461 (binary_scalar_recursive_monoblock_series
00462 <ser_separable_root_op,M,W,s,BV,t,nat> (as<series<M,W> > (f), r, init)); }
00463 };
00464
00465
00466
00467
00468 template<typename U,typename W,nat s,typename BV,nat t>
00469 struct implementation<series_slp_polynomial_regular_root,U,
00470 series_carry_monoblock<W,s,BV,t> > {
00471
00472 template<typename M, typename V, typename L>
00473 class slp_polynomial_regular_root_monoblock_series_rep: public Series_rep {
00474 protected:
00475
00476 generic f, x;
00477 M y0;
00478
00479 Monoblock_transformer blocker;
00480 bool h1_init;
00481 Series h0;
00482 series<Monoblock_modular,BV> h1;
00483
00484 public:
00485 slp_polynomial_regular_root_monoblock_series_rep (const generic& f2,
00486 const generic& x2,
00487 const M& y02):
00488 Series_rep (get_format (y02)), f(f2), x(x2), y0(y02), h1_init(false) {
00489 h0 = slp_polynomial_regular_root<M,V,L> (f, x, y0);
00490 }
00491 syntactic expression (const syntactic& z) const {
00492 return z;
00493 }
00494 virtual void Increase_order (nat l) {
00495 Series_rep::Increase_order (l);
00496 if (l <= t)
00497 increase_order (h0, l);
00498 else {
00499 if (! h1_init) {
00500 Monoblock_modular::set_modulus (binpow (L(* M::get_modulus ()), s));
00501 Monoblock_modular init= as<Lift_type(M)> (truncate (h0, s));
00502 h1= from_monoblock (slp_polynomial_regular_root<Monoblock_modular,BV,L>
00503 (f, x, init), blocker);
00504 h1_init= true;
00505 }
00506 increase_order (h1, l); } }
00507 M next () {
00508 return this->n < t ? h0[this->n] : M (* h1[this->n]); }
00509 };
00510
00511 template<typename M, typename V, typename L> static inline Series
00512 slp_pol_root (const generic& f, const generic& x, const M& y0) {
00513 typedef slp_polynomial_regular_root_monoblock_series_rep<M,W,L> Pol_rep;
00514 return as<series<M,W> > ((series_rep<M,W>*) new Pol_rep (f, x, y0));
00515 }
00516
00517 };
00518
00519
00520
00521
00522
00523
00524 template<typename U,typename W,nat s,typename BV,nat t>
00525 struct implementation<series_polynomial_regular_root,U,
00526 series_carry_monoblock<W,s,BV,t> >:
00527 public implementation<series_multiply,U>
00528 {
00529
00530 template<typename M, typename V, typename L> static inline Series
00531 pol_root (const polynomial<L>& P) {
00532 ASSERT (N(P) > 1, "Polynomial must be non constant");
00533 return as<Series>
00534 (unary_polynomial_recursive_monoblock_series
00535 <ser_carry_polynomial_regular_root_op,M,W,s,BV,t,L> (P)); }
00536
00537 template<typename M, typename V, typename L> static inline Series
00538 pol_root_init (const polynomial<L>& P, const M& init) {
00539 ASSERT (N(P) > 1, "Polynomial must be non constant");
00540 return as<Series>
00541 (unary_polynomial_recursive_monoblock_series
00542 <ser_carry_polynomial_regular_root_op,M,W,s,BV,t,L> (P, init)); }
00543
00544 };
00545
00546
00547
00548
00549
00550
00551 template<typename U,typename W,nat s,typename BV,nat t>
00552 struct implementation<series_multiply,U,series_carry_blocks<W,s,BV,t> >:
00553 public implementation<series_abstractions,U>
00554 {
00555
00556 TMPL
00557 class mul_series_rep: public Series_rep {
00558 series_carry_monoblock_transformer<M,W,s,BV> blocker;
00559 Series f, g, h;
00560
00561 public:
00562 mul_series_rep (const Series& f2, const Series& g2):
00563 Series_rep (CF(f2)), f(f2), g(g2)
00564 {
00565 VERIFY (t >= s-1, "wrong threshold");
00566 series<M,W> f0 = as<series<M,W> > (f);
00567 series<M,W> g0 = as<series<M,W> > (g);
00568 series<M,W> f1 = rshiftz (f0, (int) t);
00569 series<M,W> g1 = rshiftz (g0, (int) t);
00570 h= as<Series> (truncate_mul (f0, g0, t, t)
00571 + lshiftz (truncate_mul (f0, g1, t, Maximal (nat)) +
00572 truncate_mul (f1, g0, Maximal (nat), t),
00573 (int) t)
00574 + lshiftz (from_monoblock (to_monoblock (f1, blocker) *
00575 to_monoblock (g1, blocker),
00576 blocker),((int) t) << 1)); }
00577 syntactic expression (const syntactic& z) const {
00578 return flatten (f, z) * flatten (g, z); }
00579 virtual void Increase_order (nat l) {
00580 Series_rep::Increase_order (l);
00581 increase_order (h, l); }
00582 M next () {
00583 return h[this->n]; }
00584 };
00585
00586 TMPL static inline Series
00587 ser_mul (const Series& f, const Series& g) {
00588 if (is_exact_zero (f) || is_exact_zero (g))
00589 return Series (CF(f));
00590 return (Series_rep*) new mul_series_rep<M,V> (f, g); }
00591
00592 TMPL static inline Series
00593 ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
00594 typedef mul_series_rep<M,V> Mul_rep;
00595 if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
00596 return Series (CF(f));
00597 return (Series_rep*)
00598 new Mul_rep (piecewise (f, Series (CF(f)), nf),
00599 piecewise (g, Series (CF(g)), ng)); }
00600
00601 };
00602
00603 #undef Monoblock_modular
00604 #undef Series
00605 #undef Series_rep
00606 #undef C
00607 #undef Modulus
00608 #undef TMPL
00609 #undef Monoblock_transformer
00610 }
00611 #endif // __MMX_SERIES_CARRY_BLOCKS_HPP
