source: trunk/cppad/local/asin_op.hpp @ 2625

Last change on this file since 2625 was 2625, checked in by bradbell, 7 years ago

Change comment so end of group is included in doxygen processing
(do not know why works without this command at end).

ad_assign.hpp: fix doxygen brief description.
ad_ctor.hpp: fix doxygen brief description.

  • Property svn:keywords set to Id
File size: 5.1 KB
Line 
1/* $Id: asin_op.hpp 2625 2012-12-23 14:34:12Z bradbell $ */
2# ifndef CPPAD_ASIN_OP_INCLUDED
3# define CPPAD_ASIN_OP_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16
17CPPAD_BEGIN_NAMESPACE
18/*!
19\defgroup asin_op_hpp asin_op.hpp
20\{
21\file asin_op.hpp
22Forward and reverse mode calculations for z = asin(x).
23*/
24
25
26/*!
27Compute forward mode Taylor coefficient for result of op = AsinOp.
28
29The C++ source code corresponding to this operation is
30\verbatim
31        z = asin(x)
32\endverbatim
33The auxillary result is
34\verbatim
35        y = sqrt(1 - x * x)
36\endverbatim
37The value of y, and its derivatives, are computed along with the value
38and derivatives of z.
39
40\copydetails forward_unary2_op
41*/
42template <class Base>
43inline void forward_asin_op(
44        size_t j           ,
45        size_t i_z         ,
46        size_t i_x         ,
47        size_t nc_taylor   , 
48        Base*  taylor      )
49{       
50        // check assumptions
51        CPPAD_ASSERT_UNKNOWN( NumArg(AsinOp) == 1 );
52        CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 2 );
53        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
54        CPPAD_ASSERT_UNKNOWN( j < nc_taylor );
55
56        // Taylor coefficients corresponding to argument and result
57        Base* x = taylor + i_x * nc_taylor;
58        Base* z = taylor + i_z * nc_taylor;
59        Base* b = z      -       nc_taylor;  // called y in documentation
60
61        size_t k;
62        Base qj;
63        if( j == 0 )
64        {       z[j] = asin( x[0] );
65                qj   = Base(1) - x[0] * x[0];
66                b[j] = sqrt( qj );
67        }
68        else
69        {       qj = 0.;
70                for(k = 0; k <= j; k++)
71                        qj -= x[k] * x[j-k];
72                b[j] = Base(0);
73                z[j] = Base(0);
74                for(k = 1; k < j; k++)
75                {       b[j] -= Base(k) * b[k] * b[j-k];
76                        z[j] -= Base(k) * z[k] * b[j-k];
77                }
78                b[j] /= Base(j);
79                z[j] /= Base(j);
80                //
81                b[j] += qj / Base(2);
82                z[j] += x[j];
83                //
84                b[j] /= b[0];
85                z[j] /= b[0];
86        }
87}
88
89/*!
90Compute zero order forward mode Taylor coefficient for result of op = AsinOp.
91
92The C++ source code corresponding to this operation is
93\verbatim
94        z = asin(x)
95\endverbatim
96The auxillary result is
97\verbatim
98        y = sqrt( 1 - x * x )
99\endverbatim
100The value of y is computed along with the value of z.
101
102\copydetails forward_unary2_op_0
103*/
104template <class Base>
105inline void forward_asin_op_0(
106        size_t i_z         ,
107        size_t i_x         ,
108        size_t nc_taylor   , 
109        Base*  taylor      )
110{
111        // check assumptions
112        CPPAD_ASSERT_UNKNOWN( NumArg(AsinOp) == 1 );
113        CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 2 );
114        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
115        CPPAD_ASSERT_UNKNOWN( 0 < nc_taylor );
116
117        // Taylor coefficients corresponding to argument and result
118        Base* x = taylor + i_x * nc_taylor;
119        Base* z = taylor + i_z * nc_taylor;
120        Base* b = z      -       nc_taylor; // called y in documentation
121
122        z[0] = asin( x[0] );
123        b[0] = sqrt( Base(1) - x[0] * x[0] );
124}
125/*!
126Compute reverse mode partial derivatives for result of op = AsinOp.
127
128The C++ source code corresponding to this operation is
129\verbatim
130        z = asin(x)
131\endverbatim
132The auxillary result is
133\verbatim
134        y = sqrt( 1 - x * x )
135\endverbatim
136The value of y is computed along with the value of z.
137
138\copydetails reverse_unary2_op
139*/
140
141template <class Base>
142inline void reverse_asin_op(
143        size_t      d            ,
144        size_t      i_z          ,
145        size_t      i_x          ,
146        size_t      nc_taylor    , 
147        const Base* taylor       ,
148        size_t      nc_partial   ,
149        Base*       partial      )
150{
151        // check assumptions
152        CPPAD_ASSERT_UNKNOWN( NumArg(AsinOp) == 1 );
153        CPPAD_ASSERT_UNKNOWN( NumRes(AsinOp) == 2 );
154        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
155        CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
156        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
157
158        // Taylor coefficients and partials corresponding to argument
159        const Base* x  = taylor  + i_x * nc_taylor;
160        Base* px       = partial + i_x * nc_partial;
161
162        // Taylor coefficients and partials corresponding to first result
163        const Base* z  = taylor  + i_z * nc_taylor;
164        Base* pz       = partial + i_z * nc_partial;
165
166        // Taylor coefficients and partials corresponding to auxillary result
167        const Base* b  = z  - nc_taylor; // called y in documentation
168        Base* pb       = pz - nc_partial;
169
170        // number of indices to access
171        size_t j = d;
172        size_t k;
173        while(j)
174        {
175                // scale partials w.r.t b[j] by 1 / b[0]
176                pb[j] /= b[0];
177
178                // scale partials w.r.t z[j] by 1 / b[0]
179                pz[j] /= b[0];
180
181                // update partials w.r.t b^0
182                pb[0] -= pz[j] * z[j] + pb[j] * b[j]; 
183
184                // update partial w.r.t. x^0
185                px[0] -= pb[j] * x[j];
186
187                // update partial w.r.t. x^j
188                px[j] += pz[j] - pb[j] * x[0];
189
190                // further scale partial w.r.t. z[j] by 1 / j
191                pz[j] /= Base(j);
192
193                for(k = 1; k < j; k++)
194                {       // update partials w.r.t b^(j-k)
195                        pb[j-k] -= Base(k) * pz[j] * z[k] + pb[j] * b[k];
196
197                        // update partials w.r.t. x^k
198                        px[k]   -= pb[j] * x[j-k];
199
200                        // update partials w.r.t. z^k
201                        pz[k]   -= pz[j] * Base(k) * b[j-k];
202                }
203                --j;
204        }
205
206        // j == 0 case
207        px[0] += ( pz[0] - pb[0] * x[0]) / b[0];
208}
209
210/*! \} */
211CPPAD_END_NAMESPACE
212# endif
Note: See TracBrowser for help on using the repository browser.