# ifndef CPPAD_POW_INCLUDED # define CPPAD_POW_INCLUDED /* -------------------------------------------------------------------------- CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell CppAD is distributed under multiple licenses. This distribution is under the terms of the Common Public License Version 1.0. A copy of this license is included in the COPYING file of this distribution. Please visit http://www.coin-or.org/CppAD/ for information on other licenses. -------------------------------------------------------------------------- */ /* $begin pow$$ $spell Vec std namespace CppAD const $$ $index pow, AD$$ $index exponent, AD function$$ $section The AD Power Function$$ $head Syntax$$ $syntax%%z% = pow(%x%, %y%)%$$ $head Purpose$$ Determines the value of the power function which is defined by $latex \[ {\rm pow} (x, y) = x^y \] $$ This version of the $code pow$$ function may use logarithms and exponentiation to compute derivatives. This will not work if $italic x$$ is less than or equal zero. If the value of $italic y$$ is an integer, the $cref/pow_int/$$ function is used to compute this value using only multiplication (and division if $italic y$$ is negative). (This will work even if $italic x$$ is less than or equal zero.) $head x$$ The argument $italic x$$ has the following prototype $syntax% const %Type% &%x% %$$ where $italic Type$$ is $syntax%VecAD<%Base%>::reference%$$, $syntax%AD<%Base%>%$$, $syntax%%Base%$$, $syntax%double%$$, or $syntax%int%$$. $head y$$ The argument $italic y$$ has the following prototype $syntax% const %Type% &%y% %$$ where $italic Type$$ is $syntax%VecAD<%Base%>::reference%$$, $syntax%AD<%Base%>%$$, $syntax%%Base%$$, $syntax%double%$$, or $syntax%int%$$. $head z$$ The result $italic z$$ has prototype $syntax% AD<%Base%> %z% %$$ $head Standard Types$$ A definition for the $code pow$$ function is included in the CppAD namespace for the case where both $italic x$$ and $italic y$$ have the same type and that type is $code float$$ or $code double$$. $head Operation Sequence$$ This is an AD of $italic Base$$ $xref/glossary/Operation/Atomic/atomic operation/1/$$ and hence is part of the current AD of $italic Base$$ $xref/glossary/Operation/Sequence/operation sequence/1/$$. $head Example$$ $children% example/pow.cpp% example/pow_int.cpp %$$ The files $xref/Pow.cpp/$$, $cref/pow_int.cpp/$$ contain an examples and tests of this function. They returns true if they succeed and false otherwise. $end ------------------------------------------------------------------------------- */ // BEGIN CppAD namespace namespace CppAD { // copy of standard functions in CppAD namespace inline float pow(const float &x, const float &y) { return std::pow(x, y); } inline double pow(const double &x, const double &y) { return std::pow(x, y); } // case where x and y are AD ----------------------------------------- template AD pow(const AD &x, const AD &y) { ADTape *tape = AD::tape_ptr(); bool var_x, var_y; # ifdef NDEBUG if( tape == CPPAD_NULL ) { var_x = false; var_y = false; } else { var_x = x.id_ == tape->id_; var_y = y.id_ == tape->id_; } # else var_x = Variable(x); var_y = Variable(y); CPPAD_ASSERT_KNOWN( (! var_x) || x.id_ == tape->id_ , "pow first operand is a variable for a different thread" ); CPPAD_ASSERT_KNOWN( (! var_y) || y.id_ == tape->id_ , "pow second operand is a variable for a different thread" ); # endif AD result; result.value_ = pow(x.value_, y.value_); CPPAD_ASSERT_UNKNOWN( Parameter(result) ); if( var_x ) { if( var_y ) { // result = variable^variable CPPAD_ASSERT_UNKNOWN( NumVar(PowvvOp) == 3 ); CPPAD_ASSERT_UNKNOWN( NumInd(PowvvOp) == 2 ); // put operand addresses in tape tape->Rec_.PutInd(x.taddr_, y.taddr_); // put operator in the tape result.taddr_ = tape->Rec_.PutOp(PowvvOp); // make result a variable result.id_ = tape->id_; } else if( IdenticalZero( y.value_ ) ) { // result = variable^0 } else { // result = variable^parameter CPPAD_ASSERT_UNKNOWN( NumVar(PowvpOp) == 3 ); CPPAD_ASSERT_UNKNOWN( NumInd(PowvpOp) == 2 ); // put operand addresses in tape size_t p = tape->Rec_.PutPar(y.value_); tape->Rec_.PutInd(x.taddr_, p); // put operator in the tape result.taddr_ = tape->Rec_.PutOp(PowvpOp); // make result a variable result.id_ = tape->id_; } } else if( var_y ) { if( IdenticalZero(x.value_) ) { // result = 0^variable } else { // result = variable^parameter CPPAD_ASSERT_UNKNOWN( NumVar(PowpvOp) == 3 ); CPPAD_ASSERT_UNKNOWN( NumInd(PowpvOp) == 2 ); // put operand addresses in tape size_t p = tape->Rec_.PutPar(x.value_); tape->Rec_.PutInd(p, y.taddr_); // put operator in the tape result.taddr_ = tape->Rec_.PutOp(PowpvOp); // make result a variable result.id_ = tape->id_; } } return result; } // ========================================================================= // Fold operations in same way as CPPAD_FOLD_AD_VALUED_BINARY_OPERATION(Op) // ------------------------------------------------------------------------- // Operations with VecAD_reference and AD only template AD pow(const AD &x, const VecAD_reference &y) { return pow(x, y.ADBase()); } template AD pow(const VecAD_reference &x, const VecAD_reference &y) { return pow(x.ADBase(), y.ADBase()); } template AD pow(const VecAD_reference &x, const AD &y) { return pow(x.ADBase(), y); } // ------------------------------------------------------------------------- // Operations with Base template AD pow(const Base &x, const AD &y) { return pow(AD(x), y); } template AD pow(const Base &x, const VecAD_reference &y) { return pow(AD(x), y.ADBase()); } template AD pow(const AD &x, const Base &y) { return pow(x, AD(y)); } template AD pow(const VecAD_reference &x, const Base &y) { return pow(x.ADBase(), AD(y)); } // ------------------------------------------------------------------------- // Operations with double template AD pow(const double &x, const AD &y) { return pow(AD(x), y); } template AD pow(const double &x, const VecAD_reference &y) { return pow(AD(x), y.ADBase()); } template AD pow(const AD &x, const double &y) { return pow(x, AD(y)); } template AD pow(const VecAD_reference &x, const double &y) { return pow(x.ADBase(), AD(y)); } // ------------------------------------------------------------------------- // Special case to avoid ambuigity when Base is double inline AD pow(const double &x, const AD &y) { return pow(AD(x), y); } inline AD pow(const double &x, const VecAD_reference &y) { return pow(AD(x), y.ADBase()); } inline AD pow(const AD &x, const double &y) { return pow(x, AD(y)); } inline AD pow(const VecAD_reference &x, const double &y) { return pow(x.ADBase(), AD(y)); } // ========================================================================= // Fold operations for the cases where x is an int, // but let cppad/pow_int.hpp handle the cases where y is an int. // ------------------------------------------------------------------------- template AD pow (int x, const VecAD_reference &y) { return pow(AD(x), y.ADBase()); } template AD pow (int x, const AD &y) { return pow(AD(x), y); } } // END CppAD namespace # endif