source: trunk/cppad/local/atan_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.0 KB
Line 
1/* $Id: atan_op.hpp 3301 2014-05-24 05:20:21Z bradbell $ */
2# ifndef CPPAD_ATAN_OP_INCLUDED
3# define CPPAD_ATAN_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 atan_op.hpp
20Forward and reverse mode calculations for z = atan(x).
21*/
22
23
24/*!
25Forward mode Taylor coefficient for result of op = AtanOp.
26
27The C++ source code corresponding to this operation is
28\verbatim
29        z = atan(x)
30\endverbatim
31The auxillary result is
32\verbatim
33        y = 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_atan_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(AtanOp) == 1 );
51        CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 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        if( p == 0 )
63        {       z[0] = atan( x[0] );
64                b[0] = Base(1) + x[0] * x[0];
65                p++;
66        }
67        for(size_t j = p; j <= q; j++)
68        {
69                b[j] = Base(2) * x[0] * x[j];
70                z[j] = Base(0);
71                for(k = 1; k < j; k++)
72                {       b[j] += x[k] * x[j-k];
73                        z[j] -= Base(k) * z[k] * b[j-k];
74                }
75                z[j] /= Base(j);
76                z[j] += x[j];
77                z[j] /= b[0];
78        }
79}
80
81/*!
82Multiple direction Taylor coefficient for op = AtanOp.
83
84The C++ source code corresponding to this operation is
85\verbatim
86        z = atan(x)
87\endverbatim
88The auxillary result is
89\verbatim
90        y = 1 + x * x
91\endverbatim
92The value of y, and its derivatives, are computed along with the value
93and derivatives of z.
94
95\copydetails forward_unary2_op_dir
96*/
97template <class Base>
98inline void forward_atan_op_dir(
99        size_t q           ,
100        size_t r           ,
101        size_t i_z         ,
102        size_t i_x         ,
103        size_t cap_order   , 
104        Base*  taylor      )
105{       
106        // check assumptions
107        CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 );
108        CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 );
109        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
110        CPPAD_ASSERT_UNKNOWN( 0 < q );
111        CPPAD_ASSERT_UNKNOWN( q < cap_order );
112
113        // Taylor coefficients corresponding to argument and result
114        size_t num_taylor_per_var = (cap_order-1) * r + 1;
115        Base* x = taylor + i_x * num_taylor_per_var;
116        Base* z = taylor + i_z * num_taylor_per_var;
117        Base* b = z      -       num_taylor_per_var; // called y in documentation
118
119        size_t m = (q-1) * r + 1;
120        for(size_t ell = 0; ell < r; ell++)
121        {       b[m+ell] = Base(2) * x[m+ell] * x[0];
122                z[m+ell] = Base(q) * x[m+ell];
123                for(size_t k = 1; k < q; k++)
124                {       b[m+ell] += x[(k-1)*r+1+ell] * x[(q-k-1)*r+1+ell];
125                        z[m+ell] -= Base(k) * z[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell];
126                }
127                z[m+ell] /= ( Base(q) * b[0] );
128        }
129}
130
131/*!
132Zero order forward mode Taylor coefficient for result of op = AtanOp.
133
134The C++ source code corresponding to this operation is
135\verbatim
136        z = atan(x)
137\endverbatim
138The auxillary result is
139\verbatim
140        y = 1 + x * x
141\endverbatim
142The value of y is computed along with the value of z.
143
144\copydetails forward_unary2_op_0
145*/
146template <class Base>
147inline void forward_atan_op_0(
148        size_t i_z         ,
149        size_t i_x         ,
150        size_t cap_order   , 
151        Base*  taylor      )
152{
153        // check assumptions
154        CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 );
155        CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 );
156        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
157        CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
158
159        // Taylor coefficients corresponding to argument and result
160        Base* x = taylor + i_x * cap_order;
161        Base* z = taylor + i_z * cap_order;
162        Base* b = z      -       cap_order; // called y in documentation
163
164        z[0] = atan( x[0] );
165        b[0] = Base(1) + x[0] * x[0];
166}
167/*!
168Reverse mode partial derivatives for result of op = AtanOp.
169
170The C++ source code corresponding to this operation is
171\verbatim
172        z = atan(x)
173\endverbatim
174The auxillary result is
175\verbatim
176        y = 1 + x * x
177\endverbatim
178The value of y is computed along with the value of z.
179
180\copydetails reverse_unary2_op
181*/
182
183template <class Base>
184inline void reverse_atan_op(
185        size_t      d            ,
186        size_t      i_z          ,
187        size_t      i_x          ,
188        size_t      cap_order    , 
189        const Base* taylor       ,
190        size_t      nc_partial   ,
191        Base*       partial      )
192{
193        // check assumptions
194        CPPAD_ASSERT_UNKNOWN( NumArg(AtanOp) == 1 );
195        CPPAD_ASSERT_UNKNOWN( NumRes(AtanOp) == 2 );
196        CPPAD_ASSERT_UNKNOWN( i_x + 1 < i_z );
197        CPPAD_ASSERT_UNKNOWN( d < cap_order );
198        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
199
200        // Taylor coefficients and partials corresponding to argument
201        const Base* x  = taylor  + i_x * cap_order;
202        Base* px       = partial + i_x * nc_partial;
203
204        // Taylor coefficients and partials corresponding to first result
205        const Base* z  = taylor  + i_z * cap_order;
206        Base* pz       = partial + i_z * nc_partial;
207
208        // Taylor coefficients and partials corresponding to auxillary result
209        const Base* b  = z  - cap_order; // called y in documentation
210        Base* pb       = pz - nc_partial;
211
212        // number of indices to access
213        size_t j = d;
214        size_t k;
215        while(j)
216        {       // scale partials w.r.t z[j] and b[j]
217                pz[j] /= b[0];
218                pb[j] *= Base(2);
219
220                pb[0] -= pz[j] * z[j]; 
221                px[j] += pz[j] + pb[j] * x[0];
222                px[0] += pb[j] * x[j];
223
224                // more scaling of partials w.r.t z[j]
225                pz[j] /= Base(j);
226
227                for(k = 1; k < j; k++)
228                {       pb[j-k] -= pz[j] * Base(k) * z[k];
229                        pz[k]   -= pz[j] * Base(k) * b[j-k];
230                        px[k]   += pb[j] * x[j-k];
231                }
232                --j;
233        }
234        px[0] += pz[0] / b[0] + pb[0] * Base(2) * x[0];
235}
236
237} // END_CPPAD_NAMESPACE
238# endif
Note: See TracBrowser for help on using the repository browser.