CppAD: A C++ Algorithmic Differentiation Package  20130918
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-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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines