CppAD: A C++ Algorithmic Differentiation Package 20110419
sqrt_op.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_SQRT_OP_INCLUDED
00003 # define CPPAD_SQRT_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 sqrt_op.hpp
00020 Forward and reverse mode calculations for z = sqrt(x).
00021 */
00022 
00023 
00024 /*!
00025 Compute forward mode Taylor coefficient for result of op = SqrtOp.
00026 
00027 The C++ source code corresponding to this operation is
00028 \verbatim
00029         z = sqrt(x)
00030 \endverbatim
00031 
00032 \copydetails forward_unary1_op
00033 */
00034 template <class Base>
00035 inline void forward_sqrt_op(
00036         size_t j           ,
00037         size_t i_z         ,
00038         size_t i_x         ,
00039         size_t nc_taylor   , 
00040         Base*  taylor      )
00041 {       
00042         // check assumptions
00043         CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
00044         CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
00045         CPPAD_ASSERT_UNKNOWN( i_x < i_z );
00046         CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
00047 
00048         // Taylor coefficients corresponding to argument and result
00049         Base* x = taylor + i_x * nc_taylor;
00050         Base* z = taylor + i_z * nc_taylor;
00051 
00052         size_t k;
00053         if( j == 0 )
00054                 z[j] = sqrt( x[0] );
00055         else
00056         {
00057                 z[j] = Base(0);
00058                 for(k = 1; k < j; k++)
00059                         z[j] -= Base(k) * z[k] * z[j-k];
00060                 z[j] /= Base(j);
00061                 z[j] += x[j] / Base(2);
00062                 z[j] /= z[0];
00063         }
00064 }
00065 
00066 /*!
00067 Compute zero order forward mode Taylor coefficient for result of op = SqrtOp.
00068 
00069 The C++ source code corresponding to this operation is
00070 \verbatim
00071         z = sqrt(x)
00072 \endverbatim
00073 
00074 \copydetails forward_unary1_op_0
00075 */
00076 template <class Base>
00077 inline void forward_sqrt_op_0(
00078         size_t i_z         ,
00079         size_t i_x         ,
00080         size_t nc_taylor   , 
00081         Base*  taylor      )
00082 {
00083         // check assumptions
00084         CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
00085         CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
00086         CPPAD_ASSERT_UNKNOWN( i_x < i_z );
00087         CPPAD_ASSERT_UNKNOWN( 0 < nc_taylor );
00088 
00089         // Taylor coefficients corresponding to argument and result
00090         Base* x = taylor + i_x * nc_taylor;
00091         Base* z = taylor + i_z * nc_taylor;
00092 
00093         z[0] = sqrt( x[0] );
00094 }
00095 /*!
00096 Compute reverse mode partial derivatives for result of op = SqrtOp.
00097 
00098 The C++ source code corresponding to this operation is
00099 \verbatim
00100         z = sqrt(x)
00101 \endverbatim
00102 
00103 \copydetails reverse_unary1_op
00104 */
00105 
00106 template <class Base>
00107 inline void reverse_sqrt_op(
00108         size_t      d            ,
00109         size_t      i_z          ,
00110         size_t      i_x          ,
00111         size_t      nc_taylor    , 
00112         const Base* taylor       ,
00113         size_t      nc_partial   ,
00114         Base*       partial      )
00115 {
00116         // check assumptions
00117         CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
00118         CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
00119         CPPAD_ASSERT_UNKNOWN( i_x < i_z );
00120         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
00121         CPPAD_ASSERT_UNKNOWN( d < nc_partial );
00122 
00123         // Taylor coefficients and partials corresponding to argument
00124         Base* px       = partial + i_x * nc_partial;
00125 
00126         // Taylor coefficients and partials corresponding to result
00127         const Base* z  = taylor  + i_z * nc_taylor;
00128         Base* pz       = partial + i_z * nc_partial;
00129 
00130         // number of indices to access
00131         size_t j = d;
00132         size_t k;
00133         while(j)
00134         {       // scale partial w.r.t. z[j]
00135                 pz[j]   /= z[0];
00136 
00137                 pz[0]   -= pz[j] * z[j];
00138                 px[j]   += pz[j] / Base(2);
00139                 for(k = 1; k < j; k++)
00140                         pz[k]   -= pz[j] * z[j-k];
00141                 --j;
00142         }
00143         px[0] += pz[0] / (Base(2) * z[0]);
00144 }
00145 
00146 CPPAD_END_NAMESPACE
00147 # endif