Back to index

texmacs  1.0.7.15
function_extra.hpp
Go to the documentation of this file.
00001 
00002 /******************************************************************************
00003 * MODULE     : function_extra.hpp
00004 * DESCRIPTION: Extra routines for functions from F into T
00005 * COPYRIGHT  : (C) 2006  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license version 3 or later.
00008 * It comes WITHOUT ANY WARRANTY WHATSOEVER. For details, see the file LICENSE
00009 * in the root directory or <http://www.gnu.org/licenses/gpl-3.0.html>.
00010 ******************************************************************************/
00011 
00012 #ifndef FUNCTION_EXTRA_HPP
00013 #define FUNCTION_EXTRA_HPP
00014 #include "function.hpp"
00015 #include "vector.hpp"
00016 #include "polynomial.hpp"
00017 #define TMPL template<typename F, typename T>
00018 #define V typename properties<F>::index_type
00019 #define C typename properties<F>::scalar_type
00020 
00021 #define MATH_INFINITY HUGE_VAL
00022 
00023 /******************************************************************************
00024 * Piecewise functions
00025 ******************************************************************************/
00026 
00027 TMPL function<F,T>
00028 pw_function (vector<function<F,T> > v, bool jump_flag= false);
00029 
00030 TMPL
00031 class pw_function_rep: public function_rep<F,T> {
00032   vector<function<F,T> > v;
00033   bool jump_flag;
00034 public:
00035   inline pw_function_rep (vector<function<F,T> > v2, bool jump_flag2):
00036     v (v2), jump_flag (jump_flag2) {}
00037   int which (F x) {
00038     if (x <= 0) return 0;
00039     if (x >= 1) return N(v)-1;
00040     return (int) floor (N(v) * x);
00041   }
00042   F position (F x) { return N(v) * x - which (x); }
00043   T apply (F x) { return v[which (x)] (position (x)); }
00044   ball<T> apply (ball<F> x) {
00045     int w1= which (lower (x));
00046     int w2= which (upper (x));
00047     if (w1 == w2) return v[w1] (ball<F> (N(v)) * x - ball<F> (w1));
00048     if (jump_flag) return ball<T> (T(0), MATH_INFINITY);
00049   }
00050   function<F,T> derive (V var) {
00051     function<F,T> dummy;
00052     vector<function<F,T> > w (dummy, N(v));
00053     for (int i=0; i<N(v); i++)
00054       w[i]= function<F,T> (N(v)) * v[i]->derive (var);
00055     bool flag= jump_flag;
00056     for (int i=0; i<N(v)-1; i++)
00057       flag= flag || (v[i](F(1)) != v[i+1](F(0)));
00058     return pw_function (w, flag);
00059   }
00060   tree expression () {
00061     int i, n= N(v);
00062     array<tree> a (n+1);
00063     for (i=0; i<n; i++)
00064       a[i]= as_tree (v[i]);
00065     a[n]= (jump_flag? tree ("true"): tree ("false"));
00066     return compound ("piecewise", a);
00067   }
00068 };
00069 
00070 TMPL inline function<F,T>
00071 pw_function (vector<function<F,T> > v, bool jump_flag) {
00072   return tm_new<pw_function_rep<F,T> > (v, jump_flag);
00073 }
00074 
00075 /******************************************************************************
00076 * Vector functions
00077 ******************************************************************************/
00078 
00079 template<typename T> ball<vector<T> >
00080 as_ball (vector<ball<T> > v) {
00081   typedef typename properties<T>::norm_type radius_type;
00082   int i, n= N(v);
00083   vector<T> c (T(0), n);
00084   vector<radius_type> r (radius_type(0), n);
00085   for (i=0; i<n; i++) {
00086     c[i]= center (v[i]);
00087     r[i]= radius (v[i]);
00088   }
00089   return ball<vector<T> > (c, norm (r));
00090 }
00091 
00092 TMPL function<F,vector<T> > vector_function (vector<function<F,T> > v);
00093 
00094 TMPL
00095 class vector_function_rep: public function_rep<F,vector<T> > {
00096   vector<function<F,T> > v;
00097 public:
00098   inline vector_function_rep (vector<function<F,T> > v2): v (v2) {}
00099   vector<T> apply (F x) {
00100     int i, n= N(v);
00101     T* a= tm_new_array<T> (n);
00102     for (i=0; i<n; i++)
00103       a[i]= v[i] (x);
00104     return vector<T> (a, n);
00105   }
00106   ball<vector<T> > apply (ball<F> x) {
00107     int i, n= N(v);
00108     ball<T>* a= tm_new_array<ball<T> > (n);
00109     for (i=0; i<n; i++)
00110       a[i]= v[i] (x);
00111     return as_ball (vector<ball<T> > (a, n));
00112   }
00113   function<F,vector<T> > derive (V var) {
00114     int i, n= N(v);
00115     vector<function<F,T> > w;
00116     for (i=0; i<n; i++) w[i]= v[i]->derive (var);
00117     return vector_function (w);
00118   }
00119   tree expression () {
00120     return compound ("as_function", as_tree (v));
00121   }
00122 };
00123 
00124 TMPL inline function<F,vector<T> >
00125 vector_function (vector<function<F,T> > v) {
00126   return tm_new<vector_function_rep<F,T> > (v);
00127 }
00128 
00129 /******************************************************************************
00130 * Polynomial functions
00131 ******************************************************************************/
00132 
00133 template<typename T> function<T,T> polynomial_function (polynomial<T> p);
00134 
00135 template<typename T>
00136 class polynomial_function_rep: public function_rep<T,T> {
00137   polynomial<T> p;
00138 public:
00139   inline polynomial_function_rep (polynomial<T> p2): p (p2) {}
00140   inline T apply (T x) { return p (x); }
00141   ball<T> apply (ball<T> x) {
00142     int i, n= N(p);
00143     if (n == 0) return 0;
00144     T* a= A(p);
00145     ball<T> sum= a[n-1];
00146     for (i=n-2; i>=0; i--)
00147       sum = x * sum + ball<T> (a[i]);
00148     return sum; }
00149   function<T,T> derive (typename properties<T>::index_type var) {
00150     return polynomial_function<T> (::derive (p)); }
00151   tree expression () {
00152     return compound ("as_function", as_tree (p)); }
00153 };
00154 
00155 template<typename T> inline function<T,T>
00156 polynomial_function (polynomial<T> p) {
00157   return tm_new<polynomial_function_rep<T> > (p);
00158 }
00159 
00160 #undef TMPL
00161 #undef V
00162 #undef C
00163 #endif // FUNCTION_EXTRA_H