CppAD: A C++ Algorithmic Differentiation Package 20110419
asin_op.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_ASIN_OP_INCLUDED
00003 # define CPPAD_ASIN_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 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 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(AsinOp) == 1 );
00050         CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 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] = asin( 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 = AsinOp.
00089 
00090 The C++ source code corresponding to this operation is
00091 \verbatim
00092         z = asin(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_asin_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(AsinOp) == 1 );
00111         CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 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] = asin( x[0] );
00121         b[0] = sqrt( Base(1) - x[0] * x[0] );
00122 }
00123 /*!
00124 Compute reverse mode partial derivatives for result of op = AsinOp.
00125 
00126 The C++ source code corresponding to this operation is
00127 \verbatim
00128         z = asin(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_asin_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(AsinOp) == 1 );
00151         CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 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