source: trunk/cppad/local/acos_op.hpp @ 3301

Last change on this file since 3301 was 3301, checked in by bradbell, 6 years ago

merge in multiple forward direcitons from branches/forward_dir

  • Property svn:keywords set to Id
File size: 6.7 KB
Line 
1/* $Id: acos_op.hpp 3301 2014-05-24 05:20:21Z bradbell $ */
2# ifndef CPPAD_ACOS_OP_INCLUDED
3# define CPPAD_ACOS_OP_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 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
17namespace CppAD { // BEGIN_CPPAD_NAMESPACE
18/*!
19\file acos_op.hpp
20Forward and reverse mode calculations for z = acos(x).
21*/
22
23
24/*!
25Compute forward mode Taylor coefficient for result of op = AcosOp.
26
27The C++ source code corresponding to this operation is
28\verbatim
29        z = acos(x)
30\endverbatim
31The auxillary result is
32\verbatim
33        y = sqrt(1 - x * x)
34\endverbatim
35The value of y, and its derivatives, are computed along with the value
36and derivatives of z.
37
38\copydetails forward_unary2_op
39*/
40template <class Base>
41inline void forward_acos_op(
42        size_t p           ,
43        size_t q           ,
44        size_t i_z         ,
45        size_t i_x         ,
46        size_t cap_order   , 
47        Base*  taylor      )
48{       
49        // check assumptions
50        CPPAD_ASSERT_UNKNOWN( NumArg(AcosOp) == 1 );
51        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 2 );
52        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
53        CPPAD_ASSERT_UNKNOWN( q < cap_order );
54        CPPAD_ASSERT_UNKNOWN( p <= q );
55
56        // Taylor coefficients corresponding to argument and result
57        Base* x = taylor + i_x * cap_order;
58        Base* z = taylor + i_z * cap_order;
59        Base* b = z      -       cap_order;  // called y in documentation
60
61        size_t k;
62        Base uj;
63        if( p == 0 )
64        {       z[0] = acos( x[0] );
65                uj   = Base(1) - x[0] * x[0];
66                b[0] = sqrt( uj );
67                p++;
68        }
69        for(size_t j = p; j <= q; j++)
70        {       uj = Base(0);
71                for(k = 0; k <= j; k++)
72                        uj -= x[k] * x[j-k];
73                b[j] = Base(0);
74                z[j] = Base(0);
75                for(k = 1; k < j; k++)
76                {       b[j] -= Base(k) * b[k] * b[j-k];
77                        z[j] -= Base(k) * z[k] * b[j-k];
78                }
79                b[j] /= Base(j);
80                z[j] /= Base(j);
81                //
82                b[j] += uj / Base(2);
83                z[j] -= x[j];
84                //
85                b[j] /= b[0];
86                z[j] /= b[0];
87        }
88}
89/*!
90Multiple directions forward mode Taylor coefficient for op = AcosOp.
91
92The C++ source code corresponding to this operation is
93\verbatim
94        z = acos(x)
95\endverbatim
96The auxillary result is
97\verbatim
98        y = sqrt(1 - x * x)
99\endverbatim
100The value of y, and its derivatives, are computed along with the value
101and derivatives of z.
102
103\copydetails forward_unary2_op_dir
104*/
105template <class Base>
106inline void forward_acos_op_dir(
107        size_t q           ,
108        size_t r           ,
109        size_t i_z         ,
110        size_t i_x         ,
111        size_t cap_order   , 
112        Base*  taylor      )
113{       
114        // check assumptions
115        CPPAD_ASSERT_UNKNOWN( NumArg(AcosOp) == 1 );
116        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 2 );
117        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
118        CPPAD_ASSERT_UNKNOWN( 0 < q );
119        CPPAD_ASSERT_UNKNOWN( q < cap_order );
120
121        // Taylor coefficients corresponding to argument and result
122        size_t num_taylor_per_var = (cap_order-1) * r + 1;
123        Base* x = taylor + i_x * num_taylor_per_var;
124        Base* z = taylor + i_z * num_taylor_per_var;
125        Base* b = z - num_taylor_per_var;  // called y in documentation
126
127        size_t k, ell;
128        size_t m = (q-1) * r + 1;
129        for(ell = 0; ell < r; ell ++)
130        {       Base uq = - 2.0 * x[m + ell] * x[0];
131                for(k = 1; k < q; k++)
132                        uq -= x[(k-1)*r+1+ell] * x[(q-k-1)*r+1+ell]; 
133                b[m+ell] = Base(0);
134                z[m+ell] = Base(0);
135                for(k = 1; k < q; k++)
136                {       b[m+ell] += Base(k) * b[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell]; 
137                        z[m+ell] += Base(k) * z[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell]; 
138                }
139                b[m+ell] =  ( uq / Base(2) - b[m+ell] / Base(q) ) / b[0];
140                z[m+ell] = -( x[m+ell]     + z[m+ell] / Base(q) ) / b[0];
141        }
142}
143
144/*!
145Compute zero order forward mode Taylor coefficient for result of op = AcosOp.
146
147The C++ source code corresponding to this operation is
148\verbatim
149        z = acos(x)
150\endverbatim
151The auxillary result is
152\verbatim
153        y = sqrt( 1 - x * x )
154\endverbatim
155The value of y is computed along with the value of z.
156
157\copydetails forward_unary2_op_0
158*/
159template <class Base>
160inline void forward_acos_op_0(
161        size_t i_z         ,
162        size_t i_x         ,
163        size_t cap_order   , 
164        Base*  taylor      )
165{
166        // check assumptions
167        CPPAD_ASSERT_UNKNOWN( NumArg(AcosOp) == 1 );
168        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 2 );
169        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
170        CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
171
172        // Taylor coefficients corresponding to argument and result
173        Base* x = taylor + i_x * cap_order;
174        Base* z = taylor + i_z * cap_order;
175        Base* b = z      -       cap_order; // called y in documentation
176
177        z[0] = acos( x[0] );
178        b[0] = sqrt( Base(1) - x[0] * x[0] );
179}
180/*!
181Compute reverse mode partial derivatives for result of op = AcosOp.
182
183The C++ source code corresponding to this operation is
184\verbatim
185        z = acos(x)
186\endverbatim
187The auxillary result is
188\verbatim
189        y = sqrt( 1 - x * x )
190\endverbatim
191The value of y is computed along with the value of z.
192
193\copydetails reverse_unary2_op
194*/
195
196template <class Base>
197inline void reverse_acos_op(
198        size_t      d            ,
199        size_t      i_z          ,
200        size_t      i_x          ,
201        size_t      cap_order    , 
202        const Base* taylor       ,
203        size_t      nc_partial   ,
204        Base*       partial      )
205{
206        // check assumptions
207        CPPAD_ASSERT_UNKNOWN( NumArg(AcosOp) == 1 );
208        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 2 );
209        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
210        CPPAD_ASSERT_UNKNOWN( d < cap_order );
211        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
212
213        // Taylor coefficients and partials corresponding to argument
214        const Base* x  = taylor  + i_x * cap_order;
215        Base* px       = partial + i_x * nc_partial;
216
217        // Taylor coefficients and partials corresponding to first result
218        const Base* z  = taylor  + i_z * cap_order;
219        Base* pz       = partial + i_z * nc_partial;
220
221        // Taylor coefficients and partials corresponding to auxillary result
222        const Base* b  = z  - cap_order; // called y in documentation
223        Base* pb       = pz - nc_partial;
224
225        // number of indices to access
226        size_t j = d;
227        size_t k;
228        while(j)
229        {
230                // scale partials w.r.t b[j] by 1 / b[0]
231                pb[j] /= b[0];
232
233                // scale partials w.r.t z[j] by 1 / b[0]
234                pz[j] /= b[0];
235
236                // update partials w.r.t b^0
237                pb[0] -= pz[j] * z[j] + pb[j] * b[j]; 
238
239                // update partial w.r.t. x^0
240                px[0] -= pb[j] * x[j];
241
242                // update partial w.r.t. x^j
243                px[j] -= pz[j] + pb[j] * x[0];
244
245                // further scale partial w.r.t. z[j] by 1 / j
246                pz[j] /= Base(j);
247
248                for(k = 1; k < j; k++)
249                {       // update partials w.r.t b^(j-k)
250                        pb[j-k] -= Base(k) * pz[j] * z[k] + pb[j] * b[k];
251
252                        // update partials w.r.t. x^k
253                        px[k]   -= pb[j] * x[j-k];
254
255                        // update partials w.r.t. z^k
256                        pz[k]   -= pz[j] * Base(k) * b[j-k];
257                }
258                --j;
259        }
260
261        // j == 0 case
262        px[0] -= ( pz[0] + pb[0] * x[0]) / b[0];
263}
264
265} // END_CPPAD_NAMESPACE
266# endif
Note: See TracBrowser for help on using the repository browser.