CppAD: A C++ Algorithmic Differentiation Package 20110419
pow.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_POW_INCLUDED
00003 # define CPPAD_POW_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin pow$$
00018 $spell
00019         Vec
00020         std
00021         namespace
00022         CppAD
00023         const
00024 $$
00025 
00026 $index pow, AD$$
00027 $index exponent, AD function$$
00028 
00029 $section The AD Power Function$$
00030 
00031 $head Syntax$$
00032 $syntax%%z% = pow(%x%, %y%)%$$
00033 
00034 
00035 $head Purpose$$
00036 Determines the value of the power function which is defined by
00037 $latex \[
00038         {\rm pow} (x, y) = x^y
00039 \] $$
00040 This version of the $code pow$$ function may use
00041 logarithms and exponentiation to compute derivatives.
00042 This will not work if $italic x$$ is less than or equal zero.
00043 If the value of $italic y$$ is an integer, 
00044 the $cref/pow_int/$$ function is used to compute this value 
00045 using only multiplication (and division if $italic y$$ is negative). 
00046 (This will work even if $italic x$$ is less than or equal zero.)
00047 
00048 $head x$$
00049 The argument $italic x$$ has the following prototype
00050 $syntax%
00051         const %Type% &%x%
00052 %$$
00053 where $italic Type$$ is
00054 $syntax%VecAD<%Base%>::reference%$$,
00055 $syntax%AD<%Base%>%$$,
00056 $syntax%%Base%$$,
00057 $syntax%double%$$,
00058 or
00059 $syntax%int%$$.
00060 
00061 $head y$$
00062 The argument $italic y$$ has the following prototype
00063 $syntax%
00064         const %Type% &%y%
00065 %$$
00066 where $italic Type$$ is
00067 $syntax%VecAD<%Base%>::reference%$$,
00068 $syntax%AD<%Base%>%$$,
00069 $syntax%%Base%$$,
00070 $syntax%double%$$,
00071 or
00072 $syntax%int%$$.
00073 
00074 $head z$$
00075 The result $italic z$$ has prototype
00076 $syntax%
00077         AD<%Base%> %z%
00078 %$$
00079 
00080 $head Standard Types$$
00081 A definition for the $code pow$$ function is included
00082 in the CppAD namespace for the case where both $italic x$$
00083 and $italic y$$ have the same type and that type is
00084 $code float$$ or $code double$$.
00085  
00086 $head Operation Sequence$$
00087 This is an AD of $italic Base$$
00088 $xref/glossary/Operation/Atomic/atomic operation/1/$$
00089 and hence is part of the current
00090 AD of $italic Base$$
00091 $xref/glossary/Operation/Sequence/operation sequence/1/$$.
00092 
00093 $head Example$$
00094 $children%
00095         example/pow.cpp%
00096         example/pow_int.cpp
00097 %$$
00098 The files
00099 $xref/Pow.cpp/$$, $cref/pow_int.cpp/$$
00100 contain an examples and tests of this function.   
00101 They returns true if they succeed and false otherwise.
00102 
00103 $end
00104 -------------------------------------------------------------------------------
00105 */
00106 
00107 //  BEGIN CppAD namespace
00108 namespace CppAD {
00109  
00110 // copy of standard functions in CppAD namespace
00111 inline float pow(const float &x, const float &y)
00112 { return std::pow(x, y); }
00113 
00114 inline double pow(const double &x, const double &y)
00115 { return std::pow(x, y); }
00116 
00117 // case where x and y are AD<Base> -----------------------------------------
00118 template <class Base> AD<Base> 
00119 pow(const AD<Base> &x, const AD<Base> &y)
00120 {       ADTape<Base> *tape = AD<Base>::tape_ptr();
00121         size_t tape_id = 0;
00122         if( tape != CPPAD_NULL )
00123                 tape_id = tape->id_;
00124 
00125         // id_ setting for parameters cannot match 0
00126         bool var_x = x.id_  == tape_id;
00127         bool var_y = y.id_ == tape_id;
00128         CPPAD_ASSERT_KNOWN(
00129                 Parameter(x) || var_x ,
00130                 "pow: first argument is a variable for a different thread"
00131         );
00132         CPPAD_ASSERT_KNOWN(
00133                 Parameter(y) || var_y ,
00134                 "pow: second argument is a variable for a different thread"
00135         );
00136 
00137         AD<Base> result;
00138         result.value_  = pow(x.value_, y.value_);
00139         CPPAD_ASSERT_UNKNOWN( Parameter(result) );
00140 
00141         if( var_x )
00142         {       if( var_y )
00143                 {       // result = variable^variable
00144                         CPPAD_ASSERT_UNKNOWN( NumRes(PowvvOp) == 3 );
00145                         CPPAD_ASSERT_UNKNOWN( NumArg(PowvvOp) == 2 );
00146 
00147                         // put operand addresses in tape
00148                         tape->Rec_.PutArg(x.taddr_, y.taddr_);
00149 
00150                         // put operator in the tape
00151                         result.taddr_ = tape->Rec_.PutOp(PowvvOp);
00152 
00153                         // make result a variable
00154                         result.id_ = tape_id;
00155                 }
00156                 else if( IdenticalZero( y.value_ ) )
00157                 {       // result = variable^0
00158                 }
00159                 else
00160                 {       // result = variable^parameter 
00161                         CPPAD_ASSERT_UNKNOWN( NumRes(PowvpOp) == 3 );
00162                         CPPAD_ASSERT_UNKNOWN( NumArg(PowvpOp) == 2 );
00163 
00164                         // put operand addresses in tape
00165                         size_t p = tape->Rec_.PutPar(y.value_);
00166                         tape->Rec_.PutArg(x.taddr_, p);
00167 
00168                         // put operator in the tape
00169                         result.taddr_ = tape->Rec_.PutOp(PowvpOp);
00170 
00171                         // make result a variable
00172                         result.id_ = tape_id;
00173                 }
00174         }
00175         else if( var_y )
00176         {       if( IdenticalZero(x.value_) )
00177                 {       // result = 0^variable 
00178                 }
00179                 else
00180                 {       // result = variable^parameter 
00181                         CPPAD_ASSERT_UNKNOWN( NumRes(PowpvOp) == 3 );
00182                         CPPAD_ASSERT_UNKNOWN( NumArg(PowpvOp) == 2 );
00183 
00184                         // put operand addresses in tape
00185                         size_t p = tape->Rec_.PutPar(x.value_);
00186                         tape->Rec_.PutArg(p, y.taddr_);
00187 
00188                         // put operator in the tape
00189                         result.taddr_ = tape->Rec_.PutOp(PowpvOp);
00190 
00191                         // make result a variable
00192                         result.id_ = tape_id;
00193                 }
00194         }
00195         return result;
00196 }
00197 // =========================================================================
00198 // Fold operations in same way as CPPAD_FOLD_AD_VALUED_BINARY_OPERATION(Op)
00199 // -------------------------------------------------------------------------
00200 // Operations with VecAD_reference<Base> and AD<Base> only
00201 
00202 template <class Base> AD<Base>
00203 pow(const AD<Base> &x, const VecAD_reference<Base> &y)
00204 {       return pow(x, y.ADBase()); }
00205 
00206 template <class Base> AD<Base> 
00207 pow(const VecAD_reference<Base> &x, const VecAD_reference<Base> &y) 
00208 {       return pow(x.ADBase(), y.ADBase()); }
00209 
00210 template <class Base> AD<Base>
00211 pow(const VecAD_reference<Base> &x, const AD<Base> &y)
00212 {       return pow(x.ADBase(), y); }
00213 // -------------------------------------------------------------------------
00214 // Operations with Base
00215 
00216 template <class Base> AD<Base>
00217 pow(const Base &x, const AD<Base> &y)
00218 {       return pow(AD<Base>(x), y); }
00219 
00220 template <class Base> AD<Base>
00221 pow(const Base &x, const VecAD_reference<Base> &y)
00222 {       return pow(AD<Base>(x), y.ADBase()); }
00223 
00224 template <class Base> AD<Base>
00225 pow(const AD<Base> &x, const Base &y)
00226 {       return pow(x, AD<Base>(y)); }
00227 
00228 template <class Base> AD<Base>
00229 pow(const VecAD_reference<Base> &x, const Base &y)
00230 {       return pow(x.ADBase(), AD<Base>(y)); }
00231 // -------------------------------------------------------------------------
00232 // Operations with double
00233 
00234 template <class Base> AD<Base>
00235 pow(const double &x, const AD<Base> &y)
00236 {       return pow(AD<Base>(x), y); }
00237 
00238 template <class Base> AD<Base>
00239 pow(const double &x, const VecAD_reference<Base> &y)
00240 {       return pow(AD<Base>(x), y.ADBase()); }
00241 
00242 template <class Base> AD<Base>
00243 pow(const AD<Base> &x, const double &y)
00244 {       return pow(x, AD<Base>(y)); }
00245 
00246 template <class Base> AD<Base>
00247 pow(const VecAD_reference<Base> &x, const double &y)
00248 {       return pow(x.ADBase(), AD<Base>(y)); }
00249 // -------------------------------------------------------------------------
00250 // Special case to avoid ambuigity when Base is double
00251 
00252 inline AD<double>
00253 pow(const double &x, const AD<double> &y)
00254 {       return pow(AD<double>(x), y); }
00255 
00256 inline AD<double>
00257 pow(const double &x, const VecAD_reference<double> &y)
00258 {       return pow(AD<double>(x), y.ADBase()); }
00259 
00260 inline AD<double>
00261 pow(const AD<double> &x, const double &y)
00262 {       return pow(x, AD<double>(y)); }
00263 
00264 inline AD<double>
00265 pow(const VecAD_reference<double> &x, const double &y)
00266 {       return pow(x.ADBase(), AD<double>(y)); }
00267 
00268 // =========================================================================
00269 // Fold operations for the cases where x is an int, 
00270 // but let cppad/pow_int.hpp handle the cases where y is an int.
00271 // -------------------------------------------------------------------------
00272 template <class Base> AD<Base> pow
00273 (const int &x, const VecAD_reference<Base> &y)
00274 {       return pow(AD<Base>(x), y.ADBase()); }
00275 
00276 template <class Base> AD<Base> pow
00277 (const int &x, const AD<Base> &y)
00278 {       return pow(AD<Base>(x), y); }
00279 
00280 } // END CppAD namespace
00281 
00282 # endif