00001 # ifndef CPPAD_EXP_EPS_INCLUDED 00002 # define CPPAD_EXP_EPS_INCLUDED 00003 /* -------------------------------------------------------------------------- 00004 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell 00005 00006 CppAD is distributed under multiple licenses. This distribution is under 00007 the terms of the 00008 Common Public License Version 1.0. 00009 00010 A copy of this license is included in the COPYING file of this distribution. 00011 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00012 -------------------------------------------------------------------------- */ 00013 /* 00014 $begin exp_eps$$ 00015 $spell 00016 cppad-%yy%-%mm%-%dd% 00017 hpp 00018 Apx 00019 cpp 00020 const 00021 exp_eps 00022 bool 00023 $$ 00024 00025 $section An Epsilon Accurate Exponential Approximation$$ 00026 00027 $index exp_eps$$ 00028 $index example, algorithm$$ 00029 $index algorithm, example$$ 00030 $index exp, example$$ 00031 00032 $head Syntax$$ 00033 $syntax%# include "exp_eps.hpp"%$$ 00034 $pre 00035 $$ 00036 $syntax%%y% = exp_eps(%x%, %epsilon%)%$$ 00037 00038 00039 $head Purpose$$ 00040 This is a an example algorithm that is used to demonstrate 00041 how Algorithmic Differentiation works with loops and 00042 boolean decision variables 00043 (see $cref/exp_2/$$ for a simpler example). 00044 00045 $head Mathematical Function$$ 00046 The exponential function can be defined by 00047 $latex \[ 00048 \exp (x) = 1 + x^1 / 1 ! + x^2 / 2 ! + \cdots 00049 \] $$ 00050 We define $latex k ( x, \varepsilon ) $$ as the smallest 00051 non-negative integer such that $latex \varepsilon \geq x^k / k !$$; i.e., 00052 $latex \[ 00053 k( x, \varepsilon ) = 00054 \min \{ k \in {\rm Z}_+ \; | \; \varepsilon \geq x^k / k ! \} 00055 \] $$ 00056 The mathematical form for our approximation of the exponential function is 00057 $latex \[ 00058 \begin{array}{rcl} 00059 {\rm exp\_eps} (x , \varepsilon ) & = & \left\{ 00060 \begin{array}{ll} 00061 \frac{1}{ {\rm exp\_eps} (-x , \varepsilon ) } 00062 & {\rm if} \; x < 0 00063 \\ 00064 1 + x^1 / 1 ! + \cdots + x^{k( x, \varepsilon)} / k( x, \varepsilon ) ! 00065 & {\rm otherwise} 00066 \end{array} 00067 \right. 00068 \end{array} 00069 \] $$ 00070 00071 00072 $head include$$ 00073 The include command in the syntax is relative to 00074 $syntax% 00075 cppad-%yy%-%mm%-%dd%/introduction/exp_apx 00076 %$$ 00077 where $syntax%cppad-%yy%-%mm%-%dd%$$ is the distribution directory 00078 created during the beginning steps of the 00079 $cref%installation%Install%$$ of CppAD. 00080 00081 00082 $head x$$ 00083 The argument $italic x$$ has prototype 00084 $syntax% 00085 const %Type% &%x% 00086 %$$ 00087 (see $italic Type$$ below). 00088 It specifies the point at which to evaluate the 00089 approximation for the exponential function. 00090 00091 $head epsilon$$ 00092 The argument $italic epsilon$$ has prototype 00093 $syntax% 00094 const %Type% &%epsilon% 00095 %$$ 00096 It specifies the accuracy with which 00097 to approximate the exponential function value; i.e., 00098 it is the value of $latex \varepsilon$$ in the 00099 exponential function approximation defined above. 00100 00101 $head y$$ 00102 The result $italic y$$ has prototype 00103 $syntax% 00104 %Type% %y% 00105 %$$ 00106 It is the value of the exponential function 00107 approximation defined above. 00108 00109 $head Type$$ 00110 If $italic u$$ and $italic v$$ are $italic Type$$ objects and $italic i$$ 00111 is an $code int$$: 00112 00113 $table 00114 $bold Operation$$ $cnext $bold Result Type$$ $cnext $bold Description$$ 00115 $rnext 00116 $syntax%%Type%(%i%)%$$ 00117 $cnext $italic Type$$ 00118 $cnext object with value equal to $italic i$$ 00119 $rnext 00120 $syntax%%u% > %v%$$ 00121 $cnext $code bool$$ 00122 $cnext true, 00123 if $italic u$$ greater than $italic v$$, an false otherwise 00124 $rnext 00125 $syntax%%u% = %v%$$ 00126 $cnext $italic Type$$ 00127 $cnext new $italic u$$ (and result) is value of $italic v$$ 00128 $rnext 00129 $syntax%%u% * %v%$$ 00130 $cnext $italic Type$$ 00131 $cnext result is value of $latex u * v$$ 00132 $rnext 00133 $syntax%%u% / %v%$$ 00134 $cnext $italic Type$$ 00135 $cnext result is value of $latex u / v$$ 00136 $rnext 00137 $syntax%%u% + %v%$$ 00138 $cnext $italic Type$$ 00139 $cnext result is value of $latex u + v$$ 00140 $rnext 00141 $syntax%-%u%$$ 00142 $cnext $italic Type$$ 00143 $cnext result is value of $latex - u$$ 00144 $tend 00145 00146 $children% 00147 introduction/exp_apx/exp_eps.omh% 00148 introduction/exp_apx/exp_eps_cppad.cpp 00149 %$$ 00150 00151 $head Implementation$$ 00152 The file $xref/exp_eps.hpp/$$ 00153 contains a C++ implementation of this function. 00154 00155 $head Test$$ 00156 The file $xref/exp_eps.cpp/$$ 00157 contains a test of this implementation. 00158 It returns true for success and false for failure. 00159 00160 $head Exercises$$ 00161 $list number$$ 00162 Using the definition of $latex k( x, \varepsilon )$$ above, 00163 what is the value of 00164 $latex k(.5, 1)$$, $latex k(.5, .1)$$, and $latex k(.5, .01)$$ ? 00165 $lnext 00166 Suppose that we make the following call to $code exp_eps$$: 00167 $codep 00168 double x = 1.; 00169 double epsilon = .01; 00170 double y = exp_eps(x, epsilon); 00171 $$ 00172 What is the value assigned to 00173 $code k$$, $code temp$$, $code term$$, and $code sum$$ 00174 the first time through the $code while$$ loop in $cref/exp_eps.hpp/$$ ? 00175 $lnext 00176 Continuing the previous exercise, 00177 what is the value assigned to 00178 $code k$$, $code temp$$, $code term$$, and $code sum$$ 00179 the second time through the $code while$$ loop in $cref/exp_eps.hpp/$$ ? 00180 $lend 00181 00182 00183 $end 00184 ----------------------------------------------------------------------------- 00185 */ 00186 // BEGIN PROGRAM 00187 template <class Type> 00188 Type exp_eps(const Type &x, const Type &epsilon) 00189 { // abs_x = |x| 00190 Type abs_x = x; 00191 if( Type(0) > x ) 00192 abs_x = - x; 00193 // initialize 00194 int k = 0; // initial order 00195 Type term = 1.; // term = |x|^k / k ! 00196 Type sum = term; // initial sum 00197 while(term > epsilon) 00198 { k = k + 1; // order for next term 00199 Type temp = term * abs_x; // term = |x|^k / (k-1)! 00200 term = temp / Type(k); // term = |x|^k / k ! 00201 sum = sum + term; // sum = 1 + ... + |x|^k / k ! 00202 } 00203 // In the case where x is negative, use exp(x) = 1 / exp(-|x|) 00204 if( Type(0) > x ) 00205 sum = Type(1) / sum; 00206 return sum; 00207 } 00208 // END PROGRAM 00209 00210 # endif