CppAD: A C++ Algorithmic Differentiation Package  20130102
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-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 asin_op_hpp asin_op.hpp
00020 \{
00021 \file asin_op.hpp
00022 Forward and reverse mode calculations for z = asin(x).
00023 */
00024 
00025 
00026 /*!
00027 Compute forward mode Taylor coefficient for result of op = AsinOp.
00028 
00029 The C++ source code corresponding to this operation is
00030 \verbatim
00031      z = asin(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_asin_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(AsinOp) == 1 );
00052      CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 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] = asin( 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 = AsinOp.
00091 
00092 The C++ source code corresponding to this operation is
00093 \verbatim
00094      z = asin(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_asin_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(AsinOp) == 1 );
00113      CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 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] = asin( x[0] );
00123      b[0] = sqrt( Base(1) - x[0] * x[0] );
00124 }
00125 /*!
00126 Compute reverse mode partial derivatives for result of op = AsinOp.
00127 
00128 The C++ source code corresponding to this operation is
00129 \verbatim
00130      z = asin(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_asin_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(AsinOp) == 1 );
00153      CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines