CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_ATAN_OP_INCLUDED 00003 # define CPPAD_ATAN_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 00017 CPPAD_BEGIN_NAMESPACE 00018 /*! 00019 \file atan_op.hpp 00020 Forward and reverse mode calculations for z = atan(x). 00021 */ 00022 00023 00024 /*! 00025 Forward mode Taylor coefficient for result of op = AtanOp. 00026 00027 The C++ source code corresponding to this operation is 00028 \verbatim 00029 z = atan(x) 00030 \endverbatim 00031 The auxillary result is 00032 \verbatim 00033 y = 1 + x * x 00034 \endverbatim 00035 The value of y, and its derivatives, are computed along with the value 00036 and derivatives of z. 00037 00038 \copydetails forward_unary2_op 00039 */ 00040 template <class Base> 00041 inline void forward_atan_op( 00042 size_t j , 00043 size_t i_z , 00044 size_t i_x , 00045 size_t nc_taylor , 00046 Base* taylor ) 00047 { 00048 // check assumptions 00049 CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 ); 00050 CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 ); 00051 CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z ); 00052 CPPAD_ASSERT_UNKNOWN( j < nc_taylor ); 00053 00054 // Taylor coefficients corresponding to argument and result 00055 Base* x = taylor + i_x * nc_taylor; 00056 Base* z = taylor + i_z * nc_taylor; 00057 Base* b = z - nc_taylor; // called y in documentation 00058 00059 size_t k; 00060 if( j == 0 ) 00061 { z[j] = atan( x[0] ); 00062 b[j] = Base(1) + x[0] * x[0]; 00063 } 00064 else 00065 { 00066 b[j] = Base(2) * x[0] * x[j]; 00067 z[j] = Base(0); 00068 for(k = 1; k < j; k++) 00069 { b[j] += x[k] * x[j-k]; 00070 z[j] -= Base(k) * z[k] * b[j-k]; 00071 } 00072 z[j] /= Base(j); 00073 z[j] += x[j]; 00074 z[j] /= b[0]; 00075 } 00076 } 00077 00078 /*! 00079 Zero order forward mode Taylor coefficient for result of op = AtanOp. 00080 00081 The C++ source code corresponding to this operation is 00082 \verbatim 00083 z = atan(x) 00084 \endverbatim 00085 The auxillary result is 00086 \verbatim 00087 y = 1 + x * x 00088 \endverbatim 00089 The value of y is computed along with the value of z. 00090 00091 \copydetails forward_unary2_op_0 00092 */ 00093 template <class Base> 00094 inline void forward_atan_op_0( 00095 size_t i_z , 00096 size_t i_x , 00097 size_t nc_taylor , 00098 Base* taylor ) 00099 { 00100 // check assumptions 00101 CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 ); 00102 CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 ); 00103 CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z ); 00104 CPPAD_ASSERT_UNKNOWN( 0 < nc_taylor ); 00105 00106 // Taylor coefficients corresponding to argument and result 00107 Base* x = taylor + i_x * nc_taylor; 00108 Base* z = taylor + i_z * nc_taylor; 00109 Base* b = z - nc_taylor; // called y in documentation 00110 00111 z[0] = atan( x[0] ); 00112 b[0] = Base(1) + x[0] * x[0]; 00113 } 00114 /*! 00115 Reverse mode partial derivatives for result of op = AtanOp. 00116 00117 The C++ source code corresponding to this operation is 00118 \verbatim 00119 z = atan(x) 00120 \endverbatim 00121 The auxillary result is 00122 \verbatim 00123 y = 1 + x * x 00124 \endverbatim 00125 The value of y is computed along with the value of z. 00126 00127 \copydetails reverse_unary2_op 00128 */ 00129 00130 template <class Base> 00131 inline void reverse_atan_op( 00132 size_t d , 00133 size_t i_z , 00134 size_t i_x , 00135 size_t nc_taylor , 00136 const Base* taylor , 00137 size_t nc_partial , 00138 Base* partial ) 00139 { 00140 // check assumptions 00141 CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 ); 00142 CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 ); 00143 CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z ); 00144 CPPAD_ASSERT_UNKNOWN( d < nc_taylor ); 00145 CPPAD_ASSERT_UNKNOWN( d < nc_partial ); 00146 00147 // Taylor coefficients and partials corresponding to argument 00148 const Base* x = taylor + i_x * nc_taylor; 00149 Base* px = partial + i_x * nc_partial; 00150 00151 // Taylor coefficients and partials corresponding to first result 00152 const Base* z = taylor + i_z * nc_taylor; 00153 Base* pz = partial + i_z * nc_partial; 00154 00155 // Taylor coefficients and partials corresponding to auxillary result 00156 const Base* b = z - nc_taylor; // called y in documentation 00157 Base* pb = pz - nc_partial; 00158 00159 // number of indices to access 00160 size_t j = d; 00161 size_t k; 00162 while(j) 00163 { // scale partials w.r.t z[j] and b[j] 00164 pz[j] /= b[0]; 00165 pb[j] *= Base(2); 00166 00167 pb[0] -= pz[j] * z[j]; 00168 px[j] += pz[j] + pb[j] * x[0]; 00169 px[0] += pb[j] * x[j]; 00170 00171 // more scaling of partials w.r.t z[j] 00172 pz[j] /= Base(j); 00173 00174 for(k = 1; k < j; k++) 00175 { pb[j-k] -= pz[j] * Base(k) * z[k]; 00176 pz[k] -= pz[j] * Base(k) * b[j-k]; 00177 px[k] += pb[j] * x[j-k]; 00178 } 00179 --j; 00180 } 00181 px[0] += pz[0] / b[0] + pb[0] * Base(2) * x[0]; 00182 } 00183 00184 CPPAD_END_NAMESPACE 00185 # endif