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