source: branches/cache/cppad/local/sin_op.hpp @ 3324

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

merge trunk changes into cache

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