CppAD: A C++ Algorithmic Differentiation Package  20171217
acos_op.hpp
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
6
8 the terms of the
9  Eclipse Public License Version 1.0.
10
11 A copy of this license is included in the COPYING file of this distribution.
13 -------------------------------------------------------------------------- */
14
15
17 /*!
18 \file acos_op.hpp
19 Forward and reverse mode calculations for z = acos(x).
20 */
21
22
23 /*!
24 Compute forward mode Taylor coefficient for result of op = AcosOp.
25
26 The C++ source code corresponding to this operation is
27 \verbatim
28  z = acos(x)
29 \endverbatim
30 The auxillary result is
31 \verbatim
32  y = sqrt(1 - x * x)
33 \endverbatim
34 The value of y, and its derivatives, are computed along with the value
35 and derivatives of z.
36
38 */
39 template <class Base>
40 inline void forward_acos_op(
41  size_t p ,
42  size_t q ,
43  size_t i_z ,
44  size_t i_x ,
45  size_t cap_order ,
46  Base* taylor )
47 {
48  // check assumptions
51  CPPAD_ASSERT_UNKNOWN( q < cap_order );
52  CPPAD_ASSERT_UNKNOWN( p <= q );
53
54  // Taylor coefficients corresponding to argument and result
55  Base* x = taylor + i_x * cap_order;
56  Base* z = taylor + i_z * cap_order;
57  Base* b = z - cap_order; // called y in documentation
58
59  size_t k;
60  Base uj;
61  if( p == 0 )
62  { z[0] = acos( x[0] );
63  uj = Base(1.0) - x[0] * x[0];
64  b[0] = sqrt( uj );
65  p++;
66  }
67  for(size_t j = p; j <= q; j++)
68  { uj = Base(0.0);
69  for(k = 0; k <= j; k++)
70  uj -= x[k] * x[j-k];
71  b[j] = Base(0.0);
72  z[j] = Base(0.0);
73  for(k = 1; k < j; k++)
74  { b[j] -= Base(double(k)) * b[k] * b[j-k];
75  z[j] -= Base(double(k)) * z[k] * b[j-k];
76  }
77  b[j] /= Base(double(j));
78  z[j] /= Base(double(j));
79  //
80  b[j] += uj / Base(2.0);
81  z[j] -= x[j];
82  //
83  b[j] /= b[0];
84  z[j] /= b[0];
85  }
86 }
87 /*!
88 Multiple directions forward mode Taylor coefficient for op = AcosOp.
89
90 The C++ source code corresponding to this operation is
91 \verbatim
92  z = acos(x)
93 \endverbatim
94 The auxillary result is
95 \verbatim
96  y = sqrt(1 - x * x)
97 \endverbatim
98 The value of y, and its derivatives, are computed along with the value
99 and derivatives of z.
100
102 */
103 template <class Base>
105  size_t q ,
106  size_t r ,
107  size_t i_z ,
108  size_t i_x ,
109  size_t cap_order ,
110  Base* taylor )
111 {
112  // check assumptions
115  CPPAD_ASSERT_UNKNOWN( 0 < q );
116  CPPAD_ASSERT_UNKNOWN( q < cap_order );
117
118  // Taylor coefficients corresponding to argument and result
119  size_t num_taylor_per_var = (cap_order-1) * r + 1;
120  Base* x = taylor + i_x * num_taylor_per_var;
121  Base* z = taylor + i_z * num_taylor_per_var;
122  Base* b = z - num_taylor_per_var; // called y in documentation
123
124  size_t k, ell;
125  size_t m = (q-1) * r + 1;
126  for(ell = 0; ell < r; ell ++)
127  { Base uq = - 2.0 * x[m + ell] * x[0];
128  for(k = 1; k < q; k++)
129  uq -= x[(k-1)*r+1+ell] * x[(q-k-1)*r+1+ell];
130  b[m+ell] = Base(0.0);
131  z[m+ell] = Base(0.0);
132  for(k = 1; k < q; k++)
133  { b[m+ell] += Base(double(k)) * b[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell];
134  z[m+ell] += Base(double(k)) * z[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell];
135  }
136  b[m+ell] = ( uq / Base(2.0) - b[m+ell] / Base(double(q)) ) / b[0];
137  z[m+ell] = -( x[m+ell] + z[m+ell] / Base(double(q)) ) / b[0];
138  }
139 }
140
141 /*!
142 Compute zero order forward mode Taylor coefficient for result of op = AcosOp.
143
144 The C++ source code corresponding to this operation is
145 \verbatim
146  z = acos(x)
147 \endverbatim
148 The auxillary result is
149 \verbatim
150  y = sqrt( 1 - x * x )
151 \endverbatim
152 The value of y is computed along with the value of z.
153
155 */
156 template <class Base>
157 inline void forward_acos_op_0(
158  size_t i_z ,
159  size_t i_x ,
160  size_t cap_order ,
161  Base* taylor )
162 {
163  // check assumptions
166  CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
167
168  // Taylor coefficients corresponding to argument and result
169  Base* x = taylor + i_x * cap_order;
170  Base* z = taylor + i_z * cap_order;
171  Base* b = z - cap_order; // called y in documentation
172
173  z[0] = acos( x[0] );
174  b[0] = sqrt( Base(1.0) - x[0] * x[0] );
175 }
176 /*!
177 Compute reverse mode partial derivatives for result of op = AcosOp.
178
179 The C++ source code corresponding to this operation is
180 \verbatim
181  z = acos(x)
182 \endverbatim
183 The auxillary result is
184 \verbatim
185  y = sqrt( 1 - x * x )
186 \endverbatim
187 The value of y is computed along with the value of z.
188
190 */
191
192 template <class Base>
193 inline void reverse_acos_op(
194  size_t d ,
195  size_t i_z ,
196  size_t i_x ,
197  size_t cap_order ,
198  const Base* taylor ,
199  size_t nc_partial ,
200  Base* partial )
201 {
202  // check assumptions
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* z = taylor + i_z * cap_order;
214  Base* pz = partial + i_z * nc_partial;
215
216  // Taylor coefficients and partials corresponding to auxillary result
217  const Base* b = z - cap_order; // called y in documentation
218  Base* pb = pz - nc_partial;
219
220  Base inv_b0 = Base(1.0) / b[0];
221
222  // number of indices to access
223  size_t j = d;
224  size_t k;
225  while(j)
226  {
227  // scale partials w.r.t b[j] by 1 / b[0]
228  pb[j] = azmul(pb[j], inv_b0);
229
230  // scale partials w.r.t z[j] by 1 / b[0]
231  pz[j] = azmul(pz[j], inv_b0);
232
233  // update partials w.r.t b^0
234  pb[0] -= azmul(pz[j], z[j]) + azmul(pb[j], b[j]);
235
236  // update partial w.r.t. x^0
237  px[0] -= azmul(pb[j], x[j]);
238
239  // update partial w.r.t. x^j
240  px[j] -= pz[j] + azmul(pb[j], x[0]);
241
242  // further scale partial w.r.t. z[j] by 1 / j
243  pz[j] /= Base(double(j));
244
245  for(k = 1; k < j; k++)
246  { // update partials w.r.t b^(j-k)
247  pb[j-k] -= Base(double(k)) * azmul(pz[j], z[k]) + azmul(pb[j], b[k]);
248
249  // update partials w.r.t. x^k
250  px[k] -= azmul(pb[j], x[j-k]);
251
252  // update partials w.r.t. z^k
253  pz[k] -= Base(double(k)) * azmul(pz[j], b[j-k]);
254  }
255  --j;
256  }
257
258  // j == 0 case
259  px[0] -= azmul( pz[0] + azmul(pb[0], x[0]), inv_b0);
260 }
261
263 # endif
void forward_acos_op_dir(size_t q, size_t r, size_t i_z, size_t i_x, size_t cap_order, Base *taylor)
Multiple directions forward mode Taylor coefficient for op = AcosOp.
Definition: acos_op.hpp:104
Definition: azmul.hpp:94
size_t NumArg(OpCode op)
Number of arguments for a specified operator.
Definition: op_code.hpp:175
void forward_acos_op(size_t p, size_t q, size_t i_z, size_t i_x, size_t cap_order, Base *taylor)
Compute forward mode Taylor coefficient for result of op = AcosOp.
Definition: acos_op.hpp:40
size_t NumRes(OpCode op)
Number of variables resulting from the specified operation.
Definition: op_code.hpp:281
std::complex< double > acos(const std::complex< double > &x)
void reverse_acos_op(size_t d, size_t i_z, size_t i_x, size_t cap_order, const Base *taylor, size_t nc_partial, Base *partial)
Compute reverse mode partial derivatives for result of op = AcosOp.
Definition: acos_op.hpp:193