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