CppAD: A C++ Algorithmic Differentiation Package  20171217
div_op.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_LOCAL_DIV_OP_HPP
2 # define CPPAD_LOCAL_DIV_OP_HPP
3
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
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 namespace CppAD { namespace local { // BEGIN_CPPAD_LOCAL_NAMESPACE
16 /*!
17 \file div_op.hpp
18 Forward and reverse mode calculations for z = x / y.
19 */
20
21 // --------------------------- Divvv -----------------------------------------
22 /*!
23 Compute forward mode Taylor coefficients for result of op = DivvvOp.
24
25 The C++ source code corresponding to this operation is
26 \verbatim
27  z = x / y
28 \endverbatim
29 In the documentation below,
30 this operations is for the case where both x and y are variables
31 and the argument \a parameter is not used.
32
34 */
35
36 template <class Base>
37 inline void forward_divvv_op(
38  size_t p ,
39  size_t q ,
40  size_t i_z ,
41  const addr_t* arg ,
42  const Base* parameter ,
43  size_t cap_order ,
44  Base* taylor )
45 {
46  // check assumptions
49  CPPAD_ASSERT_UNKNOWN( q < cap_order );
50  CPPAD_ASSERT_UNKNOWN( p <= q );
51
52  // Taylor coefficients corresponding to arguments and result
53  Base* x = taylor + arg[0] * cap_order;
54  Base* y = taylor + arg[1] * cap_order;
55  Base* z = taylor + i_z * cap_order;
56
57
58  // Using CondExp, it can make sense to divide by zero,
59  // so do not make it an error.
60  size_t k;
61  for(size_t d = p; d <= q; d++)
62  { z[d] = x[d];
63  for(k = 1; k <= d; k++)
64  z[d] -= z[d-k] * y[k];
65  z[d] /= y[0];
66  }
67 }
68 /*!
69 Multiple directions forward mode Taylor coefficients for op = DivvvOp.
70
71 The C++ source code corresponding to this operation is
72 \verbatim
73  z = x / y
74 \endverbatim
75 In the documentation below,
76 this operations is for the case where both x and y are variables
77 and the argument \a parameter is not used.
78
80 */
81
82 template <class Base>
84  size_t q ,
85  size_t r ,
86  size_t i_z ,
87  const addr_t* arg ,
88  const Base* parameter ,
89  size_t cap_order ,
90  Base* taylor )
91 {
92  // check assumptions
95  CPPAD_ASSERT_UNKNOWN( 0 < q );
96  CPPAD_ASSERT_UNKNOWN( q < cap_order );
97
98  // Taylor coefficients corresponding to arguments and result
99  size_t num_taylor_per_var = (cap_order-1) * r + 1;
100  Base* x = taylor + arg[0] * num_taylor_per_var;
101  Base* y = taylor + arg[1] * num_taylor_per_var;
102  Base* z = taylor + i_z * num_taylor_per_var;
103
104
105  // Using CondExp, it can make sense to divide by zero,
106  // so do not make it an error.
107  size_t m = (q-1) * r + 1;
108  for(size_t ell = 0; ell < r; ell++)
109  { z[m+ell] = x[m+ell] - z[0] * y[m+ell];
110  for(size_t k = 1; k < q; k++)
111  z[m+ell] -= z[(q-k-1)*r+1+ell] * y[(k-1)*r+1+ell];
112  z[m+ell] /= y[0];
113  }
114 }
115
116
117 /*!
118 Compute zero order forward mode Taylor coefficients for result of op = DivvvOp.
119
120 The C++ source code corresponding to this operation is
121 \verbatim
122  z = x / y
123 \endverbatim
124 In the documentation below,
125 this operations is for the case where both x and y are variables
126 and the argument \a parameter is not used.
127
129 */
130
131 template <class Base>
132 inline void forward_divvv_op_0(
133  size_t i_z ,
134  const addr_t* arg ,
135  const Base* parameter ,
136  size_t cap_order ,
137  Base* taylor )
138 {
139  // check assumptions
142
143  // Taylor coefficients corresponding to arguments and result
144  Base* x = taylor + arg[0] * cap_order;
145  Base* y = taylor + arg[1] * cap_order;
146  Base* z = taylor + i_z * cap_order;
147
148  z[0] = x[0] / y[0];
149 }
150
151 /*!
152 Compute reverse mode partial derivatives for result of op = DivvvOp.
153
154 The C++ source code corresponding to this operation is
155 \verbatim
156  z = x / y
157 \endverbatim
158 In the documentation below,
159 this operations is for the case where both x and y are variables
160 and the argument \a parameter is not used.
161
163 */
164
165 template <class Base>
166 inline void reverse_divvv_op(
167  size_t d ,
168  size_t i_z ,
169  const addr_t* arg ,
170  const Base* parameter ,
171  size_t cap_order ,
172  const Base* taylor ,
173  size_t nc_partial ,
174  Base* partial )
175 {
176  // check assumptions
179  CPPAD_ASSERT_UNKNOWN( d < cap_order );
180  CPPAD_ASSERT_UNKNOWN( d < nc_partial );
181
182  // Arguments
183  const Base* y = taylor + arg[1] * cap_order;
184  const Base* z = taylor + i_z * cap_order;
185
186  // Partial derivatives corresponding to arguments and result
187  Base* px = partial + arg[0] * nc_partial;
188  Base* py = partial + arg[1] * nc_partial;
189  Base* pz = partial + i_z * nc_partial;
190
191  // Using CondExp, it can make sense to divide by zero
192  // so do not make it an error.
193  Base inv_y0 = Base(1.0) / y[0];
194
195  size_t k;
196  // number of indices to access
197  size_t j = d + 1;
198  while(j)
199  { --j;
200  // scale partial w.r.t. z[j]
201  pz[j] = azmul(pz[j], inv_y0);
202
203  px[j] += pz[j];
204  for(k = 1; k <= j; k++)
205  { pz[j-k] -= azmul(pz[j], y[k] );
206  py[k] -= azmul(pz[j], z[j-k]);
207  }
208  py[0] -= azmul(pz[j], z[j]);
209  }
210 }
211
212 // --------------------------- Divpv -----------------------------------------
213 /*!
214 Compute forward mode Taylor coefficients for result of op = DivpvOp.
215
216 The C++ source code corresponding to this operation is
217 \verbatim
218  z = x / y
219 \endverbatim
220 In the documentation below,
221 this operations is for the case where x is a parameter and y is a variable.
222
224 */
225
226 template <class Base>
227 inline void forward_divpv_op(
228  size_t p ,
229  size_t q ,
230  size_t i_z ,
231  const addr_t* arg ,
232  const Base* parameter ,
233  size_t cap_order ,
234  Base* taylor )
235 {
236  // check assumptions
239  CPPAD_ASSERT_UNKNOWN( q < cap_order );
240  CPPAD_ASSERT_UNKNOWN( p <= q );
241
242  // Taylor coefficients corresponding to arguments and result
243  Base* y = taylor + arg[1] * cap_order;
244  Base* z = taylor + i_z * cap_order;
245
246  // Paraemter value
247  Base x = parameter[ arg[0] ];
248
249  // Using CondExp, it can make sense to divide by zero,
250  // so do not make it an error.
251  size_t k;
252  if( p == 0 )
253  { z[0] = x / y[0];
254  p++;
255  }
256  for(size_t d = p; d <= q; d++)
257  { z[d] = Base(0.0);
258  for(k = 1; k <= d; k++)
259  z[d] -= z[d-k] * y[k];
260  z[d] /= y[0];
261  }
262 }
263 /*!
264 Multiple directions forward mode Taylor coefficients for op = DivpvOp.
265
266 The C++ source code corresponding to this operation is
267 \verbatim
268  z = x / y
269 \endverbatim
270 In the documentation below,
271 this operations is for the case where x is a parameter and y is a variable.
272
274 */
275
276 template <class Base>
278  size_t q ,
279  size_t r ,
280  size_t i_z ,
281  const addr_t* arg ,
282  const Base* parameter ,
283  size_t cap_order ,
284  Base* taylor )
285 {
286  // check assumptions
289  CPPAD_ASSERT_UNKNOWN( 0 < q );
290  CPPAD_ASSERT_UNKNOWN( q < cap_order );
291
292  // Taylor coefficients corresponding to arguments and result
293  size_t num_taylor_per_var = (cap_order-1) * r + 1;
294  Base* y = taylor + arg[1] * num_taylor_per_var;
295  Base* z = taylor + i_z * num_taylor_per_var;
296
297  // Using CondExp, it can make sense to divide by zero,
298  // so do not make it an error.
299  size_t m = (q-1) * r + 1;
300  for(size_t ell = 0; ell < r; ell++)
301  { z[m+ell] = - z[0] * y[m+ell];
302  for(size_t k = 1; k < q; k++)
303  z[m+ell] -= z[(q-k-1)*r+1+ell] * y[(k-1)*r+1+ell];
304  z[m+ell] /= y[0];
305  }
306 }
307
308 /*!
309 Compute zero order forward mode Taylor coefficient for result of op = DivpvOp.
310
311 The C++ source code corresponding to this operation is
312 \verbatim
313  z = x / y
314 \endverbatim
315 In the documentation below,
316 this operations is for the case where x is a parameter and y is a variable.
317
319 */
320
321 template <class Base>
322 inline void forward_divpv_op_0(
323  size_t i_z ,
324  const addr_t* arg ,
325  const Base* parameter ,
326  size_t cap_order ,
327  Base* taylor )
328 {
329  // check assumptions
332
333  // Paraemter value
334  Base x = parameter[ arg[0] ];
335
336  // Taylor coefficients corresponding to arguments and result
337  Base* y = taylor + arg[1] * cap_order;
338  Base* z = taylor + i_z * cap_order;
339
340  z[0] = x / y[0];
341 }
342
343 /*!
344 Compute reverse mode partial derivative for result of op = DivpvOp.
345
346 The C++ source code corresponding to this operation is
347 \verbatim
348  z = x / y
349 \endverbatim
350 In the documentation below,
351 this operations is for the case where x is a parameter and y is a variable.
352
354 */
355
356 template <class Base>
357 inline void reverse_divpv_op(
358  size_t d ,
359  size_t i_z ,
360  const addr_t* arg ,
361  const Base* parameter ,
362  size_t cap_order ,
363  const Base* taylor ,
364  size_t nc_partial ,
365  Base* partial )
366 {
367  // check assumptions
370  CPPAD_ASSERT_UNKNOWN( d < cap_order );
371  CPPAD_ASSERT_UNKNOWN( d < nc_partial );
372
373  // Arguments
374  const Base* y = taylor + arg[1] * cap_order;
375  const Base* z = taylor + i_z * cap_order;
376
377  // Partial derivatives corresponding to arguments and result
378  Base* py = partial + arg[1] * nc_partial;
379  Base* pz = partial + i_z * nc_partial;
380
381  // Using CondExp, it can make sense to divide by zero so do not
382  // make it an error.
383  Base inv_y0 = Base(1.0) / y[0];
384
385  size_t k;
386  // number of indices to access
387  size_t j = d + 1;
388  while(j)
389  { --j;
390  // scale partial w.r.t z[j]
391  pz[j] = azmul(pz[j], inv_y0);
392
393  for(k = 1; k <= j; k++)
394  { pz[j-k] -= azmul(pz[j], y[k] );
395  py[k] -= azmul(pz[j], z[j-k] );
396  }
397  py[0] -= azmul(pz[j], z[j]);
398  }
399 }
400
401
402 // --------------------------- Divvp -----------------------------------------
403 /*!
404 Compute forward mode Taylor coefficients for result of op = DivvvOp.
405
406 The C++ source code corresponding to this operation is
407 \verbatim
408  z = x / y
409 \endverbatim
410 In the documentation below,
411 this operations is for the case where x is a variable and y is a parameter.
412
414 */
415
416 template <class Base>
417 inline void forward_divvp_op(
418  size_t p ,
419  size_t q ,
420  size_t i_z ,
421  const addr_t* arg ,
422  const Base* parameter ,
423  size_t cap_order ,
424  Base* taylor )
425 {
426  // check assumptions
429  CPPAD_ASSERT_UNKNOWN( q < cap_order );
430  CPPAD_ASSERT_UNKNOWN( p <= q );
431
432  // Taylor coefficients corresponding to arguments and result
433  Base* x = taylor + arg[0] * cap_order;
434  Base* z = taylor + i_z * cap_order;
435
436  // Parameter value
437  Base y = parameter[ arg[1] ];
438
439  // Using CondExp and multiple levels of AD, it can make sense
440  // to divide by zero so do not make it an error.
441  for(size_t d = p; d <= q; d++)
442  z[d] = x[d] / y;
443 }
444 /*!
445 Multiple direction forward mode Taylor coefficients for op = DivvvOp.
446
447 The C++ source code corresponding to this operation is
448 \verbatim
449  z = x / y
450 \endverbatim
451 In the documentation below,
452 this operations is for the case where x is a variable and y is a parameter.
453
455 */
456
457 template <class Base>
459  size_t q ,
460  size_t r ,
461  size_t i_z ,
462  const addr_t* arg ,
463  const Base* parameter ,
464  size_t cap_order ,
465  Base* taylor )
466 {
467  // check assumptions
470  CPPAD_ASSERT_UNKNOWN( q < cap_order );
471  CPPAD_ASSERT_UNKNOWN( 0 < q );
472
473  // Taylor coefficients corresponding to arguments and result
474  size_t num_taylor_per_var = (cap_order-1) * r + 1;
475  Base* x = taylor + arg[0] * num_taylor_per_var;
476  Base* z = taylor + i_z * num_taylor_per_var;
477
478  // Parameter value
479  Base y = parameter[ arg[1] ];
480
481  // Using CondExp and multiple levels of AD, it can make sense
482  // to divide by zero so do not make it an error.
483  size_t m = (q-1)*r + 1;
484  for(size_t ell = 0; ell < r; ell++)
485  z[m + ell] = x[m + ell] / y;
486 }
487
488
489 /*!
490 Compute zero order forward mode Taylor coefficients for result of op = DivvvOp.
491
492 The C++ source code corresponding to this operation is
493 \verbatim
494  z = x / y
495 \endverbatim
496 In the documentation below,
497 this operations is for the case where x is a variable and y is a parameter.
498
500 */
501
502 template <class Base>
503 inline void forward_divvp_op_0(
504  size_t i_z ,
505  const addr_t* arg ,
506  const Base* parameter ,
507  size_t cap_order ,
508  Base* taylor )
509 {
510  // check assumptions
513
514  // Parameter value
515  Base y = parameter[ arg[1] ];
516
517  // Taylor coefficients corresponding to arguments and result
518  Base* x = taylor + arg[0] * cap_order;
519  Base* z = taylor + i_z * cap_order;
520
521  z[0] = x[0] / y;
522 }
523
524 /*!
525 Compute reverse mode partial derivative for result of op = DivvpOp.
526
527 The C++ source code corresponding to this operation is
528 \verbatim
529  z = x / y
530 \endverbatim
531 In the documentation below,
532 this operations is for the case where x is a variable and y is a parameter.
533
535 */
536
537 template <class Base>
538 inline void reverse_divvp_op(
539  size_t d ,
540  size_t i_z ,
541  const addr_t* arg ,
542  const Base* parameter ,
543  size_t cap_order ,
544  const Base* taylor ,
545  size_t nc_partial ,
546  Base* partial )
547 {
548  // check assumptions
551  CPPAD_ASSERT_UNKNOWN( d < cap_order );
552  CPPAD_ASSERT_UNKNOWN( d < nc_partial );
553
554  // Argument values
555  Base y = parameter[ arg[1] ];
556
557  // Partial derivatives corresponding to arguments and result
558  Base* px = partial + arg[0] * nc_partial;
559  Base* pz = partial + i_z * nc_partial;
560
561  // Using CondExp, it can make sense to divide by zero
562  // so do not make it an error.
563  Base inv_y = Base(1.0) / y;
564
565  // number of indices to access
566  size_t j = d + 1;
567  while(j)
568  { --j;
569  px[j] += azmul(pz[j], inv_y);
570  }
571 }
572
573 } } // END_CPPAD_LOCAL_NAMESPACE
574 # endif
AD< Base > azmul(const AD< Base > &x, const AD< Base > &y)
Definition: azmul.hpp:94
void forward_divpv_op_0(size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute zero order forward mode Taylor coefficient for result of op = DivpvOp.
Definition: div_op.hpp:322
void forward_divvp_op(size_t p, size_t q, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute forward mode Taylor coefficients for result of op = DivvvOp.
Definition: div_op.hpp:417
size_t NumArg(OpCode op)
Number of arguments for a specified operator.
Definition: op_code.hpp:175
void forward_divvp_op_dir(size_t q, size_t r, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Multiple direction forward mode Taylor coefficients for op = DivvvOp.
Definition: div_op.hpp:458
size_t NumRes(OpCode op)
Number of variables resulting from the specified operation.
Definition: op_code.hpp:281
void reverse_divpv_op(size_t d, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, const Base *taylor, size_t nc_partial, Base *partial)
Compute reverse mode partial derivative for result of op = DivpvOp.
Definition: div_op.hpp:357
void forward_divpv_op_dir(size_t q, size_t r, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Multiple directions forward mode Taylor coefficients for op = DivpvOp.
Definition: div_op.hpp:277
void reverse_divvp_op(size_t d, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, const Base *taylor, size_t nc_partial, Base *partial)
Compute reverse mode partial derivative for result of op = DivvpOp.
Definition: div_op.hpp:538
void forward_divvv_op_0(size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute zero order forward mode Taylor coefficients for result of op = DivvvOp.
Definition: div_op.hpp:132
void forward_divpv_op(size_t p, size_t q, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute forward mode Taylor coefficients for result of op = DivpvOp.
Definition: div_op.hpp:227