CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_LOG_OP_INCLUDED 00003 # define CPPAD_LOG_OP_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 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 CPPAD_BEGIN_NAMESPACE 00017 /*! 00018 \file log_op.hpp 00019 Forward and reverse mode calculations for z = log(x). 00020 */ 00021 00022 /*! 00023 Compute forward mode Taylor coefficient for result of op = LogOp. 00024 00025 The C++ source code corresponding to this operation is 00026 \verbatim 00027 z = log(x) 00028 \endverbatim 00029 00030 \copydetails forward_unary1_op 00031 */ 00032 template <class Base> 00033 inline void forward_log_op( 00034 size_t j , 00035 size_t i_z , 00036 size_t i_x , 00037 size_t nc_taylor , 00038 Base* taylor ) 00039 { 00040 size_t k; 00041 00042 // check assumptions 00043 CPPAD_ASSERT_UNKNOWN( NumArg(LogOp) == 1 ); 00044 CPPAD_ASSERT_UNKNOWN( NumRes(LogOp) == 1 ); 00045 CPPAD_ASSERT_UNKNOWN( i_x < i_z ); 00046 CPPAD_ASSERT_UNKNOWN( j < nc_taylor ); 00047 00048 // Taylor coefficients corresponding to argument and result 00049 Base* x = taylor + i_x * nc_taylor; 00050 Base* z = taylor + i_z * nc_taylor; 00051 00052 if( j == 0 ) 00053 z[0] = log( x[0] ); 00054 else if ( j == 1 ) 00055 z[1] = x[1] / x[0]; 00056 else 00057 { 00058 z[j] = -z[1] * x[j-1]; 00059 for(k = 2; k < j; k++) 00060 z[j] -= Base(k) * z[k] * x[j-k]; 00061 z[j] /= Base(j); 00062 z[j] += x[j]; 00063 z[j] /= x[0]; 00064 } 00065 } 00066 00067 /*! 00068 Compute zero order forward mode Taylor coefficient for result of op = LogOp. 00069 00070 The C++ source code corresponding to this operation is 00071 \verbatim 00072 z = log(x) 00073 \endverbatim 00074 00075 \copydetails forward_unary1_op_0 00076 */ 00077 template <class Base> 00078 inline void forward_log_op_0( 00079 size_t i_z , 00080 size_t i_x , 00081 size_t nc_taylor , 00082 Base* taylor ) 00083 { 00084 00085 // check assumptions 00086 CPPAD_ASSERT_UNKNOWN( NumArg(LogOp) == 1 ); 00087 CPPAD_ASSERT_UNKNOWN( NumRes(LogOp) == 1 ); 00088 CPPAD_ASSERT_UNKNOWN( i_x < i_z ); 00089 CPPAD_ASSERT_UNKNOWN( 0 < nc_taylor ); 00090 00091 // Taylor coefficients corresponding to argument and result 00092 Base* x = taylor + i_x * nc_taylor; 00093 Base* z = taylor + i_z * nc_taylor; 00094 00095 z[0] = log( x[0] ); 00096 } 00097 00098 /*! 00099 Compute reverse mode partial derivatives for result of op = LogOp. 00100 00101 The C++ source code corresponding to this operation is 00102 \verbatim 00103 z = log(x) 00104 \endverbatim 00105 00106 \copydetails reverse_unary1_op 00107 */ 00108 00109 template <class Base> 00110 inline void reverse_log_op( 00111 size_t d , 00112 size_t i_z , 00113 size_t i_x , 00114 size_t nc_taylor , 00115 const Base* taylor , 00116 size_t nc_partial , 00117 Base* partial ) 00118 { size_t j, k; 00119 00120 // check assumptions 00121 CPPAD_ASSERT_UNKNOWN( NumArg(LogOp) == 1 ); 00122 CPPAD_ASSERT_UNKNOWN( NumRes(LogOp) == 1 ); 00123 CPPAD_ASSERT_UNKNOWN( i_x < i_z ); 00124 CPPAD_ASSERT_UNKNOWN( d < nc_taylor ); 00125 CPPAD_ASSERT_UNKNOWN( d < nc_partial ); 00126 00127 // Taylor coefficients and partials corresponding to argument 00128 const Base* x = taylor + i_x * nc_taylor; 00129 Base* px = partial + i_x * nc_partial; 00130 00131 // Taylor coefficients and partials corresponding to result 00132 const Base* z = taylor + i_z * nc_taylor; 00133 Base* pz = partial + i_z * nc_partial; 00134 00135 j = d; 00136 while(j) 00137 { // scale partial w.r.t z[j] 00138 pz[j] /= x[0]; 00139 00140 px[0] -= pz[j] * z[j]; 00141 px[j] += pz[j]; 00142 00143 // further scale partial w.r.t. z[j] 00144 pz[j] /= Base(j); 00145 00146 for(k = 1; k < j; k++) 00147 { pz[k] -= pz[j] * Base(k) * x[j-k]; 00148 px[j-k] -= pz[j] * Base(k) * z[k]; 00149 } 00150 --j; 00151 } 00152 px[0] += pz[0] / x[0]; 00153 } 00154 00155 CPPAD_END_NAMESPACE 00156 # endif