00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_ANALYTIC_HPP
00014 #define __MMX_ANALYTIC_HPP
00015 #include <basix/vector.hpp>
00016 #include <basix/table.hpp>
00017 #include <basix/pair.hpp>
00018 #include <numerix/floating.hpp>
00019 #include <numerix/complex_double.hpp>
00020 #include <numerix/ball_complex.hpp>
00021 #include <algebramix/series.hpp>
00022 #include <algebramix/series_elementary.hpp>
00023 #include <analyziz/newton_polygon.hpp>
00024 #include <analyziz/polynomial_numeric.hpp>
00025
00026 namespace mmx {
00027 class std_analytic {};
00028 #define TMPL_DEF template<typename C,typename V= std_analytic>
00029 #define TMPL template<typename C, typename V>
00030 #define TMPLK template<typename C, typename V, typename K>
00031 #define Format format<C>
00032 #define Analytic analytic<C,V>
00033 #define Analytic_rep analytic_rep<C,V>
00034 #define Recursive_analytic_rep recursive_analytic_rep<C,V>
00035 #define Assume assume_bound<typename unvectorize<R >::val>
00036 #define Scalar typename unvectorize<C>::val
00037 #define Radius Abs_type(Scalar)
00038 #define Series series<C>
00039 #define Series_rep series_rep<C>
00040 #define Recursive_series_rep recursive_series_rep<C>
00041 #define Polynomial polynomial<C>
00042 #define R Abs_type(C)
00043 TMPL class analytic;
00044 TMPL Analytic lshiftz (const Analytic& f, const int& shift= 1);
00045 TMPL Analytic rshiftz (const Analytic& f, const int& shift= 1);
00046 typedef long long unsigned int llnat;
00047
00048
00049
00050
00051
00052 extern double mmx_order_ratio;
00053 extern double mmx_radius_ratio;
00054
00055
00056
00057
00058
00059 #define CACHE_EXPANSION 1
00060 #define CACHE_TAIL_BOUND 2
00061 #define CACHE_RADIUS_BOUND 4
00062 #define CACHE_MOVE 8
00063 #define CACHE_EVAL 16
00064 #define CACHE_ALL 32
00065
00066 template<typename C>
00067 struct unvectorize {
00068 typedef C val;
00069 static inline vector<val> encode (const C& x) {
00070 return vector<val> (x, 1); }
00071 static inline C decode (const vector<val>& v, const C& gauge) {
00072 (void) gauge;
00073 return v[0]; }
00074 };
00075
00076 TMPL
00077 class analytic_cache {
00078 MMX_ALLOCATORS
00079 public:
00080 Series expansion;
00081 Radius focus;
00082 table<llnat,nat> assumption;
00083 table<R,nat> head_bound;
00084 table<R,nat> tail_bound;
00085 table<R,nat> conv_bound;
00086 table<Radius,nat> radius_bound;
00087 Scalar next;
00088 Analytic* move_to;
00089 C* eval_at;
00090 };
00091
00092 template<typename C>
00093 class assume_bound {
00094 MMX_ALLOCATORS
00095 public:
00096 llnat serial;
00097 table<vector<C>,void*> assumption;
00098 table<vector<C>,void*> conclusion;
00099 inline assume_bound () {
00100 static llnat current_serial= 1;
00101 serial= current_serial++; }
00102 };
00103
00104 template<typename C> C
00105 safe_div (const C& n, const C& d) {
00106 if (d == 0) return largest_cst (get_format (n));
00107 else return n / d;
00108 }
00109
00110 template<typename C,typename T> C
00111 minimum_quotients (const table<vector<C>,T>& t, const table<vector<C>,T>& u) {
00112 C m= largest_cst (get_format1 (CF1(t)));
00113 for (iterator<T> it= entries (t); busy (it); ++it)
00114 if (contains (u, *it) && N(u[*it]) == N(t[*it]))
00115 for (nat i=0; i<N(t[*it]); i++)
00116 m= min (m, safe_div (t[*it][i], u[*it][i]));
00117 return m;
00118 }
00119
00120 template<typename C> inline bool
00121 less (const C& x, const C& y, bool sure) {
00122 (void) sure; return x <= y;
00123 }
00124
00125 template<typename C> inline bool
00126 less (const ball<C>& x, const ball<C>& y, bool sure) {
00127 if (sure) return upper (x) <= lower (y);
00128 else return lower (x) <= upper (y);
00129 }
00130
00131 template<typename C,typename T> bool
00132 less (const table<vector<C>,T>& t, const table<vector<C>,T>& u, bool sure) {
00133 for (iterator<T> it= entries (t); busy (it); ++it) {
00134 if (!contains (u, *it) || N(u[*it]) != N(t[*it])) return false;
00135 for (nat i=0; i<N(t[*it]); i++)
00136 if (!less (t[*it][i], u[*it][i], sure)) return false;
00137 }
00138 return true;
00139 }
00140
00141 template<typename C,typename T> nat
00142 total_size (const table<vector<C>,T>& t) {
00143 nat sz= 0;
00144 for (iterator<T> it= entries (t); busy (it); ++it)
00145 sz += N(t[*it]);
00146 return sz;
00147 }
00148
00149
00150
00151
00152
00153 TMPL_DEF
00154 class analytic_rep REP_STRUCT_1(C) {
00155 public:
00156 analytic_cache<C,V>* cache;
00157 void init_cache ();
00158 void init_focus (const Radius& r);
00159 void init_next (const Scalar& z);
00160 virtual void Clear_cache (nat which) const;
00161
00162 public:
00163 inline analytic_rep (const Format& fm):
00164 Format (fm), cache (NULL) {}
00165 virtual ~analytic_rep () {
00166 if (cache != NULL) {
00167 if (cache->move_to != NULL) mmx_delete_one (cache->move_to);
00168 if (cache->eval_at != NULL) mmx_delete_one (cache->eval_at);
00169 delete cache;
00170 } }
00171 virtual syntactic expression (const syntactic& z) const;
00172
00173 public:
00174 virtual Series Expand () const = 0;
00175 virtual Radius Radius_bound (nat order) const;
00176 virtual R Tail_bound (const Radius& r, nat ord, Assume& a) const = 0;
00177 virtual Analytic Move (const Scalar& z) const = 0;
00178 virtual C Eval (const Scalar& z) const;
00179 virtual Analytic Move (const list<Scalar >& z) const;
00180 virtual C Eval (const list<Scalar >& z) const;
00181 virtual Analytic Derive () const = 0;
00182 inline Analytic me () const;
00183 friend class Analytic;
00184
00185 virtual C initial () const;
00186 virtual Analytic equation () const;
00187 };
00188
00189 TMPL_DEF
00190 class analytic {
00191 INDIRECT_PROTO_2 (analytic, analytic_rep, C, V)
00192 public:
00193 analytic ();
00194 analytic (const Format& fm);
00195 analytic (int c);
00196 analytic (int c, const Format& fm);
00197 analytic (const C& c);
00198 analytic (const polynomial<C>& P);
00199 analytic (const vector<C>& coeffs);
00200 inline C operator [] (nat n) const;
00201
00202 static inline generic get_variable_name () {
00203 return Series::get_variable_name (); }
00204 static inline void set_variable_name (const generic& x) {
00205 Series::set_variable_name (x); }
00206
00207 static inline nat get_output_order () {
00208 return Series::get_output_order (); }
00209 static inline void set_output_order (const nat& x) {
00210 Series::set_output_order (x); }
00211
00212 static inline nat get_cancel_order () {
00213 return Series::get_cancel_order (); }
00214 static inline void set_cancel_order (const nat& x) {
00215 Series::set_cancel_order (x); }
00216
00217 static inline bool get_formula_output () {
00218 return Series::get_formula_output (); }
00219 static inline void set_formula_output (const bool& x) {
00220 Series::set_formula_output (x); }
00221 };
00222 INDIRECT_IMPL_2 (analytic, analytic_rep, typename C, C, typename V, V)
00223 DEFINE_UNARY_FORMAT_2 (analytic)
00224
00225 STYPE_TO_TYPE(TMPL,scalar_type,Analytic,C);
00226
00227 TMPL inline Format CF (const Analytic& f) { return f->tfm (); }
00228
00229
00230
00231
00232
00233 TMPL Analytic operator + (const Analytic& f, const Analytic& g);
00234 TMPL Analytic operator - (const Analytic& f);
00235 TMPL Analytic operator - (const Analytic& f, const Analytic& g);
00236 TMPLK Analytic operator * (const K& c, const Analytic& f);
00237 TMPL Analytic operator * (const Analytic& f, const C& c);
00238 TMPL Analytic operator / (const Analytic& f, const C& c);
00239 TMPL Analytic operator * (const Analytic& f, const Analytic& g);
00240 TMPL Analytic operator / (const Analytic& f, const Analytic& g);
00241
00242 TMPL Analytic integrate (const Analytic& f);
00243 TMPL Analytic integrate_init (const Analytic& f, const C& c);
00244
00245 TMPL Analytic sqrt (const Analytic& f);
00246 TMPL Analytic sqrt_init (const Analytic& f, const C& c);
00247 TMPL Analytic exp (const Analytic& f);
00248 TMPL Analytic log (const Analytic& f);
00249 TMPL Analytic log_init (const Analytic& f, const C& c);
00250 TMPL Analytic cos (const Analytic& f);
00251 TMPL Analytic sin (const Analytic& f);
00252 TMPL Analytic tan (const Analytic& f);
00253 TMPL Analytic acos (const Analytic& f);
00254 TMPL Analytic asin (const Analytic& f);
00255 TMPL Analytic atan (const Analytic& f);
00256
00257 TMPL Analytic fast_eval (const Analytic& f);
00258 TMPL Analytic cache_last (const Analytic& f);
00259 TMPL Analytic radial (const Analytic& f);
00260 TMPL Analytic heuristic (const polynomial<C>& f);
00261
00262
00263
00264
00265
00266 TMPL void
00267 Analytic_rep::init_cache () {
00268 if (cache == NULL) {
00269 cache= new analytic_cache<C,V> ();
00270 cache->expansion = Expand ();
00271 cache->assumption= table<llnat,nat> (0);
00272 cache->head_bound= table<R,nat> ();
00273 cache->tail_bound= table<R,nat> ();
00274 cache->conv_bound= table<R,nat> ();
00275 cache->move_to = NULL;
00276 cache->eval_at = NULL;
00277 }
00278 }
00279
00280 TMPL void
00281 Analytic_rep::init_focus (const Radius& r) {
00282 if (cache == NULL) {
00283 init_cache ();
00284 cache->focus= r;
00285 }
00286 else if (hard_neq (cache->focus, r)) {
00287 cache->focus = r;
00288 cache->assumption= table<llnat,nat> (0);
00289 cache->head_bound= table<R,nat> ();
00290 cache->tail_bound= table<R,nat> ();
00291 cache->conv_bound= table<R,nat> ();
00292 }
00293 }
00294
00295 TMPL void
00296 Analytic_rep::init_next (const Scalar& z) {
00297 if (cache == NULL) {
00298 init_cache ();
00299 cache->next= z;
00300 }
00301 else if (hard_neq (cache->next, z)) {
00302 cache->next= z;
00303 if (cache->move_to != NULL) {
00304 mmx_delete_one (cache->move_to);
00305 cache->move_to= NULL;
00306 }
00307 if (cache->eval_at != NULL) {
00308 mmx_delete_one (cache->eval_at);
00309 cache->eval_at= NULL;
00310 }
00311 }
00312 }
00313
00314 TMPL void
00315 Analytic_rep::Clear_cache (nat which) const {
00316 if (which == CACHE_ALL) {
00317 if (cache != NULL) {
00318 if (cache->move_to != NULL) mmx_delete_one (cache->move_to);
00319 if (cache->eval_at != NULL) mmx_delete_one (cache->eval_at);
00320 delete cache;
00321 const_cast<Analytic_rep*> (this) -> cache = NULL;
00322 }
00323 }
00324 else {
00325 if ((which & CACHE_EXPANSION) != 0)
00326 cache->expansion= Series ();
00327 if ((which & CACHE_TAIL_BOUND) != 0) {
00328 cache->assumption= table<llnat,nat> (0);
00329 cache->head_bound= table<R,nat> ();
00330 cache->tail_bound= table<R,nat> ();
00331 cache->conv_bound= table<R,nat> ();
00332 }
00333 if ((which & CACHE_RADIUS_BOUND) != 0)
00334 cache->radius_bound= table<Radius,nat> ();
00335 if ((which & CACHE_MOVE) != 0 && cache->move_to != NULL) {
00336 mmx_delete_one (cache->move_to);
00337 cache->move_to= NULL;
00338 }
00339 if ((which & CACHE_EVAL) != 0 && cache->eval_at != NULL) {
00340 mmx_delete_one (cache->eval_at);
00341 cache->eval_at= NULL;
00342 }
00343 }
00344 }
00345
00346 TMPL void
00347 clear_cache (const Analytic& f, nat which) {
00348 if (f->cache != NULL) inside (f) -> Clear_cache (which);
00349 }
00350
00351 TMPL Series
00352 expand (const Analytic& f) {
00353 inside (f) -> init_cache ();
00354 return f->cache->expansion;
00355 }
00356
00357 TMPL Series
00358 expand (const Analytic& f, nat order) {
00359 Series r= expand (f);
00360 if (order>0) (void) r[order-1];
00361 return r;
00362 }
00363
00364 TMPL Polynomial
00365 truncate (const Analytic& f, nat n) {
00366 return truncate (expand (f), n);
00367 }
00368
00369 TMPL Polynomial
00370 range (const Analytic& f, nat i, nat j) {
00371 return range (expand (f), i, j);
00372 }
00373
00374 TMPL Analytic
00375 derive (const Analytic& f) {
00376 return f->Derive ();
00377 }
00378
00379 TMPL Analytic
00380 xderive (const Analytic& f) {
00381 return lshiftz (f->Derive ());
00382 }
00383
00384
00385
00386
00387
00388 TMPL generic
00389 var (const Analytic& f) {
00390 return var (expand (f));
00391 }
00392
00393 TMPL inline syntactic
00394 Analytic_rep::expression (const syntactic& z) const {
00395 return flatten (expand (me ()), z);
00396 }
00397
00398 TMPL inline Analytic
00399 Analytic_rep::me () const {
00400 return Analytic (this, true);
00401 }
00402
00403 TMPL C
00404 Analytic_rep::initial () const {
00405 ERROR ("not implemented (initial)");
00406 }
00407
00408 TMPL Analytic
00409 Analytic_rep::equation () const {
00410 ERROR ("not implemented (equation)");
00411 }
00412
00413 TMPL inline C
00414 Analytic::operator [] (nat n) const {
00415 return expand (*this) [n];
00416 }
00417
00418 HARD_TO_EXACT_IDENTITY_SUGAR(TMPL,Analytic)
00419
00420 TMPL inline nat hash (const Analytic& f) {
00421 return 7002007 + hash (expand (f)); }
00422 TMPL inline bool operator == (const Analytic& f, const Analytic& g) {
00423 return expand (f) == expand (g); }
00424 TMPL inline bool operator != (const Analytic& f, const Analytic& g) {
00425 return expand (f) != expand (g); }
00426
00427 TMPL syntactic
00428 flatten (const Analytic& f, const syntactic& z) {
00429 if (Series::get_formula_output ()) return f->expression (z);
00430 else return flatten (expand (f), z);
00431 }
00432
00433 TMPL syntactic
00434 flatten (const Analytic& f) {
00435 return flatten (f, as_syntactic (var (f)));
00436 }
00437
00438
00439
00440
00441
00442
00443
00444 TMPL R
00445 tail_bound (const Analytic& f, const Radius& r, nat order, Assume& a) {
00446 Analytic_rep* rep= inside (f);
00447 rep->init_focus (r);
00448 if (read (rep->cache->assumption, order) != a.serial) {
00449 (void) expand (f, order);
00450
00451 R b= f->Tail_bound (r, order, a);
00452 rep->cache->assumption [order]= a.serial;
00453 rep->cache->tail_bound [order]= b;
00454 }
00455 return read (rep->cache->tail_bound, order);
00456 }
00457
00458 TMPL R
00459 tail_bound (const Analytic& f, const Radius& r, nat order) {
00460 Radius inc_ord= promote ((int) (order + 1), r);
00461 Radius big_one= promote (1, r) + invert (promote (4, r) * inc_ord);
00462 Assume a1;
00463 R b= tail_bound (f, r, order, a1);
00464 if (N (a1.assumption) == 0 || r == 0) {
00465 if (is_nan (b) && r == 0 && order > 0) set_zero (b);
00466 return b;
00467 }
00468
00469
00470 nat iterations;
00471 if (order < 4 ) iterations= 4;
00472 else if (order < 16 ) iterations= ((order + 1) >> 1) + 2;
00473 else if (order < 64 ) iterations= ((order + 1) >> 2) + 6;
00474 else if (order < 256) iterations= ((nat) (sqrt ((double) order))) + 14;
00475 else iterations= ((nat) (sqrt (sqrt ((double) order)))) + 26;
00476 for (nat i=0; i<iterations; i++) {
00477
00478
00479 Assume a2;
00480 a2.assumption= a1.conclusion;
00481 if (is_reliable (Radius (0))) a2.assumption= sharpen (a2.assumption);
00482 b= tail_bound (f, r, order, a2);
00483
00484 if (!(a2.assumption <= big_one * a2.conclusion)) {
00485 set_infinity (b);
00486 return b;
00487 }
00488 if (less (a2.conclusion, a2.assumption, true)) return b;
00489
00490
00491 Assume a3;
00492 table<vector<Radius >,void*> d1= a1.conclusion - a1.assumption;
00493 table<vector<Radius >,void*> d2= a2.conclusion - a2.assumption;
00494 if (d1 > d2) {
00495 Radius l = minimum_quotients (d1, d1 - d2);
00496 a3.assumption= a1.assumption + l * (a1.conclusion - a1.assumption);
00497 if (is_reliable (Radius (0))) a3.assumption= sharpen (a3.assumption);
00498 b= tail_bound (f, r, order, a3);
00499
00500 if (!(a3.assumption <= big_one * a3.conclusion)) {
00501 set_infinity (b);
00502 return b;
00503 }
00504 if (less (a3.conclusion, a3.assumption, true)) return b;
00505 }
00506 else a3= a2;
00507
00508
00509 Assume a4;
00510
00511 a4.assumption= big_one * a3.conclusion;
00512 if (is_reliable (Radius (0))) {
00513 Assume a5= a4;
00514 nat sz= total_size (a4.assumption);
00515 for (nat i=0; i<=sz; i++) {
00516 a5.assumption= sharpen (a5.assumption);
00517 b= tail_bound (f, r, order, a5);
00518
00519 if (less (a5.conclusion, a5.assumption, true)) return b;
00520 if (!less (a5.conclusion, a5.assumption, false)) break;
00521 Assume a6;
00522 a6.assumption= (a4.assumption + a5.conclusion) / promote (2, r);
00523 a5= a6;
00524 }
00525 }
00526 else {
00527 b= tail_bound (f, r, order, a4);
00528
00529 if (less (a4.conclusion, a4.assumption, true)) return b;
00530 }
00531
00532
00533 a1.assumption= a3.assumption;
00534 a1.conclusion= a3.conclusion;
00535 }
00536 set_infinity (b);
00537 return b;
00538 }
00539
00540
00541
00542
00543
00544 TMPL R
00545 head_extremum (const Analytic& f, const Radius& r, nat order, int type) {
00546 Series s= expand (f, order);
00547 return extremum (truncate (s, order), r, order, type);
00548 }
00549
00550 TMPL R
00551 head_bound (const Analytic& f, const Radius& r, nat order) {
00552 Analytic_rep* rep= inside (f);
00553 rep->init_focus (r);
00554 if (!rep->cache->head_bound->contains (order))
00555 rep->cache->head_bound [order]=
00556 head_extremum (f, r, order, MAXIMUM_DISK);
00557 return read (rep->cache->head_bound, order);
00558 }
00559
00560 TMPL R
00561 conv_bound (const Analytic& me,
00562 const Analytic& f, const Analytic& g,
00563 const Radius& r, nat order)
00564 {
00565 Analytic_rep* rep= inside (me);
00566 rep->init_focus (r);
00567 if (!rep->cache->conv_bound->contains (order)) {
00568 R sum=0, cumul= 0;
00569 Series ff= expand (f, order);
00570 Series gg= expand (g, order);
00571 for (nat i=1; i<order; i++) {
00572 cumul = r * cumul + abs (gg[order-i]);
00573 sum += abs (ff[i]) * cumul;
00574 }
00575 rep->cache->conv_bound [order]= sum * pow (r, promote ((int) order, r));
00576 }
00577 return read (rep->cache->conv_bound, order);
00578 }
00579
00580 TMPL inline nat
00581 default_order (const Analytic& f) {
00582 (void) f;
00583 return (nat) (mmx_order_ratio * ((double) precision (Radius (0))));
00584 }
00585
00586 TMPL R
00587 upper_bound (const Analytic& f, const Radius& r, nat order, Assume& a) {
00588 return head_bound (f, r, order) + tail_bound (f, r, order, a);
00589 }
00590
00591 TMPL R
00592 upper_bound (const Analytic& f, const Radius& r, nat order) {
00593 return head_bound (f, r, order) + tail_bound (f, r, order);
00594 }
00595
00596 TMPL R
00597 upper_bound (const Analytic& f, const Radius& r) {
00598 return upper_bound (f, r, default_order (f));
00599 }
00600
00601 TMPL R
00602 lower_bound (const Analytic& f, const Radius& r, nat order, Assume& a) {
00603 return max (head_extremum (f, r, order, MINIMUM_DISK) -
00604 tail_bound (f, r, order, a), R(0));
00605 }
00606
00607 TMPL R
00608 lower_bound (const Analytic& f, const Radius& r, nat order) {
00609 return max (head_extremum (f, r, order, MINIMUM_DISK) -
00610 tail_bound (f, r, order), R(0));
00611 }
00612
00613 TMPL R
00614 lower_bound (const Analytic& f, const Radius& r) {
00615 return lower_bound (f, r, default_order (f));
00616 }
00617
00618
00619
00620
00621
00622 template<typename C> inline C
00623 big_max (const C& x) {
00624 return x;
00625 }
00626
00627 template<typename C> Radius
00628 heuristic_radius (const polynomial<C>& p, nat order) {
00629 Radius* a= mmx_new<Radius > (order);
00630 for (nat i=0; i<order; i++)
00631 a[i]= big_max (abs (p[i]));
00632 nat n= order;
00633 while (n>0 && is_exact_zero (a[n-1])) n--;
00634 if (n < (order>>2) || n == 0) return Maximal (Radius );
00635 Radius r= invert (newton_scale (a + (n>>1), n - (n>>1)));
00636
00637 mmx_delete<Radius > (a, order);
00638 return r;
00639 }
00640
00641 template<typename C> inline nat
00642 default_order (const polynomial<C>& p) {
00643 (void) p;
00644 return (nat) (mmx_order_ratio * ((double) precision (Radius (0))));
00645 }
00646
00647 template<typename C> Radius
00648 heuristic_radius (const polynomial<C>& p) {
00649 return heuristic_radius (p, default_order (p));
00650 }
00651
00652 TMPL Radius
00653 Analytic_rep::Radius_bound (nat order) const {
00654 Radius r= sharpen (heuristic_radius (truncate (me (), order), order));
00655 if (r <= promote (0, r)) return promote (0, r);
00656 if (is_infinite (r)) {
00657 r= Largest (Radius);
00658 while (true) {
00659 R t= tail_bound (me (), r, order);
00660
00661 if (is_finite (t)) break;
00662 r= sqrt (r);
00663 }
00664 }
00665 Radius s= r / promote (2, r);
00666 Radius b= promote (0, r);
00667 bool ok= false;
00668 nat credit= order;
00669 while (credit != 0) {
00670 R t= tail_bound (me (), r, order);
00671
00672 if (is_finite (t)) {
00673 ok = true;
00674 b = r;
00675 r += s;
00676
00677 }
00678 else r -= s;
00679 s= s / promote (2, s);
00680 if (ok) credit= credit >> 1;
00681 else credit= credit - 1;
00682 }
00683 if (ok) return b;
00684 else return promote (0, r);
00685 }
00686
00687 TMPL Radius
00688 radius_bound (const Analytic& f, nat order) {
00689 Analytic_rep* rep= inside (f);
00690 rep->init_cache ();
00691 if (!rep->cache->radius_bound->contains (order))
00692 rep->cache->radius_bound [order]= f->Radius_bound (order);
00693 return read (rep->cache->radius_bound, order);
00694 }
00695
00696 TMPL Radius
00697 radius_bound (const Analytic& f) {
00698 return radius_bound (f, default_order (f));
00699 }
00700
00701
00702
00703
00704
00705 TMPL Analytic
00706 move (const Analytic& f, const Scalar& z) {
00707 Analytic_rep* rep= inside (f);
00708 rep->init_next (z);
00709 if (rep->cache->move_to == NULL)
00710 rep->cache->move_to= mmx_new_one<Analytic > (rep->Move (z));
00711 return *rep->cache->move_to;
00712 }
00713
00714 TMPL inline C
00715 default_eval (const Analytic& f, const Scalar& z) {
00716 nat order= max (default_order (f), (nat) 2);
00717 C sum= f[order-1];
00718 for (int i= order-2; i>=0; i--)
00719 sum= sum * z + f[i];
00720 if (is_reliable (sum))
00721 sum= blur (sum, tail_bound (f, abs(z), order));
00722 return sum;
00723 }
00724
00725 TMPL C
00726 Analytic_rep::Eval (const Scalar& z) const {
00727 return default_eval (me (), z);
00728 }
00729
00730 TMPL C
00731 eval (const Analytic& f, const Scalar& z) {
00732 Analytic_rep* rep= inside (f);
00733 rep->init_next (z);
00734 if (rep->cache->eval_at == NULL)
00735 rep->cache->eval_at= mmx_new_one<C> (rep->Eval (z));
00736 return *rep->cache->eval_at;
00737 }
00738
00739
00740
00741
00742
00743 template<typename C> list<C>
00744 simplify (const list<C>& l, const R& r) {
00745
00746
00747
00748
00749 if (is_nil (l)) return l;
00750 else if (abs (car (l)) < r) {
00751 if (is_nil (cdr (l))) return l;
00752 format<C> fm= get_format (car (l));
00753 C p1= car (l);
00754 C dp= car (cdr (l));
00755 C p2= p1 + dp;
00756 list<C> rl= cdr (cdr (l));
00757 if (abs (p2) < r)
00758 return simplify (cons<C> (p2, rl), r);
00759 else {
00760 C _2= promote (2, fm);
00761 C _4= promote (4, fm);
00762 C a = gaussian (square (Re (dp)) + square (Im (dp)), fm);
00763 C b = _2 * (Re (p1) * Im (dp) + Im (p1) * Re (dp));
00764 C c = gaussian (square (Re (p1)) + square (Im (p1)) - square (r), fm);
00765 C t = (-b + sqrt (square (b) - _4 * a * c)) / (_2 * a);
00766 C q1= p1 + t * dp;
00767 C q2= p2 - q1;
00768 return cons<C> (q1, cons<C> (q2, rl));
00769 }
00770 }
00771 else {
00772 C p = car (l);
00773 C q1= (p / abs (p)) * r;
00774 C q2= p - q1;
00775 return cons<C> (q1, cons<C> (q2, cdr (l)));
00776 }
00777 }
00778
00779 template<typename C, typename FM>
00780 struct radius_ratio_helper {
00781 static inline C op (const format<C>& fm) {
00782 return promote (4, fm) / promote (10, fm); }
00783 };
00784
00785 template<typename C>
00786 struct radius_ratio_helper<C,empty_format> {
00787 static inline C op (const format<C>& fm) {
00788 return as<C> (mmx_radius_ratio); }
00789 };
00790
00791 template<typename C> C
00792 shrink (const C& z) {
00793 typedef radius_ratio_helper<C,typename format<C>::FT> H;
00794 return H::op (get_format (z)) * z;
00795 }
00796
00797 TMPL Analytic
00798 default_move (const Analytic& f2, const list<Scalar >& l2) {
00799 Analytic f= f2;
00800 list<Scalar > l= l2;
00801 while (!is_nil (l)) {
00802 Radius r= shrink (radius_bound (f));
00803 l= simplify (l, r);
00804 Analytic g= move (f, car (l));
00805
00806 f= g;
00807 l= cdr (l);
00808 }
00809 return f;
00810 }
00811
00812 TMPL C
00813 default_eval (const Analytic& f2, const list<Scalar >& l2) {
00814 Analytic f= f2;
00815 list<Scalar > l= l2;
00816 while (!is_nil (l)) {
00817 Radius r= shrink (radius_bound (f));
00818 l= simplify (l, r);
00819 Analytic g= move (f, car (l));
00820
00821 f= g;
00822 l= cdr (l);
00823 }
00824 return f[0];
00825 }
00826
00827 TMPL Analytic
00828 Analytic_rep::Move (const list<Scalar >& z) const {
00829 return default_move (me (), z);
00830 }
00831
00832 TMPL C
00833 Analytic_rep::Eval (const list<Scalar >& z) const {
00834 return default_eval (me (), z);
00835 }
00836
00837 TMPL Analytic
00838 move (const Analytic& f, const list<Scalar >& l) {
00839 return f->Move (l);
00840 }
00841
00842 TMPL C
00843 eval (const Analytic& f, const list<Scalar >& l) {
00844 return f->Eval (l);
00845 }
00846
00847 TMPL Analytic
00848 move (const Analytic& f, const vector<Scalar >& v) {
00849 return f->Move (list<Scalar > (v));
00850 }
00851
00852 TMPL C
00853 eval (const Analytic& f, const vector<Scalar >& v) {
00854 return f->Eval (list<Scalar > (v));
00855 }
00856
00857 TMPL Analytic
00858 radial_move (const Analytic& f, const Scalar& z) {
00859 return f->Move (list<Scalar > (z));
00860 }
00861
00862 TMPL C
00863 radial_eval (const Analytic& f, const Scalar& z) {
00864 return f->Eval (list<Scalar > (z));
00865 }
00866
00867
00868
00869
00870
00871 TMPL class recursive_container_analytic_rep;
00872
00873 TMPL_DEF
00874 class recursive_analytic_rep: public Analytic_rep {
00875 protected:
00876 C init;
00877 Analytic equa;
00878 public:
00879 inline recursive_analytic_rep (const Format& fm):
00880 Analytic_rep (fm) {}
00881 virtual void Clear_cache (nat which) const {
00882 Analytic_rep::Clear_cache (which);
00883 clear_cache (equa, which); }
00884 virtual syntactic expression (const syntactic& z) const {
00885 (void) z; return "recursive"; }
00886 virtual C Initial () const = 0;
00887 virtual Analytic Equation () const = 0;
00888 virtual R Tail_bound (const Radius& r, nat order, Assume& a) const {
00889 void* code= (void*) this;
00890 R zero= promote (0, r) * abs (init);
00891 if (!a.assumption->contains (code))
00892 a.assumption [code]= unvectorize<R>::encode (zero);
00893 this->cache->assumption [order]= a.serial;
00894 this->cache->tail_bound [order]=
00895 unvectorize<R>::decode (read (a.assumption, code), zero);
00896 R bnd= tail_bound (this->equa, r, order, a);
00897 if (order == 0) bnd += abs (this->initial ());
00898 this->cache->tail_bound [order]= bnd;
00899 a.conclusion [code]= unvectorize<R>::encode (bnd);
00900 return bnd; }
00901 friend class recursive_container_analytic_rep<C,V>;
00902 };
00903
00904 TMPL
00905 class recursive_container_analytic_rep: public Analytic_rep {
00906 Analytic f;
00907 public:
00908 recursive_container_analytic_rep (const Analytic& f2):
00909 Analytic_rep (CF(f2)), f (f2) {
00910 Recursive_analytic_rep* rep= (Recursive_analytic_rep*) f.operator -> ();
00911 rep->init= rep->Initial ();
00912 rep->equa= rep->Equation (); }
00913 ~recursive_container_analytic_rep () {
00914 Recursive_analytic_rep* rep= (Recursive_analytic_rep*) f.operator -> ();
00915 rep->equa= Analytic (); }
00916 syntactic expression (const syntactic& z) const {
00917 return flatten (f, z); }
00918 Series Expand () const {
00919 return expand (f); }
00920 Radius Radius_bound (nat order) const {
00921 return radius_bound (f, order); }
00922 R Tail_bound (const Radius& r, nat ord, Assume& a) const {
00923 return tail_bound (f, r, ord, a); }
00924 Analytic Move (const Scalar& z) const {
00925 return move (f, z); }
00926 C Eval (const Scalar& z) const {
00927 return eval (f, z); }
00928 Analytic Move (const list<Scalar >& l) const {
00929 return move (f, l); }
00930 C Eval (const list<Scalar >& l) const {
00931 return eval (f, l); }
00932 Analytic Derive () const {
00933 return derive (f); }
00934 };
00935
00936 TMPL Analytic
00937 recursive (recursive_analytic_rep<C,V>* f) {
00938 return (Analytic_rep*)
00939 new recursive_container_analytic_rep<C,V> (Analytic ((Analytic_rep*) f));
00940 }
00941
00942
00943
00944
00945
00946 TMPL
00947 class zero_analytic_rep: public Analytic_rep {
00948 public:
00949 zero_analytic_rep (const Format& fm):
00950 Analytic_rep (fm) {}
00951 Series Expand () const { return Series (this->tfm ()); }
00952 Radius Radius_bound (nat order) const { return Maximal (Radius); }
00953 R Tail_bound (const Radius& r, nat order, Assume& a) const {
00954 return abs (promote (0, this->tfm ())); }
00955 Analytic Move (const Scalar& z) const {
00956 return Analytic (this->tfm ()); }
00957 C Eval (const Scalar& z) const {
00958 return promote (0, this->tfm ()); }
00959 Analytic Derive () const {
00960 return Analytic (this->tfm ()); }
00961 };
00962
00963 TMPL Analytic::analytic () {
00964 rep= new zero_analytic_rep<C,V> (format<C> (no_format ())); }
00965 TMPL Analytic::analytic (const Format& fm) {
00966 rep= new zero_analytic_rep<C,V> (fm); }
00967
00968 TMPL
00969 class polynomial_analytic_rep: public Analytic_rep {
00970 polynomial<C> p;
00971 public:
00972 polynomial_analytic_rep (const polynomial<C>& p2):
00973 Analytic_rep (CF(p2)), p (p2) {}
00974 Series Expand () const { return Series (p); }
00975 Radius Radius_bound (nat order) const { return Maximal (Radius); }
00976 R Tail_bound (const Radius& r, nat order, Assume& a) const {
00977 (void) a;
00978 (void) expand (this->me (), order);
00979 R sum= 0;
00980 for (int i= ((int) N(p)) - 1; i >= ((int) order); i--)
00981 sum= r * sum + abs (p[i]);
00982 return sum * pow (r, promote ((int) order, r)); }
00983 Analytic Move (const Scalar& z) const {
00984 return Analytic (shift (p, C (z))); }
00985 C Eval (const Scalar& z) const {
00986 return eval (p, C (z)); }
00987
00988
00989 Analytic Derive () const {
00990 return Analytic (derive (p)); }
00991 };
00992
00993 TMPL
00994 Analytic::analytic (int c) {
00995 rep= new polynomial_analytic_rep<C,V> (polynomial<C> (C (c), 0));
00996 }
00997
00998 TMPL
00999 Analytic::analytic (int c, const Format& fm) {
01000 rep= new polynomial_analytic_rep<C,V> (polynomial<C> (promote (c, fm), 0));
01001 }
01002
01003 TMPL
01004 Analytic::analytic (const C& c) {
01005 rep= new polynomial_analytic_rep<C,V> (polynomial<C> (c, 0));
01006 }
01007
01008 TMPL
01009 Analytic::analytic (const polynomial<C>& p) {
01010 rep= new polynomial_analytic_rep<C,V> (p);
01011 }
01012
01013 TMPL
01014 Analytic::analytic (const vector<C>& coeffs) {
01015 rep= new polynomial_analytic_rep<C,V> (polynomial<C> (coeffs));
01016 }
01017
01018 template<typename Real> inline analytic<Complex_type(Real) >
01019 analytic_complex (const vector<Real>& v) {
01020 return analytic<Complex_type(Real) > (as<vector<Complex_type(Real) > > (v));
01021 }
01022
01023
01024
01025
01026
01027 template<typename Op,typename C> inline analytic<C>
01028 nullary_recursive_analytic (const C& c);
01029
01030 template<typename Op,typename C,typename V= std_analytic>
01031 class nullary_recursive_analytic_rep: public Recursive_analytic_rep {
01032 protected:
01033 C c;
01034 public:
01035 inline nullary_recursive_analytic_rep (const C& c2):
01036 Recursive_analytic_rep (CF(c2)), c (c2) {}
01037 syntactic expression (const syntactic& z) const { (void) z;
01038 return Op::op_init (flatten (c)); }
01039 C Initial () const {
01040 return c; }
01041 Analytic Equation () const {
01042 return Op::def (this->me ()); }
01043 Series Expand () const {
01044 return nullary_recursive_series<Op> (c); }
01045 Analytic Move (const Scalar& z) const {
01046 return nullary_recursive_analytic<Op> (eval (this->me (), z)); }
01047 C Eval (const Scalar& z) const { (void) z;
01048 return Op::def (c); }
01049 Analytic Derive () const {
01050 return Op::diff_op (this->me ()); }
01051 };
01052
01053 template<typename Op,typename C> inline analytic<C>
01054 nullary_recursive_analytic (const C& c= Op::template op<C> ()) {
01055 return recursive (new nullary_recursive_analytic_rep<Op,C> (c));
01056 }
01057
01058 template<typename Op,typename C,typename V>
01059 class unary_analytic_rep: public Analytic_rep {
01060 protected:
01061 Analytic f;
01062 public:
01063 inline unary_analytic_rep (const Analytic& f2):
01064 Analytic_rep (CF(f2)), f(f2) {}
01065 void Clear_cache (nat which) const {
01066 Analytic_rep::Clear_cache (which);
01067 clear_cache (f, which); }
01068 Series Expand () const {
01069 return Op::op (expand (f)); }
01070 Analytic Move (const Scalar& z) const {
01071 return Op::op (move (f, z)); }
01072 C Eval (const Scalar& z) const {
01073 return Op::op (eval (f, z)); }
01074 Analytic Derive () const {
01075 return Op::diff_op (this->me (), f); }
01076 };
01077
01078 template<typename Op,typename C,typename V> inline Analytic
01079 unary_analytic (const Analytic& f) {
01080 return (Analytic_rep*) new unary_analytic_rep<Op,C,V> (f);
01081 }
01082
01083 template<typename Op,typename C,typename V>
01084 class unary_recursive_analytic_rep: public Recursive_analytic_rep {
01085 protected:
01086 Analytic f;
01087 bool with_init;
01088 C c;
01089 public:
01090 inline unary_recursive_analytic_rep (const Analytic& f2):
01091 Recursive_analytic_rep (CF(f2)), f (f2), with_init (false) {}
01092 inline unary_recursive_analytic_rep (const Analytic& f2, const C& c2):
01093 Recursive_analytic_rep (CF(f2)), f (f2), with_init (true), c (c2) {}
01094 void Clear_cache (nat which) const {
01095 Recursive_analytic_rep::Clear_cache (which);
01096 clear_cache (f, which); }
01097 syntactic expression (const syntactic& z) const {
01098 if (with_init) return Op::op_init (flatten (f, z), flatten (c));
01099 else return Op::op (flatten (f, z)); }
01100 C Initial () const {
01101 if (with_init) return c;
01102 else return Op::op (expand (f) [0]); }
01103 Analytic Equation () const {
01104 return Op::def (this->me (), f); }
01105 Series Expand () const {
01106 if (with_init) return Op::op_init (expand (f), c);
01107 else return Op::op (expand (f)); }
01108 Analytic Move (const Scalar& z) const {
01109 return Op::op_init (move (f, z), eval (this->me (), z)); }
01110
01111
01112 C Eval (const Scalar& z) const {
01113 if (with_init) return default_eval (this->me (), z);
01114
01115 else return Op::op (eval (f, z)); }
01116 Analytic Derive () const {
01117 return Op::diff_op (this->me (), f); }
01118 };
01119
01120 template<typename Op,typename C,typename V> inline Analytic
01121 unary_recursive_analytic (const Analytic& f) {
01122 return recursive (new unary_recursive_analytic_rep<Op,C,V> (f));
01123 }
01124
01125 template<typename Op,typename C,typename V> inline Analytic
01126 unary_recursive_analytic (const Analytic& f, const C& c) {
01127 return recursive (new unary_recursive_analytic_rep<Op,C,V> (f, c));
01128 }
01129
01130 template<typename Op,typename C,typename V>
01131 class binary_analytic_rep: public Analytic_rep {
01132 protected:
01133 Analytic f, g;
01134 public:
01135 inline binary_analytic_rep (const Analytic& f2, const Analytic& g2):
01136 Analytic_rep (CF(f2)), f(f2), g (g2) {}
01137 void Clear_cache (nat which) const {
01138 Analytic_rep::Clear_cache (which);
01139 clear_cache (f, which);
01140 clear_cache (g, which); }
01141 Series Expand () const {
01142 return Op::op (expand (f), expand (g)); }
01143 Analytic Move (const Scalar& z) const {
01144 return Op::op (move (f, z), move (g, z)); }
01145 C Eval (const Scalar& z) const {
01146 return Op::op (eval (f, z), eval (g, z)); }
01147 Analytic Derive () const {
01148 return Op::diff_op (this->me (), f, g); }
01149 };
01150
01151 template<typename Op,typename C,typename V> inline Analytic
01152 binary_analytic (const Analytic& f, const Analytic& g) {
01153 return (Analytic_rep*) new binary_analytic_rep<Op,C,V> (f, g);
01154 }
01155
01156 template<typename Op,typename C,typename V>
01157 class binary_recursive_analytic_rep: public Recursive_analytic_rep {
01158 protected:
01159 Analytic f, g;
01160 public:
01161 inline binary_recursive_analytic_rep (const Analytic& f2, const Analytic& g2):
01162 Recursive_analytic_rep (CF(f2)), f (f2), g (g2) {}
01163 void Clear_cache (nat which) const {
01164 Recursive_analytic_rep::Clear_cache (which);
01165 clear_cache (f, which);
01166 clear_cache (g, which); }
01167 syntactic expression (const syntactic& z) const {
01168 return Op::op (flatten (f, z), flatten (g, z)); }
01169 C Initial () const {
01170 return Op::op (expand (f), expand (g)) [0]; }
01171 Analytic Equation () const {
01172 return Op::def (this->me (), f, g); }
01173 Series Expand () const {
01174 return Op::op (expand (f), expand (g)); }
01175 Analytic Move (const Scalar& z) const {
01176 return Op::op (move (f, z), move (g, z)); }
01177 C Eval (const Scalar& z) const {
01178 return Op::op (eval (f, z), eval (g, z)); }
01179 Analytic Derive () const {
01180 return Op::diff_op (this->me (), f, g); }
01181 };
01182
01183 template<typename Op,typename C,typename V> inline Analytic
01184 binary_recursive_analytic (const Analytic& f, const Analytic& g) {
01185 return recursive (new binary_recursive_analytic_rep<Op,C,V> (f, g));
01186 }
01187
01188 template<typename Op,typename C,typename V,typename X>
01189 class binary_scalar_analytic_rep: public Analytic_rep {
01190 protected:
01191 Analytic f;
01192 X x;
01193 public:
01194 inline binary_scalar_analytic_rep (const Analytic& f2, const X& x2):
01195 Analytic_rep (CF(f2)), f(f2), x (x2) {}
01196 void Clear_cache (nat which) const {
01197 Analytic_rep::Clear_cache (which);
01198 clear_cache (f, which); }
01199 Series Expand () const {
01200 return Op::op (expand (f), x); }
01201 Analytic Move (const Scalar& z) const {
01202 return Op::op (move (f, z), x); }
01203 C eval (const Scalar& z) const {
01204 return Op::op (eval (f, z), x); }
01205 Analytic Derive () const {
01206 return Op::diff_op (this->me (), f, x); }
01207 };
01208
01209 template<typename Op,typename C,typename V,typename X> inline Analytic
01210 binary_scalar_analytic (const Analytic& f, const X& x) {
01211 return (Analytic_rep*) new binary_scalar_analytic_rep<Op,C,V,X> (f, x);
01212 }
01213
01214
01215
01216
01217
01218 TMPL
01219 class add_analytic_rep: public binary_analytic_rep<add_op,C,V> {
01220 public:
01221 add_analytic_rep (const Analytic& f, const Analytic& g):
01222 binary_analytic_rep<add_op,C,V> (f, g) {}
01223 Radius Radius_bound (nat order) const {
01224 return min (radius_bound (this->f), radius_bound (this->g)); }
01225 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01226 return tail_bound (this->f, r, order, a) +
01227 tail_bound (this->g, r, order, a); }
01228 };
01229
01230 TMPL Analytic
01231 operator + (const Analytic& f, const Analytic& g) {
01232 return (Analytic_rep*) new add_analytic_rep<C,V> (f, g);
01233 }
01234
01235 TMPL Analytic
01236 operator + (const C& f, const Analytic& g) {
01237 return (Analytic_rep*) new add_analytic_rep<C,V> (Analytic (f), g);
01238 }
01239
01240 TMPL Analytic
01241 operator + (const Analytic& f, const C& g) {
01242 return (Analytic_rep*) new add_analytic_rep<C,V> (f, Analytic (g));
01243 }
01244
01245 TMPL
01246 class neg_analytic_rep: public unary_analytic_rep<neg_op,C,V> {
01247 public:
01248 neg_analytic_rep (const Analytic& f):
01249 unary_analytic_rep<neg_op,C,V> (f) {}
01250 Radius Radius_bound (nat order) const {
01251 return radius_bound (this->f); }
01252 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01253 return tail_bound (this->f, r, order, a); }
01254 };
01255
01256 TMPL Analytic
01257 operator - (const Analytic& f) {
01258 return (Analytic_rep*) new neg_analytic_rep<C,V> (f);
01259 }
01260
01261 TMPL
01262 class sub_analytic_rep: public binary_analytic_rep<sub_op,C,V> {
01263 public:
01264 sub_analytic_rep (const Analytic& f, const Analytic& g):
01265 binary_analytic_rep<sub_op,C,V> (f, g) {}
01266 Radius Radius_bound (nat order) const {
01267 return min (radius_bound (this->f), radius_bound (this->g)); }
01268 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01269 return tail_bound (this->f, r, order, a) +
01270 tail_bound (this->g, r, order, a); }
01271 };
01272
01273 TMPL Analytic
01274 operator - (const Analytic& f, const Analytic& g) {
01275 return (Analytic_rep*) new sub_analytic_rep<C,V> (f, g);
01276 }
01277
01278 TMPL Analytic
01279 operator - (const C& f, const Analytic& g) {
01280 return (Analytic_rep*) new sub_analytic_rep<C,V> (Analytic (f), g);
01281 }
01282
01283 TMPL Analytic
01284 operator - (const Analytic& f, const C& g) {
01285 return (Analytic_rep*) new sub_analytic_rep<C,V> (f, Analytic (g));
01286 }
01287
01288 TMPL
01289 class sc_mul_analytic_rep: public binary_scalar_analytic_rep<rmul_op,C,V,C> {
01290 public:
01291 sc_mul_analytic_rep (const Analytic& f, const C& c):
01292 binary_scalar_analytic_rep<rmul_op,C,V,C> (f, c) {}
01293 Radius Radius_bound (nat order) const {
01294 return radius_bound (this->f); }
01295 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01296 return tail_bound (this->f, r, order, a) * abs (this->x); }
01297 };
01298
01299 TMPL Analytic
01300 operator * (const Analytic& f, const C& c) {
01301 return (Analytic_rep*) new sc_mul_analytic_rep<C,V> (f, c);
01302 }
01303
01304 template<typename C,typename V,typename K> Analytic
01305 operator * (const K& c, const Analytic& f) {
01306 return (Analytic_rep*) new sc_mul_analytic_rep<C,V> (f, c);
01307 }
01308
01309 TMPL Analytic
01310 operator / (const Analytic& f, const C& c) {
01311 return (Analytic_rep*) new sc_mul_analytic_rep<C,V> (f, invert (c));
01312 }
01313
01314
01315
01316
01317
01318 TMPL
01319 class mul_analytic_rep: public binary_analytic_rep<mul_op,C,V> {
01320 public:
01321 mul_analytic_rep (const Analytic& f, const Analytic& g):
01322 binary_analytic_rep<mul_op,C,V> (f, g) {}
01323 Radius Radius_bound (nat order) const {
01324 return min (radius_bound (this->f), radius_bound (this->g)); }
01325 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01326 return
01327 conv_bound (this->me (), this->f, this->g, r, order) +
01328 head_bound (this->f, r, order) * tail_bound (this->g, r, order, a) +
01329 tail_bound (this->f, r, order, a) * upper_bound (this->g, r, order, a); }
01330 };
01331
01332 TMPL Analytic
01333 operator * (const Analytic& f, const Analytic& g) {
01334 return (Analytic_rep*) new mul_analytic_rep<C,V> (f, g);
01335 }
01336
01337 TMPL
01338 class div_analytic_rep: public binary_analytic_rep<div_op,C,V> {
01339 public:
01340 div_analytic_rep (const Analytic& f, const Analytic& g):
01341 binary_analytic_rep<div_op,C,V> (f, g) {}
01342
01343
01344
01345
01346
01347 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01348 R lg= lower_bound (this->g, r, order, a);
01349 if (lg <= 0) return Maximal (R);
01350 R tf= tail_bound (this->f, r, order, a);
01351 R tg= tail_bound (this->g, r, order, a);
01352 R hh= head_bound (this->me (), r, order);
01353 R cb= conv_bound (this->me (), this->me (), this->g, r, order);
01354
01355
01356
01357
01358
01359 return (tf + cb + hh * tg) / lg; }
01360 };
01361
01362 TMPL Analytic
01363 operator / (const Analytic& f, const Analytic& g) {
01364 return (Analytic_rep*) new div_analytic_rep<C,V> (f, g);
01365 }
01366
01367 TMPL Analytic
01368 operator / (const C& f, const Analytic& g) {
01369 return (Analytic_rep*) new div_analytic_rep<C,V> (Analytic (f), g);
01370 }
01371
01372 ARITH_SCALAR_INT_SUGAR(TMPL,Analytic);
01373
01374
01375
01376
01377
01378 TMPL
01379 class lshiftz_analytic_rep: public Analytic_rep {
01380 protected:
01381 Analytic f;
01382 int shift;
01383 public:
01384 lshiftz_analytic_rep (const Analytic& f2, int shift2):
01385 Analytic_rep (CF(f2)), f (f2), shift (shift2) {}
01386 void Clear_cache (nat which) const {
01387 Analytic_rep::Clear_cache (which);
01388 clear_cache (f, which); }
01389 Series Expand () const {
01390 return lshiftz (expand (f), shift); }
01391 Radius Radius_bound (nat order) const {
01392 (void) order; return radius_bound (this->f); }
01393 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01394 Series s= expand (this->f, order);
01395 R sum= tail_bound (this->f, r, order, a);
01396 if (shift > 0) {
01397 R cum= 0;
01398 for (int i=0; i < min (shift, (int) order); i++)
01399 cum= r * cum + abs (s[order-1-i]);
01400 if (((int) order) < shift) sum += cum * pow (r, promote (shift, r));
01401 else sum += cum * pow (r, promote ((int) order, r));
01402 }
01403 else {
01404 R cum= 0;
01405 for (int i= 1-shift; i >= 0; i--)
01406 cum= r * cum + abs (s[order+i]);
01407 sum += cum * pow (r, promote ((int) order, r));
01408 }
01409 return sum; }
01410 Analytic Move (const Scalar& z) const {
01411 (void) z;
01412 ERROR ("not implemented (Move)");
01413 return this->me (); }
01414 C Eval (const Scalar& z) const {
01415 if (shift > 0) return C (pow (z, promote (shift, z))) * eval (f, z);
01416 else return default_eval (this->me (), z); }
01417 Analytic Derive () const {
01418 ERROR ("not implemented (Derive)");
01419 return this->me (); }
01420 };
01421
01422 TMPL Analytic
01423 lshiftz (const Analytic& f, const int& shift) {
01424 if (shift == 0) return f;
01425 return (Analytic_rep*) new lshiftz_analytic_rep<C,V> (f, shift);
01426 }
01427
01428 TMPL Analytic
01429 rshiftz (const Analytic& f, const int& shift) {
01430 if (shift == 0) return f;
01431 return (Analytic_rep*) new lshiftz_analytic_rep<C,V> (f, -shift);
01432 }
01433
01434
01435
01436
01437
01438 TMPL
01439 class integrate_analytic_rep: public Analytic_rep {
01440 protected:
01441 Analytic f;
01442 C c;
01443 public:
01444 integrate_analytic_rep (const Analytic& f2, const C& c2):
01445 Analytic_rep (CF(f2)), f (f2), c (c2) {}
01446 void Clear_cache (nat which) const {
01447 Analytic_rep::Clear_cache (which);
01448 clear_cache (f, which); }
01449 Series Expand () const {
01450 return integrate_init (expand (f), c); }
01451 Radius Radius_bound (nat order) const {
01452 return radius_bound (this->f); }
01453 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01454 Series s= expand (this->me (), order);
01455 R b= tail_bound (this->f, r, order, a);
01456 if (order == 0) return abs (c) + b;
01457 else {
01458 Radius or0= promote ((int) order, r);
01459 Radius or1= promote ((int) order + 1, r);
01460 return (abs (this->f [order - 1]) * pow (r, or0) / or0) + (b * r / or1);
01461 } }
01462 Analytic Move (const Scalar& z) const {
01463 return integrate_init (move (f, z), eval (this->me (), z)); }
01464 Analytic Derive () const {
01465 return f; }
01466 };
01467
01468 TMPL Analytic
01469 integrate (const Analytic& f) {
01470 return (Analytic_rep*)
01471 new integrate_analytic_rep<C,V> (f, promote (0, CF(f)));
01472 }
01473
01474 TMPL Analytic
01475 integrate_init (const Analytic& f, const C& c) {
01476 return (Analytic_rep*) new integrate_analytic_rep<C,V> (f, c);
01477 }
01478
01479
01480
01481
01482
01483 TMPL Analytic
01484 sqrt (const Analytic& f) {
01485 return unary_recursive_analytic<sqrt_op> (f);
01486 }
01487
01488 TMPL Analytic
01489 sqrt_init (const Analytic& f, const C& c) {
01490 return unary_recursive_analytic<sqrt_op> (f, c);
01491 }
01492
01493 TMPL Analytic
01494 exp (const Analytic& f) {
01495 return unary_recursive_analytic<exp_op> (f);
01496 }
01497
01498 TMPL Analytic
01499 log (const Analytic& f) {
01500 return unary_recursive_analytic<log_op> (f);
01501 }
01502
01503 TMPL Analytic
01504 log_init (const Analytic& f, const C& c) {
01505 return unary_recursive_analytic<log_op> (f, c);
01506 }
01507
01508 TMPL Analytic
01509 pow (const Analytic& f, const C& c) {
01510 return exp (c * log (f));
01511 }
01512
01513 TMPL Analytic
01514 pow (const Analytic& f, const Analytic& g) {
01515 return exp (g * log (f));
01516 }
01517
01518
01519
01520 TMPL Analytic
01521 acos (const Analytic& f) {
01522 return unary_recursive_analytic<acos_op> (f);
01523 }
01524
01525 TMPL Analytic
01526 acos_init (const Analytic& f, const C& c) {
01527 return unary_recursive_analytic<acos_op> (f, c);
01528 }
01529
01530 TMPL Analytic
01531 asin (const Analytic& f) {
01532 return unary_recursive_analytic<asin_op> (f);
01533 }
01534
01535 TMPL Analytic
01536 asin_init (const Analytic& f, const C& c) {
01537 return unary_recursive_analytic<asin_op> (f, c);
01538 }
01539
01540 TMPL Analytic
01541 atan (const Analytic& f) {
01542 return unary_recursive_analytic<atan_op> (f);
01543 }
01544
01545 TMPL Analytic
01546 atan_init (const Analytic& f, const C& c) {
01547 return unary_recursive_analytic<atan_op> (f, c);
01548 }
01549
01550 INV_TRIGO_SUGAR(TMPL,Analytic)
01551 HYPER_SUGAR(TMPL,Analytic)
01552 INV_HYPER_SUGAR(TMPL,Analytic)
01553 ARG_HYPER_SUGAR(TMPL,Analytic)
01554
01555
01556
01557
01558
01559 TMPL Analytic compose (const Analytic& f, const Analytic& g);
01560 TMPL Analytic reverse (const Analytic& f);
01561
01562 TMPL
01563 class compose_analytic_rep: public Analytic_rep {
01564 Analytic f, g;
01565 public:
01566 compose_analytic_rep (const Analytic& f2, const Analytic& g2):
01567 Analytic_rep (CF(f2)), f (f2), g (g2) {}
01568 void Clear_cache (nat which) const {
01569 Analytic_rep::Clear_cache (which);
01570 clear_cache (f, which);
01571 clear_cache (g, which); }
01572 Series Expand () const {
01573 return compose (expand (f), expand (g)); }
01574 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01575 ERROR ("not yet implemented");
01576 return 0; }
01577 Analytic Move (const Scalar& z) const {
01578 Analytic tg= move (this->g, z);
01579 Analytic tf= move (this->f, tg[0]);
01580 return compose (tf, tg - tg[0]); }
01581 C Eval (const Scalar& z) const {
01582 return default_eval (this->me (), z); }
01583 Analytic Derive () const {
01584 return derive (this->g) * compose (derive (this->f), this->g); }
01585 };
01586
01587 TMPL Analytic
01588 compose (const Analytic& f, const Analytic& g) {
01589 return (Analytic_rep*) new compose_analytic_rep<C,V> (f, g);
01590 }
01591
01592 TMPL
01593 class reverse_analytic_rep: public Analytic_rep {
01594 Analytic f;
01595 public:
01596 reverse_analytic_rep (const Analytic& f2):
01597 Analytic_rep (CF(f2)), f (f2) {}
01598 void Clear_cache (nat which) const {
01599 Analytic_rep::Clear_cache (which);
01600 clear_cache (f, which); }
01601 Series Expand () const {
01602 return reverse (expand (f)); }
01603 R Tail_bound (const Radius& r, nat order, Assume& a) const {
01604 ERROR ("not yet implemented");
01605 return 0; }
01606 Analytic Move (const Scalar& a) const {
01607 C b= eval (this->me (), a);
01608 Analytic tf= move (this->f, b);
01609 return reverse (tf - tf[0]) + b; }
01610 C Eval (const Scalar& z) const {
01611 return default_eval (this->me (), z); }
01612 Analytic Derive () const {
01613 return 1 / compose (derive (this->f), this->me ()); }
01614 };
01615
01616 TMPL Analytic
01617 reverse (const Analytic& f) {
01618 return (Analytic_rep*) new reverse_analytic_rep<C,V> (f);
01619 }
01620
01621
01622
01623
01624
01625 #define mmx_analytic(R,C) analytic<C >
01626
01627 template<typename C> inline analytic<C>
01628 make_mmx_analytic (const C& c) {
01629 return analytic<C> (c);
01630 }
01631
01632 template<typename C> inline analytic<C>
01633 make_mmx_analytic (const polynomial<C>& p) {
01634 return analytic<C> (p);
01635 }
01636
01637 template<typename C> inline analytic<C>
01638 make_mmx_analytic (const vector<C>& cs) {
01639 return analytic<C> (cs);
01640 }
01641
01642 #undef TMPL_DEF
01643 #undef TMPL
01644 #undef TMPLK
01645 #undef Format
01646 #undef Analytic
01647 #undef Analytic_rep
01648 #undef Recursive_analytic_rep
01649 #undef Assume
01650 #undef Scalar
01651 #undef Radius
01652 #undef Series
01653 #undef Series_rep
01654 #undef Recursive_series_rep
01655 #undef Polynomial
01656 #undef R
01657 }
01658 #endif // __MMX_ANALYTIC_HPP
