CppAD: A C++ Algorithmic Differentiation Package
20130102
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_POW_INCLUDED 00003 # define CPPAD_POW_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Eclipse 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 $icode%z% = pow(%x%, %y%)%$$ 00033 00034 $head See Also$$ 00035 $cref pow_int$$ 00036 00037 00038 $head Purpose$$ 00039 Determines the value of the power function which is defined by 00040 $latex \[ 00041 {\rm pow} (x, y) = x^y 00042 \] $$ 00043 This version of the $code pow$$ function may use 00044 logarithms and exponentiation to compute derivatives. 00045 This will not work if $icode x$$ is less than or equal zero. 00046 If the value of $icode y$$ is an integer, 00047 the $cref pow_int$$ function is used to compute this value 00048 using only multiplication (and division if $icode y$$ is negative). 00049 (This will work even if $icode x$$ is less than or equal zero.) 00050 00051 $head x$$ 00052 The argument $icode x$$ has one of the following prototypes 00053 $codei% 00054 const %Base%& %x% 00055 const AD<%Base%>& %x% 00056 const VecAD<%Base%>::reference& %x% 00057 %$$ 00058 00059 $head y$$ 00060 The argument $icode y$$ has one of the following prototypes 00061 $codei% 00062 const %Base%& %y% 00063 const AD<%Base%>& %y% 00064 const VecAD<%Base%>::reference& %y% 00065 %$$ 00066 00067 $head z$$ 00068 If both $icode x$$ and $icode y$$ are $icode Base$$ objects, 00069 the result $icode z$$ is also a $icode Base$$ object. 00070 Otherwise, it has prototype 00071 $codei% 00072 AD<%Base%> %z% 00073 %$$ 00074 00075 $head Operation Sequence$$ 00076 This is an AD of $icode Base$$ 00077 $cref/atomic operation/glossary/Operation/Atomic/$$ 00078 and hence is part of the current 00079 AD of $icode Base$$ 00080 $cref/operation sequence/glossary/Operation/Sequence/$$. 00081 00082 $head Example$$ 00083 $children% 00084 example/pow.cpp 00085 %$$ 00086 The file 00087 $cref pow.cpp$$ 00088 is an examples and tests of this function. 00089 It returns true if it succeeds and false otherwise. 00090 00091 $end 00092 ------------------------------------------------------------------------------- 00093 */ 00094 00095 // BEGIN CppAD namespace 00096 namespace CppAD { 00097 00098 // case where x and y are AD<Base> ----------------------------------------- 00099 template <class Base> AD<Base> 00100 pow(const AD<Base>& x, const AD<Base>& y) 00101 { 00102 // compute the Base part 00103 AD<Base> result; 00104 result.value_ = pow(x.value_, y.value_); 00105 CPPAD_ASSERT_UNKNOWN( Parameter(result) ); 00106 00107 // check if there is a recording in progress 00108 ADTape<Base>* tape = AD<Base>::tape_ptr(); 00109 if( tape == CPPAD_NULL ) 00110 return result; 00111 tape_id_t tape_id = tape->id_; 00112 00113 // tape_id cannot match the default value for tape_id_; i.e., 0 00114 CPPAD_ASSERT_UNKNOWN( tape_id > 0 ); 00115 bool var_x = x.tape_id_ == tape_id; 00116 bool var_y = y.tape_id_ == tape_id; 00117 00118 if( var_x ) 00119 { if( var_y ) 00120 { // result = variable^variable 00121 CPPAD_ASSERT_UNKNOWN( NumRes(PowvvOp) == 3 ); 00122 CPPAD_ASSERT_UNKNOWN( NumArg(PowvvOp) == 2 ); 00123 00124 // put operand addresses in tape 00125 tape->Rec_.PutArg(x.taddr_, y.taddr_); 00126 00127 // put operator in the tape 00128 result.taddr_ = tape->Rec_.PutOp(PowvvOp); 00129 00130 // make result a variable 00131 result.tape_id_ = tape_id; 00132 } 00133 else if( IdenticalZero( y.value_ ) ) 00134 { // result = variable^0 00135 } 00136 else 00137 { // result = variable^parameter 00138 CPPAD_ASSERT_UNKNOWN( NumRes(PowvpOp) == 3 ); 00139 CPPAD_ASSERT_UNKNOWN( NumArg(PowvpOp) == 2 ); 00140 00141 // put operand addresses in tape 00142 addr_t p = tape->Rec_.PutPar(y.value_); 00143 tape->Rec_.PutArg(x.taddr_, p); 00144 00145 // put operator in the tape 00146 result.taddr_ = tape->Rec_.PutOp(PowvpOp); 00147 00148 // make result a variable 00149 result.tape_id_ = tape_id; 00150 } 00151 } 00152 else if( var_y ) 00153 { if( IdenticalZero(x.value_) ) 00154 { // result = 0^variable 00155 } 00156 else 00157 { // result = variable^parameter 00158 CPPAD_ASSERT_UNKNOWN( NumRes(PowpvOp) == 3 ); 00159 CPPAD_ASSERT_UNKNOWN( NumArg(PowpvOp) == 2 ); 00160 00161 // put operand addresses in tape 00162 addr_t p = tape->Rec_.PutPar(x.value_); 00163 tape->Rec_.PutArg(p, y.taddr_); 00164 00165 // put operator in the tape 00166 result.taddr_ = tape->Rec_.PutOp(PowpvOp); 00167 00168 // make result a variable 00169 result.tape_id_ = tape_id; 00170 } 00171 } 00172 return result; 00173 } 00174 // ========================================================================= 00175 // Fold operations in same way as CPPAD_FOLD_AD_VALUED_BINARY_OPERATION(Op) 00176 // ------------------------------------------------------------------------- 00177 // Operations with VecAD_reference<Base> and AD<Base> only 00178 00179 template <class Base> AD<Base> 00180 pow(const AD<Base>& x, const VecAD_reference<Base>& y) 00181 { return pow(x, y.ADBase()); } 00182 00183 template <class Base> AD<Base> 00184 pow(const VecAD_reference<Base>& x, const VecAD_reference<Base>& y) 00185 { return pow(x.ADBase(), y.ADBase()); } 00186 00187 template <class Base> AD<Base> 00188 pow(const VecAD_reference<Base>& x, const AD<Base>& y) 00189 { return pow(x.ADBase(), y); } 00190 // ------------------------------------------------------------------------- 00191 // Operations with Base 00192 00193 template <class Base> AD<Base> 00194 pow(const Base& x, const AD<Base>& y) 00195 { return pow(AD<Base>(x), y); } 00196 00197 template <class Base> AD<Base> 00198 pow(const Base& x, const VecAD_reference<Base>& y) 00199 { return pow(AD<Base>(x), y.ADBase()); } 00200 00201 template <class Base> AD<Base> 00202 pow(const AD<Base>& x, const Base& y) 00203 { return pow(x, AD<Base>(y)); } 00204 00205 template <class Base> AD<Base> 00206 pow(const VecAD_reference<Base>& x, const Base& y) 00207 { return pow(x.ADBase(), AD<Base>(y)); } 00208 // ------------------------------------------------------------------------- 00209 // Operations with double 00210 00211 template <class Base> AD<Base> 00212 pow(const double& x, const AD<Base>& y) 00213 { return pow(AD<Base>(x), y); } 00214 00215 template <class Base> AD<Base> 00216 pow(const double& x, const VecAD_reference<Base>& y) 00217 { return pow(AD<Base>(x), y.ADBase()); } 00218 00219 template <class Base> AD<Base> 00220 pow(const AD<Base>& x, const double& y) 00221 { return pow(x, AD<Base>(y)); } 00222 00223 template <class Base> AD<Base> 00224 pow(const VecAD_reference<Base>& x, const double& y) 00225 { return pow(x.ADBase(), AD<Base>(y)); } 00226 // ------------------------------------------------------------------------- 00227 // Special case to avoid ambuigity when Base is double 00228 00229 inline AD<double> 00230 pow(const double& x, const AD<double>& y) 00231 { return pow(AD<double>(x), y); } 00232 00233 inline AD<double> 00234 pow(const double& x, const VecAD_reference<double>& y) 00235 { return pow(AD<double>(x), y.ADBase()); } 00236 00237 inline AD<double> 00238 pow(const AD<double>& x, const double& y) 00239 { return pow(x, AD<double>(y)); } 00240 00241 inline AD<double> 00242 pow(const VecAD_reference<double>& x, const double& y) 00243 { return pow(x.ADBase(), AD<double>(y)); } 00244 00245 // ========================================================================= 00246 // Fold operations for the cases where x is an int, 00247 // but let cppad/pow_int.hpp handle the cases where y is an int. 00248 // ------------------------------------------------------------------------- 00249 template <class Base> AD<Base> pow 00250 (const int& x, const VecAD_reference<Base>& y) 00251 { return pow(AD<Base>(x), y.ADBase()); } 00252 00253 template <class Base> AD<Base> pow 00254 (const int& x, const AD<Base>& y) 00255 { return pow(AD<Base>(x), y); } 00256 00257 } // END CppAD namespace 00258 00259 # endif