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