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