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