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